赵伟涛 吴九汇†

1)(西安交通大学机械学院，西安 710049)

2)(西安交通大学，机械结构强度与振动国家重点实验室，西安 710049)

(2013年4月2日收到;2013年5月4日收到修改稿)

1 引言

自从17世纪Fourier建立了导热的数学模型，Fourier定律随之被广泛应用于导热问题分析的各个领域.对于热作用时间较长的稳态传热过程以及热传播速度较快的非稳态常规传热过程，采用Fourier定律来描述热流密度与温度梯度之间的关系是可以满足精度要求的.但是，Fourier定律不涉及传热时间项，隐含了热扰动传播速度为无限大的假设.随着科学技术的进步，超短激光脉冲的出现和制冷水平的提高，存在着极高(低)温条件下的传热问题和超急速传热问题，使得Fourier定律中的准平衡条件假设不再成立[1-4].

为了克服Fourier定律的局限性，Cattaneo[5]和Vernotte[6]分别独立地提出了具有热流延迟相的非Fourier热传导模型，计及热流变化率对热传导的影响，其修正后的热传导方程为

其中，q是热流矢量，T为温度，k为热传导率，t是时间，τ0为热松弛时间.方程(1)结合能量守恒方程得到以温度T描述的双曲线型热传导方程：

其中，a=k/(ρc)为热扩散系数，ρ为密度，c为常应变比热.

为描述热以有限速度传播这一超常规热传导现象，研究者们对方程(2)在不同初始条件和边界条件下进行了求解.Tao等[7]采用数值模拟研究了固体激光器的非Fourier导热问题.李世荣等[8]采用单相延迟模型，研究了薄板在受周期热流边界条件下板内的温度响应.Sarkar和Haji-Sheikh[9]采用Laplace变换技术研究了由介电材料制成的有限平板的双曲型热传导问题.Tang和Araki[10]求解了表面周期性加热条件下有限介质的非Fourier热传导问题.Moosaie[11]在Tang等的工作基础上，运用Fourier积分表达式，求解了有限介质在周期表面热扰动条件下的非Fourier热传导问题.Barletta和Zanchini[12]分析了无限圆柱体存在内热源以及与外界流体有热对流的情况下的双曲型热传导问题.Ate fi和Talaee[13]使用分离变量法求解了边界条件不随时间变化的无限长圆筒的非Fourier温度场.Mishra和Sahai[14]运用格子Boltzmann法研究了一维圆柱和球体的非Fourier热传导问题.Jiang[15]运用Laplace变换法研究了空心球体在内外两个表面温度突然变化时的双曲型热传导问题.Shirmohammadi和Moosaie[16]采用分离变量法得到空心球体在周期表面热流条件下的解析解.

本文主要运用Duhamel积分和Fourier级数展开法得到平板前表面热流任意周期变化时双曲型热传导方程的解析解.首先采用Duhamel积分和分离变量法分析了方程在平板前表面遭受突变热流和简谐热流这两类特殊情况下解的形式，在此基础上应用Fourier级数展开法和叠加原理研究了平板前表面热流任意周期变化时解的形式.按照这些表达式，不同周期边界条件下平板的双曲线热传导行为得到分析和研究.这为工程应用和数值计算的验证提供了便利.

2 数学模型

考虑一厚度为L，边界绝缘，热物性为常数的平板，假设其初始温度T(x,0)=T0，从时间t=0时起，平板一表面x=0处(称该表面为平板前表面，另一面称平板后表面)遭受一热流为q0·q(t)的作用.在这种情况下，由方程(2)可得一维平板的非Fourier热传导方程：

为了获得方程(3)—(5)的无量纲形式，特引入以下无量纲量：

式中，Fo是 Fourier数，Ve是Vernotte数，1/Ve代表温度波传播的无量纲速度.

则方程(3)—(5)无量纲化之后为

3 分析求解

由于边界条件(8)式中的f(Fo)是任意函数，这就使得直接求解方程(7)变得不可能.因此，首先假定f(Fo)为一时间无关量f，求解此条件下的温度场;其次运用Duhamel积分，求解在任意函数f(Fo)下的温度场.

设在边界条件f下方程(7)的解为

当Ve→0时，方程(7)对应经典Fourier导热的情形.此时βn为实数，(35)式变为

