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

BubbleFoam - Source Term

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

Reply
 
LinkBack Thread Tools Display Modes
Old   September 30, 2011, 03:50
Default BubbleFoam - Source Term
  #1
Senior Member
 
Aurelien Thinat
Join Date: Jul 2010
Posts: 154
Rep Power: 7
Aurelien Thinat is on a distinguished road
Hi everybody,

I need some advice in adding a source term for the dispersed phase in the bubbleFoam solver.

At first I put a volscalarfield named "source" in my domain and I transformed the alpha equation :

fvScalarMatrix alphaEqn
(
fvm::ddt(alpha)
+ fvm::div(phi,alpha,scheme)
+ fvm::div(-fvc::flux(-phir,beta,scheme),alpha,scheme)
=
source
)

But I had a loss of "beta" in my calculations :
A column of water (z<z0)with a volumic source term of air bubbles at the bottom. The top of the column is filled with air (z>z0).
After a few time, the water column collapses on itself.

I guess that's because of the definition beta = scalar(1) - alpha. If I create some alpha, I delete some beta right ?

So I would need some way to add this source term of alpha keeping the same mass of beta in the whole domain. Is that possible without rewriting a new solver ?

Thanks.
Aurelien Thinat is offline   Reply With Quote

Old   September 30, 2011, 05:50
Default
  #2
Senior Member
 
Aurelien Thinat
Join Date: Jul 2010
Posts: 154
Rep Power: 7
Aurelien Thinat is on a distinguished road
Ok, I got the answer.

I have modified the alpha equation the right way, but I didn't do the same on the pressure equation...

Thanks Cyp.
Aurelien Thinat is offline   Reply With Quote

Old   October 20, 2011, 10:35
Default
  #3
Senior Member
 
Aurelien Thinat
Join Date: Jul 2010
Posts: 154
Rep Power: 7
Aurelien Thinat is on a distinguished road
Hello everybody,

I'm still missing something.

I have modified my phase continuity equation as written above :
fvScalarMatrix alphaEqn
(
fvm::ddt(alpha)
+ fvm::div(phi,alpha,scheme)
+ fvm::div(-fvc::flux(-phir,beta,scheme),alpha,scheme)
=
source
)

So I need to modify my pressure equation right ?
To obtain the pressure equation, one should force div(phi)=0 in the case with no source term, which becomes div(phi) = source in my case, am I still right ?
Then the final pressure equation equation should be :
pEqn
(
fvm::laplacian(Dp,p) == fvc::div(phi)-source
)

EDIT :
Well in fact I get lost in Henrik Rusche Thesis at p.112 :
He wrote that the mixture continuity equation is the equation 3.39 by combining the 2 phase continuity equation with weightning coefficients wi.
In my case the 0 should become "rhoa*source/wa".

But then he decides to solve the volumetric mixture continuity equation which is div(U)=0.
Since he got this equation by combining the 2 phase continuity equations, I guess I still should write :
div(U)= div(alpha*Ua + beta*Ub)=source
And that leads to the pressure equation :
pEqn
(
fvm::laplacian(Dp,p) == fvc::div(phi)-source
)
But when I code this, I have an infinite creation of water in my domain : the water level is rising and the mean volumic ratio of the dispersed phase is decreasing...

So I'm missing something but I don't know what. If someone has any clue...Thanks.

Last edited by Aurelien Thinat; October 21, 2011 at 03:33.
Aurelien Thinat is offline   Reply With Quote

Old   October 21, 2011, 11:38
Default
  #4
Senior Member
 
Aurelien Thinat
Join Date: Jul 2010
Posts: 154
Rep Power: 7
Aurelien Thinat is on a distinguished road
Hi again,

I continue on my problem. After few tries, I got that I need a coefficient next to my source term in the pressure equation : I have test manually few coefficients for a given value of "source" and I found one that allow mass conservation of the dispersed/continue phase.


In his thesis, Henrik Rusche wrote :
"we chose to solve for the mixture pressure by recasting the volumetric mixture continuity equation (3.3) into one for pressure. To achieve this, it is formulated at the cell faces, such that:
div(alpha_af * phi_a + alpha_bf * phi_b) = 0 "

