
[Sponsors] 
September 30, 2011, 03:50 
BubbleFoam  Source Term

#1 
Senior Member
Aurelien Thinat
Join Date: Jul 2010
Posts: 154
Rep Power: 7 
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. 

September 30, 2011, 05:50 

#2 
Senior Member
Aurelien Thinat
Join Date: Jul 2010
Posts: 154
Rep Power: 7 
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. 

October 20, 2011, 10:35 

#3 
Senior Member
Aurelien Thinat
Join Date: Jul 2010
Posts: 154
Rep Power: 7 
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. 

October 21, 2011, 11:38 

#4 
Senior Member
Aurelien Thinat
Join Date: Jul 2010
Posts: 154
Rep Power: 7 
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. 

October 24, 2011, 08:15 

#5 
Senior Member
Aurelien Thinat
Join Date: Jul 2010
Posts: 154
Rep Power: 7 
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)=04) I add a source term in alphaEqn : fvm::ddt(alpha) + fvm::div(phi,alpha,scheme) + fvm::div(fvc::flux(phir,beta,scheme),alpha,scheme) ==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. 

November 2, 2011, 11:21 

#6 
Senior Member
Alberto Passalacqua
Join Date: Mar 2009
Location: Ames, Iowa, United States
Posts: 1,895
Rep Power: 26 
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, 

November 2, 2011, 11:37 

#7 
Senior Member
Aurelien Thinat
Join Date: Jul 2010
Posts: 154
Rep Power: 7 
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. 

November 2, 2011, 11:44 

#8  
Senior Member
Alberto Passalacqua
Join Date: Mar 2009
Location: Ames, Iowa, United States
Posts: 1,895
Rep Power: 26 
Quote:


November 2, 2011, 12:00 

#9 
Senior Member
Aurelien Thinat
Join Date: Jul 2010
Posts: 154
Rep Power: 7 
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). 

November 2, 2011, 12:31 

#10 
Senior Member
Alberto Passalacqua
Join Date: Mar 2009
Location: Ames, Iowa, United States
Posts: 1,895
Rep Power: 26 
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 opensource implementation of quadraturebased moment methods Last edited by alberto; November 2, 2011 at 12:35. Reason: Added comment 

November 2, 2011, 12:36 

#11 
Senior Member
Aurelien Thinat
Join Date: Jul 2010
Posts: 154
Rep Power: 7 
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 

November 2, 2011, 13:44 

#12 
Senior Member
Alberto Passalacqua
Join Date: Mar 2009
Location: Ames, Iowa, United States
Posts: 1,895
Rep Power: 26 
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 

June 12, 2013, 08:02 

#13 
Senior Member
Aurelien Thinat
Join Date: Jul 2010
Posts: 154
Rep Power: 7 
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> (........) 

June 12, 2013, 11:16 

#14 
Senior Member
Alberto Passalacqua
Join Date: Mar 2009
Location: Ames, Iowa, United States
Posts: 1,895
Rep Power: 26 
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, 

June 12, 2013, 11:32 

#15 
Senior Member
Aurelien Thinat
Join Date: Jul 2010
Posts: 154
Rep Power: 7 
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. 

Thread Tools  
Display Modes  


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