有限元讲义弹性力学平面问题有限单元法(四结点四边形等参元八结点曲线四边形等参元问题补充)分析

/《有限元》讲义2.6四结点四边形单元(Thefour-nodequadrilateralelement)前面介绍了四结点的矩形单元其位移函数:UV

12x3y4xy56x7y8xy为双线性函数，应力，应变在单元内呈线性变化，比常应力三角形单元精度高。但它对界限要求严格。本节介绍的四结点四边形等参元，它不单拥有较高的精度，并且其网格区分也不受界限的影响。对随意四边形单元（图见下边）若仍直接采纳前面矩形单元的位移函数，在界限上它便不再是线性的（因界限不与x,y轴一致），这样会使得相邻两单元在公共界限上的位移可能会出现不连续现象（非协调元），而使收敛性遇到影响。可以考证，利用坐标变换就能解决这个问题，即可以经过坐标变换将整体坐标中的四边形（图a）变换成在局部坐标系中与四边形方向没关的边长为2的正方形。正方形四个结点i,j,m,p按反时钟次序对应四边形的四个结点ijmp。正方形的1和1二条界限，分别对应四边形的i，j界限和p,m界限；ξ=-1和ξ=+1分别对应四边形的i，p界限和j，m界限。假如用二组直线均分四边形的四个界限限段，使四边形绘成一个非正交网格，那么该非正交网格在正方形上对应着一个等距离的规则网格（见图a,b）。自然,局部坐标上的A点与整体坐标的A点对应。1《有限元》讲义一、四结点四边形等参单元的形函数及坐标变换因为可以将整体坐标下的四边形单元变换成局部坐标下的正方形单元，关于这种正方形单元，自然仍取形函数为：UV

12325678引入界限条件，即可得位移函数：UNiUiijmpVNiViijmp写成矩阵形式：f

e

UNi0Np0eNdeV0Ni0dNp式中形函数:Ni,11i1ii,j,m,p4依据等参元的定义,我们将坐标变换式亦取为:xNixiNixiNjxjNmxmNpxpijmpyNiyiNiyiNjyjNmymNpyp261ijmp式中形函数N与位移函数中的完满一致。可以考证，利用坐标变换式(2-6-1)，可以把整体坐标系中的随意四边形单元（图a）变换成在局部坐标系中与四边形对应的边长为２的正方形。所以可以将上述位移函数和形函数用于随意四边形单元，并将形函数中的ξ,η理解为随意四边形单元的局部坐标。这样由位移函数可以获取单元各点的位移。在四条界限上分别有ξ＝±和η＝±１，故界限上的位移呈线性变化，位移的连续性可获取保证。于是，我们可以理解为：随意四边形单元是从基本的正方形单元变换过来的实质单元。所以又称正方形单元为母体单元，或基本单元。例题：为了加深理解，现观察实质单元为矩形单元的坐标变换，在2.4节中，我们定义局部坐标与整体坐标的关系是：2《有限元》讲义1xx01yy0ab式中(x0,y0)为局部坐标原点。由上第一式1x0得：xaxax1xx1xx02ji2ij将其从头组合：x11x11xj2i2111x111xj111xm111xp4i444比较2.4中的形函数表达式,便知:xNixiNjxjNmxmNpxp自然同理可得:yNiyiNjyjNmymNpyp由此知,矩形单元可以看作是四结点四边形单元的特例，自然，它也是等参元。《有限元法概论》（第二版）P172中，是这样解说等参元的基本见解和推导方法的：图形变换四结点正方形(母元)图形变换四结点四边形(等参元)(ξ-η平面内)───→(x,y平面内)进行图形变换的重点是进行图形结点坐标之间的变换：正方形结点坐标坐标变换四边形结点坐标(ξi,ηi)────→(x,y)i=i,j,m,pi=i,j,m,p为了实现上述结点坐标之间的变换，可利用母元的形函数，得出(ξ,η)和(x,y)之间的坐标变换式。图形变换拥有以下性质:母元中的坐标线对应于等参元的直线；四结点正方形母元对应于四个结点可以随意部署的直边四边形等参元；变换式(2-6-1)能保证相邻等参元的界限位移相互协调。3《有限元》讲义二、几何矩阵[B]已知单元的应变与结点位移之间的关系是:x00Np00Ni262yNi0dy0Npx形函数矩阵[N]但是局部坐标ξ,η的显函数,为求形函数对整体坐标x,y的偏导数，必然用复合函数求导公式:NiNixNiyxyNiNixNiy263xyNiNi或写成：Jx263aNiNiyxy式中:Jy263bxx称为雅可比矩阵，而把它的队列式称为雅可比队列式。把式(2-6-1)代入[J]得:J

