梯度材料被广泛应用于航空航天、激光辐照等领域^[1-3]。材料性能的梯度分布有效解决了性能失配、界面剥离和热应力集中等问题^[4-5]。高熔点陶瓷和金属组合而成的梯度材料能够承受超高的温度梯度，是理想的阻热应用材料之一^[6]。研究梯度材料在高温冲击下的热传导过程，对于极端环境热分析具有重要意义。梯度材料的梯度分布形式有指数函数、幂函数以及改进的sigmod函数等，其中幂函数是梯度材料设计制备中广泛采用的形式之一^[7]。

根据傅里叶经典热传导理论^[8]：$q = - K\nabla T$，其中，$q$为热流密度矢量，$\nabla T$为温度梯度，K为导热系数。傅里叶热传导方程是对稳态热传导现象的一种唯象描述，默认满足其静态、低热流条件，抛物型控制方程对应温度扰动的传播速度趋于无穷。在瞬时高温加载工况下，传统的傅里叶热传导模型不再适用。实验研究发现，高温冲击下，靶板表面温升比傅里叶理论预测值高300 K^[9]。微秒激光脉冲加热实验发现，当热扰动足够强，在多孔材料的局部区域可以观察到明显的非傅里叶热传导现象^[10]，热流与温度响应之间存在时间差。Cattaneo和Vernotte在传统傅里叶热传导模型中引入了弛豫系数，提出了双曲型热传导方程，通过非平衡态热瞬时热力学方程来描述非傅里叶热传导现象^[11]，即$q + \tau \dfrac{{\partial q}}{{\partial t}} = - K\nabla T$，其中，热流的时间偏导项系数τ具有时间（t）的量纲，定义为热弛豫时间，用来表征温度场对扰动的响应。与传统的抛物型方程不同，双曲型热传导方程描述的温度扰动具有有限的传播速度，并引入了“温度波”这一概念。Bucă等^[12]研究了飞秒激光脉冲作用下金属中电子温度的变化，弛豫时间为皮秒量级。实验研究进一步验证了非傅里叶热传导中“温度波”的存在^[13-15]。数值计算结果表明，利用非傅里叶热传导模型预测激光加载的热传导是准确和稳定的^[16]。Peng等^[17]研究了无限长圆柱体在表面瞬时热载荷作用下的瞬态热弹性响应，利用拉普拉斯变换和解耦技术，得到了温度波波速。Jiang^[18]给出了均质材料空心圆柱体结构的双曲型热传导方程的解析解。Keles等^[19]对空心圆柱体和球体结构的非傅里叶热传导问题进行了研究，得到了指数分布的功能梯度材料的数值解。

目前，对于梯度材料双曲型热传导方程的求解主要以数值结果为主，解析解较少，同时研究的梯度结构较为单一。本文中，首先，对幂律梯度材料的一维非傅里叶双曲型温度传导方程进行求解，得到温度场关于时间和空间的理论解析解；然后，以Mo-ZrC梯度金属陶瓷材料的应用为例，讨论一般温度边界条件以及温度脉冲载荷作用下材料内部温度场的分布，研究梯度材料的热弛豫响应行为特征及其影响因素。

1.   幂律梯度材料的双曲型热传导方程

1.1   Mo-ZrC梯度材料

金属Mo具有耐高温和耐腐蚀的特性，ZrC陶瓷具有高熔点、低密度和良好的化学稳定性。两种材料制备的梯度复合材料兼具优秀的耐高温特性和抗热冲击性能，被广泛应用于发射器表面涂层、热光电辐射器等超高温耐火材料领域^[20-21]。金属Mo和ZrC陶瓷的热力学性能如表1所示。

表  1  Mo和ZrC的热力学参数

Table  1.  Thermomechanical properties of Mo and ZrC

材料	密度/(kg·m−3)	比热容/(J·kg−1·K−1)	导热系数/(W·m−1·K−1)	熔点/℃	杨氏模量/GPa	泊松比
Mo	10200	230	150	3400	279	0.32
ZrC	6510	310	10	2620	390	0.191

下载: 导出CSV 

| 显示表格

