5 交错网格与SIMPLE算法（末尾有彩蛋）

5 交错网格与SIMPLE算法（末尾有彩蛋）

提一嘴，这一节的深入讲解可在陶文铨《数值传热学 第二版》的第六章找到。

背景

对于二维流动下的动量输运，可分别写出x方向、y方向的动量方程，以及连续性方程。

\[\begin{equation} \frac{\partial(\rho\phi)}{\partial t}+\mathrm{div}(\rho\mathbf{U}\phi)=\mathrm{div}(\Gamma_\phi~\mathrm{grad}\phi)+S_\phi \end{equation} \]

x方向：

\[\begin{equation} \frac{\partial}{\partial x}(\rho uu)+\frac{\partial}{\partial y}(\rho vu)=\frac{\partial}{\partial x}\Big(\mu\frac{\partial u}{\partial x}\Big)+\frac{\partial}{\partial y}\Big(\mu\frac{\partial u}{\partial y}\Big)-\frac{\partial p}{\partial x}+S_u \end{equation} \]

y方向：

\[\begin{equation} \frac{\partial}{\partial x}(\rho uv)+\frac{\partial}{\partial y}(\rho vv)=\frac{\partial}{\partial x}\Big(\mu\frac{\partial v}{\partial x}\Big)+\frac{\partial}{\partial y}\Big(\mu\frac{\partial v}{\partial y}\Big)-\frac{\partial p}{\partial y}+S_v \end{equation} \]

连续性方程：

\[\begin{equation} \frac\partial{\partial x}(\rho u)+\frac\partial{\partial y}(\rho v)=0 \end{equation} \]

但在求解上述方程的时候，面临如下问题：

• 对流项是非线性的；
• 所有方程的速度和压力耦合：速度出现在每一个动量方程和连续性方程中，而压力梯度的贡献通常是先验未知的，等待求解。

而对于不可压缩流动，连续性方程不会给出任何压力与速度的关系。因此问题在于如何解决压力与速度之间的联系。

交错网格

若采用常规方法建立网格，把u，v和p都存在同一套网格的节点上，则在数值计算中会遇到以下问题：考虑一维流动，稳态时有

\[\begin{equation} \rho u\mathrm{~}\frac{\mathrm{d}u}{\mathrm{d}x}=-\frac{\mathrm{d}p}{\mathrm{d}x}+\eta\frac{\mathrm{d}^2u}{\mathrm{d}x^2} \end{equation} \]

​​

对于上图所示均分网格，将上式中各项取中心差分，有：

\[\begin{equation} \rho u_i\frac{u_{i+1}-u_{i-1}}{2\delta x}=-\frac{p_{i+1}-p_{i-1}}{2\delta x}+\eta_i\frac{u_{i+1}-2u_i+u_{i-1}}{(\delta x)^2} \end{equation} \]

该式表明，对中间\(i\)点的离散方程不包含在此点处的压力\(p_i\)，而是把被\(i\)点隔开的两邻点的压力联系了起来，称为\(2-\delta\)​压差。

这导致了一个问题：若在迭代求解过程的某一层次上，在压力场的当前值上，叠加一个锯齿状的压力波（隔一个网格节点来一个相同的误差值），上述方程是不能把这一不合理的分量检测出来的。

因此，在二维情况下，动量离散方程可能无法检测出如下图所示的 “棋盘型压力场”（checkerboard pressure field） ：

​​

为检测出这种每隔一个节点压力相等的不合理压力，需要让动量方程中，压力梯度的离散形式以相邻两点之间的压力差（或称为\(1-\delta\)​压差）表示。为此，我们采用交叉网格（staggered grid） 。

所谓交叉网格，就是把速度u存在压力p控制容积的东、西界面上，而把速度v存放在压力控制容积的南、北界面上，如下图所示。此时，u控制容积和压力控制容积在x方向、v控制容积和压力控制容积在y方向分别有半个网格步长的错位。

​​

在交错网格系统中，压力梯度的离散形式对\(u_e\)（东侧/右侧节点的x方向速度）为\(\frac{p_E-p_P}{(\delta x)_e}\)，对\(v_n\)为

