03-02空间问题的四面体单元解析

1、1第三章 轴对称、三维和高次单元 3-2 空间问题的四面体单元空间问题的有限单元法，和平面问题及轴对称问题的有限单元法的原理和分析过程完 全相同。由于空间问题应采用三维坐标系，因此单元的自由度、刚度矩阵的元素个数，方 程组内方程个数等要较平面问题和轴对称问题多，所以空间问题的规模一般比轴对称问题 和平面问题大得多。它要求计算机的内存大，且计算时间长，费用高。这些问题都给三维有限单元法的具体运用带来许多困难。和平面问题一样，空间有限单元法采用单元 也是多种多样的，其中最简单的是四节点四面体 单元。采用四面体单元和线性位移模式来处理空 间问题，可以看作平面问题中三角形单元的推广。在采用四面体单元离

2、散化后的空间结构物 中，一系列不相互重叠的四面体之间仅在节点处 以空间铰相互连接。四节点四面体单元仅在四个 顶点处取为节点，其编号为i,j,m,p。每个单元的 计算简图如图3-7所示。在位移法中，取节点位移为基本未知量，四 节点四面体单元共有十二个自由度（位移分量），其节点位移列阵为WiUjVjWjUm2VmWmUpVpWpT相应的节点力列阵为F-F-F FTij m p其子矩阵【iUiViw(i,j,m)3、单元法位移函数结构中各点的位移是坐标x、y、z的函数。当单元足够小时，单元内各点的位移可用 简单的线性多项式来近似描述，即u =% +ot2x +a3y+口4乙V =0( (5 +0(

3、(6X +o7y+口8乙 w=a0+a10 x+a11y+a12z式中：1，: 2，,，：12是十二个待定系数，它们可由单元的节点位移和坐标确定。假定节点i,j，m，p的坐标分别为(XjYiZj)、(XjYjZj)、(XmYmZm)、(Xpypzp)，将它们代入(3-49)式的第一式可得各个节点在X方向的位移：1，：- 2，：3，： 4，再代入(3-50)式，得1u -(aibiX ciy diz)ui-(a,b,X c.y djZ)Uj6Vj j2j/j(3-51)其子矩阵- UiViWiT(3-49)Ui=a i+ 2Xi+a3yi+ j Ziuj=a1+2Xj+cf3yj+jZj Um=

4、a1gXmgym十JZmUP=a1+ a2Xp+a3yp+ C(4Zp(3-50)解上述线性方程组，可得到4(ambmX cmy dmz)um-(apbpX Cpy dpZ)Up其中V为四面体ijmp的体积，ai,bi，,1XiyiZi1XjyjZj1XmymZm1XPYPZP(3-52)cp,dP为系数。5=(色bx Gy diZ)/6V(i,m)(冃bx yy djZ)/6/(j, p)和平面问题相似，(3-49)式中的系数,:-5,:-6代表刚性移动U0,V0,W0;系数7,12代表常量的正应变；其余6个系数反映了刚性转动Wx,Wy,Wz和常量剪应变。这就是说，12个系数充分反映了单元的

5、刚体位移和常量应变。同时，可以证明： 由于位移模式是线性的， 两个相邻单元的共同边界在变形过程中，始终是相互贴合的，使 得离散的模型变形中保持为连续体。这样，选用的位移函数满足收敛的充分必要条件，保 证了有限单元法解答收XjyjZj1yjzjai =XmymZmb =1 YmZmXpypZp1YPZPXj1ZjXjyjCi =Xm1 Zmdi =XmymXp1 ZpXPyp为了使四面体的体积V不致为负值，单元四个节点的标号 在右手坐标系中，要使得右手螺旋在按照3-1中单元那样。i,j,m,piTjTm的转向转动时,必须按照一定的顺序：向p的方向前进，象图用同样方法，可以得出其余二个位移分量:1

6、v佝bx yy diZ)wbjX y djz)Vj6V(3-54)(ambmX Cmy dmZ)Vm-(apbpX CpydpZ)Vp1w =(ai+3%+&丫9乙网一fx+Ciy+ djZ)Wj6Vj j j j j(3-55)(ambmX Cmy dmZ)Wm-(apbpX Cpy dpZ)Wp综合表达式(3-51)、(3-54)及(3-55)，可以将位移分量表示成为f=U v w=N6e=lINiINjINmINpL (3-56)其中I是三阶的单位矩阵,N为形函数矩阵，而各个形函数为NiNj(3-57)11(i,j,m,p)(3-53)16敛于精确解。、载荷移置空间问题的单元载荷

7、移置和平面问题一样，也是根据静力等效原则，将不作用在节点 上的集中力、体力、面力移置成作用在节点上的等效节点载荷。其通用公式的形式和平面 问题也是一样的，只不过多出一维空间分量。1.集中力设单元上某点(x,y,z)作用有集中力pJ巳PyPz|T则仍然得到等效节点载荷R=NTfp?(3-58)这里RF=XiYZiXjYjZjXmYmZmXpYpZpT2.分布体力单元上作用有分布体力7pf=X Y ZT，则= .NT：PdV(3-59)其中dV是单元中的微分体积，对于直角坐标系上式为二NTpdxdydz(3-60)e3.分布面力单元的某一边界面S上作用有一分布面力lP：= X Y ZTR NTP

8、dA其中dA是边界面S上的微分面积。4.常见载荷的移置1267上列公式是空间问题载荷移置的通用公式。 对于四节点四 面体单元，由于其采用线性位移模式， 采用直接计算虚功的方 法求出节点载荷比较简单。下面介绍常见的二种载荷的移置。重力四面体单元的自重为W作用在质心C处(如图3-8)。为 求得节点载荷X， 丫i,Zi,可分别假想发生u* =1,V；=1或W；=1的虚位移。在U*=1或V* = 1时，整个单元上各点的均没有z方向上8的虚位移，重力W不做功，所以X=Yi=O。(3-62)所得各节点载荷的方向和分布力的方向相同， 要求各节点载荷分量还需乘上相应的方向余 弦。当W,1时，jmp面上各点的虚

9、位移为零，即Wb= 0，又因b-bi,所以有4对于其余三个节点可得同样结论，于是有i,j,m,p)(3-61)即， 对于四节点1移置到每个节点上即可。4(2)界面压力设四面体的一个边界面ijm上受有一线性分布的压力P,共在三个节点上的强度分别为qi,0,0。很容易看出，该力向p点移置的等效节点力为零。1由水力学知，总压力P = qijm,3作用于ijm面上的d点，d点到ij边和im边的距离分别为m到ij及j到im边的距离的1/4。于是可得图 3-8 重力移置9由上述面力移置结果，可求出任意线性分布面的等效节点载荷。如在i j m面 受 有 线 性分布面力在各点强度分别为qi,qj,qm,时，在

10、i节点的等效载荷为111P(qiqjqm)-ijm(i，j，m)(3-63)6 2 2三、应力应变矩阵空间问题几何方程为点位移可表示为=BI严二Bj-BjBm一Bp L严(3-64)式中应变矩阵子矩阵为6X3矩阵:由上式可以看出，每一个单元的应变矩阵是一个常量矩阵；因此，采用线性位移模式 的四面体单元是常应变单元。这与平面问题中的三角形单元是一样的。而与平面问题的不 同之处仅在于应变矩阵的阶数不同。将表达式(3-16)代入空间问题的物理方程，即可得出用单元节点位移表示的单元应力:.:u；v-T -:y: x.:u : w : w : u十-十dzdy dx dz将四面体单元之位移表达式和(3-

11、55)代入几何方程，即得单元应变。用节_bi5 n1Bi = 6V(3-65)(3-52)、(3-54)10匕?- D - D B：F =sFe(3-66)11S= SiSjSmSp 1(3-67)1 2)显然，式(3-68)中各元素均为常量，应力矩阵 常应力单元。式中弹性矩阵D为D二1对1 -0101 -2J称2(1)0001 -2J2(1_i)00001 22(1二)应力矩阵_biACA1biCiAQA1biACdiA2GA2bi00A2diA2G.Azdi0A2bi -S=6(1E(；S(3-68)S是常量矩阵，所以，四面体单元是12四、单元刚度矩阵空间问题的单元刚度由虚功方程导出。假设

12、该单元发生某虚位移，相应节点虚位移为。此时相应的虚应变为13将上式及式(3-66)代入虚功方程，有(虽*bTFT =HJ(BbF)TDBedxdydzv通过与平面问题一样的处理，并注意到矩阵B中的元素为常量，可以得到T ?：!iBTDBdxdydzI F二BT D B*、V二Kge(3-69)v式中，Ke为单元刚度矩阵:KBTDBdxdydz =BTDBVe将式(3-64)和(3-68)式代入，可以得出其中，Krse为3X3阶方阵:(r,s=i,j,m,p)(3-72)有了单元节点力和节点位移之间的关系之后，通过分析每个节点的平衡条件可得到 KJeJ八e s=t,j ,m,pe这个矩阵形式的方

13、程实际上代表了关于r节点三个坐标轴方向的力平衡方程式。将关于结构物所有节点的线性方程式集合起来，可以得到(3-70)KHK“Kim-Kip 1Ke=-心Kjj_KjmKjpKmi- KmjKmm- Kmp_KpiKpjKpmKpp一(3-71)KrseE(1)36(1)(1一J)VAGbs+ACsAdrbs+AAdsAbrCsAGCsCrCsA2(brbsdrds)AdrCsAGdsAbrds+A2drbs1ACrds+AzdrCsdrdsA2(brbsCrCs)14K-R式中、J代表整个结构的节点的位移， 是所求之基本未知量；代表整个结构的节点载荷;K为整体刚度矩阵，其是由每个单元刚度矩阵升

14、阶后组集得到，即NEK八Ke=1其为3NP阶方阵。显然，对每一个子矩阵，应有NEKJ八心e=1和平面问题一样，K是对称、带状、稀疏矩阵，在消除刚体位移之后，它是正定的。由平衡方程组可以解出节点位移，随后即可求得所需节点和单元应力。五、形成四面体的对角线划分方法在实际计算中，用一系列的四面体来组合成一个空间物体，这个形象是很难想象的。 但是如果先用一系列较为直观的六面体（图3-9）来划分弹性体，然后由计算机来将这些六 面体及三棱柱划分为若干个四面体，则要方便得多。同时减少许多准备及输入工作，也为 将来结果分析带来方便。现在介绍一种适合计算机进行自动划分四面体的方法一一对角线划分法。1.将六面体划

15、分为四面体的方法通过连接六面体上一些四边形的对角线，可以把一个六面体划分为五个或六个四面体。为叙述方便，先将六面体的八个角点进行局部编号，编号原则是先顶面后底面，对于顶面 或底面的节点来说，则是先前后后，从左到右（见图3-9）排列。15图 3-9六面体和三棱柱(1)将一个六面体划分为五个四面体这种方法是先过六面体的一些四边形的对角线，从六面体的四个角上切下四个四面体，最后剩下中心的一个四面体，共得五个六面体单元。选择被切下的角点不同，有二种不同 的划分结果，如图3-10(a )和(b )所示。我们分别称之为A5型划分和B5型划分。A5型划分所得五个四面体为1246,1347,1467,1567

16、,4678;而B5型划分则得到1235，2348,2358,2568,3578五个四面体。以上二种划分方法的共同特点是，六面体二对面四边形的剖分对角线是交叉的。这就使得如果一个六面体按A5型划分，那么与之相邻的各个六面体必定要按B5型划分。图 3-10六面体划分为五个四面体(a) A5 型剖分； (b)B5 型剖分(2)将一个六面体划分成六个四面体616将六面体划分成六个四面体有很多种划分方法。这里介绍两种，如图3-11所示。它们的共同特点是，六面体上两对面四边形的剖分对角线是“平行”的。所不同的是在A6型剖分中取大对角线36作为划分线，而在B6型中则是取大对角线45作为划分线。为清楚起见，可

17、将A6型划分理解为先将六面体沿2367分成两个三棱柱，再将每个三棱柱分成三个四面体，分别得到1235，2356，3567和2346，3467,4678六个四面体(见图3-12。当然，A6型划分也可看成先将六面体沿3456面剖分，得到两个不同于前的三棱柱，但最后得到的六个四面体是相同的(图3-13)。对于B6型划分，六面体先以折面2457为分界面拆分成两个三棱柱，如图3-14所示。于是可见，每一个“三棱柱”被划分为三个四面体，它们分别是1235，2345，3456和2456，4567，4678。同时，也不难证明，若以3456为分界面按B6型划分将六面体拆成的二个“三棱柱”虽与前面的不同，但是划分

18、成的六个四面体和前面得到的完全相同。(a)(b)图 3-11六面体划分为六个四面体(a) A6 型剖分；(b)B6 型剖分图 3-12 A6 型划分，以折面 2376 为两个“三棱柱”的分界面17我们看到，A6型和B6型划分，由于其相对四边形的对角线“平行”，而剖分大对角线35和45并不在六面体表面上，其相邻的六面体可以全部采用A6型划分或B6型划分，两种划分也可以交替使用。一个六面体划分为六个四面体，各四面体的体积大小一般较为均 匀；但是在相等的六面体数目下，A6型和B6型划分所产生的四面体单元的总数，要比A5型和B5型产生的多六分之一。此外，在A6型与B6型划分中，如果离散体的节点整体编号

19、是按本节开头所述，从上 到下，从左到右连续进行的；同时每个六面体八个节点的整体编号的大小次序与其局部编 号的大小次序相一致(由小到大)的话，那么在划分中，对底面上任一节点，与它构成四面 体的三个节点中的最小号码，不会比其正上方那个节点的号码更小。这是由于在A6型及B6型划分中，注意到在连各个四边形对角线时，不使节点编号之差较大的三个节点出现在 同一个三角形中的结果。例如图3-11中，对于1357四边形，我们连接了35对角线而使1、7两节点分别属于两个三角形中。这样的划分能获得一个带宽较窄的刚度矩阵。(3)编号推算如果将六面体的八个顶点的节点整个编号置于数组D1:8中，而将前述图中的局图 3-1

20、3 A6型划分，以折面 3456 为两个“三棱柱”的分界面18部编号18理解为数组D1:8的下标时，于是上述问题就转化为：要在有八个元素的数组中按一定规律， 每次取四个元素构成一个四面体单元的节点编号问题。对于A6型和A5型划分所得四面体顶点编号的规律性进行一些分析之后，可以导出下列公式，分别表示 按预定规律划分成的四面体的各节点编号：Dm(1 - 3(I 1)J(1 -m)(l J)Dm(3l J一1)(1 m)(1 I J(J 1)/2)Dm(2(1I -J) J(lJ) (1m)(2 IJ(5 J)/2)(3-73)Dm(5 IJ(3 - J)/2 (1 -m)(4 IJ)(I =1,2

21、; J =0,1,2;m =0,1)通过直接代入数字检验，知道m=0对应着A6型划分所得的六个四面体，m=1则对应着A5型划分所得的五个四面体，(此时I=J=2形成的四面体应舍去)。例如m=0的情况，当I=1，J=0，1，2时，得到D1D2D3D5，D2D3D5D6，D3D5 D6D7三个四面体，而I=2，J= 0,1, 2时，得到D2D3D4D6，D3D4D6D7，D4D6D7D8三个四面体，与前述结果一致。对于B6型及B5型划分，同样可以导出一个相似的计算公式。但是也可以利用(3-73)式，只需将原来D1：8中元素位置作一定更动。更动的办法是，对于B6型划分，将8、6、4、2位置的元素置于1、2、3、4位置上，将