高可乐, 王克用

(上海工程技术大学 机械与汽车工程学院, 上海 201620)

近年来,Trefftz有限元法因兼有传统有限法和边界元法的诸多优点而备受关注,其最早可追溯到1926年Trefftz对Laplace方程的求解.1977年Jirousek等[1]在对薄板体弯曲问题的研究中正式提出杂交Trefftz有限元法(HT-FEM).与传统有限元法不同,杂交Trefftz有限元法采用双位移场插值模式,即单元域内场和辅助网线场.单元域内场采用控制方程的齐次解构建形函数,使得非常稀疏网格和相对少量自由度情况下仍能获得高精度解.此外,仅含边界积分的有限元格式显著提高了Trefftz单元的抗畸变能力.目前,Trefftz有限元法已成功应用于位势问题[2]、平面弹性问题[3]、接触问题[4-5]、轴对称问题[6-8]和复合材料问题[9-12].

汽轮机、火箭、航空发动机等装备中的机械结构常由于工作温度变化引起热荷载和机械荷载双重作用,对其进行热弹性分析十分重要.一般热弹性问题中温度场和变形场是耦合的,从数学角度讲热弹性问题是非齐次方程求解问题.借助Trefftz有限元法求解时,作为体力项的温度荷载会添加域积分到单元刚度方程中,使得杂交基本解有限元法仅含边界积分的优势消失.继区域离散法、双重互易法之后,特解法作为消除域积分的有效方法之一,由Qin等[13]引入到杂交Trefftz有限元法中.此后,Wang等[14]应用径向基函数获得特解研究正交各向异性位势问题;Wang等[15]用特解法求解极小曲面问题;刘博等[16]应用特解法分析轴对称Poisson方程问题的求解;Zhou等[6-8]利用基本解构造单元域内插值函数,提出杂交基本解有限元法(HFS-FEM),并分析轴对称位势问题以及轴对称弹性力学问题.本文基于杂交基本解有限元格式分析热弹性问题,通过应用特解法以消除其中的域积分.

1 控制方程及特解理论

1.1 控制方程和边界条件

线弹性问题的控制方程可写成

(λ+μ)uj,ji+μui,jj+bi=0

(1)

式中:u为位移;bi为非齐次项;λ和μ分别为拉梅常数和剪切弹性模量.逗号表示求导,i,j=x,y,z为笛卡尔空间直角坐标系的3个维度.考虑边界条件为

(2)

(3)

式中:t为表面力;Γ=Γu∪Γt,Γu∩Γt=∅;Γ为求解域Ω的边界;Γu和Γt分别为Dirichlet边界条件和Neumann边界条件,变量上方横线表示给定值.

在边界条件(2)和(3)下,式(1)的解可表达为齐次解和特解的叠加形式.其中,齐次解uh满足

(4)

相应的边界条件修改为

(5)

(6)

特解不受边界条件限制,只需满足关系式

(7)

1.2 离心载荷下的特解

为分析轴对称问题,首先在笛卡尔直角坐标系下给出离心力分量并求出相应位移特解,然后转化为圆柱坐标系下的特解形式.

假设轴对称物体绕z轴以角速度ω旋转,且物体密度为ρ,则体力b分量可写成

(8)

文献[14]分析了控制方程(1)的特解形式,本文给出笛卡尔空间直角坐标系下一个特解为

(9)

将式(9)转化为圆柱坐标系下形式,则位移分量特解可写成

(10)

1.3 温度荷载下的特解

单独温度荷载下,控制方程(1)中非齐次项可写成

bi=βT,i

(11)

将式(11)代入式(1)中,可得到位移分量表示的热弹性平衡微分方程为

(λ+μ)uj,ji+μui,jj=-βT,i

(12)

其中

β=α(3λ+2μ)

(13)

式中:α为热膨胀系数;T为温度.

式(12)的特解可通过热弹性位移势函数Φ(x,y,z)确定.该函数满足关系式

(14)

一旦找到符合条件的热弹性位移势,则位移特解可由式(14)计算得到.为获得上述热弹性位移势,将式(14)代入式(12),化简后得到热弹性位移势与温度场的关系式为

(15)

式(15)为标准Poisson方程形式.若温度场已知,根据式(14)和(15)即可求出位移特解,进而通过位移应变关系和本构关系求出应力特解.需要注意的是,若温度场位移势为简单函数形式,热弹性位移势可通过式(15)直接获得解析式.若温度场解析解比较复杂甚至不存在,则需采用近似函数(如径向基函数或多项式函数)求解近似热弹性位移势[16].

2 杂交基本解Trefftz有限元格式

2.1 双位移场插值模式

