Need Help w. Darcy Boussinesq Equ.
 

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

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 !

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 :

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

This answer is the same for the two others topics you have opened...

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

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

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

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:mad:)
3) solve U

ad 1:
this is easy (i think)

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
ad 3:
But something with the continuity equation went wrong. I get a pressure arround e+16 . And Openfoam need 1000 interations in this step. :mad::mad::mad:

Cheers

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

laplacianFoam
 

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

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

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).

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:

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????

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

Hi Nep
Thanks.
Can you tell me how you set up the k vector?
Cheers

Tom

Hi Cyp,

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

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 :

in OpenFOAM, "phi" is the velocity flux, defined as
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.
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

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

Hi Cyp,
i use a simple cube for testing.

T:
p:
and U:
with oben is the top and unten is the bottom, fab and seite are the side of the cube.

Cheers
Werner

Try to impose a delta P (fixedValue for inlet and fixedValue = 0 for outlet).

If you want a immobile system, I think a zeroGradient for inlet and fixedValue = 0 for outlet should work...

Note that since we solve only the pressure and the temperature in the solver, the boundary condition for the velocity field are not accounted for.

The problem we try to solve is this onset of thermal convection in a porous cylinder with conducting wall. So there is no inlet or outlet.
I also thougt that there is maybe a problem with the boundaries. In most cases when i got a floiting point exception, it was a problem with the bopundaries.

But the mainproblem is (i think), that the pressureequation is not properly solved or coded.

Cheers

Gday Cyp and Nep

Thanks so much for your input and collaboration, this has been a huge help to me.

I re-wrote my solver, based on the laplacianFoam. I wrote up some comments along the way.

I'd like to look into it to 'see' how pressure and velocity are linked. That is what I will do next, to make sure its working as expected at this stage.

I'm going to try to set up a pressure gradient with the boundary conditions as you two were discussing.

my method is
break the problem into a number of smaller steps on Cyps guidance
1. solve simple Darcy's Law, based on laplacianFoam (exclude the Temp term)
2. include a scalar transport equation
3. then look at adding temperature
4. then look at correcting T,U
5. add the (Ra*T*k) term to the U eqn
Cheers

Tom

Hi Tom,

wich equation you want to solve?
Is it necessary to handle the equations in a dimensionless form.
What i try to say is, that MAYBE (im no Pro in Openfoam) it is not so easy to handle nondimensionless equations in Openfoam.

Btw how is it going on??

Cheers
Werner

Hi Werner

Im solving for convection in a box of porous media. Sounds very similar to your problem otherwise. I do want to do it in terms of the Rayleigh number and non-dimensional temperature.

I will fudge the units if i have to ! it seems to work ok with the dimensions switched off in the /opt location. I just dont like fiddling with those settings bcoz im an amatuer.

At the moment, i am trying to make the darcyLaplacian solver with some temperature terms added. I seem to be missing a curly bracket, but I'm not sure where.

Did you get a pressure gradient working?

Cheers

Tom
edit

I added my solver and the TEqn it needs. Im having issues now with boundaryConditions. It doesnt like taking the gradient of a calculated internalField.

HI Tom,

no i dont got my continuity equation working correctly.

Send me your solver, maybe i find your problem.

Cheers
Werner

HI Werner

my solver is in the post above. I added another curly bracket and it works, but it doesnt make sense. i use bracket matching and indenting so i should be able to see where its missin..

I think my set up needs a setFields, so that the initial pressure is set before it starts to solve the fvm::laplacian(p)

Hi Tom,

which vers. of OF you use?

Hi Nep

Im using OpenFoam 2.0.1

Im having an issue with the boundary. I get the error message,

gradientInternalCoeffs cannot be called for a calculatedFvPatchField.

it suggests I migh be using a default boundary condition. I copied the boundary conditions from buoyantBoussinesqSimpleFoam

I'll work on this some more tmw ( about 8 hours from now zzzzzzzzzzzzzzzz)

I am just wondering whether this old thread has ever been solved. I think my thread may need a similar solution.

Cheers,
Marcus

DarcyFoam with Boussinesq assumption pressure diverges.
 

Hello All
I am trying to include heat transfer in DarcyFoam solver available in porousMultiphaseFoam (https://github.com/phorgue/porousMul...ties/darcyFoam.

The solver is intended to simulate Natural convective flow in porous media filled cavity. So I would have zero velocity at all walls. Two of the walls would have a temperature (different from each other) with remaining being adiabatic. In the presence of gravitational force a flow is expected to be induced.

Intended equation to be solved:
where

• u=velocity (m/s)
• K=permeability (m^2)
• mu=viscosity (kg/ms)
• p=pressure (kg/ms^2)
• g=gravity vector^(m/s^2)
• rho_0=reference density (kg/m3)
• beta=coeffcient of volume expansion (1/K)

In order to have such a solver, I modified the Main .C file as:
The createFields.H is modified as:

When I run the solver on a test case (attached herewith), the pressure diverges.
Test conditions tried and outcome

1. gravity set to zero--code works and gives a conduction profile
2. Equal temp at both walls with non zero gravity--Expected uniform temperature profile--outcome--core dumped
3. Two walls at different temperature, other wall adiabatic, gravity set to a value, coefficient of expansion beta in transport properties set to a value--outcome--core dumped
   
   Code:

   ---

   DICPCG:  Solving for p, Initial residual = 1, Final residual = 1.2601e+08, No Iterations 1000
Vel done
Temp done
DILUPBiCG:  Solving for T, Initial residual = 0.565344, Final residual = 1.85342e-12, No Iterations 1
DICPCG:  Solving for p, Initial residual = 1, Final residual = 5.48003e+10, No Iterations 1000
Vel done
Temp done
ExecutionTime = 1.51 s  ClockTime = 2 s

Time = 4

DILUPBiCG:  Solving for T, Initial residual = 0.818783, Final residual = 8.10268e-08, No Iterations 7
DICPCG:  Solving for p, Initial residual = 1, Final residual = 8.65643e+08, No Iterations 1000
Vel done
Temp done

   ---

Important parameters to be changed for testing

• constant/g
• constant/transportProperties/beta (beta is the volume expansion coeffcient, that induces a forcing term equivalent to density change induced effect.

Please note: to calculate rhok I tried rhof*(1.0 - (beta*(T - TRef))) and rhof*(0.0 - (beta*(T - TRef))) depending on whether p is treated as pressure or pressure-rho*g*h. But both of them gave similar error.


Can someone please suggest what could be the issue and how I should approach the same.

please find the solver and testcase attached.