# BubbleFoam - Source Term

 User Name Remember Me Password
 Register Blogs Members List Search Today's Posts Mark Forums Read

 September 30, 2011, 03:50 BubbleFoam - Source Term #1 Senior Member   Aurelien Thinat Join Date: Jul 2010 Posts: 165 Rep Power: 9 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 (zz0). 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: 165 Rep Power: 9 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: 165 Rep Power: 9 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: 165 Rep Power: 9 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&) [with Type = double]() + Foam::fvc::interpolate(const Foam::GeometricField&) [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: 165 Rep Power: 9 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 isfvm::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.

 November 2, 2011, 11:21 #6 Senior Member   Alberto Passalacqua Join Date: Mar 2009 Location: Ames, Iowa, United States Posts: 1,911 Rep Power: 28 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 in both physical and virtual formats (current status: http://albertopassalacqua.com/?p=1541) OpenQBMM - An open-source implementation of quadrature-based moment methods. To obtain more accurate answers, please specify the version of OpenFOAM you are using.

 November 2, 2011, 11:37 #7 Senior Member   Aurelien Thinat Join Date: Jul 2010 Posts: 165 Rep Power: 9 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,911
Rep Power: 28
Quote:
 Originally Posted by Aurelien Thinat 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 in both physical and virtual formats (current status: http://albertopassalacqua.com/?p=1541)
OpenQBMM - An open-source implementation of quadrature-based moment methods.

To obtain more accurate answers, please specify the version of OpenFOAM you are using.

 November 2, 2011, 12:00 #9 Senior Member   Aurelien Thinat Join Date: Jul 2010 Posts: 165 Rep Power: 9 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,911 Rep Power: 28 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 in both physical and virtual formats (current status: http://albertopassalacqua.com/?p=1541) OpenQBMM - An open-source implementation of quadrature-based moment methods. To obtain more accurate answers, please specify the version of OpenFOAM you are using. 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: 165 Rep Power: 9 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,911 Rep Power: 28 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 in both physical and virtual formats (current status: http://albertopassalacqua.com/?p=1541) OpenQBMM - An open-source implementation of quadrature-based moment methods. To obtain more accurate answers, please specify the version of OpenFOAM you are using.

 June 12, 2013, 08:02 #13 Senior Member   Aurelien Thinat Join Date: Jul 2010 Posts: 165 Rep Power: 9 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::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,911 Rep Power: 28 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 in both physical and virtual formats (current status: http://albertopassalacqua.com/?p=1541) OpenQBMM - An open-source implementation of quadrature-based moment methods. To obtain more accurate answers, please specify the version of OpenFOAM you are using.

 June 12, 2013, 11:32 #15 Senior Member   Aurelien Thinat Join Date: Jul 2010 Posts: 165 Rep Power: 9 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 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 zwdi FLUENT 14 June 27, 2017 15:40 [swak4Foam] swak4foam building problem GGerber OpenFOAM Community Contributions 54 April 24, 2015 16:02 hjasak OpenFOAM Post-Processing 69 April 24, 2008 01:24 jens_klostermann OpenFOAM Bugs 11 June 28, 2007 17:51 Greg Perkins FLUENT 0 October 13, 2000 23:03

All times are GMT -4. The time now is 18:17.

 Contact Us - CFD Online - Privacy Statement - Top