杂交基本解有限元法采用双位移场插值模式,如图1所示.单元域内场插值函数通过精确满足控制方程的截断完备解或基本解来构造,而辅助网线场用以保证单元间的连续性,其插值函数采用常规方法建立.

图1 轴对称体与8节点环状单元Fig.1 Axisymmetric body and 8-node annular element

由于控制方程是非齐次的,单元域内场需将特解考虑进去,可写成

(16)

式中:下标e代表单元内变量.辅助网线场不受非齐次项影响,其特解部分包含在单元节点自由度列阵de中,可写成

(17)

LDLTNe=0

(18)

式中:L为微分算子矩阵;D为弹性矩阵,其形式为

(19)

(20)

其中

(21)

根据弹性力学位移应变关系及本构关系,相应应力写成微分算子形式为

(22)

(23)

其中

(24)

辅助网线场可通过自然坐标系ξ∈[-1,1]来构造.对图1所示的8节点四边形单元,每条边布置3个节点,可构造二次网线函数为

(25)

2.2 修正变分泛函

单元域内场和辅助网线场之间的联系是通过修正变分泛函实现的.控制方程为齐次时,通过高斯散度定理可去除变分泛函中的域积分.然而,控制方程为非齐次时,变分泛函中的域积分无法直接去除.为解决这个问题,本文只构造原问题齐次情形的变分泛函,暂时排除边界条件中因非齐次项诱发的特解部分.整个求解域对应的变分泛函可写成所有单元泛函的叠加形式为

(26)

考虑式(18),对式(26)第一项应用高斯散度定理可得

(27)

将式(27)代入式(26),原泛函可简化为仅含边界积分的形式为

(28)

将式(16)、式(17)、式(22)和式(23)代入式(28),可得

(29)

其中

(30)

对式(29)应用两次驻值原理,可得待定系数列阵和单元刚度方程为

(31)

2.3 刚体运动项恢复

将上述单元刚度方程组装成总体刚度方程,采用乘大数法引入位移边界条件即可求得所有节点的齐次解,进一步与获得的特解叠加便可求得节点处的位移全解.需要注意的是,单元域内任一点处的位移还需恢复刚体位移项,具体原因和恢复方法参见文献[7-9],本文因篇幅所限不再复述.

3 数值算例

3.1 离心载荷作用下带孔圆球

为分析离心载荷下物体内部位移分布,给出带圆柱孔的球体结构,如图2(a)所示.由于对称性,只分析1/4结构,圆柱孔半径r0=1,圆球半径为r0+1.本例分析中,弹性模量E=1,泊松比v=0.3,边界条件及网格剖分如图2(b)所示.采用杂交基本解有限元法求解时,需要注意单元域外源点位置.由于这些源点位置与网格单元大小有关,为避免奇异性,子午面内的源点不得配置在z轴左侧.图3给出杂交基本解有限元和ABAQUS位移计算结果,二者吻合良好.

图2 带孔圆球及其网格剖分Fig.2 Perforated sphere and mesh configuration

图3 总位移云图Fig.3 Total displacement nephogram

3.2 温度荷载下厚壁圆筒

本节对厚壁圆筒内热应力的分布情况进行分析.厚壁圆筒尺寸如图4所示.厚壁圆筒内径a=1,外径b=a+1,高h=0.5,内外表面为自由边界条件,上下表面径向固定.弹性模量和泊松比与3.1节中算例相同.

对于热传导问题,若厚壁圆筒内表面温度变化为t1=Ta,t2=0,外表面温度变化为零,厚壁圆筒的温度场解析解可表示为

图4 厚壁圆筒及其网格划分Fig.4 Thick-walled cylinder and mesh configuration

(32)

将式(32)代入式(15),并考虑圆柱坐标系,可得方程为

(33)

对式(33)两次积分,可得热弹性位移势为

(34)

求得热弹性位移势后,位移特解可通过式(14)获得,应力特解通过位移应变关系和本构关系推得.为方便对比,给出此温度场下的应力解析解为

(35)

解析解和杂交基本解有限元解见表1和表2.两者对比,取相同精度时结果十分吻合.

表1 径向应力Table 1 Radial stress Pa

表2 周向应力Table 2 Circumferential stress Pa

4 结 语

针对轴对称热弹性问题,引入热弹性位移势得到一组位移和应力特解.若温度场解析解比较复杂甚至不存在,可根据热弹性位移势导出Poisson方程进而求得近似特解.将特解与杂交基本解有限元列式结合获得齐次解,利用线性叠加原理即可求得全解.本文方法求解过程理论清晰,编程简易,由特解法去除因温度荷载引起的域积分后,杂交基本解有限元法的固有优点得以保持.