In attempt to decrease spurious currents in VOF

akesm · June 15, 2017, 16:18

First I want to start with two questions:
The reasons for spurious currents in interFOAM is known to us, and we are familiar with the attempts that have been made to decrease them outside of OF community in multiphase flow community.

My first question is that why such spurious currents are not observed when using Level Set method. Because to me, all processes are similar to VOF (i.e. solving an advection equation, the way normal vector and curvature are calculated, and the way the surface tension force is used in the momentum balance equation) except that all are based on a distance function other than volume of fluid fraction. So if spurious currents are attributed to an imbalance between surface tension force and pressure force caused due to inaccuracies in curvature approximation and discretization of the surface tension force in the momentum equation, to me there is no reason that same thing does not occur for the level set method.

My second question is that, do you know of any successful attempts (code release) for recent OF versions that have eliminated the spurious currents?

There are two successful methods out there, one is balanced force algorithm and the other is through smoothing the alpha function using some kernels

Has anyone ever implemented either of these two, and is its code or results released somewhere to see how much they can improve the spurious velocity problem in interFOAM?

Thanks

floquation · June 16, 2017, 07:56

Warning: Wall of text.

Quote:

Originally Posted by akesm

First I want to start with two questions:
The reasons for spurious currents in interFOAM is known to us, and we are familiar with the attempts that have been made to decrease them outside of OF community in multiphase flow community.

My first question is that why such spurious currents are not observed when using Level Set method. Because to me, all processes are similar to VOF (i.e. solving an advection equation, the way normal vector and curvature are calculated, and the way the surface tension force is used in the momentum balance equation) except that all are based on a distance function other than volume of fluid fraction. So if spurious currents are attributed to an imbalance between surface tension force and pressure force caused due to inaccuracies in curvature approximation and discretization of the surface tension force in the momentum equation, to me there is no reason that same thing does not occur for the level set method.

I have never studied LS in detail, so I cannot give you a well-informed answer.
But what makes you say that spurious currents are not observed using the LS method?
If I perform a very quick search...:
This paper [1] is about reducing spurious currents using the LS method, so it certainly is not free of spurious currents. In addition, the (derivative of the) LS function suffers from discontinuities (e.g. exactly between two bubbles), which may give large errors in the curvature calculation [2], which is a disadvantage that VOF does not possess.

I think that the LS method does perform better (in terms of spurious currents) than VOF noneless... Again, I don't know exactly why, but I'm guessing that is because the VOF-field is like a Heaviside function, which does not have a nice gradient. LS on the other hand, has a monotonically increasing(decreasing) function, giving it a nice constant gradient in the radial direction. I'm guessing that this is very beneficial for calculating the normal and therefore also the curvature.

I did run across a Chalmers student report on coupling LS with VOF in OF, but I have never had time to really look at it.

Quote:

Originally Posted by akesm

My second question is that, do you know of any successful attempts (code release) for recent OF versions that have eliminated the spurious currents?

Has anyone ever implemented either of these two, and is its code or results released somewhere to see how much they can improve the spurious velocity problem in interFOAM?

Eliminated? No. Is it even possible to fully eliminate them without difficult and expensive interface reconstruction algorithms?

Reduced? Yes.
I have implemented the smoothing algorithm in the latest version: OF4.x. (Or actually, I ported the implementation of a former PhD student in our research group, Duong A. Hoang, from an old version to the latest version; and then reworked it such that it is user-friendly for others to use.)
More specifically, the alpha field is smoothed exclusively for the calculation of the curvature, while using the sharp alpha field in all other locations.
I have implemented this such that you only have to update one of OF's libraries. That is, your solver needs not be changed, which is great, as you can therefore use it with any solver (that uses VOF). This may be done without actually overwriting OF's library, but rather by overriding it (using a library with the same name in $FOAM_USER_LIBBIN).
I have uploaded my code to GitHub.
(I have also implemented the density-weighting proposed by Brackbill [3], and I am planning to try some more things, but I have not yet pushed these changes to GitHub.)

----