\(\frac{p_N-p_P}{(\delta_y)_n}\)，也就是相邻两点间的压力差构成了\(\frac{\partial p}{\partial x}\)，\(\frac{\partial p}{\partial y}\)，解决了一般网格系统遇到的困难。

特别注意：对于u和v的控制容积，边界控制容积为0，和边界相邻的内部控制容积为\(1.5\Delta x \Delta y\)（参见上图b，u的控制容积超出了东侧控制容积的中心，触碰到边界），其他内部控制容积为\(\Delta x \Delta y\)。

动量方程

速度分量的交错位置决定了其对应控制体的位置。这些控制体的面垂直于速度分量的方向，并且这些面会通过网格点。

​​

例如，对于上图所示的交错网格，对于围绕速度\(u_w\)的控制容积，它的两个垂直于x轴的面会穿过P点和W点。

考虑此时的动量方程，对于该图左侧的交错网格：

\[\begin{equation} a_wu_w=\sum a_{nb}u_{nb}+(P_W-P_P)A_w+b \end{equation} \]

该公式各项的含义为：

• 左1：西侧面w上，x方向的动量，表示净流入原控制体的动量通量
• 右1：周围网格点的速度，通过对流和扩散作用，对现中心网格点速度\(u_w\)的影响。这里的系数\(a_{nb}\)取决于采用的离散格式，可参考《4 中心差分和迎风差分》。
  这里的下标\(nb\)表示neighboring，即邻近网格点。
• 右2：压力梯度项，表示西侧点W和中间点P之间的压力差，对\(u_w\)的影响。压力差产生一个力，推动流体从高压区流向低压区。
• 右3：源项。

该公式的推导过程概述如下：

1. 从不可压缩流体的二维N-S方程，当中的x方向动量方程出发
2. 对动量方程，在围绕\(u_w\)的控制体上积分
3. 对时间和空间导数离散化，例如使用中心差分、迎风差分，然后整理公式。

同样地，对于南侧交错网格，也有如下公式：

\[\begin{equation} a_sv_s=\sum a_{nb}v_{nb}+(P_S-P_P)A_s+b \end{equation} \]

由于采用了交错网格，穿过控制体积的面的扩散导和通量需要重新计算。

如，计算\(u_w\)东侧控制面的通量（这里回顾第4讲的对流质量通量F和扩散传导系数D）：

\[\begin{equation} F_{east}=(\rho u)_{east}=\frac{\rho_wu_w+\rho_eu_e}{2} \end{equation} \]

由于采用了交错网格，此时\(u_w\)控制体积的东侧面位于P点，然后用上述分式计算。

\[\begin{equation} D_{east}=\frac{\Gamma_P}{\delta x_{wP}+\delta x_{Pe}} \end{equation} \]

这个实际上就是计算P点处的扩散传到了，注意特征长度\(\delta\)的不同。

SIMPLE算法

压力与速度的求解

SIMPLE算法的全称是：压力关联方程的半隐式方法（semi-implicit method for pressure-linked equations）。

在使用该算法之前，需要先假定一组速度\(u_0\)、\(v_0\)，把动量离散方程中的系数和常数项求出来。

该算法的第一步，是假定猜想一个压力场，记为\(p^*\)，然后用它求解上面的动量离散方程，得到一组速度\(u^*\)，\(v^*\)。

为讨论方便，我们先从已知压力后求速度开始。