NixiNiNiNjNmNpxiyiyixjyjijmpijmpNixiNiNiNjNmNpxmymyixpypijmpijmp将形函数Ni141i1ii,j,m,p代入，分别对,求偏导，即可得到四结点四边形等参元的雅可比矩阵：J

1143

A2B65A24B式中常数记为:Aiixiijmpixiijmp

Biiyiijmpiyiijmp4《有限元》讲义3ixi4iyiijmpijmp该雅可比矩阵的逆:114B2BJ4JAA31雅可比队列式:J1A2A143161B1A2A4142316

AB3可以证明，假如四结点四边形的四个内角都小于180°的话，雅可比队列式|J|大于零，其逆阵[J]-1是存在的。换句话说，为了使上述等参元能保持较好的精度，整体坐标系下所划分的随意四边形单元必然是凸四边形，即随意内角都不可以大于180°。四边形也不可以太倾斜，不然会影响其精度。利用雅可比的逆矩阵，即可求出整体坐标系下形函数的偏导数：NiNix1Ji11Ni4iNiyi

11

i266i(i,j,m,p)求出所有偏导，即代回（2-6-2）右边，即可获取几何矩阵[B],[B]是,的函数,即:Ni0Nj0Np00xxxxNiNjNpB0yN0y00yyyxNiNiNjNjNpNpyxyxyx将(2-6-6)代入即可获取[B],[B]是,的函数。三、单元刚度矩阵获取[B]后,即可由单刚的一般表达式:TKtBDBdxdy求出四结点四边形的单元刚度矩阵。在按上述公式作积分运算时,必然把面积元dxdy变换成dd,图a上的面积元abdc的面积等于矢量ab与矢量ac的矢量积的模，5《有限元》讲义即微元dAabacxy沿ξ轴对应于dξ的矢量增量是:abdidj沿η轴对应于η的矢量增量是:acxdiydj式中i,j是坐标x,y的单位矢,注意到:iijj0ij1则有:dAxdiydjxdiydjxyyxddijJdd所以刚度矩阵的积分式:11DJtddKBT26711在计算单元刚度矩阵[K]中元素时，因为被积函数中出现了雅可比队列式，使得它用分析法很难求其积分，故常采纳高斯数值积分法.6《有限元》讲义四、数值积分bd1．一维数值积分Fa基本思想:结构一个多项目式，使在i(i=1,2n)上有iFi，bd来近似原被积函数Fb今后用近似函数的积分a的积分Fd。i称为积a分点或取样点，积分点i的数值和地点决定了近似F的程度，亦即决定数值积分的精度。关于n个积分点,依据积分点地点的不一样样选择，平常采纳两种不一样样的数值积分方法，Newton-Cotes积分和高斯积分方案。两者方法基真相同，但是前者的积分点ξi是等间距散布,今后者不是等间距散布。高斯积分的积分点地点由下述方法确立：①定义n次多项式P12n②由以下条件确立n个积分点地点bia由上二式可见，

