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

Need Help w. Darcy Boussinesq Equ.

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

Like Tree2Likes

Reply
 
LinkBack Thread Tools Display Modes
Old   August 22, 2011, 10:13
Default Need Help w. Darcy Boussinesq Equ.
  #1
New Member
 
Werner
Join Date: Dec 2009
Location: Austria
Posts: 13
Rep Power: 7
nep101 is on a distinguished road
Hi everybody,

im new to Openfoam so i get stuck. I try solve this (dimensionless) Equ in a Simple Loop:

U-Ra*T*k=-grad(p)

were U is the velocity vectorfield, Ra is the Rayleighnumber, T is the temperature, k is a vector (Ra*T*k is a volScalarField) and p is the Pressure.

I want to creat a Matrix UEqn with
{
U-RaTk
}

But it dont work.

Code:
volScalarField V("V", U & mesh.C());
fvScalarMatrix UEqn
(
V-RaTk == fvc::grad(p) //Ratk is a VolScalarField
);



UEqn.solve();

Anybody got a hint?
Werner
nep101 is offline   Reply With Quote

Old   September 26, 2011, 00:56
Default
  #2
New Member
 
Tom
Join Date: Sep 2011
Posts: 11
Rep Power: 5
TomB is on a distinguished road
Hi Werner,

I am working with Darcy's eqns in OpenFoam as well.

how did you define the k vector ? Can you multiply it by a scalar the way you have done ?

Is it okay to have a UEqn with the double equals in it. I think not.

I Look fwd to working with you !





Quote:
Originally Posted by nep101 View Post
Hi everybody,

im new to Openfoam so i get stuck. I try solve this (dimensionless) Equ in a Simple Loop:

U-Ra*T*k=-grad(p)

were U is the velocity vectorfield, Ra is the Rayleighnumber, T is the temperature, k is a vector (Ra*T*k is a volScalarField) and p is the Pressure.

I want to creat a Matrix UEqn with
{
U-RaTk
}

But it dont work.

Code:
volScalarField V("V", U & mesh.C());
fvScalarMatrix UEqn
(
V-RaTk == fvc::grad(p) //Ratk is a VolScalarField
);



UEqn.solve();

Anybody got a hint?
Werner
TomB is offline   Reply With Quote

Old   October 11, 2011, 14:31
Arrow
  #3
Cyp
Senior Member
 
Cyprien
Join Date: Feb 2010
Location: Stanford University
Posts: 230
Rep Power: 9
Cyp is on a distinguished road
If you want to code a Darcy equation in OpenFOAM (or in whatever codes) it is very easy. In fact, it is much easier than a Navier-Stokes equation with PISO/SIMPLE algorithms...

The complete system is defined by :
1°) a continuity equation. Most of the time it is div (U) = 0
2°) a momentum equation, which is reduced to the Darcy's law in case of porous media and small Reynold's number. U = - K/mu*grad(P)

As in most of the (U,p) solvers, you must solve a pressure equation. It is derived by inserting the Darcy's law within the continuity equation :

div(U)=0=div(-K/mu*grad(p))=laplacian(-K/mu,p)

Consequently, in OpenFOAM you will code :

Code:
solve
(
        fvm::laplacian(K/mu,p)
);

U = -K/mu*grad(p);
U.correctBoundaryConditions();
In this example, the permeability K can be defined as :
- a dimensionedScalar (in case of isotropic homogeneous porous medium)
- a dimensionedTensor (in case of anisotropic homogeneous porous medium)
- a volScalarField (in case of isotropic heterogeneous porous medium or/and in case of velocy-dependance (Forchheimer correction for instance))
- a volTensorField (in case of anisotropic heterogeneous porous medium or/and in case of velocy-dependance)

If you have a source term in the continuity equation, for example :
div(U) = -q
Therefore the pressure equation becomes :
laplacian(-K/mu,p) = -q

I do not know your equation with your Ra*T*k term, but the method above is always available.

Regards,
Cyp
argonaut and bugmenot42 like this.
Cyp is offline   Reply With Quote

