
[Sponsors] 
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: URa*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 { URaTk } But it dont work. Code: volScalarField V("V", U & mesh.C()); fvScalarMatrix UEqn ( VRaTk == 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:


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 NavierStokes 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();  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 velocydependance (Forchheimer correction for instance))  a volTensorField (in case of anisotropic heterogeneous porous medium or/and in case of velocydependance) 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 

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 

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 temperaturedependent 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 timederivative to a pressure via the perfectgas 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: 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:
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:
Quote:
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 nondimensionally, which is adding another level of confusion. Cheers Tom 

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 rebase 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 NavierStokes equation. In fact, coding a Darcy's law can help you to understand the (U,p) solver for NavierStokes 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:
Adding rgh to the pressure is not so difficult i think look at the creatFileds in boussinesqSimpleFoam: Quote:


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

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 advectiondiffusion temperature equation : Code:
solve ( fvm::ddt(T) + fvm::div(phi,T)  fvm::laplacian(DT,T) ); Code:
phi = linearInterpolate(U) & mesh.Sf() 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) @++ 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:
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:
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:
Quote:
Quote:
Cheers Werner 

Thread Tools  
Display Modes  


Similar Threads  
Thread  Thread Starter  Forum  Replies  Last Post 
Natural Convection in a storage tank  boussinesq?  Jervds  CFX  6  October 2, 2016 02:53 
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 
Darcy + multiphase?  George Bergantz  Siemens  1  March 14, 2002 06:07 