[hiuluom]

Hi everyone,

I made a solver "buoyantPorousSimpleFoam" based on buoyantSimpleFoam OF 2.2.x. I know that buoyantSimpleFoam original can run porous with explicit.
But I want to simulate with implicit option, it will easier to reach convergence and faster than explicit option.

What I did:

In order to write buoyantPorousSimpleFoam, I review code of rhoPorousSimpleFoam and porousInterFoam. porousInterFoam solves p_rgh as the buoyantSimpleFoam but it is only explicit porous. So I prefer rhoPorousSimpleFoam to make code for buoyantPorousSimpleFoam.

My problem

I ran a test case in rhoPorousSimpleFoam and buoyantPorousSimpleFoam which gravity equal zero but the pressure loss both of solvers are difference.

Does someone can give me a hint how to set these files right?

Here is the content of code:

buoyantPorousSimpleFoam.C

#include "fvCFD.H"
#include "rhoThermo.H"
#include "RASModel.H"
#include "radiationModel.H"
#include "simpleControl.H"
#include "fvIOoptionList.H"
#include "IOporosityModelList.H"

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

int main(int argc, char *argv[])
{
    #include "setRootCase.H"
    #include "createTime.H"
    #include "createMesh.H"
    
    simpleControl simple(mesh);
    
    #include "readGravitationalAcceleration.H"
    #include "createFields.H"
    #include "createFvOptions.H"
    #include "createRadiationModel.H"
    #include "createZones.H"
    #include "initContinuityErrs.H"

    

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

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

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

        // Pressure-velocity SIMPLE corrector
        {
            #include "UEqn.H"
            #include "EEqn.H"
            #include "pEqn.H"
        }

        turbulence->correct();

        runTime.write();

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

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

    return 0;
}


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

createZones.H

    IOporosityModelList pZones(mesh);
    Switch pressureImplicitPorosity(false);

    // nUCorrectors used for pressureImplicitPorosity
    int nUCorr = 0;
    if (pZones.active())
    {
        // nUCorrectors for pressureImplicitPorosity
        simple.dict().readIfPresent("nUCorrectors", nUCorr);

        if (nUCorr > 0)
        {
            pressureImplicitPorosity = true;
            Info<< "Using pressure implicit porosity" << endl;
        }
        else
        {
            Info<< "Using pressure explicit porosity" << endl;
        }
    }

pEqn.H

{
    rho = thermo.rho();
    rho.relax();

    volScalarField rAU(1.0/UEqn().A());
    surfaceScalarField rhorAUf("Dp", fvc::interpolate(rho*rAU));

    volVectorField HbyA("HbyA", U);
    if (pressureImplicitPorosity)
    {
        HbyA = trTU() & UEqn().H();
    }
    else
    {
        HbyA = trAU()*UEqn().H();
    }
    UEqn.clear();

    surfaceScalarField phig(-rhorAUf*ghf*fvc::snGrad(rho)*mesh.magSf());

    surfaceScalarField phiHbyA
    (
        "phiHbyA",
        fvc::interpolate(rho)*(fvc::interpolate(HbyA) & mesh.Sf())
    );

    fvOptions.relativeFlux(fvc::interpolate(rho), phiHbyA);

    bool closedVolume = adjustPhi(phiHbyA, U, p_rgh);

    phiHbyA += phig;

    while (simple.correctNonOrthogonal())
    {
       tmp<fvScalarMatrix> tpEqn;

        if (pressureImplicitPorosity)
        {
            tpEqn = (fvm::laplacian(rho*trTU(), p_rgh) == fvc::div(phiHbyA));
        }
        else
        {
            tpEqn = (fvm::laplacian(rho*trAU(), p_rgh) == fvc::div(phiHbyA));
        };

        tpEqn().setReference(pRefCell, getRefCellValue(p_rgh, pRefCell));
        tpEqn().solve();
        
        if (simple.finalNonOrthogonalIter())
        {
            // Calculate the conservative fluxes
            phi = phiHbyA - tpEqn().flux();

            // Explicitly relax pressure for momentum corrector
            p_rgh.relax();

            // Correct the momentum source with the pressure gradient flux
            // calculated from the relaxed pressure
            if (pressureImplicitPorosity)
            {
                U = HbyA + (trTU() & fvc::reconstruct((phig - tpEqn().flux())/rhorAUf));
            }
            else
            {
                U = HbyA + (trAU()*fvc::reconstruct((phig - tpEqn().flux())/rhorAUf));
            }
            U.correctBoundaryConditions();
            fvOptions.correct(U);
        }
    }

    #include "continuityErrs.H"

    p = p_rgh + rho*gh;

    // For closed-volume cases adjust the pressure level
    // to obey overall mass continuity
    if (closedVolume)
    {
        p += (initialMass - fvc::domainIntegrate(psi*p))
            /fvc::domainIntegrate(psi);
        p_rgh = p - rho*gh;
    }

    rho = thermo.rho();
    rho.relax();
    Info<< "rho max/min : " << max(rho).value() << " " << min(rho).value()
        << endl;
}

