SIMPLE、PISO 、PIMPLE算法浅析

最新推荐文章于 2024-05-24 20:21:30 发布
最新推荐文章于 2024-05-24 20:21:30 发布
阅读量1.9w 收藏 114
点赞数 24
分类专栏： pimple、piso、simple 文章标签： OpenFOAM
hty
本文链接：https://blog.csdn.net/hanbingchegu/article/details/107136600
版权

为了求解Navier-Stokes方程，往往需要解耦速度场和压力场的耦合方程组。 OpenFOAM 的求解器大多使用 PISO 或者 SIMPLE 算法，或者二者的结合体 PIMPLE 算 法。这些算法一般采用预测-校正策略，通过迭代的方式将速度和压力的计算解耦。

PISO 以及 PIMPLE 算法用来处理非稳态问题，SIMPLE 算法用来处理稳态问题。

以下程序均来自于$FOAM_APP/solvers/incompressible/

文章目录

1.SIMPLE 算法

SIMPLE 算法流程图
上图为SIMPLE 算法的流程图，其求解步骤一般如下:

1.检查是否达到收敛-simple.loop()
2.使用动量预测器-UEqn.H 预测速度
3.校正压力和速度–pEqn.H
4.求解湍流模型的传输方程–湍流->校正（）
5.返回步骤1

主程序：

#include "fvCFD.H"
#include "singlePhaseTransportModel.H"
#include "turbulentTransportModel.H"
#include "simpleControl.H"
#include "fvOptions.H"

// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //

int main(int argc, char *argv[])
{
    #include "postProcess.H"

    #include "setRootCase.H"
    #include "createTime.H"
    #include "createMesh.H"
    #include "createControl.H"
    #include "createFields.H"
    #include "createFvOptions.H"
    #include "initContinuityErrs.H"
     //进行湍流模式判定
    turbulence->validate();

    // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //

    Info<< "\nStarting time loop\n" << endl;

    while (simple.loop())
    {
        Info<< "Time = " << runTime.timeName() << nl << endl;

        // --- Pressure-velocity SIMPLE corrector
        {
            #include "UEqn.H"
            #include "pEqn.H"
        }
        //重新计算湍动黏性，也就是对turbulence->divDevReff(U)需要的量进行更新。
        //比如k－e模型中下面函数包括求解k－e方程和计算有效黏性系数。
        laminarTransport.correct();
        turbulence->correct();

        runTime.write();

        Info<< "ExecutionTime = " << runTime.elapsedCpuTime() << " s"
            << "  ClockTime = " << runTime.elapsedClockTime() << " s"
            << nl << endl;
    }

    Info<< "End\n" << endl;

    return 0;
}

