如三节点三角形单元课件

第7章平面问题高阶单元

7.1位移模式阶次的选择

在前面两章中讨论了平面问题三结点三角形单元，其位移模式的最高阶是坐标x、y的一次项。这种位移模式导致单元常应变、常应力特性，单元应变矩阵、应力矩阵、刚度矩阵均为常数矩阵，因此计算非常简单。但这种单元难以反映应力梯度的迅速变化。要想提高计算精度，必须细分网格，增加单元数和点数，因而加大输入数据的工作量。

提高计算精度的另一条有效途径是采用高阶单元。由于高阶单元的应变、应力不再是常数，因此采用少量单元就第7章平面问题高阶单元7.1位移模式阶次的选择1可能达到较高的精度。图7-1悬臂梁分别采用高、低阶单元计算就是一个典型的例子。h4hPAB悬臂深梁解析解：A=1.0B=1.0常应变单元：A=0.866B=0.619高阶单元：A=0.99B=0.99••••••图7-1可能达到较高的精度。图7-1悬臂梁分别采用高、低阶单元计算2选择位移模式时，第2章提到要考虑解的收敛性，即要考虑到位移模式的完备性和协调性。实际操作中，一般应考虑位移模式的对称性。这是因为，有限元位移模式的选择实际是以帕斯卡（Pascal）三解形基础上的（如图7-2所示），由低价至高阶，顺序选取，组成多项式。多项式中的项数等于单元节点自由度数。如三节点三角形单元，位移模式取完全一次式，共3项。六节点三角形单元，位移模式取完全二次式共6项。如果某一阶次不能全取，则应按对称性原则适当选取。

1xyx2xyy2x3x2yxy2y3

x4x3yx2y2xy3y4

图7-2多项式选择的怕斯卡三角形选择位移模式时，第2章提到要考虑解的收敛性，即要考虑到位移模3

1xyx2xyy2x3x2yxy2y3

x4x3yx2y2xy3y4

