精华内容
下载资源
问答
  • 不可压缩couette流动隐式crank nicolson 技术,能够有效的解决,此程序实现《计算流体力学入门》Anderson著不可压缩Couette流动,使用Crank-Nicolson差分格式
  • 由于本章内容繁杂,篇幅较长,故分成了四部分来讲解,各部分主要内容分别为:交错网格、...这里是第一部分,主要讲解不可压缩流动问题的求解在常规网格上所碰到的问题,以及,交错网格上的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 Flows

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

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

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

    1 主要困难

    前面几章所处理的通用守恒方程,可以改写成与连续方程或动量方程相似的形式。然而,目前为止所呈现的数值技术,并不足以求解Navier-Stokes方程,因为求解一般流动的算法要能够处理压力速度的耦合问题。为理解该问题,先把连续方程和动量方程写出来
    ρt+(ρv)=0t(ρv)+(ρvv)=p+{μ[v+(v)T]}+fb\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})则基于前次迭代值做显式估算并添加到源项中去。针对上述连续方程和动量方程,从通有标量方程的数值方法(前面几章讲解的方法)中无法汲取灵感来解决的主要难题在于,出现在动量方程中的压力场无法获得其显式计算方程。

    从这俩方程中可以看出,纵然速度场可以直接使用动量方程来计算,然而出现在动量方程中的压力场无法直接从连续方程中计算。该耦合关系(压力和速度的耦合)是如此强烈却又如此隐含,为便于理解,可把方程组写成如下矩阵形式
    Au=[FBTB0][vp]=[fb0] \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)问题的特性,表明其无法通过任何迭代手法来求解压力和速度场。因此,需要推导出关于压力的方程才好。

    一种处理方法是,通过把矩阵A\bold A分解成下三角L\bold L和上三角U\bold U矩阵,将动量和连续方程系统重新改写成下述形式
    A=[FBTB0]=[F0BBF1BT][IF1BT0I]=LU \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
    其中BF1BT-\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补矩阵)来创建的。

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

    该算法可用如下的矩阵形式来描述
    [ID1BT0I][vp]=[vp] \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}
    接下来使用下式来更新速度场
    [F0BBD1BT][vp]=[fb0] \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}
    上式和上上式中的F1\bold F^{-1}用其逆对角阵D1\bold D^{-1}来近似,上标^*表示当前迭代的中间值。流程汇总如下:

    • 求解:Fv=fb\bold F \bold v^*=\bold f_b
    • 求解:BD1BTp=Bv-\bold B \bold D^{-1} \bold B^\text T \bold p^*=-\bold B \bold v^*
    • 更新:v=vD1BTp\bold v=\bold v^*-\bold D^{-1}\bold B^\text T\bold p^*
    • 更新:p=p\bold p=\bold p^*

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

    2 初步推导

    在这里插入图片描述

    针对不可压缩流动问题发展求解算法的困难之处,可通过对上图所示的一维空间均匀网格问题的离散处理来加深理解。为简便起见,假设流动是定常的。简化的连续和动量方程(写成守恒形式)为
    (ρu)x=0 (ρuu)x=x(μux)px\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 动量方程离散

    动量方程的离散,首先要对其在单元CC上做积分,得
    VC(ρuu)xdV=VCx(μux)dVVCpxdV \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
    使用散度定理,将对流项和扩散项的体积分转化成面积分
    VC(ρuu)dy=VCμuxdyVCpxdV \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)em˙eue+(ρuΔy)wm˙wuw=(μuxΔy)e(μuxΔy)wVCpxdV \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˙eue+m˙wuwConvection[(μuxΔy)e(μuxΔy)w]Diffusion=VCpxdV \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
    对流和扩散项可用之前章节中所讲的方法来离散化处理,得到如下形式的代数方程
    aCuuC+FNB(C)(aFuuF)=bCuVCpxdV 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 连续方程的离散

    连续方程的离散,也是要对其在单元CC上做积分,得
    VC(ρu)xdV=0 \int_{V_C}\frac{\partial(\rho u)}{\partial x}dV=0
    再次使用散度定理将体积分转化成面积分,然后转化成单元面通量的加和形式,连续方程的离散形式为
    fnb(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

    fnb(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积分来计算体积分
    VCpxdV=(px)CVC \int_{V_C}\frac{\partial p}{\partial x}dV=\left(\frac{\partial p}{\partial x}\right)_CV_C
    使用中心差分格式,上式的离散形式为
    VCpxdV=pEpW2ΔxVC \int_{V_C}\frac{\partial p}{\partial x}dV=\frac{p_E-p_W}{2\Delta x} V_C
    第二种方法,把压力梯度项的体积分转化成面积分
    VCpxdV=VCpdy \int_{V_C}\frac{\partial p}{\partial x}dV=\int_{\partial V_C}p dy
    把面积分写成单元面通量加和的形式,有
    VCpxdV=VCpdy=peΔyepwΔyw=(pepw)Δy=(pepw)VCΔ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}
    压力的变化选用线性插值廓线,把压力梯度项重新写成主节点上压力值的函数形式
    VCpxdV=[12(pE+pC)12(pC+pW)]VCΔx=pEpW2ΔxVC \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
    两个方法导出了同样的表达式,即,包含交替节点EEWW之间压力差值的表达式(这里的交替指的是,EEWW节点之间间隔了一个节点CC)。

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

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

    在这里插入图片描述

    针对上图所示的压力和速度值,在点WWCCEE处的压力梯度为
    VWpxdV=(pCpWW)VW2ΔxW=(1010)VW2ΔxW=0VCpxdV=(pEpW)VC2ΔxC=(100+100)VC2ΔxC=0VEpxdV=(pEEpC)VE2ΔxE=(1010)VE2ΔxE=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}
    而且连续方程看起来在每个单元上都是满足的,因为
    VWuxdV=(uCuWW)VW2ΔxW=(11)VW2ΔxW=0VCuxdV=(uEuW)VC2ΔxC=(1010)VC2ΔxC=0VEuxdV=(uEEuC)VE2ΔxE=(11)VE2ΔxE=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)。
    在这里插入图片描述

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

    ueuw=0 u_e-u_w=0
    无需再做插值了,因为速度在界面eeww处是可以获得的。此外,动量方程的积分是在单元ee上进行的,其产生如下的离散动量方程形式
    aeuue+fNB(e)afuuf=beuVe(p)e=beuVepEpCδxe 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)算法的最初实用形式。

    由上图所示网格上的离散连续方程和动量方程出发
    fnb(C)m˙f=0aeuue+fNB(e)afuuf=beuVe(px)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)}表示。在每步迭代中,先求解动量方程来获取速度场,所得解用上标^*表示,因为其并非当前迭代步的最终解(注意,由于并未考虑连续方程,所以这时候的速度场并不满足连续方程)。这样,动量方程满足
    aeuue+fNB(e)afuuf=beuVe(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
    其中压力场仍旧采用的是前次迭代值,那么,算出来的速度场uu^*满足动量方程,但是未必满足连续方程,因为压力场并不是精确值。因此,需要寻找修正措施来确保速度(或质量流量)和压力满足连续方程。

    将修正场用上标撇来表示,即,(uu'pp'),那么速度和压力修正为
    u=u+up=p+p u=u^*+u' \\ p=p^*+p'
    注意,单元面处的质量流量也修正为
    m˙f=m˙f+ρuSfx=m˙f+mf \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˙em˙w \dot m_e'+\dot m_w'=-\dot m_e^*-\dot m_w^*
    这个形式的连续方程非常有趣,其表明,一旦算得的质量流量是满足连续方程的精确解的话,那么RHS(右端项)是零,从而修正场也是零(即,如果满足连续方程话就无需再修正了)。因此,正是当前场的质量守恒误差驱动着修正场,单元面处的质量流量和质量流量修正为
    m˙e=ρveSe=ρueSex=ρueΔyem˙w=ρvwSw=ρuwSwx=ρuwΔyw\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=ρveSe=ρueSex=ρueΔyem˙w=ρvwSw=ρuwSwx=ρuwΔyw\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}
    上式中用到了Sex=ΔyeS_e^x=\Delta y_eSwx=ΔywS_w^x=-\Delta y_w

    压力场并没有出现在修正形式的连续方程m˙e+m˙w=m˙em˙w\dot m_e'+\dot m_w'=-\dot m_e^*-\dot m_w^*中,为了将其引入到式中,需要使用动量方程的离散形式。该过程先把离散动量方程aeuue+fNB(e)afuuf=beuVe(px)ea_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写成更加紧凑的形式
    ue+He(ue)=BeuDeu(px)e u_e+H_e(u_e) = B_e^u - D_e^u\left(\frac{\partial p}{\partial x}\right)_e
    其中
    He(ue)=fNB(e)afuaeuufBeu=beuaeuDeu=Veaeu 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}
    对于计算的速度场ueu_e^*,上式变为
    ue+He(ue)=BeuDeu(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
    用精确速度场表示的动量方程ue+He(ue)=BeuDeu(px)eu_e+H_e(u_e) = B_e^u - D_e^u\left(\frac{\partial p}{\partial x}\right)_e,减去,算得速度场所表示的动量方程ue+He(ue)=BeuDeu(p(n)x)eu_e^*+H_e(u_e^*) = B_e^u - D_e^u\left(\frac{\partial p^{(n)}}{\partial x}\right)_e,可得到修正场的方程为
    ue+He(u)=Deu(px)e u_e'+\underline{H_e(u')} = - D_e^u\left(\frac{\partial p'}{\partial x}\right)_e
    同样也可以推导出ww面的修正方程
    uw+Hw(u)=Dwu(px)w u_w'+\underline{H_w(u')} = - D_w^u\left(\frac{\partial p'}{\partial x}\right)_w
    m˙e=ρueΔye\dot m_e' = \rho u_e' \Delta y_em˙w=ρuwΔyw\dot m_w' = -\rho u_w' \Delta y_w代入到连续方程m˙e+m˙w=m˙em˙w\dot m_e'+\dot m_w'=-\dot m_e^*-\dot m_w^*,得
    ρeueΔye+(ρwuwΔyw)=(m˙e+m˙w) \rho_e u_e' \Delta y_e+(-\rho_w u_w' \Delta y_w)=-(\dot m_e^*+\dot m_w^*)
    然后,把上上式和上上上式的ueu_e'uwu_w'的离散形式代入上式,可得到压力修正项的方程
    ρe[He(u)Deu(px)e]Δyeρw[Hw(u)Dwu(px)w]Δyw=(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[He(u)Deu(pEpCΔx)e]Δye+ρw[Hw(u)Dwu(pCpWΔx)w](Δyw)=(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}

    ρeDeu(ΔyeΔxe)(pEpC)ρwDwu(ΔywΔxw)(pCpW)=(m˙e+m˙w)+ρeHe(u)Δye+ρwHw(u)(Δyw)\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}
    整理后的压力修正方程形式为
    aCppC+aEppE+aWppW=bCp a_C^{p'}p_C'+a_E^{p'}p_E'+a_W^{p'}p_W'=b_C^{p'}
    其中
    aEp=ρeDeuΔyeΔxeaWp=ρwDwuΔywΔxwaCp=(aEp+aWp)bCp=(m˙e+m˙w)+[ρeΔyeHe(u)ρwΔywHw(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为定值的情况,可直接把它设成11并从方程中刨掉即可。

    2.6 交错网格上的SIMPLE算法

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

    1. 从某个猜测的压力场p(n)p^{(n)}和速度场u(n)u^{(n)}开始;
    2. 求解动量方程来获得新的速度场ufu_f^*
      aeuue+fNB(e)afuuf=beuVe(p(n)x)eawuuw+fNB(e)afuuf=bwuVw(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. 使用满足动量方程的速度场ufu_f^*来更新质量流量m˙f\dot m_f^*
      m˙f=ρvfSf=ρfufSfx \dot m_f^* = \rho \bold v_f^* \cdot \bold S_f=\rho_f u_f^* S_f^x
    4. 使用新的质量流量m˙f\dot m_f^*求解压力修正方程来获取压力修正场pp'
      aCppC+aEppE+aWppW=bCp a_C^{p'}p_C'+a_E^{p'}p_E'+a_W^{p'}p_W'=b_C^{p'}
    5. 更新压力场和速度场以获取满足连续方程的场
      uf=uf+ufuf=Dfu(px)fpC=pC(n)+pCm˙f=m˙f+m˙fm˙f=ρfDfuΔyf(px)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)=uu^{(n)}=u^{**}p(n)=pp^{(n)}=p^{*}
    7. 返回第2步,并重复该流程直到收敛。

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


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

    下图展示的是水力管网系统中的一部分,管道流动中的动量方程可以写成
    m˙=ρuA=DP \dot m=\rho u A = -D\nabla P
    其中,DA=0.5D_A=0.5DB=DF=0.4D_B=D_F=0.4DC=DE=0.3D_C=D_E=0.3DD=0.19D_D=0.19DG=0.1875D_G=0.1875DH=0.35D_H=0.35。使用SIMPLE算法,计算系统中的未知质量流量和压力。

    在这里插入图片描述

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

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

    该流程汇总如下:

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

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

    step 1
    给定待求解位置处的压力猜测值,如下:
    p3(n)=300p6(n)=200p8(n)=120 p_3^{(n)}=300 \quad\quad p_6^{(n)}=200 \quad\quad p_8^{(n)}=120
    当然也可以使用其它值。

    step 2
    基于假设的压力值,使用动量方程计算不同的质量流量:
    m˙A=DA(p1p3(n))=0.5(400300)=50m˙B=DB(p3(n)p2)=0.4(300350)=20m˙C=DC(p4p3(n))=0.3(50300)=75m˙D=DD(p3(n)p6(n))=0.19(300200)=19m˙E=DE(p6(n)p5)=0.3(200300)=30m˙F=DF(p7p6(n))=0.4(80200)=48m˙G=DG(p6(n)p8(n))=0.1875(200120)=15m˙H=DH(p9p8(n))=0.35(200120)=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
    检查质量流量是否满足连续性,对所有内部节点计算km˙k\displaystyle\sum_{\sim k}\dot m_k^*,得
    Node 3:k(m˙k)=m˙B+m˙Dm˙Am˙C=20+1950+75=24Node 6:k(m˙k)=m˙E+m˙Gm˙Dm˙F=30+1519+48=14Node 8:k(m˙k)=m˙Im˙Gm˙H=501528=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)=0k(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=DA(p3)m˙B=DB(p3)m˙C=DC(p3)m˙D=DD(p3p6)m˙E=DE(p6)m˙F=DF(p6)m˙G=DG(p6p8)m˙H=DH(p8)\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}
    注意p1,2,4,5,7p'_{1,2,4,5,7}为0,因为这些压力值是已知的,本身就是精确值,所以无需修正。

    节点3683、6、8处的质量守恒方程(连续方程)为
    m˙B+m˙Dm˙Am˙C=(m˙B+m˙Dm˙Am˙C)m˙E+m˙Gm˙Dm˙F=(m˙E+m˙Gm˙Dm˙F)m˙Gm˙H=(m˙Im˙Gm˙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}
    把压力修正表示的质量流量修正代入上式,得到压力修正方程
    DB(p3)+DD(p3p6)DA(p3)DC(p3)=(m˙B+m˙Dm˙Am˙C)DE(p6)+DG(p6p8)DD(p3p6)DF(p6)=(m˙E+m˙Gm˙Dm˙F)DG(p6p8)DH(p8)=(m˙Im˙Gm˙H)\begin{aligned} D_B(p_3')+D_D(p_3'-p_6')-D_A(-p'_3)-D_C(-p_3') &= -(\dot m_B^*+\dot m_D^*-\dot m_A^*-\dot m_C^*) \\ D_E(p_6')+D_G(p_6'-p_8')-D_D(p_3'-p_6')-D_F(-p_6') &= -(\dot m_E^*+\dot m_G^*-\dot m_D^*-\dot m_F^*)\\ -D_G(p_6'-p_8')-D_H(-p_8') &= -(\dot m_I^*-\dot m_G^*-\dot m_H^*) \end{aligned}
    整理后,压力修正方程变为
    1.39(p3)0.19(p6)=(m˙B+m˙Dm˙Am˙C)1.0775(p6)0.1875(p8)0.19(p3)=(m˙E+m˙Gm˙Dm˙F)0.1875(p6)+0.5375(p8)=(m˙Im˙Gm˙H)\begin{aligned} 1.39(p_3')-0.19(p_6') &= -(\dot m_B^*+\dot m_D^*-\dot m_A^*-\dot m_C^*) \\ 1.0775(p_6')-0.1875(p_8')-0.19(p_3') &= -(\dot m_E^*+\dot m_G^*-\dot m_D^*-\dot m_F^*)\\ -0.1875(p_6')+0.5375(p_8') &= -(\dot m_I^*-\dot m_G^*-\dot m_H^*) \end{aligned}
    把算得的初始质量流量代入进去,得到修正方程为
    1.39(p3)0.19(p6)=241.0775(p6)0.1875(p8)0.19(p3)=140.1875(p6)+0.5375(p8)=7\begin{aligned} 1.39(p_3')-0.19(p_6') &= -24 \\ 1.0775(p_6')-0.1875(p_8')-0.19(p_3') &= -14 \\ -0.1875(p_6')+0.5375(p_8') &= -7 \end{aligned}
    求解压力修正方程,得
    {p3=20p6=20p8=20\begin{cases} p'_3=-20 \\ p'_6=-20 \\ p'_8=-20 \end{cases}
    使用算得的修正压力,更新速度场和压力场,以产生满足质量守恒的速度场。修正后的质量流量为:
    m˙A=m˙A+m˙A=m˙A+DA(p3)=50+0.5((20))=60m˙B=m˙B+m˙B=m˙B+DB(p3)=20+0.4(20)=28m˙C=m˙C+m˙C=m˙C+DC(p3)=75+0.3((20))=69m˙D=m˙D+m˙D=m˙D+DD(p3p6)=19+0.19(20(20))=19m˙E=m˙E+m˙E=m˙E+DE(p6)=30+0.3(20)=36m˙F=m˙F+m˙F=m˙F+DF(p6)=48+0.4((20))=40m˙G=m˙G+m˙G=m˙G+DG(p6p8)=15+0.1875(20(20))=15m˙H=m˙H+m˙H=m˙H+DH(p8)=28+0.35((20))=35\begin{aligned} & \dot m_A^{**}=\dot m_A^*+\dot m'_A = \dot m_A^*+D_A(-p'_3)=50+0.5(-(-20))=60 \\ & \dot m_B^{**}=\dot m_B^*+\dot m_B'=\dot m_B^*+D_B(p_3')=-20+0.4(-20)=-28 \\ & \dot m_C^{**}=\dot m_C^*+\dot m_C'=\dot m_C^*+D_C(-p_3')=-75+0.3(-(-20))=-69 \\ & \dot m_D^{**}=\dot m_D^*+\dot m_D'=\dot m_D^*+D_D(p_3'-p_6')=19+0.19(-20-(-20))=19 \\ & \dot m_E^{**}=\dot m_E^*+\dot m_E'=\dot m_E^*+D_E(p_6')=-30+0.3(-20)=-36 \\ & \dot m_F^{**}=\dot m_F^*+\dot m_F'=\dot m_F^*+D_F(-p_6')=-48+0.4(-(-20))=-40 \\ & \dot m_G^{**}=\dot m_G^*+\dot m_G'=\dot m_G^*+D_G(p_6'-p_8')=15+0.1875(-20-(-20))=15 \\ & \dot m_H^{**}=\dot m_H^*+\dot m_H'=\dot m_H^*+D_H(-p_8')=28+0.35(-(-20))=35 \end{aligned}
    压力修正为
    p3=p3(n)+p3=300+(20)=280p6=p6(n)+p6