UEqn.H

    // Solve the Momentum equation

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

    UEqn().relax();
    
    tmp<volScalarField> trAU;
    tmp<volTensorField> trTU;

    // fvOptions.constrain(UEqn());

    
    
    if (pressureImplicitPorosity)
    {
        tmp<volTensorField> tTU = tensor(I)*UEqn().A();
        pZones.addResistance(UEqn(), tTU());
        trTU = inv(tTU());
        trTU().rename("rAU");

        fvOptions.constrain(UEqn());

        volVectorField gradp(fvc::reconstruct((ghf*fvc::snGrad(rho) + fvc::snGrad(p_rgh))*mesh.magSf()));

        for (int UCorr=0; UCorr<nUCorr; UCorr++)
        {
            U = trTU() & (UEqn().H() - gradp);
        }
        U.correctBoundaryConditions();
        
        fvOptions.correct(U);
    }
    else
    {
        pZones.addResistance(UEqn());

        fvOptions.constrain(UEqn());

        // solve(UEqn() == -fvc::grad(p));

        // fvOptions.correct(U);
        
        if (simple.momentumPredictor())
        {
            solve
            (
                UEqn()
             ==
                fvc::reconstruct
                (
                    (
                      - ghf*fvc::snGrad(rho)
                      - fvc::snGrad(p_rgh)
                    )*mesh.magSf()
                )
            );

            fvOptions.correct(U);
        }

        trAU = 1.0/UEqn().A();
        trAU().rename("rAU");
    }

the case in attach file at

https://1drv.ms/u/s!AkLk2AeMilR1g517bxSoT0IHNg7BQA

Any suggestions would be appreciated.

Thanh

[hiuluom]

Hi all,

If I only add resistance in momentum equation, the pressure loss equivalent to rhoPorousSimpleFoam.

    // Solve the Momentum equation

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

    UEqn().relax();
    pZones.addResistance(UEqn());
    fvOptions.constrain(UEqn());

    if (simple.momentumPredictor())
    {
        solve
        (
            UEqn()
         ==
            fvc::reconstruct
            (
                (
                  - ghf*fvc::snGrad(rho)
                  - fvc::snGrad(p_rgh)
                )*mesh.magSf()
            )
        );

        fvOptions.correct(U);
    }

in porosityProperties file

porous
    {
        type            DarcyForchheimer;
       	active          yes;
        cellZone        porous;

        DarcyForchheimerCoeffs
        {
            d                       d [ 0 -2 0 0 0 0 0 ] ( 0 0 0 ); 
            f                       f [ 0 -1 0 0 0 0 0 ] ( 289.7 -1000 -1000); 
            coordinateSystem
            {
                e1                      ( 0 0 -1 );
                e2                      ( 0 1 0 );
             }    
        }       
    }

I must add porous implicit because I had a test porous fvoption which gave me bad result.
If I try my sample case as porous parameter above with rhoPorousSimpleFoam in explicit option that means nUCorrectors = 0, the pressure loss have been never reached convergence at least the mesh at porous zone must be finest.

Best,
Thanh