计算传热学基础

jedi-knight

已于 2022-04-10 11:05:23 修改

阅读量2.1k

点赞数 1

分类专栏：算法流体力学文章标签：学习经验分享

于 2022-04-09 22:07:07 首次发布

本文链接：https://blog.csdn.net/weixin_43325228/article/details/124051076

版权

流体力学同时被 2 个专栏收录

35 篇文章 15 订阅

订阅专栏

算法

30 篇文章 7 订阅

订阅专栏

绪论

计算传热学是指对描写流动与传热问题的控制方程采用数值方法通过计算机予以求解的一门传热学与数值方法相结合的交叉学科

传热学的研究方法包括理论分析（理论传热学），实验研究（实验传热学)，数值模拟（计算传热学）三个分支

计算传热学的计算流程：
（1）建立控制方程、初值条件和边界条件
（2）区域离散化（有限差分、有限体积、有限元）
（3）方程离散化（有限差分、有限体积、有限元）
（4）求解代数方程（直接法、间接法）
（5）判断解是否收敛，若未收敛则以当前值重建代数方程，返回第（3）步
（6）后处理

有限体积法是目前流动与传热问题的数值计算中应用最广的一种方法

有限体积法的特点
（1）将计算域划分为许多控制体，所有的变量存于控制体的几何中心。（空间离散的特点）
（2）将守恒型的控制方程对控制容积做积分导出离散方程，保证离散方程的守恒特性。（方程的离散的特点）
（3）物理概念清晰，强调控制体内物理量的守恒。

扩散问题

通用的控制方程

计算传热学面对的守恒型方程如下
$\frac{\partial \left( \rho \phi \right)}{\partial t}+\nabla \cdot \left( \rho \phi \boldsymbol{u} \right) =\nabla \cdot \left( \varGamma \nabla \phi \right) +S \\ \text{非稳态项}+\text{对流项}=\text{扩散项}+\text{源项}$ $\phi$ 取1可得连续方程，源项 $S$ 的单位为 $\mathrm{kg}/\left( \mathrm{m}^3\cdot \mathrm{s} \right)$ ，物理意义是单位时间内单位体积内加入的质量

$\phi$ 取 $\boldsymbol{u}$ 可得动量方程，此时 $\varGamma$ 为动力粘度 $\mu$ ，单位 $\mathrm{Pa}\cdot \mathrm{s}$ ，源项 $S$ 的单位为 $\mathrm{N}/\mathrm{m}^3$ ，物理意义是单位体积内作用的外力，或者说是单位时间内单位体积内加入的动量

$\phi$ 取 $c_{\mathrm{p}}T$ 可得能量方程，此时 $\varGamma$ 为 $k/c_{\mathrm{p}}$ ，单位 $\mathrm{kg}/\left( \mathrm{m}\cdot \mathrm{s} \right)$ ，源项 $S$ 的单位为 $\mathrm{W}/\mathrm{m}^3$ ，物理意义是单位体积内的功率（或吸热率），或者说是单位时间内单位体积内加入的能量（焓）

一维稳态导热方程

方程的离散

