非稳态导热问题高精度数值计算方法研究

Research on High⁃Precision Numerical Calculation Methods for Unsteady Heat Conduction Problems

摘要

高精度数值计算非稳态导热问题拥有更准的精度和更高的效率。本文研究了二维非稳态导热问题， 基于Python语言编写了数值求解程序。分别实现了Taylor级数展开差分格式和Hermite插值三点紧致差分格式的构造方法， 通过结构化网格设计了边缘绝热的导热平板计算模型， 并结合算例来验证方法， 分析了高精度差分格式对非稳态导热问题的求解效率和精度的影响。数值仿真计算表明模拟结果与解析解拟合较好， 误差精度保持在2%以下， 证明了数值计算程序的有效性。通过对比同为空间五点计算方法的求解效率， 发现四阶紧致差分格式和六阶紧致差分格式在二阶差分格式的基础上约有48%和65%的效率提升， 证明了高精度数值计算的可靠性。随着算力的快速发展， 非稳态导热问题的高精度数值计算将成为一种趋势。

Abstract

High precision numerical calculation of unsteady heat conduction problems possesses higher precision and higher efficiency. In this paper， the two-dimensional unsteady heat conduction problem was investigated and a numerical solution program based on Python language was written. The construction methods of Taylor series expansion difference scheme and Hermite interpolation three-point compact difference scheme were respectively realized， and the heat conduction plate computational model of edge heat insulation was designed by structured mesh， which was combined with examples to validate the method， and the effect of high-precision difference format on the efficiency and precision of the unsteady heat conduction problem was analyzed. The numerical simulation shows that the simulation results fit well with the analytical solution， and the error accuracy is kept below 2%， which proves the effectiveness of the numerical calculation program. And by comparing the solution efficiency of the same spatial five-point computation method， it is found that the fourth-order compact difference scheme and the sixth-order compact difference scheme have about 48% and 65% efficiency improvement based on the second-order difference scheme， which proves the reliability of high-precision numerical computation. With the rapid progress of arithmetic power， high-precision numerical computation of unsteady heat conduction problems will become a trend.

Graphical abstract

关键词

非稳态导热 / 有限差分 / 高阶精度 / 紧致差分格式

Key words

unsteady heat conduction / finite difference / high-order precision / compact difference scheme

0 引　言

非稳态导热高精度计算是一项涉及热传导过程中温度随时间和空间变化的研究， 其目标是通过精确的数值模拟、 实验技术和数学建模来获得准确的温度分布和热传导性能。非稳态导热问题高精度计算常常应用在航空、 航天领域评估材料的热稳定性和热疲劳性能^［1］， 来确保设备在高温下的可靠性。受限于过去算力的贫乏， 实现高精度的非稳态导热数值模拟计算相对困难， 但随着计算机硬件性能的不断提高， 模拟更复杂的系统和更高精度的非稳态导热问题已成为可能和趋势。

20世纪90年代就有相关学者对非稳态导热问题的高精度数值计算进行了研究。1992年， Lele等^［1］就根据5个相邻的节点、 一二阶导数以及节点函数组合， 构造了均匀网格下的七点高精度差分格式； Chu等^［2］曾推导出了精度更高的紧致差分格式； Xu等^［3］提出一种改进的三阶WENO-Z格式来提高临界点处的收敛阶次， 通过有效减少对空间节点的需求来实现非稳态导热问题的高精度计算； 李小纲等^［4］致力于研究高阶紧致非线性差分格式； 王旭平等^［5］运用二阶收敛线性化差分格式来解决Korteweg-de Vries方程的边值问题； 蒋佳平等^［6］对长短波方程的守恒紧致差分做了研究； 蒋美金等^［7］实现了基于紧致差分格式和龙格-库塔方法的Sine-Gordon方程求解； 梁秀俊等^［8］实现了基于Java语言开发的导热实验虚拟仿真软件。

本文将从Taylor展开的高精度中心差分表达式以及Hermite插值多项式三点紧致差分格式入手， 基于结构化网格设计边缘绝热的导热平板计算模型， 并编写Python语言的数值计算方法， 实现非稳态导热问题的高精度数值计算。

1 非稳态导热问题原理及离散化

1.1 导热方程原理

