复合材料广义热弹性问题格点型有限体积法研究

刘 琦， 肖光明， 齐文亮， 杜雁霞， 刘 磊

(1. 中国空气动力研究与发展中心空气动力学国家重点实验室， 绵阳 621000; 2. 哈尔滨工程大学动力与能源工程学院， 哈尔滨 150001； 3. 中国航空工业集团公司西安航空计算技术研究所， 西安 710068)

1 引 言

在航空航天领域，作为热障涂层而发展的功能梯度材料(FGM)目前已被广泛应用到在各种极端高温环境下工作的机器部件，如压力容器、排气管、核工程以及飞机引擎等. 宏观尺度下，FGM的传热行为和热弹性行为的研究目前已经取得了一定的研究成果. 随着纳米技术的出现，FGM的尺度达到了微/纳米级别.伴随着激光脉冲和微/纳米尺度加工技术的发展，高热流密度冲击下FGM界面挠曲脱粘断裂破坏已成为各国学者关注的热点，准确描述时间极短情形下FGM的热力学行为尤为重要.

近几十年来，广义热弹性问题的研究主要是针对均匀材料进行开展，层和材料以及FGM的广义热弹性问题的研究相对较少. Nowruzpour等[1]基于L-S理论，采用状态空间法以及Laplace变换方法研究了FGM板的广义热弹性响应. 他们假设材料力学性能沿板厚度方向变化，对比了不同热松弛时间以及FGM材料系数下的温度、位移、应力响应分布，热波波前位置处的温度没有明显阶跃. Bahtui等[2]基于L-S理论，联合FEM以及Laplace变化方法研究了热冲击下的圆柱形薄壳FGM体积分数对热弹性响应分布的影响. He等[3]基于G-L理论，联合Laplace积分变换和FEM方法，研究了无限大厚压电板的广义热弹-压电耦合问题，给出了温度、位移以及电势能的分布规律. Sur等[4]基于G-N II型、III型以及三相延迟热弹模型，采用Laplace方法研究了FGM的材料非均匀系数以及热粘对温度、位移、应力及应变的影响，温度在热波波前无显著阶跃，未能体现广义热弹性理论中的热的波动特性. Chen等采用FVM和Laplace变换方法的混合方法分析了材料温度非线性热波问题[5]，并采用相同的计算方法研究了热辐射边界条件下的广义热弹性耦合问题[6]. 以上研究大多采用Laplace变换方法进行求解，然而由于Laplace在逆变换过程中离散和截断误差造成了精度损失，导致热波或弹性波在波前位置处可能没有显著的阶跃[1-4].此外，对于FGM的广义热弹性问题，Laplace变换方法在逆变换过程中存在数值不稳定的问题[7]. 为了避免数值变换方法在变换过程中带来的计算误差以及数值求解过程的复杂性，Tian等[8]发展了时域直接求解FEM用于广义热弹性问题研究，获得了显著波动特性的温度分布，准确反映广义热弹性问题的热波动特性. 基于G-L理论，熊启林等[9]采用时域直接求解FEM对比研究了双层复合板在不同材料参数下界面处的热弹性行为，发现材料参数的不同将导致热穿过界面时界面的温度、位移、应力发生突变. 基于Green-Naghdi理论II型和III型，Xiong等[10]研究了材料热力学性能沿径向变化的FGM固体的热弹性行为，发现热冲击下FGM的材料参数以及G-N理论III型的阻尼系数对热弹性响应有明显的影响，相比于Laplace变换方法，直接FEM计算结果在热波波前有明显的阶跃现象. 基于L-S理论，Antonios等[7]采用直接FEM研究了表面受均匀热冲击下几种陶瓷/金属功能梯度材料的位移、温度以及应力分布，发现FGM材料属性从陶瓷到金属线性过渡时(P=1)，在功能梯度材料界面处的拉应力峰值最小, 但对其它广义热弹理论下FGM界面处的热弹性行为未进行讨论.

