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

Problems in understanding BuoyantBoussinesqSimpleFoam

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

Like Tree6Likes
  • 5 Post By Anne Lincke
  • 1 Post By makaveli_lcf

Reply
 
LinkBack Thread Tools Display Modes
Old   August 24, 2011, 11:41
Default Problems in understanding BuoyantBoussinesqSimpleFoam
  #1
Senior Member
 
Anne Gerdes
Join Date: Aug 2010
Location: Hamburg
Posts: 152
Rep Power: 6
Anne Lincke is on a distinguished road
Dear Foamers,

I still habe problems in understanding the source code of BuoyantBoussinesqSimpleFoam.

In the attachment you see the momentum equation which has to be solved where GT= rho g= (1-beta(T-T_0))*g

In pEqn.H the main code is settled.

First the velocity is solved without taking into account pressure or density.
HTML Code:
U = rAU*UEqn().H();
The flux is the interpolation of this velocity.
HTML Code:
phi = fvc::interpolate(U) & mesh.Sf();
Then a buoyancy flux is defined

HTML Code:
surfaceScalarField buoyancyPhi(rAUf*ghf*fvc::snGrad(rhok)*mesh.magSf());
Here my first question: Why is the normal gradient computed? In the equation there is only a term which includes the density, not the gradient of the density.

The flux is corrected with the buoyancy flux

HTML Code:
 phi -= buoyancyPhi;
After the pressure correction with the Laplace equation

HTML Code:
fvm::laplacian(rAUf, p_rgh) == fvc::div(phi)
the flux is again corrected by taking into account the pressure minus rgh.

HTML Code:
phi -= p_rghEqn.flux();
In the last step the velocity is reconstructed from the flux

HTML Code:
U -= rAU*fvc::reconstruct((buoyancyPhi + p_rghEqn.flux())/rAUf);
and p_rgh is updated

HTML Code:
 p_rgh = p - rhok*gh;
I do not understand the components of the flux. If they are added we obtain

HTML Code:
phi = fvc::interpolate(U) & mesh.Sf() - (rAUf*ghf*fvc::snGrad(rhok)*mesh.magSf()) - p_rghEqn.flux();
So the term including the normal gradient of the density is too much.
Regarding the link
http://openfoamwiki.net/index.php/Bu...sinesqPisoFoam

in the solver BuoyantBoussinesqPisoFoam the flux looks like this

HTML Code:
phi =fvc::interpolate(U) & mesh.Sf()) + rUAf*fvc::interpolate(rhok)*(g & mesh.Sf())-pEqn.flux();
which makes much more sense to me.
I can not find any explanation in the forum. Who can help me to understand the code?
Attached Images
File Type: png Bildschirmfoto.png (9.6 KB, 100 views)
Anne Lincke is offline   Reply With Quote

Old   September 5, 2011, 08:58
Default
  #2
Senior Member
 
Anne Gerdes
Join Date: Aug 2010
Location: Hamburg
Posts: 152
Rep Power: 6
Anne Lincke is on a distinguished road
No one for an answer???
Anne Lincke is offline   Reply With Quote

Old   September 5, 2011, 09:39
Default
  #3
Senior Member
 
romant's Avatar
 
Roman Thiele
Join Date: Aug 2009
Location: Stockholm, Sweden
Posts: 359
Rep Power: 11
romant is on a distinguished road
I am not sure, but you might have to look up in a numerical methods book, or somewhere else how the SIMPLE algorithm works, maybe that can give you some hints.
__________________
~roman
romant is offline   Reply With Quote

Old   September 8, 2011, 09:20
Default
  #4
Senior Member
 
Anne Gerdes
Join Date: Aug 2010
Location: Hamburg
Posts: 152
Rep Power: 6
Anne Lincke is on a distinguished road
Thank you for answering!
Now I understand the implementation. If someone is interested I could post that, too.
Anne Lincke is offline   Reply With Quote

Old   September 8, 2011, 11:24
Default
  #5
Senior Member
 
Bernhard
Join Date: Sep 2009
Location: Delft
Posts: 790
Rep Power: 12
Bernhard is on a distinguished road
If it is not on this forum yet, then it can be very interesting for future reference if you post your findings here!
Bernhard is offline   Reply With Quote

Old   September 12, 2011, 09:02
Default
  #6
Senior Member
 
Anne Gerdes
Join Date: Aug 2010
Location: Hamburg
Posts: 152
Rep Power: 6
Anne Lincke is on a distinguished road
In the momentum equation we have in z-direction (in direction of buoyancy) the terms

- dp / dz + rho*g

In OpenFOAM, g is a vector (0 0 -9,81) which ensures that the buoyancy is only valid for the right coordinate direction.