对于一维稳态导热问题，无非稳态项、无对流项，仅保留能量方程
$\frac{\mathrm{d}}{\mathrm{d}x}\left( k\frac{\mathrm{d}T}{\mathrm{d}x} \right) +S=0$ 在控制容积内积分 $\int_{\mathrm{w}}^{\mathrm{e}}{\frac{\mathrm{d}}{\mathrm{d}x}\left( k\frac{\mathrm{d}T}{\mathrm{d}x} \right)}\mathrm{d}x+\int_{\mathrm{w}}^{\mathrm{e}}{S}\mathrm{d}x=0$ 将源项写成两部分的组合（显示+隐式的形式）有
$S=S_{\mathrm{C}}+S_{\mathrm{P}}T_{\mathrm{P}}$ 并且认为源项在单元内是均匀分布的，所以积分结果为
$k_{\mathrm{e}}\left( \frac{\mathrm{d}T}{\mathrm{d}x} \right) _{\mathrm{e}}-k_{\mathrm{w}}\left( \frac{\mathrm{d}T}{\mathrm{d}x} \right) _{\mathrm{w}}+S_{\mathrm{C}}\Delta x+S_{\mathrm{P}}T_{\mathrm{P}}\Delta x=0$ 界面上的导数可以使用相邻单元中心的值差商计算，有
$k_{\mathrm{e}}\frac{T_{\mathrm{E}}-T_{\mathrm{P}}}{\delta x_{\mathrm{e}}}-k_{\mathrm{w}}\frac{T_{\mathrm{P}}-T_{\mathrm{W}}}{\delta x_{\mathrm{w}}}+S_{\mathrm{C}}\Delta x+S_{\mathrm{P}}T_{\mathrm{P}}\Delta x=0$ 整理后有
$a_{\mathrm{P}}T_{\mathrm{P}}=a_{\mathrm{W}}T_{\mathrm{W}}+a_{\mathrm{E}}T_{\mathrm{E}}+b, \\ a_{\mathrm{P}}=a_{\mathrm{W}}+a_{\mathrm{E}}-S_{\mathrm{P}}\Delta x, a_{\mathrm{W}}=\frac{k_{\mathrm{w}}}{\delta x_{\mathrm{w}}}, a_{\mathrm{E}}=\frac{k_{\mathrm{e}}}{\delta x_{\mathrm{e}}}, b=S_{\mathrm{C}}\Delta x$ 上式就是任意一个单元得到的离散形式的控制方程

总结出两大原则和四项基本法则
原则1：物理上的真实；
原则2：总的平衡；
法则1：在控制容积面上的连续性（热流密度、质量流量、动量通量）；
法则2：正系数法则，即所有的系数（ $a_{\mathrm{P}}$ 以及各相邻结点系数）必须总是正的；
法则3：源项的负斜率线性化，即 $S_{\mathrm{P}}\leqslant 0$ ；
法则4：相邻结点系数之和为 $a_{\mathrm{P}}$ ，前提是没有源项的非稳态导热

还需要注意的是，界面上的导热系数应该使用相邻节点导热系数的调和平均
$\frac{\delta x_{\mathrm{e}}}{k_{\mathrm{e}}}=\frac{\delta x_{\mathrm{Pe}}}{k_{\mathrm{P}}}+\frac{\delta x_{\mathrm{eE}}}{k_{\mathrm{E}}}, \\ \frac{\delta x_{\mathrm{w}}}{k_{\mathrm{w}}}=\frac{\delta x_{\mathrm{wP}}}{k_{\mathrm{P}}}+\frac{\delta x_{\mathrm{Ww}}}{k_{\mathrm{W}}}$

边界条件的加入

$n$ 个单元可以列写 $n$ 个方程，但是在离散边界单元时，引入了两个未知量（左侧边界一个，右侧边界还有一个）所以共有 $n + 2$ 个未知量。

为了使得方程封闭可解，需要用两个边界条件来消去多出来的两个未知量。记边界上的温度为 $T_{\mathrm{b}}$ 。

边界条件可以统一写为
$AT+B\frac{\mathrm{d}T}{\mathrm{d}x}=C$ 将该方程的导数项用差商表示，左侧的边界有
$AT_{\mathrm{b}}+B\frac{T_{\mathrm{P}}-T_{\mathrm{b}}}{\delta x_{\mathrm{w}}}=C$ 用 $T_{\mathrm{P}}$ 表示 $T_{\mathrm{b}}$ 可得
$T_{\mathrm{b}}=\frac{C\delta x_{\mathrm{w}}-BT_{\mathrm{P}}}{A\delta x_{\mathrm{w}}-B}$ 将这个方程回代到边界处单元的离散方程中即可消去 $T_{\mathrm{b}}$
同理，右侧的边界也可以这样处理，这样处理后，方程就只剩下 $n$ 个未知量，可解！
对于不同的边界条件A,B,C可以取不同的数值
对于第一类边界条件，有
$A=1, B=0, C=T_0$ 对于第二类边界条件，有
$A=0, B=-k, C=q_0$ 对于第三类边界条件，左侧有 $C=hT_{\mathrm{f}}$ 右侧有 $C=hT_{\mathrm{f}}$

