CFD有限体积4-空间离散：扩散项离散

最新推荐文章于 2023-12-08 10:20:47 发布

觅°

最新推荐文章于 2023-12-08 10:20:47 发布

阅读量623

点赞数 1

分类专栏：笔记 CFD 文章标签：学习

本文链接：https://blog.csdn.net/qq_41817886/article/details/129020271

版权

笔记同时被 2 个专栏收录

11 篇文章 1 订阅

订阅专栏

CFD

8 篇文章 0 订阅

订阅专栏

文章详细介绍了在二维笛卡尔坐标系下对定常扩散方程的离散方法，包括通用形式、离散方程的建立、边界条件（Dirichlet、Neumann、混合和对称边界）的处理，以及界面扩散因子、偏斜处理、各向异性扩散的考虑。此外，还讨论了迭代过程中的松弛处理策略，以确保计算的收敛性。

摘要由CSDN通过智能技术生成

参考资料：

F. Moukalled, L. Mangani, M. Darwish. The finite volume method in computational fluid dynamics. An advanced introduction with Openfoam and Matlab[M]. Chapter 8.

1. 在二维笛卡尔坐标系下对定常扩散方程进行离散

通用形式为

$J^{\phi,D}=-\tau ^{\phi}\nabla\phi$

经积分，可得

$\sum_{f\sim nb(C)}(J^{\phi})_{f}\cdot \vec{S}_{f}=Q^{\phi}_{C}V_{C}$

展开

$J^{\phi}_{e}\cdot \vec{S}_{e}+J^{\phi}_{w}\cdot \vec{S}_{w}+J^{\phi}_{n}\cdot \vec{S}_{n}+J^{\phi}_{s}\cdot \vec{S}_{s}=Q^{\phi}_{C}V_{C}$

已知（针对网格单元时，其各界面法向量向外，如上图所示）

$\vec{S}_{e}=S_{e}\vec{i}, \vec{S}_w=-S_{w}\vec{i}, \vec{S}_{n}=S_{n}\vec{j}, \vec{S}_{s}=-S_{s}\vec{j}$

代入上式中，有

$J^{\phi}_{e}S_{e}\cdot\vec{i}-J^{\phi}_{w}S_{w}\cdot\vec{i}+J^{\phi}_{n}S_{n}\cdot\vec{j}-J^{\phi}_{s}S_{s}\cdot\vec{j}=Q^{\phi}_{C}V_{C}$

其中

$J^{\phi}_{e}S_{e}\cdot\vec{i}=-\tau^{\phi}_{e}(\nabla\phi)_{e}S_{e}\cdot\vec{i}\\=-\tau^{\phi}_{e}S_{e}(\frac{\partial \phi}{\partial x}\cdot\vec{i}+\frac{\partial \phi}{\partial y}\cdot\vec{j})_{e}\cdot \vec{i}\\=-(\tau^{\phi}S\frac{\partial \phi}{\partial x})_{e}\\=-(\tau^{\phi}\Delta y\frac{\partial \phi}{\partial x})_{e}$

现需假设C和E之间参数phi的分布型线，假设为线性分布，则有

$(\frac{\partial \phi}{\partial x})_{e}=\frac{\phi_{E}-\phi_{C}}{(\delta x)_{e}}$

代入上式中，并将其整理为以下形式

$J^{\phi}_{e}S_{e}\cdot\vec{i}=fluxC_{e}\phi_{C}+fluxN_{e}\phi_{E}+fluxV_{e}$

整理后可得，

$J^{\phi}_{e}S_{e}\cdot\vec{i}\\=(\frac{\tau^{\phi}\Delta y}{\delta x})_{e}\phi_{C}-(\frac{\tau^{\phi}\Delta y}{\delta x})_{e}\phi_{E}+0$

令 $gdiff_{e}=(\frac{\Delta y}{\delta x})_{e}=\frac{S_{e}}{d_{EC}}$ ，有

$\\fluxC_{e}=\tau^{\phi}_{e}gdiff_{e}\\fluxN_{e}=-\tau^{\phi}_{e}gdiff_{e}\\ fluxV_{e}=0$

同理

$-J^{\phi}_{w}S_{w}\cdot\vec{i}\\=\tau^{\phi}_{w}S_{w}\frac{\phi_{C}-\phi_{W}}{(\delta x)_{w}}\\=fluxC_{w}\phi_{C}+fluxN_{w}\phi_{W}+fluxV_{w}\\=\tau^{\phi}_{w}(\frac{S}{\delta x})_{w}\phi_{C}-\tau^{\phi}_{w}(\frac{S}{\delta x})_{w}\phi_{W}+0$ ,

