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

第三章 轴对称、三维和高次单元§3-2 空间问题的四面体单元空间问题的有限单元法，和平面问题及轴对称问题的有限单元法的原理和分析过程完全相同。

由于空间问题应采用三维坐标系，因此单元的自由度、刚度矩阵的元素个数，方程组内方程个数等要较平面问题和轴对称问题多，所以空间问题的规模一般比轴对称问题和平面问题大得多。

它要求计算机的内存大，且计算时间长，费用高。

这些问题都给三维有限单元法的具体运用带来许多困难。

和平面问题一样，空间有限单元法采用单元也是多种多样的，其中最简单的是四节点四面体单元。

采用四面体单元和线性位移模式来处理空间问题，可以看作平面问题中三角形单元的推广。

在采用四面体单元离散化后的空间结构物中，一系列不相互重叠的四面体之间仅在节点处以空间铰相互连接。

四节点四面体单元仅在四个顶点处取为节点，其编号为i,j,m,p 。

每个单元的计算简图如图3-7所示。

在位移法中，取节点位移为基本未知量，四节点四面体单元共有十二个自由度(位移分量)，其节点位移列阵为{}[]Tpp p m m m j jj i i ip m j i ew v u w v u w v u w v u =⎪⎪⎭⎪⎪⎬⎫⎪⎪⎩⎪⎪⎨⎧=δδδδδ其子矩阵 {}[]i ii i w v u =δ (i,j,m)相应的节点力列阵为{}[]Tp m j ie F F F F F -图3-7 空间四面体单元其子矩阵 {}[]Ti i i i W V U F =一、单元法位移函数结构中各点的位移是坐标x 、y 、z 的函数。

当单元足够小时，单元内各点的位移可用简单的线性多项式来近似描述，即⎪⎭⎪⎬⎫+++=+++=+++=z y x w z y x v z y x u 121110087654321αααααααααααα (3-49) 式中1α，2α，…，12α是十二个待定系数，它们可由单元的节点位移和坐标确定。

