OpenFOAM基础算法 SIMPLE & PISO & PIMPLE

Created time
Jul 10, 2025 11:11 AM
type
Post
status
Published
date
Jul 10, 2025
slug
simple-piso-pimple
summary
tags
2024
ofs
ofsp2024
OpenFOAM
category
post
icon
password
Place
Last edited time
Mar 27, 2026 01:07 AM
😀
前言： 访问 https://aerosand.cc 以获取最近更新。
 
参考：
感谢原作者们的无私引路和宝贵工作。
我们首先回忆一下基本方程 （CFD理论基础00 基本方程 | 𝓐𝓮𝓻𝓸𝓼𝓪𝓷𝓭 (aerosand.cn)）
基本方程为
  1. 质量方程
  1. 动量方程
  1. 能量方程
通用形式
实际问题可能要求解多个方程，而方程组的求解，即使是数值求解也很困难。

1. 压力速度修正算法

对于不可压流动，密度恒定，压力速度耦合控制方程为（压力是除以密度的相对压力）
关于方程组
  • 对于三维问题来说，四个分量方程对应四个未知数（），看似可以求解
  • 但是压力没有自己的约束方程（不可压流动没有关于压力的状态方程）
  • 质量方程实际上是对动量方程求解后的再约束，也就是说动量方程求解后的速度场仍然要回代以满足质量方程
在阅读过 FVM 理论之后，可以知道，本地离散方程总是
系数 * 本单元量 + 系数 * 邻单元量 = 本地源项
遍历全局后，组建全局离散方程，形成待求解的线性代数系统，形如
与待求场量相关的统统的纳入系数矩阵，放在左边，也就是 LHS（Left-hand-side），涉及源项线性化之后与场不直接相关的项等等，统统放在右手边，也就是 RHS(Right-hand-side)。最终的最终，动量方程形式可以为
这里没有考虑因为离散或者因为边界条件等原因出现的 RHS
可以想见场的系数矩阵 M 是一个对角占优的稀疏方阵（或者说我们希望它能够对角占优），$M$ 矩阵的对角线都是离散方程中的本元素（有的书用字母 C 指代，有的用字母 P 指代），同一行上，对角线前后（在该行也许不紧挨）的元素都是和该单元相邻的单元。
最基本的思路是，我们最开始假设一个初始压力场，根据动量方程求出速度场。直接从动量方程求出来的速度场反过来再去满足质量方程的约束。
那，怎么从动量方程求出速度场呢？
为了更好的求解方程（3），先处理系数矩阵 $M$ 。
从 $M$ 矩阵中分解出对角矩阵 $A$ ，对角矩阵 $A$ 可以很容易求出逆矩阵。OpenFOAM 中的对角矩阵以及求解方法已经高度抽象化，也更易读。
抽离得到对角矩阵后，方程（3）的左侧可以写成
作为比较，注意以下处理是错误的。
$$ \cancel{MU = AU - HU} $$
所以，分解后的动量方程为
两边同时乘以 $A$ 的逆矩阵
可以求出速度场
其中
$H$ 的求解在 OpenFOAM 直接使用成员函数调用即可 UEqn.H() 。
到这里速度场就求解出来了，OpenFOAM 把方程求解高度抽象化，但也更加容易阅读。
动量方程求解出速度的相关代码在文件 UEqn.H 中，这一步也被称为 momentum predictor。
求出的速度场并不是最终结果，速度场需要满足质量方程，所以还要求解质量方程。 质量方程为
代入速度场的解，也就是方程（7），有
整理后得到实际用来求解的速度约束方程，
方程（9）也被称为压力方程，OpenFOAM 中的相关代码在 pEqn.H 中
求解压力方程（9），可以得到新的压力场。
压力速度修正计算的大体思路就是这样子，下面具体看不同的求解算法。

2. SIMPLE

2.1. 理论

  1. 基于上一步的压力场（启动步使用初始压力场），计算 momentum predictor ，求得速度场
  1. 求出非对角线矩阵 $H$
  1. 求解质量方程，求出新的压力场
  1. 以新的压力场，更新速度场
  1. 回到步骤 1，基于新的压力场，求解速度场，继续迭代
这一过程被称为 outer loops

2.2. 实现

我们对比 OpenFOAM 原生 simpleFoam 求解器，看看单纯的 SIMPLE 算法的实现应该是什么样子的。
  1. 检查是否继续循环——simple.loop()
  1. 使用 momentum predictor 求解速度——UEqn.H（也就是上面的 step1）
  1. 修正压力和速度——pEqn.H （也就是上面的 step3，step4）
  1. 为湍流模型求解输运方程——turbulence->correct()
  1. 返回步骤 1
Momentum predictor 也就是 UEqn.H 中的主要代码（因版本有变化只摘取主要部分）
之后进行压力修正，基于上一步的速度求解（predicted velocity）。修正后的压力被拿来求解质量方程，也就修正了速度。pEqn.H 代主要码如下（因版本有变化只摘取主要部分）
主要算法框架主要结构如下

2.3. 讨论

OpenFOAM 经过多年的更新，求解器增添了很多保证数值计算结果的方法。我们暂时不用深究，只需要把握主要算法即可。

3. PISO

3.1. 理论

  1. 基于上一步的压力场（启动步使用初始压力场），计算 momentum predictor ，求得速度场
  1. 求出非对角线参与矩阵
  1. 求解质量方程，求出新的压力场
  1. 以新的压力场，更新速度场
  1. 回到步骤 2（不再重新求解原始动量方程），继续迭代
  1. 求解满足要求后，时间推进
Loading...