Old   October 11, 2011, 15:00
Default
  #4
Cyp
Senior Member
 
Cyprien
Join Date: Feb 2010
Location: Stanford University
Posts: 230
Rep Power: 9
Cyp is on a distinguished road
This answer is the same for the two others topics you have opened...
Cyp is offline   Reply With Quote

Old   October 11, 2011, 20:26
Default
  #5
New Member
 
Tom
Join Date: Sep 2011
Posts: 11
Rep Power: 5
TomB is on a distinguished road
Hi CyP

thanks very much for your guidance. I'll put some thought into it and reply again later.

The Ra*T*k term is there as we are looking at flow in porous media that is heated (from below) so the convection is the interesting thing.

Ra is rayliegh number
T is the temperature (im unsure if it will be Tref, or a 'local' Temp)
and k is just the unit vector in the vertical direction.

Cheers
TomB is offline   Reply With Quote

Old   October 11, 2011, 23:03
Default
  #6
New Member
 
Tom
Join Date: Sep 2011
Posts: 11
Rep Power: 5
TomB is on a distinguished road
Cyp,
I compiled a darcySimpleFoam based on your guidance.
i attached the solver code.

I based it on the simpleFoam as the name suggests, so it still has some of the simpleFoam headers etc.

I need to test it, see if it makes sense. say in a box of porous media with a pressure difference equal to rho*g*h from top to bottom,

From the way I have set up the pEqn, it doesnt seem to take any inputs, so how does it 'know' what the pressure difference is.. ?

Also, I understood from your post that I dont need a PISO/SIMPLE method to solver these equations.

What I do want though is a transient solver? Does this complicate things, or will the same code work?

Cheers
Attached Files
File Type: c darcySimpleFoam.C (2.4 KB, 34 views)
TomB is offline   Reply With Quote

Old   October 12, 2011, 01:58
Default
  #7
Cyp
Senior Member
 
Cyprien
Join Date: Feb 2010
Location: Stanford University
Posts: 230
Rep Power: 9
Cyp is on a distinguished road
Hi TomB/neph101,

Why are you persisting with PISO/SIMPLE loop ? If you want to base your darcy solver on an existing solver, you should start from laplacianFoam.

I suggest you to work in several step:
1°) Fisrt you code and compile a simple Darcy solver (as the one I wrote below)

1°bis) Add a simple scalar transport equation (see example in scalarTransportFoam)

2°) Then, think to the way to add the temperature-dependent term in your solver. if U = Ra*k*T -grad(p) therefore, div(U) = Ra*k&grad(T) - laplacian(p) = 0 ...

3°) Consider transient effect by starting from a compressible continuity equation. Most of the time we transform the density of the time-derivative to a pressure via the perfect-gas law. There is probably other way. From example, you can solve an steadyState darcy at each time step while you temperature is solved in transient manner..

Regards,
Cyp
Cyp is offline   Reply With Quote

Old   October 12, 2011, 06:17
Default
  #8
New Member
 
Werner
Join Date: Dec 2009
Location: Austria
Posts: 13
Rep Power: 7
nep101 is on a distinguished road
Hi Cyp, TomB!

First thanks for your help and tips.
But i dont get it.

My Equations for the thermal convection in a porous cylinder (dimensionless) is this:
ddt(T)+U&grad(T)-laplace(T)=0
Make it a little easier for testing:
U&grad(T)-laplace(T)=0,
div(U)=0 and
U+grad(p)-Ra*T*k=0

The Way i want to solve it is this:
1) solve the Tempeqn.
2) solve the continuity equation (this is what not work)
3) solve U

ad 1:
this is easy (i think)

Quote:
solve
(
(U&fvc::grad(T))
- fvm::laplacian(T)
);
ad 2:
with div(U)=0=div(Ra*k*T)-laplacian(p)
or like Cyp already wrote
div(U) = Ra*k&grad(T) - laplacian(p) = 0
Quote:
fvScalarMatrix pEqn
(
fvm::laplacian(p) == fvc::div(Ra*k*T)
);
pEqn.solve();
ad 3:
Quote:
U=Ra*k*T-fvc::grad(p);
But something with the continuity equation went wrong. I get a pressure arround e+16 . And Openfoam need 1000 interations in this step.

