CFD Online Logo CFD Online URL
www.cfd-online.com
[Sponsors]
Home > Forums > OpenFOAM Programming & Development

Continuity equation in implicit form

Register Blogs Members List Search Today's Posts Mark Forums Read

Reply
 
LinkBack Thread Tools Display Modes
Old   November 27, 2012, 12:43
Default Continuity equation in implicit form
  #1
Senior Member
 
HECKMANN Frédéric
Join Date: Jul 2010
Posts: 236
Rep Power: 7
fredo490 is on a distinguished road
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

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 ?
Thanks in advance.
fredo490 is offline   Reply With Quote

Old   November 27, 2012, 13:21
Default
  #2
Member
 
Meindert de Groot
Join Date: Jun 2012
Location: Netherlands
Posts: 34
Rep Power: 5
meindert is on a distinguished road
Did you already try fvm::div(phi,rho)? I think that should do the trick for you.
meindert is offline   Reply With Quote

Old   November 27, 2012, 20:54
Default
  #3
Senior Member
 
HECKMANN Frédéric
Join Date: Jul 2010
Posts: 236
Rep Power: 7
fredo490 is on a distinguished road
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})
fredo490 is offline   Reply With Quote

Old   November 28, 2012, 06:02
Default
  #4
Member
 
Meindert de Groot
Join Date: Jun 2012
Location: Netherlands
Posts: 34
Rep Power: 5
meindert is on a distinguished road
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.

Last edited by meindert; November 28, 2012 at 06:23.
meindert is offline   Reply With Quote

Old   November 28, 2012, 06:48
Default
  #5
Senior Member
 
HECKMANN Frédéric
Join Date: Jul 2010
Posts: 236
Rep Power: 7
fredo490 is on a distinguished road
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
fredo490 is offline   Reply With Quote

Old   November 28, 2012, 06:53
Default
  #6
Member
 
Meindert de Groot
Join Date: Jun 2012
Location: Netherlands
Posts: 34
Rep Power: 5
meindert is on a distinguished road
What happens if you try

Code:
    phiR = fvc::interpolate(rho*U) & mesh.Sf()
?
meindert is offline   Reply With Quote

Old   November 28, 2012, 07:18
Default
  #7
Senior Member
 
HECKMANN Frédéric
Join Date: Jul 2010
Posts: 236
Rep Power: 7
fredo490 is on a distinguished road
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 :
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 ?
fredo490 is offline   Reply With Quote

Reply

Thread Tools
Display Modes

Posting Rules
You may not post new threads
You may not post replies
You may not post attachments
You may not edit your posts

BB code is On
Smilies are On
[IMG] code is On
HTML code is Off
Trackbacks are On
Pingbacks are On
Refbacks are On


Similar Threads
Thread Thread Starter Forum Replies Last Post
Conservative form of Navier Stokes equation. balkrishna OpenFOAM Running, Solving & CFD 2 January 25, 2012 09:33
Derivation of Momentum Equation in Integral Form Demonwolf Main CFD Forum 2 October 29, 2009 20:53
IcoFoam parallel woes msrinath80 OpenFOAM Running, Solving & CFD 9 July 22, 2007 02:58
Could anybody help me see this error and give help liugx212 OpenFOAM Running, Solving & CFD 3 January 4, 2006 19:07
continuity equation for the mixture mateus FLUENT 0 April 24, 2003 07:53


All times are GMT -4. The time now is 01:27.