Pd0i0,1n1有以下性质：①在积分点上Pi0；②多项式P与0,1,2,...n1在（a,b）域内正交。所以可知n个积分点的地点0,1,2,...n1正交的n次多项式P构成方程ξi是在求积域（a,b）内与biPd0的解。a于是,被积函数F可由2n-1次多项式来近似。nn1lin1FiiPi1i0llii1i2ni是n-1阶Lagrange插值函数1nbb用ad近似aFd,并考虑上边确立n个积分点地点的商定条件得:7《有限元》讲义bndHiFiRFai1bn1式中HilidaHi称为积分的权系数(加权系数,权)，可见加权系数H与被积函数F没关，只与积分i点个数和地点相关。为便于计算积分点地点ξi和加权系数Hi,常将上式中的积分限规格化，即令:a1b1由此计算出的i,Hi对应于原积分域(a,b)的关系为:ababibaHi222关于多重积分，按重积分规则，计算内层时，保持外层为常数，逐层计算即可。相关数值积分更详细的资料可参阅其教科书。．等参元计算中数值积分阶次的选择数值积分中，怎样选择积分阶次将直接影响计算精度、工作量，选择不妥则有可能使计算失败。选择积分阶次的原则第一是要保证积分的精度能知足所求问题的要求。如一维问题：设：插值函数中的多项式阶数为p，微分算子中导数的最高阶为m，则有限元获取的近似能量是2(p-m)次多项式。若被积函数是2(p-m)次多项式，应选高斯积分的阶次为n=p-m+1。三结点三角形单元（线性单元）刚度矩阵中被积函数是常数，故只要一个积分点。双线性单元22p-m+1=1,但插值函数将包括有二次项（ξ,ξη,η）中的ξη项，所以要达到精准积分应采纳2×2阶的高斯积分。有的有限元书中列出了不一样样插值函数的高斯积分点数n及相应积分点坐标i和权系数Hi的值，可供编写有限元程序时参照。五、四边形等参元的荷载等效变换四边形等参元等效结点荷载的计算，仍旧利用局部坐标系统。1．集中力PPx设在单元上随意一点C作用有一集中力,依据荷载等效变换的一般公式，其等Py效结点荷载计算公式是：RNCTP式中，NTC是形函数[N]在集中力作用点C上的值。2．体积力Wx设在单元上作用有单位体积力WWy，其等效结点荷载的计算公式是：8《有限元》讲义11{R}[N]T{W}tJdξdη113．散布力qx设单元的某一界限上承受的散布表面力是qqy，其等效结点荷载的计算公式是，当散布力作用在单元ij,mp界限上时：{R}[N]T{q}t(x)2(y)2d11当散布力作用在单元ij,mp界限上时：{R}[N]T{q}t(x)2(y)2d119《有限元》讲义2.7八结点曲线四边形等参元在常例有限元程序的单元库内，当前工程上对平面问题最合用的是八结点等参元。下边作一简要介绍。单元及结点编号如图a所示，位移和力列向量仍采纳与前面近似的摆列方式。为了结构其位移函数，依据前面等参元的见解，可先结构八结点正方形母体单元的位移函数（图b），并取为双二次插值位移函数：U222212345678V2222910111213141516代入界限条件，获取用结点位移{d}表示的位移场:8Ui1NiUi2718NiVi1式中的形函数:Ni11i1iii1i1,2,3,4427212222Nj11ij5,6,7,82jji依据等参元的定义，实质单元照耀成正方的母体单元的坐标变换式与式2-7-2等同,得:8xNixii12738yNiyii110《有限元》讲义前面我们在推导单元刚度矩阵时，都是假设单元结点编号从左下角开始，逆时针编号。从有限元理论上说，其编号应当是可以随意的。下边给出一种不一样样结点编号的例子。实质单元基本单元11《有限元》讲义12《有限元》讲义13《有限元》讲义14《有限元》讲义15《有限元》讲义16《有限元》讲义17《有限元》讲义18《有限元》讲义19《有限元》讲义20《有限元》讲义2.8几个问题的增补一、变厚问题各单元取不一样样t。二、不一样样资料问题各单元取不一样样E,μ三、平面应力与平面应变问题在前面三角形单元的推导中，我们假设其为平面应力问题。工程中，还有另一类情况─平面应变状态。比方，在对坝体或遂洞等长柱体进行分析时，假如取xoy坐标平面与其横截平面行,而Z轴与其长度方向一致（如图）那么，因为所观察物体在Z方向的尺寸很大，且又遇到平行于xoy平面，且不沿长度方向变化的荷载作用，即可以为各个横截面应处于相同的情况，即近似以为Z方向的位移重量W=0,(位移与Z没关)于是，由弹性力学知，在六个应力重量中也仅有三个独立重量σx,σy和τxy,而zxy不独立。并可获取平面应变问题的物理方程。1221xxy，yyxE1E1xy2xy1E112比较平面应力问题：11xy2xy(1)y，yyx，ExExEE得悉只要把应力问题中的E换成1，换成1即得应变问题。所以在这种问题的程序设计中，平常可以同时求解应力和应变问题，只要设置一开关变量便可以实现。四、各向异性资料在弹性矩阵[D]中反响。如正交各向异性时的弹性矩阵[D]为：21Ex

xEy

《有限元》讲义01xy

1xyD

xEy

Ey01xy0

1xyGxy（见凯维奇《有限元法》，中P102～104，英P98～109）五、设置不一样样种类的单元在工程中，同一构件常常要用到一种以上资料。如利用角、槽钢等在开洞板内作增强筋用。另一种情况是钢筋混凝土构件的全过程分析中，纵向受拉筋的单元区分，如图。混凝土被划成若干平面三角形单元，而将纵向受力钢筋（或箍筋）看作线单元（图中红线所示），也可把钢筋等效成与混凝土叠合的三角形单元，但此时均为两种不一样样资料。有限元分析时，可以为结构是由若干（面）单元（三角形或矩形）和线（杆）单元（二力杆）共同构成，区分网格时，若遇到线单元