[fredo490]

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 :

	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:

 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 ?
Thanks in advance.

[meindert]

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

[fredo490]

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

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]

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

What happens if you try

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

?

[fredo490]

It also works and it simplify my code


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; }
else
{ fRer = d; }

However, this code cannot be applied to a Field. Is there any method to avoid the cell loop ?