对称结构有限元分析

1、对称结构有限元分析3节点三角形单元的分析一问题分析（对称框架线弹性实体的静力平衡问题）图是一个方形弹性实体，单位边长、单位厚度、承受等效竖向压力1KN m2，其中边界条件暗示着存在两组相对称的平面，因此现考虑的仅是问题的J., -：-o每个节点上的自由度号码代表了各自在x和y方向上可能的位移。结构和单元信息NELSNCENNNIP8291BBAA约束节点自由度信息NR5K , NF(:,K),1=1, NR101 4 017 00 8 1910载荷信息LOADED_NODES3(K, LOADS(NF(:,K),1=1,LOADED_NODES)1 .0-.252.0-.53.0-.25.5.

2、551.E6.33333节点三角形单元网络的总体节点和单元编号3节三角形单元局部坐标系中节点和自由度编号二理论基础（有限元方法原理）通过弹性力学变分原理建立弹性力学问题有限元方法表达格式的基本步骤。最小位能原理的未知场变量是位移，以结点位移为基本未知量，并以最小位能原理为基础建立的有限元 为位移元。它是有限元方法中应用最为普遍的单元，也是本书主要讨论的单元。对于一个力学或无力问题，在建立其数学模型以后，用有限元方法对它进行分析的首要 步骤是选择单元形式。平面问题 3结点三角形单元是有限元方法最早采用，而且至今仍经常 采用的单元形式。我们将以它作为典型，讨论如何应用广义坐标建立单元位移模式与位移

3、插 值函数，以及如何根据最小位能原理建立有限元求解方程的原理、方法与步骤，并进而引出弹性力学问题有限元方法的一般表达格式。对于前一问题，着重讨论选择广义坐标和有限元位移模式的一般原则和建立其位移插值函数的一般步骤。对于后一问题，着重讨论单元刚度矩阵和单元载荷向量的形式，总体刚度矩阵和总体载荷向量集成的原理和方法，以及它们各自的特性。作为一种数值方法，有限元解的收敛性无疑是十分重要的问题，以后将讨论解的收敛准则及其物理意义，所阐明的原则在以后还将得到进一步的应用和具体化。在建立了有限元的一般表达格式以后，原则上可以将它推广到平面问题以外的其他弹性力学问题和采用任何形式的单元。轴对称问题具有很广泛

4、的应用领域，轴对称问题3结点三角形单元的表达格式可以看作是平面问题此种单元表达格式的直接推广。）弹性力学平面问题的有限元格式结点三角形单元是有限元方法中最早提出，并且至今仍广泛应用的单元，由于三角形单元对复杂边界有较强的适应能力，因此很容易将一个二维离散成有限个三角形单元，如图1所示。在边界上以若干段直线近似原来的曲线边界，随着单元增多，这种拟合将趋于精确。 我们在讨论如何应用有限元方法分析各类具体问题的开始，将以平面问题3结点三角形单元为例来阐明弹性力学问题有限元分析的表达格式和一般步 1.1 ）单元位移模式及插值函数的构造典型的3节点三角形单元节点编码i,j,m ，以逆时针方向编码为正向。

5、每个节点有位移分量如图所示。每个单元有6个节点位移即6个节点自由度，亦即ajUiViu j VjUmVm1.2） 单元的位移模式和广义坐标在有限元方法中单元的位移模式或称位移函数一般采用多项式作为近似函数，因为 多项式运算简便，并且随着项数的增多，可以逼近任何一段光滑的函数曲线。多项式的选 取由低次到高次。(2)13结点三角形单元位移模式选取一次多项式V =：4：5X 7它的矩阵表达式是其中1 x y I - 1-1：2称为位移模式，它表示位移作为坐标x，y的函数中所包含的项次。对于现在的情况,单元内的位移是坐标x，y的线性函数；-6是待定系数，称之为广义坐标。6个广义坐标可由单元的6个结点位

6、移来表示。在（1 ）的1式中代入结点i的坐标（x，y）可得到结点i 在x方向的位移叫，同理可得和。它们表示为Ui 二;-；1 jXi - ,yi TOC o 1-5 h z Uj = /Xj ij(3)Um 八1 FmMm解（2）式可以得到广义坐标有结点位移表示的表达式。上式的系数行列式是XiyiXj yj=2A(4)Xm Xm其实A是三角形单元的面积。广义坐标爲为UiXiyiXjyj2! aiUi ajUj amUmUmym1UiyiUjyj1云(biUi bjUjbmUm)(5)UmymXiUiXjUj1云(CiUi CjUjCmUm)XmUm同理，利用3个结点y方向的位移,即（a）式的第

7、2式可求得56在（c）式和（d）式中aibiCi1(ay2A1 ,(b Vi2A1(C Vi2AXjyj-ajVjCjVj-amVm)-bmVm)XmymyjymXjXm上式（i, j， m表示下标轮换，如二 XjymXm yj-ym ( i, j,m)-Xmj t m, nri。以下同此。1.3）位移插值函数 将求得的广义坐标代入（1）式，可将位移函数表示成结点位移的函数，即U = N iUiN jU j N mUmV 二 NjVjN jVjN mVm（8）其中1 Ni佝 biX y） （i,j,m）（9）2ANi,Nj,Nm称为单元的插值函数或形函数，对于当前情况，他的坐标x、y的一次函数

8、。其中的 ai,bi ,ci, cm是常熟，取决于单元的3个结点坐标。g）式中的单元面积 A可通过（7）的系数表示为1（ai2aj1 am )(bi C j2_ bjCi)(10)（f ）式的矩阵形式是Ui!J】NiNj00 NmNm3iUjVjUmN称为插值函数矩阵或形函数矩阵,jN iIN jIN mNiN j Nmaj丿e=Na（ii）称为单元结点位移列阵。（1）插值函数具有如下性质j -ij - i在结点上插值函数的值有N,Xj,yj)二(12)即有 N i(Xj, yj), Ni(Xj,yj)二 N j(xm，ym)也就是说在i结点上结点上N i =0。由（8）式可见，当x = Xj