方程求解

TDMA直接算法

针对三对角线性方程组，将 $T_{i-1}$ 用 $T_i$ 线性表出
$T_{i-1}=P_{i-1}T_i+Q_{i-1}$
然后代入方程组，可以得到
$T_{i}=P_{i}T_{i+1}+Q_{i}$ 进而得到 $P_{i}$ 和 $Q_{i}$ 的递推关系式
注意到 $T_{n}=Q_{n}$ ,则可以回代得到所有 $T$ 值

G-S迭代算法

$T_{\mathrm{P}}=\frac{a_{\mathrm{W}}T_{\mathrm{W}}+a_{\mathrm{E}}T_{\mathrm{E}}+b}{a_{\mathrm{P}}}$ 注意这里的等号是赋值的含义，G-S迭代法收敛的充分条件是：主对角线严格占优

一维非稳态导热方程

一维非稳态导热的控制方程为
$\frac{\partial \left( \rho cT \right)}{\partial t}=\frac{\partial}{\partial x}\left( k\frac{\partial T}{\partial x} \right) +S$ 与稳态问题类似，需要对 $t$ 和 $x$ 进行双重积分
$\int_t^{t+\Delta t}{\int_{\mathrm{w}}^{\mathrm{e}}{\frac{\partial \left( \rho cT \right)}{\partial t}\mathrm{d}x}\mathrm{d}t}=\int_t^{t+\Delta t}{\int_{\mathrm{w}}^{\mathrm{e}}{\frac{\partial}{\partial x}(k\frac{\partial T}{\partial x})\mathrm{d}x}\mathrm{d}t}+\int_t^{t+\Delta t}{\int_{\mathrm{w}}^{\mathrm{e}}{\left( S_{\mathrm{C}}+S_{\mathrm{P}}T \right) \mathrm{d}x}\mathrm{d}t}$ 积分后的结果为
$\left[ \left( \rho cT_{\mathrm{P}} \right) ^1-\left( \rho cT_{\mathrm{P}} \right) ^0 \right] \Delta x=\int_t^{t+\Delta t}{\left[ \left( k\frac{\partial T}{\partial x} \right) _{\mathrm{e}}-\left( k\frac{\partial T}{\partial x} \right) _{\mathrm{w}}+S_{\mathrm{C}}\Delta x+S_{\mathrm{P}}T_{\mathrm{P}}\Delta x \right] \mathrm{d}t}$ 再使用差商替代界面上的导数，可得
$\left[ \left( \rho cT_{\mathrm{P}} \right) ^1-\left( \rho cT_{\mathrm{P}} \right) ^0 \right] \Delta x=\int_t^{t+\Delta t}{\left[ \frac{k_{\mathrm{e}}\left( T_{\mathrm{E}}-T_{\mathrm{P}} \right)}{\delta x_{\mathrm{e}}}-\frac{k_{\mathrm{w}}\left( T_{\mathrm{P}}-T_{\mathrm{W}} \right)}{\delta x_{\mathrm{w}}}+S_{\mathrm{C}}\Delta x+S_{\mathrm{P}}T_{\mathrm{P}}\Delta x \right] \mathrm{d}t}$
关于时间的积分采用显式与隐式混合的形式，比如对 $T_\mathrm{P}$ 的积分为
$\int_t^{t+\Delta t}{T_{\mathrm{P}}\mathrm{d}t}=\left( fT_{\mathrm{P}}^{1}+\left( 1-f \right) T_{\mathrm{P}}^{0} \right) \Delta t$ 所以有
$\left[ \left( \rho cT_{\mathrm{P}} \right) ^1-\left( \rho cT_{\mathrm{P}} \right) ^0 \right] \Delta x=f\left( \frac{k_{\mathrm{e}}\left( T_{\mathrm{E}}^{1}-T_{\mathrm{P}}^{1} \right)}{\delta x_{\mathrm{e}}}-\frac{k_{\mathrm{w}}\left( T_{\mathrm{P}}^{1}-T_{\mathrm{W}}^{1} \right)}{\delta x_{\mathrm{w}}}+S_{\mathrm{P}}T_{\mathrm{P}}^{1}\Delta x \right) \Delta t \\ +\left( 1-f \right) \left( \frac{k_{\mathrm{e}}\left( T_{\mathrm{E}}^{0}-T_{\mathrm{P}}^{0} \right)}{\delta x_{\mathrm{e}}}-\frac{k_{\mathrm{w}}\left( T_{\mathrm{P}}^{0}-T_{\mathrm{W}}^{0} \right)}{\delta x_{\mathrm{w}}}+S_{\mathrm{P}}T_{\mathrm{P}}^{0}\Delta x \right) \Delta t+S_{\mathrm{C}}\Delta x\Delta t$ 当 $f$ 为1时，为全隐式格式，当 $f$ 为0时，为显示格式，当 $f$ 为0.5时，为Crank-Nicolson格式

