CFD Online Logo CFD Online URL
www.cfd-online.com
[Sponsors]
Home > Forums > Software User Forums > OpenFOAM > OpenFOAM Running, Solving & CFD

strange curvature with interFoam (comparison with Brackbill work)

Register Blogs Members List Search Today's Posts Mark Forums Read

Like Tree1Likes

Reply
 
LinkBack Thread Tools Search this Thread Display Modes
Old   March 8, 2011, 10:36
Default strange curvature with interFoam (comparison with Brackbill work)
  #1
Member
 
Duong A. Hoang
Join Date: Apr 2009
Location: Delft, Netherlands
Posts: 93
Rep Power: 16
duongquaphim is on a distinguished road
Send a message via Yahoo to 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
Attached Images
File Type: jpg initial.jpg (44.7 KB, 376 views)
File Type: jpg K.jpg (37.2 KB, 381 views)
File Type: jpg p.jpg (43.6 KB, 362 views)
Attached Files
File Type: zip Brackbill_test_case.zip (12.2 KB, 158 views)
duongquaphim is offline   Reply With Quote

Old   March 9, 2011, 03:37
Default
  #2
Senior Member
 
sega's Avatar
 
Sebastian Gatzka
Join Date: Mar 2009
Location: Frankfurt, Germany
Posts: 729
Rep Power: 20
sega is on a distinguished road
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
__________________
Schrödingers wife: "What did you do to the cat? It's half dead!"
sega is offline   Reply With Quote

Old   March 9, 2011, 04:56
Default
  #3
Senior Member
 
Gijsbert Wierink
Join Date: Mar 2009
Posts: 383
Rep Power: 18
gwierink is on a distinguished road
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 ...



__________________
Regards, Gijs
gwierink is offline   Reply With Quote

Old   March 9, 2011, 05:19
Default
  #4
Senior Member
 
sega's Avatar
 
Sebastian Gatzka
Join Date: Mar 2009
Location: Frankfurt, Germany
Posts: 729
Rep Power: 20
sega is on a distinguished road
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!"
sega is offline   Reply With Quote

Old   March 9, 2011, 05:29
Default
  #5
Member
 
Duong A. Hoang
Join Date: Apr 2009
Location: Delft, Netherlands
Posts: 93
Rep Power: 16
duongquaphim is on a distinguished road
Send a message via Yahoo to 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
Attached Images
File Type: jpg alpha1_smooth.jpg (44.2 KB, 202 views)
File Type: jpg K_smooth.jpg (43.3 KB, 179 views)
duongquaphim is offline   Reply With Quote

Old   March 9, 2011, 05:39
Default
  #6
Member
 
Duong A. Hoang
Join Date: Apr 2009
Location: Delft, Netherlands
Posts: 93
Rep Power: 16
duongquaphim is on a distinguished road
Send a message via Yahoo to 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
duongquaphim is offline   Reply With Quote

Old   March 9, 2011, 06:39
Default
  #7
Senior Member
 
sega's Avatar
 
Sebastian Gatzka
Join Date: Mar 2009
Location: Frankfurt, Germany
Posts: 729
Rep Power: 20
sega is on a distinguished road
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.
__________________
Schrödingers wife: "What did you do to the cat? It's half dead!"
sega is offline   Reply With Quote

Old   March 9, 2011, 09:00
Default
  #8
Senior Member
 
Gijsbert Wierink
Join Date: Mar 2009
Posts: 383
Rep Power: 18
gwierink is on a distinguished road
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
gwierink is offline   Reply With Quote

Old   March 9, 2011, 15:20
Default
  #9
Member
 
Duong A. Hoang
Join Date: Apr 2009
Location: Delft, Netherlands
Posts: 93
Rep Power: 16
duongquaphim is on a distinguished road
Send a message via Yahoo to 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
Attached Images
File Type: jpg K.jpg (66.0 KB, 164 views)
tom_flint2012 likes this.
duongquaphim is offline   Reply With Quote

Old   March 10, 2011, 03:21
Default
  #10
Senior Member
 
Gijsbert Wierink
Join Date: Mar 2009
Posts: 383
Rep Power: 18
gwierink is on a distinguished road
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.
Quote:
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 ...
Attached Images
File Type: png Brackbill_kappa-bar_small.png (48.2 KB, 188 views)
__________________
Regards, Gijs
gwierink is offline   Reply With Quote

