向DPMFoam或MPPICFoam求解器中添加源項

OpenFOAM-2.3.x版本中的DPMFoam與MPPICFoam求解器無法對使用者自定義源操作（即fvOptions）進行處理，需要修改DPMFoam以及MPPICFoam求解器的源碼，以實作其對fvOptions的處理。

在對OpenFOAM中給定的算例進行分析時發現，

$FOAM_TUTORIALS/incompressible/pimpleFoam

算例的

而DPMFoam及MPPICFoam中流體求解部分采取的是PIMPLE方法，是以可以借鑒pimpleFoam求解器對于fvOptions的處理，并對照pimpleFoam求解器的源碼将fvOptions加入到DPMFoam以及MPPICFoam求解器中。

求解器分析

實際上，MPPICFoam求解器與DPMFoam求解器基本一緻，隻是對顆粒的部分處理方式不一樣，這從MPPICFoam求解器源碼也可以看出：

// MPPICFoam.C
#define MPPIC

#include "DPMFoam.C"

是以，對DPMFoam求解器的修改也就自然包含了對MPPICFoam求解器的修改。DPMFoam求解器的主程式如下：

// DPMFoam.C

#include "fvCFD.H"
#include "singlePhaseTransportModel.H"
#include "PhaseIncompressibleTurbulenceModel.H"
#include "pimpleControl.H"
#include "fixedFluxPressureFvPatchScalarField.H"

#ifdef MPPIC
    #include "basicKinematicMPPICCloud.H"
    #define basicKinematicTypeCloud basicKinematicMPPICCloud
#else
    #include "basicKinematicCollidingCloud.H"
    #define basicKinematicTypeCloud basicKinematicCollidingCloud
#endif