可以进一步分析稳定性：全隐式格式无条件稳定，显示格式需要满足CFL条件
$\text{库朗数}=\frac{\varGamma \Delta t}{\rho \Delta x^2}=\frac{\alpha \Delta t}{\Delta x^2}\leqslant 0.5$ 非稳态问题的求解方法：
（1）根据初始条件得到方程组的所有系数
（2）求解方程组（隐式）或直接沿时间推进（显式）
（3）根据求解得到的结果和边界条件更新方程组
（4）返回第(2)步迭代，直到收敛

多维问题的方程离散同一维问题，这里不再推导。需要注意的是由于多维问题的方程组的系数矩阵不再是三对角阵，TDMA算法不再适用。可以使用逐列法、逐行法或ADI方法将TDMA算法与迭代法相结合，提高收敛速度。

对于逐列扫描，将 $T_{\mathrm{W}}$ 和 $T_{\mathrm{E}}$ 用上一步的值固定可以得到三对角方程;同理，对于逐行扫描，将 $T_{\mathrm{N}}$ 和 $T_{\mathrm{S}}$ 用上一步的值固定可以得到三对角方程。逐列扫描与逐行扫描交替进行即可实现交替方向隐式迭代法（ADI迭代法）

对流-扩散问题

一维稳态对流-扩散方程

$\frac{\mathrm{d}\left( \rho \phi u \right)}{\mathrm{d}x}=\frac{\mathrm{d}}{\mathrm{d}x}\left( \varGamma \frac{\mathrm{d}\phi}{\mathrm{d}x} \right)$ 对 $x$ 积分，有
$\int_{\mathrm{w}}^{\mathrm{e}}{\frac{\mathrm{d}\left( \rho \phi u \right)}{\mathrm{d}x}\mathrm{d}x}=\int_{\mathrm{w}}^{\mathrm{e}}{\frac{\mathrm{d}}{\mathrm{d}x}\left( \varGamma \frac{\mathrm{d}\phi}{\mathrm{d}x} \right) \mathrm{d}x}$ 积分结果为
$\left( \rho \phi u \right) _{\mathrm{e}}-\left( \rho \phi u \right) _{\mathrm{w}}=\left( \varGamma \frac{\mathrm{d}\phi}{\mathrm{d}x} \right) _{\mathrm{e}}-\left( \varGamma \frac{\mathrm{d}\phi}{\mathrm{d}x} \right) _{\mathrm{w}}$
为了简化方程，引入 $D$ 和 $F$ ， $D$ 表示界面处由扩散引起的质量通量， $F$ 表示界面处由对流引起的质量通量，它们的单位均为 $\mathrm{kg}/\left( \mathrm{m}^2\cdot \mathrm{s} \right)$
$F_{\mathrm{e}}\phi _{\mathrm{e}}-F_{\mathrm{w}}\phi _{\mathrm{w}}=D_{\mathrm{e}}\left( \phi _{\mathrm{E}}-\phi _{\mathrm{P}} \right) -D_{\mathrm{w}}\left( \phi _{\mathrm{P}}-\phi _{\mathrm{W}} \right) \\ D=\frac{\varGamma}{\delta x}, F=\rho u$

中心差分格式