Cheers

Last edited by nep101; October 12, 2011 at 06:17. Reason: forget to say Cheers :)
nep101 is offline   Reply With Quote

Old   October 12, 2011, 06:49
Default
  #9
New Member
 
Tom
Join Date: Sep 2011
Posts: 11
Rep Power: 5
TomB is on a distinguished road
Hi Nep,

Why solve the Temp Eqn first? From what I understand, the scheme would be to solve first the pEqn, then UEqn and then TEqn.

I have been working on the scheme today, guided by Cyp's advice.

I have attached my solver.

I have issues with it,
1) I need to add the p_rgh (by that I mean pressure should increase with depth)
2) I need to add a TEqn
3) The UEqn needs the Ra*T term, which I still have to make a dictionary for and implement in the UEqn

I also want to do this all non-dimensionally, which is adding another level of confusion.

Cheers

Tom
Attached Files
File Type: c darcySimpleFoam.C (2.4 KB, 27 views)
TomB is offline   Reply With Quote

Old   October 12, 2011, 07:02
Default laplacianFoam
  #10
New Member
 
Tom
Join Date: Sep 2011
Posts: 11
Rep Power: 5
TomB is on a distinguished road
Hi Cyp,

Thanks for your input. I was confused by the using the piso/simple loop. As you guessed, it was just a convenience to base my solver on an existing one. I'll re-base it on laplacianFoam and then follow your advice.

In regards to the temperature, would you recommend taking guidance from bouyantBoussinesq(Simple or Pimple)Foam ?

I do want to include the boussinesq approx, and the next thing I need to figure out is how to include the change in pressure with depth.

Cheers

Tom
TomB is offline   Reply With Quote

Old   October 12, 2011, 07:52
Default
  #11
Cyp
Senior Member
 
Cyprien
Join Date: Feb 2010
Location: Stanford University
Posts: 230
Rep Power: 9
Cyp is on a distinguished road
Your problem is very easy to code ! much easier than Navier-Stokes equation. In fact, coding a Darcy's law can help you to understand the (U,p) solver for Navier-Stokes equations...

For the moment, we forget the temperature scalar transport equation. So your system is defined by
1°) a continuity equation
2°) a momentum equation

You have to notice that the unknow variable in the momentum equation is the velocity. At this point, there is no differential equation that govern the pressure field. However, the knowledge of this pressure field is very important since it appears as the main source term in the momentum equation.

You must also notice that the continuity equation is never directly solve in a (U,p) solver. In fact, we use this equation as a constraint to built the equation that govern p. The pressure equation laplacian(-K/mu,p) will always satisfy the continuity equation, even if we do not directly solve div(U) = 0

The first step of your development is to code (and understand how to code) a simple Darcy's law ! Then you could the pressure effect.

Another remark : the divergence operator in OpenFOAM require a velocity flux (something like div(phi) where phi = linearInterpolate(U) & mesh.Sf() )

I can help you, but you must proceed step by step.

@++
Cyp
Cyp is offline   Reply With Quote

Old   October 12, 2011, 07:56
Default
  #12
New Member
 
Werner
Join Date: Dec 2009
Location: Austria
Posts: 13
Rep Power: 7
nep101 is on a distinguished road
Hi Tom,

i think it is not so important in wich order you solve the equatios. My Problem is that something is going wrong with my continuity equation. After every interation my pressure is increases until Openfoam interrupts (at a point the pressure goes inf).

Quote:
solve(fvm::laplacian(p) - fvc::div(Ra*k*T));
My Equat. are also nondimensional. What you must do is set dimensionSet to 0 in opt/openfoamXXX/etc/controlDict. Or you will get problems with the Units.

