[duongquaphim]

Dear all,

I am working with surface tension dominant flow using interFoam. In my simulation, the calculation of surface tension force or curvature is important. I did verify interFoam with Brackbill work (Brackbill et al. "A continuum method for modeling surface tension") with the case described in page 342.

The case is: a static liquid drop with R = 2cm is placed in the middle of a tank 6x6 cm. The mesh size is 60x60. The boundary condition around is atmospheric boundary condition. The test case and initial condition are in the attachment.

I ran the case with interFoam till 0.1 to make sure that the simulation is converged. The courant number is set to 0.05 to make sure there is no instability. The obtained normalized curvature is very different from what described in Brackbill work (Figure 3). The maximum and minimum curvature obtained in OF is 30 times larger than the correct curvature. For pressure distribution, I also observed some fluctuations around the edge of the droplet. Please look at the pictures in the attachment.

I can not find any explanation why there is such a big difference between OF and that paper. If anyone have same experiences or any ideas about this, it would be great.

Regards,

Duong

[sega]

Hi Duong.

I did exactly the same calcluation in my bachelor thesis and the results were equally bad.

Furhtermore I did calculations following the works of

S. Vincent and J.-P. Caltagirone. Test-case No 10: Parasitic currents induced by surface tension (PC).

D.J.E. Harvie and M.R. Davidson and M. Rudmann. An analysis of parsitic current generation in Volume of Fluid simulations.

D. Gerlach and G. Tomar and G. Biswas and F. Durst. Comparision of volume-of-fluid methods for surface tension-dominant two-phase flows.

In every case OpenFOAM (1.5 at this time) was not able to solve this problem satisfactory. The reason for this behavious was in my option the creation of destabilizing parasitic currents. Especiall Harvie et al. give you a rather complete overview about the problem.

Sebastian

[gwierink]

Hi guys,

I'm kind of in the same boat, or bubble actually.

I downloaded your case, Duong, many thanks for that! One thing that I noticed, when checking the transportProperties and the Brackbill paper, is that the viscosity ratio seems off. The kinematic viscosity for phase1 is 1e-6 and the density for phase1 is 1000. For phase2 there is density 500 (as in the Brackbill paper), but shouldn't nu be 2e-6 instead of 1.56e-5 then? I tried it but (although the difference is small) the results seems worse ...

[sega]

Of course it gets worse.
If there is such a small viscosity there is even less damping to decrease the parasitic currents which will corrupt the solution.

[duongquaphim]

Hi Sebastian,

Thank your for quick respond. So as I understood, you face similar problem with interFoam. I found out that Unno Ubbink can obtain quite nice curvature with one of the first version of OpenFOAM (I guest it was call FOAM at that time). Then I do not know why we obtain worse curvature than before with new and better OF.

Actually, I tried different way to figured out what is the reason for this strange behavior. I tested with a smoother proposed in Ubbink's thesis (equation 3.51) but it did not improve much. I also tried to used the same discretization schemes proposed in Brackbill's paper (MAC) but it is sad that I got no improvement as well. Since I working with low-capillary (surface tension dominant flow), it is crucial for me to get good curvature. Do you have any ideas to improve it or still live with it?

Also I attached the result I got by using smooth function in Ubbink's paper if you are interested in.

Regards,

Duong

[duongquaphim]

Hi all,

I also tried to switch off the convection and viscosity term in UEqn.H and ran the case. It give incorrect result also. And Sebastian is right, a high viscosity would damp out the parasitic current then the bubble/droplet become more spherical. In Ubbink's thesis, there is a test case to confim it as well.

Duong

[sega]

You can try to play around with the compression term in the alpha-equation.
This term is intended to dampen numerical smearing of the interface and may as well have "improvement character" to parasitic currents.

Swithing to interDyMFoam looks like a good idea in the first place, but as many authors report parasitic currents tend to increase with deceasing cell size.