$J^{\phi}_{n}S_{n}\cdot\vec{j}\\=-\tau^{\phi}_{n}S_{n}\frac{\phi_{N}-\phi_{C}}{(\delta y)_{n}}\\=fluxC_{n}\phi_{C}+fluxN_{n}\phi_{N}+fluxV_{n}\\=\tau^{\phi}_{n}(\frac{S}{\delta y})_{n}\phi_{C}-\tau^{\phi}_{n}(\frac{S}{\delta y})_{n}\phi_{N}+0$ ,

$-J^{\phi}_{s}S_{s}\cdot\vec{j}\\=\tau^{\phi}_{s}S_{s}\frac{\phi_{C}-\phi_{S}}{(\delta y)_{s}}\\=fluxC_{s}\phi_{C}+fluxN_{s}\phi_{S}+fluxV_{s}\\=\tau^{\phi}_{s}(\frac{S}{\delta y})_{s}\phi_{C}-\tau^{\phi}_{s}(\frac{S}{\delta y})_{s}\phi_{S}+0$

将上述结果代入

$J^{\phi}_{e}S_{e}\cdot\vec{i}-J^{\phi}_{w}S_{w}\cdot\vec{i}+J^{\phi}_{n}S_{n}\cdot\vec{j}-J^{\phi}_{s}S_{s}\cdot\vec{j}=Q^{\phi}_{C}V_{C}$

并整理为

$a_{C}\phi_{C}+a_{E}\phi_{E}+a_{W}\phi_{W}+a_{N}\phi_{N}+a_{S}\phi_{S}=b_{C}$

得，

$\\a_{E}=-\tau^{\phi}_{e}(\frac{S}{\delta x})_{e}\\a_{W}=-\tau^{\phi}_{w}(\frac{S}{\delta x})_{w}\\a_{N}=-\tau^{\phi}_{n}(\frac{S}{\delta y})_{n}\\a_{S}=-\tau^{\phi}_{s}(\frac{S}{\delta y})_{s}\\a_{C}=-(a_{E}+a_{W}+a_{N}+a_{S})\\b_{C}=Q^{\phi}_{C}V_{C}-(fluxV_{e}+fluxV_{w}+fluxV_{n}+fluxV_{s})$

可整理为通用形式

$\\a_{NEIGHBOR}=-\tau^{\phi}_{f}gdiff_{f}=-\tau^{\phi}_{f}\frac{S_f}{abs(d_{CN})}\\a_{OWNER}=-(a_{E}+a_{W}+a_{N}+a_{S})\\b_{C}=Q^{\phi}_{C}V_{C}-(fluxV_{e}+fluxV_{w}+fluxV_n+fluxV_s)$

合理的离散方程应反映原始守恒方程的特性。离散方程的系数应满足以下规则：

1）零和法则

单元C的参数值由源项及其相邻单元参数值决定，当无源项时，单元C的值仅由相邻单元参数值加权决定，因此单元C的参数值应在相邻单元参数值的范围内，否则参数分布是非物理的。这一思想可用零和法则描述，如下式所示

2）符号相反法则

当相邻单元的参数值增大/减小时，单元C的参数值也应增大/减小，因此aC与aF的符号应相反。

2. 边界条件

对于对流/扩散问题，边界条件类型包括Dirichlet, Neumann, 混合边条，对称边条。边条被应用于边界网格单元上，离散的参数Phi值存储于边界单元中心以及边界面上。

如上图所示，对边界单元C，其离散方程可写作

对内部界面n, s, w，通量计算方式与内部界面一致。对于边界界面b，其通量为

$\\J^{\phi, D}_b\cdot \vec{S}_{b}\\=-\tau^{\phi}_{b}\nabla \phi_{b}\cdot S_b\vec{i}\\=-\tau^{\phi}_bS_b(\frac{\partial \phi}{\partial x}\vec{i}+\frac{\partial \phi}{\partial y}\vec{j})\cdot \vec{i}\\=-\tau^{\phi}_bS_b\frac{\phi_{b}-\phi_{C}}{d_{bC}}\\=fluxC_b\phi_{C}+fluxN_b\phi_b+fluxV_b$

1) Dirichlet边界条件

