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

1第三章轴对称、三维和高次单元§3-2空间问题的四面体单元空间问题的有限单元法，和平面问题及轴对称问题的有限单元法的原理和分析过程完全相同。由于空间问题应采用三维坐标系，因此单元的自由度、刚度矩阵的元素个数，方程组内方程个数等要较平面问题和轴对称问题多，所以空间问题的规模一般比轴对称问题和平面问题大得多。它要求计算机的内存大，且计算时间长，费用高。这些问题都给三维和平面问题一样，空间有限单元法采用单元也是多种多样的，其中最简单的是四节点四面体单元。采用四面体单元和线性位移模式来处理空间问题，可以看作平面问题中三角形单元的推广。在采用四面体单元离散化后的空间结构物中，一系列不相互重叠的四面体之间仅在节点处以空间铰相互连接。四节点四面体单元仅在四个顶点处取为节点，其编号为i,j,m,p。每个单元的计算简图如图3-7所示。在位移法中，取节点位移为基本未知量，四节点四面体单元共有十二个自由度(位移分量)，其节点位移列阵为Tpppmmmjjjiiipmjiewvuwvuwvuwvu其子矩阵iiiiwvu(i,j,m)相应的节点力列阵为TpmjieFFFFF图3-7空间四面体单元2其子矩阵TiiiiWVUF一、单元法位移函数结构中各点的位移是坐标x、y、z的函数。当单元足够小时，单元内各点的位移可用简单的线性多项式来近似描述，即zyxwzyxvzyxu121110087654321(3-49)式中1，2，…，12是十二个待定系数，它们可由单元的节点位移和坐标确定。假定节点i,j,m,p的坐标分别为(ixiyiz)、(jxjyjz)、(mxmymz)、(pxpypz)，将它们代入(3-49)式的第一式可得各个节点在x方向的位移ppppmmmmjjjjiiiizyxuzyxuzyxuzyxu4321432143214321(3-50)解上述线性方程组，可得到1，2，3，4，再代入(3-50)式，得])()()()[(61pppppmmmmmjjjjjiiiiiuzdycxbauzdycxbauzdycxbauzdycxbaVu(3-51)其中V为四面体ijmp的体积，ai,bi，…，cp,dp为系数。pppmmmjjjiiizyxzyxzyxzyxV1111(3-52)3pppmmmjjjizyxzyxzyxa111jjimmppyzbyzyzppmmjjizxzxzxc111111ppmmjjiyxyxyxd(i,j,m,p)(3-53)为了使四面体的体积V不致为负值，单元四个节点的标号i,j,m,p必须按照一定的顺序：在右手坐标系中，要使得右手螺旋在按照i→j→m的转向转动时，向p的方向前进，象图3-1中单元那样。用同样方法，可以得出其余二个位移分量：])()()()[(61pppppmmmmmjjjjjiiiiivzdycxbavzdycxbavzdycxbavzdycxbaVv(3-54)])()()()[(61pppppmmmmmjjjjjiiiiiwzdycxbawzdycxbawzdycxbawzdycxbaVw(3-55)综合表达式(3-51)、(3-54)及(3-55)，可以将位移分量表示成为epmjieTININININNwvuf][(3-56)其中I是三阶的单位矩阵，[N]为形函数矩阵，而各个形函数为),(6/)(),(6/)(pjVzdycxbaNmiVzdycxbaNiiiijiiiii(3-57)和平面问题相似，(3-49)式中的系数1，5，6代表刚性移动0u，0v，0w；系数2，7，12代表常量的正应变；其余6个系数反映了刚性转动xw，yw，zw和常量剪应变。这就是说，12个系数充分反映了单元的刚体位移和常量应变。同时，可以证明：由于位移模式是线性的，两个相邻单元的共同边界在变形过程中，始终是相互贴合的，使4得离散的模型变形中保持为连续体。这样，选用的位移函数满足收敛的充分必要条件，保证了有限单元法解答收敛于精确解。二、载荷移置空间问题的单元载荷移置和平面问题一样，也是根据静力等效原则，将不作用在节点上的集中力、体力、面力移置成作用在节点上的等效节点载荷。其通用公式的形式和平面问题也是一样的，只不过多出一维空间分量。1.集中力设单元上某点(x,y,z)作用有集中力TxyzPPPP则仍然得到等效节点载荷PNRT][(3-58)这里TpppmmmjjjiiieZYXZYXZYXZYXR][2.分布体力单元上作用有分布体力TZYXP][，则dVPNRTe][(3-59)其中dV是单元中的微分体积，对于直角坐标系上式为dxdydzpNReTe][(3-60)3.分布面力单元的某一边界面S上作用有一分布面力TZYXP则dAPNRTe][其中dA是边界面S上的微分面积。4.常见载荷的移置5上列公式是空间问题载荷移置的通用公式。对于四节点四面体单元，由于其采用线性位移模式，采用直接计算虚功的方法求出节点载荷比较简单。下面介绍常见的二种载荷的移置。(1)重力四面体单元的自重为W，作用在质心C处(如图3-8)。为求得节点载荷Xi,Yi,Zi,可分别假想发生1*iu,1*iv或1*iw的虚位移。在1*iu或1*iv时，整个单元上各点的均没有z方向上的虚位移，重力W不做功，所以Xi=Yi=0。当1*iw时，jmp面上各点的虚位移为零，即0*bw，又因bibc41,所以有41*cw，4WZi对于其余三个节点可得同样结论，于是有TeiWR400(i,j,m,p)(3-61)即，对于四节点四面体单元承受的重力载荷，只需要把共41移置到每个节点上即可。(2)界面压力设四面体的一个边界面ijm上受有一线性分布的压力P，共在三个节点上的强度分别为qi,0,0。很容易看出，该力向p点移置的等效节点力为零。由水力学知，总压力ijmqPi31,作用于ijm面上的d点，d点到ij边和im边的距离分别为m到ij及j到im边的距离的1/4。于是可得TijmiTeiqPPPR021211610442(3-62)所得各节点载荷的方向和分布力的方向相同，要求各节点载荷分量还需乘上相应的方向余弦。图3-8重力移置6由上述面力移置结果，可求出任意线性分布面的等效节点载荷。如在ijm面受有线性分布面力在各点强度分别为qi，qj，qm，时，在i节点的等效载荷为ijmmjiiqqqP)2121(61(i，j，m)(3-63)三、应力应变矩阵空间问题几何方程为Tzyxzyxzuxwywzuxvyuzwyvxu将四面体单元之位移表达式(3-52)、(3-54)和(3-55)代入几何方程，即得单元应变。用节点位移可表示为epmjieBBBBB][(3-64)式中应变矩阵子矩阵为6×3矩阵：iiiiiiiiiibdcdbcdcbVB00000000061][(i,j,m,p)(3-65)由上式可以看出，每一个单元的应变矩阵是一个常量矩阵；因此，采用线性位移模式的四面体单元是常应变单元。这与平面问题中的三角形单元是一样的。而与平面问题的不同之处仅在于应变矩阵的阶数不同。将表达式(3-16)代入空间问题的物理方程，即可得出用单元节点位移表示的单元应力：eeSBDD][][][(3-66)7式中弹性矩阵[]D为)1(22100000)1(2210000)1(221000111111][称对D应力矩阵pmjiSSSS][S(3-67)令11A，)1(2212A则iiiiiibAdAcAdAbAVES22222i2ii1i1i1ii1i1i1i000cAdcAbAdAcbAdAcAb)21)(1(6)1(][(i,j,m,p)(3-68)显然，式(3-68)中各元素均为常量，应力矩阵[Ｓ]是常量矩阵，所以，四面体单元是常应力单元。四、单元刚度矩阵空间问题的单元刚度由虚功方程导出。假设该单元发生某虚位移，相应节点虚位移为*e。此时相应的虚应变为8eB**][将上式及式(3-66)代入虚功方程，有dxdydzBDBFeTeveTe]][[)]([)(**通过与平面问题一样的处理，并注意到矩阵[B]中的元素为常量，可以得到eeeTevTeKVBDBdxdydzBDBF][]][[][]][[][(3-69)式中，eK][为单元刚度矩阵：eTTeVBDBdxdydzBDBK]][[][]][[][][(3-70)将式(3-64)和(3-68)式代入，可以得出pppmpjpimpmmmjmijpjmjjjiipimijiieKKKKKKKKKKKKKKKKK][(3-71)其中，ersK][为3×3阶方阵：212121221212122()(1)[]()36(1)(1)()rsrsrsrsrsrsrsersrsrsrsrsrsrsrsrsrsrsrsrsrsrsbbAccddAbcAccAbdAdbEKAcbAbcccAbbddAcdAdcVAdbAbdAdcAcdddAbbcc(r,s=i,j,m,p)(3-72)有了单元节点力和节点位移之间的关系之后，通过分析每个节点的平衡条件可得到erepmjissersRK,,,][这个矩阵形式的方程实际上代表了关于r节点三个坐标轴方向的力平衡方程式。将关于结构物所有节点的线性方程式集合起来，可以得到9RK][式中代表整个结构的节点的位移，是所求之基本未知量；R代表整个结构的节点载荷；][K为整体刚度矩阵，其是由每个单元刚度矩阵升阶后组集得到，即NEeeKK1][][其为3NP阶方阵。显然，对每一个子矩阵，应有NEeersrsKK1][][和平面问题一样，][K是对称、带状、稀疏矩阵，在消除刚体位移之后，它是正定的。由平衡方程组可以解出节点位移，随后即可求得所需节点和单元应力。五、形成四面体的对角线划分方法在实际计算中，用一系列的四面体来组合成一个空间物体，这个形象是很难想象的。但是如果先用一系列较为直观的六面体(图3-9)来划分弹性体，然后由计算机来将这些六面体及三棱柱划分为若干个四面体，则要方便得多。同时减少许多准备及输入工作，也为将来结果分析带来方便。现在介绍一种适合计算机进行自动划分四面体的方法——对角线划分法。1.将六面体划分为四面体的方法通过连接六面体上一些四边形的对角线，可以把一个六面体划分为五个或六个四面体。为叙述方便，先将六面体的八个角点进行局部编号，编号原则是先顶面后底面，对于顶面或底面的节点来说，则是先前后后，从左到右(见图3-9)排列。10图3-9六面体和三棱柱(1)将一个六面体划分为五个四面体这种方法是先过六面体的一些四边形的对角线，从六面体的四个角上切下四个四面体，最后剩下中心的一个四面体