FVM in CFD 学习笔记_第15章_流动计算:不可压缩流动_1_交错网格上的SIMPLE算法

2020年06月10日 阅读数:143
这篇文章主要向大家介绍FVM in CFD 学习笔记_第15章_流动计算:不可压缩流动_1_交错网格上的SIMPLE算法,主要内容包括基础应用、实用技巧、原理机制等方面,希望对大家有所帮助。

学习自F. Moukalled, L. Mangani, M. Darwish所著The Finite Volume Method in Computational Fluid Dynamics - An Advanced Introduction with OpenFOAM and Matlab
Chapter 15 Fluid Flow Computation: Incompressible Flowshtml

前面几章讲解的关于变量 ϕ \phi 的通常输运方程的离散和求解流程,均是创建在事先已知速度场的前提下。可是通常状况下,速度场是未知的,且须要经过求解Navier-Stokes方程组来获取。对于不可压缩流动而言,该项工做尤其复杂,由于压力和速度是强烈耦合在一块儿的,并且压力并不以主变量的形式出如今动量或是连续方程中。本章的重点在于展现解决上述两个问题的方法,以及不可压缩流动问题的流场计算方法。先讲解一维交错网格,而后是一维同位网格,最后是同位三维非结构网格。除了阐明SIMPLE、SIMPLEC、PRIME和PISO算法的前因后果,还清晰阐明了Rhie-Chow插值方法,以及将其扩展到瞬态、松弛和体积力项的方法。最后,展示了一些常见边界条件的添加细节。web

因为本章内容繁杂,篇幅较长,故分红了四部分来说解,各部分主要内容分别为:交错网格、同位网格、边界条件、SIMPLE家族算法。算法

这里是第一部分,主要讲解不可压缩流动问题的求解在常规网格上所碰到的问题,以及,交错网格上的SIMPLE算法是如何解决该问题的。网络

1 主要困难

前面几章所处理的通用守恒方程,能够改写成与连续方程或动量方程类似的形式。然而,目前为止所呈现的数值技术,并不足以求解Navier-Stokes方程,由于求解通常流动的算法要可以处理压力速度的耦合问题。为理解该问题,先把连续方程和动量方程写出来
ρ t + ( ρ v ) = 0 t ( ρ v ) + ( ρ v v ) = p + { μ [ v + ( v ) T ] } + f b \begin{aligned} \frac{\partial\rho}{\partial t}+\nabla\cdot(\rho\bold v)&=0 \\ \frac{\partial}{\partial t}(\rho \bold v)+\nabla\cdot(\rho\bold v\bold v)&= -\nabla p+\nabla\cdot\left\{ \mu\left[ \nabla\bold v+(\nabla\bold v)^\text{T} \right] \right\} + \bold f_b \end{aligned}
这俩方程是非线性的,但也并非不可克服的困难,由于这种问题一般可采用迭代方法来应对。此外,第二式是矢量方程,若把其展开成份量形式的话,则是个标量方程系统,须要顺序求解。再者,应力张量可改写成相似扩散项的形式,并可隐式处理,其第2项(即,速度梯度的转置 ( v ) T (\nabla\bold v)^\text{T} )则基于前次迭代值作显式估算并添加到源项中去。针对上述连续方程和动量方程,从通有标量方程的数值方法(前面几章讲解的方法)中没法汲取灵感来解决的主要难题在于,出如今动量方程中的压力场没法得到其显式计算方程。app

从这俩方程中能够看出,纵然速度场能够直接使用动量方程来计算,然而出如今动量方程中的压力场没法直接从连续方程中计算。该耦合关系(压力和速度的耦合)是如此强烈却又如此隐含,为便于理解,可把方程组写成以下矩阵形式
A u = [ F B T B 0 ] [ v p ] = [ f b 0 ] \bold A \bold u=\begin{bmatrix} \bold F & \bold B^\text T \\ \bold B & \bold 0 \end{bmatrix} \begin{bmatrix} \bold v \\ \bold p \end{bmatrix}= \begin{bmatrix} \bold f_b \\ \bold 0 \end{bmatrix}
在该形式中,上式中出现了一个零对角块,这是鞍点(saddle point)问题的特性,代表其没法经过任何迭代手法来求解压力和速度场。所以,须要推导出关于压力的方程才好。svg

一种处理方法是,经过把矩阵 A \bold A 分解成下三角 L \bold L 和上三角 U \bold U 矩阵,将动量和连续方程系统从新改写成下述形式
A = [ F B T B 0 ] = [ F 0 B B F 1 B T ] [ I F 1 B T 0 I ] = L U \bold A = \begin{bmatrix} \bold F & \bold B^\text T \\ \bold B & \bold 0 \end{bmatrix}= \begin{bmatrix} \bold F & \bold 0 \\ \bold B & -\bold B\bold F^{-1}\bold B^\text T \end{bmatrix} \begin{bmatrix} \bold I & \bold F^{-1}\bold B^\text T \\ \bold 0 & \bold I \end{bmatrix}= \bold L \bold U
其中 B F 1 B T -\bold B\bold F^{-1}\bold B^\text T 为Schur补矩阵(Schur complement matrix)。函数