I have not published the results of this implementation anywhere, as it is not novel: it is merely implementing existing methods into OF. I have, however, performed two rising bubble simulations to evaluate its performance. I have attached these mesh-convergence figures for you. I have only looked at the magnitude of the spurious currents and the terminal rise velocity.
The first simulation is at $\mathrm{Re}\approx11.4$ with $\mathrm{Eo}=0.900$ and $\mathrm{Mo}=1.00e-5$ , which corresponds to $\mathrm{Ca}=3.8e-2$ .
Attachment [spurious_Re10_tcxpi_10.pdf] shows the magnitude of the spurious currents extracted from a zero-gravity simulation as a function of mesh resolution (number of cells per bubble diameter). "normal" is OF's implementation, and "vofsmooth" is the smoothed-curvature implementation.
Note that the spurious currents are reduced by an order of magnitude.
Also note that it is important to satisfy the capillary time step constraint! This is not enforced by OF, but it is absolutely crucial! It is given by [3]:
$\Delta t \leq \tau_c = \sqrt{\frac{\left(\rho_a+\rho_b\right)(\Delta x)^3}{4\pi\sigma}}.$
This is not new... This was even proposed by Brackbill [3] in his original paper! Yet I feel like many (especially OF users) do not know this. This also had me struggling with spurious currents for months.
That criterion is, in fact, the reason that spurious currents get worse with finer meshes: for a coarse mesh, the Courant criterion is more strict than the Capillary constraint... but as you refine your mesh, the Capillary time step constraint becomes the limiting factor. (Look at scaling with $\Delta x$ !)

Then, attachment [terminal_Re10_tcxpi_10.pdf] shows the effect on the terminal rising velocity. Note that the deviation is of the order of the magnitude of the spurious currents, resulting in a very decent mesh independence (note the scale of the y-axis).
You might wonder why the experimental result (Eotvos/Morton or Clift map) is off... This is because I was performing 2D simulations in a too small domain. Attachment [slip_Re10.pdf] shows what happens after correcting a 3D sphere to an infinite 3D cylinder (which is a circle in 2D) and using a larger domain. Then we obtain a good match.

The second simulation is for an air-water system at $\mathrm{Re}\approx160$ with $\mathrm{Eo}=0.14$ and $\mathrm{Mo}=2.86e-11$ , which corresponds to $\mathrm{Ca}=2.3e-3$ . This is a more difficult case due to its low capillary number.
Without attaching all figures comparing "normal" with "vofsmooth", etc.: the conclusion is that vofsmooth combined with the capillary timestep constraint could also drastically reduce the spurious currents and mesh convergence could be achieved. (See attachment [terminalRe160_lowerdt.pdf].)

However, I could not get these simulations to match the Eotvos-Morton/Clift map, despite them showing a nice mesh convergence when using 32 cells per diameter with $\Delta t_{\mathrm{max}} = \tau_c$ . (See attachment [terminalRe160_lowerdt.pdf].) Hypothesis: The experiments in the Clift map probably didn't use perfectly pure water. The present surfactants will increase the drag and thereby reduce the terminal rise velocity. If these experiments are therefore underpredicting the terminal rise velocity, they could still be consistent with my perfectly pure water simulations.
This is something I am looking into right now. If you have any further ideas, please have me know.

Overall, you can see that the spurious currents are reduced by two orders of magnitude when using vofsmooth and satisfying the Capillary timestep constraint.

Quote:

Originally Posted by akesm

There are two successful methods out there, one is balanced force algorithm and the other is through smoothing the alpha function using some kernels

As said, I have been using the smoothing method.
I should look into the balanced force algorithm as well... I do not know the method.
I also see that a height function is used by these authors [4], which would be nice to try. Given that OF has a general polyhedral mesh data structure rather than a nice orthogonal mesh data structure this method might be slightly difficult to implement though...

-----------------------------

[1] R. Meland, I.R. Gran, R. Olsen and S.T. Munkejord, Reduction of parasitic currents in level-set calculations with a consistent discretization of the surface-tension force for the CSF model, 16th Australasian Fluid Mechanics Conference (2007)
http://people.eng.unimelb.edu.au/ima.../16/Meland.pdf
[2] K.Y. Lervag, Calculation of interface curvature with the level-set method, Comp. Phys (2014)
[3] Brackbill, J. U., Kothe, D. B., & Zemach, C. (1992). A continuum method for modeling surface tension. Journal of computational physics, 100(2), 335-354.
researchgate-URL
arXiv:1407.7340 [physics.comp-ph]
[4] M.M. Francois et. al., A balanced-force algorithm for continuous and sharp interfacial surface tension models within a volume tracking framework, J. Comp. Phys. (2006) vol. 213
https://doi.org/10.1016/j.jcp.2005.08.004

akidess · June 16, 2017, 08:49

Fantastic post Kevin. Two small notes:

