CFD Online Discussion Forums

CFD Online Discussion Forums (
-   OpenFOAM Programming & Development (
-   -   Continuity equation in implicit form (

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 :
\frac{\partial \alpha}{\partial t} + {\nabla \alpha \vec{u}} = 0


        Info<< "Solve Continuity" << endl;
        tmp<fvScalarMatrix> CEqn
          + fvc::div(phi)

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

I'm a little bit confused by the OpenFoam documentation :

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:

word scheme("div(phi)");

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

Does anybody has an idea ?
Thanks in advance.

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

Actually, your advice rise more questions than answers ^^

For more details, my model is defined as followed:
Nomenclature :
- \alpha volume fraction of the particles
- \vec{u} velocity of the particles
- f(Re) function of the relative Reynolds number
- \tau time respond of the particle
- \vec{V} velocity of the fluid (constant)

Continuity equation :
\frac{\partial \alpha}{\partial t} + {\nabla \alpha \vec{u}} = 0

Momentum equation
\frac{\partial \alpha \vec{u}}{\partial t} + {\nabla \alpha \vec{u} \vec{u}} = \frac{f(Re)}{\tau } \alpha (\vec{V} -\vec{u})

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 (\alpha and \alpha\vec{u}). 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 \alpha 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

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

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

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

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

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


    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 :

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

if (Rer <= A)
{ fRer = a; }
else if ((Rer > A) && (Rer <= B))
{ fRer = b; }
else if ((Rer > B) && (Rer <= C))
{ fRer = c; }
{ 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.