[AmirBaqa1987]

Hi everyone.

the momentum equation of porousSimpleFoam is :



the last two terms are darcy and Forchheimer resistance in order.UEqn.h and PEqn.h in solver are:

Ueqn.H

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

    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());
        trTU = inv(tTU());
        trTU().rename("rAU");

        volVectorField gradp = fvc::grad(p);

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

        eqnResidual = solve
        (
            UEqn() == -fvc::grad(p)
        ). initialResidual();

        maxResidual = max(eqnResidual, maxResidual);

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

PEqn.H

if (pressureImplicitPorosity)
{
    U = trTU()&UEqn().H();
}
else
{
    U = trAU()*UEqn().H();
}

UEqn.clear();
phi = fvc::interpolate(U) & mesh.Sf();
adjustPhi(phi, U, p);

for (int nonOrth=0; nonOrth<=nNonOrthCorr; nonOrth++)
{
    tmp<fvScalarMatrix> tpEqn;

    if (pressureImplicitPorosity)
    {
        tpEqn = (fvm::laplacian(trTU(), p) == fvc::div(phi));
    }
    else
    {
        tpEqn = (fvm::laplacian(trAU(), p) == fvc::div(phi));
    }

    tpEqn().setReference(pRefCell, pRefValue);
    // retain the residual from the first iteration
    if (nonOrth == 0)
    {
        eqnResidual = tpEqn().solve().initialResidual();
        maxResidual = max(eqnResidual, maxResidual);
    }
    else
    {
        tpEqn().solve();
    }

    if (nonOrth == nNonOrthCorr)
    {
        phi -= tpEqn().flux();
    }
}

#include "continuityErrs.H"

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

if (pressureImplicitPorosity)
{
    U -= trTU()&fvc::grad(p);
}
else
{
    U -= trAU()*fvc::grad(p);
}

U.correctBoundaryConditions();

now I want to add Brinkman Term to momentum equation, I mean:



for this I've just define in createFields.H and modify UEqn.H like below:

UEqn.H

 tmp<fvVectorMatrix> UEqn
    (
            fvm::div(phi, U)
         -  fvm::laplacian(mue,U)
         + turbulence->divDevReff(U)
    );

am I right?
shouldn't I modify PEqn.H?

thanks a lot.
your hints are appreciated .
Best Regards,
Arjang

[AmirBaqa1987]

I'm still waiting for HELP.

[Caio Martins]

Were you able to find a solution? Were you right?