假定梯度材料性能满足幂律分布，即其任意位置的热力学性能P（如密度、比热容、导热系数等）分布满足：$P=P_0\left(x/l_0\right)^{\beta\mathrm{_{\mathit{p}}}n}$，其中，P₀表示ZrC的热力学性能，l₀为材料厚度，x、β_p、n分别表示位置坐标、材料组合系数和梯度指数，x/l₀∈[1,2]。x=l₀对应ZrC材料，x=2l₀对应Mo材料，中间为梯度部分。材料组合系数通过表1中金属Mo和ZrC陶瓷的热力学参数计算得到，β_p=(1/n) log₂(P₁/P₀)，其中P₁表示Mo的热力学性能。令y=(P−P₀)/(P₁−P₀)，简化指数m=β_pn。归一化处理后的幂函数形式为y=[(x/l₀)^m−1]/(2^m−1)，其分布如图1所示，可以根据不同的材料性能组合得到相应参数进行后续计算。

图  1  归一化幂函数分布

Figure  1.  The distribution of normalized power function

下载: 全尺寸图片 幻灯片

1.2   幂律梯度双曲型热传导方程

线性双曲型热传导方程^[19]如下：

式中：ρ为密度，c为比热容，Q为内部热源，T为变温。对于一维问题，假设内部热源Q=0，则热传导方程为：

对其进行无量纲处理，得到：

式中：$\bar q = \dfrac{q}{{{q_0}}}，\bar T = \dfrac{T}{{{T_0}}}，\bar t = \dfrac{t}{{{t_0}}}，\bar \tau = \dfrac{\tau }{{{\tau _0}}}，\bar \rho = \dfrac{\rho }{{{\rho _0}}}，c = \dfrac{c}{{{c_0}}}，\bar x = \dfrac{x}{{{l_0}}}，\bar K = \dfrac{K}{{{K_0}}}$；q₀为参考热流，T₀为环境温度，ρ₀、c₀、K₀分别为陶瓷材料的密度、比热容和导热系数；t₀、τ₀分别为参考时间以及参考热弛豫时间，且满足t₀=τ₀，$t_0=\rho_0 c_0 l_0^2/K_0$，q₀=K₀T₀/l₀。

假设边界条件和初值条件满足如下形式：

式中：$h(\bar t)$表示温度随无量纲时间的变化函数。对热传导方程时间项进行拉普拉斯变换，得到：

式中：“^～”表示与时间相关的无量纲化参数经过拉普拉斯变换后的象函数，s为频域变量。式(5)消去$\tilde q$得到：

式中：下标x表示对x的导数，下标xx表示对x的二阶导数。假设弛豫时间系数为常数，即$\bar \tau = {k_\tau }$。梯度材料的热力学性能满足幂律分布，即$\bar \rho \propto {x^{{\beta _\rho }n}},\bar c \propto {x^{{\beta _c}n}},\bar K \propto {x^{{\beta _k}n}}$，其中β_ρ、β_c、β_k分别由表1中的密度、比热容和热导系数计算得到，代入式(6)化简得到：

转化为一般形式的Lommel方程^[22]：

对应的系数如下：

式(7)所示的双曲型非傅里叶热传导方程解为：

式中：$b = {b_1}g(s)，g(s) = \sqrt { - s(1 + {k_\tau }s)}，{b_1} = \dfrac{2}{{n({\beta _\rho } + {\beta _c} - {\beta _k}) + 2}}$；系数C₁(s)、C₂(s)由边界条件确定；J为第一类贝塞尔函数；Y为第二类贝塞尔函数。当$n({\beta _\rho } + {\beta _c} - {\beta _k}) + 2 \ne 0$时，将得到的解析解定义为第一类解；当$n({\beta _\rho } + {\beta _c} - {\beta _k}) + 2 = 0$时，将得到的解析解定义为第二类解。

2.   第一类解析解求解

2.1   一般温度边界条件求解

对于如下形式的温度边界条件（$\bar t {\,\text{＞}} 0$）：

进行拉普拉斯变换后得到：

结合式(10)和式(12)求解得到${C_2} = - {C_1} \dfrac{{{J_v}({2^{{c_1}}}b)}}{{{Y_v}({2^{{c_1}}}b)}}$，${C_1} = \dfrac{{{T_1}}}{s}\dfrac{{{Y_v}({2^{{c_1}}}b)}}{{{J_v}(b){Y_v}({2^{{c_1}}}b) - {Y_v}(b){J_v}({2^{{c_1}}}b)}}$。代入式(10)得到：