迄今为止，在计算流体力学领域发展的有限体积法(FVM)，由于其实施简单，离散方程具有强守恒性，以及离散后方程的各项具有明确物理意义等优点，被国内外学者不断发展并成功地应用到计算固体力学领域[11-15]. 在计算复合材料热弹性问题上，相比于FEM，FVM获得的应力分布更为真实连续，而FEM无论是基于高阶单元或是基于低阶单元, 均存在明显的应力不连续问题[12]. 时域直接求解FEM已经被成功应用于广义热弹性问题的求解，然而，直接时域求解FVM在此方面上的研究工作还相对较少，Liu等[13]基于L-S理论建立了适用于FGM广义热弹问题时域直接求解的FVM，结果表明FVM可以和直接FEM一样有效捕捉热波动效应下温度在波前处的阶跃现象，准确反映广义热弹性问题的波动特性. 本文基于G-L及T-C理论对FVM进一步扩展，讨论不同热弹耦合理论下FGM材料系数对界面处的温度、位移和应力响应的影响.

2 控制方程

基于任意控制体Ω建立热传导方程和平衡方程，其中，控制体Ω的体积为V，控制体的边界面为S.研究稳态下各向同性线弹性复合材料热应力问题.

2.1 热传导方程

Ω中能量守恒方程为：

(1)

式中，q为热流密度矢量;Q为体积热源生成热;ST为体积熵;T0为参考温度.

熵平衡方程为[14]：

(2)

热传导方程为[14]：

(3)

式中，t0为CV模型引入的热松弛时间;k为导热系数矩阵.

(4)

式中，ni(nx,ny)为微元面外法矢量.

边界条件：

θ=θB,边界SD

(5)

q=qB,边界SN

(6)

式中，下角标“B”表示边界;θB为Dirichlet边界上给定的温度;qB表示Neumann边界上给定的热流.

温度场初始条件为：

(7)

2.2 热弹性方程

Ω中运动平衡方程为：

(8)

式中，u为位移矢量;σ为应力张量;f为单位质量体积力矢量.

线弹性体本够关系为[14]：

(9)

针对各向同性体，弹性模量张量C可简化为：

(10)

应变-位移关系为：

(11)

式(9)和式(10)中，t1为G-L理论引入的热松弛因子；C为四阶弹性模量张量；μ、λ和Γ分别为材料的拉梅系数和热弹性系数，可采用材料弹性模量E、泊松比υ以及热膨胀系数αt进行表示.

(12)

将式(9)～式(11)带入式(8)，可得：

(13)

将式(13)写为张量形式并在控制体Ω内进行积分，得T-C、L-S和G-L理论热弹性问题待解平衡方程：

(14)

边界条件：

ui=uiB

(15)

σijnj=σnB

(16)

式中，uiB为Dirichlet边界上给定的位移;σnB为Neumann边界上给定的力载荷.

位移场初始条件为：

(17)

3 数值方法

本文采用双线性四边形单元划分计算域. 图1为交错网格示意图，其中实心点为单元节点，空心点为单元中心. 我们将待解变量(温度)存储在单元节点，并假设其在单元内分布一致，其余物理量(密度、比热、导热系数、热松弛时间、弹性模量等)存在单元中心，并假设在单元内分布一致；以内部节点N1为例，连接相临边中点以及单元中心建立控制体Ω1～6.

图1 单元网格及控制体示意图[15]

3.1 热传导方程离散

借鉴FEM，单元内任意变量φ采用形函数进行插值近似：

(18)

式中Ni为形函数;ns为单元内的节点数，对于4点四边形单元ns=4. 四边形单元的形函数表达式为：

(19)

式中ξi,ηi为四边形单元所对应的局部坐标.

为了更清晰地表达能量方程式，将其改写成如下格式：

(20)

将式(19)带入式(20)，将其离散为：

(21)

