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

《有限元》讲义12.6四结点四边形单元(Thefour-nodequadrilateralelement)前面介绍了四结点的矩形单元其位移函数:xyyxU4321xyyxV8765为双线性函数，应力，应变在单元内呈线性变化，比常应力三角形单元精度高。但它对边界要求严格。本节介绍的四结点四边形等参元，它不但具有较高的精度，而且其网格划分也不受边界的影响。对任意四边形单元（图见下面）若仍直接采用前面矩形单元的位移函数，在边界上它便不再是线性的（因边界不与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点对应。《有限元》讲义2一、四结点四边形等参单元的形函数及坐标变换由于可以将整体坐标下的四边形单元变换成局部坐标下的正方形单元，对于这种正方形单元，自然仍取形函数为：2321U8765V引入边界条件，即可得位移函数：ijmpiiUNUiijmpiVNV写成矩阵形式：eepipiedNdNNNNVUf0000式中形函数:iiiN1141,pmji,,,按照等参元的定义,我们将坐标变换式亦取为:ppmmjjiiiijmpixNxNxNxNxNxppmmjjiiiijmpiyNyNyNyNyNy162式中形函数N与位移函数中的完全一致。可以验证，利用坐标变换式(2-6-1)，可以把整体坐标系中的任意四边形单元（图a）变换成在局部坐标系中与四边形对应的边长为２的正方形。因此可以将上述位移函数和形函数用于任意四边形单元，并将形函数中的ξ,η理解为任意四边形单元的局部坐标。这样由位移函数可以得到单元各点的位移。在四条边界上分别有ξ＝±和η＝±１，故边界上的位移呈线性变化，位移的连续性可得到保证。于是，我们可以理解为：任意四边形单元是从基本的正方形单元变换过来的实际单元。因此又称正方形单元为母体单元，或基本单元。例题：为了加深理解，现考察实际单元为矩形单元的坐标变换，在2.4节中，我们定义局部坐标与整体坐标的关系是：《有限元》讲义301xxa01yyb式中(x0,y0)为局部坐标原点。由上第一式01xxa得：jiijxxxxxax21210将其重新组合：jixxx121121pmjixxxx1141114111411141对照2.4中的形函数表达式,便知:ppmmjjiixNxNxNxNx自然同理可得:ppmmjjiiyNyNyNyNy由此知,矩形单元可以看作是四结点四边形单元的特例，自然，它也是等参元。《有限元法概论》（第二版）P172中，是这样解释等参元的基本概念和推导方法的：图形变换四结点正方形(母元)图形变换四结点四边形(等参元)(ξ-η平面内)───→(x,y平面内)进行图形变换的关键是进行图形结点坐标之间的变换：正方形结点坐标坐标变换四边形结点坐标(ξi,ηi)────→(x,y)i=i,j,m,pi=i,j,m,p为了实现上述结点坐标之间的变换，可利用母元的形函数，得出(ξ,η)和(x,y)之间的坐标变换式。图形变换具有如下性质:1.母元中的坐标线对应于等参元的直线；2.四结点正方形母元对应于四个结点可以任意布置的直边四边形等参元；3.变换式(2-6-1)能保证相邻等参元的边界位移彼此协调。《有限元》讲义4二、几何矩阵[B]已知单元的应变与结点位移之间的关系是:dNNNNxyyxpipi000000262形函数矩阵[N]只是局部坐标ξ,η的显函数,为求形函数对整体坐标x,y的偏导数，必须用复合函数求导公式:yyNxxNNiiiyyNxxNNiii362或写成：yNxNJNNiiiia362式中:xyxyxJb362称为雅可比矩阵，而把它的行列式称为雅可比行列式。把式(2-6-1)代入[J]得:ppmmjjiipmjipmjiijmpiiijmpiiijmpiiijmpiiyxyxyxyxNNNNNNNNyNxNyNxNJ将形函数iiiN1141pmji,,,代入，分别对,求偏导，即可得到四结点四边形等参元的雅可比矩阵：BABAJ432141562式中常数记为:ijmpiiixAijmpiiiyBijmpiix1ijmpiiy2《有限元》讲义5ijmpiix3ijmpiiy4该雅可比矩阵的逆:AABBJJ1324141雅可比行列式:AAAJ324116134213241161BAAB可以证明，如果四结点四边形的四个内角都小于180°的话，雅可比行列式|J|大于零，其逆阵[J]-1是存在的。换句话说，为了使上述等参元能保持较好的精度，整体坐标系下所划分的任意四边形单元必须是凸四边形，即任意内角都不能大于180°。四边形也不能太歪斜，否则会影响其精度。利用雅可比的逆矩阵，即可求出整体坐标系下形函数的偏导数：),,,(111141pmjiiiiiiiiiiJNNJyNxN662求出全部偏导，即代回（2-6-2）右侧，即可得到几何矩阵[B],[B]是,的函数,即:xNyNxNyNxNyNyNyNyNxNxNxNNxyyxBppjjiipjipji00000000将(2-6-6)代入即可获得[B],[B]是,的函数。三、单元刚度矩阵获得[B]后,便可由单刚的一般表达式:dxdyBDBtKT求出四结点四边形的单元刚度矩阵。在按上述公式作积分运算时,必须把面积元dxdy变换成dd,图a上的面积元abdc的面积等于矢量ab与矢量ac的矢量积的模，《有限元》讲义6即微元acabdA沿ξ轴对应于dξ的矢量增量是:jdyidxab沿η轴对应于η的矢量增量是:jdyidxac式中ji,是坐标x,y的单位矢,注意到:0jjii1ji则有:jdyidxjdyidxdAddJjiddxyyx因此刚度矩阵的积分式:1111dtdJDBKT762在计算单元刚度矩阵[K]中元素时，由于被积函数中出现了雅可比行列式，使得它用解析法很难求其积分，故常采用高斯数值积分法.《有限元》讲义7四、数值积分1．一维数值积分badF基本思想:构造一个多项目式，使在i(i=1,2……n)上有iiF，然后用近似函数的积分bad来近似原被积函数F的积分badF。i称为积分点或取样点，积分点i的数值和位置决定了近似F的程度，亦即决定数值积分的精度。对于n个积分点,按照积分点位置的不同选择，通常采用两种不同的数值积分方法，Newton-Cotes积分和高斯积分方案。二者方法基本相同，只是前者的积分点ξi是等间距分布,而后者不是等间距分布。高斯积分的积分点位置由下述方法确定：①定义n次多项式nP21②由下列条件确定n个积分点位置baidP011,0ni由上二式可见，P有以下性质：①在积分点上0iP；②多项式P与1210...,,,n在（a,b）域内正交。由此可见n个积分点的位置ξi是在求积域（a,b）内与1210...,,,n正交的n次多项式P构成方程baidP0的解。于是,被积函数F可由2n-1次多项式来近似。niniiiniPFl1101il是n-1阶Lagrange插值函数nniiil121用bad近似badF,并考虑上面确定n个积分点位置的约定条件得:baniiiRFHdF1《有限元》讲义8式中baniidlH1Hi称为积分的权系数(加权系数,权)，可见加权系数Hi与被积函数F无关，只与积分点个数和位置有关。为便于计算积分点位置ξi和加权系数Hi,常将上式中的积分限规格化，即令:1a1b由此计算出的i,Hi对应于原积分域(a,b)的关系为:ibaba22iHab2对于多重积分，按重积分规则，计算内层时，保持外层为常数，逐层计算即可。有关数值积分更详尽的资料可参阅其教科书。2．等参元计算中数值积分阶次的选择数值积分中，如何选择积分阶次将直接影响计算精度、工作量，选择不当则有可能使计算失败。选择积分阶次的原则首先是要保证积分的精度能满足所求问题的要求。如一维问题：设：插值函数中的多项式阶数为p，微分算子中导数的最高阶为m，则有限元得到的近似能量是2(p-m)次多项式。若被积函数是2(p-m)次多项式，应选高斯积分的阶次为n=p-m+1。三结点三角形单元（线性单元）刚度矩阵中被积函数是常数，故只需一个积分点。双线性单元p-m+1=1,但插值函数将包含有二次项（ξ2,ξη,η2）中的ξη项，所以要达到精确积分应采用2×2阶的高斯积分。有的有限元书中列出了不同插值函数的高斯积分点数n及相应积分点坐标i和权系数Hi的值，可供编写有限元程序时参考。五、四边形等参元的荷载等效变换四边形等参元等效结点荷载的计算，仍然利用局部坐标体系。1．集中力设在单元上任意一点C作用有一集中力yxPPP,根据荷载等效变换的一般公式，其等效结点荷载计算公式是：PNRTC式中，TCN是形函数[N]在集中力作用点C上的值。2．体积力设在单元上作用有单位体积力yx，其等效结点荷载的计算公式是：1111{}[]{}TRN