这正是下面要讲的迭代求解Navier-Stokes方程的方法的本质,Patankar和Spalding将该技术嵌入在经典的分离式SIMPLE(Semi Implicit Method for Pressure Linked Equations)算法中。学习

求解流程是,将Navier-Stokes方程从新整理成动量方程和压力方程的形式,而后作离散和顺序求解。压力方程是经过联合半离散动量方程和连续方程(Schur补矩阵)来建立的。ui

算法是由Picard类型的迭代流程驱动的,其率先求解动量方程,使用的是前次迭代的压力场。求出的速度场是知足动量方程的可是未必知足质量方程。接下来用该速度场来构建压力方程,并用所得解来修正压力和速度场以便知足质量守恒。而后开始新的迭代过程,并重复上述流程,直到速度和压力场既知足质量守恒又知足动量守恒为止。spa

该算法可用以下的矩阵形式来描述
[ I D 1 B T 0 I ] [ v p ] = [ v p ] \begin{bmatrix} \bold I & \bold D^{-1}\bold B^\text T \\ \bold 0 & \bold I \end{bmatrix}\begin{bmatrix} \bold v \\ \bold p \end{bmatrix}= \begin{bmatrix} \bold v^* \\ \bold p^* \end{bmatrix}
接下来使用下式来更新速度场
[ F 0 B B D 1 B T ] [ v p ] = [ f b 0 ] \begin{bmatrix} \bold F & \bold 0 \\ \bold B & -\bold B\bold D^{-1}\bold B^\text T \end{bmatrix} \begin{bmatrix} \bold v^* \\ \bold p^* \end{bmatrix}= \begin{bmatrix} \bold f_b \\ \bold 0 \end{bmatrix}
上式和上上式中的 F 1 \bold F^{-1} 用其逆对角阵 D 1 \bold D^{-1} 来近似,上标 ^* 表示当前迭代的中间值。流程汇总以下:

  • 求解: F v = f b \bold F \bold v^*=\bold f_b
  • 求解: B D 1 B T p = B v -\bold B \bold D^{-1} \bold B^\text T \bold p^*=-\bold B \bold v^*
  • 更新: v = v D 1 B T p \bold v=\bold v^*-\bold D^{-1}\bold B^\text T\bold p^*
  • 更新: p = p \bold p=\bold p^*

这种类型的分裂与SIMPLE系列算法中所用的分裂是相同的,这正是本章的主题。

2 初步推导

在这里插入图片描述

针对不可压缩流动问题发展求解算法的困难之处,可经过对上图所示的一维空间均匀网格问题的离散处理来加深理解。为简便起见,假设流动是定常的。简化的连续和动量方程(写成守恒形式)为
( ρ u ) x = 0   ( ρ u u ) x = x ( μ u x ) p x \begin{aligned} \frac{\partial (\rho u)}{\partial x}&=0 \\ ~ \\ \frac{\partial(\rho u u)}{\partial x}&=\frac{\partial}{\partial x} \left(\mu \frac{\partial u}{\partial x} \right) -\frac{\partial p}{\partial x} \end{aligned}

2.1 动量方程离散

动量方程的离散,首先要对其在单元 C C 上作积分,得
V C ( ρ u u ) x d V = V C x ( μ u x ) d V V C p x d V \int_{V_C}\frac{\partial(\rho u u)}{\partial x}dV= \int_{V_C}\frac{\partial}{\partial x} \left(\mu \frac{\partial u}{\partial x} \right)dV- \int_{V_C}\frac{\partial p}{\partial x}dV
使用散度定理,将对流项和扩散项的体积分转化成面积分
V C ( ρ u u ) d y = V C μ u x d y V C p x d V \int_{\partial V_C}(\rho u u)dy=\int_{\partial V_C}\mu \frac{\partial u}{\partial x}dy - \int_{V_C}\frac{\partial p}{\partial x}dV
把面积分用流过单元面的通量加和来替换,并使用单Gauss点来作面积分,获得半离散形式为
( ρ u Δ y ) e m ˙ e u e + ( ρ u Δ y ) w m ˙ w u w = ( μ u x Δ y ) e ( μ u x Δ y ) w V C p x d V \underbrace{(\rho u \Delta y)_e}_{\dot m_e} u_e+ \underbrace{-(\rho u \Delta y)_w}_{\dot m_w} u_w= \left( \mu \frac{\partial u}{\partial x}\Delta y \right)_e- \left( \mu \frac{\partial u}{\partial x}\Delta y \right)_w- \int_{V_C}\frac{\partial p}{\partial x}dV
其可重写为
m ˙ e u e + m ˙ w u w C o n v e c t i o n [ ( μ u x Δ y ) e ( μ u x Δ y ) w ] D i f f u s i o n = V C p x d V \underbrace{\dot m_e u_e+\dot m_w u_w}_{Convection}- \underbrace{\left[\left( \mu \frac{\partial u}{\partial x}\Delta y \right)_e- \left( \mu \frac{\partial u}{\partial x}\Delta y \right)_w\right]}_{Diffusion}=- \int_{V_C}\frac{\partial p}{\partial x}dV
对流和扩散项可用以前章节中所讲的方法来离散化处理,获得以下形式的代数方程
a C u u C + F N B ( C ) ( a F u u F ) = b C u V C p x d V a_C^u u_C+\sum_{F\sim NB(C)}(a_F^u u_F)=b_C^u-\int_{V_C}\frac{\partial p}{\partial x}dV
压力项的离散待我们讲完连续方程的离散后再说。