给定边界上参数值，

$\phi_{b}=\phi_{specified}$

此时

$\\fluxC_b=\tau^{\phi}_b\frac{S_b}{d_{bC}}=a_b \\fluxN_b=-\tau^{\phi}_b\frac{S_b}{d_{bC}}=-a_b\\ fluxV_b=0$

将其代入通用积分代数式中

$a_{C}\phi_{C}+a_{E}\phi_{E}+a_{W}\phi_{W}+a_{N}\phi_{N}+a_{S}\phi_{S}=b_{C}$

有

$\\a_{E}=a_b=fluxN_b\\ a_W=fluxN_w=-\tau^{\phi}_w\frac{S_w}{d_{CW}}\\ a_N=fluxN_n=-\tau^{\phi}_n\frac{S_n}{d_{CN}}\\ a_S=fluxN_s=-\tau^{\phi}_s\frac{S_s}{d_{CS}}\\ a_C=fluxC_b+fluxC_w+fluxC_n+fluxC_s\\=-(a_b+a_W+a_N+a_S)\\b_{C}=Q^{\phi}_CV_C-fluxV_b-\sum_{f\sim w,s,n}fluxV_f$

由于 $a_E\phi_E=a_b\phi_b$ 已知，因此可以将其移至右项bc中，有

$b_{C}=Q^{\phi}_CV_C-\sum_{f\sim b,w,s,n}fluxV_f-a_b\phi_b$

因此代数方程组为

$a_{C}\phi_{C}+a_{W}\phi_{W}+a_{N}\phi_{N}+a_{S}\phi_{S}=b_{C}$

$\\ a_W=fluxN_w=-\tau^{\phi}_w\frac{S_w}{d_{CW}}\\ a_N=fluxN_n=-\tau^{\phi}_n\frac{S_n}{d_{CN}}\\ a_S=fluxN_s=-\tau^{\phi}_s\frac{S_s}{d_{CS}}\\ a_C=fluxC_b+fluxC_w+fluxC_n+fluxC_s\\=-(a_b+a_W+a_N+a_S)\\b_{C}=Q^{\phi}_CV_C-\sum_{f\sim b,w,s,n}fluxV_f-a_b\phi_b$

当其他条件相同时，由于C与界面b之间的距离小于C与W, S, N点之间的距离，因此ab大于aW等系数，即界面b处参数对C处参数的影响更大。

2）Von Neumann边界条件

直接给定边界处phi的通量：

$J^{\phi, D}_b\cdot \vec{i}=q_b$

因此

$J^{\phi, D}_b\cdot \vec{S}_b=J^{\phi, D}_b\cdot S_b\vec {i}=q_bS_b\\=flux_C\phi_C+flux_N\phi_b+fluxV_b$

其中fluxC=fluxN=0, fluxVb=qbSb

得代数方程组为

$a_{C}\phi_{C}+a_{W}\phi_{W}+a_{N}\phi_{N}+a_{S}\phi_{S}=b_{C}$

$\\ a_W=fluxN_w=-\tau^{\phi}_w\frac{S_w}{d_{CW}}\\ a_N=fluxN_n=-\tau^{\phi}_n\frac{S_n}{d_{CN}}\\ a_S=fluxN_s=-\tau^{\phi}_s\frac{S_s}{d_{CS}}\\ a_C=fluxC_w+fluxC_n+fluxC_s\\=-(a_W+a_N+a_S)\\b_{C}=Q^{\phi}_CV_C-\sum_{f\sim b,w,s,n}fluxV_f$

在Von Neumann条件下，单元中心C处参数值可以超过/低于其相邻单元参数值，这是符合物理实际的（有通量流入/流出单元）。

3）混合边界条件

混合边界条件如下图所示

此时边界处通量为

其中 $h_\infty$ 为传热系数， $\phi_\infty$ 为该边界外部参数值，将上式写作（边界处通量守恒）

$-(\tau^{\phi}\nabla \phi)_bS_b\cdot \vec{i}=-h_{\infty }(\phi_{\infty}-\phi_b)\Delta y_C\\ -\tau^{\phi}_b(\frac{\phi_b-\phi_C}{d_{bC}})S_b=-h_{\infty }(\phi_{\infty}-\phi_b)\Delta y_C$

可计算得

代入 $J^{\phi,D}_b\vec {S}_b=-h_{\infty }(\phi_{\infty }-\phi_b)S_b$ 中，有