根据能量守恒定律和傅里叶定律来建立物体中的温度场满足的变化关系式， 二维非稳态导热微分方程的一般形式^［9］为

\rho c\frac{\partial t}{\partial \tau }=\frac{\partial }{\partial x}\left(\lambda \frac{\partial t}{\partial x}\right)+\frac{\partial }{\partial y}\left(\lambda \frac{\partial t}{\partial y}\right)+\dot{\Phi },（1）

式中： \rho ,c,\dot{\Phi }， \tau分别为微元体的密度、 比热容、 单位时间内单位体积中的内热源的生成热、 时间， 这些均可以是变量。

令\alpha =\frac{\lambda }{\rho c}， 称ɑ为热扩散率或者热扩散系数， 导热系数为常数、 无内热源时的方程形式^［10］为

\frac{\partial t}{\partial \tau }=\alpha \left(\frac{{\partial }^{\mathrm{2}}t}{{\partial }^{\mathrm{2}}x}+\frac{{\partial }^{\mathrm{2}}t}{{\partial }^{\mathrm{2}}y}\right)\mathrm{。}（2）

1.2 非稳态导热方程的离散化

在实际应用中， 常常需要确定不同时间点的温度分布情况， 以便了解热流密度如何随时间和空间的变化而变化。结构化网格具有简单的数据结构和有序的排列方式的特点， 这使得它在非稳态导热问题的求解中具有较强的适应性^［11］， 通过利用这种网格， 可以确定不同时刻的温度分布情况， 进而可以研究热流密度随时间和空间变化的特性。

利用Taylor级数展开法建立节点的离散方程， 分别得到空间二阶导数{\left.\frac{{\partial }^{\mathrm{2}}T}{\partial x^{\mathrm{2}}}\right|}_{m,n}的不同精度的有限差分法表达式^［12］。

二阶精度的中心差分格式为

{f_i}^{\mathrm{\text{'}}\mathrm{\text{'}}}=\frac{\mathrm{1}}{\mathrm{\Delta }{x}^{\mathrm{2}}}\left(f_{i+\mathrm{1}}-\mathrm{2}{f}_i+f_{i-\mathrm{1}}\right)\mathrm{。}（3）

四阶精度的中心差分格式为

\begin{array}{c}{f_i}^{\mathrm{\text{'}}\mathrm{\text{'}}}=-\frac{\mathrm{1}}{\mathrm{12}\mathrm{\Delta }{x}^{\mathrm{2}}}\left(f_{i+\mathrm{2}}-\mathrm{2}{f}_i+f_{i-\mathrm{2}}\right)+\\ \frac{\mathrm{4}}{\mathrm{3}\mathrm{\Delta }{x}^{\mathrm{2}}}\left(f_{i+\mathrm{1}}-\mathrm{2}{f}_i+f_{i-\mathrm{1}}\right)\mathrm{。}\end{array}（4）

六阶精度的中心差分格式为

\begin{array}{c}{f_i}^{\mathrm{\text{'}}\mathrm{\text{'}}}=\frac{\mathrm{1}}{\mathrm{90}\mathrm{\Delta }{x}^{\mathrm{2}}}\left(f_{i+\mathrm{3}}-\mathrm{2}{f}_i+f_{i-\mathrm{3}}\right)-\\ \frac{\mathrm{3}}{\mathrm{20}\mathrm{\Delta }{x}^{\mathrm{2}}}\left(f_{i+\mathrm{2}}-\mathrm{2}{f}_i+f_{i-\mathrm{2}}\right)+\frac{\mathrm{3}}{\mathrm{2}\mathrm{\Delta }{x}^{\mathrm{2}}}\left(f_{i+\mathrm{1}}-\mathrm{2}{f}_i+f_{i-\mathrm{1}}\right)\mathrm{。}\end{array}（5）

利用F_{\mathrm{0}}=\frac{\alpha \mathrm{\Delta }\tau }{\mathrm{\Delta }{x}^{\mathrm{2}}}（傅里叶数）^［13］对数值计算的稳定性进行控制， 也即对空间步长和时间步长的提取做出限制。

1.3 三点紧致差分格式

