CFD Online Discussion Forums - Meaning of "fvc::div(phi)"

CFD Online Discussion Forums (https://www.cfd-online.com/Forums/)

- OpenFOAM Programming & Development (https://www.cfd-online.com/Forums/openfoam-programming-development/)

- - Meaning of "fvc::div(phi)" (https://www.cfd-online.com/Forums/openfoam-programming-development/97922-meaning-fvc-div-phi.html)

enoch

February 28, 2012 10:51

Meaning of "fvc::div(phi)"

Hi Foamers,

When I look at pEqn.C in the "twoPhaseEulerFoam" solver, I don't really understand the following codes.

.......
for(int nonOrth=0; nonOrth<=nNonOrthCorr; nonOrth++)
{
fvScalarMatrix pEqn
(
fvm::laplacian(Dp, p) == fvc::div(phi)
);
.....

Speaking of div(phi), divergence of any scalar is zero mathematically.

Divergence reduce the order rank of maxtrix.
For an example, let's say, velcotiy vector u=(ux, uy)

div(u)=ux/dx + uy/dy ---> scalar quantity

But, scalar phi=phi(x,y)
div(phi)=0

So what is the meaning of div(phi) and what's a mathematical formulation for the term in openfoam?

sharonyue

July 16, 2013 04:01

Hi Kim,

This problem happens to me in the same time!http://www.cfd-online.com/Forums/ope...-equation.html Did you find a solution?

fportela

July 16, 2013 06:01

Isn't phi simply the mass flux rho*U*A ? check this

HTML Code:

http://openfoamwiki.net/index.php/Uguide_table_of_fields

and this

HTML Code:

http://openfoamwiki.net/index.php/Main_FAQ#What_is_the_field_phi_that_the_solver_is_writing

I have not used this particular solver, so I'm not sure what's going on, but from continuity: div(phi) = 0 implies that ddt(rho) is zero, if this is not the case, then div(phi) is not zero.

Cheers,
Felipe

sharonyue

July 16, 2013 08:31

Hi Felipre,

Thanks for your prompt reply,

Quote:

Originally Posted by fportela (Post 440008)

but from continuity: div(phi) = 0

regarding this, isnt it should be div(u)=0?

Bernhard

July 16, 2013 08:35

No, because U is the cell centre value, and the divergence is obtained from the face values, i.e. phi.

sharonyue

July 16, 2013 08:38

Hi Bernhard,
Could you explain more?

Bernhard

July 16, 2013 08:39

Can you be more specific?

sharonyue

July 16, 2013 08:44

Quote:

Originally Posted by Bernhard (Post 440031)

No, because U is the cell centre value, and the divergence is obtained from the face values, i.e. phi.

I mean all the books say from continuity we have: $\nabla \cdot u=0$ without mentioning whether its cell centre value or the face value. But Im sure I confused about this. So?

Bernhard

July 16, 2013 08:49

The book are correct, but http://www.cfd-online.com/Forums/vbL...5ac33fb3-1.gif is valid without a mesh. If you integrate this equation over a control volume, this converts to a summation over the faces of the velocity times the area ( http://en.wikipedia.org/wiki/Divergence_theorem ). phi is nothing less then the velocity at the face (times the density and the face area).

fportela

July 16, 2013 08:51

Quote:

Originally Posted by sharonyue (Post 440034)

I mean all the books say from continuity we have: $\nabla \cdot u=0$ without mentioning weather its cell centre value or the face value. But Im sure I confused about this. So?

This is only true for incompressible flow.

For incompressible flow, you have $\frac{D\rho}{Dt} = \frac{\partial\rho}{\partial t} + \vec{u} \cdot \nabla (\rho) = 0$

Plug this into continuity

$\frac{\partial\rho}{\partial t} + \nabla (\rho \vec{u}) = 0$

And you get

$\nabla \cdot \vec{u} = 0$

sharonyue

July 16, 2013 09:24

Quote:

Originally Posted by Bernhard (Post 440037)

$\int\limits_V^{} {\nabla \cdot udV} = \sum\limits_f {{u_f} \cdot {S_f} = 0}$

I know this, but phi is a scalar, what is div(scalar).....

Felipe. Yeah, you are right~

Bernhard

July 16, 2013 09:27

Quote:

Originally Posted by sharonyue (Post 440049)

I know this, but phi is a scalar, what is div(scalar).....

Be careful here. phi is a surfaceScalarField, so there is always a direction vector defined by the face area vector.

sharonyue

July 16, 2013 09:56

$\int\limits_V^{} {\nabla \cdot udV} = \sum\limits_f {{u_f} \cdot {S_f} = 0}$

Regarding this, if div(phi)=0 Do you mean:
$\int\limits_V^{} {u_f \cdot S_fdV} = 0$

Quote:

phi is a surfaceScalarField, so there is always a direction vector defined by the face area vector

Is there any difference between "surfaceScalarField" and "volScaklaField" expect where they are stored?

Bernhard

July 16, 2013 10:01

You apply the divergence theorem.
$\int\limits_V^{} {\nabla \cdot udV} = \int\limits_{\partial V}u\cdot\hat{n}dS =\sum\limits_f {{u_f} \cdot {S_f} = 0}$
You can't integrate face values over a volume as you are posting.

The difference between a volScalarField and a surfaceScalarField, is that for a volScalarField, there is a value stored per control volume or cell. For a surfaceScalarField there is a value stored per face.

You can of course interpolate from the one to the other, but this is only accurate to some order.

sharonyue

July 16, 2013 10:23

Quote:

Originally Posted by Bernhard (Post 440063)

Bernhard, Thanks very very much for your consistent help, But I still cannot understand this "div(u)=0" turns into "div(phi)=0" in OpenFOAM, even mathematicly div(vector) works but div(scalar) not.
Does it mean fvc::div(phi)=0 equals to sum(phi)=0 in OpenFOAM? Why doesnt it use sum(phi)=0 instead.
At last,I think I need to clear my head. Thanks for you patience. Really thankful.:)

sharonyue

July 17, 2013 23:52

1 Attachment(s)

Im still confused!!!!!!

div(u)=0 I know its meaning.

1) If via continuity, we have div(phi)=0, then we make an intergration on this equation like this:

[LaTeX Error: Syntax error]

Is this weird?

Bernhard

July 18, 2013 01:59

Yes, this is extremely weird. phi only lives in the discretized domain on faces of cells. The is no way you can do a volume integral over such a variable.

Do you know how to derive the equation $\nabla\cdot\vec{u}$ ? You can do this using a mass-micro-balanse. This has not yet anything to do with a mesh, but if you understand this, it is easily translated.

sharonyue

July 18, 2013 18:46

Quote:

Originally Posted by Bernhard (Post 440449)

Yeah, Im very clear about this $\nabla\cdot\vec{u}$ , and I know after intergral it can be a sum. I dont know whether you know what I dont know.

Code:

fvScalarMatrix pEqn

                (

                    fvm::laplacian(rAU, p) == SUM(phiHbyA)//If there existed "SUM"

                );

Even I can understand this. But I cannot understand fvc::div(phiHbyA).
Anyway thanks bro.

BTW, the unit of fvc::div(u) and fvc::div(phi) is the same:

Code:

fvc::div(HbyA):dimensions      [0 0 -1 0 0 0 0];



internalField   nonuniform List<scalar> 9(0.157469 -0.0674757 -0.487085 -0.960691 0.789192 0.715649 30.1762 0.258981 -30.5822);



fvc::div(phiHbyA):dimensions      [0 0 -1 0 0 0 0];



internalField   nonuniform List<scalar> 9(-0.435688 -0.0982547 0.152804 -0.0448476 0.720155 -0.152663 28.2072 0.236062 -28.5848);

ARTem

July 22, 2013 04:01

Take into account that mathematical and OpenFOAM's languages are different.
Originally, mathematical divergence is vector operation that gives you source or sink at a point. When you use finite volume method to discretize differential equation, you get linear form of diff. equation and volumes to store discrete variables. That's why in order to get divergence you should firstly compute fluxes at volume faces.
In OF there are two ways to compute divergence:
1) take fvc::div( of volVectorField ). You can see in sources, it calls for fvc::surfaceIntegrate ( volVectorField & mesh.Sf() ).
2) make first step manually (i.e. phi = U&mesh.Sf()) and call for fvc::div( surfaceScalarField phi ). Again, OF will make fvc::surfaceIntegrate ( surfaceScalarField ).
Mathematically div(scalar phi) has no sense. Because you think of variables as continuous fields, but in OF variables are discrete fields.

