铁球在土壤中的散热1
拉普拉斯变换法
文章目录
1. 问题描述
考虑一个半径为 R R R 的铁球，初始温度为 T i T_i Ti，放置在初始温度为 T 0 T_0 T0 的均匀土壤中。铁球的热物性参数为密度 ρ s \rho_s ρs、比热容 c s c_s cs；土壤的热物性参数为导热系数 k k k、热扩散率 α \alpha α。如下图：
由于铁球导热性能良好，可视为集总系统；土壤被视为半无限大介质，热传导满足球对称条件。
注：以前采用分离变量法，发现最后会走进死胡同，无法求解。改为拉普拉斯变换法后，仍然无法直接得到理论解，但是可以得到比较实用的级数解。
2. 数学模型
2.1 控制方程
土壤区域的瞬态热传导方程在球坐标系下表示为：
∂ T ( r , t ) ∂ t = α ( ∂ 2 T ∂ r 2 + 2 r ∂ T ∂ r ) , r ≥ R , t > 0 (1) \frac{\partial T(r,t)}{\partial t} = \alpha \left( \frac{\partial^2 T}{\partial r^2} + \frac{2}{r} \frac{\partial T}{\partial r} \right), \quad r \geq R, t > 0 \tag{1} ∂t∂T(r,t)=α(∂r2∂2T+r2∂r∂T),r≥R,t>0(1)
铁球的能量守恒方程为：
ρ s c s 4 3 π R 3 d T s ( t ) d t = 4 π R 2 k ∂ T ∂ r ∣ r = R (2) \rho_s c_s \frac{4}{3}\pi R^3 \frac{dT_s(t)}{dt} = 4\pi R^2 k \left. \frac{\partial T}{\partial r} \right|_{r=R} \tag{2} ρscs34πR3dtdTs(t)=4πR2k∂r∂T r=R(2)
2.2 初始条件和边界条件
初始条件：
T
(
r
,
0
)
=
T
0
,
T
s
(
0
)
=
T
i
(3)
T(r,0) = T_0, \quad T_s(0) = T_i \tag{3}
T(r,0)=T0,Ts(0)=Ti(3)
边界条件：
T
(
R
,
t
)
=
T
s
(
t
)
,
lim
r
→
∞
T
(
r
,
t
)
=
T
0
(4)
T(R,t) = T_s(t), \quad \lim_{r\to\infty} T(r,t) = T_0 \tag{4}
T(R,t)=Ts(t),r→∞limT(r,t)=T0(4)
3. 拉普拉斯变换求解
3.1 变量变换
引入过余温度：
θ
(
r
,
t
)
=
T
(
r
,
t
)
−
T
0
,
θ
s
(
t
)
=
T
s
(
t
)
−
T
0
(5)
\theta(r,t) = T(r,t) - T_0, \quad \theta_s(t) = T_s(t) - T_0 \tag{5}
θ(r,t)=T(r,t)−T0,θs(t)=Ts(t)−T0(5)
问题变为：
∂
θ
∂
t
=
α
(
∂
2
θ
∂
r
2
+
2
r
∂
θ
∂
r
)
(6)
\frac{\partial \theta}{\partial t} = \alpha \left( \frac{\partial^2 \theta}{\partial r^2} + \frac{2}{r} \frac{\partial \theta}{\partial r} \right) \tag{6}
∂t∂θ=α(∂r2∂2θ+r2∂r∂θ)(6)
边界条件：
θ
(
R
,
t
)
=
θ
s
(
t
)
,
θ
(
∞
,
t
)
=
0
(7)
\theta(R,t) = \theta_s(t), \quad \theta(\infty,t) = 0 \tag{7}
θ(R,t)=θs(t),θ(∞,t)=0(7)
初始条件：
θ
(
r
,
0
)
=
0
,
θ
s
(
0
)
=
T
i
−
T
0
=
θ
i
(8)
\theta(r,0) = 0, \quad \theta_s(0) = T_i - T_0 = \theta_i \tag{8}
θ(r,0)=0,θs(0)=Ti−T0=θi(8)
铁球能量方程：
d
θ
s
d
t
=
3
k
ρ
s
c
s
R
∂
θ
∂
r
∣
r
=
R
(9)
\frac{d\theta_s}{dt} = \frac{3k}{\rho_s c_s R} \left. \frac{\partial \theta}{\partial r} \right|_{r=R} \tag{9}
dtdθs=ρscsR3k∂r∂θ
r=R(9)
3.2 拉普拉斯变换
定义拉普拉斯变换：
Θ
(
r
,
s
)
=
L
[
θ
(
r
,
t
)
]
=
∫
0
∞
e
−
s
t
θ
(
r
,
t
)
d
t
(10)
\Theta(r,s) = \mathcal{L}[\theta(r,t)] = \int_0^\infty e^{-st} \theta(r,t) dt \tag{10}
Θ(r,s)=L[θ(r,t)]=∫0∞e−stθ(r,t)dt(10)
Θ s ( s ) = L [ θ s ( t ) ] = ∫ 0 ∞ e − s t θ s ( t ) d t (11) \Theta_s(s) = \mathcal{L}[\theta_s(t)] = \int_0^\infty e^{-st} \theta_s(t) dt \tag{11} Θs(s)=L[θs(t)]=∫0∞e−stθs(t)dt(11)
对控制方程 (6) 进行拉普拉斯变换：
L
[
∂
θ
∂
t
]
=
s
Θ
(
r
,
s
)
−
θ
(
r
,
0
)
=
s
Θ
(
r
,
s
)
(12)
\mathcal{L}\left[\frac{\partial \theta}{\partial t}\right] = s\Theta(r,s) - \theta(r,0) = s\Theta(r,s) \tag{12}
L[∂t∂θ]=sΘ(r,s)−θ(r,0)=sΘ(r,s)(12)
L [ α ( ∂ 2 θ ∂ r 2 + 2 r ∂ θ ∂ r ) ] = α ( d 2 Θ d r 2 + 2 r d Θ d r ) (13) \mathcal{L}\left[\alpha \left( \frac{\partial^2 \theta}{\partial r^2} + \frac{2}{r} \frac{\partial \theta}{\partial r} \right)\right] = \alpha \left( \frac{d^2 \Theta}{dr^2} + \frac{2}{r} \frac{d\Theta}{dr} \right) \tag{13} L[α(∂r2∂2θ+r2∂r∂θ)]=α(dr2d2Θ+r2drdΘ)(13)
因此，变换后的方程为：
s
Θ
=
α
(
d
2
Θ
d
r
2
+
2
r
d
Θ
d
r
)
(14)
s\Theta = \alpha \left( \frac{d^2 \Theta}{dr^2} + \frac{2}{r} \frac{d\Theta}{dr} \right) \tag{14}
sΘ=α(dr2d2Θ+r2drdΘ)(14)
整理得：
d
2
Θ
d
r
2
+
2
r
d
Θ
d
r
−
s
α
Θ
=
0
(15)
\frac{d^2 \Theta}{dr^2} + \frac{2}{r} \frac{d\Theta}{dr} - \frac{s}{\alpha} \Theta = 0 \tag{15}
dr2d2Θ+r2drdΘ−αsΘ=0(15)
3.3 求解变换后的方程
方程 (15) 是球贝塞尔方程。令 Θ ( r , s ) = U ( r , s ) r \Theta(r,s) = \frac{U(r,s)}{r} Θ(r,s)=rU(r,s)，则：
d Θ d r = 1 r d U d r − U r 2 (16) \frac{d\Theta}{dr} = \frac{1}{r} \frac{dU}{dr} - \frac{U}{r^2} \tag{16} drdΘ=r1drdU−r2U(16)
d 2 Θ d r 2 = 1 r d 2 U d r 2 − 2 r 2 d U d r + 2 U r 3 (17) \frac{d^2 \Theta}{dr^2} = \frac{1}{r} \frac{d^2 U}{dr^2} - \frac{2}{r^2} \frac{dU}{dr} + \frac{2U}{r^3} \tag{17} dr2d2Θ=r1dr2d2U−r22drdU+r32U(17)
代入方程 (15)：
( 1 r d 2 U d r 2 − 2 r 2 d U d r + 2 U r 3 ) + 2 r ( 1 r d U d r − U r 2 ) − s α U r = 0 (18) \left( \frac{1}{r} \frac{d^2 U}{dr^2} - \frac{2}{r^2} \frac{dU}{dr} + \frac{2U}{r^3} \right) + \frac{2}{r} \left( \frac{1}{r} \frac{dU}{dr} - \frac{U}{r^2} \right) - \frac{s}{\alpha} \frac{U}{r} = 0 \tag{18} (r1dr2d2U−r22drdU+r32U)+r2(r1drdU−r2U)−αsrU=0(18)
简化得：
1
r
d
2
U
d
r
2
−
s
α
U
r
=
0
(19)
\frac{1}{r} \frac{d^2 U}{dr^2} - \frac{s}{\alpha} \frac{U}{r} = 0 \tag{19}
r1dr2d2U−αsrU=0(19)
即：
d
2
U
d
r
2
−
s
α
U
=
0
(20)
\frac{d^2 U}{dr^2} - \frac{s}{\alpha} U = 0 \tag{20}
dr2d2U−αsU=0(20)
这个方程的通解为：
U
(
r
,
s
)
=
A
(
s
)
e
−
r
s
/
α
+
B
(
s
)
e
r
s
/
α
(21)
U(r,s) = A(s) e^{-r\sqrt{s/\alpha}} + B(s) e^{r\sqrt{s/\alpha}} \tag{21}
U(r,s)=A(s)e−rs/α+B(s)ers/α(21)
因此：
Θ
(
r
,
s
)
=
A
(
s
)
r
e
−
r
s
/
α
+
B
(
s
)
r
e
r
s
/
α
(22)
\Theta(r,s) = \frac{A(s)}{r} e^{-r\sqrt{s/\alpha}} + \frac{B(s)}{r} e^{r\sqrt{s/\alpha}} \tag{22}
Θ(r,s)=rA(s)e−rs/α+rB(s)ers/α(22)
由边界条件
Θ
(
∞
,
s
)
=
0
\Theta(\infty,s) = 0
Θ(∞,s)=0 得
B
(
s
)
=
0
B(s) = 0
B(s)=0，所以：
Θ
(
r
,
s
)
=
A
(
s
)
r
e
−
r
s
/
α
(23)
\Theta(r,s) = \frac{A(s)}{r} e^{-r\sqrt{s/\alpha}} \tag{23}
Θ(r,s)=rA(s)e−rs/α(23)
由边界条件
Θ
(
R
,
s
)
=
Θ
s
(
s
)
\Theta(R,s) = \Theta_s(s)
Θ(R,s)=Θs(s) 得：
Θ
s
(
s
)
=
A
(
s
)
R
e
−
R
s
/
α
⇒
A
(
s
)
=
R
Θ
s
(
s
)
e
R
s
/
α
(24)
\Theta_s(s) = \frac{A(s)}{R} e^{-R\sqrt{s/\alpha}} \quad \Rightarrow \quad A(s) = R\Theta_s(s) e^{R\sqrt{s/\alpha}} \tag{24}
Θs(s)=RA(s)e−Rs/α⇒A(s)=RΘs(s)eRs/α(24)
因此：
Θ
(
r
,
s
)
=
R
r
Θ
s
(
s
)
e
−
(
r
−
R
)
s
/
α
(25)
\Theta(r,s) = \frac{R}{r} \Theta_s(s) e^{-(r-R)\sqrt{s/\alpha}} \tag{25}
Θ(r,s)=rRΘs(s)e−(r−R)s/α(25)
3.4 耦合条件求解
计算在
r
=
R
r = R
r=R 处的温度梯度：
d
Θ
d
r
=
R
r
Θ
s
(
s
)
e
−
(
r
−
R
)
s
/
α
(
−
s
α
)
+
Θ
s
(
s
)
e
−
(
r
−
R
)
s
/
α
(
−
R
r
2
)
(26)
\frac{d\Theta}{dr} = \frac{R}{r} \Theta_s(s) e^{-(r-R)\sqrt{s/\alpha}} \left( -\sqrt{\frac{s}{\alpha}} \right) + \Theta_s(s) e^{-(r-R)\sqrt{s/\alpha}} \left( -\frac{R}{r^2} \right) \tag{26}
drdΘ=rRΘs(s)e−(r−R)s/α(−αs)+Θs(s)e−(r−R)s/α(−r2R)(26)
在
r
=
R
r = R
r=R 处：
d
Θ
d
r
∣
r
=
R
=
Θ
s
(
s
)
[
−
s
α
−
1
R
]
(27)
\left. \frac{d\Theta}{dr} \right|_{r=R} = \Theta_s(s) \left[ -\sqrt{\frac{s}{\alpha}} - \frac{1}{R} \right] \tag{27}
drdΘ
r=R=Θs(s)[−αs−R1](27)
铁球能量方程 (9) 的拉普拉斯变换为：
L
[
d
θ
s
d
t
]
=
s
Θ
s
(
s
)
−
θ
s
(
0
)
=
s
Θ
s
(
s
)
−
θ
i
(28)
\mathcal{L}\left[ \frac{d\theta_s}{dt} \right] = s\Theta_s(s) - \theta_s(0) = s\Theta_s(s) - \theta_i \tag{28}
L[dtdθs]=sΘs(s)−θs(0)=sΘs(s)−θi(28)
L [ 3 k ρ s c s R ∂ θ ∂ r ∣ r = R ] = 3 k ρ s c s R d Θ d r ∣ r = R (29) \mathcal{L}\left[ \frac{3k}{\rho_s c_s R} \left. \frac{\partial \theta}{\partial r} \right|_{r=R} \right] = \frac{3k}{\rho_s c_s R} \left. \frac{d\Theta}{dr} \right|_{r=R} \tag{29} L[ρscsR3k∂r∂θ r=R]=ρscsR3kdrdΘ r=R(29)
代入梯度表达式 (27)：
s
Θ
s
(
s
)
−
θ
i
=
3
k
ρ
s
c
s
R
Θ
s
(
s
)
[
−
s
α
−
1
R
]
(30)
s\Theta_s(s) - \theta_i = \frac{3k}{\rho_s c_s R} \Theta_s(s) \left[ -\sqrt{\frac{s}{\alpha}} - \frac{1}{R} \right] \tag{30}
sΘs(s)−θi=ρscsR3kΘs(s)[−αs−R1](30)
整理得：
s
Θ
s
(
s
)
−
θ
i
=
−
3
k
ρ
s
c
s
R
2
Θ
s
(
s
)
−
3
k
ρ
s
c
s
R
α
s
Θ
s
(
s
)
(31)
s\Theta_s(s) - \theta_i = -\frac{3k}{\rho_s c_s R^2} \Theta_s(s) - \frac{3k}{\rho_s c_s R\sqrt{\alpha}} \sqrt{s} \Theta_s(s) \tag{31}
sΘs(s)−θi=−ρscsR23kΘs(s)−ρscsRα3ksΘs(s)(31)
令：
a
=
3
k
ρ
s
c
s
R
2
,
b
=
3
k
ρ
s
c
s
R
α
(32)
a = \frac{3k}{\rho_s c_s R^2}, \quad b = \frac{3k}{\rho_s c_s R\sqrt{\alpha}} \tag{32}
a=ρscsR23k,b=ρscsRα3k(32)
则：
s
Θ
s
(
s
)
−
θ
i
=
−
a
Θ
s
(
s
)
−
b
s
Θ
s
(
s
)
(33)
s\Theta_s(s) - \theta_i = -a \Theta_s(s) - b \sqrt{s} \Theta_s(s) \tag{33}
sΘs(s)−θi=−aΘs(s)−bsΘs(s)(33)
整理得：
Θ
s
(
s
)
(
s
+
a
+
b
s
)
=
θ
i
(34)
\Theta_s(s) \left( s + a + b\sqrt{s} \right) = \theta_i \tag{34}
Θs(s)(s+a+bs)=θi(34)
所以：
Θ
s
(
s
)
=
θ
i
s
+
a
+
b
s
(35)
\Theta_s(s) = \frac{\theta_i}{s + a + b\sqrt{s}} \tag{35}
Θs(s)=s+a+bsθi(35)
3.5 拉普拉斯逆变换
3.5.1 铁球温度的逆变换
将
Θ
s
(
s
)
\Theta_s(s)
Θs(s) 展开为级数：
Θ
s
(
s
)
=
θ
i
∑
n
=
0
∞
(
−
1
)
n
(
a
+
b
s
)
n
s
−
n
−
1
(36)
\Theta_s(s) = \theta_i \sum_{n=0}^{\infty} (-1)^n (a + b\sqrt{s})^n s^{-n-1} \tag{36}
Θs(s)=θin=0∑∞(−1)n(a+bs)ns−n−1(36)
应用二项式定理：
(
a
+
b
s
)
n
=
∑
k
=
0
n
(
n
k
)
a
n
−
k
b
k
s
k
/
2
(37)
(a + b\sqrt{s})^n = \sum_{k=0}^{n} \binom{n}{k} a^{n-k} b^k s^{k/2} \tag{37}
(a+bs)n=k=0∑n(kn)an−kbksk/2(37)
代入得：
Θ
s
(
s
)
=
θ
i
∑
n
=
0
∞
∑
k
=
0
n
(
−
1
)
n
(
n
k
)
a
n
−
k
b
k
s
k
/
2
−
n
−
1
(38)
\Theta_s(s) = \theta_i \sum_{n=0}^{\infty} \sum_{k=0}^{n} (-1)^n \binom{n}{k} a^{n-k} b^k s^{k/2 - n - 1} \tag{38}
Θs(s)=θin=0∑∞k=0∑n(−1)n(kn)an−kbksk/2−n−1(38)
利用拉普拉斯逆变换公式：
L
−
1
[
s
−
ν
]
=
t
ν
−
1
Γ
(
ν
)
,
ν
>
0
(39)
\mathcal{L}^{-1}[s^{-\nu}] = \frac{t^{\nu-1}}{\Gamma(\nu)}, \quad \nu > 0 \tag{39}
L−1[s−ν]=Γ(ν)tν−1,ν>0(39)
对每一项进行逆变换：
L
−
1
[
s
k
/
2
−
n
−
1
]
=
t
n
+
1
−
k
/
2
−
1
Γ
(
n
+
1
−
k
/
2
)
=
t
n
−
k
/
2
Γ
(
n
+
1
−
k
/
2
)
(40)
\mathcal{L}^{-1}[s^{k/2 - n - 1}] = \frac{t^{n + 1 - k/2 - 1}}{\Gamma(n + 1 - k/2)} = \frac{t^{n - k/2}}{\Gamma(n + 1 - k/2)} \tag{40}
L−1[sk/2−n−1]=Γ(n+1−k/2)tn+1−k/2−1=Γ(n+1−k/2)tn−k/2(40)
因此，铁球过余温度的时域解为：
θ
s
(
t
)
=
θ
i
∑
n
=
0
∞
∑
k
=
0
n
(
−
1
)
n
(
n
k
)
a
n
−
k
b
k
t
n
−
k
/
2
Γ
(
n
+
1
−
k
/
2
)
(41)
\theta_s(t) = \theta_i \sum_{n=0}^{\infty} \sum_{k=0}^{n} (-1)^n \binom{n}{k} a^{n-k} b^k \frac{t^{n - k/2}}{\Gamma(n + 1 - k/2)} \tag{41}
θs(t)=θin=0∑∞k=0∑n(−1)n(kn)an−kbkΓ(n+1−k/2)tn−k/2(41)
铁球实际温度为：
T
s
(
t
)
=
T
0
+
θ
s
(
t
)
(42)
T_s(t) = T_0 + \theta_s(t) \tag{42}
Ts(t)=T0+θs(t)(42)
3.5.2 土壤温度场的逆变换
由式 (25)：
Θ
(
r
,
s
)
=
R
r
Θ
s
(
s
)
e
−
(
r
−
R
)
s
/
α
(43)
\Theta(r,s) = \frac{R}{r} \Theta_s(s) e^{-(r-R)\sqrt{s/\alpha}} \tag{43}
Θ(r,s)=rRΘs(s)e−(r−R)s/α(43)
代入
Θ
s
(
s
)
\Theta_s(s)
Θs(s) 的级数表达式：
Θ
(
r
,
s
)
=
θ
i
R
r
∑
n
=
0
∞
∑
k
=
0
n
(
−
1
)
n
(
n
k
)
a
n
−
k
b
k
s
k
/
2
−
n
−
1
e
−
(
r
−
R
)
s
/
α
(44)
\Theta(r,s) = \theta_i \frac{R}{r} \sum_{n=0}^{\infty} \sum_{k=0}^{n} (-1)^n \binom{n}{k} a^{n-k} b^k s^{k/2 - n - 1} e^{-(r-R)\sqrt{s/\alpha}} \tag{44}
Θ(r,s)=θirRn=0∑∞k=0∑n(−1)n(kn)an−kbksk/2−n−1e−(r−R)s/α(44)
利用拉普拉斯逆变换公式：
L
−
1
[
s
−
ν
e
−
c
s
]
=
1
Γ
(
ν
)
t
ν
−
1
exp
(
−
c
2
4
t
)
,
c
=
r
−
R
α
(45)
\mathcal{L}^{-1}[s^{-\nu} e^{-c\sqrt{s}}] = \frac{1}{\Gamma(\nu)} t^{\nu-1} \exp\left(-\frac{c^2}{4t}\right), \quad c = \frac{r-R}{\sqrt{\alpha}} \tag{45}
L−1[s−νe−cs]=Γ(ν)1tν−1exp(−4tc2),c=αr−R(45)
对每一项进行逆变换：
L
−
1
[
s
k
/
2
−
n
−
1
e
−
(
r
−
R
)
s
/
α
]
=
t
n
−
k
/
2
Γ
(
n
+
1
−
k
/
2
)
exp
(
−
(
r
−
R
)
2
4
α
t
)
(46)
\mathcal{L}^{-1}[s^{k/2 - n - 1} e^{-(r-R)\sqrt{s/\alpha}}] = \frac{t^{n - k/2}}{\Gamma(n + 1 - k/2)} \exp\left(-\frac{(r-R)^2}{4\alpha t}\right) \tag{46}
L−1[sk/2−n−1e−(r−R)s/α]=Γ(n+1−k/2)tn−k/2exp(−4αt(r−R)2)(46)
因此，土壤过余温度的时域解为：
θ
(
r
,
t
)
=
θ
i
R
r
∑
n
=
0
∞
∑
k
=
0
n
(
−
1
)
n
(
n
k
)
a
n
−
k
b
k
t
n
−
k
/
2
Γ
(
n
+
1
−
k
/
2
)
exp
(
−
(
r
−
R
)
2
4
α
t
)
(47)
\theta(r,t) = \theta_i \frac{R}{r} \sum_{n=0}^{\infty} \sum_{k=0}^{n} (-1)^n \binom{n}{k} a^{n-k} b^k \frac{t^{n - k/2}}{\Gamma(n + 1 - k/2)} \exp\left(-\frac{(r-R)^2}{4\alpha t}\right) \tag{47}
θ(r,t)=θirRn=0∑∞k=0∑n(−1)n(kn)an−kbkΓ(n+1−k/2)tn−k/2exp(−4αt(r−R)2)(47)
土壤实际温度为：
T
(
r
,
t
)
=
T
0
+
θ
(
r
,
t
)
(48)
T(r,t) = T_0 + \theta(r,t) \tag{48}
T(r,t)=T0+θ(r,t)(48)
4. 最终解
4.1 铁球温度
T s ( t ) = T 0 + ( T i − T 0 ) ∑ n = 0 ∞ ∑ k = 0 n ( − 1 ) n ( n k ) a n − k b k t n − k / 2 Γ ( n + 1 − k / 2 ) (49) T_s(t) = T_0 + (T_i - T_0) \sum_{n=0}^{\infty} \sum_{k=0}^{n} (-1)^n \binom{n}{k} a^{n-k} b^k \frac{t^{n - k/2}}{\Gamma(n + 1 - k/2)} \tag{49} Ts(t)=T0+(Ti−T0)n=0∑∞k=0∑n(−1)n(kn)an−kbkΓ(n+1−k/2)tn−k/2(49)
其中：
a
=
3
k
ρ
s
c
s
R
2
,
b
=
3
k
ρ
s
c
s
R
α
a = \frac{3k}{\rho_s c_s R^2}, \quad b = \frac{3k}{\rho_s c_s R\sqrt{\alpha}}
a=ρscsR23k,b=ρscsRα3k
4.2 土壤温度场
T ( r , t ) = T 0 + ( T i − T 0 ) R r ∑ n = 0 ∞ ∑ k = 0 n ( − 1 ) n ( n k ) a n − k b k t n − k / 2 Γ ( n + 1 − k / 2 ) exp ( − ( r − R ) 2 4 α t ) (50) T(r,t) = T_0 + (T_i - T_0) \frac{R}{r} \sum_{n=0}^{\infty} \sum_{k=0}^{n} (-1)^n \binom{n}{k} a^{n-k} b^k \frac{t^{n - k/2}}{\Gamma(n + 1 - k/2)} \exp\left(-\frac{(r-R)^2}{4\alpha t}\right) \tag{50} T(r,t)=T0+(Ti−T0)rRn=0∑∞k=0∑n(−1)n(kn)an−kbkΓ(n+1−k/2)tn−k/2exp(−4αt(r−R)2)(50)
5. 单重级数表达式
5.1 铁球温度的单重级数解
通过重新组织双重级数，铁球温度可表示为：
T s ( t ) = T 0 + θ i ∑ m = 0 ∞ C m t m / 2 T_s(t) = T_0 + \theta_i \sum_{m=0}^{\infty} C_m t^{m/2} Ts(t)=T0+θim=0∑∞Cmtm/2
其中系数 C m C_m Cm 由下式给出：
C m = ∑ n , k 2 n − k = m ( − 1 ) n ( n k ) a n − k b k 1 Γ ( n + 1 − k / 2 ) C_m = \sum_{\substack{n,k \\ 2n-k=m}} (-1)^n \binom{n}{k} a^{n-k} b^k \frac{1}{\Gamma(n + 1 - k/2)} Cm=n,k2n−k=m∑(−1)n(kn)an−kbkΓ(n+1−k/2)1
更具体地，对于每个 m m m，我们需要找到所有满足 2 n − k = m 2n - k = m 2n−k=m 的 ( n , k ) (n,k) (n,k) 对，其中 n ≥ ⌈ m / 2 ⌉ n \geq \lceil m/2 \rceil n≥⌈m/2⌉， k = 2 n − m k = 2n - m k=2n−m，且 0 ≤ k ≤ n 0 \leq k \leq n 0≤k≤n。
5.2 土壤温度场的单重级数解
土壤温度场可表示为：
T ( r , t ) = T 0 + θ i R r ∑ m = 0 ∞ C m t m / 2 exp ( − ( r − R ) 2 4 α t ) T(r,t) = T_0 + \theta_i \frac{R}{r} \sum_{m=0}^{\infty} C_m t^{m/2} \exp\left(-\frac{(r-R)^2}{4\alpha t}\right) T(r,t)=T0+θirRm=0∑∞Cmtm/2exp(−4αt(r−R)2)
其中系数 C m C_m Cm 与铁球温度解中的系数完全相同。
5.3 算法步骤
- 确定最大项数 M M M
- 预计算系数 C m C_m Cm， m = 0 , 1 , 2 , … , M m = 0, 1, 2, \ldots, M m=0,1,2,…,M
- 对于每个
m
m
m：
- 初始化 C m = 0 C_m = 0 Cm=0
- 对于
n
=
⌈
m
/
2
⌉
n = \lceil m/2 \rceil
n=⌈m/2⌉ 到
N
max
N_{\text{max}}
Nmax：
- 计算 k = 2 n − m k = 2n - m k=2n−m
- 如果
0
≤
k
≤
n
0 \leq k \leq n
0≤k≤n：
C m ← C m + ( − 1 ) n ( n k ) a n − k b k 1 Γ ( n + 1 − k / 2 ) C_m \leftarrow C_m + (-1)^n \binom{n}{k} a^{n-k} b^k \frac{1}{\Gamma(n + 1 - k/2)} Cm←Cm+(−1)n(kn)an−kbkΓ(n+1−k/2)1
- 存储系数 C m C_m Cm 供后续使用
6. 算法实现
6.1 Python实现
如上所述，单重级数的表达式较双重级数的计算量有了明显降低，基于Python实现了该理论解的求解：
import numpy as np
from scipy.special import gamma, binom
import matplotlib.pyplot as plt
class SingleSeriesSphereSolution:
def __init__(self, T_i, T_0, R, k, alpha, rho_s, c_s, M_max=100):
"""
单重级数解实现
参数:
T_i: 铁球初始温度 [K]
T_0: 土壤初始温度 [K]
R: 铁球半径 [m]
k: 土壤导热系数 [W/m·K]
alpha: 土壤热扩散率 [m²/s]
rho_s: 铁球密度 [kg/m³]
c_s: 铁球比热容 [J/kg·K]
M_max: 最大项数
"""
self.T_i = T_i
self.T_0 = T_0
self.theta_i = T_i - T_0
self.R = R
self.k = k
self.alpha = alpha
self.rho_s = rho_s
self.c_s = c_s
# 计算参数
self.a = 3 * k / (rho_s * c_s * R**2)
self.b = 3 * k / (rho_s * c_s * R * np.sqrt(alpha))
# 预计算系数
self.M_max = M_max
self.C_coeffs = self._compute_coefficients(M_max)
print(f"单重级数解初始化完成，预计算了 {M_max} 项系数")
print(f"参数: a = {self.a:.6f}, b = {self.b:.6f}")
def _compute_coefficients(self, M_max):
"""计算单重级数系数 C_m"""
C_coeffs = np.zeros(M_max)
N_max = 2 * M_max # 确保覆盖所有可能的(n,k)对
for m in range(M_max):
C_m = 0.0
# 寻找满足 2n - k = m 的(n,k)对
for n in range((m + 1) // 2, N_max):
k = 2 * n - m
if 0 <= k <= n:
term = ((-1)**n * binom(n, k) *
self.a**(n - k) * self.b**k /
gamma(n + 1 - k/2))
C_m += term
C_coeffs[m] = C_m
return C_coeffs
def sphere_temperature(self, t, M_terms=None):
"""计算铁球温度"""
if t <= 0:
return self.T_i
if M_terms is None:
M_terms = self.M_max
else:
M_terms = min(M_terms, self.M_max)
theta_sum = 0.0
for m in range(M_terms):
theta_sum += self.C_coeffs[m] * t**(m/2)
return self.T_0 + self.theta_i * theta_sum
def soil_temperature(self, r, t, M_terms=None):
"""计算土壤温度场"""
if t <= 0:
if r < self.R:
return self.T_i
else:
return self.T_0
if M_terms is None:
M_terms = self.M_max
else:
M_terms = min(M_terms, self.M_max)
if r < self.R:
# 铁球内部，近似为铁球表面温度
return self.sphere_temperature(t, M_terms)
theta_sum = 0.0
for m in range(M_terms):
theta_sum += self.C_coeffs[m] * t**(m/2)
# 添加径向衰减项
radial_factor = self.R / r
exponential_factor = np.exp(-(r - self.R)**2 / (4 * self.alpha * t))
return self.T_0 + self.theta_i * radial_factor * theta_sum * exponential_factor
def get_coefficient_info(self, max_display=10):
"""显示系数信息"""
print("\n单重级数系数信息:")
print("m\tC_m\t\t|C_m|")
print("-" * 40)
for m in range(min(max_display, self.M_max)):
print(f"{m}\t{self.C_coeffs[m]:.6e}\t{abs(self.C_coeffs[m]):.6e}")
# 分析系数衰减
abs_coeffs = np.abs(self.C_coeffs[:max_display])
decay_ratio = abs_coeffs[1:] / abs_coeffs[:-1]
print(f"\n系数衰减比 (平均): {np.mean(decay_ratio):.4f}")
def plot_temperature_profiles(self, t_values, r_locations=None, M_terms=100):
"""绘制温度分布图"""
if r_locations is None:
r_locations = [self.R, self.R*1.1, self.R*1.5, self.R*2, self.R*3, self.R*4]
print(f"r_locations: {r_locations}")
plt.figure(figsize=(12, 10))
# 铁球温度随时间变化
plt.subplot(2, 2, 1)
T_sphere = [self.sphere_temperature(t, M_terms) for t in t_values]
plt.plot(t_values/3600, T_sphere, 'r-', linewidth=2, label='Sphere Temperature')
plt.xlabel('Time (hours)')
plt.ylabel('Temperature (K)')
plt.title('Sphere Temperature vs Time')
plt.grid(True, alpha=0.3)
# 土壤不同位置温度
plt.subplot(2, 2, 2)
colors = plt.cm.viridis(np.linspace(0, 1, len(r_locations)))
for i, r in enumerate(r_locations):
T_profile = [self.soil_temperature(r, t, M_terms) for t in t_values]
label = f'r = {r:.2f} m' if r > self.R else 'Sphere Surface'
plt.plot(t_values/3600, T_profile, color=colors[i], linewidth=2, label=label)
plt.xlabel('Time (hours)')
plt.ylabel('Temperature (K)')
plt.title('Soil Temperature at Different Locations')
plt.legend()
plt.grid(True, alpha=0.3)
# 径向温度分布
plt.subplot(2, 2, 3)
r_array = np.linspace(self.R, self.R*3, 100)
time_points = [3600, 3600*6, 3600*12, 3600*24] # 1, 6, 12, 24小时
colors = plt.cm.plasma(np.linspace(0, 1, len(time_points)))
for i, t in enumerate(time_points):
T_radial = [self.soil_temperature(r, t, M_terms) for r in r_array]
plt.plot(r_array, T_radial, color=colors[i], linewidth=2, label=f't = {t/3600:.1f} h')
plt.axvline(x=self.R, color='k', linestyle='--', alpha=0.5, label='Sphere Surface')
plt.xlabel('Radial Distance (m)')
plt.ylabel('Temperature (K)')
plt.title('Radial Temperature Distribution')
plt.legend()
plt.grid(True, alpha=0.3)
# 系数衰减
plt.subplot(2, 2, 4)
m_values = np.arange(min(20, self.M_max))
abs_coeffs = np.abs(self.C_coeffs[m_values])
plt.plot(m_values, abs_coeffs, 'bo-', linewidth=2, markersize=6)
plt.xlabel('Coefficient Index m')
plt.ylabel('|C_m|')
plt.title('Coefficient Decay')
plt.grid(True, alpha=0.3)
plt.tight_layout()
plt.show()
# 使用示例
if __name__ == "__main__":
# 参数设置
T_i, T_0 = 373.15, 293.15 # 100°C -> 20°C
R, k, alpha = 0.1, 0.5, 1e-6
rho_s, c_s = 7800, 500
# 创建求解器
solver = SingleSeriesSphereSolution(T_i, T_0, R, k, alpha, rho_s, c_s, M_max=100)
# 显示系数信息
solver.get_coefficient_info()
# 计算温度
t_test = 3600 # 1小时
T_sphere = solver.sphere_temperature(t_test)
T_soil = solver.soil_temperature(R*1.5, t_test)
print(f"\n计算结果:")
print(f"铁球温度 (t={t_test/3600}h): {T_sphere:.2f} K")
print(f"土壤温度 (r={R*1.5:.2f}m, t={t_test/3600}h): {T_soil:.2f} K")
# 绘制温度分布
t_array = np.logspace(0, 5, 50) # 1秒到约1天
solver.plot_temperature_profiles(t_array)
# 收敛性测试
print("\n收敛性测试 (t=1h):")
for M in [5, 10, 15, 20, 25, 30]:
T_val = solver.sphere_temperature(3600, M)
print(f"M = {M:2d}: T = {T_val:.4f} K")
仿真结果如下：
由此可知，单重级数的表达式在 M M M 较大时，能够很好地逼近铁球的温度分布。
6.2 FDM方法实现
为了验证理论解的正确性，基于显式有限差分法（FDM）实现了数值解的求解：
import numpy as np
import matplotlib.pyplot as plt
from scipy.special import gamma, binom
import time
class SingleSeriesSphereSolution:
def __init__(self, T_i, T_0, R, k, alpha, rho_s, c_s, M_max=30):
"""单重级数解实现"""
self.T_i = T_i
self.T_0 = T_0
self.theta_i = T_i - T_0
self.R = R
self.k = k
self.alpha = alpha
self.rho_s = rho_s
self.c_s = c_s
# 计算参数
self.a = 3 * k / (rho_s * c_s * R**2)
self.b = 3 * k / (rho_s * c_s * R * np.sqrt(alpha))
# 预计算系数
self.M_max = M_max
self.C_coeffs = self._compute_coefficients(M_max)
def _compute_coefficients(self, M_max):
"""计算单重级数系数 C_m"""
C_coeffs = np.zeros(M_max)
N_max = 2 * M_max # 确保覆盖所有可能的(n,k)对
for m in range(M_max):
C_m = 0.0
# 寻找满足 2n - k = m 的(n,k)对
for n in range((m + 1) // 2, N_max):
k = 2 * n - m
if 0 <= k <= n:
term = ((-1)**n * binom(n, k) *
self.a**(n - k) * self.b**k /
gamma(n + 1 - k/2))
C_m += term
C_coeffs[m] = C_m
return C_coeffs
def sphere_temperature(self, t, M_terms=None):
"""计算铁球温度"""
if t <= 0:
return self.T_i
if M_terms is None:
M_terms = self.M_max
else:
M_terms = min(M_terms, self.M_max)
theta_sum = 0.0
for m in range(M_terms):
theta_sum += self.C_coeffs[m] * t**(m/2)
return self.T_0 + self.theta_i * theta_sum
class FDMSphereSolution:
def __init__(self, T_i, T_0, R, k, alpha, rho_s, c_s, r_max=1.0, Nr=200):
"""有限差分法实现"""
self.T_i = T_i
self.T_0 = T_0
self.R = R
self.k = k
self.alpha = alpha
self.rho_s = rho_s
self.c_s = c_s
self.r_max = r_max
self.Nr = Nr
# FDM网格设置
self.dr = (r_max - R) / (Nr - 1)
self.r_nodes = np.linspace(R, r_max, Nr)
# 铁球热容和表面积
self.C_sphere = rho_s * c_s * (4/3) * np.pi * R**3
self.A_sphere = 4 * np.pi * R**2
def solve_fdm(self, t_max, dt, save_times=None):
"""有限差分法求解"""
if save_times is None:
save_times = np.linspace(0, t_max, 100)
# 初始化
T_soil = np.full(self.Nr, self.T_0)
T_sphere = self.T_i
# 存储结果
results = {
'times': [],
'T_sphere': []
}
# 时间循环
t_current = 0
save_index = 0
while t_current <= t_max:
# 保存结果
if save_index < len(save_times) and t_current >= save_times[save_index] - dt/2:
results['times'].append(t_current)
results['T_sphere'].append(T_sphere)
save_index += 1
# 创建新数组
T_new = T_soil.copy()
# 更新土壤内部节点
for i in range(1, self.Nr-1):
r = self.r_nodes[i]
term1 = (T_soil[i+1] - 2*T_soil[i] + T_soil[i-1]) / self.dr**2
term2 = (2/r) * (T_soil[i+1] - T_soil[i-1]) / (2*self.dr)
T_new[i] = T_soil[i] + self.alpha * dt * (term1 + term2)
# 边界条件
# 计算热流
dT_dr = (T_soil[1] - T_soil[0]) / self.dr
heat_flux = -self.k * dT_dr
# 更新铁球温度
delta_T = - (heat_flux * self.A_sphere * dt) / self.C_sphere
T_sphere_new = T_sphere + delta_T
T_sphere_new = max(T_sphere_new, self.T_0)
# 更新边界
T_new[0] = T_sphere_new
T_new[-1] = self.T_0
# 更新
T_soil = T_new
T_sphere = T_sphere_new
t_current += dt
return results
# 参数设置
T_i, T_0 = 373.15, 293.15 # 100°C -> 20°C
R, k, alpha = 0.1, 0.5, 1e-6
rho_s, c_s = 7800, 500
print("计算参数:")
print(f"铁球初始温度: {T_i} K")
print(f"土壤初始温度: {T_0} K")
print(f"铁球半径: {R} m")
print(f"土壤导热系数: {k} W/(m·K)")
print(f"土壤热扩散率: {alpha} m²/s")
# 创建单重级数求解器
print("\n初始化单重级数求解器...")
series_solver = SingleSeriesSphereSolution(T_i, T_0, R, k, alpha, rho_s, c_s, M_max=30)
# 创建FDM求解器
print("初始化FDM求解器...")
fdm_solver = FDMSphereSolution(T_i, T_0, R, k, alpha, rho_s, c_s, r_max=0.5, Nr=100)
# 计算时间步长
dt_max = 0.5 * fdm_solver.dr**2 / alpha
dt_use = 0.3 * dt_max
print(f"FDM时间步长: {dt_use:.2f} s")
# 时间范围
t_max = 3600 * 24 # 24小时
t_array = np.logspace(0, np.log10(t_max), 100) # 对数时间点
# 计算单重级数解
print("计算单重级数解...")
start_time = time.time()
T_series = [series_solver.sphere_temperature(t) for t in t_array]
series_time = time.time() - start_time
print(f"单重级数计算时间: {series_time:.2f} s")
# 计算FDM解
print("计算FDM解...")
start_time = time.time()
fdm_results = fdm_solver.solve_fdm(t_max, dt_use, save_times=t_array)
fdm_time = time.time() - start_time
print(f"FDM计算时间: {fdm_time:.2f} s")
# 提取FDM结果
t_fdm = np.array(fdm_results['times'])
T_fdm = np.array(fdm_results['T_sphere'])
# 绘制对比图
plt.figure(figsize=(10, 6))
plt.semilogx(t_array/3600, T_series, 'b-', linewidth=2, label='Single Series Solution')
plt.semilogx(t_fdm/3600, T_fdm, 'r--', linewidth=2, label='FDM Solution')
plt.axhline(y=T_0, color='k', linestyle=':', alpha=0.5, label='Initial Soil Temperature')
plt.xlabel('Time (hours)')
plt.ylabel('Sphere Temperature (K)')
plt.title('Comparison of Sphere Temperature: Single Series vs FDM')
plt.legend()
plt.grid(True, alpha=0.3)
plt.ylim(T_0 - 5, T_i + 5)
# 添加计算时间信息
plt.figtext(0.02, 0.02, f"Single Series: {series_time:.2f}s, FDM: {fdm_time:.2f}s",
fontsize=10, bbox=dict(boxstyle="round,pad=0.3", facecolor="white", alpha=0.8))
plt.tight_layout()
plt.show()
# 打印最终温度对比
print(f"\n最终温度对比 (t = {t_max/3600:.1f} 小时):")
print(f"单重级数解: {T_series[-1]:.2f} K")
print(f"FDM解: {T_fdm[-1]:.2f} K")
print(f"温度差: {abs(T_series[-1] - T_fdm[-1]):.2f} K")
# 计算平均误差
# 在相同时间点上比较
T_series_interp = np.interp(t_fdm, t_array, T_series)
mean_error = np.mean(np.abs(T_series_interp - T_fdm))
max_error = np.max(np.abs(T_series_interp - T_fdm))
print(f"\n误差分析:")
print(f"平均绝对误差: {mean_error:.3f} K")
print(f"最大绝对误差: {max_error:.3f} K")
运行结果如下：
最终温度对比 (t = 24.0 小时):
单重级数解: 297.40 K
FDM解: 297.93 K
温度差: 0.53 K
误差分析:
平均绝对误差: 0.753 K
最大绝对误差: 1.089 K
铁球温度降低对比图如下：
如上对比证实了解析结果的正确性。
扫一扫
77
被折叠的 条评论
为什么被折叠?
到【灌水乐园】发言
