CFD Online Discussion Forums

CFD Online Discussion Forums (https://www.cfd-online.com/Forums/)
-   OpenFOAM Running, Solving & CFD (https://www.cfd-online.com/Forums/openfoam-solving/)
-   -   interpolation for pressure (https://www.cfd-online.com/Forums/openfoam-solving/90435-interpolation-pressure.html)

harry July 10, 2011 08:33
interpolation for pressure
 
Does Openfoam implement body-force-weighted scheme for pressure interpolation?

alberto July 11, 2011 02:36
What solver are you referring to?

harry July 11, 2011 02:44
Quote:
Originally Posted by alberto (Post 315587)
What solver are you referring to?
I am new to OpenFoam. this type of interpolations is necessary in the Eulerian model based on the SIMPLE velocity and pressure coupling method in the solution of gas-solid/liquid flows where large source terms present due to the strong interaction between phases. it is necessary associated with the intrinsic flaw in Rhie and Chow momentum interpolation practice. A big CFD issue but very limited publications devoted to the problem. It has puzzled me for a long time.

alberto July 11, 2011 03:09
Yes, bubbleFoam and twoPhaseEulerFoam implement an improved interpolation approach to treat gravity and drag. In OpenFOAM such a practise is implemented in some of the solvers (buoyant, multiphase, ...). You can see how this is managed in bubbleFoam here: http://openfoamwiki.net/index.php/BubbleFoam


Basically the force terms are treated as a source term in the momentum equation (well in bubbleFoam they are literally dropped, see interFoam for a complete example), after face-flux reconstruction, and then are moved to the flux used to construct the pressure equation.

You might also want to take a look at how OpenFOAM does Rhie-Chow:

http://web.student.chalmers.se/~f98faka/rhiechow.pdf

P.S. My area of research is multiphase flow too, and I absolutely agree there is not enough good literature on the topic, and what is available is usually hidden in some conference proceedings!

harry July 11, 2011 06:57
Quote:
Originally Posted by alberto (Post 315595)
Yes, bubbleFoam and twoPhaseEulerFoam implement an improved interpolation approach to treat gravity and drag. In OpenFOAM such a practise is implemented in some of the solvers (buoyant, multiphase, ...). You can see how this is managed in bubbleFoam here: http://openfoamwiki.net/index.php/BubbleFoam


Basically the force terms are treated as a source term in the momentum equation (well in bubbleFoam they are literally dropped, see interFoam for a complete example), after face-flux reconstruction, and then are moved to the flux used to construct the pressure equation.

You might also want to take a look at how OpenFOAM does Rhie-Chow:

http://web.student.chalmers.se/~f98faka/rhiechow.pdf

P.S. My area of research is multiphase flow too, and I absolutely agree there is not enough good literature on the topic, and what is available is usually hidden in some conference proceedings!
thank you so much, alberto! this is a great help to me. Would you please check the link to download the pdf regarding Rhie-chow interpolation. it is not working. If possible, you can also send it to the following email address s.b.kuang@gmail.com.

l_r_mcglashan July 11, 2011 07:01
Try this:

http://www.tfd.chalmers.se/~hani/kur...7/rhiechow.pdf

alberto July 11, 2011 10:41
Sorry, the link on the wiki (where I took it) to the Rhie-Chow document is broken. Laurence gave you the correct one.

harry July 12, 2011 00:14
I have read through the relevant documents but cannot get the main idea about how OpenFoam implement the momentum interpolation under the condition with/without body force as source terms. I appreciate that some guy would like to outline this method or redirect me to some further publications.
Best regards,
Shibo

alberto July 12, 2011 01:34
Quote:
Originally Posted by harry (Post 315783)
I have read through the relevant documents but cannot get the main idea about how OpenFoam implement the momentum interpolation under the condition with/without body force as source terms. I appreciate that some guy would like to outline this method or redirect me to some further publications.
Best regards,
Shibo
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):

Code:
    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:

Code:
- 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:

Code:
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 :D

Best,

harry July 12, 2011 07:12
thank you so much, alberto! this is a quite interesting treatment which I am not quite sure that I have completely understood it.
Below is my understanding of the entire procedure, which seems to be achieved through two steps.
1) OpenFoam solves the momentum twice.
Firstly, the momentum is solved with pressure gradient and body force.
Secondly, the momentum is solved without pressure gradient and body force.
2) the velocity from the second solution of momentum is used to establish pressure correction equation, considering the contributions of pressure gradient and body force.

If my understanding is correct, it seems to me that it is not necessary to interpolate of gradient and body force to the cell face and add them to the momentum.

Did I missed something? I need to get confirmed before I can raise other questions.

BTW: I am completely new to OpenFoam, and cannot thus far discuss this using Openfoam style language.

Andrea_85 July 12, 2011 12:05
Hi all,
i'm writing my problem here, since it is similar.
i have a question concerning the solution procedure in OF. I'm struggling with problem of spurious velocities at the interface using interFoam and i've read that basically these spurious currents could come from a sort of errors in the splitting of the equations, because there is a force unbalance between the pressure gradient and the surface force term in the 2 steps predictor and corrector. They have suggested to add another correction in the flux calculation (before solving the pressure), which basically has to take into account the difference between these two terms.