假定节点i,j,m,p 的坐标分别为(i x i y i z )、(j x j y j z )、(m x m y m z )、 (p x p y p z )，将它们代入(3-49)式的第一式可得各个节点在x 方向的位移⎪⎪⎭⎪⎪⎬⎫⎪⎪⎩⎪⎪⎨⎧+++=+++=+++=+++=p p p p m m m m j j j j i i i i z y x u z y x u z y x u z y x u 4321432143214321αααααααααααααααα (3-50)解上述线性方程组，可得到1α，2α，3α，4α，再代入(3-50)式，得])()()()[(61p p p p p m m m m m jj j j j i i i i i u z d y c x b a u z d y c x b a u z d y c x b a u z d y c x b a Vu +++-+++++++-+++=(3-51) 其中V 为四面体ijmp 的体积，a i ,b i ，…，c p ,d p 为系数。

pppm m m j j j i i i z y x z y x z y x z y x V 1111=(3-52)p p pm m mj j ji z y x z y x z y x a = 111j ji m m p py z b y z y z = p p m m j j i z x z x z x c 111= 111p pm mj ji y x y x y x d = (i,j,m,p) (3-53) 为了使四面体的体积V 不致为负值，单元四个节点的标号i,j,m,p 必须按照一定的顺序：在右手坐标系中，要使得右手螺旋在按照i →j →m 的转向转动时，向p 的方向前进，象图3-1中单元那样。

用同样方法，可以得出其余二个位移分量：])()()()[(61p p p p p m m m m m jj j j j i i i i i v z d y c x b a v z d y c x b a v z d y c x b a v z d y c x b a Vv +++-+++++++-+++=(3-54) ])()()()[(61p p p p p m m m m m jj j j j i i i i i w z d y c x b a w z d y c x b a w z d y c x b a w z d y c x b a Vw +++-+++++++-+++=(3-55) 综合表达式(3-51)、(3-54)及(3-55)，可以将位移分量表示成为{}[]{}[]{}ep mj ieTIN IN IN IN N w v uf δδ===][ (3-56)其中I 是三阶的单位矩阵，[N]为形函数矩阵，而各个形函数为⎭⎬⎫+++-=+++=),(6/)(),(6/)(p j Vz d y c x b a N m i V z d y c x b a N i i i i j i i i i i (3-57) 和平面问题相似，(3-49)式中的系数1α，5α，6α代表刚性移动0u ，0v ，0w ；系数2α，7α，12α代表常量的正应变；其余6个系数反映了刚性转动x w ，y w ，z w 和常量剪应变。

这就是说，12个系数充分反映了单元的刚体位移和常量应变。

同时，可以证明：由于位移模式是线性的，两个相邻单元的共同边界在变形过程中 ，始终是相互贴合的，使得离散的模型变形中保持为连续体。

这样，选用的位移函数满足收敛的充分必要条件，保证了有限单元法解答收敛于精确解。

二、载荷移置空间问题的单元载荷移置和平面问题一样，也是根据静力等效原则，将不作用在节点上的集中力、体力、面力移置成作用在节点上的等效节点载荷。

其通用公式的形式和平面问题也是一样的，只不过多出一维空间分量。

1. 集中力设单元上某点(x,y,z)作用有集中力{}Txyz P P P P ⎡⎤=⎣⎦则仍然得到等效节点载荷{}{}P N R T ][= (3-58)这里 {}T p pp m m m j j j i i ieZ Y X Z Y X Z Y X Z Y X R ][=2. 分布体力单元上作用有分布体力{}T Z YXP ][=，则{}{}dV P N R T e ⎰=][ (3-59)其中dV 是单元中的微分体积，对于直角坐标系上式为{}{}dxdydz p N R eT e ⎰⎰⎰=][ (3-60)3. 分布面力单元的某一边界面 S 上作用有一分布面力{}[]TZ YX P =则 {}{}dA P N R T e⎰=][其中dA 是边界面S 上的微分面积。

4. 常见载荷的移置上列公式是空间问题载荷移置的通用公式。

对于四节点四面体单元，由于其采用线性位移模式，采用直接计算虚功的方法求出节点载荷比较简单。

下面介绍常见的二种载荷的移置。

(1) 重力四面体单元的自重为W ，作用在质心C 处(如图3-8)。

为求得节点载荷X i ,Y i ,Z i ,可分别假想发生1*=i u ,1*=i v 或1*=i w 的虚位移。

在1*=i u 或1*=i v 时，整个单元上各点的均没有z 方向上的虚位移，重力W 不做功，所以X i =Y i =0。

当1*=i w 时，jmp 面上各点的虚位移为零，即0*=b w ，又因bi bc 41=,所以有 41*=c w ， 4WZ i -= 对于其余三个节点可得同样结论，于是有{}Tei W R ⎥⎦⎤⎢⎣⎡-=400(i,j,m,p) (3-61)即，对于四节点四面体单元承受的重力载荷，只需要把共41移置到每个节点上即可。

(2) 界面压力设四面体的一个边界面ijm 上受有一线性分布的压力P ，共在三个节点上的强度分别为q i ,0,0。

很容易看出，该力向p 点移置的等效节点力为零。

由水力学知，总压力ijm q P i ∆=31,作用于ijm 面上的d 点，d 点到ij 边和im 边的距离分别为m 到ij 及j 到im 边的距离的1/4。

于是可得{}Tijm i Tei q P P P R ⎥⎦⎤⎢⎣⎡∆=⎥⎦⎤⎢⎣⎡=021211610442(3-62) 所得各节点载荷的方向和分布力的方向相同，要求各节点载荷分量还需乘上相应的方向余弦。

图3-8 重力移置由上述面力移置结果，可求出任意线性分布面的等效节点载荷。

如在ijm 面受有线性分布面力在各点强度分别为q i ，q j ，q m ，时，在i 节点的等效载荷为ijm m j i i q q q P ∆++=)2121(61 (i ，j ，m) (3-63)三、应力应变矩阵空间问题几何方程为{}Tz y x z y x z u x w yw z u x v y u z w y vx u ⎥⎦⎤⎢⎣⎡∂∂+∂∂∂∂+∂∂∂∂+∂∂∂∂∂∂∂∂=⎪⎪⎪⎪⎭⎪⎪⎪⎪⎬⎫⎪⎪⎪⎪⎩⎪⎪⎪⎪⎨⎧=γγγεεεε 将四面体单元之位移表达式(3-52)、(3-54)和(3-55)代入几何方程，即得单元应变。

用节点位移可表示为{}{}[]{}ep mj ie B B B B B δδε--==][ (3-64)式中应变矩阵子矩阵为6×3矩阵：⎥⎥⎥⎥⎥⎥⎥⎥⎦⎤⎢⎢⎢⎢⎢⎢⎢⎢⎣⎡=i ii i i ii i ii b d c d b c d c b V B 00000000061][ (i,j,m,p) (3-65) 由上式可以看出，每一个单元的应变矩阵是一个常量矩阵；因此，采用线性位移模式的四面体单元是常应变单元。

这与平面问题中的三角形单元是一样的。

而与平面问题的不同之处仅在于应变矩阵的阶数不同。

将表达式(3-16)代入空间问题的物理方程，即可得出用单元节点位移表示的单元应力：{}{}[]{}{}eeS B D D δδεσ][][][=== (3-66)式中弹性矩阵[]D 为⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎦⎤⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎣⎡---------=)1(221000)1(2210000)1(221000111111][μμμμμμμμμμμμ称对D 应力矩阵 []p m j iS S S S ][=S (3-67)令 μμ-=11A ， )1(2212μμ--=A则⎥⎥⎥⎥⎥⎥⎥⎥⎦⎤⎢⎢⎢⎢⎢⎢⎢⎢⎣⎡-+-=i ii i i i b A d A c A d A b A V E S 22222i2i i 1i 1i 1i i 1i 1i 1i 000c A d c A b A d A c b A d A c A b )21)(1(6)1(][μμμ (i,j,m,p) (3-68) 显然，式(3-68)中各元素均为常量，应力矩阵[Ｓ]是常量矩阵，所以，四面体单元是常应力单元。