In order to guarantee, that in the pressure correction of the simple algorithm also the buoyancy term rho*g*z is taken into account, the pressure and the buoyancy are melted together in one term

p_rgh = p - rho*g*z

Instead of the normal gradient of p, the normal gradient of p_rgh is on the RHS.

- d/dn (p_rgh) which equals in OF to

HTML Code:
 - fvc :: snGrad (p_rgh)
Computing the derivative of this term with product rule and assuming that rho may change with respect to the z-direction we obtain

- d/dz [ p - rho * g *z] = - dp/dz + rho*g + g*z* d rho/dz

So the third term g*z* d rho/dz is "too much" and is therefore substracted via

HTML Code:
- ghf *fvc::snGrad(rhok)
So the derivative of the density is a correction term which is nonzero for changes in density.
Anne Lincke is offline   Reply With Quote

Old   September 14, 2011, 02:08
Default
  #7
Senior Member
 
Dr. Alexander Vakhrushev
Join Date: Mar 2009
Posts: 213
Rep Power: 10
makaveli_lcf is on a distinguished road
Send a message via ICQ to makaveli_lcf
Hallo Anne!

Would you please clarify following points:

1. What do you mean by term "two much" (is it phisical meaning or regards difference in some mathematical formulation).

2. Did you use some reference to understand SIMPLE/PISO treatment of the buoyancy?

Thank you in advance!
Mfg,
Alexander
Tobi likes this.
__________________
Best regards,

Dr. Alexander VAKHRUSHEV

Christian Doppler Laboratory for "Advanced Process Simulation of
Solidification and Melting"

Simulation and Modelling of Metallurgical Processes
Department of Metallurgy
University of Leoben

Franz-Josef-Str. 18
A - 8700 Leoben
Österreich / Austria
Tel.: +43 3842 - 402 - 3125
http://smmp.unileoben.ac.at
makaveli_lcf is offline   Reply With Quote

Old   September 14, 2011, 03:47
Default
  #8
Senior Member
 
Anne Gerdes
Join Date: Aug 2010
Location: Hamburg
Posts: 152
Rep Power: 6
Anne Lincke is on a distinguished road
Hallo Alexander,

I will try to answer your questions:

1. With "too much" I mean the term g*z* d rho/dz which is the result of taking the derivative of - [ p - rho * g *z] in z-direction. It is an additional term due to the product rule. Basically we have the terms
- dp / dz + rho*g

So by taking the derivative of both, p and rho*g*z we obtain together with the correction term - g*z* d rho/ dz:

- d/dz [ p - rho * g *z] - g*z* d rho/ dz=
- dp /dz + rho*g + g*z* d rho/ dz - g*z* d rho/ dz =
- dp/ dz + rho*g

and this is the orignial term in the buoyancy driven momentum equation.
In OpenFOAM- code this equals to the right-hand side of the UEqn

HTML Code:
 - fvc::snGrad(p_rgh) - ghf*fvc::snGrad(rhok)
2. I have a lecture script where this method is described, but it is not open for everyone. I also found a notice for this method in Peric and Ferziger's "Computational Methods for Fluid Dynamics" at the very beginning of the book (page 11 in the version I have). Search for "hydrostatic pressure" and you will find a short explanation there.
Anne Lincke is offline   Reply With Quote

Old   September 14, 2011, 04:16
Default
  #9
Senior Member
 
Dr. Alexander Vakhrushev
Join Date: Mar 2009
Posts: 213
Rep Power: 10
makaveli_lcf is on a distinguished road
Send a message via ICQ to makaveli_lcf
Thank you Anne for the detailed explanation! I will search for the algorithm description in Peric's book.

Got it, thanx)
__________________
Best regards,

Dr. Alexander VAKHRUSHEV

Christian Doppler Laboratory for "Advanced Process Simulation of
Solidification and Melting"

Simulation and Modelling of Metallurgical Processes
Department of Metallurgy
University of Leoben

Franz-Josef-Str. 18
A - 8700 Leoben
Österreich / Austria
Tel.: +43 3842 - 402 - 3125
http://smmp.unileoben.ac.at
makaveli_lcf 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
[ICEM] Problems with coedge curves and surfaces tommymoose ANSYS Meshing & Geometry 0 August 5, 2011 16:02
Needed Benchmark Problems for FSI Mechstud Main CFD Forum 4 July 26, 2011 12:13
Two-phase air water flow problems by activating Wall Lubrication Force challenger85 CFX 5 November 5, 2009 06:44
Problems understanding some piso details tehache OpenFOAM Running, Solving & CFD 3 July 27, 2007 06:02
Some problems with Star CD Micha CD-adapco 0 August 6, 2003 13:55


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