9、, y = y j即在结点i，应有u =比，因此也必然要N m =0。其他两个形函数也具有同样的性质。此性质称为kronecker delta性质。(13)（2）在单元中任一点各插值函数之和应等于1，即Ni N j Nm 1因为若单元发生刚体位移，如X方向有刚体位移 J0，则单元内（包括结点上）到处应有位移巴，即ui =Uj =Um =u0，又由（g）式得到UrNjUj NjUj NmUm=(Ni N j Nm)Uo=Uo(17)因此必然要求NiNjNm=1。若插值函数不满足此要求，则不能反映单元的刚体位移,用以求解必然得不到的正确的结果。 单元的各个结点位移插值函数之和等于1的性质称为规一性

10、。（3）对于现在的单元，插值函数是线性的，在单元内部及单元的边界上位移也是线性的，可由结点上的位移值唯一地确定。由于相邻单元公共结点的结点位移是相等的，因此保证了相邻单元在公共边界上位移的连续性。1 . 4）应变矩阵和应力矩阵确定了单元位移后，可以很方便地利用几何方程和物理方程求得单元的应变和应力。 在几何方程中，位移用（11）式代入，得到单元应变为Bi二 LNa e 二 L N if eea BaNjI( 14)B称为应变矩阵，L是平面问题的微分算子。 应变矩阵B的分块子矩阵是Bi = LN i =00NiI cNijNicyo1jNi;:N i(15)e i1cN1=bi-cex2Acy2

11、 A代入（15）式得到fbi0 1Bi10Ci(i,j,m)2ACibi 一对（9）式求导可得3结点单元的应变矩阵是(16)(18)v(23)biB - BiB j Bm 丨2Ao b.Jo 5 oo bm oCj o Cmbi Cjb j Cm b m式中 ,bj,bm,Ci ,Cj ,cm 由(7)式确定，它们是单元形状的参数。当单元的结点坐标确定后，这些参数都是常量（与坐标变量X，y无关），因此B是常量阵。当单元的结点位移ae确定后，由B转换求得的单元应变都是常量，也就是说在载荷作用下单元中各点具有同样的值、 ；y值及xy 值。因此3结点三角形单元称为常应变单元。在应变梯度较大（也xy即

12、应力梯度较大）的部位，单元的划分应适当密集，否则将不能反映应变的真实变化而导致 较大的误差。单元应力可以根据物理方程求得，即在物理方程中代入（14）式可以得到-xy其中e二 D ；二 DBa Sa(19)S = DB 二D Bj Bm -Si Sj Sm 1(20)S称为应力矩阵。将平面应力或平面应变的弹性矩阵及（ 得到计算平面应力或平面应变问题的单元应力矩阵。rbi18）式代入（2o）S的分块矩阵为1ivo C i式，可以Si二 DB i2(1Eo2-vo ) AVobi-v oCi2Ci1 - V0(i, j, m)(21)bi其中、o，讥为材料常数。对于平面应力问题EoV。(22)对于平