2.2 连续方程的离散

连续方程的离散,也是要对其在单元 C C 上作积分,得
V C ( ρ u ) x d V = 0 \int_{V_C}\frac{\partial(\rho u)}{\partial x}dV=0
再次使用散度定理将体积分转化成面积分,而后转化成单元面通量的加和形式,连续方程的离散形式为
f n b ( C ) ( ρ u Δ y ) f = ( ρ u Δ y ) e ( ρ u Δ y ) w = 0 \sum_{f\sim nb(C)}(\rho u \Delta y)_f=(\rho u \Delta y)_e - (\rho u \Delta y)_w=0

f n b ( C ) m ˙ f = m ˙ e + m ˙ w = 0 \sum_{f\sim nb(C)}\dot m_f=\dot m_e + \dot m_w=0

2.3 棋盘问题(Checkerboard Problem)

压力项离散可由下面两种方法中的任意一个来实现。第一种方法,经过单点Gauss积分来计算体积分
V C p x d V = ( p x ) C V C \int_{V_C}\frac{\partial p}{\partial x}dV=\left(\frac{\partial p}{\partial x}\right)_CV_C
使用中心差分格式,上式的离散形式为
V C p x d V = p E p W 2 Δ x V C \int_{V_C}\frac{\partial p}{\partial x}dV=\frac{p_E-p_W}{2\Delta x} V_C
第二种方法,把压力梯度项的体积分转化成面积分
V C p x d V = V C p d y \int_{V_C}\frac{\partial p}{\partial x}dV=\int_{\partial V_C}p dy
把面积分写成单元面通量加和的形式,有
V C p x d V = V C p d y = p e Δ y e p w Δ y w = ( p e p w ) Δ y = ( p e p w ) V C Δ x \int_{V_C}\frac{\partial p}{\partial x}dV=\int_{\partial V_C}p dy= p_e\Delta y_e - p_w \Delta y_w=(p_e-p_w)\Delta y= (p_e-p_w)\frac{V_C}{\Delta x}
压力的变化选用线性插值廓线,把压力梯度项从新写成主节点上压力值的函数形式
V C p x d V = [ 1 2 ( p E + p C ) 1 2 ( p C + p W ) ] V C Δ x = p E p W 2 Δ x V C \int_{V_C}\frac{\partial p}{\partial x}dV= \left[\frac12(p_E+p_C)-\frac12(p_C+p_W)\right]\frac{V_C}{\Delta x}= \frac{p_E-p_W}{2\Delta x} V_C
两个方法导出了一样的表达式,即,包含交替节点 E E W W 之间压力差值的表达式(这里的交替指的是, E E W W 节点之间间隔了一个节点 C C )。

采用一样的方法,使用线性插值廓线,并注意到密度是常数,且有 ( Δ y ) e = ( Δ y ) w = ( Δ y ) C (\Delta y)_e=(\Delta y)_w=(\Delta y)_C ,则连续方程为
u E u W = 0 u_E-u_W=0
用的也是两个交替网格节点的速度值。

在上上式中,单元 C C 处的压力梯度项是由两个交替的(而非连续的)横跨单元两侧的网格节点上的压力值所决定的。连续方程亦是如此,这样一来就只是对交替的速度单元施加了守恒特性。这意味着不符合物理意义的之字形(或棋盘型,国际象棋中的黑白格子间隔排列)的压力场或速度场(以下图所示的那样),在该数值方法中将会被感知成均匀场。换句话说,像这种…-A-B-A-B-…排列的场,本来并不是是均匀场,可是因为数值算法中的压力梯度和连续方程均考虑的是间隔单元间的特性,因此反却是把其错误滴认定为是均匀场了。

在这里插入图片描述