图7-2多项式选择的怕斯卡三角形例如在下节将要讨论的四结点矩形单元中，位移模式不能取1，x，y，x2四项，也不能取1，x，y，y2四项，而应取1，x，y，xy四项。7.2四节点矩形单元图7-3示出的矩形单元，边长分别为2a和2b。取4个角点为节点，编号为i，j，l，m。将x轴和y轴置于单元的对称轴上。单元的位移函数可取为：1、位移函数1图7-24在上式表示的位移模式中，a1,a2,a3,a5,a6,a7，a8反映了单元的刚体位移和常应变。在单元的边界(x=±a或y=±a)上（或），位移是按线性分布的。因此，相邻单元在公共边上的位移是连续的。这样，位移模式满足了解答收敛性的充分条件。ijlmxyaabb图7-3在式（7-1）中代入节点位移和节点坐标后，可解出（7-1）在上式表示的位移模式中，a1,a2,5式中形函数为：（7-3）（7-2）各待定系数(a1…a8)。将这些系数再代入式（7-1），可得：式中形函数为：（7-3）（7-2）各待定系数(a1…a86则式（7-3）可简写为（7-4）将位移函数写成矩阵形式，即有与式（2-20）相同的形式（7-5）式中（7-6）令在节点上的值为：则式（7-3）可简写为（7-4）将位移函数写成矩阵形式，即有7（7-7）其中，I为二阶单位矩阵。2、应变矩阵根据几何方程，可得与式（2-25）同样的形式（7-8）把应变矩阵[B]写成子矩阵形式（7-9）其中（7-10）（7-7）其中，I为二阶单位矩阵。2、应变矩阵根据几何方程，8由此可见，[B]是、的函数，即是x、y的函数。因此单元中的应变不再是常数。3、应力矩阵根据应力-应变关系，可以计算单元中的应力，得到式（2-28）同样形式（7-11）应力矩阵[S]具有与式（2-29）同样形式（7-12）将[S]写成子矩阵形式（7-13）由此可见，[B]是、的函数，即是x、y9其中（7-14）上式对应平面应力情形。对于平面应变情形，只需将其中的E，作相应的改变即可。4、单元刚度矩阵单元刚度矩阵可采用式（2-33a）进行计算（2-33a）其中（7-14）上式对应平面应力情形。对于平面应变情形，10在四节点矩形单元中，[k]是一个8×8的矩阵。将[k]写成分块形式：（7-16）其中的子矩阵[krs]2×2可由下式计算在四节点矩形单元中，[k]是一个8×8的矩阵。将[k]写成分11（7-17）上式对应平面应力情形。对于平面应变情形，只须将上式中的E、作相应的改变。5、等价节点力单元体积力和表面力引起的节点力仍可用式（2-45）和（2-46）进行计算。（7-17）上式对应平面应力情形。对于平面应变情形，只须将12对本问题给定的位移函数，若体积力是重力的情形（设重度为），单元等价节点载荷列阵为：（2-45）（2-46）（7-18）有了对单元的上述结果，便可应用第5章的方法组集结构刚度矩阵和节点荷载向量；求解节点位移；计算内力和应力。

四节点矩形单元采用较高阶的位移模式，具有比三节点三角形单元较高的计算精度。但矩形单元也有缺点，对本问题给定的位移函数，若体积力是重力的情形（设13在三角形单元i,j,m的各边中点增设一个节点，使每个单元具有6个节点，得到图7-4所示的六节点三角形单元。这种单元具有12个自由度，可以采用完全二次多项式的位移模式：一是不能适应斜线及曲线边界，二是不便于采用大小不同的单元。7.3六节点三角形单元1、位移模式

•••ijmijmxy图7-4在三角形单元i,j,m的各边中点增设一个14（7-20）所取位移模式反映了单元的刚体位移和常应变；单元内部是连续的；在单元边界上位移分量按抛物线变化，而每条公共边界上有3个公共结点，可以保证相邻两单元位移的连续性。因此，上述位移模式满足收敛的必要和充分条件。上述位移模式确定之后，可以用分析三节点三角形单元和四节点矩形单元相同的方法进行分析。得到形函数、应变矩阵、应力矩阵、单元刚度矩阵、等价节点力向量。但其过程十分繁复，采用面积坐标可以大大简化计算。（7-20）所取位移模式反映了单元的刚体位移和常应变；单元内152、面积坐标

对于一个三角形ijm（图7-5），三角形内任一点P（x,y）的位置，可以用如下的三个比值来确定：ijmxy图7-5·P（7-21）AiAjAm（1）定义2、面积坐标对于一个三角形ijm（图7-5），三角形内任一16其中A为三角形ijm的面积，Ai,Aj,Am分别为三角形的Pjm,Pmi,Pijd的面积。这三个比值Li,Lj,Lm称为P点的面积坐标。由于则（7-22）由此可见，P点的三个面积坐标不是独立的。同时，面积坐标只是用以确定三角形内部某点的位置，因而是一种局部坐标。下面进一步给出面积坐标的几个性质。（2）面积坐标与直角坐标的关系在图7-5中，三角形Pjm的面积为其中A为三角形ijm的面积，Ai,Aj,Am分别为三角形17（7-23）由式（7-23），式（7-21）化为（7-24）将式(7-24)、(7-23a)和式(2-18)、(2-17)对比，可知，面积坐标就是三节点三角形单元的形函数Ni、Nj、Nm。（7-23a）（7-23）由式（7-23），式（7-21）化为（7-2418将式（7-24）的3个式子分别乘以xi,xj,xm，然后相加，并利用关系式（7-23a），有同理（7-25）（3）面积坐标的导数公式根据面积坐标与直角坐标的关系，由复合函数的求导公式，有将式（7-24）的3个式子分别乘以xi,xj,19（7-26）（4）面积坐标的积分公式下面给出面积坐标的幂函数积分公式。它们在计算单元刚度矩阵和等效结点载荷时有用。（7-26）（4）面积坐标的积分公式下面给出面积坐标的幂函数20在三角形单元上进行积分时，有（7-27）在三角形某一边（设ij边，边长为l）上进行积分时，有（7-28）3、用面积坐标表示六节点三角形单元计算公式

对应如图7-4所示的六节点三角形单元，形函数可用面积坐标表示为

（1）形函数和位移表达式在三角形单元上进行积分时，有（7-27）在三角形某一边（设i21•••ijmijmxy图7-6现利用形函数的性质检验式（7-29）的正确性。先考虑三角形的角点，例如图7-6中的i点，有由式（7-21）（P16），有代入式（7-29），有（7-29）•••ijmijmxy图7-6现利用形函数的性质检验式22再考虑三角形的边中点，例如i点，面积划分如图7-7所示。显然有：•••ijmijmxy图7-7由式（7-21）(P16)，有代入式（7-29）(P16)，进一步说明式（7-29）所表示的形函数的正确性。说明形函数Ni在i点等于1，在其它节点等于0，因此是正确的。再考虑三角形的边中点，例如i点，面积划分如23形函数确定后，单元中任意一点的位移可以表示为：（7-30）其中（7-31）（7-32）其中I为二阶单位阵，形函数由式（7-29）确定。（2）应变矩阵单元中的应变仍可表示为：（7-33）形函数确定后，单元中任意一点的位移可以表示为：（7-30）其24式中应变矩阵[B]为：（7-34）其中（7-35）式中应变矩阵[B]为：（7-34）其中（7-35）25单元中的应力仍可表示为：（3）应力矩阵（7-36）式中[D]是弹性矩阵，由式（2-9）确定；应变矩阵由式（7-34）、（7-35）确定。根据矩阵乘法，可以给出用面积坐标表示的应力矩阵[S]（4）单元刚度矩阵单元刚度矩阵仍可表示为：（7-37）根据[B]、[D]的表达式以及面积坐标的积分公式（7-27），可以求出[k]中元素的显式表示。由于较为繁复，这里就不列出详细结果。

单元中的应力仍可表示为：（3）应力矩阵（7-36）26（5）等价节点力向量由于位移模式是非线性的，因此体积力和表面力引起的节点力向量不能采用静力等效原理进行分配，而应采用相应公式进行计算。单元体积力引起的等价节点力计算公式仍为：（7-38）

将由式（7-29）、（7-32）表示的[N]代入，并应用积分式（7-27），可以计算FVe。例如对于重力引起的FVe，有

它表示各边中点承担单元重力的1/3。（5）等价节点力向量由于位移模式是非线性的，27单元表面力引起的结点力计算公式仍为：

（7-39）设在ij边上受有x方向的均匀分布力ps，对应的等价节点力向量为（图7-8）pslh/6pslh/64pslh/6•••ijmijmxy图7-9lps•••ijmijmxy图7-8lpspslh/62pslh/6单元表面力引起的结点力计算公式仍为：（7-39）设在ij边28如在ij边上受到x方向的三角形分布面力，其集度在i点为ps，在j点为0。对应的节点力向量为（图3-9）它表示边中点承担载荷的2/3，载荷集度大的角节点承担1/3。六结点三角形单元中的应变、应力不为常量，因此可以应用于应力梯度较大的地方，精度较高。显然，其计算也较复杂。

7.4四节点四边形等参数单元1、等参数单元的概念

如在ij边上受到x方向的三角形分布面力，其集度在i点为ps，29现在，我们从任意四边形单元着手，介绍等参数单元的概念。1234xy图7-10任意四边形单元前面讲到的四节点矩形单元虽然比较简单，但难以应用于斜线边界。图7-10所示四节点任意四边形单元容易适应这种边界，但要在整体坐标系内，写出它的统一的形函数又是相当复杂和困难的。现在，我们从任意四边形单元着手，介绍等参数单30但是若能找到它与一个规则正方形的关系，就能写出它的统一的位移模式，这可以通过坐标变换来解决。在图7-10所示四边形单元上，用等分四边的两族直线分割该四边形，以两族曲线的中心（=0、=0）为原点，沿、增大的方向作轴和轴，并令四边的=±1、=±1，就得出一组新坐标系（图7-11）。1234xy图11实际单元=-1=1=1=-1这里，、是一种局部（单元）坐标，它只应用于单元范围内。而x，y是整体（结构）坐标，它适用于所有的单元。图中的任意四边形单元是研究对象，称为实际单元。但是若能找到它与一个规则正方形的关系，就能写出它的31参照式（7-2）和（7-3）P6，此基本单元位移函数可写为：（7-40）

1234=-1=1=-1=1图7-12基本单元为了得出实际单元的位移模式和局部坐标与整体坐标之间的变换关系，引入一个四节点的正方形单元，称基本单元（图7-12）。参照式（7-2）和（7-3）P6，此基本单元位移函数可写为：32其中，形函数应为：引入新变量

i、i（i=1,2,3,4)基本单元的形函数被写成：其中，形函数应为：引入新变量基本单元的形函数被写成：33（7-41）