方程左侧的界面处的对流项使用中心差分
$\begin{array}{l} F_{\mathrm{e}}\phi _{\mathrm{e}}=F_{\mathrm{e}}\frac{\phi _{\mathrm{P}}+\phi _{\mathrm{E}}}{2}\\ F_{\mathrm{w}}\phi _{\mathrm{w}}=F_{\mathrm{w}}\frac{\phi _{\mathrm{P}}+\phi _{\mathrm{W}}}{2}\\ \end{array}$ 代回方程可得
$a_{\mathrm{P}}\phi _{\mathrm{P}}=a_{\mathrm{W}}\phi _{\mathrm{W}}+a_{\mathrm{E}}\phi _{\mathrm{E}}, \\ a_{\mathrm{P}}=a_{\mathrm{W}}+a_{\mathrm{E}}+\left( F_{\mathrm{e}}-F_{\mathrm{w}} \right) ,a_{\mathrm{W}}=D_w-0.5F_w,a_{\mathrm{E}}=D_e-0.5F_e$ 但是，对于中心差分格式，为了保证稳定性，由正系数法则可得稳定性条件
$\left| Pe \right|=\frac{\left| F \right|}{D}<2$

一阶迎风格式

考虑到 $F$ 的方向性，计算界面上的物理量时使用上游节点的值，有
$\begin{array}{l} F_{\mathrm{e}}\phi _{\mathrm{e}}=\phi _{\mathrm{P}}\max \left( F_{\mathrm{e}},0 \right) -\phi _{\mathrm{E}}\max \left( -F_{\mathrm{e}},0 \right)\\ F_{\mathrm{w}}\phi _{\mathrm{w}}=\phi _{\mathrm{W}}\max \left( F_{\mathrm{w}},0 \right) -\phi _{\mathrm{P}}\max \left( -F_{\mathrm{w}},0 \right)\\ \end{array}$
这样导出的方程不会产生负的系数，解在物理上真实，求解过程是稳定的

指数格式

考虑一维稳态对流-扩散方程的精确解
$\frac{\phi -\phi _{\mathrm{P}}}{\phi _{\mathrm{E}}-\phi _{\mathrm{P}}}=\frac{e^{Pe\frac{x}{\delta x_{\mathrm{e}}}}-1}{e^{Pe}-1}, Pe=\frac{\rho u\delta x_{\mathrm{e}}}{\varGamma}$ 根据该精确解，可以精确地推算出 $\phi _{\mathrm{e}}$ 和 $\phi _{\mathrm{w}}$ 以及 $\left( \frac{\mathrm{d}\phi}{\mathrm{d}x} \right) _{\mathrm{e}}$ 和 $\left( \frac{\mathrm{d}\phi}{\mathrm{d}x} \right) _{\mathrm{w}}$ ，由此得到的格式被称为指数格式

混合格式

由于指数格式计算耗时所以又提出了混合格式来近似指数格式得到的 $a_\mathrm{E}$
$a_{\mathrm{E}}=D_{\mathrm{e}}\max \left( -Pe,1-\frac{Pe}{2},0 \right)$

幂函数格式

相比混合格式更精确地近似指数格式
$a_{\mathrm{E}}=D_{\mathrm{e}}\max \left( 0,\left( 1-0.1\frac{\left| F_{\mathrm{e}} \right|}{D_{\mathrm{e}}} \right) ^5 \right) +\max \left( 0,-F_{\mathrm{e}} \right)$
与一维对流-扩散方程的离散方式类似，可以得到二维或三维的方程

流动问题

SIMPLE算法

SIMPLE算法可以分离地求解连续方程和动量方程，其步骤如图
在这里插入图片描述
（1）假设一个初始 $u^*$ 和 $v^*$ ，用于定义动量方程的系数
（2）假设一个压力场 $p^*$
（3）根据 $p^*$ 使用动量方程解出新的 $u^*$ 和 $v^*$
（4）根据 $u^*$ 和 $v^*$ 使用压力修正方程接触 $p^{'}$
（5）根据 $p^{'}$ 进行速度修正，得到 $u^{'}$ 和 $v^{'}$
（6）更新 $u^*$ 和 $v^*$ ，用于重新定义动量方程的系数
（7）更新 $p^*$ ，返回第（3）步迭代直到收敛