针对上图所示的压力和速度值,在点 W W C C E E 处的压力梯度为
V W p x d V = ( p C p W W ) V W 2 Δ x W = ( 10 10 ) V W 2 Δ x W = 0 V C p x d V = ( p E p W ) V C 2 Δ x C = ( 100 + 100 ) V C 2 Δ x C = 0 V E p x d V = ( p E E p C ) V E 2 Δ x E = ( 10 10 ) V E 2 Δ x E = 0 \begin{aligned} &\int_{V_W}\frac{\partial p}{\partial x}dV=(p_C-p_{WW})\frac{V_W}{2\Delta x_W}=(10-10)\frac{V_W}{2\Delta x_W}=0 \\ &\int_{V_C}\frac{\partial p}{\partial x}dV=(p_E-p_{W})\frac{V_C}{2\Delta x_C}=(-100+100)\frac{V_C}{2\Delta x_C}=0 \\ &\int_{V_E}\frac{\partial p}{\partial x}dV=(p_{EE}-p_{C})\frac{V_E}{2\Delta x_E}=(10-10)\frac{V_E}{2\Delta x_E}=0 \\ \end{aligned}
并且连续方程看起来在每一个单元上都是知足的,由于
V W u x d V = ( u C u W W ) V W 2 Δ x W = ( 1 1 ) V W 2 Δ x W = 0 V C u x d V = ( u E u W ) V C 2 Δ x C = ( 10 10 ) V C 2 Δ x C = 0 V E u x d V = ( u E E u C ) V E 2 Δ x E = ( 1 1 ) V E 2 Δ x E = 0 \begin{aligned} & \int_{V_W}\frac{\partial u}{\partial x}dV=(u_C-u_{WW})\frac{V_W}{2\Delta x_W}=(1-1)\frac{V_W}{2\Delta x_W}=0 \\ & \int_{V_C}\frac{\partial u}{\partial x}dV=(u_E-u_{W})\frac{V_C}{2\Delta x_C}=(10-10)\frac{V_C}{2\Delta x_C}=0 \\ & \int_{V_E}\frac{\partial u}{\partial x}dV=(u_{EE}-u_{C})\frac{V_E}{2\Delta x_E}=(1-1)\frac{V_E}{2\Delta x_E}=0 \\ \end{aligned}
在多维状况下,也会产生类似的非物理特性,即便很难具象化展现(其实若是是二维结构网格的话,这种状况跟国际象棋的黑白格子就很像了,即,黑格子是一个值,而白格子是另外一个值,可是按照上面这种数值方法处理起来,整场各个单元反却是知足了连续方程,并且压力梯度也都是零了)。这为下面展现一种处理该问题的方法作好了铺垫。

2.4 交错网格

前面公式中所呈现的问题是压力和速度场的解耦问题。能够经过把不一样的变量存储在交错位置上来增强它们的耦合特性,由于这样一来,就不须要用插值手段来计算动量方程中的压力梯度以及连续方程中的速度场了。这样的交错网格以下图所示,在交错网格中速度场是存储在单元面上的(图a),而压力场则是存储在单元形心上的(图b)。
在这里插入图片描述

单元 C C 的连续方程离散为
f n b ( C ) m ˙ f = m ˙ e + m ˙ w = 0 \sum_{f\sim nb(C)}\dot m_f=\dot m_e+\dot m_w=0

u e u w = 0 u_e-u_w=0
无需再作插值了,由于速度在界面 e e w w 处是能够得到的。此外,动量方程的积分是在单元 e e 上进行的,其产生以下的离散动量方程形式
a e u u e + f N B ( e ) a f u u f = b e u V e ( p ) e = b e u V e p E p C δ x e a_e^u u_e+\sum_{f\sim NB(e)}a_f^u u_f=b_e^u-V_e(\nabla p)_e= b_e^u-V_e\frac{p_E-p_C}{\delta x_e}
压力梯度的计算用的是横跨单元面两侧的连续单元形心节点值,故不须要作插值。所以棋盘压力场和速度场的状况将不会出现,由于这些状况将很容易地被数值方法探测并衰减掉。

2.5 压力修正方程

接下来要呈现的推导是基于Patankar和Spalding的工做,他们发展了SIMPLE(Semi Implicit Method for Pressure Linked Equations)算法的最初实用形式。

由上图所示网格上的离散连续方程和动量方程出发
f n b ( C ) m ˙ f = 0 a e u u e + f N B ( e ) a f u u f = b e u V e ( p x ) e \begin{aligned} \sum_{f\sim nb(C)}\dot m_f &= 0 \\ a_e^u u_e+\sum_{f\sim NB(e)}a_f^u u_f&= b_e^u-V_e\left(\frac{\partial p}{\partial x}\right)_e \end{aligned}
求解过程首先须要提供初始的猜想速度和压力场。把初始猜想的或是迭代开始时的解用上标 ( n ) ^{(n)} 表示,那么这些个暂定的(初始猜想的或是未达到收敛的不肯定的)速度和压力场用 u ( n ) u^{(n)} p ( n ) p^{(n)} 表示。在每步迭代中,先求解动量方程来获取速度场,所得解用上标 ^* 表示,由于其并不是当前迭代步的最终解(注意,因为并未考虑连续方程,因此这时候的速度场并不知足连续方程)。这样,动量方程知足
a e u u e + f N B ( e ) a f u u f = b e u V e ( p ( n ) x ) e a_e^u u_e^*+\sum_{f\sim NB(e)}a_f^u u_f^*= b_e^u-V_e\left(\frac{\partial p^{(n)}}{\partial x}\right)_e
其中压力场仍旧采用的是前次迭代值,那么,算出来的速度场 u u^* 知足动量方程,可是未必知足连续方程,由于压力场并非精确值。所以,须要寻找修正措施来确保速度(或质量流量)和压力知足连续方程。