现在，把基本单元的位移模式（7-40）和形函数式（7-41）移用于图（7-11）所示的实际单元，则实际单元的位移模式取为：

（7-40）在结点处：在其它结点处：（7-41）现在，把基本单元的位移模式（7-34且，式中的形函数Ni仍由式（7-41）确定。而把式（7-41）中的、理解为图7-11所示实际单元的局部坐标，i、i便是实际单元中节点i的局部坐标。（7-42）利用形函数的上述性质，可以将任意四边形的整体坐标写成：任意四边形单元中结点的整体坐标，如果它已知，那么(7-42)表示了局部坐标与整体坐标的变换且，式中的形函数Ni仍由式（7-41）确定。而把式（7-4135另一方面，式（7-42）表明了实际单元中局部坐标（、）与整体坐标（x、y）的一一对应关系，是一个坐标变换式。实际单元是任意四边形四节点单元，基本单元是正方形单元，可以认为：实际单元是对基本单元通过变换得来的。由于实际单元的位移模式中采用了基本单元等同的形函数，这个实际单元就称为等参数单元。类似于本章3.2节进行的四结点矩形单元的特性分析，可以建立等参单元的应变矩阵、应力矩阵、刚度矩阵、节点力向量等的计算公式。与前面不同之处在于，在等参数单元法中，要将对整体坐标x、y的导数计算和积分计算转换为对局部坐标、的导数计算和积分计算。另一方面，式（7-42）表明了实际单元中局部36例：实际单元的结点整体坐标如图(a)中括号内数字所示，基本单元的结点局部坐标如图(b)中括号内数字所示。图(a)实际单元1(0,0)(1)试验证基本单元上的结点局部坐标与实际单元上对应点的整体坐标的对应关系。(2)求基本单元的局部坐标原点(),在实际单元上的整体坐标(x,y)是多少？4(0,1)3(1,2)2(2,0)xy0(3/4,3/4)01(－1,－1)2(1,－1)3(1,1)4(－1,1)图(b)基本单元例：实际单元的结点整体坐标如图(a)中括号内数字所示，基本单37解（1）以3结点为例，根据形函数的性质：在3结点处：应用式（7-42）说明了由基本单元上结点的局部坐标可映射出实际单元上对应的结点整体坐标（2）由于，由（7-40）得：解（1）以3结点为例，根据形函数的性质：在3结点处：应用式38代入（7-42）由上例可知：利用（7-42）在基本单元上任意一点，都可以在实际单元上找到一个对应点的坐标(x,y),这样就把实际单元与基本单元紧密地联系起来。反之，则比较困难，这是因为形函数是一个二次函数。为了避开这个困难，一般都假定基本单元上已知点去求实际单元上的对应点。代入（7-42）由上例可知：利用（7-42）在基本单元上任意39