here the post of my problem http://www.cfd-online.com/Forums/ope...tml#post315660
(in the post can be also found the title of the paper)

i would like to know what do you think about that and if it would be possibile to implement in OF?

thanks
andrea

alberto July 12, 2011 14:50
Quote:
Originally Posted by harry (Post 315814)
thank you so much, alberto! this is a quite interesting treatment which I am not quite sure that I have completely understood it.
Below is my understanding of the entire procedure, which seems to be achieved through two steps.
1) OpenFoam solves the momentum twice.
Firstly, the momentum is solved with pressure gradient and body force.
Secondly, the momentum is solved without pressure gradient and body force.
The full PDE for the momentum equation is solved once, with all the terms (well, bubbleFoam is an exception, since the momentum equation does not contain some of the terms, which are directly moved to the pressure equation).
The estimated value of U from this equation is used to update what OF calls H (see my previous post). Then, at this point the flux is defined consistently with the semi-discrete form of the momentum equation, in order to build the pressure equation. This means you need to compute U = H/A, and add the terms formally left out in the predictor.

Quote:
2) the velocity from the second solution of momentum is used to establish pressure correction equation, considering the contributions of pressure gradient and body force.

If my understanding is correct, it seems to me that it is not necessary to interpolate of gradient and body force to the cell face and add them to the momentum.
It depends on how accurately you want to estimate H. If you do not have particularly large terms, I agree. In some cases, like the case of Euler-Euler codes for granular materials for example, where the granular phase pressure is very large, then including it makes quite some difference. Similarly for strongly changing body forces.

Best,

harry July 12, 2011 21:14
Thanks, alberto! this helps a lot! I have almost got the picture of the method in OF, but still with one unclear point. I appreciate your further helps. My questions are given below:



Quote:
Originally Posted by alberto (Post 315788)

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:

Code:
- 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:

Code:
solve
(
      UEqn == fvc::reconstruct
      (
          - fvc::snGrad(p)*mesh.magSf()/rhoa
          + fvc::interpolate(F) & mesh.Sf()
  )
)
I cannot understand that why OF needs to interpolate pressure gradient and body force to cell face. How can I understand the improved method? what is its difference from the ordinary method where pressure gradient and body force are treated as source terms?

alberto July 12, 2011 21:56
Quote:
Originally Posted by harry (Post 315900)
I cannot understand that why OF needs to interpolate pressure gradient and body force to cell face. How can I understand the improved method? what is its difference from the ordinary method where pressure gradient and body force are treated as source terms?
Hi, this is not a requirement. If you look at many solvers (incompressible codes for example), the momentum equation is solved directly, treating the pressure term as source term:

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

This does not work so well in the case of strong density gradients or changes in porosity, for example. It is a known problem, and one of the solutions was proposed in Zhang S., Zhao, X., General formulations for Rhie-Chow interpolation, ASME Heat Transfer/Fluids Engineering Summer Conference, HT-FED04, Charlotte, USA, 2004. Such a solution suggests to average the body force term (and other strongly varying terms) on cell faces, and then reconstruct the source term from the face-averaged value. This is very similar to what OpenFOAM does in the lines of codes you are referring to.

Best,

harry July 13, 2011 01:10
Dear alberto,
In the Zhang's paper, the treatment of body force seems to be similar to that done by OF. It is achieved through threes steps:
1. First define body forces primitively at cell centers, FP
2. Average these to cell faces, Ff, using weight factors of one half.
3. Redefine F at cell centers in a manner consistent with the second order formula for nodal pressure gradient.

I cannot understand the third step. this is also the reason that I cannot understand the treatment in OF
Would you please detail the theory for the third step? In OF, this is done by fvc: reconstruct.

Andrea_85 July 14, 2011 04:30
Hi,
if the momentum predictor is turned off, where the solution "can see" the effect of the old pressure?
(I'm talking about interFoam) in the flux definition only the surface force and the gravity are explicitly added to U_f=(H/A)_f (where both H and A do not depend on the old pressure...correct me if i'm wrong)
So if the momentum predictor is turned off, does the new pressure depend on the old pressure?

thanks
andrea

alberto July 14, 2011 11:09
Quote:
Originally Posted by harry (Post 315913)
Dear alberto,
In the Zhang's paper, the treatment of body force seems to be similar to that done by OF. It is achieved through threes steps:
1. First define body forces primitively at cell centers, FP
2. Average these to cell faces, Ff, using weight factors of one half.
3. Redefine F at cell centers in a manner consistent with the second order formula for nodal pressure gradient.

I cannot understand the third step. this is also the reason that I cannot understand the treatment in OF
Would you please detail the theory for the third step? In OF, this is done by fvc: reconstruct.
Hi,

@harry: each face flux will provide you a contribution to the cell centered value, which can be re-constructed then based on the face values you defined. In other words you integrate the contribution of the face values to recover the cell value by means of integration.

@Andrea_85: the velocity U sees the pressure immediately after the solution of the pressure equation, when you correct it based on the new pressure gradient (you do this by flux reconstruction in interFoam).