把(37)式代入(34)式可得经典Fourier导热时平板内部的温度场.

4 结果分析与讨论

为了得到任意周期边界条件下平板温度场的解析解，将做如下递进式推导.4.1给出了边界条件为常数的情况;4.2的边界条件为简谐周期函数;4.3在前两步的基础上，计算与讨论了边界条件可用Fourier级数展开的任意周期函数的情况.

4.1 突变边界条件

当厚度为L的平板前表面遭受大小为q0的突变热扰动(q(t)=1)时，无量纲化后

此时(30)式就是方程(7)在边界条件为常数下的解

对于Fourier导热，当Ve→0时，方程的解为

根据(39)式，可以计算得到厚度为L的平板前面遭受一突变热扰动时的温度分布.为了验证结果的可靠性，直接对方程(7)—(9)进行数值求解.图1给出了平板前表面遭受突变热扰动时，在Ve=0.8下，平板后表面温度响应的数值解与解析解的对比.可以看出，数值解和解析解的误差非常小，几乎看不到差异，从而验证了解析解的正确性.

图1 平板前表面遭受突变热扰动时后表面温度响应的数值解与解析解的比较

图2 平板前表面遭受突变热扰动时不同热松弛时间Ve下的表面温度响应 (a)前表面温度响应;(b)后表面温度响应

图2给出了平板前表面遭受突变热扰动时不同无量纲热松弛时间Ve下，平板两表面无量纲温度随无量纲时间的变化曲线.可以看出当Ve=0.3，0.6，0.8时，其温度分布曲线与典型的Fourier温度分布曲线已不一致，尤其是当Ve=0.6，0.8时，其曲线有了凸凹点，热传导的非Fourier效应变得更加明显.从图2(b)可以看出其温度响应曲线均存在一明显的延迟，由于热波的传播速度为1/Ve，所以在平板后表面(无量纲化后X=1)，其延迟时间为Fo0=Ve.

4.2 简谐变化边界条件

当厚度为L的平板前表面遭受一幅值为q0，角频率为α的简谐热扰动(q(t)=cos(αt))时

式中，Fo1=a/αL2.

把(41)式代入(34)和(35)式得平板前表面热扰动简谐变化时平板内部的温度场分布

式中

对于Fourier导热，当Ve→0时，把(41)式代入(34)和(37)式可得

使用(42)式，可以模拟得到厚度为L的平板在其前表面受到简谐热扰动时，平板前后两表面的温度响应.

图3给出了平板前表面在简谐热流作用下，当Fo1=0.25，Ve=0.8时，平板前后两表面的温度响应曲线，反映出了非Fourier热传导的瞬态温度特性.从图3可以看出两条曲线都存在一系列阶跃点，它们表示温度波波前经过边界的传播和反射到达的时刻.由于热波的传播速度为1/Ve，所以在两表面X(X=0,1)处，阶跃点位置为Fo=Ve·X，Ve(2+X)，Ve(4+X),···.由于扩散的作用，前表面温度波的阶跃值和幅值要比后表面的大.从图3还可以看出温度波曲线的阶跃点随着时间的增加而减小并逐渐消失，曲线最终变成光滑的周期曲线.

图3 平板前表面遭受简谐热扰动时前后两表面的温度响应(Fo1=0.25，Ve=0.8)

图4 平板前表面遭受简谐热扰动时不同热松弛时间Ve下前后两表面的温度响应(Fo1=0.25) (a)前表面温度响应;(b)后表面温度响应

图4给出了在前表面遭受简谐热扰动和Fo1=0.25的情况下，不同无量纲热松弛时间Ve下，平板前后两表面处无量纲温度随无量纲时间的变化曲线，以及采用(44)式得到的Fourier热传导模型下的温度响应曲线.从图中可以看出，随着热松弛时间Ve的减小，热波的传播速度增大，热量能够迅速地传递，从而导致平板内温度响应的幅值也随之减小，非Fourier效应逐渐减弱.Ve越小，非Fourier和Fourier的温度响应曲线越接近.

4.3 任意周期边界条件

当厚度为L的平板前表面遭受任意周期热流变化时，假设该周期函数q0·q(t)的周期为2l.在实际问题中，q(t)只能在t≥0上有定义，故我们可以在t＜0的区间内将函数q(t)进行偶延拓，使q(t)=q(-t)，即延拓后的函数为偶函数.这样周期函数q(t)可以展开成余弦级数的形式