令$F(\bar x,s) = {J_v}(b{\bar x^{{c_1}}}){Y_v}({2^{{c_1}}}b) - {Y_v}(b{\bar x^{{c_1}}}){J_v}({2^{{c_1}}}b)$，则频域下温度解析表达式可以写为：

式(14)存在一阶极点s=0和普通极点s_k。其中，s_k通过下式求得：

通过极点留数法，对$\tilde T$进行拉普拉斯逆变换，可以解得：

2.2   幂律梯度Mo-ZrC材料的一般温度边界求解

取梯度指数n=1，根据表1所示的金属Mo和ZrC陶瓷的热力学参数，可以计算得到：β_ρ=0.65，β_c=−0.43，β_k=3.91，a=−1.45，b₁=−1.18，c₁=−0.85，v=1.72。

图2为v=1.72时的贝塞尔级数及其渐近展开式以及函数F(1,X)，横坐标X表示自变量，纵坐标F_X表示对应的函数值。可以看出，在第一个零点位置，贝塞尔级数接近其渐近展开式，贝塞尔级数渐近展开式如下：

图  2  贝塞尔级数渐近展开可行性验证

Figure  2.  Feasibility verification of asymptotic expansion of Bessel equation

下载: 全尺寸图片 幻灯片

两者相乘得到：

因此：

式中：$A = {\bar x^{{c_1}}}b$，$B = {2^{{c_1}}}b$。代入式(19)和式(15)，化简得到：

代入式(14)得到频域温度的表达式：

式(21)存在一阶极点s=0和普通极点s_k，s_k满足方程$b({2^{{c_1}}} - 1) = k\pi$$(k = 0,1,2,3, \cdots)$。解得：

利用极点留数法求解式(21)的逆变换，对于一阶极点$s = 0$，有：

对于普通极点${s_k}$，有：

整理可得：

2.3   幂律梯度Mo-ZrC材料的温度脉冲边界求解

对于温度脉冲边界条件，假设脉冲幅值为T₁，持续时间为t_d，有：

式中：H为赫维赛德阶跃函数。H满足：

利用线性叠加原理，可以得到温度脉冲下无量纲温度$\bar{T}_{\mathrm{d}}$的表达式为：

3.   温度函数分布及弛豫时间系数讨论

3.1   数值解对比及解析解分析

假设函数f(t)满足拉普拉斯变换存在的条件，F(s)是它的象函数，则当t>0时^[23]：

式中：i为虚数单位；积分是沿着任一直线$\mathrm{Re}(s)=\beta$进行积分，积分值不依赖β。可以选择合适的β值使得F(s)所有的极点在Re(s)=β左侧。令s=β+ωi，其中ω为替换变量，得到：

利用梯形积分公式进行数值积分运算，数值结果和解析结果如图3所示，其中：红色虚线为数值计算结果，黑色实线为解析式计算结果。弛豫时间系数k_τ=1，$\bar x = 1.5$。数值计算结果在间断位置存在较为明显的数值振荡。数值结果和解析结果的整体趋势一致，证明两种计算方式均可靠。

图  3  拉普拉斯逆变换数值计算结果和第一类解析结果对比

Figure  3.  Comparison of numerical results of inverse Laplace transformation and the results of the first analytic solution

下载: 全尺寸图片 幻灯片

如图3绿色虚线框所示，当弛豫时间系数k_τ=1时，温度扰动以梯形波的形式传播，该结果与Keles等^[19]的分析一致。框内为一个波形单元，AB为线性下降段，BC为间断上升段，CD为线性上升段，DE为间断下降段。因为式(1)中温度对时间的一阶偏导项引入了阻尼项，所以随着时间的发展，温度-时间曲线的振幅逐渐衰减。

3.2   解析函数形貌分析及系数讨论

取弛豫时间系数k_τ=1，温度边界T₁=1，得到温度随时间和位置的二维分布如图4所示。一侧（$\bar x = 1$）的温度扰动以温度波的形式向另一侧（$\bar x = 2$）传播，温度场表现出波动和传导的双重特性。

图  4  第一类温度函数的二维分布图像

Figure  4.  Two-dimensional distribution image of temperature function of the first analytic solution

下载: 全尺寸图片 幻灯片