Best,

Andrea_85 July 14, 2011 12:16
Hi Alberto and thanks for reply.
What i cannot figure out is the case in which the momentum predictor is turned off (set off in the PISO subdictionary in fvSolution). Which is in this case the predicted velocity? In Ueqn.H OF defines A and H, but H still depends on the velocity (that is unknown in my understanding if the momentum preditor is not solved). The system is:
A[U]=[SourceTerm],

then you can split in diagonal and off-diagonal part

A[U] + Aoffd[U] = [S]

and then H is defined as

H=[S]-aoffd[U]

If the momentum predictor is turned to off, H is something like:

H= - sum_neigbours(a_n*U_n)+U^0/deltaT
(U^0 is the velocity from the last time step)

the velocities of the "neighboring" are unknown (or not?...are from the previous time step?)
so i dont understand how can you calculate the velocity using something that depends on the velocity itself (U=H(U)/A)...


Best Regards

andrea

makaveli_lcf September 6, 2011 10:09
Hi all!

Thanx for a very interesting discussion! Would you please clarify what does fvc::reconstruct operator actually do? Is it:

fvc::reconstruct (f_flux) = Surf_integral (f_face_interpolated . Sf) = Sum_over_all_faces(f_face_interpolated . Sf) ?????

Thank you in advance!

alberto September 6, 2011 21:56
From ~/OpenFOAM/OpenFOAM-2.0.x/src/finiteVolume/finiteVolume/fvc/fvcReconstruct.C, you find that the reconstruct method is defined as:

Code:
template<class Type>
tmp
<
    GeometricField
    <
        typename outerProduct<vector,Type>::type, fvPatchField, volMesh
    >
>
reconstruct
(
    const GeometricField<Type, fvsPatchField, surfaceMesh>& ssf
)
{
    typedef typename outerProduct<vector, Type>::type GradType;

    const fvMesh& mesh = ssf.mesh();

    tmp<GeometricField<GradType, fvPatchField, volMesh> > treconField
    (
        new GeometricField<GradType, fvPatchField, volMesh>
        (
            IOobject
            (
                "volIntegrate("+ssf.name()+')',
                ssf.instance(),
                mesh,
                IOobject::NO_READ,
                IOobject::NO_WRITE
            ),
            inv(surfaceSum(sqr(mesh.Sf())/mesh.magSf()))
          & surfaceSum((mesh.Sf()/mesh.magSf())*ssf),
            zeroGradientFvPatchField<GradType>::typeName
        )
    );

    treconField().correctBoundaryConditions();

    return treconField;
}
Best,

makaveli_lcf September 7, 2011 08:27
1 Attachment(s)
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?

alberto September 7, 2011 21:16
I do not have a reference for this, but I think you got it.

akidess September 8, 2011 03:37
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 September 8, 2011 03:39
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 September 13, 2011 01:20
Alexander,

Quote:
Originally Posted by makaveli_lcf (Post 323277)
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 July 26, 2012 04:28
Could you be so kind to write down the corresponding code for the equation for p?
 
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


Quote:
Originally Posted by alberto (Post 315788)
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):

Code:
    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:

Code:
- 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:

Code:
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 :D

Best,

dl6tud July 29, 2012 15:20
Hallo Vladimir,

I don't understand everything yet, but based on interFoam it should be sth like that:
Quote:
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:

Quote:
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

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

vahid.najafi July 30, 2012 01:33
plz help me
 
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 June 4, 2014 10:51
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 February 23, 2015 09:16
Quote:
Originally Posted by harry (Post 315532)
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 May 12, 2017 09:26
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:

Code:
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:

Code:
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:

Code:
U = HbyA + rAU*fvc::reconstruct((phil - pEqn.flux())/rAUf);
But when I run the solver I get an continuity error:
Code:
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 May 18, 2017 05:15
The problem is solved. I had to add the flux from the force after adjusting phi.

Code:
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 September 19, 2017 04:42
Quote:
Originally Posted by Novel (Post 649373)
The problem is solved. I had to add the flux from the force after adjusting phi.

Code:
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
Code:
            // Calculate the conservative fluxes
            phi = phiHbyA - p_rghEqn.flux();
??

gaza September 21, 2017 04:29
high residuals for p_rgh
 
1 Attachment(s)
Quote:
Originally Posted by alberto (Post 315788)
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):

Code:
    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:

Code:
- 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:

Code:
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 :D

Best,
Hi Alberto,
I am trying to add body force to momentum equation like this
\frac{\partial \mathbf{U}}{\partial t}+\frac{\nabla\cdot(\mathbf{UU})}{\rho}-\frac{\nabla\cdot\mathbf{\tau}}{\rho} =-\frac{\nabla p}{\rho}+\rho_k \mathbf{g} - \mathbf{M}

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
Code:
    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
Code:
{
    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 November 3, 2017 04:34
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.
Quote:
Originally Posted by gaza (Post 665120)
hi alberto,
peqn.h
Code:
{
    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 November 10, 2017 06:35
Hi Novel,
Thank you for your tip however in my case I do not see the difference after adjustPhi.


All times are GMT -4. The time now is 08:48.