Old   March 10, 2011, 05:23
Default
  #11
Member
 
Duong A. Hoang
Join Date: Apr 2009
Location: Delft, Netherlands
Posts: 93
Rep Power: 16
duongquaphim is on a distinguished road
Send a message via Yahoo to 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 is offline   Reply With Quote

Old   March 10, 2011, 15:43
Default
  #12
Member
 
Duong A. Hoang
Join Date: Apr 2009
Location: Delft, Netherlands
Posts: 93
Rep Power: 16
duongquaphim is on a distinguished road
Send a message via Yahoo to 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
Attached Images
File Type: jpg K_filter.jpg (39.5 KB, 112 views)
File Type: jpg pressure.jpg (72.7 KB, 135 views)
Attached Files
File Type: zip Brackbill_square.zip (16.8 KB, 62 views)
duongquaphim is offline   Reply With Quote

Old   March 10, 2011, 16:09
Default
  #13
Senior Member
 
sega's Avatar
 
Sebastian Gatzka
Join Date: Mar 2009
Location: Frankfurt, Germany
Posts: 729
Rep Power: 20
sega is on a distinguished road
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!"
sega is offline   Reply With Quote

Old   March 10, 2011, 18:15
Default
  #14
Member
 
Duong A. Hoang
Join Date: Apr 2009
Location: Delft, Netherlands
Posts: 93
Rep Power: 16
duongquaphim is on a distinguished road
Send a message via Yahoo to 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
duongquaphim is offline   Reply With Quote

Old   March 11, 2011, 01:29
Default
  #15
Senior Member
 
Gijsbert Wierink
Join Date: Mar 2009
Posts: 383
Rep Power: 18
gwierink is on a distinguished road
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
gwierink is offline   Reply With Quote

Old   March 11, 2011, 03:39
Default
  #16
Senior Member
 
sega's Avatar
 
Sebastian Gatzka
Join Date: Mar 2009
Location: Frankfurt, Germany
Posts: 729
Rep Power: 20
sega is on a distinguished road
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.
Attached Files
File Type: zip sufaceTension.zip (6.5 KB, 119 views)
File Type: zip waterAirReal.zip (6.5 KB, 86 views)
__________________
Schrödingers wife: "What did you do to the cat? It's half dead!"
sega is offline   Reply With Quote

Old   March 11, 2011, 09:00
Default
  #17
Senior Member
 
Andrea Ferrari
Join Date: Dec 2010
Posts: 319
Rep Power: 16
Andrea_85 is on a distinguished road
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
Andrea_85 is offline   Reply With Quote

Old   March 14, 2011, 03:10
Default Smoother Implementation
  #18
ata
Senior Member
 
ata's Avatar
 
ata kamyabi
Join Date: Aug 2009
Location: Kerman
Posts: 323
Rep Power: 17
ata is on a distinguished road
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
ata is offline   Reply With Quote

Old   March 14, 2011, 04:52
Default
  #19
Member
 
Duong A. Hoang
Join Date: Apr 2009
Location: Delft, Netherlands
Posts: 93
Rep Power: 16
duongquaphim is on a distinguished road
Send a message via Yahoo to 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
duongquaphim is offline   Reply With Quote

Old   March 14, 2011, 05:20
Default
  #20
Senior Member
 
Andrea Ferrari
Join Date: Dec 2010
Posts: 319
Rep Power: 16
Andrea_85 is on a distinguished road
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
Andrea_85 is offline   Reply With Quote

Reply

Thread Tools Search this Thread
Search this Thread:

Advanced Search
Display Modes

Posting Rules
You may not post new threads
You may not post replies
You may not post attachments
You may not edit your posts

BB code is On
Smilies are On
[IMG] code is On
HTML code is Off
Trackbacks are Off
Pingbacks are On
Refbacks are On


Similar Threads
Thread Thread Starter Forum Replies Last Post
new curvature model with interFoam DaChris OpenFOAM Programming & Development 15 July 31, 2017 08:57
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 09:41
ATTENTION! Reliability problems in CFX 5.7 Joseph CFX 14 April 20, 2010 16: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 03:01


All times are GMT -4. The time now is 07:31.