[gaza]

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]

Originally Posted by sharonyue

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:

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]

Originally Posted by sharonyue

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:

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]

Originally Posted by kc95

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]

Originally Posted by fportela

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]

Originally Posted by sharonyue

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:

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]

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:

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

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

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]

add this

div(HbyA)      Gauss upwind;

ps

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

[anhkenyt]

Originally Posted by gaza

add this

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,

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:

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]

Hi
I want to write below equation to foam-extend 4.0

I thought to use

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]

Originally Posted by Hgholami

Hi
I want to write below equation to foam-extend 4.0

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:

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]

Originally Posted by anhkenyt

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

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:

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]

Originally Posted by anhkenyt

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.

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.

[saiguruprasad]

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]

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]

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.