1.动量预测：（该代码来自`simpleFoam 求解器）

    MRF.correctBoundaryVelocity(U);
//利用tmp构建速度方程
    tmp<fvVectorMatrix> tUEqn
    (
        fvm::div(phi, U)
      + MRF.DDt(U)
      + turbulence->divDevReff(U)
     ==
        fvOptions(U)
    );
    fvVectorMatrix& UEqn = tUEqn.ref();

    UEqn.relax();

    fvOptions.constrain(UEqn);

    if (simple.momentumPredictor())
    {
        solve(UEqn == -fvc::grad(p));

        fvOptions.correct(U);
    }

2.校正

根据求得的预测速度校正压力场，再根据连续性方程，校正获得新的速度场。

{
    volScalarField rAU(1.0/UEqn.A());
    volVectorField HbyA(constrainHbyA(rAU*UEqn.H(), U, p));
    surfaceScalarField phiHbyA("phiHbyA", fvc::flux(HbyA));
    MRF.makeRelative(phiHbyA);
    adjustPhi(phiHbyA, U, p);
//引用tmp指针
    tmp<volScalarField> rAtU(rAU);

    if (simple.consistent())
    {
        rAtU = 1.0/(1.0/rAU - UEqn.H1());
        phiHbyA +=
            fvc::interpolate(rAtU() - rAU)*fvc::snGrad(p)*mesh.magSf();
        HbyA -= (rAU - rAtU())*fvc::grad(p);
    }
//清理tmp类型，减少峰值内存
    tUEqn.clear();

    // Update the pressure BCs to ensure flux consistency
    constrainPressure(p, U, phiHbyA, rAtU(), MRF);

    // Non-orthogonal pressure corrector loop
    while (simple.correctNonOrthogonal())
    {
        fvScalarMatrix pEqn
        (
            fvm::laplacian(rAtU(), p) == fvc::div(phiHbyA)
        );

        pEqn.setReference(pRefCell, pRefValue);

        pEqn.solve();//亚松弛

        if (simple.finalNonOrthogonalIter())
        {
            phi = phiHbyA - pEqn.flux();
        }
    }

    #include "continuityErrs.H"
//压力亚松驰
    // Explicitly relax pressure for momentum corrector
    p.relax();

    // Momentum corrector
    U = HbyA - rAtU()*fvc::grad(p);
    U.correctBoundaryConditions();
    fvOptions.correct(U);
}

2.PISO算法

PISO算法也采用预测—校正策略，并对速度场进行预测。然后根据事先设定的迭代次数进行压力场和速度场的修正。最后，求解湍流模型的传递函数。
PISO算法流程图

  • SIMPLE算法对压力方程进行松弛迭代p.relax()，而PISO算法则没有
  • SIMPLE算法解决的是稳态问题，因此每个内循环压力只修正了一次，而PISO算法压力方程需要多次修正。
  • SIMPLE算法速度方程中含有时间离散项，并且使用了tmp类型，而PISO算法则没有，这主要是因为在PISO算法中，在每个内循环中，压力方程进行了多次修正，每一次都需要更新速度矩阵系数，因此UEqn不能被释放。但是在SIMPLE算法中，在一个循环内，压力只修正了一次，因此在构造HbyA 项后，UEqn 即可以被释放以减少峰值内存（peak memory）。
#include "fvCFD.H"
#include "singlePhaseTransportModel.H"
#include "turbulentTransportModel.H"
#include "pisoControl.H"
#include "fvOptions.H"

// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //

int main(int argc, char *argv[])
{
    #include "postProcess.H"

    #include "setRootCase.H"
    #include "createTime.H"
    #include "createMesh.H"
    #include "createControl.H"
    #include "createFields.H"
    #include "createFvOptions.H"
    #include "initContinuityErrs.H"

    turbulence->validate();

    // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //

    Info<< "\nStarting time loop\n" << endl;

    while (runTime.loop())
    {
        Info<< "Time = " << runTime.timeName() << nl << endl;

        #include "CourantNo.H"

        // Pressure-velocity PISO corrector
        {
            #include "UEqn.H"

            // --- PISO loop
            while (piso.correct())
            {
                #include "pEqn.H"
            }
        }

        laminarTransport.correct();
        turbulence->correct();

        runTime.write();

        Info<< "ExecutionTime = " << runTime.elapsedCpuTime() << " s"
            << "  ClockTime = " << runTime.elapsedClockTime() << " s"
            << nl << endl;
    }

    Info<< "End\n" << endl;

    return 0;
}

UEqn.H

MRF.correctBoundaryVelocity(U);

fvVectorMatrix UEqn
(
    fvm::ddt(U) + fvm::div(phi, U)
  + MRF.DDt(U)
  + turbulence->divDevReff(U)
 ==
    fvOptions(U)
);

UEqn.relax();

fvOptions.constrain(UEqn);

if (piso.momentumPredictor())
{
    solve(UEqn == -fvc::grad(p));

    fvOptions.correct(U);
}

pEqn.H

volScalarField rAU(1.0/UEqn.A());
volVectorField HbyA(constrainHbyA(rAU*UEqn.H(), U, p));
surfaceScalarField phiHbyA
(
    "phiHbyA",
    fvc::flux(HbyA)
  + fvc::interpolate(rAU)*fvc::ddtCorr(U, phi)
);

MRF.makeRelative(phiHbyA);

adjustPhi(phiHbyA, U, p);

// Update the pressure BCs to ensure flux consistency
constrainPressure(p, U, phiHbyA, rAU, MRF);

// Non-orthogonal pressure corrector loop
while (piso.correctNonOrthogonal())
{
    // Pressure corrector

    fvScalarMatrix pEqn
    (
        fvm::laplacian(rAU, p) == fvc::div(phiHbyA)
    );

    pEqn.setReference(pRefCell, pRefValue);

    pEqn.solve(mesh.solver(p.select(piso.finalInnerIter())));

    if (piso.finalNonOrthogonalIter())
    {
        phi = phiHbyA - pEqn.flux();
    }
}

#include "continuityErrs.H"

U = HbyA - rAU*fvc::grad(p);
U.correctBoundaryConditions();
fvOptions.correct(U);

3.PIMPLE算法

在有限容积离散中，时间项的离散仍采用的差分格式，这样做可以得到某个时间点的流场信息，而非某个时间步长的内的平均值。采用传统的piso算法求解变化较快的流动的时候，需要的时间步长较小（因为相邻两个时间点的流场不能差别太大，否则会发散），这样会造成求解的某种流动需要的耗时过长。
pimple算法将每个时间步长内看成一种稳态的流动（采用亚松驰来解决相邻两个时间段变化大的情况），当按照稳态的求解器求解到一定的时候，则采用标准的piso做最后一步求解。(本段摘录自苏老师的博客，即参考内容1)

PIMPLE算法流程图
由上述流程图，我们可以看出：

  • PIMPLE算法和PISO算法大体相同，但是在时间步内增加了速度压力耦合求解的循环while (pimple.1oop()), 注意PISO算法速度方程组建一次， 多次求解压力。
  • PIMPLE算法可以使用松弛技术。值得注意的是在源代码中pisoFoam存在UEqn.relax(),但是不存在p.relax()。这说明PISO算法中的动量预测完全可以删掉。
  • 时间步内对速度矩阵方程系数进行更新。

主程序代码:

#include "fvCFD.H" //有限体积库头文件的集合
//单相流流动粘性选择器，通过粘性模型指针（私有变量），根据选择的粘性模型更新粘性（viscosity）
//在`correct`内完成更新，通过nu()返回更新后的粘性
#include "singlePhaseTransportModel.H" 
#include "turbulentTransportModel.H"//湍流传递模型
#include "pimpleControl.H"//pimple算法
#include "fvOptions.H"//控制方程的源项或约束

// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //

int main(int argc, char *argv[])
{
    #include "postProcess.H"

    #include "setRootCase.H"
    #include "createTime.H"
    #include "createMesh.H"
    #include "createControl.H"
    #include "createTimeControls.H"
    #include "createFields.H" //一般需要自定义
    #include "createFvOptions.H"
    #include "initContinuityErrs.H"
   //进行湍流模式验证
    turbulence->validate();

    // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //

    Info<< "\nStarting time loop\n" << endl;
    while (runTime.run()) 
    {
        #include "readTimeControls.H"
        #include "CourantNo.H"  //计算并输出库朗数
        #include "setDeltaT.H" //通过库朗数设置时间步长`dealtT`

        runTime++;//时间步进
		//显示当前运行时间
        Info<< "Time = " << runTime.timeName() << nl << endl;
//simple修正——外循环
        // --- Pressure-velocity PIMPLE corrector loop
        while (pimple.loop())      
        {
            #include "UEqn.H"
//piso修正——内循环
            // --- Pressure corrector loop
            while (pimple.correct())
            {
                #include "pEqn.H"
            }

            if (pimple.turbCorr())  //湍流修正
            {
                laminarTransport.correct();
                //重新计算湍动黏性，对turbulence->divDevReff(U)需要的量进行更新。
                turbulence->correct();
            }
        }

        runTime.write();

        Info<< "ExecutionTime = " << runTime.elapsedCpuTime() << " s"
            << "  ClockTime = " << runTime.elapsedClockTime() << " s"
            << nl << endl;
    }

    Info<< "End\n" << endl;

    return 0;
}

速度方程：

MRF.correctBoundaryVelocity(U);
//tmp瞬态对象模板类，详见(https://editor.csdn.net/md/?articleId=107122685)
tmp<fvVectorMatrix> tUEqn
(
    fvm::ddt(U) + fvm::div(phi, U)
  + MRF.DDt(U)
  + turbulence->divDevReff(U)
 ==
    fvOptions(U)
);
fvVectorMatrix& UEqn = tUEqn.ref();

UEqn.relax();  //亚松弛，提高稳定性，加快收敛

fvOptions.constrain(UEqn);

if (pimple.momentumPredictor())//动量预测开关，默认为关闭，这里的p由上一时间步求得
{
    solve(UEqn == -fvc::grad(p));

    fvOptions.correct(U);
}

压力方程：

volScalarField rAU(1.0/UEqn.A());
volVectorField HbyA(constrainHbyA(rAU*UEqn.H(), U, p));
surfaceScalarField phiHbyA
(
    "phiHbyA",
    fvc::flux(HbyA)
  + fvc::interpolate(rAU)*fvc::ddtCorr(U, phi)
);

MRF.makeRelative(phiHbyA);

adjustPhi(phiHbyA, U, p);

tmp<volScalarField> rAtU(rAU);

if (pimple.consistent())
{
    rAtU = 1.0/max(1.0/rAU - UEqn.H1(), 0.1/rAU);
    phiHbyA +=
        fvc::interpolate(rAtU() - rAU)*fvc::snGrad(p)*mesh.magSf();
    HbyA -= (rAU - rAtU())*fvc::grad(p);
}
//nCorrPISO<= 1，压力方程只需修正一次，不需要反复更新速度矩阵系数，可直接释放
if (pimple.nCorrPISO() <= 1)
{
    tUEqn.clear();   //清空场对象
}

// Update the pressure BCs to ensure flux consistency
constrainPressure(p, U, phiHbyA, rAtU(), MRF);

// Non-orthogonal pressure corrector loop
while (pimple.correctNonOrthogonal())
{
    // Pressure corrector
    fvScalarMatrix pEqn
    (
        fvm::laplacian(rAtU(), p) == fvc::div(phiHbyA)
    );

    pEqn.setReference(pRefCell, pRefValue);
	/通过查询system/fvSolution下的PIMPLE更新压力参考网格并重新设定参考值
    pEqn.solve(mesh.solver(p.select(pimple.finalInnerIter())));

    if (pimple.finalNonOrthogonalIter())
    {
        phi = phiHbyA - pEqn.flux();
    }
}
//计算连续性方程误差
#include "continuityErrs.H"

// Explicitly relax pressure for momentum corrector
p.relax();//此处外迭代次数为n,则压力场的松弛仅在n-1次外循环前进行，若外迭代次数为1，这里使用piso算法，压力不进行亚松驰

U = HbyA - rAtU()*fvc::grad(p);
//校正速度，满足边界条件(主要针对第二类边界条件）
U.correctBoundaryConditions();
fvOptions.correct(U);

注：最后一次外循环对压力方程不采用亚迭代方式，以便采用标准的PISO循环做最后一步求解！

参考内容：
1.http://blog.sina.com.cn/s/blog_5fdfa7e60100fbb1.html
2.http://openfoamwiki.net/index.php/SimpleFoam
3.http://openfoamwiki.net/index.php/OpenFOAM_guide/The_PIMPLE_algorithm_in_OpenFOAM

寒冰彻骨
关注
  • 24
    点赞
  • 踩
  • 114
    收藏
    觉得还不错? 一键收藏
  • 3
    评论
专栏目录
PHP实现PIMPLE算法（附完整源码）
希望我的博客，能帮上你解决学习中工作中所遇到的问题
05-02 48
PHP实现PIMPLE算法（附完整源码）
PISO、SIMPLE算法对比分析
寒冰彻骨的博客
07-13 6221
本文对对PIMPLE算法进行具体分析。
2.2 SIMPLE系列算法 | 2.3 PISO算法（OpenFOAM理论笔记系列）
热门推荐
CloudBird07的博客
07-31 1万+
2.2 SIMPLE系列算法 2.2.1标准SIMPLE算法 SIMPLE算法(Semi-Implicit Method for PressureLinked Equations)1最初被设计用来求解稳态问题，即控制方程中不包含瞬态项的计算。按照1.3.3节的约定，我们假设计算开始的时候有初始的压力和速度值Po,Uo⃗P^o,\vec{U^o}Po,Uo,待求的真值为Pn,Un⃗P^n,\vec{U^n}Pn,Un。实际上，速度的初始值可以完全随意给定，因为我们完全首先可以将初始压力带入式2.1解出一个速度
icoFoam源码解析（附测试用的源代码文件）
最新发布
BeanChai的博客
05-24 326
本文基于多篇icoFoam解析的文章，做了一些补充以供参考。 本文给出的源代码中包含了大量测试命令，便于大家边运行边学习，提高学习效率。 为不污染OpenFOAM源代码，单独将icoFoam源码提取出来并重命名为icoFoam_learn，方便大家使用。
非定常流动计算的PISO算法_王彤.pdf
12-04
与SIMPLE 算法比较, PISO 算法有其计算非定常流动的优势。本文针对非交错网格下PISO 算法提出改进, 简化其求解过程, 使其计算步骤统一。斯托克斯第二问题计算结果比较, 说明在非交错网格下改进后PI SO 算法比SIMPLE 算法更适宜于非定常流动的计算, 并可应用于叶轮机械内部动静耦合非定常流动的数值求解中。
OpenFOAM的PISO算法【转载】
asd165654的博客
07-19 407
转载自：http://openfoam.blog.sohu.com/94234375.html 流体力学的控制方程是耦合方程组，形式上体现为连续方程和运动方程的耦合，变量上体现为速度和压强的耦合。在数值计算过程中，同时求解速度和压强的方法称为"同时算法"，这通常比较困难繁琐，为此"分离算法"得到了广泛的应用。 分离算法的思路是建立速度和压强之间的迭代关系，在迭代中去寻求稳定...
2.4 PIMPLE算法 | 2.5 附加显式力的压力速度耦合（OpenFOAM理论笔记系列）
CloudBird07的博客
08-04 6017
2.4 PIMPLE算法 2.4.1速度的非线性耦合 在2.1节中，我们曾对流动系统的压力速度耦合问题进行过讨论，实际上，在流动系统中除了压力和速度的耦合关系，我们还存在着另一种重要的耦合关系，速度的非线性耦合关系。当我们对速度进行求解时（例如求解动量预测方程），速度是未知量，因此对流项∇⋅(U⃗U⃗)\nabla\cdot(\vec U\vec U)∇⋅(UU)是一个未知量乘以未知量的结果，这样会使得原来呈线性的方程变为非线性方程（即未知数的最高次数从1变为2）。直接求解非线性系统的矩阵方程需要调用非线性
Simple算法算例
06-02
Simple算法3个典型算例，详细讲解。二维、结构网格单相流SIMPLE 算法的C语言实现。
SIMPLE算法讲解课程PPT课件
10-14
PISO算法针对SIMPLE系列算法中动量方程和连续性方程修正不同步的问题，通过增加速度修正步，同时满足质量和动量守恒。 在实际应用中，SIMPLE算法因其广泛适用性和简单性而被广泛应用，SIMPLEC则适合于对速度场和...
力学流体simple算法完成.zip
07-03
然而，由于迭代过程，SIMPLE算法的计算效率相对较慢，因此，对于大规模的复杂问题，人们通常会寻求更高效的算法，如SIMPLEC或PISO等改进版本。尽管如此，SIMPLE仍然是学习流体仿真和理解流体动力学基础的宝贵工具。
SIMPLE算法-MATLAB.zip_simple_simple matlab_simple 算法_simple-matlab
07-14
SIMPLE算法虽然在计算效率上不如其他高级算法（如SIMPLEC或PISO），但其概念清晰，易于理解和实现，适合初学者学习CFD的基础知识。通过阅读和分析"SIMPLE算法-MATLAB.txt"，你可以更深入地了解如何在实际编程中应用...
pimple的例子代码，仅供参考
12-05
pimple的例子代码，仅供参考,用C++写的，简单易懂
simple算法
05-18
simple算法，介绍一种改进后的simple算法，计算流体
simple 算法程序
08-25
运行软件： Microsoft PowerStation 4.0. 注意：本程序如果采用 CVF（Compaq Visual FORTRAN）会发生编译出错。这是因为在程序中产生数组越界问题，大概是CVF编译检查比较严格，所以编译无法通过。但是如果在CVF中采用Release 模式编译，编译可以通过。Microsoft PowerStation 不会出现此问题。
SIMPLE算法Matlab代码
01-25
%--------------------------------------------------------% % SIMPLE Algorithm - staggered mesh % Incompressible Lid-driven Flow in Cavity %--------------------------------------------------------%
openFOAM编程日志_单相不可压缩流动2020edition2.5.pdf
03-23
本文档从有限体积离散出发，介绍了CFD算法中最基础的PISO，SIMPLE和PIMPLE算法，并将对应各个算法的求解器代码进行了简单解析。本文档适合openFOAM编程的初学者掌握基础的理论知识。
C++ Pimpl编程技法
上善若水
07-06 5500
C++ Pimple 编程技法
Pimple相关的源码
weixin_34204057的博客
10-20 188
已经有了非常好的Pimple的相关解析，建议先看下：Pimple - 一个简单的 PHP 依赖注入容器读 PHP - Pimple 源码笔记（上）读 PHP - Pimple 源码笔记（下） 这里通过例子补充下核心方法的说明：相关的类型： 服务(类似单例)： $container['session'] = function ($c) ...
OpenFOAM——SIMPLE,PISO,PIMPLE
qq_45682183的博客
02-24 352
OpenFOAM-SIMPLE & PISO&PIMPLE
fluent里PISO方法和simple方法有什么不同
04-20
PISO方法和simple方法是求解流体动力学方程的两种数值方法，其中PISO方法在求解速度场和压力场的过程中引入了一个瞬时压力修正，可以提高数值解的精度和收敛速度，而simple方法则没有这个修正步骤。具体来说，PISO方法使用预估-校正的迭代格式求解流体动力学方程，而simple方法则是一个单步迭代方法。因此，PISO方法的求解速度一般会比simple方法慢，但是PISO方法的数值解更加精确。

“相关推荐”对你有帮助么？

  • 非常没帮助
  • 没帮助
  • 一般
  • 有帮助
  • 非常有帮助
提交
评论 3
添加红包

请填写红包祝福语或标题

个

红包个数最小为10个

元

红包金额最低5元

当前余额3.43元 前往充值 >
需支付：10.00元
成就一亿技术人!
领取后你会自动成为博主和红包主的粉丝 规则
hope_wisdom
发出的红包
实付元
使用余额支付
点击重新获取
扫码支付
钱包余额 0

抵扣说明：

1.余额是钱包充值的虚拟货币，按照1:1的比例进行支付金额的抵扣。
2.余额无法直接购买下载，可以购买VIP、付费专栏及课程。

余额充值