图5为不同位置的温度-时间曲线。对于边界温度扰动，温度场需要一定的响应时间。将温度场的响应时间t_r定义为温度-时间曲线的第一个起跳位置。可以发现，随着坐标位置逐渐远离扰动源，整体对应的扰动响应时间变长，$\bar{x}$为1.3、1.5、1.8对应的响应时间t_r分别为0.18、0.34、0.42。如图5绿色虚线框所示，随着坐标位置远离扰动源，线性下降段AB变长，线性上升段CD变短，振动幅值减小。CD段中心大致处于同一时间，不同位置的温度场振荡周期趋于一致。

图  5  第一类解在不同位置处的温度-时间曲线

Figure  5.  Temperature-time curves at different positions of the first analytic solution

下载: 全尺寸图片 幻灯片

图6为不同时刻温度的空间分布曲线。如图6(a)所示，在初始阶段，温度波从一侧传播至另一侧，0.6时刻前在边界位置发生反射，从而反向传播，波速随坐标位置发生变化。如图6(b)所示，随着时间发展，在1～6的时间内，温度波振动衰减，温度场的空间分布趋于平衡，平衡态与边界条件以及材料梯度分布相关。

图  6  第一类解在不同时刻的温度分布曲线

Figure  6.  Temperature-position curves at different times of the first analytic solution

下载: 全尺寸图片 幻灯片

图7(a)给出了弛豫时间系数k_τ分别为0.25、1、4时，相同位置（$\bar x = 1.5$）处的温度-时间曲线。随着弛豫时间系数的增加，同一位置处的温度响应滞后更加明显，温度波动幅值增大，脉宽增加。单位波形从梯形向矩形发展，由图7(a)中线性上升段C₁D₁向C₂D₂发展。统计响应时间和弛豫时间系数k_τ（图7(b)）可以发现，随着弛豫时间系数的增大，响应时间增加，同时温度极值（$\bar T_{\mathrm{max}}$）也随之增大。

图  7  第一类解在不同弛豫系数下的温度-时间曲线以及响应时间、温度极值与弛豫系数的关系

Figure  7.  The temperature-time curves corresponding to different relaxation coefficients and the relationships between response time, maximum temperature and relaxation coefficient of the first analytic solution

下载: 全尺寸图片 幻灯片

双曲型热传导行为具有波动和传导的双重特性。随着时间的发展，波动逐渐衰减，最终趋于平衡态。图8为相同位置（$\bar x = 1.5$）、相同弛豫时间系数（k_τ=1）时，3种脉冲（t_d=0.5，1.0，2.0）的温度-时间曲线。可以发现，温度脉冲作用下，温度时间曲线表现为波动特性，且波动逐渐衰减，温度趋于稳定。随着脉冲持续时间的不同，温度场的波动形式有所改变。

图  8  不同温度脉冲对应的温度-时间曲线

Figure  8.  Temperature-time curves corresponding to the temperature pulses of different widths

下载: 全尺寸图片 幻灯片

4.   第二类解析解求解

4.1   简化欧拉方程求解

当$({\beta _\rho } + {\beta _c} - {\beta _k})n + 2 = 0$时，式(7)可以简化为欧拉方程形式。由一般Lommel方程得到的解析解为第一类解，由欧拉方程得到的解为第二类解。化简后的欧拉方程形式为：

式(31)具有$T = {\bar x^r}$形式的特征解。代入特征解求得特征根满足：

解得：

式(31)的通解为：

假设β_ρ=β_c=β_k=1，n=−2，简化特征方程，得到${r_{1,2}} = \dfrac{{3 \pm \sqrt {9 + 4s(1 + {k_\tau }s)} }}{2}$，则：

式中：$f(s) = \sqrt {2.25 + s(1 + {k_\tau }s)}$。另一种形式通解为：

式中：双曲函数解的系数D₁(s)、D₂(s)由边界条件确定。对于温度边界条件（式(11)～(12)），求解系数得到${D_2} = - {D_1}\tanh (f(s)\ln 2),{D_2} = \dfrac{{{T_1}}}{s}$。

频域内温度的表达式为：

式(37)存在一阶极点s=0和$f(s)\ln2=\mathrm{i}k{\text{π}}\; (k=0,\pm1,\pm2,\cdots)$。解得：

利用极点留数法求解式(37)的逆变换，对于一阶极点$s = 0$，有：

对于普通极点${s_k}$，有：

解得第二类解析解：

4.2   数值逆变换和解析式对比