2、应变矩阵单元的几何方程与式（7-8）、（7-9）相同，即：（7-8）（7-9）（7-43）式中（7-44）?2、应变矩阵单元的几何方程与40这里采用记号由于形函数式（3-41）是用局部坐标、给出的，将、看作x、y的函数，根据复合函数的求导规则，有：上式可记为：（7-45）这里采用记号由于形函数式（3-41）是用局部坐标41上式右边第一个矩阵称为雅可比（Jacobi）矩阵：其逆矩阵为：式中|J|为雅可比行列式（7-48）由式（7-42），有(7-46)(7-47)上式右边第一个矩阵称为雅可比（Jacobi）矩阵：其逆矩阵42（7-49）由式（7-41），有（7-50）由式（7-45），有（7-51）（7-49）由式（7-41），有（7-50）由式（7-45）43式中分别由式（7-47）和（7-50）确定。从而由式（3-43）、（3-44）确定出应变矩阵[B]。和3、应力矩阵应力矩阵仍由下式得到（7-52）

4、单元刚度矩阵单元刚度矩阵是一个8×8的矩阵，仍为式中分别由式（7-47）和（7-50）确定。从而由式（44（7-53）由于[B]是用局部坐标系、给出的，坐标变换时有面积微元公式因此，[k]可由下式计算5、等价节点力向量（1）体积力设单元的体积力是，则等价节点力公式为（7-53）由于[B]是用局部坐标系、给出的，坐标变换时45（7-54）(2)表面力设单元的某边（如对应的=±1）上作用有表面力ps=[psxpsy]T，则该边上节点的节点力为（7-55）（7-54）(2)表面力设单元的某边（如对应的46

上式中，用到坐标变换时线积分微元公式(当=常量时)对于=±1的边界表面力，等价节点力的计算过程完全一样。6、关于高斯（Gauss）积分在计算单元刚度矩阵式（7-53）和等价节点载荷向量式（7-54）、（7-55）时，由于被积函数比较复杂，通常可采用数值积分。即在单元上选择某些点，称为上式中，用到坐标变换时线积分微元公式(当=常量时)对于47积分点，求出被积函数在这些积分点上的数值，再用一些权函数乘这些函数值后求和，就可得到近似积分值。高斯积分法是数值积分法中具有较高精度的方法。现简要介绍如下：（1）一维高斯积分公式

（7-56）式中Hi对应于积分点i的权函数，对于n=2~4个积分点，其坐标i和权函数Hi列于表7-1中。积分点，求出被积函数在这些积分点上的数值，再用一些权函数乘这48n±iHi20.57735026921.000000000030.77459666920.00000000000.55555555560.888888888940.86113631160.33998104360.34785484510.6521451549表3-1高斯求积公式中的积分点坐标和权系数（2）二维高斯积分公式

（7-57）n±iHi20.57735026921.00000000049积分点j、i和权函数Hi、Hj同样可按表3-1进行取值。对式（7-53）进行高斯积分，有

