
[Sponsors] 
June 5, 2015, 07:33 
A rising bubble in a stagnant liquid with twofluidmodel

#1 
Senior Member
Dr. Fabian Schlegel
Join Date: Apr 2009
Location: Dresden, Germany
Posts: 222
Rep Power: 11 
Dear Foamers,
after a couple of month without CFD I started again, but with a new topic: multiphase flows. My aim is to include the so called GENTOP concept into OpenFoam (Haensch et al. 2012, International Journal of Multiphase Flow, 47, pp. 171  182). This concept is intended to simulate a dispersed and segregated fluid simultaneously by using a twofluidstrategy. Something similar is done by Kent Wardle in Wardle et al. 2013, International Journal of Chemical Engineering. As a mechanism for interface capturing is cruical for this strategy, multiphaseEulerFoam is the only appropriate solver. As a starting point I decided to deal with a single bubble in a stagnant liquid water column. The bubble has a diameter of about 9mm and the dimensions of the solution domain are 50 x 50 x 100mm. The mesh resolution is 1mm. The bubble is initialized as a sphere and is expected to become capshaped during its way up to the upper boundary. The horizontal boundaries are freeslip, the bottom a fixed wall and the upper boundary an outlet. The bad thing is, that in OpenFoam the bubble starts rising and forming the capshape. But then the bubble breaks into a vortex ring and a little bit later the bubblering explodes and disappears. As a note, CFX is capable to simulate this problem and gives reasonable results. I have to mention that CFX was configured to use a so called freesurface drag formulation: F = 1/2 * rho * cd * A * U_slip**2, where rho is the liquid density, cd a drag coefficient of 0.44, A the interfacial area density A = grad(alpha.liquid) and U_slip the slip velocity. I implemented a new drag model using this formulation. For a simple test case I get reasonable values for K (function that is used to calculate the drag), but I am not able to figure out how the drag is finally implemented or used in pEqn.H. This is very complex code and hard to understand. My problem is, that the bubble "explodes". I tested several configurations: interface compression on/off, surface tension on/off, time discretization Euler/backward, various limiter functions for the convective term in the momentum equations (linear, limitedLinearV, vanLeerV, gammaV), Schiller/Nauman drag model, higher mesh resolution, alpha subcycling (2,4). Additionally, I tried several simplified test cases, e.g. without gravitation, where everything is fine (bubble keeps its spherical shape and stays in place). However with gravitation, the only way to get a rising bubble without "explosion" is to increase the surface tension by a factor of 7 (sigma = 0.5 instead of 0.07). The CourantNumber is well below 0.5, either for explizit or implizit time discretization. Now I try to understand why this happens. Did anybody tried such a case before, or does anybody has an idea how to stabilize the bubble? I am running out of new ideas at the moment. I attached an example case (forum.zip). If there are further questions concerning the setup, please ask. I am happy for every comment. An important questions for me is, why the backward scheme does not work. With implicit time discetization, the bubble does not break into a vortex ring, it explodes directly. kind regards, Fabian Last edited by fs82; June 5, 2015 at 08:48. Reason: Correction of sigma value for running case 

June 5, 2015, 09:26 

#2 
Senior Member
Dr. Fabian Schlegel
Join Date: Apr 2009
Location: Dresden, Germany
Posts: 222
Rep Power: 11 
Allright, I continued testing (PCG/GMAG, various values for sigma, interface compression 0.5). I have a feeling that the breakage into this ring structure is caused by the interface compression. If I turn it off, the bubble starts rising, forming the capshape, getting four tentacles (might be caused by the edges), getting back its capshape and than exploding into small bubbles. But still I have no clue how to fix this behaviour.
kind regards, Fabian P.S: Is it possible to figure out where and how the interface compression acts on my interface? 

June 11, 2015, 14:02 