当i=0时，方程(7)的解可以由(39)式得到：

把(50)与(52)式进行叠加，可得平板表面温度任意周期变化时的非Fourier温度场：

当Ve→0时，根据叠加原理与(40)，(44)和(45)式，可得在任意周期边界条件下的非Fourier导热温度场：

4.4 应用举例

设q(t)是周期为2π的三角波函数(如图5所示)，它在[0,2π)上的表达式为

将函数q(t)进行偶延拓后，展开成余弦级数

将(57)式分别代入(53)式和(54)式得三角波条件下的非Fourier和Fourier导热温度场.

图5 三角波热扰动

图6给出了在前表面遭受三角波热扰动和a/L2=0.25的情况下，不同无量纲热松弛时间Ve下，平板前后两表面处无量纲温度随无量纲时间的变化曲线.从图6可以看出，该温度响应曲线除具有前面介绍的规律外，还具有随着Ve的增大，特别是Ve=0.50时，温度响应曲线形状越接近前表面输入的三角波热扰动信号这样的现象.这是因为，Ve越大，热波的传播速度越小，其延迟效应越明显，随着时间的增加，当阶跃消失时，其温度响应曲线形状将越接近输入信号.

图6 平板前表面遭受三角波热扰动时不同热松弛时间Ve下前后两表面的温度响应(a/L2=0.25) (a)前表面温度响应;(b)后表面温度响应

5 结论

本文首先通过分离变量法和Duhamel积分原理，对平板前表面遭受突变热扰动和简谐热扰动这两类边界条件下的非Fourier热传导问题进行了分析求解，在此基础上，利用Fourier级数展开法和叠加原理，得到了双曲线型热传导方程在平板前表面遭受任意周期变化热流时这个最一般情况下的解析解.通过理论计算与数值模拟，展示了非Fourier热传导模型所给出的温度响应与Fourier热传导模型的如下差别：

1)温度响应曲线存在一系列的阶跃点，在两表面X(X=0,1)处，它们依次出现的时间为Fo=Ve·X，Ve(2+X)，Ve(4+X)，···，但随着时间的延长，温度响应曲线逐渐变得光滑;

2)同一位置处，温度响应的幅值随着热松弛时间Ve的减小而减小，最终趋近于Fourier热传导时的温度响应幅值;

3)具有有限的传播速度，其大小与1/Ve成线性关系;

4)热松弛时间Ve越大，温度响应曲线的形状越和前表面遭受的热流形状相似.

[1]Wang Y Z，Song X N 2012Acta Phys.Sin.61 234601(in Chinese)[王颖泽，宋新南2012物理学报61 234601]

[2]Guo Z Y，Cao B Y 2008Acta Phys.Sin.57 4273(in Chinese)[过增元，曹炳阳2008物理学报57 4273]

[3]Tian X G，Shen Y P 2012Adv.Mech.42 18(in Chinese)[田晓耕，沈亚鹏2012力学进展42 18]

[4]Tung T L，Fong E 2011Int.J.Heat Mass Transfer54 4796

[5]Cattaneo C 1948Atti.Sem.Mat.Fis.Univ.Modena3 83

[6]Vernotte P 1958C.R.Acad.Sci.246 3154

[7]Tao Y J，Huai X L，Li Z G 2006Chin.Phys.Lett.23 2487

[8]Li S R，Zhou F X，Wu H M 2007Engineer.Mech.24 48(in Chinese)[李世荣，周凤玺，吴红梅2007工程力学24 48]

[9]Sarkar D，Haji-Sheikh A 2012Int.Commun.Heat Mass Transfer39 1009

[10]Tang D W，Araki N 1996Int.J.Heat Mass Transfer39 1585

[11]Moosaie A 2007Int.Commun.Heat Mass Transfer34 996

[12]Barletta A，Zanchini E 1997Heat and Mass Transfer32 285

[13]Ate fiG，Talaee M R 2011Arch.Appl.Mech.81 569

[14]Mishra S C，Sahai H 2012Int.J.Heat Mass Transfer55 7015

[15]Jiang F M 2006Heat and Mass Transfer42 1083

[16]Shirmohammadi R，Moosaie A 2009Int.Commun.Heat Mass Transfer36 827