[fs82]

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 two-fluid-strategy. 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 cap-shaped during its way up to the upper boundary. The horizontal boundaries are free-slip, 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 cap-shape. But then the bubble breaks into a vortex ring and a little bit later the bubble-ring 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 free-surface 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 Courant-Number 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 set-up, 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

[fs82]

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 cap-shape, getting four tentacles (might be caused by the edges), getting back its cap-shape 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?

[kwardle]

Fabian,
I will take a look at your case and see if I can offer some suggestions. One quick thing you have not mentioned--you have hit upon one subtlety of the multi-fluid--VOF 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

[kwardle]

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 (multi-fluid, 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 two-fluid model only or with GENTOP?

[fs82]

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 Bubble-Reynolds-Number and the Eötvös-Number suggest the sperical-cap shape.


Best regards
Fabian

[kwardle]

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 multi-fluid part and drag (which are not in interFoam).

[fs82]

Not yet, but this is a good idea. I will try it on Monday.

Best regards
Fabian

[Frenk_T]

Dear Fabian,
Do you have any updates regarding your problem?

Best regards
Francesco

[fs82]

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 front-capturing volume-of-fluid methods. It would be better to use, e.g., PLIC methods. However, as I am interested in Euler-Euler-Methods 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

[Andrea_85]

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

[fs82]

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 cylinder-wall 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. Level-set simulations of buoyancy-driven 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

[Andrea_85]

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.

[fs82]

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

[ssss]

Dear Fabian,

Which is your max Co? This is one of the most important parameters in compressive VOF

[fs82]

Co is less than 0.5.

Best regards
Fabian

[Andrea_85]

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

[ssss]

As stated by Andrea, use a max Co of 0.05 and you will hopefully see better results

[fs82]

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 two-fluid model in multiphaseEulerFoam. But it might be a step forward towards a solution. Thanks a lot for this help.

Fabian

[Hoyt]

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

[akidess]

Originally Posted by Hoyt

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.

Hey Nathan, this sounds like something that should be reported to OpenCFD as a bug report.