为得到满足连续性方程的速度场，上述星标速度需要被修正，这是通过对估计压力加上修正值\(P'\)达成的：

\[\begin{equation} \begin{aligned} P = P^* + P^{\prime} \\ u = u^* + u^{\prime} \\ v = v^* + v^{\prime} \end{aligned} \end{equation} \]

则动量方程变成：

\[\begin{equation} \begin{aligned} a_e(u_e^*+u_e^{\prime})=\sum a_{nb}(u_{nb}^*+u_{nb}^{\prime})+b\\+[(p_P^*+p_P^{\prime})-(p_E^*+p_E^{\prime})]A_e \end{aligned} \end{equation} \]

（注意压力的角标，这里我参考的是陶文铨院士的《数值传热学》，考虑的是东侧、北侧的下游节点，流向还是从西向东、从南向北，总之是流向上游的减去流向下游的）

我们假定由源项构成的b的值保持不变，则有：

\[\begin{equation} a_{e}u_{e}^{\prime}=\sum a_{nb}u_{nb}^{\prime}+(p_{P}^{\prime}-p_{E}^{\prime})A_{e} \end{equation} \]

这里，我们近似地忽略四周邻点速度修正值的影响，即令系数\(a_{nb}=0\)（这就是半隐式方法的含义，保留这部分就是全隐了），则有：

\[\begin{equation} \begin{aligned} &a_{e}u_{e}^{^{\prime}}=(p_{P}^{^{\prime}}-p_{E}^{^{\prime}})A_{e}\\ &u_{e}^{^{\prime}}=\left(\frac{A_{e}}{a_{e}}\right)(p_{P}^{^{\prime}}-p_{E}^{^{\prime}})=d_{e}(p_{P}^{^{\prime}}-p_{E}^{^{\prime}})\\ &v_{n}^{\prime}=d_{n}(p_{P}^{\prime}-p_{N}^{\prime}),~d_{n}=\frac{A_{n}}{a_{n}} \end{aligned} \end{equation} \]

这样获得改进后的速度：

\[\begin{equation} u_{e}=u_{e}^{\star}+d_{e}(p_{P}^{\prime}-p_{E}^{\prime}),~v_{n}=v_{n}^{\star}+d_{n}(p_{P}^{\prime}-p_{N}^{\prime}) \end{equation} \]

下一步，我们讨论确定压力修正值的方法。

回顾连续性方程（式4），在时间间隔\(\Delta t\)内，对主控制体（没交错网格的）做积分，采用全隐格式，得：

\[\begin{equation} [(\rho u)_{e}-(\rho u)_{w}]\Delta y+[(\rho v)_{n}-(\rho v)_{s}]\Delta x=0 \end{equation} \]

对于一个控制体积，连续性方程变成：

\[\begin{equation} (\rho uA)_w-(\rho uA)_e+(\rho vA)_s-(\rho vA)_n=0 \end{equation} \]

把改进后的速度回代进去，整理得到：

\[\begin{equation} a_Pp_P^{\prime}=a_Ep_E^{\prime}+a_Wp_W^{\prime}+a_Np_N^{\prime}+a_Sp_S^{\prime}+b_P^{\prime} \end{equation} \]

其中：

\[a_E = (\rho u A)_e，a_W = (\rho u A)_w，a_N = (\rho u A)_n，a_S = (\rho u A)_s \]

\[\begin{equation} \begin{aligned} a_{P}&=a_{E}+a_{W}+a_{N}+a_{S}\\ b_{P}^{\prime}&=(\rho u^{*}A)_{w}-(\rho u^{*}A)_{e}+(\rho v^{*}A)_{s}-(\rho v^{*}A)_{n} \end{aligned} \end{equation} \]

这里，源项\(b_P'\)表示由不正确的速度场\(u^*\)，\(v^*\)导致的连续性不平衡，或者说，是一个控制容积不满足连续性的，剩余质量的大小。

当速度场迭代收敛时，各控制容积b的绝对值最大值，以及控制容积的b值代数和都应为最小量。

此时，修正值u'和v'将为0，压力修正为0，\(p^{*}=p,~u^{*}=u,~v^{*}=v\)。

例子

​​

已知：\(p_W = 60\)，\(p_S = 50\)，\(u_e=20\)，\(v_n=7\)。

又给定：\(u_w=0.7(p_W-p_P)\)，\(v_s=0.6(p_S-p_P)\)。

求解\(p_P\)，\(u_w\)和\(v_s\)。

解：假设\(p_P=20\)，则：

\(u_w^*=0.7(60-20)=28\)，\(v_s^*=12\)。

设连续性方程：\(u_w+v_s=u_e+v_n\)，则使用SIMPLE算法，\(u_w\)和\(v_s\)表示为：

\(u_{w}=u_{w}^{\star}+d_{w}(p_{W}^{\prime}-p_{P}^{\prime})\)

\(v_s=v_s^*+d_s(p'_S-p'_P)\)​

按已知条件，由于\(p_W\)和\(p_S\)已知，有\(p'_W=p'_S=0\)，则：

\(u_w=28-0.7p'_P，v_s=12-0.6p'_P\)

代入连续性方程，得到关于\(p'_P\)的方程：

\(40-1.3p'_P=27\)

解得\(p'_P=10，p_P=30\)。

回代入速度的方程，得：

\(u_w=u_w^* -0.7p'_P=28-7=21\)，\(v_s=6\)

亚松弛

对于压力，由于在速度修正值公式中，略去了邻点的影响，用p'解得的修正速度是合适的，但压力修正值p'本身被夸大了，因而要进行亚松弛处理。

\[\begin{equation} p=p^{\star}+\alpha_{p}p^{\prime} \end{equation} \]

这里，\(\alpha_p\)（alpha，p小写!!!!!）是压力亚松弛因子。

对于速度的亚松弛，是把本层次计算结果，与上一层次结果的差值做适当缩减，避免由于差值过大，引起非线性迭代过程中的发散。

当然，根据杨教授的讲义，我们在这里考虑的不是“上一层次的结果”，而是同一迭代中星标估计值和实际值的差异。但数值传热学书中的定义可能不太一样，欢迎讨论。

根据讲义，我们对速度的动量方程略作修改：

\[\begin{equation} u_w=u_w^*+\alpha_u\left[\frac{\sum a_{nb}u_{nb}+(p_W-p_P)A_w+b}{a_w}-u_w^*\right] \end{equation} \]

其中\(\alpha_u\)是速度的亚松弛因子。

重排一下得：

\[\begin{equation} \frac{a_w}{\alpha_u}u_w=\sum a_{nb}u_{nb}+(p_W-p_P)A_w+b+\frac{1-\alpha_u}{\alpha_u}a_wu_w^* \end{equation} \]

改进SIMPLE算法

SIMPLER

这个算法的英文叫SIMPLE Revised。其思想为：

1. 假定了速度分布之后，和这一速度分布相协调的压力场，可由动量方程计算得到，无需另行单独假定一个压力场；
2. 只把p'用于修正速度，另外寻找更适合修正压力场的方法。

上述思想的结合就构成了SIMPLER算法。

我们从已知速度分布，计算压力场开始。动量离散方程可以写作：

\[\begin{equation} u_{e}=\underbrace{\frac{\sum a_{nb}u_{nb}+b}{a_{e}}}_{\hat{u}_{e}}+d_{e}(p_{P}-p_{E}) \end{equation} \]

在事先假定了速度分布以后，上式右边第一项即可算出，这一项也具有速度的量纲，称为假拟速度（pseudo velocity），记作\(\hat{u}_{e}\)。

这样，动量离散方程可简记作：

\[\begin{equation} u_{e}=\hat{u}_{e}+d_{e}(p_{P}-p_{E}) \end{equation} \]

\[\begin{equation} v_n=\hat{v}_n+d_n(p_P-p_N) \end{equation} \]

回代到离散连续性方程中，得：

\[\begin{equation} a_Pp_P=a_Ep_E+a_Wp_W+a_Np_N+a_Sp_S+b_P \end{equation} \]

这里，\(a_E = (\rho u A)_e，a_W = (\rho u A)_w，a_N = (\rho u A)_n，a_S = (\rho u A)_s\)，和SIMPLE是一样的。

\(a_p\)等于上述四个a的和也是一样的。

不同的地方在于\(b'_p\)的计算（当然这里换成了\(b_P\)）：

\[\begin{equation} b_P=(\rho\widehat{u}A)_w-(\rho\widehat{u}A)_e+(\rho\widehat{u}A)_s-(\rho\widehat{u}A)_n \end{equation} \]

SIMPLEC

在求解速度修正值时，我们忽略了\(\sum a_{nb}u_{nb}^{\prime}\)项，但这使得方程不协调。

为了能略去该项，同时使得方程基本协调，在速度修正值的两边同时减去\(\sum a_{nb}u_{e}^{\prime}\)，即：

\[\begin{equation} (a_e-\sum a_{nb})u_e^{\prime}=\sum a_{nb}(u_{nb}^{\prime}-u_e^{\prime})+A_e(p_P^{\prime}-p_E^{\prime}) \end{equation} \]

这里，预期\(u'_e\)和其邻点的修正值\(u'_{nb}\)量级相同，因此现在略去右一项产生的影响比原来要小得多。于是得：

\[\begin{equation} \begin{aligned} u'_{e}&=d_e(p^{_{P}^{\prime}}-p^{_{E}^{\prime}})=\frac{A_{e}}{a_{e}-\sum a_{nb}}(p^{_{P}^{\prime}}-p^{_{E}^{\prime}})\\ v'_n&=d_n(p^{_{P}^{\prime}}-p^{_{N}^{\prime}})=\frac{A_{n}}{a_{n}-\sum a_{nb}}(p^{_{P}^{\prime}}-p^{_{N}^{\prime}}) \end{aligned} \end{equation} \]

SIMPLEC和SIMPLE的主要区别：

• 用\(\frac{A_{e}}{a_{e}-\sum a_{nb}}\)代替SIMPLE中的\(\frac{A_e}{a_e}\)；
• 压力修正项p'不再亚松弛，即\(a_p=1\)。

总结

压力修正算法的共同特征有：

• 非线性：动量方程中，未知数以非线性的方式出现
• 耦合：动量、能量、组分方程等不同的输运方程之间相互耦合，一个方程的解会影响其他方程的解，如速度场影响温度场、温度场通过密度变化影响速度场
• 迭代求解：由于非线性和耦合，上述方程组无法直接求解。
• 交错网格：速度分量定义在交错的网格上，以避免棋盘压力场问题
• 速度存储在标量（压力、温度等）的控制体上
• 离散的动量方程在交错的控制台（速度控制体）上求解

SIMPLE算法：

• 一种用于计算压力场和速度场的迭代算法

• 首先猜测一个初始压力场\(p^*\)，可以不准确，迭代过程中会逐步修正

• 使用猜测的压力场，求解离散化的动量方程，得到中间速度场\(u^*\)，\(v^*\)。由于压力场是猜测的，速度场通常不满足连续性方程（即质量守恒）

• 把连续性方程转化为压力修正方程，用压力修正值修正速度场：

  \[\begin{equation} \begin{aligned}&p_{P}=p_{P}^{*}+p_{P}^{\prime}\\&u_{w}=u_{w}^{*}+d_{w}(p_{W}^{\prime}-p_{P}^{\prime})\\&v_{s}=v_{s}^{*}+d_{s}(p_{S}^{\prime}-p_{P}^{\prime})\end{aligned} \end{equation} \]

• 求解其他标量的离散输运方程

• 重复直到压力场、速度场、其他标量场都收敛，即迭代结果不再发生明显变化。

SIMPLE的改进：

• SIMPLER假定了速度分布之后，由动量方程计算得和这一速度分布相协调的压力场，只把p'用于修正速度；
• SIMPLEC用\(\frac{A_{e}}{a_{e}-\sum a_{nb}}\)代替SIMPLE中的\(\frac{A_e}{a_e}\)；
• 亚松弛缩减变量。

写完啦！

到此，2024-25学年的Computer Modeling Techniques课程的笔记就基本完成了。

数值方法和湍流是9月10月上的，那时候我还是尽量紧跟每一堂课的进度。

到中心差分就不行了，赶毕设进度，还有不少作业要做，而且上课能听懂1个小时，后面俩小时懵了。

所以先更新的结构力学部分，做做那部分的作业；然后拖到期末考试之前再写4和5。

由于二维平面应力是选做，我完成这份笔记的日期距离考试只有3天了，因此不再考虑。有兴趣的同学可以自行写一写，然后贴上来。

这个学期真的挺难的，毕设、大作业、各种杂务的压力都很大，同时还要追课堂进度，可以说是高三模式了0v0。

坚持下来的每一位朋侑都不容易，请允许我在此向我们送上敬意和祝福。

祝我们未来前程似锦，好运连连。

​​