So, the pressure equation is written at cell's faces only. Does it mean I have to use some kind of "phi_source", a projection of the source term on each face ?


EDIT :

I'm trying to interpolate my "volScalarField source" at the cell's face. The volScalarField alpha is interpolate at the cell's face by the command :
"surfaceScalarField alphaf(fvc::interpolate(alpha));"

So I defined :
"surfaceScalarField sourcef(fvc::interpolate(source));"


But then, when I try to include sourcef in the pEqn equation, I have an issue about their type with the operator "+" :

pEqn.H:54: error: no match for 'operator+' in
'Foam::fvc::div(const Foam::GeometricField<Type, Foam::fvsPatchField, Foam::surfaceMesh>&) [with Type = double]()
+
Foam::fvc::interpolate(const Foam::GeometricField<Type, Foam::fvPatchField, Foam::volMesh>&) [with Type = double]()'


If anyone has an idea...

Last edited by Aurelien Thinat; October 24, 2011 at 04:30.
Aurelien Thinat is offline   Reply With Quote

Old   October 24, 2011, 08:15
Default
  #5
Senior Member
 
Aurelien Thinat
Join Date: Jul 2010
Posts: 154
Rep Power: 7
Aurelien Thinat is on a distinguished road
No one has ever try to add a source term of alpha in this solver ?

A little summary (and a up of the post) :

1) fvc::div(phi) is a volScalarField,
2) fvm::laplacian(Dp,p) is a volScalarField too,
3) The initial pressure equation is
fvm::laplacian(Dp,p)=fvc::div(phi)
as the continuity equation at the face center is :
div(alpha_f*phia + beta_f*phib)=0
4) I add a source term in alphaEqn :
fvm::ddt(alpha)
+ fvm::div(phi,alpha,scheme)
+ fvm::div(-fvc::flux(-phir,beta,scheme),alpha,scheme)
==
source
5) My continuity equation at the face is so :
div(alpha_f*phia + beta_f*phib)=source_f (if anyone can confirm I am right)
I defined sourcef as fvc::interpolate(source). It is a surfaceScalarField.

6) How can I add it in the fvScalarMatrix pEqn ?

Last edited by Aurelien Thinat; November 2, 2011 at 05:21.
Aurelien Thinat is offline   Reply With Quote

Old   November 2, 2011, 11:21
Default
  #6
Senior Member
 
Alberto Passalacqua
Join Date: Mar 2009
Location: Ames, Iowa, United States
Posts: 1,894
Rep Power: 26
alberto will become famous soon enoughalberto will become famous soon enough
Hi,

consider the two continuity equations:

ddt(alpha) + div(alpha*Ua) = Sa
ddt(beta) + div(beta*Ub) = 0

Since they are volumetric already (we assumed constant density, so the equations were divided by rhoa and rhob respectively), there is no need of weighting again.

Sum the equations to obtain:

div(alpha*Ua) + div(beta*Ub) = Sa

Replace the velocities with face fluxes, and you have what you are looking for.

Best,
__________________
Alberto Passalacqua

GeekoCFD - A free distribution based on openSUSE 64 bit with CFD tools, including OpenFOAM. Available as live DVD/USB, hard drive image and virtual image.
OpenQBMM - An open-source implementation of quadrature-based moment methods
alberto is offline   Reply With Quote

Old   November 2, 2011, 11:37
Default
  #7
Senior Member
 
Aurelien Thinat
Join Date: Jul 2010
Posts: 154
Rep Power: 7
Aurelien Thinat is on a distinguished road
Thank you for answering Alberto.

So if I understand what you are saying I will have :
div(alpha_f*phia + beta_f*phib) - laplacian(Dp,p)=Sa (I take your notation)

That's the equation I have coded at first but with this one, the mean volume fraction in the whole domain is not stable. It leads to an endless creation of continuous phase.

Last edited by Aurelien Thinat; November 2, 2011 at 12:33.
Aurelien Thinat is offline   Reply With Quote

Old   November 2, 2011, 11:44
Default
  #8
Senior Member
 
