[shock77]

Hi,


I have noticed that rhoCentralFoam has no porosity treatment unlike sonicFoam or rhoPimpleFoam. Since I have the best results with rhoCentralFoam I would like to add the porosity treatment.


I want to use the Darcy-Forchheimer model explicitly by using fvOptions, which works well for other solvers except rhoCentralFoam.


I have noticed, that in rhoPimpleFoam UEqn.H fvOptions is mentioned:

// Solve the Momentum equation

MRF.correctBoundaryVelocity(U);

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

UEqn.relax();

fvOptions.constrain(UEqn);

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

    fvOptions.correct(U);
    K = 0.5*magSqr(U);
}

while in rhoCentralFoam.C it is not:

        // --- Solve momentum
        solve(fvm::ddt(rhoU) + fvc::div(phiUp));

        U.ref() =
            rhoU()
           /rho();
        U.correctBoundaryConditions();
        rhoU.boundaryFieldRef() == rho.boundaryField()*U.boundaryField();

        if (!inviscid)
        {
            solve
            (
                fvm::ddt(rho, U) - fvc::ddt(rho, U)
              - fvm::laplacian(muEff, U)
              - fvc::div(tauMC)
            );
            rhoU = rho*U;
        }

Since my flow is viscous, I should concentrate on the if statement. Replacing it by

        if (!inviscid)
        {
            solve
            (
                fvm::ddt(rho, U) - fvc::ddt(rho, U)
              - fvm::laplacian(muEff, U)
              - fvc::div(tauMC) ==    fvOptions(rho, U)
            );
            rhoU = rho*U;
        }

doesn't seem to work.




I would be glad to hear your advice!




Regards,
shock77