[makaveli_lcf]

Thank you Alberto, I also found it and now trying to understand this expression.
Please see attached file Volumefield_reconstruction.pdf. I am looking for the original reference where this formulation is established. Could you please give some tips?

[alberto]

I do not have a reference for this, but I think you got it.

[akidess]

Alexander, I think it would be nice to have this on the wiki. If you agree, do you have the time to copy it to an article there?

[makaveli_lcf]

Hi Anton!

Ok, I will put it there. And if Alberto doesn't mind, it would look better with his description of the body forces and the pressure gradient treatment technique.

[santiagomarquezd]

Alexander,

Originally Posted by makaveli_lcf

Thank you Alberto, I also found it and now trying to understand this expression.
Please see attached file Attachment 9124. I am looking for the original reference where this formulation is established. Could you please give some tips?

we were chatting about fvc::reconstruct() with Alberto few months ago and I took similar notes of the method, writing out the algorithm in mathematical form. I was looking for a "continuum form" of this operator, I think it's based in a generalized form of Gauss Theorem [like (3.9) to (3.16) in Hrv. thesis] due to a surface defined field is transformed in cell defined one, but I couldn't hack it, maybe these ideas could help you to find the basis.

Regards.

[vgalindo]

Dear Alberto,

thank you very much for this very detailed explanation about the issue how to add correctly a body force to the (incompressible) momentum equation!
I would be grateful if you would still specify the corresponding code for the p equation.

Best regards, Vladimir

Originally Posted by alberto

I will try to give you the basic idea in the case of a generic body force term F, so that your momentum equation reads:

(1) ddt(U) + div(UU) = div(tau) - grad(p)/rho + F

This equation can be written in semi-discrete form as:

(2) A*U = H - grad(p)/rho + F

where the pressure gradient and the force term we want to include in the momentum interpolation is left explicitly out at this point. In other words, only the first three terms (time derivative, unsteady, divergence of the stress tensor) of Eq. 1 are used to define A and H. For example, in an incompressible code, this could read (it might be different, depending on the solver):

     fvm::ddt(U)
  + fvm::div(phi, U)
  + turbulence->divDevReff(U)

Normally the pressure gradient would be treated directly as a source term in an incompressible code (see pisoFoam for example). Here we assume we want to treat it with the "improved approach".

From Eq. 2, interpolating on faces and dotting with the surface area vector S the pressure gradient and the force term, we have:

- snGrad(p)*|S|/rho + F_f . S

where F_f is F interpolated on faces. This in OF corresponds to:

- fvc::snGrad(p)*mesh.magSf()/rhoa + fvc::interpolate(F) & mesh.Sf()

This term can be used as argument of fvc::reconstruct() to add the contribution of the pressure gradient and of the force term to the momentum predictor:

solve
(
      UEqn == fvc::reconstruct
      (
          - fvc::snGrad(p)*mesh.magSf()/rhoa 
          + fvc::interpolate(F) & mesh.Sf()
   )
)

Then OpenFOAM re-computes U as

U = H/A

where A is updated with the predicted value of U. However, keep in mind that H does not directly include the effect of grad(p) and F.

From Eq. 2 we can derive the flux to construct the pressure equation:

phi = (H/A)_f . S - (1/A)_f * snGrad(p) |S|/rho + (1/A)_f*F_f . S

and imposing div(phi) = 0, you obtain the pressure equation, which will include for the first time all the effects.

The flux is then corrected based on the solution of the pressure equation, and the velocity correction is reconstructed from the flux. Only at this point U will "see" the effect of p and F.

If you take a look at VOF solvers (i.e. interFoam), you will notice that the code does not solve for p, but for p_rgh = p + rho*gh. This is another way to treat the gravity to address the weakness of the Rhie-Chow interpolation, but the idea is similar.

I hope this helps, but please let us know if you have more questions

Best,

[dl6tud]

Hallo Vladimir,

I don't understand everything yet, but based on interFoam it should be sth like that:

ScalarField rUA = 1.0/UEqn.A();
surfaceScalarField rUAf = fvc::interpolate(rUA);

U = rUA * UEqn.H();

surfaceScalarField phi = fvc::interpolate(U) & mesh.Sf();
surfaceScalarField Ff = fvc::interpolate(F) & mesh.Sf();