图9为第二类解的数值计算结果和解析解的结果（k_τ=1，$\bar x = 1.5$），两者的整体变化趋势一致，数值解存在一定数值振荡，所需的计算时间较长，相较而言，解析解是一种较优的计算方式。对比图3和图9，可以发现，第二类温度波和第一类温度波的单位波形不同，其差异主要表现在AB段和CD段。第一类解的上升和下降段均为线性，对应梯形波形单元。第二类解的AB段为曲线上升段，伴随斜率的减小，CD段为曲线下降段，单位波形与梯度结构相关。出现这一差异主要是因为第一类解和第二类解对应的材料梯度结构不同。

图  9  拉普拉斯逆变换数值计算结果和第二类解析结果对比

Figure  9.  Comparison of numerical results of inverse Laplace transformation and the results of the second analytic solution

下载: 全尺寸图片 幻灯片

图10为第二类解的时间空间二维分布，其分布表现出波动和传导双重特性。随着时间的发展，波动衰减，温度根据边界条件趋于平衡，第二类解的波动形式和幅值与第一类解不同。

图  10  第二类温度函数的二维分布

Figure  10.  Two-dimensional distribution of temperature function of the second analytic solution

下载: 全尺寸图片 幻灯片

图11显示了$\bar{x}$为1.3、1.5、1.8时对应的温度-时间曲线。可以发现，接近温度扰动源处的响应时间更短，波动幅度更小，但是衰减速度更快。当$\bar x = 1.3$时，AB段和BC段的温度变化幅度相当，而$\bar x = 1.8$时，温度变化以BC段为主，说明温度的振荡和衰减具有一定的偏向性。

图  11  第二类解在不同位置处的温度-时间曲线

Figure  11.  Temperature-time curves at different positions of the second analytic solution

下载: 全尺寸图片 幻灯片

图12为不同时刻的温度分布曲线。如图12(a)所示，初始阶段，温度从扰动侧向另一侧传播，温度幅值随时间的发展而增大。图12(b)为第二阶段的温度分布曲线。可以发现，随着温度波在两侧温度边界之间来回传播，温度分布具有向平衡态发展的趋势，温度间断界面的温度差值逐渐减小。温度曲线与梯度结构相关，分析图12(b)和图6(b)可以发现，温度曲线的凸性不同。

图  12  第二类解在不同时刻的温度分布曲线

Figure  12.  Temperature-position curves at different times of the second analytic solution

下载: 全尺寸图片 幻灯片

图13(a)为第二类解在不同弛豫系数下的温度-时间曲线。可以看出，随着弛豫时间系数的增大，曲线下降阶段CD趋向线性，振荡脉宽增加。图13(b)显示了第二类解的弛豫系数与响应时间、温度极值的关系。可以看出，随着弛豫时间系数的增大，温度响应时间增大，温度幅值增加，与图7(b)的趋势一致。

图  13  第二类解在不同弛豫系数下的温度时间曲线及响应时间、温度极值与弛豫系数的关系

Figure  13.  The temperature-time curves corresponding to different relaxation coefficients and the relationships between response time, maximum temperature and relaxation coefficient of the second analytic solution

下载: 全尺寸图片 幻灯片

5.   结　论

结合Cattaneo-Vernotte线性双曲型热传导方程，得到了幂律梯度材料的一维双曲型热传导方程；利用积分变换法和极点留数法，得到了控制方程在一般温度边界条件下的第一类解析解和第二类解析解；基于线性叠加原理，得到了温度脉冲载荷下的温度解，验证了数值结果和解析结果的可靠性；结合Mo-ZrC梯度材料的热力学性能，对特定梯度材料组合下温度场的时间和空间分布进行了分析，主要结论如下。

(1) 温度扰动以温度波的形式传播，表现出波动和传导双重特性，波动随时间逐渐衰减，最终趋于平衡态。在一般温度边界条件下，边界温度和梯度结构决定平衡态温度分布曲线的凸性。

(2) 第一类解对应的温度-时间曲线表现为随时间衰减的周期性梯形波动；第二类解对应的温度-时间曲线表现为随时间衰减的周期性异形波动，其波形单元存在曲线上升段和曲线下降段。两类解析解反映材料的梯度结构差异。

(3) 不同位置处的温度-时间曲线振荡具有偏向性，远离扰动源时，扰动响应时间增长。响应时间和温度极值随着弛豫时间系数的增大而增大，并趋向稳定值。单元波形随着弛豫系数的增大而变化。