[fs82]

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.

Yes, I could confirm this error. I discovered it a few days ago, and I put it on my list to submitt this bug to OpenFoam. However, I did not have time yet. For me it seems to be a mistake, in the new solver reactingMultiphaseEulerFoam in release 3.0.x it was added to the reconstruction of the velocity.

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.

This is very interesting, as this is the way to a possible solution I figured out as well. I tried a simple height function method for a better representation of the surface curvature, but I did not solve my problem. The idea of a smoothing I discovered a few days ago and gave it a quick try. However, It did not work so far, but after quickly going through the paper you mentioned I made a mistake, as I smoothed the phase fraction before calculating the surface normals and the curvature. I will give it a try again. Thanks for your reply. I will report as soon as I have new results.

Fabian

[fs82]

Dear Hoyt

is it possible to share a working set-up? Or do you use my set-up attached to my first post?

kind regards,
Fabian

[Hoyt]

Hi Fabian,

I've been using the 2-D bubble rise benchmark scenario that can be found here: http://www.featflow.de/en/benchmarks...ng/bubble.html

The case entails comparison of the interface shape after allowing an initially spherical bubble to rise for a fixed amount of time. I've attached a plot to this post showing the interface shape that is predicted by interFoam, by the original multiphaseEulerFoam, and by the fixed version of multiphaseEulerFoam.

Give me a day or two and I'll post the case and the updated solver.

I agree with Anton that a bug report should definitely be made at some point.

Best regards,
Nathan

bubbleRiseInterface.png

[Hoyt]

See the attached files for the bubble rise case and an updated version of multiphaseEulerFoam.

In updating the solver to OpenFoam 3.0.0, I've realized that the interfacial smoothing that I talked about earlier isn't quite as critical as I'd remembered. The drag between the phases is actually what seems to cause the bubbles to “explode” and distort non-physically. Along with smoothing and the fixed velocity reconstruction, the version I've attached also multiplies the interphase drag for a given phase pair by the respective value of (1-cAlpha). This seems to make the rising bubble case and other high capillary number cases work fine, but additional investigation probably needs to be done into this issue.

It is possible that the free-surface drag formulation that Fabian was attempting may be a better way to go. Alternatively, cAlpha switching (as in Shonibare, O., Wardle, K., “Numerical Investigation of Vertical Plunging Jet Using a Hybrid Multifluid–VOF Multiphase CFD Solver” International Journal of Chemical Engineering, Volume 2015 (2015)) is another approach that could also work.

[fs82]

Thanks a lot for sharing the set-up and the solver. I will give it a try and post my results.

The problem is to my opinion that the rising or more simplified the stagnant (no gravitation) bubble case is a limiting case for the two-fluid model. As alpha is always multiplied in all terms of the momentum equation one will have the result 0=0 except in the vicinity of the interface. To avoid this the drag is set to a small but always present value (residualPhaseFraction and residualSlip). I have the feeling that this destabalizes the numerics, and, henceforth, the bubble due to the small velocity which is always present. Unfortunatly, I do not have a proper solution yet.

kind regards
Fabian

[fs82]

Allright, after some quick and dirty tests, it seems to work with OF 3.0.x. However, with OF 2.4.x it still do not work. I tried Nathans solver with the his bubble test case and this works fine. Nathans solver with my bubble test case (scheme and solution settings taken from Nathans case) works well, but gives wrong results in terms of bubble shape, but this might be due to the simple drag model. Now, backporting Nathans changes (nHatfvSmooth and K function) to my solver in OF 2.4.x (with fixed surface tension) does not work. As in 2.4.x the p_rgh pressure does not exist I tried to implement it as well. First I thought this should be easy, but now after spending some hours on it and it still not work, it seems to be more challenging

kind regards
Fabian

Edit: I found the problem, but currently I do not understand the reason for it. In Nathans solver is another modification in the dragCoeffs function. The line of interest is

max
(
  dm.phase1()*dm.phase2()*(1-cAlpha()),
  dm.residualPhaseFraction()
)

the term (1-cAlpha()) seems to be important to stabilize the bubble. If I understand it right, it reduces the drag as 1-cAlpha should be zero, and, hence, the result of the K function is always multiplied by residualPhaseFraction().

[Hoyt]

Hi Fabian,

Sorry for the confusion. As i mentioned in my earlier message, the drag term is more important than I had thought in my original post. In the solver that I posted the (1-cAlpha) term does serve to reduce the drag for the phase pair with interface compression turned on.

This is perhaps too drastic a step for cases where the phases might be dispersed in some regions of the simulation; for those cases treating cAlpha as a field variable and setting it to 1 or 0 locally depending on the ability of the solver to resolve interfacial features is an approach that I can confirm does work for more general problems.

Best regards,
Nathan

[Frenk_T]

Hi Nathan,

Thanks a lot for sharing the solver. I have one question: treating cAlpha as a field variable is not possible in the release version of the solver?

Best regards,
Francesco

[akidess]

Originally Posted by Frenk_T

Hi Nathan,

Thanks a lot for sharing the solver. I have one question: treating cAlpha as a field variable is not possible in the release version of the solver?

Best regards,
Francesco

It's not (it's defined as a scalar), but the code change is trivial: define a volScalarField IOobject cAlpha, read it from a file.

[ddovizio]

Originally Posted by Frenk_T

Hi Nathan,

Thanks a lot for sharing the solver. I have one question: treating cAlpha as a field variable is not possible in the release version of the solver?

Best regards,
Francesco

Hi Francesco,
I am interested in treating cAlpha as a non constant value. Did you try that as well?
Thanks,

Daniele

[kwardle]

Anton, you are mostly correct in saying that the changes needed are straightforward. One challenging point is that cAlpha is not used as a volScalarField in most places--it is actually a surfaceScalarField--though there are places that it needs to be used as a volScalarField and so some switching back and forth through interpolation/averaging is required. Additionally, cAlpha is defined for each phase pair so it needs to be constructed as a PtrDictionary lookup similar to some other properties in transportProperties. This perhaps makes it less trivial but still straightforward.

I am hopeful that we can open-source release our modified solver in some form at some point. Nathan will continue working on this front.