You are partially right, there is no "SUM", but "surfaceIntegrate" called by OF to compute div.

And again, phi = U & mesh.Sf(). div(phi) = div(U) = div(U&mesh.Sf()). That's why you get same dimensions.

sharonyue

August 23, 2013 00:08

Quote:

Originally Posted by ARTem (Post 441149)

U make it very very clear, Thanks very much. So I just make it a example.:

1. volVectorField U;
fvc::div(U);

2. volVectorField U;
phi = U & mesh.Sf();
fvc::div(phi);

So these two functions are the same rite? fvc::div(u)=fvc::div(phi).

Looks like div() has been overloaded.

I make a testify:

Code:

fvc::div(HbyA)dimensions      [0 0 -1 0 0 0 0];



internalField   nonuniform List<scalar> 9(0.000569046 -8.589e-06 -0.000648869 -0.00128281 0.000183935 0.00124553 0.00416702 3.71583e-06 -0.00422898);



fvc::div(phiHbyA):dimensions      [0 0 -1 0 0 0 0];



internalField   nonuniform List<scalar> 9(0.000569046 -8.589e-06 -0.000648869 -0.00128281 0.000183935 0.00124553 0.00416702 3.71583e-06 -0.00422898);

Its total the same. Good job.

gaza	March 25, 2015 08:22

