CFD Online Discussion Forums (https://www.cfd-online.com/Forums/)
-   OpenFOAM Programming & Development (https://www.cfd-online.com/Forums/openfoam-programming-development/)
-   -   Continuity equation in implicit form (https://www.cfd-online.com/Forums/openfoam-programming-development/109819-continuity-equation-implicit-form.html)

 fredo490 November 27, 2012 12:43

Continuity equation in implicit form

Dear everyone,
I'm looking for a way of converting my continuity equation in an implicit form. This equation correspond to a particle tracking using an Eulerian approach. "rho" is not a gas, it is my particle concentration.

Right now I have :
Continuity equation :

Code:

Info<< "Solve Continuity" << endl;
tmp<fvScalarMatrix> CEqn
(
fvm::ddt(rho)
+ fvc::div(phi)
);
CEqn().solve();

I want to convert "fvc::div(phi)" in an implicit form.

I'm a little bit confused by the OpenFoam documentation : http://www.foamcfd.org/Nabla/guides/...sGuidese9.html

It says that the convection term can be Imp/Exp with the following form: div(psi,scheme)* where psi is a surfaceScalarField (just like my "phi").

However, if I try the following code, I get a mistake while compiling:
Code:

word scheme("div(phi)");

Info<< "Solve Continuity" << endl;
tmp<fvScalarMatrix> CEqn
(
fvm::ddt(rho)
+ fvm::div(phi, scheme)
);
CEqn().solve();

Does anybody has an idea ?

 meindert November 27, 2012 13:21

Did you already try fvm::div(phi,rho)? I think that should do the trick for you.

 fredo490 November 27, 2012 20:54

Well, if I use fvm::div(phi,rho) , it would imply that phi is an interpolation of U only. However, in my solver, phi is an interpolation of U and rho.

I'm not sure, can I use 2 different phi ?
- phiU = (interpolation of U) for the continuity equation
- phiRhoU = (phiU * interpolation of Rho) for the momentum equation

For more details, my model is defined as followed:
Nomenclature :
- volume fraction of the particles
- velocity of the particles
- function of the relative Reynolds number
- time respond of the particle
- velocity of the fluid (constant)

Continuity equation :

Momentum equation

 meindert November 28, 2012 06:02

Hi Frédéric,

Defining two different phi's would indeed be the way to go. Otherwise, you would end up with two unknowns in the continuity equation ( and ). You would not be able to solve for both.
I would be careful with the definition of phiRhoU. It might be better to do an interpolation of U * Rho, but this is something you could test of course.

How do you solve the momentum equation? Is given from solving the continuity equation?

I have one small remark. You seem to be using both rho and alpha for the volume fraction of the particles. Correct me if I am wrong. I think it would be better to be consistent to avoid confusion.

Good luck.

 fredo490 November 28, 2012 06:48

I have applied your method and it seems to work perfectly. I made a short run and the results look similar to my explicit solver.

To gives some details, here is roughly what happen in my code
Code:

//------------------------------------------------------------//
phi = fvc::interpolate(U) & mesh.Sf();

//------------------------------------------------------------//
Info<< "Solve Continuity" << endl;
tmp<fvScalarMatrix> CEqn
(
fvm::ddt(rho)
+ fvm::div(phi,rho)
);
CEqn().solve();

//------------------------------------------------------------//
phiR = fvc::interpolate(rho)*(fvc::interpolate(U) & mesh.Sf());

//------------------------------------------------------------//
Info<< "Solve Momentum" << endl;
tmp<fvVectorMatrix> UEqn
(
fvm::ddt(rho, U)
+ fvm::div(phiR, U)
==
fRer/tau*rho*Uf
- fvm::Sp(fRer/tau*rho, U)
);
UEqn().solve();

I can now use much bigger timestep and the computation remains stable.
Thanks again ;)

 meindert November 28, 2012 06:53

What happens if you try

Code:

phiR = fvc::interpolate(rho*U) & mesh.Sf()
?

 fredo490 November 28, 2012 07:18

It also works and it simplify my code :D

I have another question related to the field management. Right now I have a RelativeReynolds number stored in a volScalarField and a drag value in another volScalarField.

Until now, I use the following code :
Code:

// Compute the Relative Reynolds number field
Rer = mag( Uf - U )

// Compute the drag value field
fRer = 1.0 + 0.15 * pow(Rer, 0.687) + ...

But my function of drag value is limited to a small range of Rer... Therefore I want to use a code like this:
Quote:
 if (Rer <= A) { fRer = a; } else if ((Rer > A) && (Rer <= B)) { fRer = b; } else if ((Rer > B) && (Rer <= C)) { fRer = c; } else { fRer = d; }
However, this code cannot be applied to a Field. Is there any method to avoid the cell loop ?

 All times are GMT -4. The time now is 05:02.