So you are dealing with capillary flow.
Be aware that the problem with parasitic currents is dominant in resting forms, meaning when higher velocities are infolved the influence of parasitic currents decreases.

Thats what I have found out in my bachelor thesis. Cases with moving bubbles showed no alteration due to parasitic currents but resting bubbles do.

[gwierink]

Hi Duong,

Ite helps also to switch off dynamic time stepping, especially if you want to do the Brackbill test. Also, if you want to do comparison of the method with the Lapace solution with the Brackbill case, I would not use interDyMFoam because there the mesh changes and messes up the comparison. Probably the result is better, but the quality of the comparison worse.

I tried to replicate the Brackbill case, how he compared his solution to the Laplace solution and interFoam is actually not doing so bad. In the attached graph the blue line ("p_jump") should approach the Laplace line ...

[duongquaphim]

Hi Gijs,

Indeed, the obtained pressure value is quite good. I also see that behaviour. The only thing I worry is that the calculated curvature is very different from the correct value. The curvature obtained from my simulation is 30 times bigger than the correct value while in Brackbill's work, the calculated curvature is only 1.5-2 times bigger (Onno Ubbink obtained similar value). I tried to run my test case with fixed deltaT = 1E-6 but the result after 0.1s is still the same. And I think that oscillation of curvature (or in another way the surface force) makes the parasite current bigger in OF. If we can find some way to reproduce Brackbill case, I think we can have a lot of benefit.

I did check the computation process of the curvature by constructing an analytical smooth function. I choose to compute the distance function from every cell point to the curve (x-0.03)^2+(y-0.03)^2 = R^2 like:

distance[celli] = mag(mesh.C()[celli]-(0.03,0.03,0.0025)) - R

here R = 0.02 is the radius of the drop.

The dimensionless curvature obtained from that distance function is very nice (see the figure attached). The curvature varies from 1.2-1.5K_correct. Therefore, there is no bug in curvature calculation. Then I do not know why the curvature varies that much in normal case.

One possible reason might be the shape of the drop is distorted by parasite currents during simulation. Even though that distortion is small to be recognized, it affects strongly on curvature calculation (curvature is quite sensitive to the smoothness of the volume fraction field). This is the only explanation I can think about at this moment.

Regards,

Duong

[gwierink]

Hi Duong,

Your curvature calculation looks very good! Hmm, so there must be something else wrong. Just to check, I made a field where curvature kappa is calculated from sigmaK, i.e. kappa-bar = sigmaK()/sigma/(1/0.02). A bit of a primitive trick, but maybe it gives some idea about numbers. The result after 0.1 s look like the pic attached.

The curvature obtained from my simulation is 30 times bigger than the correct value while in Brackbill's work

... seems to be the same problem here ...

[duongquaphim]

Hi Gijs,

Yes it is. And this is the problem I am addressing. I can not find any explanations for that behaviour. I thought initially that OF must give similar result as Brackbill or Ubbink but it is not the case.

I am thinking about it ....

Regards,

Duong

[duongquaphim]

Hi all,

