双层球对称热传导问题的解及数值验证

期刊菜单

双层球对称热传导问题的解及数值验证
The Solution and Numerical Verification of the Two-Layer Spherical Symmetric Heat Conduction Problem

DOI: 10.12677/mos.2024.134425, PDF, HTML, XML, 下载: 11 浏览: 22
作者: 白宇龙：西安建筑科技大学理学院，陕西西安
关键词: 间断的函数正交系；双层瞬态热传导；间断的解析解；不光滑的解；Orthogonal System of Discontinuous Series； Double-Layer Transient Heat Conduction； Discontinuous Analytic Solution； The Solution Is Not Smooth

摘要: 具有两种材料双层热传导问题，由于具有两种传热性质不同的材料，标准的三角函数正交系难以满足两种材料界面处温度和热流量的连续条件，问题难以求解。本文基于三角函数在正交性上良好的性质，根据函数系正交性的定义，用连续条件构造了一个间断的三角函数正交系，并且该正交函数系的系数中含有待定函数。利用该正交函数系级数展开的方法可以满足连续条件、边界条件和任意的初始条件，并且其系数中的待定函数由连续条件和初始条件共同确定，而在一般情况下系数中的待定函数是由初始条件决定的。该解析解与数值模拟对比一致性较好。文中提出的求解方法构造了一种新的间断的正交函数系，为多相介质的热传导问题提供了解析解。也为其他多相介质求解析解或者求解不光滑的解析解提供了一种新的思路。

Abstract: Because there are two kinds of materials with different heat transfer properties, the standard trigonometric orthogonal system makes it difficult to meet the continuous conditions of temperature and heat flow at the interface of the two materials, and the problem is difficult to solve. In this paper, a discontinuous orthogonal system of trigonometric functions is constructed based on the good properties of orthogonality of trigonometric functions and the definition of orthogonality of function systems, and its coefficients contain undetermined functions. The continuous condition, boundary condition and any initial condition can be satisfied by the method of series expansion of the orthogonal system of functions, and the undetermined function in the coefficient is determined by the continuous condition and the initial condition together. In contrast, the initial condition, in general, determines the undetermined function in the coefficient. The analytical solution is in good agreement with the numerical simulation. The method proposed in this paper constructs a new discontinuous orthogonal function system, which provides an analytical solution for the heat conduction problem in multiphase media. It also provides a new way of thinking for solving analytic solutions of other polyphase media or solving non-smooth analytical solutions.

文章引用：白宇龙. 双层球对称热传导问题的解及数值验证[J]. 建模与仿真, 2024, 13(4): 4694-4701. https://doi.org/10.12677/mos.2024.134425

1. 引言

热传导方程，是偏微分方程中最重要的基本方程，是空气动力学、弹性力学、电动力学、流体动力学、近场动力学、量子力学、等离子体物理、数字信号处理等领域的基本理论方程[1] [2]；多相介质热传导问题的研究方法，对上述领域相关问题研究大多有借鉴意义[3]；也即，热传导问题的研究成果可为上述力学领域的研究提供理论支撑。

对于经典的一维有限区间热传导问题，将材料的热参数假定为已知常数，标准的三角函数正交系[4]就可以满足边界条件和任给的初始条件。对于无限区域上的热传导问题，可利用Fourier变换和Laplace变换等积分变换方法[5]求解，这也是这类问题的基本解。

从1822年第一个热传导模型被解出以来[6]，热传导方程一直是数学物理方法研究的基础内容之一[7] [8]。但至今为止还有一些比较常见的问题难以直接用现有的正交函数系求解，如多相介质的热传导问题[9] [1 0 ]，或非线性方程[11]等问题。在求解两种介质中的温度场分布的过程中，正交函数系需要满足界面处的温度和热流量连续条件，一般连续的正交函数系难以满足，因此无法直接求解。

正因为界面处连续条件的特殊性，导致这方面的理论研究较少；但两相热传导问题在自然界中十分普遍[12]-[14]，所以他的解析解在相应的研究中十分重要[15] [16]。针对上述问题，无法根据现有的文献求解，故本文基于分离变量法[17]，根据界面处的连续条件[18]构造了间断的正交函数系，该函数系可以满足边界条件和界面处的连续条件，用级数叠加的方式满足了任給的初始条件，求出了该问题的解。这种求解方法，可为类似的多相热传导问题[19] [20]或其他多相的偏微分方程问题提供一种求得解析解的方法。