（7-58）对式（7-54）进行高斯积分，有（7-59）对式（7-55）进行高斯积分，有（7-60）积分点j、i和权函数Hi、Hj同样可按表3-1进行507.5八节点四边形等参数单元上节讨论的四结点四边形等参数单元有时仍然不够理想。一是其实际单元为直线边界，不能准确拟合物体的曲线边界；二是位移模式的阶次还不够高，影响计算精度。为此，本节介绍一种精度更高、应用广泛的八节点四边形等参数单元。其实际单元和基本单元如图7-13、7-14所示。图7-14基本单元=-1=1=-1=113••567824••图7-13实际单元xy12345678••••••••=-1=1=-1=17.5八节点四边形等参数单元上节讨论的四结点51基本单元的位移模式可取为：（7-61）采用形函数表示，将位移模式写成（7-62）式中基本单元的位移模式可取为：（7-61）采用形函数表示，将位52（7-63）仿照位移模式，将坐标变换式取为（7-64）（7-63）仿照位移模式，将坐标变换式取为（7-64）53显然，该坐标变换式将平面上的正方形映射为xy平面上的曲边四边形。xy平面上每一条边都是一条二次曲线，它完全由对应边上3个结点的坐标唯一确定。因此，单元是协调的，同时也可证明，单元的位移函数反映刚体位移和常应变，具有完备性。有关八节点四边形单元的特性分析和等价节点力计算过程与四节点四边形单元完全相同，具体公式形式也一致。区别仅在于两种单元有关矩阵的维数不同，具体是：（1）四节点单元的e是一个8×1列阵；而八结点单元的e是一个16×1列阵。显然，该坐标变换式将平面上的正方形映射为54（2）四节点单元的[B]是一个3×8矩阵，由B1～B44个子块阵组成；而八节点单元的[B]是一个3×16矩阵，由B1～B88个个子块阵组成。（3）四节点单元的[k]是一个8×8矩阵，而八节点单元的[k]是一个16×16矩阵。（4）四结点单元某边表面力引起的等价节点力Fsi公式对应于该边上的2个节点；而八节点单元对应于该边上的3个节点。（2）四节点单元的[B]是一个3×8矩阵，由B1～B44个55第7章平面问题高阶单元

7.1位移模式阶次的选择

在前面两章中讨论了平面问题三结点三角形单元，其位移模式的最高阶是坐标x、y的一次项。这种位移模式导致单元常应变、常应力特性，单元应变矩阵、应力矩阵、刚度矩阵均为常数矩阵，因此计算非常简单。但这种单元难以反映应力梯度的迅速变化。要想提高计算精度，必须细分网格，增加单元数和点数，因而加大输入数据的工作量。

提高计算精度的另一条有效途径是采用高阶单元。由于高阶单元的应变、应力不再是常数，因此采用少量单元就第7章平面问题高阶单元7.1位移模式阶次的选择56可能达到较高的精度。图7-1悬臂梁分别采用高、低阶单元计算就是一个典型的例子。h4hPAB悬臂深梁解析解：A=1.0B=1.0常应变单元：A=0.866B=0.619高阶单元：A=0.99B=0.99••••••图7-1可能达到较高的精度。图7-1悬臂梁分别采用高、低阶单元计算57选择位移模式时，第2章提到要考虑解的收敛性，即要考虑到位移模式的完备性和协调性。实际操作中，一般应考虑位移模式的对称性。这是因为，有限元位移模式的选择实际是以帕斯卡（Pascal）三解形基础上的（如图7-2所示），由低价至高阶，顺序选取，组成多项式。多项式中的项数等于单元节点自由度数。如三节点三角形单元，位移模式取完全一次式，共3项。六节点三角形单元，位移模式取完全二次式共6项。如果某一阶次不能全取，则应按对称性原则适当选取。

1xyx2xyy2x3x2yxy2y3

x4x3yx2y2xy3y4

图7-2多项式选择的怕斯卡三角形选择位移模式时，第2章提到要考虑解的收敛性，即要考虑到位移模58

1xyx2xyy2x3x2yxy2y3

x4x3yx2y2xy3y4