13、面应变问题E与应变矩阵B相同，应力矩阵S也是常量阵，即3结点单元中各点的应力是相同的。 情况下，不单独定义应力矩阵 S,而直接用DB进行应力计算。EoVo 二1 -V在很多(31)）利用最小位能原理建立有限元方程最小位能原理的泛函数总位能的表达式，在平面问题中的矩阵表达形式为lipD ；tdxdyu ftdxdy-uTTtdSV s T(24)其中，t是二维体厚度；f是作用在二维体内的体积力；T是作用在二维体边界上的面积力。对于离散模型，系统位能是各个单元位能的和，利用（24）式并代入（11）及（14）式，及得到离散模型的总位能为eeT-p _p = (a.ee-eTT(a N ftdxdy

14、) (aT_B2eTeDBtdxdya )-TTtdS )(25)将结构总位能的各项矩阵表达成各个单元总位能的各对应项矩阵之和，隐含着要求单元各 项矩阵的阶系数（即单元的结点自由度数） 和结果各项矩阵的阶数 （即结构的结点自由度数） 相同。为此需要引入单元结点自由度和结构结点自由度的转换矩阵G,从而将单元结点位移列阵G用结构结点列阵ae表示，即a二 Ga(26)其中a =U1V1u2 v2 Ui Vi unVn T( 27)其中n为结构结点数。令Ke=btDBtdxdyPfe = N T ftdxdyeTeeePseNTtdSP二 PfPs（28）CJ0)read(10,*)(k,nf(:,k

15、),i=1,nr)call formnf (nf);neq=maxval(nf)nband=0! this is a plane stress analysisdee=.0; dee(1,1)=e/(1.-v*v);dee(2,2)=dee(1,1);dee(3,3)=.5*e/(1.+v) dee(1,2)=v*dee(1,1);dee(2,1)=dee(1,2); call sample(element,points,weights)elements_1: do iel=1,nelscall geometry_3tx(iel,nce,aa,bb,coord,num);call num_to

16、_g(num,nf,g) g_num(:,iel)=num;g_coord(:,num)=transpose(coord);g_g(:,iel)=g if(nbandbandwidth(g)nband=bandwidth(g)end do elements_1write(11, (a) “Global coordinates “do k=1,nn;write(11, (a,i5,a,2e12.4) )”Node”,k,” “,g_coord(:,k);end do write(11, (a) “Global node numbers “do k=1, nels; write(11, (a,15

17、,a,3i5) )&“Element “,k,” “,g_num(:,k); end dowrite(11, (2(a,i5) )“There are “,neq,”eduations and the half-bandwidth is “,nband allocate(kv(neq*(nband+1),loads(0:neq); kv=.0elements_2: do iel=1,nelsnum=g_num(:,iel); g=g_g( : ,iel )coord=transpose(g_coordd(:, num); km=0.0gauss_pts_1: do i=1, nipcall s

18、hape_der(der,points,i); jac=matmul(der,coord) det=determinant(jac); call invert(jac) deriv=matmul(jac,der); call beemat(bee,deriv) km=km+matmul(matmul(transpose(bee),dee),bee) *det* weights(i) end do gauss_pts_1call formkv (kv,km,g,neq)end do elements_2load=.0; read(10,*)loaded_nodes,(k,loads(nf(:,k

19、),i=1,loaded_nodes)call banred(kv,neq);call bacsub(kv,loads)write(11, (a) “The nodal displacements Are :”write(11, (a)”Node Displacement ”do k=1,nn; write(11, (i5,a,2e12.4) ) k,” “,loads(nf(:,k); end donip=1; deallocate(points,weights);allocate(points(nip,ndim),weights(nip) call sample(element,point

20、s,weights)write(11, (a) “The central point stresses are: ” elements_3:do iel=1,nelswrite(11, (a,i5) “Element No. “,ielnum=g_num(:,iel); coord=transpose(g_coord(: ,num)g=g_g( :,iel); eld=loads(g) gauss_pts_2: do i=1,nipcall shape_der(der,points,i);jac=matmul(der,coord)call invert(jac); deriv=matmul(j

21、ac,der)call beemat(bee,deriv); sigma=matmul(dee,matmul(bee,eld) write(11, (a,i5) “Point “,i ; write(11, (3e12.4) ) sigma end do gauss_pts_2end do elements_3end program p50五程序运行结果GlobalcoordinatesNode10.0000E+000.0000E+00Node20.5000E000.0000E+00Node30.1000E+010.0000E+00Node40.0000E+00-0.5000E00Node50

22、.5000E00-0.5000E00Node60.1000E+01-0.5000E00Node70.0000E+00-0.1000E+01Node80.5000E00-0.1000E+01Node90.1000E+01-0.1000E+01Global node numbersElement1124Element2542Element3235Element4453Element5457Element6875Element7568Element8986There are 12 equations and the half-bandwidth is 6 The nodal displacements are :NodeDisplacement10.0000E+00-0.1000E-0520.1500E-06-0.1000E-0530.3000E-06-0.1000E-0540.0000E+00-0.50