Hi all
I have only the question why divU = fvc::div(phi)
When we see in the prof. Jasak's PhD in Eqn. 3.15 we have that: divU = sum(S.Uf)/Vp
where Vp is a volume of the cell.
So why in the code we don't have to do division by the cell volume?

kc95	May 26, 2018 03:23

Quote:

Originally Posted by sharonyue (Post 447595)

Code:

fvc::div(HbyA)dimensions      [0 0 -1 0 0 0 0];



internalField   nonuniform List<scalar> 9(0.000569046 -8.589e-06 -0.000648869 -0.00128281 0.000183935 0.00124553 0.00416702 3.71583e-06 -0.00422898);



fvc::div(phiHbyA):dimensions      [0 0 -1 0 0 0 0];



internalField   nonuniform List<scalar> 9(0.000569046 -8.589e-06 -0.000648869 -0.00128281 0.000183935 0.00124553 0.00416702 3.71583e-06 -0.00422898);

Its total the same. Good job.

How did you find out these numerical values? i am very new to OF. so my question might appear as silly. But, please do reply. Thanks in advance!!

kc95	May 26, 2018 03:28

Quote:

Originally Posted by sharonyue (Post 447595)

Code:

fvc::div(HbyA)dimensions      [0 0 -1 0 0 0 0];



internalField   nonuniform List<scalar> 9(0.000569046 -8.589e-06 -0.000648869 -0.00128281 0.000183935 0.00124553 0.00416702 3.71583e-06 -0.00422898);



fvc::div(phiHbyA):dimensions      [0 0 -1 0 0 0 0];



internalField   nonuniform List<scalar> 9(0.000569046 -8.589e-06 -0.000648869 -0.00128281 0.000183935 0.00124553 0.00416702 3.71583e-06 -0.00422898);

Its total the same. Good job.

How did you get this numerical values? I am very new to OpenFOAM, so the question might appear as very silly. Please do reply. Thanks in Advance!!

fportela

May 26, 2018 09:23

Quote:

Originally Posted by kc95 (Post 693637)

How did you get this numerical values? I am very new to OpenFOAM, so the question might appear as very silly. Please do reply. Thanks in Advance!!

What exactly do you want to know?

kc95	May 26, 2018 09:30

Quote:

Originally Posted by fportela (Post 693658)

What exactly do you want to know?

How did you get these values of fvc::div(phiHbyA) as wellas fvc::(HbyA)?
I mean, how can we calculate and display those values?

randolph

July 1, 2018 12:44

Quote:

Originally Posted by sharonyue (Post 447595)

Code:

fvc::div(HbyA)dimensions      [0 0 -1 0 0 0 0];



internalField   nonuniform List<scalar> 9(0.000569046 -8.589e-06 -0.000648869 -0.00128281 0.000183935 0.00124553 0.00416702 3.71583e-06 -0.00422898);



fvc::div(phiHbyA):dimensions      [0 0 -1 0 0 0 0];



