第43卷第2期力学与实践2021 年4月
刘承斌*，*12)龚顺风*王惠明t
*(浙江大学土木工程学院，杭州310058)
t(浙江大学工程力学系，杭州310〇27)
摘要基于非傅里叶热传导理论，分别从频域和时域分析了球对称瞬态热传导问题。频域求解时，采用 L a p l a c e变换结合状态空间法以及层合近似模型，过程清晰统一。时域求解时采用了分离变量法和叠加法，能 减小数值反变换所带来的误差。数值算例比较了两种方法的计算效果，结果表明时域法可清晰捕捉热波传播过
程中的前沿波阵面。
关键词非傅里叶理论，球壳，瞬态热传导，频域法，时域法
中图分类号：0343.6, 0343.8, T B381文献标识码：A doi: 10.6052/1000-0879-20-396
ANALYSIS OF SPHERICALLY SYMMETRIC TRA NSIENT HEAT CO NDUCTION BASED ON NON-FOURIER THEORY1)
Liu C h e n g b i n*’2)G o n g Shunfeng*W a n g H u i m i n g卞
* (D ep artm en t of Civil E ngineering, Zhejiang University, H angzhou 310058, C hina)
"^(Departm ent of E ngineering M echanics, Zhejiang University, H angzhou 310027, C hina)
A b s t r a c t
B a s e d o n the non-Fourier heat conduction theory,the spherically s y m m e t r i c transient heat c o n d u ction p r o b l e m is analyzed in b o t h frequency d o m a i n a n d time d o m a i n.B y using the Laplace transform a n d the state-space formula in the frequency d o m a i n along with the laminated approximation m o d e l,the time responses are obtained numerically through the inverse Laplace transform.This appr o a c h is straightforward a n d with a unified f o r m.For the analysis in the time d o m a i n,a combination of the m e t h o d of separation of variables a n d the superposition m e t h o d is a d o p t e d,without requiring the numerical inverse Laplace transform.N
umerical results indicate that the analysis in the time d o m a i n can clearly capture the w a v e front of the thermal w a v e s. K e y w o r d s non-Fourier theory,spherical shell,transient heat conduction,frequency d o m a i n m e t h o d,time d o m a i n m e t h o d
球壳结构因其受力合理在工程中有广泛的应用, 大至储油罐，小至压电传感器等。当球壳结构工作于 温差突然有较大变化的环境中时，急剧的温度变化 使得结构内部产生瞬态热应力。通常将这种短促时 间内出现的热应力数值上非常大的现象称为热冲击, 热冲击有可能产生很大的拉压应力，甚至可能达到脆性材料的破坏应力，从而造成材料破坏。以往常规 热分析大多基于经典傅里叶热传导定律，这种经典 傅里叶热传导的描述可以满足绝大部分工程理论分 析的精度要求，但存在热传播速度无穷大这一与实 际物理规律不符合的弊端。随着研宄的不断深入和 应用范围的不断扩大，特别是在极端环境下服役时，
2020-09-15收到第1稿，2021-11-16收到修改稿。
1)国家自然科学基金（11772296)和浙江省公益技术应用研究项目（L G G21E030007)资助。
2)刘承斌，闻级工程师，从事压电弹性力学方面的研究工作。E-m ail: lcb@zju.edu
引用格式：刘承斌，龚顺风，王惠明.球对称瞬态非傅里叶热传导分析.力学与实践，2021, 43(2): 212-218
Liu C hengbin, G ong Shunfeng, W ang Huim ing. Analysis of spherically sy m m etric tra n sie n t heat conduction based on non-Fourier theory. M echanics in Engineering,2021, 43(2): 212-218
第2期刘承斌等：球对称瞬态非傅里叶热传导分析213
经典理论己不能满足分析要求W。为克服经典傅里叶
热传导定律的局限性，Cattaneo[2l和Vernotte丨3丨给出了能够描述热以有限速度传播的双曲型热传导方
程，被称为C-V方程或非傅里叶热传导方程。这种 广义热传导问题可采用积分变换、差分法、有限元以 及解析法等多种方法进行求解。积分变换方法将时 间域的偏微分方程变换为频域的常微分方程，再进 行求解，然后进行积分反变换，求得时域响应。Jiang 等W利用拉普拉斯变换，针对某些边界条件，得到了 空心球壳受热冲击的解析解。B a b a e i等W利用拉普 拉斯变换和常规数值反变换得到功能梯度空心球壳 的广义热传导响应。C o n k e r等[6]将拉普拉斯变换和 修正的数值反变换方法相结合，探讨了材料特性按 指数变化的功能梯度厚壁圆柱壳和球壳双曲型非傅 里叶热传导问题。W a n g等W利用积分变换研究了考 虑界面裂纹的层状复合材料的非傅里叶热传导现象。有限差分方法根据不同的边界条件和初始条件选择 合适的差分格式，将问题离散化后求解，所采用的差 分格式决定了求解精度，且必须满足收敛性和稳定
性要求M。广义热弹性有限元方法，可避免数值积分 反变换造成的截断误差，直接在时间域求解，但推导 相对复杂A11】。分离变量法I12-15】在球壳表面受热冲 击作用的热传导分析中也非常有效。此外，R a h i d e h 等1161利用微分求积法研究了叠层功能梯度平面域非
傅里叶热传导问题。张浙等对非傅里叶热传导的研宄进展做了评述，从实质、模型的建立和求解及应 用于实验等多个方面进行了论述，并指出了今后的 研宄方向。蒋方明等1181进一步介绍了非傅里叶导热 的研究进展。近年来，M〇〇saie[19l研宄了空心球的球 对称瞬态热传导问题，为便于求解，将解分为两部分 之和：一部分稳态解对应非齐次边界条件，另一部分 瞬态解对应齐次边界条件。Z h a o等[2C)!通过分离变 量法结合杜哈梅原理，分析了实心球外受任意表面 热冲击条件下的非傅里叶导热响应。许光映等[21]利 用积分变换方法探讨了非傅里叶导热边界条件对激 光辐射生物组织热传导的影响。刘承斌等[22]从三维 压电弹性理论出发，利用拉普拉斯变换结合状态空 间法，分析了压电球壳受表面球对称热冲击作用下 的热弹性问题。张彦博等[231基于双曲型单相延迟非 傅里叶热传导方程，采用有限元法研宄了含裂纹厚
壁圆筒结构在热冲击载荷下的热力学数值解。陈豪 龙等[241提出微分转换双重互易边界元法求解功能梯度材料非傅里叶瞬态热传导问题。结果表明时间 步长对计算结果的影响较小，该方法可以有效求解 功能梯度材料非傅里叶瞬态热传导问题。郭松林
基于双相延迟非傅里叶热传导模型对工程中常见的 圆柱、板、涂层以及夹芯板等结构在热冲击载荷作用 下的裂纹问题进行了分析，系统分析了结构热弹性 响应和裂纹尖端应力强度因子。
综上所述，尽管广义热传导问题己有很多有效 的求解方法，但随着新材料技术的不断发展和服役 环境的更加严苛，需要寻求计算效率高且准确的方 法。本文针对球对称瞬态非傅里叶热传导问题，给出 了两种不同的求解方法，并通过具体算例分析比较 了两种方法的计算效果。
1球壳瞬态热传导基本方程
考虑内外半径分别为a和6的空心球壳受球对 称热冲击作用丨图1)，球对称情形的L-S型广义热 传导方程为
^33
d22d
.dr2
+
r d r,
T=pcv I
dT d2T、
m+T~d^.
⑴
式中，为径向热传导系数，7■为热松弛时间，r为 温度变化值，p为质量密度，&为比热容。由非傅里 叶定理可知
f c33V2T= -(^0r+T0rj(2)式中，V2 =r3/5r，6>r =rgr，qr为径向热流密度，符 号上面的.表示对时间i求导。
◎
图1球坐标及球壳示意图
存在三种不同的温度边界条件：给定边界处的 温度、热流密度或者两者的线性组合。为便于推导，本文考虑温度边界条件，其他边界条件采用同样的 步骤可以类推。球壳内外表面的温度边界条件和初 始条件可表示为
工业反哺农业T(r,t)\r=a =i p〇H(t),T(r,t)\r=b=(3
)
214
力学与实践2021年第43卷
=〇
⑷
t =0
式中，(A )和^1为某特定值，丑⑴为阶跃函数。2频域求解
6 (r , 0) = 0, k33
86 (r , t )
d t
r = a + (2i - 1)0 — a )/2p ，由矩阵理论得式（10)的
解为
V
(r ) = e x p (r — rj _i ) M
(!)j V (rj _i )
(r ,-! ^ r < z = 1.2,
,p )
(11)
2.1
状态方程的建立
将式（2)代入式（1)，整理后得到
相邻子层间的连续条件要求状态变量连续，由式（11) 出发递推可得
Q Q V
26>r +
+ @r -h t
a t
2 f d T d T \
-r p c v
+
d O r ~df
(5)
张武力为消除时间变量，将偏微分方程组转化为常微分方 程组，引入拉普拉斯变换
F (s ) = ^[f (t )}= r e
~s t f (t )dt ,
J o
R e (s ) >0, s = a -\-i c j
(6)
将式（2)和式（5)利用式（6)变换后，联合得到状态 方程
▽2
T
1 + TS
^33
(7)
V (b ) = T V {a )
(12)
其中
T
H
exp
(6-a )
M (
i )为二阶方阵。引入球壳内外温度边界条件式（8)，代
入式（12)，得从而可求得球壳内表面的热流密度状态量
〇r {〇)
(f l — (f 〇Tn
T \2S
(13)
(14)
球壳边界条件式（3)经过拉普拉斯变换后变为
T \r
S
’
T |r =b
s
2.2
状态方程的求解
为求解方便，对变量进行无量纲化
(8)
(9)
其中c 33为材料的弹性参数，c 为弹性波速，r /二
p c …/f c 33为热黏度系数。为简便考虑，省略后续公式
符号中的上标符号将式（9)代入式（7)，可得
t
Y :
=M V
式中V "=[T .〇r ]0
—
(1 + T s )-1
〇
—
r z s
V
(10)
T 。由于矩阵M
中的元素含有与
半径r ■相关的量，可采用层合模型将其转化为常系 数矩阵，当分层的模型足够细时，这种近似方法有 足够的精度。将球壳沿厚度方向等分为P 层，矩 阵中的坐标变量r 取为子层中间厚度处坐标值，即
式中为矩阵了的元素。利用式（11)可以求得沿 径向任意位置的状态量温度值。
对于大部分边界条件，无法直接将频域解解析 为时域解，必须采用近似方法进行拉普拉斯数值反 变换丨26_29)。本文采用D u r b i n ^26丨的傅里叶级数法
/⑴=
-全 R e (/(o 〇) +
N U S M
k =0
R e
/ o ； + i
kn
~U
I m
(15)
其中[/为反变换的周期，合理的参数取值为5彡
a U < 10, 50 ^ N U S M ^ 5000〇
3时域求解 3.1分解
将式（9)代入式（5)化简后为
(^ + l t )T =
(l +r^)T
(16)
第2期刘承斌等：球对称瞬态非傅里叶热传导分析215
由于边界条件式（3)非齐次，不能直接利用分离变量 法来求解方程（16)，因而需要进行边界条件的齐次 化，将温度场解分解成两部分：对应非齐次边界条件 的解r s (r )和对应齐次边界条件的解T D (r ，$
即T (r ,t ) = Ts (r ) + TD (r ,t )
(17)
3.2稳态解r s (r )
非齐次边界条件下的热传导边值问题可以表示为
d2r s (r)
2dT s (r)
dr2 r
dr
Ts (r )
式（18)的解为
^>o , Ts (r )
<P i
Ts (r ) = A
〇 +
B 〇
(18)
(19)
(20)
其中砒和说为待定常数，将式（20)代入边界条 件式（19)，得到
Ts (r ) =
-------1-----^
(21)
b — a
r b — a
3.3瞬态解(r )
满足齐次边界条件及初始条件的球对称热传导 问题为
+
2-l )TD (r ,t )=r d2T^ '
\d r 2
r dr J u 、’’ dt 2
dTD (r, t )
dt t d (r, t)|r=a = 0, T d (r, t )\r =b = 0
(22)
(23)T d (r , 0) ^ i
,dTD (r ^)U ( f c 33^^ = 0
a ^p 〇 —
b <f \ a b ((p 〇 — (fi ) a r (b — a )假设r D (r ,《）=i ?T (r )I /⑴，将其代入式（22)，从而 得到
1 (d2R T 2d R T 、
r \ dr2 r dr y
1 ( d2L d L 、
—L J
(25)
考虑式（25)两端只能等于常数-a ;2，分解成两个2 阶常微分方程
d2凡! d r2 2 c I j R t 2 r >H --------h u 2 R t = 0r dr d2L dL 2r n r d ^+d ^+w L = 0
(26)
(27)
方白羽智枭相应的边界条件为
W |r=a = 〇, -R
t
(r -)|r=6 = 〇 (28)
式（26)为贝塞尔方程，其解析解为
R t (r )
C 〇\ — sin (ur ) +
熏洗仪、’ 7TW
D 〇\l
cos (cjr)
n u j
(29)
其中C 。和
A
)为待定系数。将式（29)代入边界条
件式（28)，可以得到
C 〇 sin (〇ja ) +
D 〇 cos (c o a ) — 0 C 〇 sin (u j b ) + D
〇 cos (t u b ) = 0
(3〇)
为保证式（30)存在非零解，系数矩阵行列式必须为 零，即
sin (c o a ) cos (c u b ) — cos (u i a ) sin (c o b ) = 0
(31)
求解上述超越方程，可以得到无穷个特征根(f c = 1，2,…），对于任意叫，式（29)的解为
(24)得到
Rrk (r) = \-------[cos (wf e 6) sin (uikr )-V r
sin (cjfc6) cos (cjfcr)] (/c = l ,2, •••)
(32)
二阶常微分方程（27)的解为
L k (t) = e ~^ (E 〇k cos K ,kt + F 〇k sin nkt) (33)
其中 M - l /2r ，联合式（32)和式（33)
T d (r, t) = e_</(2r) (E 〇k cosK,kt-\-k =l
F 〇k sin K ,kt) R Tk (r )
由初始条件式（24)第1式得
智血
(34)
T d
(r , 0) = ^ E 〇kR Tk (r )=
k =l
a (f 〇 -
b <f i — ab ((p 〇 - (p i )
b — a (35)r (b — a)
对应不同本征值的本征函数(r )和(r )正交
性证明推导过程及模值详见文献[30]，这里直接给出
模值
216力学与实践2021年第43卷
求解时，频域法采用层合近似法将球壳沿厚度方向 细分60等分，时域法采用前70阶特征根，和文献 [31]得到的球壳内表面温度随时间变化的结果吻合 得非常好，如图2所示，表明本文理论推导和数值
计算的正确性。
t = 0.4时域法
t = 0.4 文献[31]< =0.4频域法
1.0
0.8 -
Q l = [ r X l d r =
[ R r kr2dr
2^f
r 3
(^) +
r 2R TkdRTk
dr
dr
^lr：i R\k
利用级数正交展幵技术，由式（35)可得
(36)
E 〇k = —
I r2
Q l J a Rrk (r ) dr
a (p 〇 —烟酒伴侣
b (f \ ab ((f 〇 — ^p \) b — a
r (b — a )
(37)
由初始条件式（24)第2式得
s r D (r ,〇)
^i E 〇k + K k F \RTk (r ) = 〇 (38)
d t
k =l
求得系数F o f c
2rnk
从而得到瞬态温度场
TD (r ,t )
=e ~^J 2
k =1 Ql V 2^f c
r b
sin Hkt +
c o s Kkt j x R k (r ) J r2R k (r ) x
a (f 〇 -
b (f i
a 6 (p
〇 —
d
厂
(39)
b — a r (b — a )
3.4全解
将稳态解和瞬态解两部分叠加，得到总的温度场
T (r , t ) = Ts (r ) + T D (r , t )=
b i p i - atpo + 1 ab (i p 〇 - i p x) + b — a
r
b — a
e ^E
Ql \2TKk
sin tikt + cos ] x
R k {r ) j r2R k (r ) x
a(p 〇 - 6(^i ab (ipQ - ifO
dr
(40)
b — a
r (b — a )
4算例
为验证本文解的正确性和有效性，将两种方法 与文献[31]计算结果进行对比。球壳无量纲参数为： 内半径a = 1，外半径6 = 2，波速ci = 1，C 2二 0.535, e = 0.02。温度边界条件为
T {rA)\r
1 - (1 + 100t )e -loot
dT
0 (41)
0.2 -
0 - 1.0
1.2
1.4
1.6
1.8
2.0
r 图2球壳内表面温度随时间变化
下面以晒化镉球壳受内外温差变化为例，材料
参数为：c 33 = 83.6 G P a , /9 = 5.68 t /m 3，A ：33 =
12.9 W /(m .K )，= 0.46 x 103 J .k g /K 。无量纲 热弛豫时间t = 0.1，球壳内外径之比a /6取为 0.5，内/外表面分别同时突加温升2f )0°C /400°C 。
图3为采用频域法所求得对应不同时刻内部温度变 化沿着球壳径向分布情况，计算时，将球壳沿厚度方 向80等分。观察后可以发现，内外表面热冲击所产 生的热波，以无量纲速度向球壳对向传播。在
温度 波前尚未到达的区域内，温度保持初始温度不变，当 热波波前到达某位置，该处形成一个明显的温度阶 跃现象，当双向的热波交汇后，球体内部的温度会超 过边界处设定的温度。产生的原因在于热传导时间 与热驰豫时间相当时，材料的局部热平衡己经不能 维持W 。
U
图3对应不同时间f 温度沿球壳厚度方向变化（
频域法）
豫公网安备41160202000603
站长QQ:729038198
关于我们
投诉建议