2. 基本模型

2.1. 控制方程的通解

首先确定一个半径为r₁和r₂的同心圆球域(见图1)。原点O位于圆域的中心(坐标(0,0))。其中， $I_{1} \subset [0, r_{1}]$ ， $I_{2} \subset (r_{1}, r_{2}]$ 。

Figure 1. Spherical coordinates of the double layer

图1. 双层球面坐标

蓝色部分为区域I₁，黄色部分为区域I₂。控制方程分别如下：

$\frac{\partial u_{1}}{\partial t} = \frac{k_{1}}{c_{1} ρ_{1}} Δ u_{1} \begin{matrix} \end{matrix} 0 \leq r \leq r_{1}, t \geq 0$

$\frac{\partial u_{2}}{\partial t} = \frac{k_{2}}{c_{2} ρ_{2}} Δ u_{2} \begin{matrix} \end{matrix} 0 \leq r \leq r_{2}, t \geq 0$

u为温度分布函数，k，ρ，c分别是导热系数、密度、比热容。下标1和2分别代表区域I₁和I₂。

球对称的热传导方程式(1)如下

$\frac{\partial u_{1}}{\partial t} = 2 \frac{k_{1}}{c_{1} ρ_{1}} \frac{1}{r} \frac{\partial u_{1}}{\partial r} + \frac{k_{1}}{c_{1} ρ_{1}} \frac{\partial^{2} u_{1}}{\partial r^{2}}$

用变量分离法[17]可得，通解为

$u_{1} = \frac{1}{r} C_{1} \cos (β r) e^{- β^{2} \frac{k_{1}}{c_{1} ρ_{1}} t} + \frac{1}{r} C_{2} \sin (β r) e^{- β^{2} \frac{k_{1}}{c_{1} ρ_{1}} t}$ (1)

相同的，区域I₂的通解为

$u_{2} = \frac{1}{r} D_{1} \cos (ω (r - r_{2})) e^{- ω^{2} \frac{k_{2}}{c_{2} ρ_{2}} t} + \frac{1}{r} D_{2} \sin (ω (r - r_{2})) e^{- ω^{2} \frac{k_{2}}{c_{2} ρ_{2}} t}$ (2)

其中 $β$ 和 $ω$ 是特征值，C₁、C₂、D₁和D₂是待定常数。

2.2. 边界条件，初始条件和连续条件

球对称的状态下， $r = 0$ 处的边界条件为

${\frac{\partial u_{1}}{\partial r} |}_{r = 0} = 0$ (3)

$r = r_{2}$ 处的边界条件为

${(α_{3} \frac{\partial u_{2}}{\partial r} + α_{4} u_{2}) |}_{r = r_{2}} = 0$ (4)

其中 $α_{3}$ 和 $α_{4}$ 是常数。 $r = r_{1}$ 处的连续条件为

${u_{1} |}_{r = r_{1}} = {u_{2} |}_{r = r_{1}}$ (5)

${k_{1} \frac{\partial u_{1}}{\partial r} |}_{r = r_{1}} = {k_{2} \frac{\partial u_{2}}{\partial r} |}_{r = r_{1}}$ (6)

初始条件分别为

${u_{1} |}_{t = 0} = φ_{1} (r)$ (7)

${u_{2} |}_{t = 0} = φ_{2} (r)$ (8)

3. 两相介质热传导问题的解

3.1. 满足边界条件

将通解(1)式代入边界条件(3)式，得待定系数 $C_{1} = 0$ ，将通解(2)式代入边界条件(4)式，得待定系数， $(α_{3} / r_{2} - α_{4}) D_{1} = ω α_{3} D_{2}$ 。简写为 $D_{1} = B_{2} D$ ， $D_{2} = B_{1} D$ ，其中， $B_{1} = α_{3} / r_{2} - α_{4}$ ， $B_{2} = ω α_{3}$ ，得到满足边界条件的通解如下

$u_{1} = \frac{1}{r} C_{2} \sin (β r) e^{- β^{2} \frac{k_{1}}{c_{1} ρ_{1}} t}$ (9)

