[kaifu]

Hi Foamers,

I create a new steady state solver based on the standard buoyantSimpleFoam to study natural convection. Part of solid-fluid zones are modelled with porous media. Thus I add some features from the standard solver rhoPorousSimpleFoam. The major equation is coded as,
1. UEqn.H

    // Construct the Momentum equation

    MRF.correctBoundaryVelocity(U);

    tmp<fvVectorMatrix> tUEqn
    (
        fvm::div(phi, U)
      + MRF.DDt(rho, U)
      + turbulence->divDevRhoReff(U)
     ==
        fvOptions(rho, U)
    );
    fvVectorMatrix& UEqn = tUEqn.ref();

    UEqn.relax();
    

    // Include the porous media resistance and solve the momentum equation
    // either implicit in the tensorial resistance or transport using by
    // including the spherical part of the resistance in the momentum diagonal

    tmp<volScalarField> trAU;
    tmp<volTensorField> trTU;

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

        fvOptions.constrain(UEqn);

        volVectorField gradp(fvc::grad(p));

        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())  // modify coeffs in rAU
        {
            solve
            (
                UEqn
             ==
                fvc::reconstruct
                (
                    (
            // 2019-09-16
            //          - ghf*fvc::snGrad(rhok)
                      - ghf*fvc::snGrad(rho)
                      - fvc::snGrad(p_rgh)
                    )*mesh.magSf()
                )
            );

            fvOptions.correct(U);
        }

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

2. pEqn.H