Alberto Passalacqua
Join Date: Mar 2009
Location: Ames, Iowa, United States
Posts: 1,894
Rep Power: 26
alberto will become famous soon enoughalberto will become famous soon enough
Quote:
Originally Posted by Aurelien Thinat View Post
Thank you for answering Alberto.

So if I understand what you are saying I will have :
div(alpha_f*phia + beta_f*phib) - laplacian(Dp,p)=Sa (I take your notation)

That's the equation I have coded at first but with this one, the mean volume fraction in the whole domain is not stable. It leads to an endless creation of continuous phase.
That's the equation that comes out from what you said. However the formulation of the model seems a bit inconsistent. If the mass leaves phase a, where does it go? In deriving the pressure equation you are imposing that the total mass (volume, since rho = const) is conserved.
__________________
Alberto Passalacqua

GeekoCFD - A free distribution based on openSUSE 64 bit with CFD tools, including OpenFOAM. Available as live DVD/USB, hard drive image and virtual image.
OpenQBMM - An open-source implementation of quadrature-based moment methods
alberto is offline   Reply With Quote

Old   November 2, 2011, 12:00
Default
  #9
Senior Member
 
Aurelien Thinat
Join Date: Jul 2010
Posts: 154
Rep Power: 7
Aurelien Thinat is on a distinguished road
I'm trying to simulate the case of bubble plume. So I have a sharp interface between the dispersed phase a and the continuous phase b.
At t= 0 :
- In the top part of the domain I have alpha = 1,
- In the bottom part, alpha = 0.
The dispersed phase a can go through an outlet at the top of the domain. The creation should be ensured by the source term I'm trying to implement (I'd like to simulate an air bubble creation from a chemical reaction).

Normally, my interface should be moved to the top and my volume fraction a should be constant after few seconds (well, less than seconds).
Aurelien Thinat is offline   Reply With Quote

Old   November 2, 2011, 12:31
Default
  #10
Senior Member
 
Alberto Passalacqua
Join Date: Mar 2009
Location: Ames, Iowa, United States
Posts: 1,894
Rep Power: 26
alberto will become famous soon enoughalberto will become famous soon enough
Then I think your model is missing a term. If you transfer mass from a to b, then b has a sink of mass in its continuity equation.

Also, keep in mind that transferring mass, you transfer momentum. As a consequence you will have also a source term in the momentum equations.
__________________
Alberto Passalacqua

GeekoCFD - A free distribution based on openSUSE 64 bit with CFD tools, including OpenFOAM. Available as live DVD/USB, hard drive image and virtual image.
OpenQBMM - An open-source implementation of quadrature-based moment methods

Last edited by alberto; November 2, 2011 at 12:35. Reason: Added comment
alberto is offline   Reply With Quote

Old   November 2, 2011, 12:36
Default
  #11
Senior Member
 
Aurelien Thinat
Join Date: Jul 2010
Posts: 154
Rep Power: 7
Aurelien Thinat is on a distinguished road
I have edited the wrong post :

"I (and you) may have solved my problem. I wasn't confident in my pressure equation but it seems the problem came from the boundary condition in pressure I fixed : buoyantPressure, which I'm used for VOF models.
I have defined a zeroGradient boundary condition to the walls and the volume fraction is now constant."

Thank you again for your time.

Best regards,

Aurélien
Aurelien Thinat is offline   Reply With Quote

Old   November 2, 2011, 13:44
Default
  #12
Senior Member
 
Alberto Passalacqua
Join Date: Mar 2009
Location: Ames, Iowa, United States
Posts: 1,894
Rep Power: 26
alberto will become famous soon enoughalberto will become famous soon enough
Hi,

buoyantPressure is identical to zeroGradient if p is used instead of p_rgh, so it is a bit strange it was the source of your problems.

Best,
Alberto
__________________
Alberto Passalacqua

GeekoCFD - A free distribution based on openSUSE 64 bit with CFD tools, including OpenFOAM. Available as live DVD/USB, hard drive image and virtual image.
OpenQBMM - An open-source implementation of quadrature-based moment methods
alberto is offline   Reply With Quote

Old   June 12, 2013, 08:02
Default
  #13
Senior Member
 
Aurelien Thinat
Join Date: Jul 2010
Posts: 154
Rep Power: 7
Aurelien Thinat is on a distinguished road
Hi everyone,

