谷小强， 常利武， 孙玉周， 张旺玺

(1.中原工学院 建筑工程学院， 河南 郑州 450007； 2.中原工学院 材料与化工学院， 河南 郑州 450007)

随着科技的进步和问题研究的深入，许多原来能满足精度要求的理论由于限制条件的增加而不再能满足，进而发展为新的理论[1]。如伴随着热以有限速度传播的热机耦合理论的提出，经典热弹性理论发展为广义热弹性理论[2-3]；考虑材料的微结构特征对热弹性介质的影响，在控制方程中引入偶应力张量，广义热弹性理论发展为热弹性偶应力理论[4]。新近大量的研究表明，热弹性介质材料在超高温度梯度、超高温升速率等特殊环境中表现出复杂的热力耦合特征，温度梯度的存在必然对温度梯度固体的应力场产生影响[5]。由于温度梯度会使物质结构产生偶应力，故可将温度梯度产生的偶应力引入到热弹性本构关系中,以研究其对结构力学性能的影响。

热弹性偶应力理论虽然考虑了应变梯度对结构的影响，但并没有考虑温度梯度对应力场的影响。本文将温度梯度引入到偶应力热弹性理论的本构关系中，建立了考虑温度梯度影响的无网格数值计算方法，并编写了相应的Fortran计算程序，对温度梯度影响下的结构性能进行了研究。

1 无网格法计算框架

1.1 本构方程

若只考虑稳态热传导问题，可以认为温度梯度只与偶应力有关。将温度梯度引入到偶应力理论的本构关系中，即可得到考虑温度梯度影响的偶应力理论的本构方程。对于各向同性的线弹性热传导固体，Helmholtz自由能函数Φ可表示为应变张量、应变梯度张量、温度以及温度梯度的函数[4-5]

(1)

式中：λ、μ是拉梅常数；η、η,i为偶应力介质的材料常数；β是热力耦合系数，β=(3λ+2μ)α；α是材料的线性膨胀系数；εji为应变张量；χji为应变梯度张量；ΔT表示温度的变化量，ΔT=T-T0；c是材料的比热容。

参考张晓敏等提出的方法[5]，考虑温度梯度影响的偶应力理论的本构方程可表示为：

(2)

(3)

式中,γ为介质常量。本文取γ=1.68×10-4N/K。

1.2 无网格法计算框架

与传统数值计算方法如有限元方法、边界元法相比，无网格方法具有许多突出的优点，如移动最小二乘近似的形函数具有高阶连续的特征[6-8]。利用形函数的高阶连续特征，温度梯度、结构的应变和旋转曲率可直接由节点的温度或位移的插值得到，这样，问题大大简化。

在无热源的情况下，由于温度发生改变，弹性结构体内部产生的势能由结构的弹性应变能提供。系统的弹性应变能为：

(4)

其中:

[εT]T=[αΔTαΔT0αT,iαT,i]

(5)

[εE]T=[ε11ε22γ12χ1χ2]

(6)

[ε]T=[εE]+[εT]

(7)

式中，[D]为刚度矩阵。

问题域上的位移函数u(x)可表示为

u(x)=uh(x)=φ(x)u

(8)

其中:φ(x)为形函数，u为节点的位移向量。

应变分量ε可写成：

(9)

uT=[u1v1u2v2…unvn]

(10)

利用移动最小二乘近似，温度场可表示为：

T(x)=Th(x)=φ(x)T

(11)

温度梯度可表示为：

(12)

TT=[T1xT1yT2xT2y…TnxTny]

(13)

将式(13)、式(9)和式(5)代入式(4),可得

(15)

其中：

(16)

(17)

(18)

式中，C只与温度有关，不产生势能，可略去。式(15)可简化为

(19)

令δΠ=0，得无网格法的离散控制方程为

[K][u]=[F]

(20)

式中：

(21)

(22)

2 算例及结果分析

图1所示为矩形薄板，长度L=20 cm，高度H=5 cm，厚度t=0.01 cm，材料的弹性模量E=2.6 GPa，泊松比ν=0.32，线性膨胀系数α=1.0×10-5/K。对薄板施加的连续温度场的边界条件分别取TL=0，TU=10°和TL=0，TU=40°两种工况。数值计算中，薄板模型的总节点分布为6×21，背景网格采用3×3的高斯积分方案，且温度场和位移场采用相同的节点分布和背景网格积分方案。

图1 薄板模型

考虑温度梯度影响的偶应力热弹性力学问题的研究是在偶应力的基础上研究温度梯度对结构力学性能的影响，故还需要考虑微结构的尺度效应问题。本文通过含节理特征的层状岩体的尺度因子(l=0.4t)来表述微结构的影响[9-10]。

图2和图3分别为偶应力热弹性力学下的薄板弯曲变形图和考虑温度梯度影响的偶应力理论下的薄板弯曲变形图。从图中可以看出，两种情况下结构内部的弯曲程度从边缘到中间均不断加大，且考虑温度梯度影响的薄板的弯曲变形得到了加强。

图2 偶应力理论的热弹性位移云图(TL=0，TU=10 ℃)

表1所示为边界条件TL=0，TU=10°和TL=0，TU=40°两种工况下薄板上边缘5个节点的Y方向位移值(括号里的数值表示边界条件为TL=0，TU=40°时的位移值)。由表1可知，薄板弯曲变形从边缘部位向薄板中间隆起并在薄板中间达到峰值，考虑温度梯度影响的偶应力理论下得到的位移峰值较偶应力热弹性理论下的结果，在两种温度梯度下分别提高了2.5%和5.8%。这说明温度梯度的存在增加了薄板的弯曲变形值，其提高幅度随着温度梯度的增大而增大。

图3 温度梯度偶应力位移云图(TL=0，TU=10 ℃)

表1 Y方向位移值

3 结论

在假设温度梯度只产生偶应力的基础上，将其引入到偶应力热弹性理论的本构方程中，得到了一个简化的考虑温度梯度影响的偶应力本构方程，建立了无网格数值计算框架，并编写了相应的Fortran计算程序，对温度梯度影响下薄板的弯曲变形规律进行了研究。研究结果表明，利用移动最小二乘近似建立起来的无网格方法可以较好地模拟考虑温度梯度和应变梯度的高阶连续问题，薄板在考虑温度梯度影响时有显著的变形，且变形值随着温度梯度的增大而增大。本文结论为考虑应变梯度场和温度梯度场影响的结构的数值模拟提供了一种新的思路。