int main(int argc, char *argv[])
{
    argList::addOption
    (
        "cloudName",
        "name",
        "specify alternative cloud name. default is 'kinematicCloud'"
    );

    #include "setRootCase.H"
    #include "createTime.H"
    #include "createMesh.H"
    #include "readGravitationalAcceleration.H"
    #include "createFields.H"
    #include "initContinuityErrs.H"

    pimpleControl pimple(mesh);

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

    while (runTime.run())
    {
        #include "readTimeControls.H"
        #include "CourantNo.H"
        #include "setDeltaT.H"

        runTime++;

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

        continuousPhaseTransport.correct();
        muc = rhoc*continuousPhaseTransport.nu();

        Info<< "Evolving " << kinematicCloud.name() << endl;
        kinematicCloud.evolve();

        // Update continuous phase volume fraction field
        alphac = max(1.0 - kinematicCloud.theta(), alphacMin);
        alphac.correctBoundaryConditions();
        alphacf = fvc::interpolate(alphac);
        alphaPhic = alphacf*phic;

        fvVectorMatrix cloudSU(kinematicCloud.SU(Uc));
        volVectorField cloudVolSUSu
        (
            IOobject
            (
                "cloudVolSUSu",
                runTime.timeName(),
                mesh
            ),
            mesh,
            dimensionedVector
            (
                "0",
                cloudSU.dimensions()/dimVolume,
                vector::zero
            ),
            zeroGradientFvPatchVectorField::typeName
        );

        cloudVolSUSu.internalField() = -cloudSU.source()/mesh.V();
        cloudVolSUSu.correctBoundaryConditions();
        cloudSU.source() = vector::zero;

        // --- Pressure-velocity PIMPLE corrector loop
        while (pimple.loop())
        {
            #include "UcEqn.H"

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

            if (pimple.turbCorr())
            {
                continuousPhaseTurbulence->correct();
            }
        }

        runTime.write();

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

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

    return 0;
}

可以發現其中并不包含對fvOptions的處理。此外可以發現DPMFoam求解器在對流場進行求解時使用的是PIMPLE算法，而OpenFOAM的不可壓縮流體求解器pimpleFoam是可以對fvOptions進行處理的，因而可以借鑒pimpleFoam源碼中對fvOptions的處理方式，對DPMFoam源碼進行修改。

pimpleFoam求解器的主程式如下，程式中為了處理fvOptions所添加的頭檔案已認證注釋标示出：

// pimpleFoam.C

#include "fvCFD.H"
#include "singlePhaseTransportModel.H"
#include "turbulenceModel.H"
#include "pimpleControl.H"

// !!! Here
#include "fvIOoptionList.H"
#include "IOporosityModelList.H"
#include "IOMRFZoneList.H"

#include "fixedFluxPressureFvPatchScalarField.H"

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

int main(int argc, char *argv[])
{
    #include "setRootCase.H"
    #include "createTime.H"
    #include "createMesh.H"
    #include "createFields.H"

    // !!! Here
    #include "createFvOptions.H"

    #include "initContinuityErrs.H"

    pimpleControl pimple(mesh);

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

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

    while (runTime.run())
    {
        #include "readTimeControls.H"
        #include "CourantNo.H"
        #include "setDeltaT.H"

        runTime++;

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

        // --- Pressure-velocity PIMPLE corrector loop
        while (pimple.loop())
        {
            #include "UEqn.H"

            // --- Pressure corrector loop
            while (pimple.correct())
            {
                #include "pEqn.H"
            }

            if (pimple.turbCorr())
            {
                turbulence->correct();
            }
        }

        runTime.write();

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

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

    return 0;
}

此外，PIMPLE算法循環中的速度方程

UEqn.H

以及壓力方程

pEqn.H

中均對fvOptions進行了處理，是以在對DPMFoam求解器進行修改的時候，也需要考慮速度方程以及壓力方程的修改。

求解器修改

首先需要将OpenFOAM初始的DPMFoam求解器備份：

cd $FOAM_APP/solvers/lagrangian
cp -rfv DPMFoam DPMFoam-bak

然後再修改DPMFoam求解器。

`DPMFoam.C`

對照

pimpleFoam.C

，可以初步将fvOptions相關的頭檔案加入到

DPMFoam.C

中，如下：

// DPMFoam.C

#include "fvCFD.H"
#include "singlePhaseTransportModel.H"
#include "PhaseIncompressibleTurbulenceModel.H"
#include "pimpleControl.H"
#include "fvIOoptionList.H"
#include "IOporosityModelList.H"
#include "IOMRFZoneList.H"
#include "fixedFluxPressureFvPatchScalarField.H"

#ifdef MPPIC
    #include "basicKinematicMPPICCloud.H"
    #define basicKinematicTypeCloud basicKinematicMPPICCloud
#else
    #include "basicKinematicCollidingCloud.H"
    #define basicKinematicTypeCloud basicKinematicCollidingCloud
#endif

int main(int argc, char *argv[])
{
    argList::addOption
    (
        "cloudName",
        "name",
        "specify alternative cloud name. default is 'kinematicCloud'"
    );

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

    pimpleControl pimple(mesh);

    // ......
}

`UEqn.H`

對照pimpleFoam中包含的

UEqn.H

檔案：

// Solve the Momentum equation

tmp<fvVectorMatrix> UEqn
(
    fvm::ddt(U)
  + fvm::div(phi, U)
  + turbulence->divDevReff(U)
 ==
    fvOptions(U)
);

UEqn().relax();

fvOptions.constrain(UEqn());

volScalarField rAU(1.0/UEqn().A());

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

    fvOptions.correct(U);
}

可以将DPMFoam中包含的

UEqn.H

檔案修改為：

// fvOptions is added in the UcEqn.
fvVectorMatrix UcEqn
(
    fvm::ddt(alphac, Uc) + fvm::div(alphaPhic, Uc)
  - fvm::Sp(fvc::ddt(alphac) + fvc::div(alphaPhic), Uc)
  + continuousPhaseTurbulence->divDevRhoReff(Uc)
 ==
    (1.0/rhoc)*cloudSU
  + fvOptions(Uc)
);

UcEqn.relax();

// fvOptions is added here.
fvOptions.constrain(UcEqn);

volScalarField rAUc(1.0/UcEqn.A());
surfaceScalarField rAUcf("Dp", fvc::interpolate(rAUc));

surfaceScalarField phicForces
(
   (fvc::interpolate(rAUc*cloudVolSUSu/rhoc) & mesh.Sf())
 + rAUcf*(g & mesh.Sf())
);

if (pimple.momentumPredictor())
{
    solve
    (
        UcEqn
     ==
        fvc::reconstruct
        (
            phicForces/rAUcf - fvc::snGrad(p)*mesh.magSf()
        )
    );

    // fvOptions is added here.
    fvOptions.correct(Uc);
}

`pEqn.H`

對照pimpleFoam中包含的

pEqn.H

檔案：

surfaceScalarField rAUf("rAUf", fvc::interpolate(rAU));

volVectorField HbyA("HbyA", U);
HbyA = rAU*UEqn().H();

if (pimple.nCorrPISO() <= 1)
{
    UEqn.clear();
}

surfaceScalarField phiHbyA
(
    "phiHbyA",
    (fvc::interpolate(HbyA) & mesh.Sf())
  + rAUf*fvc::ddtCorr(U, phi)
);

fvOptions.makeRelative(phiHbyA);

adjustPhi(phiHbyA, U, p);

// Update the fixedFluxPressure BCs to ensure flux consistency
setSnGrad<fixedFluxPressureFvPatchScalarField>
(
    p.boundaryField(),
    (
        phiHbyA.boundaryField()
      - fvOptions.relative(mesh.Sf().boundaryField() & U.boundaryField())
    )/(mesh.magSf().boundaryField()*rAUf.boundaryField())
);

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

    pEqn.setReference(pRefCell, pRefValue);

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

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

可以将DPMFoam中包含的

pEqn.H

檔案修改為：

{
    volVectorField HbyA("HbyA", Uc);
    HbyA = rAUc*UcEqn.H();

    surfaceScalarField phiHbyA
    (
        "phiHbyA",
        (
           (fvc::interpolate(HbyA) & mesh.Sf())
         + alphacf*rAUcf*fvc::ddtCorr(Uc, phic)
         + phicForces
        )
    );

    // fvOptions is added here.
    fvOptions.makeRelative(phiHbyA);

    // fvOptions is added here.
    // Update the fixedFluxPressure BCs to ensure flux consistency
    setSnGrad<fixedFluxPressureFvPatchScalarField>
    (
        p.boundaryField(),
        (
            phiHbyA.boundaryField()
          - fvOptions.relative(mesh.Sf().boundaryField() & Uc.boundaryField())
        )/(mesh.magSf().boundaryField()*rAUcf.boundaryField())
    );

    // Non-orthogonal pressure corrector loop
    while (pimple.correctNonOrthogonal())
    {
        fvScalarMatrix pEqn
        (
            fvm::laplacian(alphacf*rAUcf, p)
         ==
            fvc::ddt(alphac) + fvc::div(alphacf*phiHbyA)
        );

        pEqn.setReference(pRefCell, pRefValue);

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

        if (pimple.finalNonOrthogonalIter())
        {
            phic = phiHbyA - pEqn.flux()/alphacf;

            p.relax();

            Uc = HbyA
              + rAUc*fvc::reconstruct((phicForces - pEqn.flux()/alphacf)/rAUcf);
            Uc.correctBoundaryConditions();

            // fvOptions is added here.
            fvOptions.correct(Uc);
        }
    }
}

#include "continuityErrs.H"

其他檔案修改

`Make/options`

在該檔案中的

EXE_INC

選項中添加

-I$(LIB_SRC)/fvOptions/lnInclude

及

-I$(LIB_SRC)/sampling/lnInclude

。在該檔案中的

EXE_LIBS

選項中添加

-lfvOptions

及

-lsamping

。這裡要注意每行行尾的轉義字元

。

添加這些選項的原因是使用OpenFOAM自帶的編譯工具wmake進行編譯時，可以通過

EXE_INC

及

EXE_LIBS

來尋找需要包含的源檔案以及需要連結的庫檔案。增加fvOptions相關的源檔案及庫檔案路徑可以保證編譯的順利進行。

`MPPICFoam/make/options`

該檔案的修改方式與

Make/options

檔案的修改方式相同。

編譯

在DPMFoam目錄中執行如下指令：

./Allwclean
./Allwmake

即可以實作修改後求解器的編譯。

向DPMFoam或MPPICFoam求解器中添加源項

求解器分析

求解器修改

`DPMFoam.C`

`UEqn.H`

`pEqn.H`

其他檔案修改

`Make/options`

`MPPICFoam/make/options`

編譯

繼續閱讀

Regarding modified drag force formulation for application in dense flows

ANSYS FLUENT三維非結構網格數值計算及後處理——教室全新風送風與排風

openFOAM中的tensor

在OpenFOAM中使用foamCloneCase快速建立參數化研究所需的各算例

一維熱傳導問題與C++描述1一維熱傳導問題與C++描述

OPEN AI-Chat GPT& CFD逆天回答（一）

生物流體力學及血流動力學模組化仿真技術實戰

openfoam v8 波浪算例學習日記： 5.邊界條件及初始條件

CentOS6.5 安裝 OpenFOAM-2.4.0

【計算流體力學CFD】Fluent軟體模拟：方腔熱對流&圓柱繞流（卡門渦街）｜Matlab

OpenFOAM類庫介紹（一）網格和場Part1 網格Part2 場

openfoam學習心得-dynamic mesh_day1

openFoam4.0 學習記錄（七）：Tutorials之Lid-driven Cavity Flow (後續3-高雷諾數RAS計算）

openFoam4.0 學習記錄（六）：Tutorials之Lid-driven Cavity Flow (後續2）

openFoam4.0 學習記錄（四）：Tutorials 之 Stress analysis of a plate with a hole （後續）

OpenFOAM大渦模拟湍流模型之Smagorinsky模型代碼詳解