Me again. But with another "source term" problem. I have to implement a turbulent dispersion force in Euler/Euler solvers. So let's start with bubbleFoam.

Starting with Lopez de Bertodano Model :
Force of turbulent dispersion, phase a = Ftd,a = -Ftd,b = C * rhob * k * grad(alpha)
with k the kinetic energy of the continuous phase.

Since the momemtum equations are divided by rhoa * alpha (rhob * beta for the second one), I have got :

Mtd,a = C * rhob * k * grad(alpha) / ( rhoa * alpha)

So I'm adding this term in the UaEqn and it's opposite in UbEqn (divided by rhob * beta instead). And...at the first pressure equation it blows up because of some rules in the operator "divide".

Am I missing something ?

Edit :

This is the UaEqn :

UaEqn =
(
(scalar(1) + Cvm*rhob*beta/rhoa)*
(
fvm::ddt(Ua)
+ fvm::div(phia,Ua,"div(phia,Ua)")
- fvm::Sp(fvc::div(phia),Ua)
)
- fvm::laplacian(nuEffa,Ua)
+fvm::Sp(fvc::div(phiRa),Ua)
+ (fvc::grad(alpha)/(fvc::average(alpha)+scalar(0.001)) & Rca)
==
-fvm::Sp(beta/rhoa*dragCoef,Ua)
- beta/rhoa*(liftCoef - Cvm*rhob*DDtUb)
-rhob*k*fvc::grad(alpha)/(alpha*rhoa) //this is the turbulent dispersion force
);

This is the error message :

#0 Foam::error:rintStack(Foam::Ostream&) in "......../libOpenFOAM.so"
#1 Foam::sigFpe::sigHandler..../libOpenFOAM.so
#2 ?? in "lib64/libc.so.6"
#3 void Foam::divide<Foam::Vector<double>,Foam::fvPatchFie ld,Foam::volMesh>
(........)
Aurelien Thinat is offline   Reply With Quote

Old   June 12, 2013, 11:16
Default
  #14
Senior Member
 
Alberto Passalacqua
Join Date: Mar 2009
Location: Ames, Iowa, United States
Posts: 1,894
Rep Power: 26
alberto will become famous soon enoughalberto will become famous soon enough
Hi Aurelien,

if your term is

Mtd,a = C * rhob * k * grad(alpha)

you have to add

C*rhob/rhoa*k*grad(alpha)/(alpha + 0.0001) to UaEqn

and

C*k*grad(alpha)/(beta + 0.0001)

since alpha and beta can be zero.

If this term is large, it's worth moving it to the flux term.

Best,
__________________
Alberto Passalacqua

GeekoCFD - A free distribution based on openSUSE 64 bit with CFD tools, including OpenFOAM. Available as live DVD/USB, hard drive image and virtual image.
OpenQBMM - An open-source implementation of quadrature-based moment methods
alberto is offline   Reply With Quote

Old   June 12, 2013, 11:32
Default
  #15
Senior Member
 
Aurelien Thinat
Join Date: Jul 2010
Posts: 154
Rep Power: 7
Aurelien Thinat is on a distinguished road
Thank you Alberto for your answer.

It's correct, alpha = 0 was the problem. For information, I'm testing right now this implementation :

UaEqn : C*rhob/rhoa*k*fvc::grad(alpha)/(fvc::average(alpha) + scalar(0.0001))
UbEqn : C*k*fvc::grad(alpha)/(fvc::average(beta) + scalar(0.0001))

Next step is the Simonin's model and working on a degassing boundary condition.

Thank you again.
Aurelien Thinat 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
swak4foam building problem GGerber OpenFOAM Installation 54 April 24, 2015 16:02
momentum source term zwdi FLUENT 13 December 5, 2013 18:58
DxFoam reader update hjasak OpenFOAM Post-Processing 69 April 24, 2008 01:24
DecomposePar links against liblamso0 with OpenMPI jens_klostermann OpenFOAM Bugs 11 June 28, 2007 17:51
UDFs for Scalar Eqn - Fluid/Solid HT Greg Perkins FLUENT 0 October 13, 2000 23:03


All times are GMT -4. The time now is 15:32.