{
// 1. p --> p_rgh
// 2. comment implicit part

//    const volScalarField& psi = thermo.psi();

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

//    Info<< "Pre rho max/min : " << max(rho).value() << " " << min(rho).value()
//        << endl;

    tmp<volVectorField> tHbyA;
    if (pressureImplicitPorosity)
    {
//        tHbyA = constrainHbyA(trTU()&UEqn.H(), U, p);
    }
    else
    {
        tHbyA = constrainHbyA(trAU()*UEqn.H(), U, p_rgh);
    }
    volVectorField& HbyA = tHbyA.ref();

    tUEqn.clear();

    surfaceScalarField rhorAUf("rhorAUf", fvc::interpolate(rho*trAU()));
// 2019-09-16
//    surfaceScalarField phig(-rhorAUf*ghf*fvc::snGrad(rhok)*mesh.magSf());
    surfaceScalarField phig(-rhorAUf*ghf*fvc::snGrad(rho)*mesh.magSf());

    bool closedVolume = false;

    surfaceScalarField phiHbyA("phiHbyA", fvc::flux(rho*HbyA));
    MRF.makeRelative(fvc::interpolate(rho), phiHbyA);

    closedVolume = adjustPhi(phiHbyA, U, p_rgh);

    phiHbyA += phig;

    // Update the pressure BCs to ensure flux consistency
    constrainPressure(p_rgh, rho, U, phiHbyA, rhorAUf, MRF);

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

        if (pressureImplicitPorosity)
        {
//            tpEqn =
//            (
//                fvm::laplacian(rho*trTU(), p)
//              + fvOptions(psi, p, rho.name())
//             ==
//                fvc::div(phiHbyA)
//            );
        }
        else
        {
            tp_rghEqn =
            (
                fvm::laplacian(rhorAUf, p_rgh)
//              + fvOptions(psi, p, rho.name())
             ==
                fvc::div(phiHbyA)
            );
        }

        fvScalarMatrix& p_rghEqn = tp_rghEqn.ref();

        p_rghEqn.setReference
        (
            pRefCell, 
            getRefCellValue(p_rgh, pRefCell)
        );

        p_rghEqn.solve();

        if (simple.finalNonOrthogonalIter())
        {
            // Calculate the conservative fluxes
            phi = phiHbyA - p_rghEqn.flux();

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

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

    #include "continuityErrs.H"

// 2019-09-16
//    p = p_rgh + rhok*gh;
    p = p_rgh + rho*gh;
/****
    if (pressureImplicitPorosity)
    {
        U = HbyA - (trTU() & fvc::grad(p));
    }
    else
    {
        U = HbyA - trAU()*fvc::grad(p);
    }

    U.correctBoundaryConditions();
    fvOptions.correct(U);
****/
//    pressureControl.limit(p);

    // For closed-volume cases adjust the pressure and density levels
    // to obey overall mass continuity
/****
    if (!thermo.incompressible() && closedVolume)
    {
        p += (initialMass - fvc::domainIntegrate(psi*p))
            /fvc::domainIntegrate(psi);
        p_rgh = p - rho*gh;
    }
****/
    rho = thermo.rho();     // depends on p and T
    rho.relax();
    Info<< "rho max/min : " << max(rho).value() << " " << min(rho).value()
        << endl;

    // update display var
//    gradP_rgh = fvc::grad(p_rgh);
//    gradRho =  fvc::grad(rho);
//    driveSource = -gradP_rgh - gh_display*gradRho;

}

3.EEqn.H

{
    volScalarField& he = thermo.he();

/****
    dimensionedScalar Cp1
    (
        "liquidHeatCapacity",
        dimEnergy/dimMass/dimTemperature, 
        thermo.Cp()
    );

    volScalarField Cp = thermo.Cp();
****/

    gammaPhi = gammaf*phi;

    volScalarField alphaMixEff
    (
        "alphaMixEff", 
        gamma*turbulence->alphaEff() + (1.0-gamma)*lambda2/Cp1
    );
// h1 -> fluid
// h2 -> solid
// h2 = Cp2 *T = Cp2 * h1/Cp1
// laplacian(alpha1*gamma, h1) + laplacian(alpha2*(1-gamma), h2)
//=laplacian(alpha1*gamma, h1) + laplacian(alpha2*Cp2/Cp1*(1-gamma), h1)
// alpha = conductivity/heat capacity [kg/m/s]
// Cp  [J/kg/K]
// lambda: conductivity  [kg m /s^3 /K]
// alpha = nu / prandtl
// the following he refers to enthalpy of fluid
    fvScalarMatrix EEqn
    (
        fvm::div(gammaPhi, he)
      + (
            he.name() == "e"
          ? fvc::div(gammaPhi, volScalarField("Ekp", 0.5*magSqr(U) + p/rho))
          : fvc::div(gammaPhi, volScalarField("K", 0.5*magSqr(U)))
        )
      - fvm::laplacian(alphaMixEff, he)
     ==
        rho*(U&g)*gamma
      + radiation->Sh(thermo, he)
      + fvOptions(rho, he)
    );

    fvScalarMatrix EEqnFvO
    (
      - fvOptions(rho, he)
    );

    const scalarField& V = mesh.V();
    const scalarField& eqnSource = EEqnFvO.source();

    scalarField& heSource = volHeatSource;
    heSource = eqnSource/V;

//    volHeatSource = EEqnFvO.source()/mesh.V();

//    forAll(mesh.C(), celli)
//    {
//        volHeatSource[celli] = eqnSource[celli]/V[celli];
//    }

    EEqn.relax();

    fvOptions.constrain(EEqn);

    EEqn.solve();

    fvOptions.correct(he);

    thermo.correct();
    radiation->correct();

    Cp1 = thermo.Cp();

// 2019-09-16
//    rhok = rho*(1.0 - beta*(T - TRef));
//    rhok.relax();
}

In brief, the porous model is implemented in an explicit way.

However I found the velocity along the channel axis shows checkerboard feature, which is located in porous zones. However, there is no checkerboard out of the porous zones. And the temperature is good everywhere.

In details, there are a lot of channel (in which the cross section is square) hang in the space. The fluid will flow into the channel from bottom and flow out the channel from top, as shown in the figure. There are spaces between two channels, and we call it as interchannel. The fluid is flow from top to bottom in interchannels, which is also shown in the figure. Besides, there are holes in the channel wall so that the fluid may flow through the holes aside. The channel wall is modelled as baffles. The major difference between channel and interchannel lies in: 1. the fuild in channels is heated by volumetric source and interchannel no vol heat sources; 2. there are solid structure in channels so the channel is modelled as porous media, whereas interchannel is pure fluid.
Uz:

p_rgh:

T:

I also investigated the reason of checkerboard feature in the forum. The problem should be solved by Rhie-Chow interpolation,
make a "pseudo-staggered" version of the pressure equation in co-located grid, as discussed in,
Differences in solution method for pisoFoam and buoyantBoussinesqPisoFoam
Causes of checkerboarding

I swithed the divSchemes from upwind to bounded Gauss limitedLinear but it doesnt help.

I dont understand why velocity field has checkerboard feature only in porous zones. The source of porous resistance is add as

pZones.addResistance(UEqn);

in my case, I selected Darcy Forchheimer submodel, so the source is added in an explicit way into the matrix in DarcyForchheimerTemplates.C as,

        forAll(cells, i)
        {
            const label celli = cells[i];
            const label j = this->fieldIndex(i);
            const tensor D = dZones[j];
            const tensor F = fZones[j];

            AU[celli] += mu[celli]*D + (rho[celli]*mag(U[celli]))*F;
        }

Is there any trouble so that the velocity shows checkerboard feature?