图7-2多项式选择的怕斯卡三角形例如在下节将要讨论的四结点矩形单元中，位移模式不能取1，x，y，x2四项，也不能取1，x，y，y2四项，而应取1，x，y，xy四项。7.2四节点矩形单元图7-3示出的矩形单元，边长分别为2a和2b。取4个角点为节点，编号为i，j，l，m。将x轴和y轴置于单元的对称轴上。单元的位移函数可取为：1、位移函数1图7-259在上式表示的位移模式中，a1,a2,a3,a5,a6,a7，a8反映了单元的刚体位移和常应变。在单元的边界(x=±a或y=±a)上（或），位移是按线性分布的。因此，相邻单元在公共边上的位移是连续的。这样，位移模式满足了解答收敛性的充分条件。ijlmxyaabb图7-3在式（7-1）中代入节点位移和节点坐标后，可解出（7-1）在上式表示的位移模式中，a1,a2,60式中形函数为：（7-3）（7-2）各待定系数(a1…a8)。将这些系数再代入式（7-1），可得：式中形函数为：（7-3）（7-2）各待定系数(a1…a861则式（7-3）可简写为（7-4）将位移函数写成矩阵形式，即有与式（2-20）相同的形式（7-5）式中（7-6）令在节点上的值为：则式（7-3）可简写为（7-4）将位移函数写成矩阵形式，即有62（7-7）其中，I为二阶单位矩阵。2、应变矩阵根据几何方程，可得与式（2-25）同样的形式（7-8）把应变矩阵[B]写成子矩阵形式（7-9）其中（7-10）（7-7）其中，I为二阶单位矩阵。2、应变矩阵根据几何方程，63由此可见，[B]是、的函数，即是x、y的函数。因此单元中的应变不再是常数。3、应力矩阵根据应力-应变关系，可以计算单元中的应力，得到式（2-28）同样形式（7-11）应力矩阵[S]具有与式（2-29）同样形式（7-12）将[S]写成子矩阵形式（7-13）由此可见，[B]是、的函数，即是x、y64其中（7-14）上式对应平面应力情形。对于平面应变情形，只需将其中的E，作相应的改变即可。4、单元刚度矩阵单元刚度矩阵可采用式（2-33a）进行计算（2-33a）其中（7-14）上式对应平面应力情形。对于平面应变情形，65在四节点矩形单元中，[k]是一个8×8的矩阵。将[k]写成分块形式：（7-16）其中的子矩阵[krs]2×2可由下式计算在四节点矩形单元中，[k]是一个8×8的矩阵。将[k]写成分66（7-17）上式对应平面应力情形。对于平面应变情形，只须将上式中的E、作相应的改变。5、等价节点力单元体积力和表面力引起的节点力仍可用式（2-45）和（2-46）进行计算。（7-17）上式对应平面应力情形。对于平面应变情形，只须将67对本问题给定的位移函数，若体积力是重力的情形（设重度为），单元等价节点载荷列阵为：（2-45）（2-46）（7-18）有了对单元的上述结果，便可应用第5章的方法组集结构刚度矩阵和节点荷载向量；求解节点位移；计算内力和应力。

四节点矩形单元采用较高阶的位移模式，具有比三节点三角形单元较高的计算精度。但矩形单元也有缺点，对本问题给定的位移函数，若体积力是重力的情形（设68在三角形单元i,j,m的各边中点增设一个节点，使每个单元具有6个节点，得到图7-4所示的六节点三角形单元。这种单元具有12个自由度，可以采用完全二次多项式的位移模式：一是不能适应斜线及曲线边界，二是不便于采用大小不同的单元。7.3六节点三角形单元1、位移模式

•••ijmijmxy图7-4在三角形单元i,j,m的各边中点增设一个69（7-20）所取位移模式反映了单元的刚体位移和常应变；单元内部是连续的；在单元边界上位移分量按抛物线变化，而每条公共边界上有3个公共结点，可以保证相邻两单元位移的连续性。因此，上述位移模式满足收敛的必要和充分条件。上述位移模式确定之后，可以用分析三节点三角形单元和四节点矩形单元相同的方法进行分析。得到形函数、应变矩阵、应力矩阵、单元刚度矩阵、等价节点力向量。但其过程十分繁复，采用面积坐标可以大大简化计算。（7-20）所取位移模式反映了单元的刚体位移和常应变；单元内702、面积坐标

对于一个三角形ijm（图7-5），三角形内任一点P（x,y）的位置，可以用如下的三个比值来确定：ijmxy图7-5·P（7-21）AiAjAm（1）定义2、面积坐标对于一个三角形ijm（图7-5），三角形内任一71其中A为三角形ijm的面积，Ai,Aj,Am分别为三角形的Pjm,Pmi,Pijd的面积。这三个比值Li,Lj,Lm称为P点的面积坐标。由于则（7-22）由此可见，P点的三个面积坐标不是独立的。同时，面积坐标只是用以确定三角形内部某点的位置，因而是一种局部坐标。下面进一步给出面积坐标的几个性质。（2）面积坐标与直角坐标的关系在图7-5中，三角形Pjm的面积为其中A为三角形ijm的面积，Ai,Aj,Am分别为三角形72（7-23）由式（7-23），式（7-21）化为（7-24）将式(7-24)、(7-23a)和式(2-18)、(2-17)对比，可知，面积坐标就是三节点三角形单元的形函数Ni、Nj、Nm。（7-23a）（7-23）由式（7-23），式（7-21）化为（7-2473将式（7-24）的3个式子分别乘以xi,xj,xm，然后相加，并利用关系式（7-23a），有同理（7-25）（3）面积坐标的导数公式根据面积坐标与直角坐标的关系，由复合函数的求导公式，有将式（7-24）的3个式子分别乘以xi,xj,74（7-26）（4）面积坐标的积分公式下面给出面积坐标的幂函数积分公式。它们在计算单元刚度矩阵和等效结点载荷时有用。（7-26）（4）面积坐标的积分公式下面给出面积坐标的幂函数75在三角形单元上进行积分时，有（7-27）在三角形某一边（设ij边，边长为l）上进行积分时，有（7-28）3、用面积坐标表示六节点三角形单元计算公式

对应如图7-4所示的六节点三角形单元，形函数可用面积坐标表示为

（1）形函数和位移表达式在三角形单元上进行积分时，有（7-27）在三角形某一边（设i76•••ijmijmxy图7-6现利用形函数的性质检验式（7-29）的正确性。先考虑三角形的角点，例如图7-6中的i点，有由式（7-21）（P16），有代入式（7-29），有（7-29）•••ijmijmxy图7-6现利用形函数的性质检验式77再考虑三角形的边中点，例如i点，面积划分如图7-7所示。显然有：•••ijmijmxy图7-7由式（7-21）(P16)，有代入式（7-29）(P16)，进一步说明式（7-29）所表示的形函数的正确性。说明形函数Ni在i点等于1，在其它节点等于0，因此是正确的。再考虑三角形的边中点，例如i点，面积划分如78形函数确定后，单元中任意一点的位移可以表示为：（7-30）其中（7-31）（7-32）其中I为二阶单位阵，形函数由式（7-29）确定。（2）应变矩阵单元中的应变仍可表示为：（7-33）形函数确定后，单元中任意一点的位移可以表示为：（7-30）其79式中应变矩阵[B]为：（7-34）其中（7-35）式中应变矩阵[B]为：（7-34）其中（7-35）80单元中的应力仍可表示为：（3）应力矩阵（7-36）式中[D]是弹性矩阵，由式（2-9）确定；应变矩阵由式（7-34）、（7-35）确定。根据矩阵乘法，可以给出用面积坐标表示的应力矩阵[S]（4）单元刚度矩阵单元刚度矩阵仍可表示为：（7-37）根据[B]、[D]的表达式以及面积坐标的积分公式（7-27），可以求出[k]中元素的显式表示。由于较为繁复，这里就不列出详细结果。

单元中的应力仍可表示为：（3）应力矩阵（7-36）81（5）等价节点力向量由于位移模式是非线性的，因此体积力和表面力引起的节点力向量不能采用静力等效原理进行分配，而应采用相应公式进行计算。单元体积力引起的等价节点力计算公式仍为：（7-38）

将由式（7-29）、（7-32）表示的[N]代入，并应用积分式（7-27），可以计算FVe。例如对于重力引起的FVe，有

它表示各边中点承担单元重力的1/3。（5）等价节点力向量由于位移模式是非线性的，82单元表面力引起的结点力计算公式仍为：

（7-39）设在ij边上受有x方向的均匀分布力ps，对应的等价节点力向量为（图7-8）pslh/6pslh/64pslh/6•••ijmijmxy图7-9lps•••ijmijmxy图7-8lpspslh/62pslh/6单元表面力引起的结点力计算公式仍为：（7-39）设在ij边83如在ij边上受到x方向的三角形分布面力，其集度在i点为ps，在j点为0。对应的节点力向量为（图3-9）它表示边中点承担载荷的2/3，载荷集度大的角节点承担1/3。六结点三角形单元中的应变、应力不为常量，因此可以应用于应力梯度较大的地方，精度较高。显然，其计算也较复杂。

7.4四节点四边形等参数单元1、等参数单元的概念

如在ij边上受到x方向的三角形分布面力，其集度在i点为ps，84现在，我们从任意四边形单元着手，介绍等参数单元的概念。1234xy图7-10任意四边形单元前面讲到的四节点矩形单元虽然比较简单，但难以应用于斜线边界。图7-10所示四节点任意四边形单元容易适应这种边界，但要在整体坐标系内，写出它的统一的形函数又是相当复杂和困难的。现在，我们从任意四边形单元着手，介绍等参数单85但是若能找到它与一个规则正方形的关系，就能写出它的统一的位移模式，这可以通过坐标变换来解决。在图7-10所示四边形单元上，用等分四边的两族直线分割该四边形，以两族曲线的中心（=0、=0）为原点，沿、增大的方向作轴和轴，并令四边的=±1、=±1，就得出一组新坐标系（图7-11）。1234xy图11实际单元=-1=1=1=-1这里，、是一种局部（单元）坐标，它只应用于单元范围内。而x，y是整体（结构）坐标，它适用于所有的单元。图中的任意四边形单元是研究对象，称为实际单元。但是若能找到它与一个规则正方形的关系，就能写出它的86参照式（7-2）和（7-3）P6，此基本单元位移函数可写为：（7-40）

1234=-1=1=-1=1图7-12基本单元为了得出实际单元的位移模式和局部坐标与整体坐标之间的变换关系，引入一个四节点的正方形单元，称基本单元（图7-12）。参照式（7-2）和（7-3）P6，此基本单元位移函数可写为：87其中，形函数应为：引入新变量

i、i（i=1,2,3,4)基本单元的形函数被写成：其中，形函数应为：引入新变量基本单元的形函数被写成：88（7-41）