[...]

fvScalarMatrix pEqn
(
fvm::div(rUAf*fvm::snGrad(p)*mesh.magSf()/rho) == fvc::div(phi) + fvc::div(rUAf*Ff);
);

[...]

U += rUA*fvc::reconstruct(Ff/rUAf);

I am not sure how to transform div.(snGrad), but probably we have to use:

fvScalarMatrix pEqn
(
fvm::laplacian(rUAf,p*mesh.magSf()/rho) == fvc::div(phi) + fvc::div(rUAf*Ff) ;
);

In interFoam the velocity is corrected only for the force. But they put phi and Ff together (and name it phi) -> I am not sure why.
Working on pisoFoam, I think rho has to be removed. I do not know if phi has to contain only the velocity flux, or the 'force flux', too. Any idea?

For my understanding: Can someone give me a hint, what 'reconstruct' does? It 'adds' terms to the equation and changes A and H? But there must be a difference between

solve ( UEqn == fvc::reconstruct ( - fvc::snGrad(p)*mesh.magSf()/rhoa + fvc::interpolate(F) & mesh.Sf() ) )

and

solve ( UEqn == - fvc::snGrad(p)*mesh.magSf()/rhoa + fvc::interpolate(F) & mesh.Sf() )

[vahid.najafi]

Hi Dear alberto again:

I want to add surface tension(sigma) in one solver,for this reason I added :
#include ''fvCFD.H''
fvc::interpolate(interface.sigma())

in this code:
Foam::tmp<Foam::volScalarField>
Foam:haseChangeTwoPhaseMixtures::SchnerrSauer: Coeff
(
const volScalarField& p
) const
{
volScalarField limitedAlpha1(min(max(alpha1_, scalar(0)), scalar(1)));
volScalarField rho
(
limitedAlpha1*rho1() + (scalar(1) - limitedAlpha1)*rho2()
);
return

//......I want to change it( <<sigma>> surface tension multiple in it):
(3*rho1()*rho2())*sqrt(2/(3*rho1()))*(fvc::interpolate(interface.sigma()))
*rRb(limitedAlpha1)/(rho*sqrt(mag(p - pSat()) + 0.01*pSat()));
//.................................................. ......
}
dont successful wmake, and seen(was not declared ):
phaseChangeTwoPhaseMixtures/SchnerrSauer/SchnerrSauer.C:113: error: 'interface' was not declared in this scope
make: *** [Make/linux64GccDPOpt/SchnerrSauer.o] Error 1

please help me,and tell me ,How to correct this problem???

[be_inspired]

HI all,

How could it be done in case it is not a body force but a vectorField?

All is in relation to the rotorDiskSource and how introducing this source in the momentum equation can create wiggles in the pressure field and velocity field.

Thank you

[chinaduck]

Originally Posted by harry

Does Openfoam implement body-force-weighted scheme for pressure interpolation?

Hello Harry,
Did you finally solve your problem? I have met a similar problem regarding the body-force-weighted scheme for pressure interpolation.
Do you have any ideas on this? Thanks a lot for your help and time!

Best,

Peter

[Novel]

Hi,
I know this treat is old but I'm struggling with the implementation of the force term by the use of the Rhi-Chow interpolation. I'm using pimpleFoam.

So far I added in the momentum predictor my force term like this:

surfaceScalarField B_fs = fvc::interpolate(lorentz/rho) & mesh.Sf()

if (pimple.momentumPredictor())
{
   solve
   (
      UEqn   
      ==
      fvc::reconstruct
      (   
         - fvc::snGrad(p)* mesh.magSf()
         + B_fs    
      )
    );
    fvOptions.correct(U);
}

In pEqn.H

I added a flux term in the phiHbyA equation:

surfaceScalarField rAUf("rAUf", fvc::interpolate(rAU));
surfaceScalarField phil(rAUf*B_fs);

surfaceScalarField phiHbyA
(
    "phiHbyA",
    fvc::flux(HbyA)
  + fvc::interpolate(rAU)*fvc::ddtCorr(U, phi)
  + phil
);

and at the end I correct the momentum source with the pressure gradient flux:

U = HbyA + rAU*fvc::reconstruct((phil - pEqn.flux())/rAUf);