将修正场用上标撇来表示,即,( u u' p p' ),那么速度和压力修正为
u = u + u p = p + p u=u^*+u' \\ p=p^*+p'
注意,单元面处的质量流量也修正为
m ˙ f = m ˙ f + ρ u S f x = m ˙ f + m f \dot m_f=\dot m_f^*+\rho u' S_f^x=\dot m_f^*+m_f'
这样,精确质量流量应知足连续方程,即
m ˙ e + m ˙ w = m ˙ e + m ˙ e + m ˙ w + m ˙ w = 0 \dot m_e+\dot m_w=\dot m_e^*+\dot m_e'+\dot m_w^*+\dot m_w'=0
能够写为
m ˙ e + m ˙ w = m ˙ e m ˙ w \dot m_e'+\dot m_w'=-\dot m_e^*-\dot m_w^*
这个形式的连续方程很是有趣,其代表,一旦算得的质量流量是知足连续方程的精确解的话,那么RHS(右端项)是零,从而修正场也是零(即,若是知足连续方程话就无需再修正了)。所以,正是当前场的质量守恒偏差驱动着修正场,单元面处的质量流量和质量流量修正为
m ˙ e = ρ v e S e = ρ u e S e x = ρ u e Δ y e m ˙ w = ρ v w S w = ρ u w S w x = ρ u w Δ y w \begin{aligned} \dot m_e^* &= \rho \bold v_e^* \cdot \bold S_e=\rho u_e^* S_e^x = \rho u_e^* \Delta y_e \\ \dot m_w^* &= \rho \bold v_w^* \cdot \bold S_w=\rho u_w^* S_w^x = -\rho u_w^* \Delta y_w \end{aligned}

m ˙ e = ρ v e S e = ρ u e S e x = ρ u e Δ y e m ˙ w = ρ v w S w = ρ u w S w x = ρ u w Δ y w \begin{aligned} \dot m_e' &= \rho \bold v_e' \cdot \bold S_e=\rho u_e' S_e^x = \rho u_e' \Delta y_e \\ \dot m_w' &= \rho \bold v_w' \cdot \bold S_w=\rho u_w' S_w^x = -\rho u_w' \Delta y_w \end{aligned}
上式中用到了 S e x = Δ y e S_e^x=\Delta y_e S w x = Δ y w S_w^x=-\Delta y_w