令上式等于 $flux_C\phi_C+flux_N\phi_b+fluxV_b$ ，则有

$flux_C=R_{eq}, flux_N=0, fluxV_b=-R_{eq}\phi_{\infty }$

代数方程组为

$a_{C}\phi_{C}+a_{W}\phi_{W}+a_{N}\phi_{N}+a_{S}\phi_{S}=b_{C}$

4）对称边界条件

当边界为对称边界时，其法向参数通量为0，因此对称边界等价于一种Von Neumann边界条件的在通量为0时的特例（fluxCb=fluxNb=fluxVb=0）。令Von Neumann中qb=0，即可得对称边界条件下代数方程组

$a_{C}\phi_{C}+a_{W}\phi_{W}+a_{N}\phi_{N}+a_{S}\phi_{S}=b_{C}$

$\\ a_W=fluxN_w=-\tau^{\phi}_w\frac{S_w}{d_{CW}}\\ a_N=fluxN_n=-\tau^{\phi}_n\frac{S_n}{d_{CN}}\\ a_S=fluxN_s=-\tau^{\phi}_s\frac{S_s}{d_{CS}}\\ a_C=fluxC_w+fluxC_n+fluxC_s\\=-(a_W+a_N+a_S)\\b_{C}=Q^{\phi}_CV_C-\sum_{f\sim w,s,n}fluxV_f$

3. 界面扩散因子

当扩散因子随空间坐标变化时，其存储于单元中心E, W, ...处，在这种情况下需要预估扩散因子在界面上的值。

假设在单元中心C与其相邻单元中心之间扩散因子呈线性分布，则有

$\tau^{\phi}_{e}=(1-g_e)\tau^{\phi}_C+g_{e}\tau^{\phi}_E$

这种方法不适用于扩散系数急剧变化的情况。另一个更好方法的基本思想是：界面处扩散系数的精度影响较小，而界面处扩散通量的精度影响要大得多，该方法即为调和平均法

虽然调和平均仅对一维问题是精确的，但也可将其应用至多维情况。

4. 偏斜（Skewness）

界面上的参数值应等于该参数在该界面上的平均值。通常的做法是采用线性插值，预估参数在界面与两个单元中心点连线的交点处的值。当网格扭曲时，该交点f '不一定与界面中心点f重合，如下图所示。

为了保证离散方法的二阶精度，界面上积分点应选取在界面中心点f，因此需要定义从f '处参数值计算f处参数值的方法。泰勒展开法：

$\vec{d}_{f'f}$ 是由f '点指向f的向量。

5. 各项异性扩散

上述方法假设扩散是各向同性的，当扩散系数大小与方向相关时，该问题属于各向异性扩散。此时半离散化扩散方程为

此时扩散系数为二阶对称张量，有

等价于

此时只要令各向同性扩散系数为1，并将 $\vec{S}_f$ 替换为 ${\vec{S}_f}'$ 即可采用与各向同性扩散问题相同代码计算各向异性扩散问题。

6. 迭代过程松弛处理

由于问题的非线性，源项受当前迭代步参数值phi的影响，当迭代步之间phi值变化剧烈时，计算很可能发散，为了提高计算的收敛性，可通过亚松弛（under-relaxation）法限制迭代步之间phi值的变化大小。这里介绍亚松弛的一种方法。

已知通用代数方程为

$\phi^{*}_C$ 为上一步迭代结果，有

括号内计算结果代表两个相邻迭代步之间phi的变化大小，通过给该变化乘以一个松弛系数控制phi的变化

或

当系数 $\lambda ^{\phi}>1$ 时称为过松弛，当 $0<\lambda ^{\phi}<1$ 时称为亚松弛。在CFD计算中通常采用的是亚松弛。最优松弛系数的选取与算例相关，需要依据经验给定。影响最优松弛系数的因素包括：所求解问题，网格数大小，网格尺寸及膨胀比，所采用的迭代方法等。在整个计算域中不一定都采用相等的松弛系数，松弛系数也可能随着迭代的进行而变化。

7. uFVM

在uFVM中，内部单元界面上的扩散项离散是在方程cfdAssembleDiffusionTermInterior中进行的。扩散项是在cfdProcessOpenFoamMesh.m中设置的。并对边界条件进行处理，具体在文件cfdAssembleDiffusionTerm.m中。最后在cfdAssembleIntoGlobalMatrixFaceFluxes函数中生成求解矩阵。