Adding rgh to the pressure is not so difficult i think look at the creatFileds in boussinesqSimpleFoam:


Quote:
Info<< "Calculating field g.h\n" << endl;
volScalarField gh("gh", g & mesh.C());
surfaceScalarField ghf("ghf", g & mesh.Cf());

volScalarField p
(
IOobject
(
"p",
runTime.timeName(),
mesh,
IOobject::NO_READ,
IOobject::AUTO_WRITE
),
p_rgh + rhok*gh
);

label pRefCell = 0;
scalar pRefValue = 0.0;
setRefCell
(
p,
p_rgh,
mesh.solutionDict().subDict("PIMPLE"),
pRefCell,
pRefValue
);
nep101 is offline   Reply With Quote

Old   October 12, 2011, 09:21
Default
  #13
New Member
 
Werner
Join Date: Dec 2009
Location: Austria
Posts: 13
Rep Power: 7
nep101 is on a distinguished road
Hi Cyp,

Thanks for your input. But I really dont now how to code it.

Ok forget about the Temp so i got the 2 Equations:
div(U)=0 and U+grad(p)-Ra*k*T=0.
Ra is the Rayleighnumber and k is a Vector. I start with U=0 and heat a side of my object, so I can say my U is known and =0. So I must solve the pressure and satisfy the continuity equation.

I cannot solve div(U) the way i tried before.
Maybe with Rhie-Chow interpolation?
Like:
Ueqn = U - Ra*k*T
AU = UEqn().A();
U = UEqn().H()/AU
create phi
solve(fvm::laplacian(1.0/AU, p) == fvc::div(phi))
U=-grad(p)+Ra*k*T

What do you think Cyp or is this totally wrong????
nep101 is offline   Reply With Quote

Old   October 12, 2011, 09:27
Default
  #14
Cyp
Senior Member
 
Cyprien
Join Date: Feb 2010
Location: Stanford University
Posts: 230
Rep Power: 9
Cyp is on a distinguished road
Forget the Temperature, so your velocity is U = -grad p (according to your dimensionless analysis).

Forget also the predictor/corrector algorithm. Just think about mathematic with pen and paper.

The first step is to code :
div(U) = 0
U = -grad(p)

I mentioned all the stages in my previous posts. At this point you may also forget the "phi" part.

let me know of your progression.

@++
Cyp
Cyp is offline   Reply With Quote

Old   October 12, 2011, 09:43
Default
  #15
New Member
 
Tom
Join Date: Sep 2011
Posts: 11
Rep Power: 5
TomB is on a distinguished road
Hi Nep
Thanks.
Can you tell me how you set up the k vector?
Cheers

Tom
TomB is offline   Reply With Quote

Old   October 12, 2011, 11:01
Default
  #16
New Member
 
Werner
Join Date: Dec 2009
Location: Austria
Posts: 13
Rep Power: 7
nep101 is on a distinguished road
Hi Cyp,

Quote:
solve
(
fvm::laplacian(p)
);

U=-fvc::grad(p);
I dont get it.
You said U=-grad(p) so with div(U)=0 you get laplacian(p)=0 and U=-grad(p).
But in my case i got from div(U)=0 laplacian(p)=div(Ra*k*T)=Ra*d(T)/dz.

Hi Tom,

I defined k as a volVectorField [0 0 1]

Cheers
Werner
nep101 is offline   Reply With Quote

Old   October 12, 2011, 11:10
Default
  #17
Cyp
Senior Member
 
Cyprien
Join Date: Feb 2010
Location: Stanford University
Posts: 230
Rep Power: 9
Cyp is on a distinguished road
Indeed in your case you have another term in your equation. However it is easy to add it when you have a pressure/velocity solver that works fine !

Now that you have your velocity/pressure solver (if we neglect the temperature term) you can add the advection-diffusion temperature equation :

Code:
solve
(
    fvm::ddt(T) + fvm::div(phi,T) - fvm::laplacian(DT,T)
);
in OpenFOAM, "phi" is the velocity flux, defined as
Code:
phi = linearInterpolate(U) & mesh.Sf()
You must add this line just after the calculation of U