压力场并无出如今修正形式的连续方程 m ˙ e + m ˙ w = m ˙ e m ˙ w \dot m_e'+\dot m_w'=-\dot m_e^*-\dot m_w^* 中,为了将其引入到式中,须要使用动量方程的离散形式。该过程先把离散动量方程 a e u u e + f N B ( e ) a f u u f = b e u V e ( p x ) e a_e^u u_e+\displaystyle\sum_{f\sim NB(e)}a_f^u u_f=b_e^u-V_e\left(\dfrac{\partial p}{\partial x}\right)_e 写成更加紧凑的形式
u e + H e ( u e ) = B e u D e u ( p x ) e u_e+H_e(u_e) = B_e^u - D_e^u\left(\frac{\partial p}{\partial x}\right)_e
其中
H e ( u e ) = f N B ( e ) a f u a e u u f B e u = b e u a e u D e u = V e a e u H_e(u_e) =\sum_{f\sim NB(e)}\frac{a_f^u}{a_e^u} u_f \quad\quad B_e^u = \frac{b_e^u}{a_e^u} \quad\quad D_e^u = \frac{V_e}{a_e^u}
对于计算的速度场 u e u_e^* ,上式变为
u e + H e ( u e ) = B e u D e u ( p ( n ) x ) e u_e^*+H_e(u_e^*) = B_e^u - D_e^u\left(\frac{\partial p^{(n)}}{\partial x}\right)_e
用精确速度场表示的动量方程 u e + H e ( u e ) = B e u D e u ( p x ) e u_e+H_e(u_e) = B_e^u - D_e^u\left(\frac{\partial p}{\partial x}\right)_e ,减去,算得速度场所表示的动量方程 u e + H e ( u e ) = B e u D e u ( p ( n ) x ) e u_e^*+H_e(u_e^*) = B_e^u - D_e^u\left(\frac{\partial p^{(n)}}{\partial x}\right)_e ,可获得修正场的方程为
u e + H e ( u ) = D e u ( p x ) e u_e'+\underline{H_e(u')} = - D_e^u\left(\frac{\partial p'}{\partial x}\right)_e
一样也能够推导出 w w 面的修正方程
u w + H w ( u ) = D w u ( p x ) w u_w'+\underline{H_w(u')} = - D_w^u\left(\frac{\partial p'}{\partial x}\right)_w
m ˙ e = ρ u e Δ y e \dot m_e' = \rho u_e' \Delta y_e m ˙ w = ρ u w Δ y w \dot m_w' = -\rho u_w' \Delta y_w 代入到连续方程 m ˙ e + m ˙ w = m ˙ e m ˙ w \dot m_e'+\dot m_w'=-\dot m_e^*-\dot m_w^* ,得
ρ e u e Δ y e + ( ρ w u w Δ y w ) = ( m ˙ e + m ˙ w ) \rho_e u_e' \Delta y_e+(-\rho_w u_w' \Delta y_w)=-(\dot m_e^*+\dot m_w^*)
而后,把上上式和上上上式的 u e u_e' u w u_w' 的离散形式代入上式,可获得压力修正项的方程
ρ e [ H e ( u ) D e u ( p x ) e ] Δ y e ρ w [ H w ( u ) D w u ( p x ) w ] Δ y w = ( m ˙ e + m ˙ w ) \begin{aligned} &\rho_e \left[ -\underline{H_e(u')} - D_e^u\left(\frac{\partial p'}{\partial x}\right)_e \right] \Delta y_e \\ &-\rho_w \left[ -\underline{H_w(u')} - D_w^u\left(\frac{\partial p'}{\partial x}\right)_w \right] \Delta y_w=-(\dot m_e^*+\dot m_w^*) \end{aligned}
在该方程中,压力场以扩散项的形式出现,离散后变为
ρ e [ H e ( u ) D e u ( p E p C Δ x ) e ] Δ y e + ρ w [ H w ( u ) D w u ( p C p W Δ x ) w ] ( Δ y w ) = ( m ˙ e + m ˙ w ) \begin{aligned} &\rho_e \left[ -\underline{H_e(u')} - D_e^u\left(\frac{p_E'-p_C'}{\Delta x}\right)_e \right] \Delta y_e \\ &+\rho_w \left[ -\underline{H_w(u')} - D_w^u\left(\frac{p_C'-p_W'}{\Delta x}\right)_w \right] (-\Delta y_w)=-(\dot m_e^*+\dot m_w^*) \end{aligned}

ρ e D e u ( Δ y e Δ x e ) ( p E p C ) ρ w D w u ( Δ y w Δ x w ) ( p C p W ) = ( m ˙ e + m ˙ w ) + ρ e H e ( u ) Δ y e + ρ w H w ( u ) ( Δ y w ) \begin{aligned} & -\rho_e D_e^u\left(\frac{\Delta y_e}{\Delta x_e}\right)(p_E'-p_C') -\rho_w D_w^u\left(-\frac{\Delta y_w}{\Delta x_w}\right) (p_C'-p_W')\\ &=-(\dot m_e^*+\dot m_w^*)+\rho_e \underline{H_e(u')}\Delta y_e+\rho_w \underline{H_w(u')}(-\Delta y_w) \end{aligned}
整理后的压力修正方程形式为
a C p p C + a E p p E + a W p p W = b C p a_C^{p'}p_C'+a_E^{p'}p_E'+a_W^{p'}p_W'=b_C^{p'}
其中
a E p = ρ e D e u Δ y e Δ x e a W p = ρ w D w u Δ y w Δ x w a C p = ( a E p + a W p ) b C p = ( m ˙ e + m ˙ w ) + [ ρ e Δ y e H e ( u ) ρ w Δ y w H w ( u ) ] \begin{aligned} a_E^{p'} &= -\frac{\rho_e D_e^u\Delta y_e}{\Delta x_e} \\ a_W^{p'} &= -\frac{\rho_w D_w^u\Delta y_w}{\Delta x_w}\\ a_C^{p'} &= -(a_E^{p'}+a_W^{p'}) \\ b_C^{p'} &= -(\dot m_e^*+\dot m_w^*)+ \underline{[\rho_e \Delta y_e H_e(u')-\rho_w \Delta y_w H_w(u')]} \end{aligned}
注意,对于上式下划线所标记的项,在达到稳态收敛的状况下,其值将变为零。所以,它们对于最终的解将没有影响。对于这些项的不一样近似方法将会衍生不一样的算法,这将在下面详细展开。在最初的SIMPLE算法中,这些项是简单地直接忽略掉了。此外,对于一维问题,且面积 Δ y \Delta y 为定值的状况,可直接把它设成 1 1 并从方程中刨掉便可。

2.6 交错网格上的SIMPLE算法

使用动量方程和压力修正方程,能够得到流动问题的解。在SIMPLE算法中的求解过程是,在每步迭代中,交替地生成压力和速度场,并前后知足动量方程和连续方程,而后再逐步逼近最终解(连续方程和动量方程这两个方程都知足的最终解)。该流程并不是是同步进行的,这种求解方程的法子常被称为分离式方法(segregated approach)(实际上还有一种把压力和速度直接耦合起来同步求解的法子呢,不过比较耗内存了)。该分离式的SIMPLE算法的计算流程汇总以下:

  1. 从某个猜想的压力场 p ( n ) p^{(n)} 和速度场 u ( n ) u^{(n)} 开始;
  2. 求解动量方程来得到新的速度场 u f u_f^*
    a e u u e + f N B ( e ) a f u u f = b e u V e ( p ( n ) x ) e a w u u w + f N B ( e ) a f u u f = b w u V w ( p ( n ) x ) w a_e^u u_e^*+\sum_{f\sim NB(e)}a_f^u u_f^*= b_e^u-V_e\left(\frac{\partial p^{(n)}}{\partial x}\right)_e \\ a_w^u u_w^*+\sum_{f\sim NB(e)}a_f^u u_f^*= b_w^u-V_w\left(\frac{\partial p^{(n)}}{\partial x}\right)_w
  3. 使用知足动量方程的速度场 u f u_f^* 来更新质量流量 m ˙ f \dot m_f^*
    m ˙ f = ρ v f S f = ρ f u f S f x \dot m_f^* = \rho \bold v_f^* \cdot \bold S_f=\rho_f u_f^* S_f^x
  4. 使用新的质量流量 m ˙ f \dot m_f^* 求解压力修正方程来获取压力修正场 p p'
    a C p p C + a E p p E + a W p p W = b C p a_C^{p'}p_C'+a_E^{p'}p_E'+a_W^{p'}p_W'=b_C^{p'}
  5. 更新压力场和速度场以获取知足连续方程的场
    u f = u f + u f u f = D f u ( p x ) f p C = p C ( n ) + p C m ˙ f = m ˙ f + m ˙ f m ˙ f = ρ f D f u Δ y f ( p x ) f \begin{aligned} u_f^{**} &= u_f^* + u_f' \quad\quad u_f'=-D_f^u\left(\frac{\partial p'}{\partial x}\right)_f \\ p_C^* &=p_C^{(n)} + p_C' \\ \dot m_f^{**} &= \dot m_f^* + \dot m_f' \quad\quad \dot m_f'=-\rho_f D_f^u \Delta y_f \left(\frac{\partial p'}{\partial x}\right)_f \end{aligned}
  6. u ( n ) = u u^{(n)}=u^{**} p ( n ) = p p^{(n)}=p^{*}
  7. 返回第2步,并重复该流程直到收敛。

理论是晦涩的,而实例则是生动的,用个简单的小例子来展现下SIMPLE的强大是再好不过的了。


例1 管道网络(管网)中的流动

下图展现的是水力管网系统中的一部分,管道流动中的动量方程能够写成
m ˙ = ρ u A = D P \dot m=\rho u A = -D\nabla P
其中, D A = 0.5 D_A=0.5 D B = D F = 0.4 D_B=D_F=0.4 D C = D E = 0.3 D_C=D_E=0.3 D D = 0.19 D_D=0.19 D G = 0.1875 D_G=0.1875 D H = 0.35 D_H=0.35 。使用SIMPLE算法,计算系统中的未知质量流量和压力。

在这里插入图片描述

在该系统中,使用质量流量场做为变量,来替代速度场。这并不会产生什么问题,由于该例中的动量方程已经忽略掉了对流项和扩散项,仅保留了压力梯度项,即认为对流和扩散效应,与压差驱动力相比很是小,直接忽略了。

使用SIMPLE算法的求解要先用假设的压力场来求解一个初始速度场,而后再根据这个计算的速度场预测一个压力场,以知足连续方程。

该流程汇总以下:

  1. 由一个猜想的压力场开始;
  2. 使用给定的动量方程计算质量流量;
  3. 构造压力修正方程来知足连续性(质量守恒),并用其修正压力和速度场。

该例中无需迭代,由于其不含对流项所带来的非线性效应(若是是用SIMPLE来求解不可压缩流动,那是要迭代屡次才能收敛的)。

step 1
给定待求解位置处的压力猜想值,以下:
p 3 ( n ) = 300 p 6 ( n ) = 200 p 8 ( n ) = 120 p_3^{(n)}=300 \quad\quad p_6^{(n)}=200 \quad\quad p_8^{(n)}=120
固然也可使用其它值。

step 2
基于假设的压力值,使用动量方程计算不一样的质量流量:
m ˙ A = D A ( p 1 p 3 ( n ) ) = 0.5 ( 400 300 ) = 50 m ˙ B = D B ( p 3 ( n ) p 2 ) = 0.4 ( 300 350 ) = 20 m ˙ C = D C ( p 4 p 3 ( n ) ) = 0.3 ( 50 300 ) = 75 m ˙ D = D D ( p 3 ( n ) p 6 ( n ) ) = 0.19 ( 300 200 ) = 19 m ˙ E = D E ( p 6 ( n ) p 5 ) = 0.3 ( 200 300 ) = 30 m ˙ F = D F ( p 7 p 6 ( n ) ) = 0.4 ( 80 200 ) = 48 m ˙ G = D G ( p 6 ( n ) p 8 ( n ) ) = 0.1875 ( 200 120 ) = 15 m ˙ H = D H ( p 9 p 8 ( n ) ) = 0.35 ( 200 120 ) = 28 \begin{aligned} & \dot m_A^*=D_A(p_1-p_3^{(n)})=0.5*(400-300)=50 \\ & \dot m_B^*=D_B(p_3^{(n)}-p_2)=0.4*(300-350)=-20 \\ & \dot m_C^*=D_C(p_4-p_3^{(n)})=0.3*(50-300)=-75 \\ & \dot m_D^*=D_D(p_3^{(n)}-p_6^{(n)})=0.19*(300-200)=19 \\ & \dot m_E^*=D_E(p_6^{(n)}-p_5)=0.3*(200-300)=-30 \\ & \dot m_F^*=D_F(p_7-p_6^{(n)})=0.4*(80-200)=-48 \\ & \dot m_G^*=D_G(p_6^{(n)}-p_8^{(n)})=0.1875*(200-120)=15 \\ & \dot m_H^*=D_H(p_9-p_8^{(n)})=0.35*(200-120)=28 \end{aligned}

step 3
检查质量流量是否知足连续性,对全部内部节点计算 k m ˙ k \displaystyle\sum_{\sim k}\dot m_k^* ,得
Node 3: k ( m ˙ k ) = m ˙ B + m ˙ D m ˙ A m ˙ C = 20 + 19 50 + 75 = 24 Node 6: k ( m ˙ k ) = m ˙ E + m ˙ G m ˙ D m ˙ F = 30 + 15 19 + 48 = 14 Node 8: k ( m ˙ k ) = m ˙ I m ˙ G m ˙ H = 50 15 28 = 7 \begin{aligned} & \text{Node 3:\quad}\sum_{\sim k}(\dot m_k^*)=\dot m_B^*+\dot m_D^*-\dot m_A^*-\dot m_C^*=-20+19-50+75=24 \\ & \text{Node 6:\quad}\sum_{\sim k}(\dot m_k^*)=\dot m_E^*+\dot m_G^*-\dot m_D^*-\dot m_F^*=-30+15-19+48=14 \\ & \text{Node 8:\quad}\sum_{\sim k}(\dot m_k^*)=\dot m_I^*-\dot m_G^*-\dot m_H^*=50-15-28=7 \end{aligned}
因为质量守恒并不知足,须要修正处理,压力修正方程推导以下:
k ( m ˙ k + m ˙ k ) = 0 k ( m ˙ k ) = k ( m ˙ k ) \sum_{\sim k}(\dot m_k^*+\dot m_k')=0 \Rightarrow \sum_{\sim k}(\dot m_k')=-\sum_{\sim k}(\dot m_k^*)
基于动量方程,把质量流量的修正用压力修正表示为
m ˙ A = D A ( p 3 ) m ˙ B = D B ( p 3 ) m ˙ C = D C ( p 3 ) m ˙ D = D D ( p 3 p 6 ) m ˙ E = D E ( p 6 ) m ˙ F = D F ( p 6 ) m ˙ G = D G ( p 6 p 8 ) m ˙ H = D H ( p 8 ) \begin{aligned} & \dot m'_A = D_A(-p'_3) \\ & \dot m_B'=D_B(p_3') \\ & \dot m_C'=D_C(-p_3') \\ & \dot m_D'=D_D(p_3'-p_6') \\ & \dot m_E'=D_E(p_6') \\ & \dot m_F'=D_F(-p_6') \\ & \dot m_G'=D_G(p_6'-p_8') \\ & \dot m_H'=D_H(-p_8') \end{aligned}
注意 p 1 , 2 , 4 , 5 , 7 p'_{1,2,4,5,7} 为0,由于这些压力值是已知的,自己就是精确值,因此无需修正。

节点 3 6 8 三、六、8 处的质量守恒方程(连续方程)为
m ˙ B + m ˙ D m ˙ A m ˙ C = ( m ˙ B + m ˙ D m ˙ A m ˙ C ) m ˙ E + m ˙ G m ˙ D m ˙ F = ( m ˙ E + m ˙ G m ˙ D m ˙ F ) m ˙ G m ˙ H = ( m ˙ I m ˙ G m ˙ H ) \begin{aligned} \dot m_B'+\dot m_D'-\dot m_A'-\dot m_C' &= -(\dot m_B^*+\dot m_D^*-\dot m_A^*-\dot m_C^*) \\ \dot m_E'+\dot m_G'-\dot m_D'-\dot m_F' &= -(\dot m_E^*+\dot m_G^*-\dot m_D^*-\dot m_F^*)\\ -\dot m_G'-\dot m_H' &= -(\dot m_I^*-\dot m_G^*-\dot m_H^*) \end{aligned}
把压力修正表示的质量流量修正代入上式,获得压力修正方程
D B ( p 3 ) + D D ( p 3 p 6 ) D A ( p 3 ) D C ( p 3 ) = ( m ˙ B + m ˙ D m ˙ A m ˙ C ) D E ( p 6 ) + D G ( p 6 p 8 ) D