#3 
Senior Member
Kent Wardle
Join Date: Mar 2009
Location: Illinois, USA
Posts: 218
Rep Power: 14 
Fabian,
I will take a look at your case and see if I can offer some suggestions. One quick thing you have not mentionedyou have hit upon one subtlety of the multifluidVOF coupling and that is that the drag model for the multifluid part is active everywhere, even in interface sharpened regions. At least this is true of the standard multiphaseEulerFoam. You will find that the results are also dependent on what you specify for the dispersed phase diameter of your gas phase and how you specify the drag model pairing in transportProperties (which should be in the order (dispersedPhase continuousPhase) ). You might try to specify a very large diameter for the air bubble phase. Even in the case where you apply interface sharpening everywhere you will see an impact. For a vertical plunging liquid jet case I ran recently, we found that even in the case of interface sharpening everywhere (using C_alpha = 1), the bubble size specified for the gas phase made a big difference on the entrainment behavior of the jet. This is obviously a problem. You are correct that the drag implementation in pEqns.H is a little messy, but that is because it was meant to be generic so that you could specify different drag models for different phase pairs and also have different blending schemes if you want. Kent 

June 11, 2015, 16:15 

#4 
Senior Member
Kent Wardle
Join Date: Mar 2009
Location: Illinois, USA
Posts: 218
Rep Power: 14 
I have played with your case a little bit now. With the hybrid scheme, keep in mind that you inherit the capabilities of each method (multifluid, VOF) but you also get the drawbacks. If you are trying to resolve the interface curvature with interfaceCompression then you need enough mesh to do so or your bubble will breakup with a size on the order of the mesh. I ran it with as much as 4x the mesh in the central region and it still breaks in the middle. Have you tried this same case with the same mesh using interFoam? When you say CFX can do this, you are meaning for twofluid model only or with GENTOP?


June 12, 2015, 02:32 

#5 
Senior Member
Dr. Fabian Schlegel
Join Date: Apr 2009
Location: Dresden, Germany
Posts: 222
Rep Power: 11 
Dear Kent,
thank you for your reply and for your help. Currently I am playing around with the interface compression as I have the same feeling, that it causes the bubble breakup. Unfortunately I am not able to get a good understanding of the compression mechanism only from the equations, and, hence, I started to play with it in 1D and with finite differences. But this is ongoing work and just a sidemark. CFX was not configured to use GENTOP. It uses the free surface drag formulation and its compressive scheme, which should be something similar (from the numerical point of view) as the interface compression used in OF. Concerning the drag I use my own dragmodel, which is hopefully attached to my first post. The drag formulation does not use a bubble diameter any more and is based on the gradient of the phase fraction. I borrowed this from the CFX manual. Hence, it should not generate drag outside the interface region. This is a quick image of the drag in one of my test simulations. I also tried a finer grid as well with a resolution of about half a millimetre in the whole domain. The behaviour is almost the same, the bubble breaks up in the middle. I plan to do, after finishing playing around with the compression mechanism, to try different bubble sizes. May be the bubble has a "wrong" diameter and is in a "bad" regime. However, the BubbleReynoldsNumber and the EötvösNumber suggest the spericalcap shape. Best regards Fabian 

June 12, 2015, 09:14 

#6 
Senior Member
Kent Wardle
Join Date: Mar 2009
Location: Illinois, USA
Posts: 218
Rep Power: 14 
Have you tried it with interFoam? It would be good to narrow down whether it is simply the challenge to resolve the thin interface in the center of the cap or whether it has something to do with the multifluid part and drag (which are not in interFoam).


June 12, 2015, 09:23 

#7 
Senior Member
Dr. Fabian Schlegel
Join Date: Apr 2009
Location: Dresden, Germany
Posts: 222
Rep Power: 11 
Not yet, but this is a good idea. I will try it on Monday.
Best regards Fabian 

October 25, 2015, 12:03 

#8 
New Member
Join Date: Sep 2015
Location: Germany
Posts: 11
Rep Power: 4 
Dear Fabian,
Do you have any updates regarding your problem? Best regards Francesco 

October 26, 2015, 05:20 

#9 
Senior Member
Dr. Fabian Schlegel
Join Date: Apr 2009
Location: Dresden, Germany
Posts: 222
Rep Power: 11 
I do not have a solution yet, but I am still working on that issue.
It took me some time to define a test case for a rising bubble with interFoam (boundary conditions are a petty) and I am still a bit disappointed about the results (errors in rising velocity between 10 and 30 percent). However, the results are reproducable and seems to be a problem of such frontcapturing volumeoffluid methods. It would be better to use, e.g., PLIC methods. However, as I am interested in EulerEulerMethods I decided to accept the errors (for reference it should be good enough and most important thing is: the bubble does not explode) and currently I am still trying to figure out why the bubble explodes in multiphaseEulerFoam. I do have an idea but I am still testing and implementing various possible solutions and may be I have to adopt them somehow. That is all I am able to say at the moment. kind regards Fabian 