But when I run the solver I get an continuity error:

Continuity error cannot be removed by adjusting the outflow.
Please check the velocity boundary conditions and/or run potentialFoam to initialise the outflow.
Total flux              : 8.226416123979263e-05
Specified mass inflow   : 0.0005613679754126387
Specified mass outflow  : 6.154738583391114e-17
Adjustable mass outflow : 0

What do I miss here?

Thanks in advance

[Novel]

The problem is solved. I had to add the flux from the force after adjusting phi.

surfaceScalarField phil(rAUf*B_fs);
surfaceScalarField phiHbyA
(
    "phiHbyA",
    fvc::flux(HbyA)
  + fvc::interpolate(rAU)*fvc::ddtCorr(U, phi)
);

MRF.makeRelative(phiHbyA);

adjustPhi(phiHbyA, U, p);

phiHbyA += phil;

[gaza]

Originally Posted by Novel

The problem is solved. I had to add the flux from the force after adjusting phi.

surfaceScalarField phil(rAUf*B_fs);
surfaceScalarField phiHbyA
(
    "phiHbyA",
    fvc::flux(HbyA)
  + fvc::interpolate(rAU)*fvc::ddtCorr(U, phi)
);

MRF.makeRelative(phiHbyA);

adjustPhi(phiHbyA, U, p);

phiHbyA += phil;

Hi Novel,
Have you also modified the line

            // Calculate the conservative fluxes
            phi = phiHbyA - p_rghEqn.flux();

??

[gaza]

Originally Posted by alberto

I will try to give you the basic idea in the case of a generic body force term F, so that your momentum equation reads:

(1) ddt(U) + div(UU) = div(tau) - grad(p)/rho + F

This equation can be written in semi-discrete form as:

(2) A*U = H - grad(p)/rho + F

where the pressure gradient and the force term we want to include in the momentum interpolation is left explicitly out at this point. In other words, only the first three terms (time derivative, unsteady, divergence of the stress tensor) of Eq. 1 are used to define A and H. For example, in an incompressible code, this could read (it might be different, depending on the solver):

     fvm::ddt(U)
  + fvm::div(phi, U)
  + turbulence->divDevReff(U)

Normally the pressure gradient would be treated directly as a source term in an incompressible code (see pisoFoam for example). Here we assume we want to treat it with the "improved approach".

From Eq. 2, interpolating on faces and dotting with the surface area vector S the pressure gradient and the force term, we have:

- snGrad(p)*|S|/rho + F_f . S

where F_f is F interpolated on faces. This in OF corresponds to:

- fvc::snGrad(p)*mesh.magSf()/rhoa + fvc::interpolate(F) & mesh.Sf()

This term can be used as argument of fvc::reconstruct() to add the contribution of the pressure gradient and of the force term to the momentum predictor:

solve
(
      UEqn == fvc::reconstruct
      (
          - fvc::snGrad(p)*mesh.magSf()/rhoa 
          + fvc::interpolate(F) & mesh.Sf()
   )
)

Then OpenFOAM re-computes U as

U = H/A

where A is updated with the predicted value of U. However, keep in mind that H does not directly include the effect of grad(p) and F.

From Eq. 2 we can derive the flux to construct the pressure equation:

phi = (H/A)_f . S - (1/A)_f * snGrad(p) |S|/rho + (1/A)_f*F_f . S

and imposing div(phi) = 0, you obtain the pressure equation, which will include for the first time all the effects.

The flux is then corrected based on the solution of the pressure equation, and the velocity correction is reconstructed from the flux. Only at this point U will "see" the effect of p and F.

If you take a look at VOF solvers (i.e. interFoam), you will notice that the code does not solve for p, but for p_rgh = p + rho*gh. This is another way to treat the gravity to address the weakness of the Rhie-Chow interpolation, but the idea is similar.

I hope this helps, but please let us know if you have more questions

Best,

Hi Alberto,
I am trying to add body force to momentum equation like this


