
[Sponsors] 
strange curvature with interFoam (comparison with Brackbill work) 

LinkBack  Thread Tools  Display Modes 
March 8, 2011, 10:36 
strange curvature with interFoam (comparison with Brackbill work)

#1 
Member

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 

March 9, 2011, 03:37 

#2 
Senior Member
Sebastian Gatzka
Join Date: Mar 2009
Location: Frankfurt, Germany
Posts: 729
Rep Power: 11 
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. Testcase 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 volumeoffluid methods for surface tensiondominant twophase 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
__________________
Schrödingers wife: "What did you do to the cat? It's half dead!" 

March 9, 2011, 04:56 

#3 
Senior Member
Gijsbert Wierink
Join Date: Mar 2009
Posts: 382
Rep Power: 9 
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 1e6 and the density for phase1 is 1000. For phase2 there is density 500 (as in the Brackbill paper), but shouldn't nu be 2e6 instead of 1.56e5 then? I tried it but (although the difference is small) the results seems worse ...
__________________
Regards, Gijs 

March 9, 2011, 05:19 

#4 
Senior Member
Sebastian Gatzka
Join Date: Mar 2009
Location: Frankfurt, Germany
Posts: 729
Rep Power: 11 
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.
__________________
Schrödingers wife: "What did you do to the cat? It's half dead!" 

March 9, 2011, 05:29 

#5 
Member

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 lowcapillary (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 

March 9, 2011, 05:39 

#6 
Member

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 

March 9, 2011, 06:39 

#7 
Senior Member
Sebastian Gatzka
Join Date: Mar 2009
Location: Frankfurt, Germany
Posts: 729
Rep Power: 11 
You can try to play around with the compression term in the alphaequation.
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.
__________________
Schrödingers wife: "What did you do to the cat? It's half dead!" 

March 9, 2011, 09:00 

#8 
Senior Member
Gijsbert Wierink
Join Date: Mar 2009
Posts: 382
Rep Power: 9 
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 ...
__________________
Regards, Gijs 

March 9, 2011, 15:20 

#9 
Member

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.52 times bigger (Onno Ubbink obtained similar value). I tried to run my test case with fixed deltaT = 1E6 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 (x0.03)^2+(y0.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.21.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 

March 10, 2011, 03:21 

#10  
Senior Member
Gijsbert Wierink
Join Date: Mar 2009
Posts: 382
Rep Power: 9 
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. kappabar = 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. Quote:
__________________
Regards, Gijs 

March 10, 2011, 05:23 

#11 
Member

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 

March 10, 2011, 15:43 

#12 
Member

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 

March 10, 2011, 16:09 

#13 
Senior Member
Sebastian Gatzka
Join Date: Mar 2009
Location: Frankfurt, Germany
Posts: 729
Rep Power: 11 
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?
__________________
Schrödingers wife: "What did you do to the cat? It's half dead!" 

March 10, 2011, 18:15 

#14 
Member

Hi Sebastian,
In my case, the droplet eventually move out of the center after a while (12 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 45 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 

March 11, 2011, 01:29 

#15 
Senior Member
Gijsbert Wierink
Join Date: Mar 2009
Posts: 382
Rep Power: 9 
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?
__________________
Regards, Gijs 

March 11, 2011, 03:39 

#16 
Senior Member
Sebastian Gatzka
Join Date: Mar 2009
Location: Frankfurt, Germany
Posts: 729
Rep Power: 11 
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 postprocess the pressure drop.
__________________
Schrödingers wife: "What did you do to the cat? It's half dead!" 

March 11, 2011, 09:00 

#17 
Senior Member
Andrea Ferrari
Join Date: Dec 2010
Posts: 275
Rep Power: 7 
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 

March 14, 2011, 03:10 
Smoother Implementation

#18 
Senior Member
ata kamyabi
Join Date: Aug 2009
Location: Kerman
Posts: 322
Rep Power: 8 
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 

March 14, 2011, 04:52 

#19 
Member

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 

March 14, 2011, 05:20 

#20 
Senior Member
Andrea Ferrari
Join Date: Dec 2010
Posts: 275
Rep Power: 7 
Hi Duong,
basically i set to zero alpha1 where alpha1<1e3 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 

Thread Tools  
Display Modes  


Similar Threads  
Thread  Thread Starter  Forum  Replies  Last Post 
new curvature model with interFoam  DaChris  OpenFOAM Programming & Development  14  February 6, 2015 08:19 
Curvature of the interface using interFoam  Andrea_85  OpenFOAM  1  March 2, 2011 05:19 
strange turbulent BCs work  marine  OpenFOAM  3  June 14, 2010 08:41 
ATTENTION! Reliability problems in CFX 5.7  Joseph  CFX  14  April 20, 2010 15:45 
Strange results from interFoam solution converges but sum of all forces not equal to zero  nicasch  OpenFOAM Running, Solving & CFD  0  April 15, 2008 02:01 