October 26, 2015, 09:18 

#10 
Senior Member
Andrea Ferrari
Join Date: Dec 2010
Posts: 313
Rep Power: 9 
Hi All,
just a comment on interfoam. I recently simulated with interfoam bubble rising for various combination of the dimensionless numbers (Mo, Bo, Ca and Re) and got reasonably good results (error with experimental data almost always below 5%). Simulations were 2D axysimmetric with domain size 4D width and 12D lenght (D original bubble diameter). I used fixed pressure bc at top and bottom (pressureInletOutlet for velocity, zeroGradient for alpha) and slip condition on the side. i compared the results in terms of rising velocity and bubble shape with those from this pubblication "Numerical simulation of bubble rising in viscous liquid", by Jinsong Hua and Jing Lou, 2006 I did also some 2D simulations and compared the results with "Numerical simulation of a single rising bubble by VOF with surface compression" (http://onlinelibrary.wiley.com/doi/1....3692/abstract) in which the Authors also used openFoam. Results are particulary dependent on time step size and mesh size and it took me some time to get acceptable errors. For all the cases where surface tension were not that important i got good results compared with experiments. On the other hand if surface tension becomes important the results get worse due to large spurios currents on the bubble surface (that also affect the rising velocity). I do not know exactly which combination of dimensionless numbers you are using but since your bubble is expected to deform during its way up to the boundary i would argue that surface tension is smaller compared to the other forces and you should be able to get reasonably good results (at least using interfoam). Hope this help. Andrea 

October 26, 2015, 09:36 

#11 
Senior Member
Dr. Fabian Schlegel
Join Date: Apr 2009
Location: Dresden, Germany
Posts: 222
Rep Power: 11 
This is very interesting, I have errors between 10% to 30% in terms of the rising velocity. I use a cylindrical domain with a diameter of 8D and a height of 12D (D equals the bubble diameter). The mesh resolution was 16 and 32 cells per bubble diameter. The funny thing is that the results on a finer mesh are always worst than the ones on the coarser grid. My boundary conditions are a fixed wall at the bottom (fixedValue for velocity, fixedFluxPressure for p_rgh and fixedValue for alpha), slip at the cylinderwall and a air buffer at the outlet (pressureInletOutlet for velocity, totalPressure for p_rgh and inletOutlet for alpha). The air buffer was added because I experienced that the bubble was destroyed while approaching the outlet.
I used a density ratio of 100, Morton numbers of 41, 1 and 0.001, and Eötvös numbers of 116, 7 and 10, respectively. My reference is Balcázar, N.; Lehmkuhl, O.; Jofre, L. & Oliva, A. Levelset simulations of buoyancydriven motion of single and multiple bubbles International Journal of Heat and Fluid Flow, 2015, 56, 91  107 I guess that the mesh resolution is very important. Do you use the "default" fvSchemes form the tutorial cases or did you changed something? Best regards Fabian 

October 26, 2015, 10:12 

#12 
Senior Member
Andrea Ferrari
Join Date: Dec 2010
Posts: 313
Rep Power: 9 
I took fvSchemes and fvSolution from the damBreak tutorial (i have just changed div(phirb,alpha) to Gauss vanLeer, but no big difference with respect to interfaceCompression). Momentum predictor on (important in my case for the finest grid) and Co = 0.05.
For some cases i also got "worst" results on a finer grid than the ones on the coaser grid, in the sense that i am not converging to the same value as in the experiments but to a slightly different value (but still with relative error less than 5% with respect to the experiment). I think that match exactly the condition of the experiment is not easy and some difference are somehow expected. I have the same simulations in fluent with the same BC and interFoam simulation always converge to fluent results. 

October 26, 2015, 10:31 

#13 
Senior Member
Dr. Fabian Schlegel
Join Date: Apr 2009
Location: Dresden, Germany
Posts: 222
Rep Power: 11 
OK, I did not use the momentum predictor so far. May be I should give it a try. Do you have a test case for me? I just want to try it by myself.
Best regards Fabian 

October 26, 2015, 11:34 

#14 
Senior Member
anonymous
Join Date: Aug 2014
Posts: 204
Rep Power: 6 
Dear Fabian,
Which is your max Co? This is one of the most important parameters in compressive VOF 

October 26, 2015, 11:38 

#15 
Senior Member
Dr. Fabian Schlegel
Join Date: Apr 2009
Location: Dresden, Germany
Posts: 222
Rep Power: 11 
Co is less than 0.5.
Best regards Fabian 

October 26, 2015, 11:57 

#16 
Senior Member
Andrea Ferrari
Join Date: Dec 2010
Posts: 313
Rep Power: 9 
Here is a test case
M = 848; E0 = 116, reference solution from the paper in my previos post in terms of reynolds number is Re = 2.47. Re calculated as Re = rho_liquid*U_inf*D/mu_liquid U_inf rising velocity. run blockmesh. Then initialize the bubble centered at (0, 2D) with radius D/2 (D=3mm) using funkysetfield. run the simulation without gravity for some time to let the bubble relax. Then use the last saved time as initial condition for the bubble rising simulation. I got Re around 2.1, less than 3% error. ps C0 = 0.5 is too high, i would set it at least 1 order of magnitude smaller. Andrea edit sorry, i got Re = 2.41 

October 26, 2015, 12:13 

#17 
Senior Member
anonymous
Join Date: Aug 2014
Posts: 204
Rep Power: 6 
As stated by Andrea, use a max Co of 0.05 and you will hopefully see better results


December 3, 2015, 12:29 

#18 
Senior Member
Dr. Fabian Schlegel
Join Date: Apr 2009
Location: Dresden, Germany
Posts: 222
Rep Power: 11 
Allright, a Co number of less than 0.05 reduces my error for Morton number 41 and Eötvös number 116 on a coarser grid (16 cells per diameter) down to 2.75%. What a shame that I do not noticed that, but I allways thought the limit is 1 for explicit time scheme. Nevertheless, after a quick test this do not solve my problem with the twofluid model in multiphaseEulerFoam. But it might be a step forward towards a solution. Thanks a lot for this help.
Fabian 

December 3, 2015, 14:09 

#19 
New Member
Join Date: Dec 2010
Posts: 4
Rep Power: 8 
Hi Fabian,
Along with Kent Wardle I’ve been using multiphaseEulerFoam to simulate some rising bubble cases and we have noticed the same problems as you (exploding bubbles, invalid reproduction of the Laplace pressure, etc.). We’ve been able to address the issue with two additions to the code. First, the surface tension force needs to be added to the fvc::reconstruct() at the end of the pEqn.H file. Surface tension is included earlier in the construction of the pressure equation itself, but for some reason it has been omitted from the reconstruction of the velocity field. Even with the addition of this term though, the rising bubble cases still do not perform as expected and the interface is wildly distorted. We have found that the addition of interfacial curvature smoothing (as in Raeini, et al. Journal of Computational Physics 231 (2012) 5653–5668) rectifies this issue. With these changes we’ve been able to run rising bubble cases with Co=0.5 that behave as expected. Best regards, Nathan 

December 4, 2015, 03:02 

#20  
Senior Member
Anton Kidess
Join Date: May 2009
Location: Germany
Posts: 1,262
Rep Power: 23 
Quote:
__________________
*On twitter @akidTwit *Spend as much time formulating your questions as you expect people to spend on their answer. 

Thread Tools  
Display Modes  


Similar Threads  
Thread  Thread Starter  Forum  Replies  Last Post 
bubble rising problem  swamysrikanth  ANSYS Meshing & Geometry  3  May 31, 2016 11:09 
Rising bubble with interFoam  tayo  OpenFOAM Running, Solving & CFD  15  December 5, 2014 04:49 
Simulation of a single bubble with a VOFmethod  Suzzn  CFX  18  October 2, 2009 04:18 
mass flow in is not equal to mass flow out  saii  CFX  2  September 18, 2009 08:07 
Terrible Mistake In Fluid Dynamics History  Abhi  Main CFD Forum  12  July 8, 2002 09:11 