平面问题高阶单元

平面问题高阶单元第1页，课件共54页，创作于2023年2月可能达到较高的精度。图7-1悬臂梁分别采用高、低阶单元计算就是一个典型的例子。h4hPAB悬臂深梁解析解：

A=1.0

B=1.0常应变单元：

A=0.866

B=0.619高阶单元：

A=0.99

B=0.99••••••图7-1第2页，课件共54页，创作于2023年2月选择位移模式时，第2章提到要考虑解的收敛性，即要考虑到位移模式的完备性和协调性。实际操作中，一般应考虑位移模式的对称性。这是因为，有限元位移模式的选择实际是以帕斯卡（Pascal）三解形基础上的（如图7-2所示），由低价至高阶，顺序选取，组成多项式。多项式中的项数等于单元节点自由度数。如三节点三角形单元，位移模式取完全一次式，共3项。六节点三角形单元，位移模式取完全二次式共6项。如果某一阶次不能全取，则应按对称性原则适当选取。

1xyx2xyy2x3x2yxy2y3

x4x3yx2y2xy3y4

图7-2多项式选择的怕斯卡三角形第3页，课件共54页，创作于2023年2月

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、位移函数第4页，课件共54页，创作于2023年2月在上式表示的位移模式中，a1,a2,a3,a5,a6,a7，a8反映了单元的刚体位移和常应变。在单元的边界(x=±a或y=±a)上（或），位移是按线性分布的。因此，相邻单元在公共边上的位移是连续的。这样，位移模式满足了解答收敛性的充分条件。ijlmxyaabb图7-3在式（7-1）中代入节点位移和节点坐标后，可解出（7-1）第5页，课件共54页，创作于2023年2月式中形函数为：（7-3）（7-2）各待定系数(a1…a8)。将这些系数再代入式（7-1），可得：第6页，课件共54页，创作于2023年2月则式（7-3）可简写为（7-4）将位移函数写成矩阵形式，即有与式（2-20）相同的形式（7-5）式中（7-6）令在节点上的值为：第7页，课件共54页，创作于2023年2月（7-7）其中，I为二阶单位矩阵。2、应变矩阵根据几何方程，可得与式（2-25）同样的形式（7-8）把应变矩阵[B]写成子矩阵形式（7-9）其中（7-10）第8页，课件共54页，创作于2023年2月由此可见，[B]是

、

的函数，即是x、y的函数。因此单元中的应变不再是常数。3、应力矩阵根据应力-应变关系，可以计算单元中的应力，得到式（2-28）同样形式（7-11）应力矩阵[S]具有与式（2-29）同样形式（7-12）将[S]写成子矩阵形式（7-13）第9页，课件共54页，创作于2023年2月其中（7-14）上式对应平面应力情形。对于平面应变情形，只需将其中的E，

作相应的改变即可。4、单元刚度矩阵单元刚度矩阵可采用式（2-33a）进行计算（2-33a）第10页，课件共54页，创作于2023年2月在四节点矩形单元中，[k]是一个8×8的矩阵。将[k]写成分块形式：（7-16）其中的子矩阵[krs]2×2可由下式计算第11页，课件共54页，创作于2023年2月（7-17）上式对应平面应力情形。对于平面应变情形，只须将上式中的E、

作相应的改变。5、等价节点力单元体积力和表面力引起的节点力仍可用式（2-45）和（2-46）进行计算。第12页，课件共54页，创作于2023年2月对本问题给定的位移函数，若体积力是重力的情形（设重度为

），单元等价节点载荷列阵为：（2-45）（2-46）（7-18）有了对单元的上述结果，便可应用第5章的方法组集结构刚度矩阵和节点荷载向量；求解节点位移；计算内力和应力。

四节点矩形单元采用较高阶的位移模式，具有比三节点三角形单元较高的计算精度。但矩形单元也有缺点，第13页，课件共54页，创作于2023年2月在三角形单元i,j,m的各边中点增设一个节点，使每个单元具有6个节点，得到图7-4所示的六节点三角形单元。这种单元具有12个自由度，可以采用完全二次多项式的位移模式：一是不能适应斜线及曲线边界，二是不便于采用大小不同的单元。7.3六节点三角形单元1、位移模式

•••ijmi

j

m

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

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

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

（1）形函数和位移表达式第21页，课件共54页，创作于2023年2月•••ijmi

j

m

xy图7-6现利用形函数的性质检验式（7-29）的正确性。先考虑三角形的角点，例如图7-6中的i点，有由式（7-21）（P16），有代入式（7-29），有（7-29）第22页，课件共54页，创作于2023年2月再考虑三角形的边中点，例如i