现在，把基本单元的位移模式（7-40）和形函数式（7-41）移用于图（7-11）所示的实际单元，则实际单元的位移模式取为：

（7-40）在结点处：在其它结点处：（7-41）现在，把基本单元的位移模式（7-89且，式中的形函数Ni仍由式（7-41）确定。而把式（7-41）中的、理解为图7-11所示实际单元的局部坐标，i、i便是实际单元中节点i的局部坐标。（7-42）利用形函数的上述性质，可以将任意四边形的整体坐标写成：任意四边形单元中结点的整体坐标，如果它已知，那么(7-42)表示了局部坐标与整体坐标的变换且，式中的形函数Ni仍由式（7-41）确定。而把式（7-4190另一方面，式（7-42）表明了实际单元中局部坐标（、）与整体坐标（x、y）的一一对应关系，是一个坐标变换式。实际单元是任意四边形四节点单元，基本单元是正方形单元，可以认为：实际单元是对基本单元通过变换得来的。由于实际单元的位移模式中采用了基本单元等同的形函数，这个实际单元就称为等参数单元。类似于本章3.2节进行的四结点矩形单元的特性分析，可以建立等参单元的应变矩阵、应力矩阵、刚度矩阵、节点力向量等的计算公式。与前面不同之处在于，在等参数单元法中，要将对整体坐标x、y的导数计算和积分计算转换为对局部坐标、的导数计算和积分计算。另一方面，式（7-42）表明了实际单元中局部91例：实际单元的结点整体坐标如图(a)中括号内数字所示，基本单元的结点局部坐标如图(b)中括号内数字所示。图(a)实际单元1(0,0)(1)试验证基本单元上的结点局部坐标与实际单元上对应点的整体坐标的对应关系。(2)求基本单元的局部坐标原点(),在实际单元上的整体坐标(x,y)是多少？4(0,1)3(1,2)2(2,0)xy0(3/4,3/4)01(－1,－1)2(1,－1)3(1,1)4(－1,1)图(b)基本单元例：实际单元的结点整体坐标如图(a)中括号内数字所示，基本单92解（1）以3结点为例，根据形函数的性质：在3结点处：应用式（7-42）说明了由基本单元上结点的局部坐标可映射出实际单元上对应的结点整体坐标（2）由于，由（7-40）得：解（1）以3结点为例，根据形函数的性质：在3结点处：应用式93代入（7-42）由上例可知：利用（7-42）在基本单元上任意一点，都可以在实际单元上找到一个对应点的坐标(x,y),这样就把实际单元与基本单元紧密地联系起来。反之，则比较困难，这是因为形函数是一个二次函数。为了避开这个困难，一般都假定基本单元上已知点去求实际单元上的对应点。代入（7-42）由上例可知：利用（7-42）在基本单元上任意94

2、应变矩阵单元的几何方程与式（7-8）、（7-9）相同，即：（7-8）（7-9）（7-43）式中（7-44）?2、应变矩阵单元的几何方程与95这里采用记号由于形函数式（3-41）是用局部坐标、给出的，将、看作x、y的函数，根据复合函数的求导规则，有：上式可记为：（7-45）这里采用记号由于形函数式（3-41）是用局部坐标96上式右边第一个矩阵称为雅可比（Jacobi）矩阵：其逆矩阵为：式中|J|为雅可比行列式（7-48）由式（7-42），有(7-46)(7-47)上式右边第一个矩阵称为雅可比（Jacobi）矩阵：其逆矩阵97（7-49）由式（7-41），有（7-50）由式（7-45），有（7-51）（7-49）由式（7-41），有（7-50）由式（7-45）98式中分别由式（7-47）和（7-50）确定。从而由式（3-43）、（3-44）确定出应变矩阵[B]。和3、应力矩阵应力矩阵仍由下式得到（7-52）

4、单元刚度矩阵单元刚度矩阵是一个