internalField   nonuniform List<scalar> 9(0.000569046 -8.589e-06 -0.000648869 -0.00128281 0.000183935 0.00124553 0.00416702 3.71583e-06 -0.00422898);

Its total the same. Good job.

1. From the continuity equation, the divergence-free condition is strictly applied to the face flux in the FVM. So it makes sense that cell value is not strictly divergence free.

2. However, I am still confused that why fvc::div(U) = fvc::div(phi) ? Isn't the phi = U*A?

anhkenyt

October 13, 2020 00:39

Display value of fvc::div(HbyA)

Quote:

Originally Posted by sharonyue View Post
U make it very very clear, Thanks very much. So I just make it a example.:

1. volVectorField U;
fvc::div(U);

2. volVectorField U;
phi = U & mesh.Sf();
fvc::div(phi);

So these two functions are the same rite? fvc::div(u)=fvc::div(phi).

Looks like div() has been overloaded.

I make a testify:

Code:
fvc::div(HbyA)dimensions [0 0 -1 0 0 0 0];

internalField nonuniform List<scalar> 9(0.000569046 -8.589e-06 -0.000648869 -0.00128281 0.000183935 0.00124553 0.00416702 3.71583e-06 -0.00422898);

fvc::div(phiHbyA):dimensions [0 0 -1 0 0 0 0];

internalField nonuniform List<scalar> 9(0.000569046 -8.589e-06 -0.000648869 -0.00128281 0.000183935 0.00124553 0.00416702 3.71583e-06 -0.00422898);
Its total the same. Good job.

I tried it in peqn.H in pisoFoam like that:

Code:

volScalarField rAU(1.0/UEqn.A());

volVectorField HbyA(constrainHbyA(rAU*UEqn.H(), U, p));

surfaceScalarField phiHbyA

(

    "phiHbyA",

    fvc::flux(HbyA)

  + MRF.zeroFilter(fvc::interpolate(rAU)*fvc::ddtCorr(U, phi))

);

Info<<"fvc::div(phiHbyA): "<<fvc::div(phiHbyA)<<endl;

Info<<"fvc::div(HbyA): "<<fvc::div(HbyA)<<endl;

But I faced to error

Code:

FOAM FATAL IO ERROR: 

keyword div(HbyA) is undefined in dictionary "/home/anh/OpenFOAM/anh-6/applications/solvers/src/cyclicSimpleFoam/system/fvSchemes.divSchemes"

And this is my divSchemes

Code:

divSchemes

{

    default         none;

        div(HybA)  Gauss upwind; //vda add    

//        div(U)           Gauss upwind; //vda add 

        div(phi,U)      bounded Gauss linearUpwind grad(U);

    div(phi,k)      bounded Gauss limitedLinear 1;

    div(phi,epsilon) bounded Gauss limitedLinear 1;

    div(phi,omega)  bounded Gauss limitedLinear 1;

    div(phi,v2)     bounded Gauss limitedLinear 1;

    div((nuEff*dev2(T(grad(U))))) Gauss linear;

    div(nonlinearStress) Gauss linear;

}

Anyone could help me display this value ?

gaza	October 13, 2020 08:35

add this

Code:

div(HbyA) Gauss upwind;

ps

you have the typo "HybA" instead of "HbyA"

anhkenyt

October 13, 2020 09:13

Quote:

Originally Posted by gaza (Post 785152)

add this

Code:

div(HbyA) Gauss upwind;

ps

you have the typo "HybA" instead of "HbyA"

Thanks for the quick answer, sorry about my stupid mistake. However when I fixed it,

Code:

divSchemes

{

    default         none;

        div(HbyA)  Gauss upwind; //vda add    

//        div(U)           Gauss upwind; //vda add 

        div(phi,U)      bounded Gauss linearUpwind grad(U);

    div(phi,k)      bounded Gauss limitedLinear 1;

    div(phi,epsilon) bounded Gauss limitedLinear 1;

    div(phi,omega)  bounded Gauss limitedLinear 1;

    div(phi,v2)     bounded Gauss limitedLinear 1;

    div((nuEff*dev2(T(grad(U))))) Gauss linear;

    div(nonlinearStress) Gauss linear;

}

I have another error:

Code:

fvc::div(HbyA): 



--> FOAM FATAL IO ERROR: 

attempt to read beyond EOF



file: /home/anh/OpenFOAM/anh-6/applications/solvers/src/vda_meanVelocityForce/testcode/cyclicSimpleFoam/system/fvSchemes.divSchemes.div(HbyA) at line 31.

