有限体积法（1）——一维扩散方程的推导

薛定谔的alpha狗

已于 2022-09-16 23:06:37 修改

阅读量8.9k

点赞数 16

分类专栏： CFD 文章标签： cfd

于 2020-06-15 07:39:54 首次发布

本文链接：https://blog.csdn.net/weixin_42562856/article/details/106750975

版权

CFD 专栏收录该内容

16 篇文章 104 订阅

订阅专栏

对于变量 $\phi$ 的输运方程为：

$\frac{\partial ( \rho \phi)}{\partial t} + \nabla \cdot (\rho \phi \bold u) = \nabla \cdot (\Gamma \nabla \phi) + S_\phi \tag{1}$

其中， $\Gamma$ 为扩散系数。
方程 $(1)$ 从左到右的各项分别是时间项、对流项、扩散项和源项。将方程 $(1)$ 中略去时间项和对流项就是稳态扩散方程：
$\nabla \cdot (\Gamma \nabla \phi) + S_\phi = 0 \tag{2}$

对方程 $(2)$ 在有限控制体内积分，并根据高斯散度定理有，
$\begin{aligned} &\int_{CV} \nabla \cdot (\Gamma \nabla \phi) dV + \int_{CV} S_\phi dV \\ \\ &= \int_{\tilde A} \bold n \cdot (\Gamma \nabla \phi) dA + \int_{CV} S_\phi dV = 0 \end{aligned} \tag{3}$

其中， $\bold n$ 为边界面 $\tilde A$ 的法向量。
考虑一维模型，扩散方程为
$\frac{d}{dx}\left( \Gamma \frac{d\phi}{dx} \right) + S = 0 \tag{4}$

要离散方程 $(4)$ 首先要生成网格，如下图，将一维计算域分成5段（即5个网格单元），每个网格单元就是一个有限控制体。两端边界为A和B，每个网格单元有一个节点，如W、P和E，节点一般为网格单元的中心点，也是变量保存的位置。两个相邻节点的中点是网格单元的边界。
网格划分
网格与边界之间的距离关系如下图，
在这里插入图片描述
扩散方程的离散主要包括两个部分：散度的离散和梯度的离散。

散度的离散

根据方程式 $(3)$ ，一维扩散方程 $(4)$ 在控制体单元 $\Delta V$ （边界为 $\tilde A$ ）上的积分有，
$\begin{aligned} &\int_{\Delta V} \frac{d}{dx} \left( \Gamma \frac{d \phi}{dx} \right)dV + \int_{\Delta V} SdV \\ \\ &= \int_{\tilde A} \bold n \cdot \left( \Gamma \frac{d \phi}{dx} \bold i \right) dA + \int_{\Delta V}SdV \\ \\ &= \left( \Gamma A \frac{d \phi}{dx} \right)_e - \left( \Gamma A \frac{d \phi}{dx} \right)_w + \bar{S} \Delta V = 0 \end{aligned} \tag{5}$

其中， $A$ 是边界面的面积， $\bar{S}$ 为控制体单元内的平均值。从上图中可见，在边界 $w$ 处，边界的法向量 $\bold n$ 是指向x轴负向的，所以 $\bold n \cdot \bold i = -1$ 。

梯度的离散

方程式 $(5)$ 仅是半离散方程，需要进一步对梯度项 $\frac{d \phi}{dx}$ 进行离散。从式中看见，梯度项都是位于单元边界面位置的，而边界处是不保存变量值，边界面两边的节点会保存变量值，因此可以用有限差分法求解边界面的导数，比如用二阶中心差分格式，
$\left( \Gamma \frac{d \phi}{dx} \right)_w = \Gamma_w \frac{\phi_P - \phi_W}{\delta x_{WP}} \tag{6a}$

$\left( \Gamma \frac{d \phi}{dx} \right)_e = \Gamma_e \frac{\phi_E - \phi_P}{\delta x_{PE}} \tag{6b}$

源项 $\bar S \Delta V$ 通常是变量 $\phi_P$ 的函数，这里假设为线性关系，即
$\bar S \Delta V = S_u + S_P \phi_P \tag{7}$

将式 $(6 a)$ 、 $(6 b)$ 和 $(7)$ 带入到式 $(5)$ 中，
$\Gamma_e A_e \frac{\phi_E - \phi_P}{\delta x_{PE}} - \Gamma_w A_w \frac{\phi_P - \phi_W}{\delta x_{WP}} + (S_u + S_P \phi_P) = 0 \tag{8}$
重新整理方程式 $(8)$ 得，
$\left( \frac{\Gamma_e}{\delta x_{PE}} + \frac{\Gamma_w}{\delta x_{WP}} - S_P\right) \phi_P = \left( \frac{\Gamma_w}{\delta x_{WP}} \right) \phi_W + \left( \frac{\Gamma_e}{\delta x_{PE}} \right) \phi_E + S_u \tag{9}$
简化一下，
$a_P \phi_P = a_W \phi_W + a_E \phi_E + S_u \tag{10}$
其中各系数为
$\begin{cases} a_W &= \frac{\Gamma_w}{\delta x_{WP}} \\ \\ a_E &= \frac{\Gamma_e}{\delta x_{PE}} \\ \\ a_P &= a_W + a_E - S_P \end{cases}\tag{11}$

方程式 $(5)$ 中的扩散系数 $\Gamma$ 有两种方法来计算：算术平均法（相当于线性插值）和调和平均法（详见文献[2]）。这里先采用简单的线性插值来表示，即
$\Gamma_w = \frac{\Gamma_W + \Gamma_P}{2} \tag{12a}$

$\Gamma_e = \frac{\Gamma_P + \Gamma_E}{2} \tag{12b}$

参考资料：

Versteeg H K , Malalasekera W . An introduction to computational fluid dynamics : the finite volume method = 计算流体动力学导论[M]. 世界图书出版公司, 2010.
陶文铨. 数值传热学-第2版[M]. 西安交通大学出版社, 2001.