式中，nc为节点n周围的单元总数;θij为第i个单元的第j个节点的温度变化量;Si为第i个单元内控制体边界面的长度.

式(20)中左端一阶时间项采用向后差分法进行离散、二阶时间项采用向中心差分法进行离散、其余项采用节点n处的值进行近似，同时引入边界条件，则式(4)离散为：

(22)

其中，

(23)

其中，|ABi|为处于边界上第i个单元的面积矢量的模;nN为处于Neumann边界上的单元个数.

3.2 热弹性方程离散

采用类似处理方法，热弹性方程离散为：

(24)

式中

(25)

式(23)和式(25)中有关形函数及其空间导数的积分项在标准坐标系下进行计算，具体计算流程可见文献[15].

为方便数值计算，引入以下无量纲参数：

(26)

式中，CE表示弹性波波速;下角标“m”表示金属基体.特征长度ξm及特征速度cm为：

(27)

为方便表示，下文不显示无量纲角标“*”.

采用FORTRAN语言实现本文计算方法，开始计算之前，将与网格信息相关的形函数导数、积分点位置、形函数积分作为几何常数进行一次性存储，减少计算量．对于临近单元存在Dirichlet边界的节点采用置大数法将对其导热系数矩阵或刚度矩阵相应的行进行处理.离散后的热传导及热弹性方程采用直接求解法进行求解.

4.1 热弹性耦合理论CV-FVM算法验证

针对一半无限大方板热冲击问题[16-18]进行数值验证，采用三种耦合理论进行求解. 计算模型如图2所示，计算过程中不考虑热波和弹性波在方板边界处的反射. 方板左侧受θ=H(t)θ0，θ0=1的热冲击作用，其余三个边界绝热，方板右侧施加固支约束. 为节省计算资源，将计算域沿x轴方向分为两部分，仅将区域1内(0≤x≤2)的网格进行加密，采用均匀四边形网格进行网格划分，网格密度为60×30. 区域2内采用非均匀四边形网格进行划分，网格增长率为1.1，网格密度为60×30，网格模型如图2b所示. 材料选取为不锈钢，其物性参数见表1. 根据波动方程性质，无量纲热波和弹性波波速分别为CT=1/(t0+t2)1/2，CE=1. 针对T-C理论下热弹性问题，选取δe=0和δe=1进行算法验证；针对L-S和G-L理论下热弹性问题选取δe=0，δe=0.05和δe=0.36进行算法验证.

(a) 计算模型

(b) 网格模型

计算初始条件为：

(28)

表1 不锈钢物性参数

4.1.1 T-C理论CV-FVM算例验证 图3给出了方板x=1处(弹性波t=1时刻波前位置)的无量纲温度、位移以及应力随时间的变化规律. 从图3可以看出CV-FVM结果与FEM和解析解[16]得到的结果吻合良好. 相比于无耦合计算结果，热弹性耦合作用导致温度、位移以及应力响应峰值明显降低. 从图3a可以看出，耦合作用加速了弹性波波前位置处温度扩散速度，在弹性波波前产生了一个“负”的温度梯度，导致区域(x=1)内温度降低.在弹性波经过该区域后温度扩散速度逐渐降低，温度逐渐恢复并趋近无耦合作用的温度结果. 发生上述现象的原因是热能和机械能转换主要发生在波前位置处. 从图3b, 3c可以看出，相比于CV-FVM结果，FEM得到的位移和应力幅值明显偏低，同时应力在弹性波波前位置处存在明显的振荡现象.

4.1.2 L-S、G-L理论CV-FVM算例验证 图4为方板沿x轴方向不同时刻温度分布. 从图4可以看出，L-S和G-L理论的热波在波前均有明显突跳现象，二者预测的温度形式无明显区别，随着时间的推进，热波强度逐渐降低. 不同耦合系数下L-S和G-L理论下的得到温度曲线中均出现了两个激烈的温度波动，分别由弹性波(靠前)和热波(靠后)引起. 随着耦合系数的增大，弹性波引起的温度扰动越提前，热扰动越强烈，影响范围越大，同时热波引起的温度扰动幅值降低.