When you will have check the validity of your solver, you can add the temperature correction in the calculation of U and p.
Code:
Ra*k&fvc::grad(T)
Don't forget that you always have to work step by step when you code a solver. It is too difficult (and source of error) to program all your problem in one time...

@++
Cyp
Cyp is offline   Reply With Quote

Old   October 12, 2011, 12:05
Default
  #18
New Member
 
Werner
Join Date: Dec 2009
Location: Austria
Posts: 13
Rep Power: 7
nep101 is on a distinguished road
Hi Cyp,

I realy thank you for your help, but I think I give up. It is allwasy the same and I dont know were the problem is.

My Code (another attempt):
Quote:
volVectorField kT("kt",Ra*k*T);
volScalarField divT = fvc::grad(kT)().component(vector::Z);
solve
(
fvm::laplacian(p)- divT
);
U=-fvc::grad(p)+Ra*k*T;
phi = linearInterpolate(U) & mesh.Sf();
solve
(
fvm::ddt(T) + fvm::div(phi,T) - fvm::laplacian(T)
);
But it is alltime the same. After one or two interations the pressure explode (e+12) and thats it.

I cant compute the pressure right. So once again Cyp thanks for your time and your help.

Cheers
Werner
nep101 is offline   Reply With Quote

Old   October 12, 2011, 12:11
Default
  #19
Cyp
Senior Member
 
Cyprien
Join Date: Feb 2010
Location: Stanford University
Posts: 230
Rep Power: 9
Cyp is on a distinguished road
Quote:
Originally Posted by nep101 View Post
Hi Cyp,

I realy thank you for your help, but I think I give up. It is allwasy the same and I dont know were the problem is.

My Code (another attempt):


But it is alltime the same. After one or two interations the pressure explode (e+12) and thats it.

I cant compute the pressure right. So once again Cyp thanks for your time and your help.

Cheers
Werner

What are the boundary condition you applied ?

And I am not sure about the way you code divT.. I think you should use the snippet I previously proposed.

And you should also use U.correctBoundaryConditions() after the calculation of U since U is evaluated from another field
Cyp is offline   Reply With Quote

Old   October 12, 2011, 13:16
Default
  #20
New Member
 
Werner
Join Date: Dec 2009
Location: Austria
Posts: 13
Rep Power: 7
nep101 is on a distinguished road
Hi Cyp,
i use a simple cube for testing.

T:
Quote:
internalField uniform 450;

boundaryField
{

seite
{
type zeroGradient;
}
fab
{
type zeroGradient;
}


unten
{
type fixedValue;
value uniform 500;
}
oben
{
type fixedValue;
value uniform 450;
}
p:
Quote:
internalField uniform 0;

boundaryField
{

seite
{
type zeroGradient;
}
fab
{
type zeroGradient;
}


unten
{
type zeroGradient;
}
oben
{
type zeroGradient;
}
}
and U:
Quote:
internalField uniform (0 0 0);

boundaryField
{

seite
{
type fixedValue;
value uniform (0 0 0);
}
fab
{
type fixedValue;
value uniform (0 0 0);
}


unten
{
type fixedValue;
value uniform (0 0 0);
}
oben
{
type fixedValue;
value uniform (0 0 0);
}
}
with oben is the top and unten is the bottom, fab and seite are the side of the cube.

Cheers
Werner
nep101 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
Boussinesq Assumption tstorm FLUENT 3 October 26, 2009 06:11
Boussinesq lucifer FLUENT 0 September 26, 2009 09:31
Outflow boundary condition for a boussinesq fluid lin OpenFOAM Running, Solving & CFD 19 July 31, 2009 08:50
Natural Convection in a storage tank - boussinesq? Jervds CFX 5 November 19, 2008 07:01
Darcy + multiphase? George Bergantz CD-adapco 1 March 14, 2002 06:07


All times are GMT -4. The time now is 23:51.