点，面积划分如图7-7所示。显然有：•••ijmi

j

m

xy图7-7由式（7-21）(P16)，有代入式（7-29）(P16)，进一步说明式（7-29）所表示的形函数的正确性。说明形函数Ni在i点等于1，在其它节点等于0，因此是正确的。第23页，课件共54页，创作于2023年2月形函数确定后，单元中任意一点的位移可以表示为：（7-30）其中（7-31）（7-32）其中I为二阶单位阵，形函数由式（7-29）确定。（2）应变矩阵单元中的应变仍可表示为：（7-33）第24页，课件共54页，创作于2023年2月式中应变矩阵[B]为：（7-34）其中（7-35）第25页，课件共54页，创作于2023年2月单元中的应力仍可表示为：（3）应力矩阵（7-36）式中[D]是弹性矩阵，由式（2-9）确定；应变矩阵由式（7-34）、（7-35）确定。根据矩阵乘法，可以给出用面积坐标表示的应力矩阵[S]（4）单元刚度矩阵单元刚度矩阵仍可表示为：（7-37）根据[B]、[D]的表达式以及面积坐标的积分公式（7-27），可以求出[k]中元素的显式表示。由于较为繁复，这里就不列出详细结果。

第26页，课件共54页，创作于2023年2月（5）等价节点力向量由于位移模式是非线性的，因此体积力和表面力引起的节点力向量不能采用静力等效原理进行分配，而应采用相应公式进行计算。单元体积力引起的等价节点力计算公式仍为：（7-38）

将由式（7-29）、（7-32）表示的[N]代入，并应用积分式（7-27），可以计算FV

e。例如对于重力引起的FV

e，有

它表示各边中点承担单元重力的1/3。第27页，课件共54页，创作于2023年2月单元表面力引起的结点力计算公式仍为：

（7-39）设在ij边上受有x方向的均匀分布力ps，对应的等价节点力向量为（图7-8）pslh/6pslh/64pslh/6•••ijmi

j

m

xy图7-9lps•••ijmi

j

m

xy图7-8lpspslh/62pslh/6第28页，课件共54页，创作于2023年2月如在ij边上受到x方向的三角形分布面力，其集度在i点为ps，在j点为0。对应的节点力向量为（图3-9）它表示边中点承担载荷的2/3，载荷集度大的角节点承担1/3。六结点三角形单元中的应变、应力不为常量，因此可以应用于应力梯度较大的地方，精度较高。显然，其计算也较复杂。

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

第29页，课件共54页，创作于2023年2月现在，我们从任意四边形单元着手，介绍等参数单元的概念。1234xy图7-10任意四边形单元前面讲到的四节点矩形单元虽然比较简单，但难以应用于斜线边界。图7-10所示四节点任意四边形单元容易适应这种边界，但要在整体坐标系内，写出它的统一的形函数又是相当复杂和困难的。第30页，课件共54页，创作于2023年2月但是若能找到它与一个规则正方形的关系，就能写出它的统一的位移模式，这可以通过坐标变换来解决。在图7-10所示四边形单元上，用等分四边的两族直线分割该四边形，以两族曲线的中心（=0、=0）为原点，沿、增大的方向作轴和轴，并令四边的=±1、=±1，就得出一组新坐标系（图7-11）。1234xy图11实际单元=-1=1=1=-1

这里，、是一种局部（单元）坐标，它只应用于单元范围内。而x，y是整体（结构）坐标，它适用于所有的单元。图中的任意四边形单元是研究对象，称为实际单元。第31页，课件共54页，创作于2023年2月参照式（7-2）和（7-3）P6，此基本单元位移函数可写为：（7-40）

1234=-1=1=-1=1图7-12基本单元为了得出实际单元的位移模式和局部坐标与整体坐标之间的变换关系，引入一个四节点的正方形单元，称基本单元（图7-12）。第32页，课件共54页，创作于2023年2月其中，形函数应为：引入新变量

i、i（i=1,2,3,4)基本单元的形函数被写成：第33页，课件共54页，创作于2023年2月（7-41）

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

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

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

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

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

（7-56）式中Hi对应于积分点

i的权函数，对于n=2~4个积分点，其坐标

i和权函数Hi列于表7-1中。第48页，课件共54页，创作于2023年2月n±

iHi20.57735026921.000000000030.77459666920.00000000000.55555555560.888888888940.86113631160.3399810436