图5 x=1处温度T随时间的变化规律(无量纲)

图6和图7为基于L-S和G-L理论下沿x轴方向和x=1位置处不同时刻位移分布. 从图6和图7可以看出，L-S和G-L理论下的位移响应幅值均随着耦合系数的增大而降低，但二者位移响应分布规律明显不同. G-L理论下的位移响应较L-S理论下的结果更为剧烈，对应位移峰值更高. G-L理论下应力松弛时间t1的引入明显改变了位移响应的分布形式，使得温度波动对位移和应力的作用更加明显. CV-FVM预测的位移响应与BEM结果吻合较好.

图6 沿x轴方向不同时刻位移u的变化规律(无量纲)

图7 x=1处位移u随时间的变化规律(无量纲)

图8和图9给出了方板沿x轴方向和x=1处的应力分布规律. 从图8和图9可以看出，耦合系数对L-S和G-L理论下的应力分布的影响与对位移分布的影响相一致，相比于L-S理论计算结果，G-L理论下的应力响应更为剧烈，峰值更高. 图9b给出了Prevosr[16]及Chen[18]耦合系数δe=0时计算结果，对比发现，基于G-L理论采用CV-FVM预测的应力结果与文献预测结果分布规律相同，仅在应力峰值上存在差别.

图8 沿x轴方向应力σxx随时间的变化规律(无量纲)

图9 x=1处应力σxx随时间的变化规律(无量纲)

4.2 FGM复合板热弹性行为的研究

考虑文献[7]陶瓷/金属功能复合板上端面(x=0)受热冲击下瞬态热冲击问题(图10)，本文采用G-L、T-C理论对其进行进一步研究. 冲击强度θ=1. 计算时我们考虑了热波和弹性波在矩形板下端面(x=1)的反射. 模型无量纲计算域为H×L=1×10，我们采用均匀四边形单元对其进行网格划分. 我们采用平面应变假设进行求解，计算过程中不考虑界面热松弛时间t0的空间梯度对热弹性分布的影响. 计算初始条件为：

(29)

图10 陶瓷/金属复合板热冲击作用下热弹性问题

表2 涂层、基体物性参数[13](参考温度为300 K)

物性参数采用Voigt混合律模型进行计算，其有效物性参数Peff(x)表示为：

Peff(x)=PmVm(x)+Pc(1-Vm(x))

(30)

式中，Pm为金属材料物性;Pc为陶瓷材料物性;Vm为沿x方向金属材料体积分数.Vm的表达式采用Sigmoid函数给出：

(31)

式中，p为描述材料空间变化的指数系数.数值模拟选取p=0.1, 1, 10.

任意层内的弹性波波速CE和热波波速CT为：

(32)

根据式(32)可以得到双层金属/涂层复合板内的弹性波CE和热波波速CT(表3).

表3 全陶瓷/金属双层复合板金属域和陶瓷域内特征长度、特征速度、热波速度和弹性波速度

基于L-S理论，图12给出了p=0.1,p=1和p=10时Ti-6Al-4V/ZrO2复合功能型材料交界面x=0.5处的温度、位移以及应力随时间的变化曲线. 从图14a和14b可以看出，随着材料指数系数p的增大，热波传播速度明显降低，热强度峰值明显降低，同时热波达到稳定时的温度降低. 产生这种现象的原因主要是随着p的增大，在0x0.5区域内导热系数降低，沿x方向温度梯度增大，导致x=0.5处温度幅值降低(如图15). 从图14b可以看出，随着p的增大，x=0.5处的位移响应降低.相比于位移响应，温度扰动更快达到稳定. 从图14c应力响应分布可以看出，p=1时复合材料界面处的拉应力峰值最小，这和文献[7]中应力变化规律一致.

