2.2SIMPLE 系列算法2.3PISO 算法（OpenFOAM 理论笔记系列）
2.2 SIMPLE 系列算法
2.2.1标准SIMPLE 算法
SIMPLE算法(Semi-Implicit Method for PressureLinked Equations)最初被设计⽤来求解稳态问题，即控制⽅程中不包含瞬态项的计算。按照1.3.3节的约定，我们假设计算开始的时候有初始的压⼒和速度值，待求的真值为。实际上，速度的初始值可以完全随意给定，因为我们完全⾸先可以将初始压⼒带⼊式2.1解出⼀个速度:
上述步骤称为动量预测，解得的速度称为预测速度。
预测速度并不是待求的真实速度，根据式上⼀节的推导，待求速度和待求压⼒之间满⾜压⼒泊松⽅程：
显然，是的函数，知道我们就能算出待求的压⼒，但⽬前我们只有求解式(2.13)得到的,SIMPLE算法在这⾥做了⼀个妥协与假设，其⽤代替计算，即：
其本质上是将⽅程(2.4)转化为
再带⼊连续性⽅程后进⾏离散得到的。式(2.15)当然是不精确的，在下⼀节中我们将分析式(2.15)中做出的简化的背后含义。同时，式(2.15)还可以⽤来根据式2.14解出的压⼒来更新速度以得到。这样，我们先通过求解式(2.13)得到预测速度，再根据式2.13得到待求压⼒，最后根据式2.15更新现有的速度得到。由于式(2.14)是不准确地，因此我们计算得到的和也是不准确的，我们将和作为新的和重新进⾏上述过程，直到计算收敛为⽌。SIMPLE算法的流程图如下：
图2.1 稳态SIMPLE算法框图
2.2.2 SIMPLE 算法的亚松弛
P ,o U o P ,n U n a +P U p r
a =N
∑
N U N
r −S ∇P
o
(2.13)
U r ∇⋅∇P =(a P 1
n )∇⋅(HbyA n )
ml和g(2.6)
HbyA n U n U n P n U r U r U n P n ∇⋅∇P =(a P 1
n )∇⋅(HbyA r )
(2.14)
=U p
n −HbyA r ∇P a P 1
n
(2.15)
U n U r P n U n P n U n P n U n P o U o
在上⼀节中，我们指出SIMPLE算法唯⼀的简化假设是将式(2.4)简化为了式(2.15)，在这⼀节我们详细分析⼀下这⼀简化背后的逻辑。为了⽅便表述，我们将式(2.13)中的压⼒改写为，表⽰为预测压⼒。
对于待求的压⼒和速度，同样存在着上述关系，即式2.1：
使⽤式(2.1)减式(2.16)可得：
其中，，，分别表⽰修正速度和修正压⼒，它们是真实值与当前计算值的差。SIMPLE算法认为，在计算过程中可以略去相邻点修正带来的影响，即将式(2.17)简化为：
将式2.17与式2.16相加即可得到SIMPLE算法中使⽤的式(2.15)。
这⼀略去相邻点修正的假设通常称为略去邻点假设，这⼀假设是⾮常激进(粗暴)的，在有些情况会导致速度和压⼒的修正过于剧烈，使模拟发散。解决这⼀问题最好的办法是使⽤亚松弛技术(under-relax)。
假设⼀次SIMPLE循环开始时候的初始值为和，这次循环结束后得到的计算值为，，我们并不使⽤计算值作为最终得到的
值，⽽是使⽤下⾯的值作为最终计算得到的值:
式(2.18)中为亚松弛因⼦，其值在0和1之间进⾏选择。显然，当时，我们将完全采⽤计算得到的值作为最终的值，但正如前⾯所提到的，这样做往往会导致修正过⼤，计算发散。如果，我们则将完全将旧值赋值给新值，这样相当于没做计算，也是不可取的。整体来说，愈趋近于1，表⽰速度和压⼒修正得越激进，越容易导致发散但如果可以收敛其会使⽤较少的迭代步，也就是计算效率较⾼。当趋近于0时，计算不容易发散，但是会需要更多的迭代步来达到收敛，计算效率会降低。合适的选择应该兼顾结果的收敛性和计算的效率。
2.2.3 SIMPLEC 算法
SIMPLEC(SIMPLE-Consistent)算法是⼀类修正的SIMPLE算法，其也认为SIMPLE算法中的略去邻点假设过于激进，为了避免这种激进的假设，其没有采⽤亚松弛技术，⽽是做了以下修正：
P o P r a +P U p r
a =N
∑
N U N
r −S ∇P
r
(2.16)
a +P U p n
a =N
∑
N U N
n −S ∇P
n
(2.1)
a +P U p ∗
a =N
∑
N U N
∗−∇P
∗
(2.17)
=U ∗−U n U r P =∗P −n P r a =P U
p ∗−∇P ∗
(2.17)
车身广告效果图
U o P o U tem
P tem =+α−U n U o U (U tem U o )P =P +αP −
P n
o
P
(tem o )
(2.18)
αα=1α=0ααα
⾸先，SIMPLEC算法认为当前点的修正速度应该是相邻点修正速度的平均，即：
将式2.19带⼊修正量的⽅程2.17可得：
整理可得：
在有限体积法中，为负值，因此，在相同的U_p^*时，使⽤2.21得到的修正压⼒要⽐式(2.17)⼩。将式(2.21)与式(2.16)相加可得：
整理可得：
将式2.23带⼊连续性⽅程即可得到压⼒泊松⽅程：
SIMPLEC算法与SIMPLE算法的流程完全⼀致，仅是将式(2.14)替换为式(2.24)，将式(2.15)替换为式(2.23)。
SIMPLE算法还有另⼀种修正版本，称为SIMPLER算法，另外，SIMPLE算法还可以拓展到瞬态计算中，称为瞬态SIMPLE算法，由于OpenFOAM中不涉及这两类算法，我们在这⾥不做介绍，有兴趣的读者可以参考相关⽂献。关于SIMPLE和SIMPLEC算法的内容除了参考原始⽂献也可以参考⽂献.
稻米油好吗
2.2.4 OpenFOAM 中的SIMPLE 算法
2.2.4.1 simpleFoam 代码解析
=U p ∗
a ∑N
a ∑N U N ∗
(2.19)
a +
P U p ∗a =(
N
∑
N )U p ∗
−∇P ∗
(2.20)
=U p ∗脸下垂怎么办
−∇P a +a P ∑N N 1
∗
(2.21)
学习麻辣烫
a ∑N N a +P U p n
a =N
∑
N U N r
−S ∇P −r ∇P a +a P ∑N N a P
∗
(2.22)
=
−
−
∇P −
∇P U p
n a P
S a P
a ∑N N U N
r
a P 1
r a +a P ∑N N 1∗=−∇P −
∇P HbyA r
a P 1
r a +a P ∑N N 1∗=−
∇P −∇P +∇P −∇P HbyA r
a P 1
r a +a P
∑N N 1∗a +a P ∑N N 1r a +a P ∑N N 1r =+−
∇P −∇P HbyA r (a +a P ∑N N 1
a P 1
)r a +a P ∑N N 1n
(2.23)
∇⋅∇P =(a +a P ∑N N 1n )∇⋅+−∇P [HbyA r
(a +a P ∑N N 1a P 1
)r ]
(2.24)
simpleFoam是OpenFOAM中使⽤SIMPLE算法的最经典的求解器，在这⼀节我们对其代码进⾏解析。
在代码解析中，我们对OpenFOAM⾃带的注释进⾏翻译，翻译内容以“译：”开头，对于我们⾃⼰的解释内容则以"释："开头。
代码1 simpleFoam代码解读
simpleFoam.C:
#include"fvCFD.H"
//释：OpenFOAM最基础的库⽂件，必须包含
#include"singlePhaTransportModel.H"
//释：导⼊单相流体输运模型库
#include"turbulentTransportModel.H"
//释：导⼊湍流模型库
#include"simpleControl.H"
//释：导⼊SIMPLE算法的控制库
#include"fvOptions.H"
//释：导⼊处理源项的库
//释：上述⽂件为程序运⾏的必要库⽂件，每个库⽂件的作⽤及详细信息⽤户可以查找对应的原代码查看
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
int main(int argc,char*argv[])
{
#include"postProcess.H"
//释：导⼊⽤于后处理相关功能的库
#include"tRootCaLists.H"
//释：导⼊与算例⽬录设置相关的库
#include"createTime.H"
//释：创建时间相关对象
免费资料下载#include"createMesh.H"
//释：创建⽹格相关对象
#include"createControl.H"
装修基础知识//释：创建算法控制相关对象
#include"createFields.H"
//释：创建并初始化速度场，压⼒场等变量。
#include"initContinuityErrs.H"
//释：初始化计算连续性⽅程模拟误差的相关代码
turbulence->validate();
/
/湍流模型相关功能
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
Info<<"\nStarting time loop\n"<< endl;
怪诞心理学while(simple.loop(runTime))
//SIMPLE算法主循环
{
Info<<"Time = "<< runTime.timeName()<< nl << endl;
// --- Pressure-velocity SIMPLE corrector
//译：SIMPLE压⼒速度修正代码
{
#include"UEqn.H"
/
/求解速度⽅程，这部分代码之后详解
#include"pEqn.H"
//求解压⼒⽅程，这部分代码之后详解
}
turbulence->correct();
//输运模型与湍流模型的计算修正代码
//输运模型与湍流模型的计算修正代码
runTime.write();
//写出需要输出的场
Info<<"ExecutionTime = "<< runTime.elapdCpuTime()<<" s"
<<" ClockTime = "<< runTime.elapdClockTime()<<" s"
<< nl << endl;
}
Info<<"End\n"<< endl;
return0;
}
UEqn.H:
// Momentum predictor
//译：动量预测代码
//⽤来处理旋转坐标系下的速度修正
tmp<fvVectorMatrix> tUEqn
(
fvm::div(phi, U)
+ MRF.DDt(U)
+ turbulence->divDevReff(U)
//释： turbulence->divDevReff(U)即扩散项，因为湍流模型，流变模型等代码实际上是对扩散项进⾏修正，因此在OpenFOAM中扩散项的计算归于湍流模型的代码中。
==
fvOptions(U)
);
fvVectorMatrix& UEqn = f();
/
/释：上述代码离散了不包含压⼒的动量⽅程，本质上是组建了式(2.8)的矩阵系统，为之后通过式(2.9)和式(2.10)计算HbyA做准备。
// 释：速度⽅程进⾏亚松弛
umPredictor())
{
solve(UEqn ==-fvc::grad(p));
}
//释：求解动量预测⽅程(2.13),umPredictor()将读取算法配置字典中的momentumPredictor选项，如果不设置该项，默认为true，即进⾏动量预测步骤
pEqn.H
{
volScalarField rAU(1.0/UEqn.A());
//释：计算1/aP，其中UEqn.A()获得Ap，从矩阵的⾓度看即是由系数矩阵的对⾓线元素组成的矩阵
volVectorField HbyA(constrainHbyA(rAU*UEqn.H(), U, p));
surfaceScalarField phiHbyA("phiHbyA", fvc::flux(HbyA));
//释：计算HbyA，即式(2.9)⾄(2.10)的过程
MRF.makeRelative(phiHbyA);
//释： MRF为旋转坐标系相关内容
adjustPhi(phiHbyA, U, p);
tmp<volScalarField>rAtU(rAU);