Could you give me some suggestions?

Hgholami

October 13, 2020 11:10

Hi
I want to write below equation to foam-extend 4.0
https://uupload.ir/files/ybl_untitled.jpg
I thought to use

Quote:

fvScalarMatrix pEqn
(
fvm::grad(p) == fvc::div(U)/runTime().deltaT().value()
);
pEqn.solve();

but I know it's wrong. Do you have any tip for this?
Thanks.

anhkenyt

October 13, 2020 12:03

Quote:

Originally Posted by Hgholami (Post 785171)

Hi
I want to write below equation to foam-extend 4.0
https://uupload.ir/files/ybl_untitled.jpg
I thought to use but I know it's wrong. Do you have any tip for this?
Thanks.

Hi Hgholami,
In your equation, does not grad p, It should be lapacian(p). You can try:

Quote:

fvScalarMatrix pEqn
(
fvm::laplacian(p) == fvc::div(U)*rhoF/dT
);

However, I don't know to extract the meaning of index "p" and "n+1" in your eqn.

gaza	October 13, 2020 12:40

Quote:

Originally Posted by anhkenyt (Post 785160)

Thanks for the quick answer, sorry about my stupid mistake. However when I fixed it,

Code:

divSchemes

{

    default         none;

    div(HbyA)  Gauss upwind; //vda add    

//    div(U)       Gauss upwind; //vda add 

    div(phi,U)      bounded Gauss linearUpwind grad(U);

    div(phi,k)      bounded Gauss limitedLinear 1;

    div(phi,epsilon) bounded Gauss limitedLinear 1;

    div(phi,omega)  bounded Gauss limitedLinear 1;

    div(phi,v2)     bounded Gauss limitedLinear 1;

    div((nuEff*dev2(T(grad(U))))) Gauss linear;

    div(nonlinearStress) Gauss linear;

}

I have another error:

Code:

fvc::div(HbyA): 



--> FOAM FATAL IO ERROR: 

attempt to read beyond EOF



file: /home/anh/OpenFOAM/anh-6/applications/solvers/src/vda_meanVelocityForce/testcode/cyclicSimpleFoam/system/fvSchemes.divSchemes.div(HbyA) at line 31.

Could you give me some suggestions?

try Guass linear maybe
you wrote something in fvSchemes

Hgholami

October 13, 2020 22:22

Quote:

Originally Posted by anhkenyt (Post 785177)

Hi Hgholami,
In your equation, does not grad p, It should be lapacian(p). You can try:

However, I don't know to extract the meaning of index "p" and "n+1" in your eqn.

Hi anhkenyt
I want to use fractional-step projection method.
https://uupload.ir/files/rij7_untitled.jpg
the index "p" and n+1 is not matter. they are predicted velocity and next time step, respectively.
my mistake may be comparison delta_t with delta_p in Eq. 10.:rolleyes:

saiguruprasad

May 3, 2022 14:43

I am trying to write a new Particle Dispersion Model in OpenFOAM. The model requires the divergence of the Reynolds stress tensor. I am able to read the Reynolds stress tensor into the model and am able to perform the fvc::div operation.

My doubt is if this approach is correct. I tried fvc::div(U) in the same way and it gave really high values (O(1)) which are not correct because the simulation has converged and the fvc::div(U) should be a small value (O(1e-3)). So now I have doubts if fvc::div(R) is giving the right values.

I have the phi variable too which I think is U*surfaceArea. I computed div(phi, U) and div(phi), and both gave small values in O(1e-3).

Should I be using phi to compute the divergence of R? I am getting an error if I use fvc::div(phi,R). What is the difference between using div(R) and div(phi,R)?

kagen816

July 29, 2022 14:05

Hi Sai,
Have you found any clue, I found a similar problem. fvc::div(U) has a much greater value than fvc::div(phi) but they have the same dimension!
Thanks

saiguruprasad

July 30, 2022 08:40

Hi Kagen816.

It seems that div(phi) is correct if you want to verify the continuity equation. phi corresponds to the flux, while U is the cell centre value of velocity.

Divergence needs to be calculated at the faces.

But I am not sure how this applies to other variables. I am still looking at div(R), where R is the Reynolds Stress.

All times are GMT -4. The time now is 10:13.