where M is a body force. I based on buoyantBoussinesqPimpleFoam solver. I followed according to the instructions you provided and others on the forum and I changed the UEqn.H and pEqn.H as follows
UEqn.H

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

    UEqn.relax();

    fvOptions.constrain(UEqn);

    if (pimple.momentumPredictor())
    {
        solve
        (
            UEqn
         ==
            fvc::reconstruct
            (
                (
                  - ghf*fvc::snGrad(rhok)
                  - fvc::snGrad(p_rgh)
                )*mesh.magSf()
              - (fvc::interpolate(M) & mesh.Sf())
            )
        );

        fvOptions.correct(U);
    }

pEqn.H

{
    volScalarField rAU("rAU", 1.0/UEqn.A());
    surfaceScalarField rAUf("rAUf", fvc::interpolate(rAU));
    volVectorField HbyA(constrainHbyA(rAU*UEqn.H(), U, p_rgh));

    surfaceScalarField phig(-rAUf*ghf*fvc::snGrad(rhok)*mesh.magSf());
    surfaceScalarField phiM(-rAUf*(fvc::interpolate(M) & mesh.Sf()));

    surfaceScalarField phiHbyA
    (
        "phiHbyA",
        fvc::flux(HbyA)
      + rAUf*fvc::ddtCorr(U, phi)
      + phig
      + phiM
    );

    MRF.makeRelative(phiHbyA);

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

                                                

    while (pimple.correctNonOrthogonal())
    {
        fvScalarMatrix p_rghEqn
        (
            fvm::laplacian(rAUf, p_rgh) == fvc::div(phiHbyA)
        );

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

        p_rghEqn.solve(mesh.solver(p_rgh.select(pimple.finalInnerIter())));

        if (pimple.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 + rAU*fvc::reconstruct((phiM + phig - p_rghEqn.flux())/rAUf); //5   (+)
            U.correctBoundaryConditions();
            fvOptions.correct(U);
        }
    }

    #include "continuityErrs.H"

    p = p_rgh + rhok*gh;

    if (p_rgh.needReference())
    {
        p += dimensionedScalar
        (
            "p",
            p.dimensions(),
            pRefValue - getRefCellValue(p, pRefCell)
        );
        p_rgh = p - rhok*gh;
    }
}

However I got high residuals for p_rgh as in the image attached. Have I missed something or done something wrong?

[Novel]

Hi gaza,
I had problems when adding a source term without adjusting the flux to ensure mass conservation. Try to add the source terms after adjusting phi to your flux phiHybA.

Originally Posted by gaza

hi alberto,
peqn.h

{
    volscalarfield rau("rau", 1.0/ueqn.a());
    surfacescalarfield rauf("rauf", fvc::interpolate(rau));
    volvectorfield hbya(constrainhbya(rau*ueqn.h(), u, p_rgh));

    surfacescalarfield phig(-rauf*ghf*fvc::sngrad(rhok)*mesh.magsf());
    surfacescalarfield phim(-rauf*(fvc::interpolate(m) & mesh.sf()));

    surfacescalarfield phihbya
    (
        "phihbya",
        fvc::flux(hbya)
      + rauf*fvc::ddtcorr(u, phi)
    );
    mrf.makerelative(phihbya);
    adjustPhi(phihbya, U, p_rgh);
    phihyba += phig + phim;


    // update the pressure bcs to ensure flux consistency
    constrainpressure(p_rgh, u, phihbya, rauf, mrf);

                                                

    while (pimple.correctnonorthogonal())
    {
        fvscalarmatrix p_rgheqn
        (
            fvm::laplacian(rauf, p_rgh) == fvc::div(phihbya)
        );

        p_rgheqn.setreference(prefcell, getrefcellvalue(p_rgh, prefcell));

        p_rgheqn.solve(mesh.solver(p_rgh.select(pimple.finalinneriter())));

        if (pimple.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 + rau*fvc::reconstruct((phim + phig - p_rgheqn.flux())/rauf); //5   (+)
            u.correctboundaryconditions();
            fvoptions.correct(u);
        }
    }

    #include "continuityerrs.h"

    p = p_rgh + rhok*gh;

    if (p_rgh.needreference())
    {
        p += dimensionedscalar
        (
            "p",
            p.dimensions(),
            prefvalue - getrefcellvalue(p, prefcell)
        );
        p_rgh = p - rhok*gh;
    }
}

Please let me know if this solved your problem.

[gaza]

Hi Novel,
Thank you for your tip however in my case I do not see the difference after adjustPhi.