$u_{2} = \frac{1}{r} D (B_{2} \cos (ω (r - r_{2}))) e^{- ω^{2} \frac{k_{2}}{c_{2} ρ_{2}} t} + \frac{1}{r} D (B_{1} \sin (ω (r - r_{2}))) e^{- ω^{2} \frac{k_{2}}{c_{2} ρ_{2}} t}$ (10)

3.2. 满足连续条件的间断正交函数系

将满足边界条件的解式(9)和(10)代入连续条件(5)式和(6)式，分别得

$C_{2} \sin (β r_{1}) e^{- β^{2} \frac{k_{1}}{c_{1} ρ_{1}} t} = D (B_{2} \cos (ω (r_{1} - r_{2}))) e^{- ω^{2} \frac{k_{2}}{c_{2} ρ_{2}} t} + D (B_{1} \sin (ω (r_{1} - r_{2}))) e^{- ω^{2} \frac{k_{2}}{c_{2} ρ_{2}} t}$

$\begin{array}{l} k_{1} C_{2} [- \frac{1}{r_{1}} \sin (β r_{1}) + β \cos (β r_{1})] e^{- β^{2} \frac{k_{1}}{c_{1} ρ_{1}} t} \\ = k_{2} D [(ω B_{1} - \frac{B_{2}}{r_{1}}) \cos (ω (r_{1} - r_{2})) + (B_{1} - ω B_{2}) \sin (ω (r_{1} - r_{2}))] e^{- ω^{2} \frac{k_{2}}{c_{2} ρ_{2}} t} \end{array}$

从上两式中可以看出，要满足所有 $t > 0$ 的连续条件，前提是变量分离后的时间项相等，如下

$e^{- β^{2} \frac{k_{1}}{c_{1} ρ_{1}} t} = e^{- ω^{2} \frac{k_{2}}{c_{2} ρ_{2}} t}$

由此可得 $ω^{2} = β^{2} (k_{1} c_{2} ρ_{2}) / (k_{2} c_{1} ρ_{1})$ ，简写 $ω = β K$ ，为了方便证明函数系的正交性，将通解简写为

$u_{1 n} = \frac{1}{r} C_{2 n} ϕ_{1 n} T_{1 n}$ (11)

$u_{2 n} = \frac{1}{r} D_{n} ϕ_{2 n} T_{1 n}$ (12)

其中 $ϕ_{1 n} = \sin (β_{n} r)$ ， $ϕ_{2 n} = B_{2 n} \cos (K β_{n} (r - r_{2})) + B_{1} \sin (K β_{n} (r - r_{2}))$ ， $T_{1 n} = e^{- β_{n}^{2} \frac{k_{1}}{c_{1} ρ_{1}} t}$ ，整理可得，

$\int_{0}^{r_{1}} \frac{ϕ_{1 n} ϕ_{1 m}}{ϕ_{1 n} (r_{1}) ϕ_{1 m} (r_{1})} d r = \frac{1}{β_{n}^{2} - β_{m}^{2}} [(\frac{{ϕ^{'}}_{1 m} (r_{1})}{ϕ_{1 m} (r_{1})} - \frac{1}{r_{1}}) - (\frac{{ϕ^{'}}_{1 n} (r_{1})}{ϕ_{1 n} (r_{1})} - \frac{1}{r_{1}})]$ (13)

$\int_{r_{1}}^{r_{2}} \frac{K^{2} ϕ_{2 n} ϕ_{2 m}}{ϕ_{2 n} (r_{1}) ϕ_{2 m} (r_{1})} d r = - \frac{1}{β_{n}^{2} - β_{m}^{2}} [(\frac{{ϕ^{'}}_{2 m} (r_{1})}{ϕ_{2 m} (r_{1})} - \frac{1}{r_{1}}) - (\frac{{ϕ^{'}}_{2 n} (r_{1})}{ϕ_{2 n} (r_{1})} - \frac{1}{r_{1}})]$ (14)

其中 $ϕ_{1 n} (r_{1})$ ， $ϕ_{2 n} (r_{1})$ ， $ϕ_{1 m} (r_{1})$ 和 $ϕ_{2 m} (r_{1})$ 分别为函数 $ϕ_{1 n}$ $ϕ_{2 n}$ $ϕ_{1 n}$ 和 $ϕ_{2 m}$ 在 $r = r_{1}$ 处的取值。将通解(11)式和(12)式代入连续条件(5)式和(6)式，联立可得