As you already have assumed, one reason LS behaves better surely is the continuous nature of the LS function and thus the more accurate computation of the gradient. You try to get the same advantage by smoothing the alpha field.

Second, the Brackbill criterion is not the only stability criterion (you are ignoring viscosity effects). Check the Deshpande paper [1] for the full story.

[1] http://dx.doi.org/10.1088/1749-4699/5/1/014016

floquation · June 16, 2017, 09:11

Quote:

Originally Posted by akidess

Fantastic post Kevin. Two small notes:

As you already have assumed, one reason LS behaves better surely is the continuous nature of the LS function and thus the more accurate computation of the gradient. You try to get the same advantage by smoothing the alpha field.

Second, the Brackbill criterion is not the only stability criterion (you are ignoring viscosity effects). Check the Deshpande paper [1] for the full story.

[1] http://dx.doi.org/10.1088/1749-4699/5/1/014016

Thank you, Anton. I'll have a look at that.

Z.Q. Niu · July 31, 2017, 08:45

Dear Anton,
Do you have good solution for reducing the parasitic current presently?

Saideep · August 1, 2017, 09:28

Hi guys;

Found this post interesting.

One thing about the LSM is the mass conservation issue. Though, there are few reinitialization concepts available for the signed distance function, this is still a concern. Hence, came out the hybrid methods which couple the VOF with LSM. Ultimately the idea is to be mass conservative (positive aspect of VOF) and better estimate the interface curvature (from the signed LS distance function).

-> Coming to few other VOF formulations, Raeinis work is a good case where they discuss about filtering the spurious currents by catching capillary forces and fluxes that potentially cause spurious currents. Their group also made their code public so you can check that (solver name: poreFoam).
https://workspace.imperial.ac.uk/ear...hys%202012.pdf