I have quite interesting thing poped up and want to share. I tried a very similar test case with initializing a square bubble (let's call like that) instead of a circle and ran the case for a while (11s in this case) till the it reaches an equilibrium shape which is a circle. I found out that by doing that all of the calculations, curvature, pressure and surface force, is much better.

In previous pictures, I showed the curvature in the entire domain. However, if the mag(gradient(alpha1)) is small, the curvature is not important any more. Therefore, I filtered out the curvature with the criteria that mag(gradient(alpha1)) > 0.1 in this case. Attached figures show the filtered curvature and the pressure variation. It seems that in most of interfacial area, the curvature is approximately about 2K_correct. That is a good news. If I compare the pressure drop like Brackbill did, the difference <p>/p_correct = 0.9822.

So at lease, we have some good news about the case. The only thing that make difference between two cases might be the parasite currents. I guess there are tiny deformations of the droplet or alpha1 field in the previous case (initialize the circular drop) due to parasite currents and those deformations are not damped out. The curvature is very sensitive to the smoothness of alpha1 then there is error in computing curvature and pressure. In the current case, since we ran the case till equilibrium, there are much smaller deformation in alpha1 field. If that is the case, it will make sense.

If you guys have any other ideas to make this more clear, it is great to hear.

Regards,

Duong

[sega]

I did the same calculations with an initiated "square" drop which should reach its equilibrum state after a while.

With "real" values for viscosity I observed a motion of the whole drop.
http://therealsega.th.funpic.de/open...thanolReal.mpg

Did you observe this too?

[duongquaphim]

Hi Sebastian,

In my case, the droplet eventually move out of the center after a while (1-2 seconds after getting equilibrium). However, the shape of the droplet is still in equilibrium when it moves. I saw in your case, the interface is oscillating during the movement.

I tried that case with the smoother as proposed in Ubbink's thesis (equation 3.51)
alpha[celli] = sum(alpha[facei]*Sf[facei])/sum(Sf[facei])
and obtain the equilibrium without movement in 4-5 seconds. I did not ran the simulation longer so I did not know what happend latter yet.

You use ethanol which has nearly the same viscosity with fluid 1 in my simulation. How is the fluid 2? The explanation might be many things: the grid size (which affects the magnitude of the parasite currents) or the ratio between viscous force and surface tension force or the time step.

If you can provide me the operating parameters of your case, I can test it.

Regards,

Duong

[gwierink]

Hi guys,

Nice video, Sebastian. Could it be the "atmospheric" BC that causes instabilities and movement? In some of my simulations with rising bubbles a lot of air can be "pulled in" from the top, causing waves ... What kind of BCs did Brackbill use actually? Somthing like zeroGradient?

[sega]

At least I used zeroGradient.
You may have a look at my case, but I think you will learn nothing new.

I acchieved some results using aritificially increased viscosity (factor 10000) on this small mesh I see in the vide.

The "best" results were acchieved with real viscosities (water and air) and a far bigger mesh in comparison to the drop size.
There was a drifting of the drop, but at least it was in constant shape so I could post-process the pressure drop.

[Andrea_85]

Hi all,
I also work on similar problems. I'm trying to simulate flow in more or less complex geometries using interFoam. i want to calculate the curvature at the interface but i get strange results as you. The problem in my opinion lies on the definition of curvature itself. If you look at surface forces instead of the curvature you will find good results. Alpha1 will never be exactly zero or one far from the interface and when you calculate grad(alpha1)/mag(grada(alpha1) you will find curvature more or less everywhere!. At the moment the only way I found to get acceptable results is to manually apply a kind of filter on grad(alpha1) where there is no interface and then calculate the curvature, but of course the results depend quite on where i put the cut off.

any suggestion is welcome

andrea

[ata]

Hi Hoang Anh Duong
Would you please explain how you implement that smoother? Could you please release your code?
Thank you very much
Best regards

[duongquaphim]

Hi Andrea,

So you applied a filter on grad(alpha1) such as grad(alpha1) = 0 when 0<alpha1<1. Do you find much improvement on that? I tried to filter in that way but it did not work very well. Can you please post some of your result to make comparison?

Actually, I found some better results when I apply the smoother for alpha1 like the previous post. I am thinking of implementing Level Set or other methods to test.

Regards,

Duong

PS ata: send me an email to h.a.duong@tudelft.nl please

[Andrea_85]

Hi Duong,
basically i set to zero alpha1 where alpha1<1e-3 and to one where alpha1>0.999. So i can get curvature only where there is interface because the gradient far from the interface are exactly zero.
But I do not think it's a good way to solve the problem. For my particular problem I need to know the curvature of the interface and the interfacial area and seems that the results depends on where i put the cur off. (for example if I choose to put to zero alpha1 when alpha1<0.1 and to one when alpha1>0.9,the solution changes). I had no time to make some test but maybe the results could be better by increasing calpha and let the interface more compressed. I dont know if you try.

Andrea