$k_{1} (- \frac{1}{r_{1}} + \frac{{φ^{'}}_{1 n} (r_{1})}{φ_{1 n} (r_{1})}) = k_{2} (- \frac{1}{r_{1}} + \frac{{φ^{'}}_{2 n} (r_{1})}{φ_{2 n} (r_{1})})$ (15)

联立(13)式(14)式和(15)式

$k_{1} \int_{0}^{r_{1}} \frac{ϕ_{1 n} ϕ_{1 m}}{ϕ_{1 n} (r_{1}) ϕ_{1 m} (r_{1})} d r + k_{2} \int_{r_{1}}^{r_{2}} \frac{K^{2} ϕ_{2 n} ϕ_{2 m}}{ϕ_{2 n} (r_{1}) ϕ_{2 m} (r_{1})} d r = 0$ (16)

(16)式说明，在 $r \in (0, r_{2})$ 上有间断的正交函数系 $ϕ_{n}^{} (r)$ ，满足正交性如下

$\int_{0}^{r_{2}} ϕ_{n} ϕ_{m} d r = {\begin{array}{l} \int_{0}^{r_{2}} ϕ_{n}^{2} d r & n = m \\ 0 & n \neq m \end{array}$

根据函数系的正交性确定待定系数的表达式如下，

$C_{2 n} = \frac{\int_{0}^{r_{1}} r ϕ_{1 n} φ_{1} d r + \int_{r_{1}}^{r_{2}} r ϕ_{2 n} \frac{K \sqrt{k_{2}} ϕ_{1 n} (r_{1})}{ϕ_{2 n} (r_{1}) \sqrt{k_{1}}} φ_{3} d r}{\int_{0}^{r_{1}} ϕ_{1 n}^{2} d r + \int_{r_{1}}^{r_{2}} ϕ_{2 n}^{2} {[\frac{K \sqrt{k_{2}} ϕ_{1 n} (r_{1})}{ϕ_{2 n} (r_{1}) \sqrt{k_{1}}}]}^{2} d r}$ ， $D_{n} = \frac{\int_{0}^{r_{1}} r ϕ_{1 n} \frac{ϕ_{2 n} (r_{1}) \sqrt{k_{1}}}{K \sqrt{k_{2}} ϕ_{1 n} (r_{1})} φ_{4} d r + \int_{r_{1}}^{r_{2}} r ϕ_{2 n} φ_{2} d r}{\int_{0}^{r_{1}} ϕ_{1 n}^{2} {[\frac{ϕ_{2 n} (r_{1}) \sqrt{k_{1}}}{K \sqrt{k_{2}} ϕ_{1 n} (r_{1})}]}^{2} d r + \int_{r_{1}}^{r_{2}} ϕ_{2 n}^{2} d r}$

其中含有待定函数 $φ_{3}$ 和 $φ_{4}$ ，将待定系数表达式带入连续条件(5)式，求出待定函数 $φ_{3}$ 和 $φ_{4}$ ，最终得待定系数如下

$C_{2 n} = \frac{\int_{0}^{r_{1}} r ϕ_{1 n} φ_{1} d r + \int_{r_{1}}^{r_{2}} r ϕ_{2 n} \frac{K \sqrt{k_{2}} ϕ_{1 n} (r_{1})}{ϕ_{2 n} (r_{1}) \sqrt{k_{1}}} \frac{K \sqrt{k_{2}}}{\sqrt{k_{1}}} φ_{2} d r}{\int_{0}^{r_{1}} ϕ_{1 n}^{2} d r + \int_{r_{1}}^{r_{2}} ϕ_{2 n}^{2} {[\frac{K \sqrt{k_{2}} ϕ_{1 n} (r_{1})}{ϕ_{2 n} (r_{1}) \sqrt{k_{1}}}]}^{2} d r}$ (17)

$D_{n} = \frac{\int_{0}^{r_{1}} r ϕ_{1 n} \frac{ϕ_{2 n} (r_{1}) \sqrt{k_{1}}}{K \sqrt{k_{2}} ϕ_{1 n} (r_{1})} \frac{\sqrt{k_{1}}}{K \sqrt{k_{2}}} φ_{1} d r + \int_{r_{1}}^{r_{2}} r ϕ_{2 n} φ_{2} d r}{\int_{0}^{r_{1}} ϕ_{1 n}^{2} {[\frac{ϕ_{2 n} (r_{1}) \sqrt{k_{1}}}{K \sqrt{k_{2}} ϕ_{1 n} (r_{1})}]}^{2} d r + \int_{r_{1}}^{r_{2}} ϕ_{2 n}^{2} d r}$ (18)

4. 结果与分析

就前几节阐述的问题而言，本文总共分析了4种不同的参数和初始条件组合的例子(见表1)。

本文用matlab自编程序实现了理论解中特征值的取值，计算了函数系的系数，实现了理论解和数值解的可视化。自编程序代码与详细的理论推导过程见附件1，分析结果见图2~图5。

Table 1. Thermal properties

表1. 材料的热属性

例	$k_{1}$	$k_{2}$	$c_{1}$	$c_{2}$	$ρ_{1}$	$ρ_{2}$	$r_{1}$	$r_{2}$	$α_{3}$	$α_{4}$	${u_{2} \|}_{t = 0}$	${u_{1} \|}_{t = 0}$
1	1	11	2	12	3	13	4	20	1	3	20	20
2	21	51	32	42	3	33	8	10	1	1	$e^{- x}$	x
3	2	51	37	4	7	33	4	10	1	1	$e^{\sin x}$	$\sin x$
4	47	3	38	5	7	28	5	10	1	1	20	-20

(a) (b)

Figure 2. Temperature distribution (a) and error (b) of theoretical solution and numerical solution at 200 s, 400 s, 600 s, 800 s and 1000 s respectively under the parameter combination of Example 1

图2. 例1的参数组合下，理论解和数值解分别在200 s、400 s、600 s、800 s、1000 s时的温度分布(a)和误差(b)

(a) (b)

Figure 3. Temperature distribution (a) and error (b) of theoretical solution and numerical solution at 200 s, 400 s, 600 s, 800 s and 1000 s respectively under the parameter combination of Example 2

图3. 例2的参数组合下，理论解和数值解分别在200 s、400 s、600 s、800 s、1000 s时的温度分布(a)和误差(b)

(a) (b)

Figure 4. Temperature distribution (a) and error (b) of theoretical solution and numerical solution at 200 s, 400 s, 600 s, 800 s and 1000 s respectively under the parameter combination of Example 3

图4. 例3的参数组合下，理论解和数值解分别在200 s、400 s、600 s、800 s、1000 s时的温度分布(a)和误差(b)

(a) (b)

Figure 5. Temperature distribution (a) and error (b) of theoretical solution and numerical solution at 200 s, 400 s, 600 s, 800 s and 1000 s respectively under the parameter combination of Example 4

图5. 例4的参数组合下，理论解和数值解分别在200 s、400 s、600 s、800 s、1000 s时的温度分布(a)和误差(b)

图2~图5的两相材料中的温度分布结果显示，温度在两种材料的界面处连续，满足连续条件式(5)，而温度梯度在界面处会产生突变。温度梯度突变的比例满足连续条件式(6)，与两种材料对应的热传导系数k₁和k₂成反比。解析解与数值解的误差在各种参数与初始条件的组合情况下，始终小于5%，具体精度受matlab自编程序代码影响。数值解的网格密度与解析解计算系数时的积分精度，会影响解析解与数值解的误差。在以上所有情况下，解析解与数值解的温度结果高度一致。这说明该解析解是有效且准确的。可以为双层球体中的温度分布和其他不光滑的解析解求解提供参考。

5. 结论

本文通过理论分析结合数值模拟，得出结论如下。

1) 用两种介质间温度的连续条件，构造了间断的正交函数系，利用该间断的正交函数系首次推导出了任给初始条件下双层热传导问题的解析解。说明间断的正交函数系可以以级数的形式逼近初始条件，从理论求解过程可见，间断的正交函数系在两相介质中的与连续的正交函数系在单相介质中的地位相同。