-> The other software is Gerris (http://gfs.sourceforge.net/wiki/index.php/Main_Page) where the interface is constructed using line segments and the interface is advected geometrically rather than numerically, thereby controlling numerical diffusion. Though they solve for the same capillary force equation, the power of this solver, lies in the way the interface curvature is computed using height functions (a higher order method). Gerris treats the Brackbill time step constraint on a serious note to capture the shortest amplitude of a capillary wave. There are additional test cases (apart from droplet relaxation in their website, note the developer 'Popinet' moved from Gerris to Basilisk). The larger stencil for height functions make the computation time larger (around 5 times openfoam run time). However, Gerris uses Adaptive meshing if required and used only structured grids.

HTH.

akidess · August 1, 2017, 09:43

Gerris hasn't been actively developed since 2013 I believe. The authors work on a new code called PARIS: http://www.ida.upmc.fr/~zaleski/paris/index.html

gaza · February 2, 2019, 09:23

Quote:

Originally Posted by floquation

Warning: Wall of text.
(I have also implemented the density-weighting proposed by Brackbill [3], and I am planning to try some more things, but I have not yet pushed these changes to GitHub.)

Hi floquation,
What do you mean by "density-weighting"?
Which formula of Brackbill's article is it?
In your code we have

Code:

        return fvc::interpolate(sigmaK()*rho)*fvc::snGrad(alpha1_) * 2/(rho1_+rho2_);

but I do not see such formula in the article. Can you elaborate?

floquation · February 4, 2019, 05:41

Quote:

Originally Posted by gaza

Hi floquation,
What do you mean by "density-weighting"?
Which formula of Brackbill's article is it?
In your code we have

Code:

        return fvc::interpolate(sigmaK()*rho)*fvc::snGrad(alpha1_) * 2/(rho1_+rho2_);

but I do not see such formula in the article. Can you elaborate?

The correction factor is $g(\vec{x})$ in Eq. 32 of the article I linked. The equation you copied from my code is OpenFOAM's original surface tension force, with that term added. That's in essence Eq. 33 from Brackbill's article (but using the VOF alpha-field, rather than using a density field).
The idea of adding this term made sense to me: using this term, the two fluids will have the same acceleration at the interface, which sounds to me like a desirable property. However, in practice, I obtained worse results with this term activated for my work, hence I've left it off.

gaza · February 4, 2019, 06:06

thank you for your reply
I also did not have good results with this formula
but for me this is not the same as Eq. (33)
in the implementation I see only the term
$<\rho> = \frac{\rho_1+\rho_2}{2}$
and in Eq. (33) there is also
$\frac{\rho(x)}{[\rho]}$
I do not see the last term implemented

floquation · February 5, 2019, 05:24

Quote:

Originally Posted by gaza

thank you for your reply
I also did not have good results with this formula
but for me this is not the same as Eq. (33)
in the implementation I see only the term
$<\rho> = \frac{\rho_1+\rho_2}{2}$
and in Eq. (33) there is also
$\frac{\rho(x)}{[\rho]}$
I do not see the last term implemented

The term originally implemented by OpenFOAM is essentially Eq. 28, where

is the color-function called $\alpha$ in OpenFOAM's code. Since $\alpha$ is between 0 and 1, $[\alpha]=1$ , and you have what is written in OpenFOAM's code (numerical details like interpolation aside: $\sigma\kappa\nabla\alpha$ ):

Code:

return fvc::interpolate(sigmaK())*fvc::snGrad(alpha1_);

In Brackbill's paper in Sec. III.A he says that you could use the density field as color function in Eq. 31: $c=\rho$ . In that case, the term that you are missing in the code comes in: $\frac{\nabla\rho(x)}{[\rho]}$ ... but in OpenFOAM we do not use that field: we use the original $c=\alpha$ , giving

Code:

fvc::snGrad(alpha1_)

in the code, which is "missing" in Eq. 33 of Brackbill. That is, they are the exact same terms in a different notation - in disguise.

Now, why's it the same? In OpenFOAM, $\rho=\rho_1\alpha + \rho_2\left(1-\alpha\right)$ , so $\nabla\rho=\rho_1\nabla\alpha + \rho_2\nabla\left(1-\alpha\right) = \left(\rho_1-\rho_2\right)\nabla\alpha$ .
Now, using that

represents the "jump" in the color field:

(Eq. 21), you can see that:
$\frac{\nabla\rho(x)}{[\rho]} = \frac{\left(\rho_1-\rho_2\right)\nabla\alpha}{\rho_2-\rho_1} = -\nabla\alpha = \frac{\nabla\alpha}{0-1} = \frac{\nabla\alpha}{[\alpha]}$ , where I used that $\alpha_1=1$ and $\alpha_2=0$ , which is consistent with $\rho=\rho_1\alpha + \rho_2\left(1-\alpha\right)$ .

gaza · February 5, 2019, 06:50

thank you Kevin for a detailed explanation
what a pity it does not work

tom_flint2012 · June 27, 2020, 15:47

This thread has been very helpful to me.

Kevin, would you be willing to share your changes to the setdeltat.H file to implement the timestep restriction due to the capillary time step constraint?

I plan to implement this along with the level set approach for the curvature calculation this week, and compare with your vofsmoothed approach.

EDIT: I've just realized my time-steps are way below the capillary time step constraint anyway. Thanks again for the information Kevin

All the best,
Tom

June 15, 2017, 16:18	In attempt to decrease spurious currents in VOF	#1
akesm New Member Lieh Join Date: Mar 2017 Posts: 26 Rep Power: 9	First I want to start with two questions: The reasons for spurious currents in interFOAM is known to us, and we are familiar with the attempts that have been made to decrease them outside of OF community in multiphase flow community. My first question is that why such spurious currents are not observed when using Level Set method. Because to me, all processes are similar to VOF (i.e. solving an advection equation, the way normal vector and curvature are calculated, and the way the surface tension force is used in the momentum balance equation) except that all are based on a distance function other than volume of fluid fraction. So if spurious currents are attributed to an imbalance between surface tension force and pressure force caused due to inaccuracies in curvature approximation and discretization of the surface tension force in the momentum equation, to me there is no reason that same thing does not occur for the level set method. My second question is that, do you know of any successful attempts (code release) for recent OF versions that have eliminated the spurious currents? There are two successful methods out there, one is balanced force algorithm and the other is through smoothing the alpha function using some kernels Has anyone ever implemented either of these two, and is its code or results released somewhere to see how much they can improve the spurious velocity problem in interFOAM? Thanks

June 16, 2017, 08:49		#3
akidess Senior Member Anton Kidess Join Date: May 2009 Location: Germany Posts: 1,377 Rep Power: 30	Fantastic post Kevin. Two small notes: As you already have assumed, one reason LS behaves better surely is the continuous nature of the LS function and thus the more accurate computation of the gradient. You try to get the same advantage by smoothing the alpha field. Second, the Brackbill criterion is not the only stability criterion (you are ignoring viscosity effects). Check the Deshpande paper [1] for the full story. [1] http://dx.doi.org/10.1088/1749-4699/5/1/014016 floquation likes this. __________________ On twitter @akidTwit Spend as much time formulating your questions as you expect people to spend on their answer.

August 1, 2017, 09:28		#6
Saideep Senior Member Saideep Join Date: Apr 2015 Location: INDIA Posts: 203 Rep Power: 12	Hi guys; Found this post interesting. One thing about the LSM is the mass conservation issue. Though, there are few reinitialization concepts available for the signed distance function, this is still a concern. Hence, came out the hybrid methods which couple the VOF with LSM. Ultimately the idea is to be mass conservative (positive aspect of VOF) and better estimate the interface curvature (from the signed LS distance function). -> Coming to few other VOF formulations, Raeinis work is a good case where they discuss about filtering the spurious currents by catching capillary forces and fluxes that potentially cause spurious currents. Their group also made their code public so you can check that (solver name: poreFoam). https://workspace.imperial.ac.uk/ear...hys%202012.pdf -> The other software is Gerris (http://gfs.sourceforge.net/wiki/index.php/Main_Page) where the interface is constructed using line segments and the interface is advected geometrically rather than numerically, thereby controlling numerical diffusion. Though they solve for the same capillary force equation, the power of this solver, lies in the way the interface curvature is computed using height functions (a higher order method). Gerris treats the Brackbill time step constraint on a serious note to capture the shortest amplitude of a capillary wave. There are additional test cases (apart from droplet relaxation in their website, note the developer 'Popinet' moved from Gerris to Basilisk). The larger stencil for height functions make the computation time larger (around 5 times openfoam run time). However, Gerris uses Adaptive meshing if required and used only structured grids. HTH. rasool_soofi likes this.

August 1, 2017, 09:43		#7
akidess Senior Member Anton Kidess Join Date: May 2009 Location: Germany Posts: 1,377 Rep Power: 30	Gerris hasn't been actively developed since 2013 I believe. The authors work on a new code called PARIS: http://www.ida.upmc.fr/~zaleski/paris/index.html arjun and arsalan.dryi like this. __________________ On twitter @akidTwit Spend as much time formulating your questions as you expect people to spend on their answer.

February 4, 2019, 06:06		#10
gaza Senior Member Przemek Join Date: Jun 2011 Posts: 249 Rep Power: 16	thank you for your reply I also did not have good results with this formula but for me this is not the same as Eq. (33) in the implementation I see only the term $<\rho> = \frac{\rho_1+\rho_2}{2}$ and in Eq. (33) there is also $\frac{\rho(x)}{[\rho]}$ I do not see the last term implemented __________________ best regards pblasiak

July 31, 2017, 08:45		#5
Z.Q. Niu Member Niu Join Date: Apr 2014 Posts: 55 Rep Power: 12	Dear Anton, Do you have good solution for reducing the parasitic current presently?

February 5, 2019, 06:50		#12
gaza Senior Member Przemek Join Date: Jun 2011 Posts: 249 Rep Power: 16	thank you Kevin for a detailed explanation what a pity it does not work __________________ best regards pblasiak

June 27, 2020, 15:47		#13
tom_flint2012 Member Thomas Flint Join Date: Jan 2016 Posts: 60 Rep Power: 10	This thread has been very helpful to me. Kevin, would you be willing to share your changes to the setdeltat.H file to implement the timestep restriction due to the capillary time step constraint? I plan to implement this along with the level set approach for the curvature calculation this week, and compare with your vofsmoothed approach. EDIT: I've just realized my time-steps are way below the capillary time step constraint anyway. Thanks again for the information Kevin All the best, Tom Last edited by tom_flint2012; June 28, 2020 at 07:24.

Similar Threads
Thread	Thread Starter	Forum	Replies	Last Post
whats the cause of error?	immortality	OpenFOAM Running, Solving & CFD	13	March 24, 2021 08:15
interFoam : presence of strong spurious currents in static drop in equilibrium test	swap_9068	OpenFOAM Running, Solving & CFD	8	July 17, 2018 11:08
Parasitic currents in VOF	josephmp	Fluent Multiphase	2	February 10, 2016 16:23
is internalField(U) equivalent to zeroGradient?	immortality	OpenFOAM Running, Solving & CFD	7	March 29, 2013 02:27
Spurious currents using interFoam (OF 1.6)	peterW	OpenFOAM	2	October 5, 2009 19:39