# Need Help w. Darcy Boussinesq Equ.

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

 August 22, 2011, 10:13 Need Help w. Darcy Boussinesq Equ. #1 New Member   Werner Join Date: Dec 2009 Location: Austria Posts: 13 Rep Power: 9 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

September 26, 2011, 00:56
#2
New Member

Tom
Join Date: Sep 2011
Posts: 11
Rep Power: 8
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 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

 October 11, 2011, 14:31 #3 Senior Member   Cyprien Join Date: Feb 2010 Location: Stanford University Posts: 260 Rep Power: 11 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, bugmenot42 and tomdylan like this.

 October 11, 2011, 15:00 #4 Senior Member   Cyprien Join Date: Feb 2010 Location: Stanford University Posts: 260 Rep Power: 11 This answer is the same for the two others topics you have opened...

 October 11, 2011, 20:26 #5 New Member   Tom Join Date: Sep 2011 Posts: 11 Rep Power: 8 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

October 11, 2011, 23:03
#6
New Member

Tom
Join Date: Sep 2011
Posts: 11
Rep Power: 8
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
 darcySimpleFoam.C (2.4 KB, 36 views)

 October 12, 2011, 01:58 #7 Senior Member   Cyprien Join Date: Feb 2010 Location: Stanford University Posts: 260 Rep Power: 11 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

October 12, 2011, 06:17
#8
New Member

Werner
Join Date: Dec 2009
Location: Austria
Posts: 13
Rep Power: 9
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:
Make it a little easier for testing:
div(U)=0 and

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

this is easy (i think)

Quote:
 solve ( (U&fvc::grad(T)) - fvm::laplacian(T) );
with div(U)=0=div(Ra*k*T)-laplacian(p)
div(U) = Ra*k&grad(T) - laplacian(p) = 0
Quote:
 fvScalarMatrix pEqn ( fvm::laplacian(p) == fvc::div(Ra*k*T) ); pEqn.solve();
Quote:
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 :)

October 12, 2011, 06:49
#9
New Member

Tom
Join Date: Sep 2011
Posts: 11
Rep Power: 8
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
 darcySimpleFoam.C (2.4 KB, 30 views)

 October 12, 2011, 07:02 laplacianFoam #10 New Member   Tom Join Date: Sep 2011 Posts: 11 Rep Power: 8 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

 October 12, 2011, 07:52 #11 Senior Member   Cyprien Join Date: Feb 2010 Location: Stanford University Posts: 260 Rep Power: 11 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

October 12, 2011, 07:56
#12
New Member

Werner
Join Date: Dec 2009
Location: Austria
Posts: 13
Rep Power: 9
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 );

 October 12, 2011, 09:21 #13 New Member   Werner Join Date: Dec 2009 Location: Austria Posts: 13 Rep Power: 9 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????

 October 12, 2011, 09:27 #14 Senior Member   Cyprien Join Date: Feb 2010 Location: Stanford University Posts: 260 Rep Power: 11 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

 October 12, 2011, 09:43 #15 New Member   Tom Join Date: Sep 2011 Posts: 11 Rep Power: 8 Hi Nep Thanks. Can you tell me how you set up the k vector? Cheers Tom

October 12, 2011, 11:01
#16
New Member

Werner
Join Date: Dec 2009
Location: Austria
Posts: 13
Rep Power: 9
Hi Cyp,

Quote:
I dont get it.
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

 October 12, 2011, 11:10 #17 Senior Member   Cyprien Join Date: Feb 2010 Location: Stanford University Posts: 260 Rep Power: 11 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

October 12, 2011, 12:05
#18
New Member

Werner
Join Date: Dec 2009
Location: Austria
Posts: 13
Rep Power: 9
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

October 12, 2011, 12:11
#19
Senior Member

Cyprien
Join Date: Feb 2010
Location: Stanford University
Posts: 260
Rep Power: 11
Quote:
 Originally Posted by nep101 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

October 12, 2011, 13:16
#20
New Member

Werner
Join Date: Dec 2009
Location: Austria
Posts: 13
Rep Power: 9
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

 Thread Tools Display Modes Linear Mode

 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 OffTrackbacks are On Pingbacks are On Refbacks are On Forum Rules

 Similar Threads Thread Thread Starter Forum Replies Last Post Jervds CFX 6 October 2, 2016 02:53 tstorm FLUENT 3 October 26, 2009 06:11 lucifer FLUENT 0 September 26, 2009 09:31 lin OpenFOAM Running, Solving & CFD 19 July 31, 2009 08:50 George Bergantz Siemens 1 March 14, 2002 06:07

All times are GMT -4. The time now is 19:45.