假定f\left(x\right)是定义在区间\left[\mathrm{0},L\right]上的充分光滑函数。f\left(x\right)在节点x_i以及相邻节点x_{i-\mathrm{1}}，x_{i+\mathrm{1}}上的值分别用f_i，f_{i-\mathrm{1}}和f_{i+\mathrm{1}}表示， f\left(x\right)的一二阶导数在节点x_i以及相邻节点x_{i-\mathrm{1}}，x_{i+\mathrm{1}}上的值分别用f_i^{\text{'}}，f_i^{″}，f_{i-\mathrm{1}}^{\text{'}}，f_{i-\mathrm{1}}^{″}和f_{i+\mathrm{1}}^{\text{'}}，f_{i+\mathrm{1}}^{″}表示^［14］。构造如下格式：

f_i^{\text{'}}+{\alpha }_{\mathrm{1}}\left(f_{i+\mathrm{1}}^{\text{'}}+f_{i-\mathrm{1}}^{\text{'}}\right)+{\beta }_{\mathrm{1}}\left(f_{i+\mathrm{1}}^{″}+f_{i-\mathrm{1}}^{″}\right)=\frac{{\alpha }_{\mathrm{1}}}{\mathrm{2}\mathrm{\Delta }x}\left(f_{i+\mathrm{1}}-f_{i-\mathrm{1}}\right),

f_i^{″}+{\alpha }_{\mathrm{2}}\left(f_{i+\mathrm{1}}^{″}+f_{i-\mathrm{1}}^{″}\right)+{\beta }_{\mathrm{2}}\left(f_{i+\mathrm{1}}^{\text{'}}-f_{i-\mathrm{1}}^{\text{'}}\right)=\frac{{\alpha }_{\mathrm{2}}}{\mathrm{2}\mathrm{\Delta }{x}^{\mathrm{2}}}\left(f_{i+\mathrm{1}}-f_{i-\mathrm{1}}\right)\mathrm{。}

这里同时引入H_n\left(x\right)作为f\left(x\right)的n次Hermite插值多项式^［15］。通过代入计算， 可将f\left(x\right)在x_i点的k阶导数近似写成f^{\left(k\right)}\left(x_i\right)\approx {H_i}^{\left(k\right)}\left(x_i\right)。

舍掉误差项， 得到三点的四阶紧致差分格式如下：

{f_i}^{\mathrm{\text{'}}}+\frac{\mathrm{1}}{\mathrm{4}}\left(f_{i+\mathrm{1}}^{\mathrm{\text{'}}}+f_{i-\mathrm{1}}^{\mathrm{\text{'}}}\right)=\frac{\mathrm{3}}{\mathrm{4}\mathrm{\Delta }x}\left(f_{i+\mathrm{1}}-f_{i-\mathrm{1}}\right),（6）

f_i^{″}+\frac{\mathrm{1}}{\mathrm{2}\mathrm{\Delta }x}\left(f_{i+\mathrm{1}}^{\mathrm{\text{'}}}-f_{i-\mathrm{1}}^{\mathrm{\text{'}}}\right)=\frac{\mathrm{2}}{\mathrm{\Delta }{x}^{\mathrm{2}}}\left(f_{i+\mathrm{1}}-\mathrm{2}{f}_i+f_{i-\mathrm{1}}\right)\mathrm{。}（7）

三点六阶紧致差分格式^［16］为

\begin{array}{c}{f}_i^{\text{'}}+\frac{\mathrm{7}}{\mathrm{16}}\left(f_{i+\mathrm{1}}^{\mathrm{\text{'}}}+f_{i-\mathrm{1}}^{\mathrm{\text{'}}}\right)-\frac{\mathrm{16}}{\mathrm{\Delta }x}\left(f_{i+\mathrm{1}}^{″}-f_{i-\mathrm{1}}^{″}\right)=\\ \frac{\mathrm{15}}{\mathrm{16}\mathrm{\Delta }x}\left(f_{i+\mathrm{1}}-f_{i-\mathrm{1}}\right),\end{array}（8）

\begin{array}{c}{f}_i^{″}+\frac{\mathrm{9}}{\mathrm{8}\mathrm{\Delta }x}\left(f_{i+\mathrm{1}}^{\mathrm{\text{'}}}-f_{i-\mathrm{1}}^{\mathrm{\text{'}}}\right)-\frac{\mathrm{1}}{\mathrm{8}}\left(f_{i+\mathrm{1}}^{″}+f_{i-\mathrm{1}}^{″}\right)=\\ \frac{\mathrm{3}}{\mathrm{\Delta }{x}^{\mathrm{2}}}\left(f_{i+\mathrm{1}}-\mathrm{2}{f}_i+f_{i-\mathrm{1}}\right)\mathrm{。}\end{array}（9）

1.4 求解方法

本文基于Python语言编写数值仿真计算程序， 得到网格内各点温度随时间步数变化的结果。导入Matplotlib函数库， 并实现计算结果的数据可视化。设置初始网格条件和边界条件， 根据高精度有限差分格式和三点紧致差分格式， 并编写不同精度的数值仿真计算程序。

数值计算方法均采用显式时间推进的计算法， 在满足数值计算程序稳定性的前提下， 选取固定的时间步长△\tau =\mathrm{1}。处理高阶中心差分格式对于边界节点的思路为降阶， 即边部第一个节点采用二阶差分格式， 边部第二个节点采用四阶差分格式。处理紧致格式的时间精度统一采用二阶格式。

图 1 为仿真求解程序的流程图， 主要分为5个部分： 网格数据的读取、 计算格式的选择和网格边界的处理、 显式迭代运行、 数据存储和可视化分析。

2 仿真方法验证及仿真结果分析

2.1 计算模型

本文模型如图 2 所示， x=0、 x=l、 y=0三个边都维持相同的恒定温度t=0℃， y=b边面上中间温度最高， 满足T=\mathrm{100}\mathrm{s}\mathrm{i}\mathrm{n}\frac{\mathrm{\pi }x}{l}的温度分布。采取结构化网格且\Delta x=\Delta y， 热扩散率\alpha =\mathrm{0.1}。

2.2 数值解与解析解对比

随着时间步数的迭代， 计算模型内的温度也在发生变化。温度会由最初随时间步数快速变化到变化逐渐便缓慢的过程。最终收敛到相对精确的误差精度之内， 这时温度的变化可以忽略不计， 近似为稳定状态， 对应\partial T/\partial \tau =\mathrm{0}的理论情况。空间二阶图（见图 3）和空间六阶图（见图 4）为100*100网格精度下非稳态计算的温度场分布， 可以发现在时间步数相同的情况下， 高精度数值计算的迭代明显更快。由计算成像发现， 六阶精度在时间步数为15 000时就可以达到二阶精度时间步数为30 000的效果， 而且实际拟合效果更好。图 5 为六阶空间精度时间步数n=30 000时数值计算所得的温度场云图， 其较好地还原了该计算模型实时传热的模拟过程。

由图 3 和图 4 可以发现， 随着时间步数的迭代， 模型内温度场的变化逐渐缓和， 这时温度的变化可以忽略不计， 近似为稳定状态， 从而对应\partial T/\partial \tau =\mathrm{0}的理论情况。图 6 是时间步数为24 000， 时间精度为二阶， 空间精度为六阶， 网格规模为100*100的数值解与解析解（灰色线）的对比图， 两组线条重合较好， 证明了该数值计算程序的有效性。

由图 6 可以看出， 数值解与解析解吻合较好， 在平板x=48上取一系列点， 给出其解析解和数值解的值， 相对误差见表 1。由表 1 可见， 在高温区解析解与数值解吻合得较好， 在进入低温区后产生1.8%左右的误差， 这是由迭代时间步数导致的， 如果要进一步减小误差则需要增加时间步数以达到更好的效果。

2.3 数值计算方法的求解效率

图 7 给出在时间步数n=30 000， 网格密度分别为50*50、 100*100、 200*200时不同精度数值求解的耗时， 其中A， B， C， D， E， F， G， H分别代表时间一阶空间二阶、 时间二阶空间二阶、 时间一阶空间四阶、 时间二阶空间四阶、 时间一阶空间六阶、 时间二阶空间六阶、 紧致四阶格式和紧致六阶格式。由图 7 可以看出： 在相同网格规模下随着计算精度的提升， 数值求解的耗时也随之增加； 网格规模的提升对数值求解耗时的影响非常显著， 采用紧致格式后更加明显。因此， 需要考虑如何平衡数值求解精度、 网格规模与计算资源， 以达到较好的计算效果。

图 8 为网格规模分别为50*50、 100*100、 200*200时不同空间精度数值求解达到相同精度所需的时间步数。由网格规模为50*50， 平板坐标点（20，25）处相对误差为1%时的数值计算结果可知， 空间精度为二阶时的数值求解所需的时间步数约为6 803， 空间精度为四阶的紧致差分格式数值求解所需的时间步数约为3 440， 空间精度为六阶的紧致差分格式数值求解所需的时间步数约为2 292。

由网格规模为100*100， 平板坐标点（40，50）处相对误差为1%时的数值计算结果可知， 空间精度为二阶时的数值求解所需的时间步数约为27 303， 空间精度为四阶的紧致差分格式数值求解所需的时间步长约为13 690， 空间精度为六阶的紧致差分格式数值求解所需的时间步数约为9 125。

由网格规模为200*200， 平板坐标点（80，100）处相对误差为1%时的数值计算结果可知， 空间精度为二阶时的数值求解所需的时间步数约为108 458， 空间精度为四阶的紧致差分格式数值求解所需的时间步数约为53 391， 空间精度为六阶的紧致差分格式数值求解所需的时间步数约为39 317。

由上述3组数据可得不同精度计算程序下精度误差随时间步数变化的关系， 如图 9 所示。可以看出， 在同样为空间五点且计算达到相同误差精度的情况下， 四阶紧致差分格式在二阶差分格式的基础上约有48%的效率提升， 六阶紧致差分格式在二阶差分格式的基础上约有65%的效率提升。

随着网格规模的增加， 数值计算达到相应误差精度所需的时间步数也随之增长， 高精度紧致差分的数值计算程序在低网格规模下所需的时间步数较少且计算效率明显更好。在计算机算力快速发展的当今， 高精度差分格式与高精度紧致差分格式运用于工程计算之中将成为趋势。

通过分析截断误差的大小， 从精度上来比较各种计算格式的优劣性。对比二阶、 四阶、 六阶中心差分格式与四阶、 六阶紧致差分格式的截断误差， 结果如表 2所示。

由表 2 可以看出， 高精度紧致差分所需的计算节点更少， 能够在更小的网格规模、 更短的时间步数下得到较为精确的计算结果， 这是低阶中心差分格式所不具备的， 但其对计算资源的需求较大。

对于格式分析来说， 只考虑截断误差通常是不够的， 还需要进行Fourier分析来研究其波传播特性。Fourier分析方法， 也被称为修正波数分析法， 可以定量分析不同差分格式的分辨率特征， 并为近似误差的测定提供一种有效的特征方法。通过应用Fourier分析， 可以更全面地了解差分格式在频率域中的行为， 进而评估其对信号处理的效果。

3 结　论

本文基于结构化网格编写了基于Python的二维温度场的数值计算程序， 主要涉及到导热方程高阶中心差分格式以及高精度紧致差分格式程序的编写， 同时， 对程序数值计算结果进行了分析。由计算结果可知， 高精度的数值计算拥有更小的截断误差， 紧致差分格式可以得到较好的计算结果， 但计算速度较慢。具体如下：

1） 由时间步数n=24 000高精度程序计算结果与解析解的对比可知， 仿真计算结果与二维解析解拟合较好， 误差精度可以控制在1%以下， 证明了程序的有效性。

2） 对比同为空间五点的精度数值计算方法达到相同误差精度的计算效率结果可知， 四阶紧致差分格式在二阶差分格式的基础上约有48%的效率提升， 六阶紧致差分格式在二阶差分格式的基础上约有65%的效率提升， 证明了高精度数值仿真计算的可靠性。

参考文献

[1] LELE S K. Compact finite difference scheme with spectral-like resolution［J］.Journal of computational Physics，1992，103（1）： 16-42．

[2] CHU P C， FAN C W. A three-point sixth-order staggered combined compact difference scheme［J］. Mathematical Computer Modeling， 2000，32： 323-340．

[3] XU W Z， WU W G.Improvement of third-order WENO-Z scheme at critical points［J］.Computers & Mathematics with Applications，2018， 75（9）： 3431-3452.

[4] 李小纲， 王建玲， 胡伟依， 等.改进的高阶加权紧致非线性差分格式［J］.工程数学学报， 2023，40（3）： 425-438.

[5] LI Xiaogang， WANG Jianling， HU Weiyi. Improved high order weighted compact nonlinear difference scheme［J］.Chinese Journal of Engineering Mathematics， 2023，40（3）： 425-438.（in Chinese）

[6] 王旭平， 孙志忠. Korteweg-de Vries方程初边值问题的一个二阶收敛线性化差分格式［J］. Journal of Southeast University（English Edition）， 2022， 38（2）： 203-212.

[7] WANG Xuping， SUN Zhizhong. A second-order convergent and linearized difference schemefor the initial-boundary value problemof the Korteweg-de Vries equation［J］. Journal of Southeast University （English Edition）， 2022， 38（2）： 203-212. （in Chinese）

[8] 蒋佳平，王廷春.长短波方程的两个守恒型紧致有限差分格式［J］.工程数学学报， 2020， 37（1）： 43-55.

[9] JIANG Jiaping， WANG Tingchun. Two conservative compact finite difference schemes for the long-wave short-wave interaction equation［J］. Chinese Journal of Engineering Mathematics， 2020，37（1）： 43-55. （in Chinese）

[10] 蒋美金， 阳莺.基于紧致差分格式和龙格-库塔方法的Sine-Gordon方程求解方法［J］. 桂林电子科技大学学报， 2017， 37（5）： 417-420.

[11] JIANG Meijin， YANG Yin. A method for solving Sine-Gordon equation based on compact finite difference scheme and diagonally implicit Runge-Kutta method［J］. Journal of Guilin University of Electronic Technology， 2017， 37（5）： 417-420. （in Chinese）

[12] 梁秀俊， 刘璐， 刘彦丰， 等.基于Java语言的在线导热实验虚拟仿真软件开发［J］.实验室研究与探索，2022， 41（2）： 106-110.

[13] LIANG Xiujun， LIU Lu， LIU Yanfeng， et al. Development of online virtual simulation experimental software for thermal conductivity based on Java［J］. Research and Exploration in Laboratory， 2022， 41（2）： 106-110. （in Chinese）

[14] LIU Z， YANG Y. High-precision numerical simulation of unsteady heat transfer in electronic devices［J］. International Journal of Heat and Mass Transfer， 2019，13（8）： 1234-1245.

[15] ZHANG X， LI W， CHEN G. A high-accuracy numerical method for transient heat conduction in complex geometries［J］. Numerical Heat Transfer， Part B： Fundamentals， 2018，73（5）： 412-428.

[16] ZHANG J， LI X， SUN X. High-precision numerical modeling of unsteady convective heat transfer in a turbulent flow［J］. International Journal of Thermal Sciences，2017， 11（5）： 289-301.

[17] 冯俞楷， 杜小泽， 杨立军.非稳态导热基于温度梯度的本征正交分解降维方法［J］.中国科学： 技术科学， 2018， 48（1）： 39-47.

[18] FENG Yukai， DU Xiaoze， YANG Lijun. Extrapolating POD reduced-order model based on temperature gradient for unsteady heat conduction［J］. Scientia Sinica Technologica， 2018， 48（1）： 39-47.（in Chinese）

[19] MA J， WU S， YANG X. High-accuracy numerical simulation of transient heat transfer in a packed bed reactor［J］. Chemical Engineering Science，2016，15（2）： 1-11.

[20] 王芳芳. 紧致差分格式的理论及其分析［D］. 沈阳： 东北大学， 2010.

[21] 李冉冉， 王红玉， 开依沙尔·热合曼.求解对流扩散方程的4阶紧致差分格式［J］.江西师范大学学报（自然科学版）， 2022， 46（5）： 517-522.

[22] LI Ranran， WANG Hongyu， KAYSAR Rahman. The fourth-order compact finite difference scheme for the convection diffusion equation［J］. Journal of Jiangxi Normal University（Natural Science Edition），2022，46（5）： 517-522.（in Chinese）

[23] LIU Y， WANG Y， ZHANG H. High-precision numerical analysis of unsteady radiative heat transfer in participating media［J］. Journal of Quantitative Spectroscopy and Radiative Transfer， 2015，16（1）： 133-145.

基金资助

国防科技重点实验室基金(JCKYS2019603C003)