2) 本研究提出的间断正交函数系，可以有效地满足两相介质的初始条件，对多相介质的问题的理论求解具有一定意义。这将为多相介质的解析解求解提供一种新的思路。并对以其他正交函数系为解答的，多相介质问题提供一种构造解析解的思路。该间断正交函数系本身也可以作为一种间断的解析解，满足相关的间断解的问题。

3) 当两相介质的“热传导系数/比热容/密度”的差异较大时，界面处的温度梯度变化较大。理论解与数值模拟间的误差在界面两侧的波动较大。不同的界面位置和不同的初始温度分布则对上述规律的影响较小。

参考文献

[1]	王彩云, 姜冬菊, 黄丹. 基于常规态型近场动力学的热力耦合变形破坏分析[J]. 应用力学学报, 2020, 37(3): 938-945.
[2]	熊辉, 杨光. 复变量热传导方程的动力系统[J]. 应用数学和力学2014, 31(9): 1655-1602.
[3]	陶月赞, 姚梅. 地下水渗流力学发展进程与动向[J]. 吉林大学学报(地球科学版), 2007, 37(2): 221-230.
[4]	韦兰英. 傅里叶级数在级数求和中的应用[J]. 高师理科学刊, 2018, 38(2): 67-69.
[5]	丁波, 欧志华, 奉瑞萍. 热传导方程的积分变换求解与MATLAB实现[J]. 湖南科技学院学报, 2021, 42(3): 2-7.
[6]	Ozisik, M.N. (1989) Boundary Value Problems of Heat Conduction. Dover Publications, New York.
[7]	姜玉山, 徐延钦, 王晓敏, 等. 数学物理方程[M]. 北京: 清华大学出版社, 2014.
[8]	苏变萍, 陈东立. 复变函数与积分变换[M]. 北京: 高等教育出版社, 2018.
[9]	张品戈, 于泳博, 兰天然. 多层热防护服的温度分布与最优厚度研究[J]. 工业技术, 2020, 18(16): 67-68.
[10]	汪易平, 丁颖, 于志财, 王震, 徐丽慧. 热防护服散热性的数学模型仿真[J]. 现代纺织技术, 2024, 32(3): 102-109.
[11]	刘建军, 叶礼友, 薛强, 何翔. 基于N-S方程的泡沫陶瓷渗流微观数值模拟[J]. 动能材料, 2007, 38(A7): 2770-2772.
[12]	周航, 邵珠山, 乔汝佳, 郭银波, 席慧慧. 工业微波照射下混凝土砂浆-骨料界面力学性能与多场效应数值分析[J]. 材料科学与工程学报, 2023, 41(1): 30-38.
[13]	田苗, 苏云, 李俊. “火灾环境-防护服-人体”CFD传热模拟及皮肤烧伤分布预测[J]. 清华大学学报(自然科学版), 2024, 64(6): 1032-1038.
[14]	齐兵, 杨松, 曹振生, 张少强, 隆能增. 高地温隧道隔热层方案优化设计及应用——以尼格隧道为例[J]. 科学技术与工程, 2023, 23(9): 4004-4010.
[15]	唐松花, 罗迎社, 彭相华. 火灾情况下混凝土板温度场的解析解[J]. 应用力学学报, 2013, 30(4): 544-549.
[16]	唐松花, 崔宇鹏. 火灾高温下混凝土柱温度场的解析解[J]. 应用力学学报, 2018, 35(6): 1361-1366.
[17]	王元明. 数学物理方程与特殊函数[M]. 北京: 高等教育出版社, 2004.
[18]	李长玉, 方彦奎, 刘福旭, 阮宇航. 热防护服-空气-皮肤热传导模型及其解析解[J]. 应用数学和力学, 2021, 42(2): 162-169.
[19]	胡亮, 马兰荣, 谷磊, 李丹丹, 韩艳浓. 高温高压对微波破岩效果的影响模拟研究[J]. 石油钻探技术, 2019, 47(2): 50-55.
[20]	李想, 赵先航, 钟华, 谢宇, 茹佳胜, 刘鑫, 李旭东, 李悦芳. 多层复合材料筒状结构在温度载荷作用下的层间应力建模与试验验证[J]. 复合材料学报, 2024, 41(8): 4371-4380.

为你推荐

友情链接