图10 不同时刻位移u的分布(无量纲)

图11 不同时刻温度T的分布(无量纲)

为了体现研究的广义热弹性问题的热波动以及热-机耦合特性，在t=0时刻陶瓷端施加压力F=1的机械冲击并与热冲击载荷计算结果进行对比，对比结果如图12和图13所示.

从图12中可以看出，机械冲击仅在弹性波波前产生应力阶跃，同时机械冲击由于热弹性耦合作用而引起的温度波动十分微小，但仍然存在着由弹性波和热波引起的温度阶跃现象. 从图13中可以看出热冲击载荷下应力波和热波波前存在两次阶跃，温度响应在热波波前位置处仅存在一次阶跃. 热冲击作用下，冲击作用面(y=0)由于温度的升高产生膨胀作用，同时由热波诱导产生的机械波也会产生压缩作用，导致热冲击作用面具有最大的膨胀位移，并随着x的增加不断减小，变成压缩位移，最终在弹性波波前位置处消失.

图16给出了G-L理论下的计算结果. 假设G-L理论下的热松弛时间t1=t2=t0，从图15a可以看出G-L理论下的温度响应和L-S理论下的温度响应(图14a)基本一致，当p=0.1时的温度响应峰值最大. 从图16b可以看出，G-L理论下的位移响应结果明显高于L-S理论计算结果(图14b)，同时，当p=1时的位移的响应峰值最大. G-L理论下的应力响应峰值的预测规律与L-S理论明显不同，这是由应力松弛时间t1的引入所致. 当p=10时G-L理论下的拉应力响应峰值最小.

图17给出了T-C理论下温度、位移以及应力响应计算结果. 由于T-C理论假设热波的传播速度为无限大，x=0.5处的温度很快趋于稳定，与L-S理论和G-L理论预测的温度响应规律一致. 在0

图15 导热系数k沿x轴方向的变化规律

图17 x=0.5处温度T、位移u和应力σxx随时间的变化规律(无量纲)

表4 Ti-6Al-4V/ZrO2复合板x=0.5处温度T、位移u和应力σxx的分布(无量纲)

从图14、图16和图17的对比结果中可以看出，基于三种耦合理论的计算结果存在明显区别. 当材料的梯度系数p增大时，三种耦合理论下的温度峰值变化一致，随着陶瓷层内陶瓷体积分数p增加，温度响应峰值降低，同时温度扰动达到稳定的速度要快于位移扰动达到稳定的速度. 不同指数系数p下，G-L理论下的位移响应和应力响应峰值最大. L-S理论下p=1时交界面处的拉应力峰值最小；G-L理论和T-C理论下p=10时交界面处的拉应力峰值最小(表4). 相比与L-S和T-C理论下的计算结果，G-L理论由于热松弛时间t1的引入导致位移和应力响应幅值显著增强，幅值变化更为剧烈.

5 结 论

本文基于L-S、G-L、T-C耦合理论, 发展了一种适用于复合材料广义热弹性问题时域直接求解的CV-FVM方法. 基于双线性四边形单元，对控制方程的离散过程进行了详细推导. 针对均质无限大方板热冲击问题对CV-FVM的求解过程进行了数值验证. 计算结果表明，本文发展的数值方法可以很好地捕捉热波波前和弹性波前的阶跃特性，以及热弹耦合特性. 相比于FEM，本文的计算应力分布在弹性波波前位置处更为稳定. 本文采用三种热弹性耦合理论对不同梯度系数p下钛合金/氮化硅复合板的热冲击问题进行了对比研究，针对本文算例，我们发现L-S理论复合板p=1时预测的界面应力最小；T-C、G-L理论下，复合板p=10时预测的界面应力最小. G-L理论由于热松弛时间t1的引入导致位移和应力响应幅值显著增强，应力变化更为剧烈.