devbox
从前慢现在也慢
这个屌丝很懒，什么也没留下！
关注作者
热门标签
热门文章
当前位置:   article > 正文

SIMPLE、PISO 、PIMPLE算法浅析_piso算法

作者：从前慢现在也慢 | 2024-05-03 11:14:13

赞

踩

piso算法

为了求解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

  • 1
  • 2
  • 3
  • 4
  • 5

主程序：

#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
  • 2
  • 3
  • 4
  • 5
  • 6
  • 7
  • 8
  • 9
  • 10
  • 11
  • 12
  • 13
  • 14
  • 15
  • 16
  • 17
  • 18
  • 19
  • 20
  • 21
  • 22
  • 23
  • 24
  • 25
  • 26
  • 27
  • 28
  • 29
  • 30
  • 31
  • 32
  • 33
  • 34
  • 35
  • 36
  • 37
  • 38
  • 39
  • 40
  • 41
  • 42
  • 43
  • 44
  • 45
  • 46
  • 47
  • 48
  • 49
  • 50
  • 51

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);
    }

  • 1
  • 2
  • 3
  • 4
  • 5
  • 6
  • 7
  • 8
  • 9
  • 10
  • 11
  • 12
  • 13
  • 14
  • 15
  • 16
  • 17
  • 18
  • 19
  • 20
  • 21
  • 22

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);
}


  • 1
  • 2
  • 3
  • 4
  • 5
  • 6
  • 7
  • 8
  • 9
  • 10
  • 11
  • 12
  • 13
  • 14
  • 15
  • 16
  • 17
  • 18
  • 19
  • 20
  • 21
  • 22
  • 23
  • 24
  • 25
  • 26
  • 27
  • 28
  • 29
  • 30
  • 31
  • 32
  • 33
  • 34
  • 35
  • 36
  • 37
  • 38
  • 39
  • 40
  • 41
  • 42
  • 43
  • 44
  • 45
  • 46
  • 47
  • 48
  • 49
  • 50
  • 51

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;
}

  • 1
  • 2
  • 3
  • 4
  • 5
  • 6
  • 7
  • 8
  • 9
  • 10
  • 11
  • 12
  • 13
  • 14
  • 15
  • 16
  • 17
  • 18
  • 19
  • 20
  • 21
  • 22
  • 23
  • 24
  • 25
  • 26
  • 27
  • 28
  • 29
  • 30
  • 31
  • 32
  • 33
  • 34
  • 35
  • 36
  • 37
  • 38
  • 39
  • 40
  • 41
  • 42
  • 43
  • 44
  • 45
  • 46
  • 47
  • 48
  • 49
  • 50
  • 51
  • 52
  • 53
  • 54
  • 55
  • 56
  • 57

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);
}

  • 1
  • 2
  • 3
  • 4
  • 5
  • 6
  • 7
  • 8
  • 9
  • 10
  • 11
  • 12
  • 13
  • 14
  • 15
  • 16
  • 17
  • 18
  • 19
  • 20
  • 21

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);

  • 1
  • 2
  • 3
  • 4
  • 5
  • 6
  • 7
  • 8
  • 9
  • 10
  • 11
  • 12
  • 13
  • 14
  • 15
  • 16
  • 17
  • 18
  • 19
  • 20
  • 21
  • 22
  • 23
  • 24
  • 25
  • 26
  • 27
  • 28
  • 29
  • 30
  • 31
  • 32
  • 33
  • 34
  • 35
  • 36
  • 37
  • 38
  • 39
  • 40
  • 41

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;
}

  • 1
  • 2
  • 3
  • 4
  • 5
  • 6
  • 7
  • 8
  • 9
  • 10
  • 11
  • 12
  • 13
  • 14
  • 15
  • 16
  • 17
  • 18
  • 19
  • 20
  • 21
  • 22
  • 23
  • 24
  • 25
  • 26
  • 27
  • 28
  • 29
  • 30
  • 31
  • 32
  • 33
  • 34
  • 35
  • 36
  • 37
  • 38
  • 39
  • 40
  • 41
  • 42
  • 43
  • 44
  • 45
  • 46
  • 47
  • 48
  • 49
  • 50
  • 51
  • 52
  • 53
  • 54
  • 55
  • 56
  • 57
  • 58
  • 59
  • 60
  • 61
  • 62
  • 63
  • 64
  • 65
  • 66
  • 67
  • 68

速度方程：

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);
}

  • 1
  • 2
  • 3
  • 4
  • 5
  • 6
  • 7
  • 8
  • 9
  • 10
  • 11
  • 12
  • 13
  • 14
  • 15
  • 16
  • 17
  • 18
  • 19
  • 20
  • 21
  • 22

压力方程：

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);

  • 1
  • 2
  • 3
  • 4
  • 5
  • 6
  • 7
  • 8
  • 9
  • 10
  • 11
  • 12
  • 13
  • 14
  • 15
  • 16
  • 17
  • 18
  • 19
  • 20
  • 21
  • 22
  • 23
  • 24
  • 25
  • 26
  • 27
  • 28
  • 29
  • 30
  • 31
  • 32
  • 33
  • 34
  • 35
  • 36
  • 37
  • 38
  • 39
  • 40
  • 41
  • 42
  • 43
  • 44
  • 45
  • 46
  • 47
  • 48
  • 49
  • 50
  • 51
  • 52
  • 53
  • 54
  • 55
  • 56
  • 57
  • 58
  • 59

注：最后一次外循环对压力方程不采用亚迭代方式，以便采用标准的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

声明：本文内容由网友自发贡献，不代表【wpsshop博客】立场，版权归原作者所有，本站不承担相应法律责任。如您发现有侵权的内容，请联系我们。转载请注明出处：https://www.wpsshop.cn/w/从前慢现在也慢/article/detail/529345
推荐阅读
相关标签

Copyright © 2003-2013 www.wpsshop.cn 版权所有，并保留所有权利。

  

闽ICP备14008679号