[Ralinus]

Hey!

I have a Problem with a interFoam (OpenFOAM-2.2.1 on ubuntu 13.10) simulation regarding stability. I was working at this problems the last 5 weeks, but did not find a solution. Now I don't know what I can still do. Any help or tip is very appreciated.
I tried the best I could do to pose my question the right way (http://www.cfd-online.com/Forums/ope...-get-help.html) If there is any bit of information missing, don't hesitate to ask.

My simplified Scenario is as follows:
I have a 3D cube ca. with 0.05m scale, which is filled up to approx half with a fluid with is like water, but has over 1000 times higher viscosity. (nu = 0.002, rho = 1000). The other medium is air at 20 degrees. The right wall of my box is moving downwards with approx. 2m/s, at top there is open atmosphere and all other walls are fixed walls (no-slip).

You can see a video of my result here:
http://pizzard.de/forumQuestion/CavityInstableVideo.mkv

In this scenario I experience 2 problems:
- stability issue regarding surface sharpening terms

Due to the computing-cost-expensiveness I analyzed this problems considering a much simpler (and faster computable) scenario, which I will describe now in short terms and which you can find in the http://pizzard.de/forumQuestion/aspectRatioTest.tar . (don't get carried away by the includes, I use an include pattern like described in http://olesenm.github.io/2009/11/17/...ludeIfPresent/)

I took a rectangular domain, x = 0.1, y = 0.5 , z = 0.1 and placed a rectangular fluid block at [0.02;0.08]*[0.375;0.475]*[0.02;0.08]. All boundaries are wall, except top, which is atmosphere. As simulation starts, drop starts falling and impacts surface and spreads on it.
You can see a video of this here:
http://pizzard.de/forumQuestion/AspectRatioBase.mkv
(sry its a bit broken, I messed up the fps rate, i will update this to a better one later)

The problems also appear here, the surface of the block osillates in an unphysical way along the y-axis.

I attached a picture to illustrate my first (and bigger) problem, see this picture short before Droplet impact, this cut in half and rendered with volume rendering and contour filter:


I searched the forums for a solution or for anyone having similar problems (and read most of the results).
Regarding my second problem I tried the following things:
- smaller Courant number, up to 0.01 - did not help
- other schemes, according to: http://www.cfd-online.com/Forums/ope...wing-up-2.html - did not help

Then I found this thread: http://www.cfd-online.com/Forums/ope...-behavior.html
Which gaves me the first hint that surface sharpening could be the problem. I turned cAlpha=0 and the oscillations vanish but the
solution gets too diffuse. So I know the bad guy know, but I still need him.
Second hint was the impact of the aspect ratio, which I investgated further in.
But this lead to absolutely unexplainable results for me.

I ran the scenario multiple times with different aspect ratios. The following table lists them, basically I took cubic cells and streched them in all three axis by 2,4 and approx 8. (arx4 is four times longer in x than y and z, az2 is two times longer in z than x and y, and so on.)

Due to the fact that I only changed parameters of my solver, but do solve the same physical problem all ten times, a stable and converged solver should show the same result, right?
I plotted the alpha1 value over a line, the line is showed in the picture below.

This line is placed exactly at the surface of the block that it can show the oscillations and goes to the ground that the wetting is visible. The paraview state file of the visualization can be found here: http://pizzard.de/forumQuestion/ThousandLinesGraph.pvsm
Now I plotted all 10 lines into the same graph, using a black for the base, red for x, green for y and blue for y stretching and made the lines less solid the higher there aspect ratio is.(solid = 2, dashed=4, dash-dot-dot=8).

The following video shows what happens:
http://pizzard.de/forumQuestion/Aspe...ghResVideo.mkv
Here the log of the solver: http://pizzard.de/forumQuestion/interFoamBaseLog
First, the x and z cases show the same results, which is due to the symmetry of the scenario. But most cases oscillate and yield different results from each other.
How do I get rid of this oscillations? How do I get a stable result, independent of what aspect ratio I have?

For comparison, here the same results with cAlpha=0, without oscillations:
http://pizzard.de/forumQuestion/Aspe...GraphVideo.mkv

Side note:
This also shows different drop impact times depending on aspect ratio. I asked myself, which impact time is right.
According to Newtons Law neglecting air friction I got 0.2756s. So I repeated all simulations with a near vacuum instead of air (rho=1e-6, nu=1e-6) and
I got the correct impact times for the block between 0.27 and 0.28 for the base and arxyz2 cases, but over 0.30 for the other cases.
Also the block deformation with the peek at the lower end is caused by the air flow and vanished in this case.

[Ralinus]

I reworked and improved my explanation in the question, due to the fact that it is still unanswered and I still haven't found out a solution yet.

[kwardle]

So you are using a hex mesh? In your original video from your first post it looked like tets for you original case--these will not do so good for you. Anyway, you have most likely found the correct problem in Calpha. I don't know what value you were using originally, but it shouldn't be greater than 1 or you will have these interfacial currents.
For your simplified case, I am not sure I see the instabilities in the video. The snapshot you have included makes them look very regular and very symmetric which makes me wonder why you are so certain they are unphysical in this particular case. Just because they disappear when Calpha=0 does not mean they should not exist. The Calpha=0 case is in itself unphysical since it has interface governed by numerical diffusion.
Also, I would avoid using high aspect ratio cells (>~2) with interFoam or its variants.
Just a few thoughts.
-Kent

[Ralinus]

Originally Posted by kwardle

So you are using a hex mesh? In your original video from your first post it looked like tets for you original case--these will not do so good for you.

in the Demo Cases I'm working with cubes of aspect ratio 1. This is of course not realizeable with more complex geometries.

Originally Posted by kwardle

Anyway, you have most likely found the correct problem in Calpha. I don't know what value you were using originally, but it shouldn't be greater than 1 or you will have these interfacial currents.

I also tried cAlpha=0.001 and still had these currents. I know for the formulation to be conservative, it is necessary to have cAlpha<=1.

Originally Posted by kwardle

For your simplified case, I am not sure I see the instabilities in the video. The snapshot you have included makes them look very regular and very symmetric which makes me wonder why you are so certain they are unphysical in this particular case.

They reproduce the symmetry of the given scenario, but there is one thing that really makes me sure this is not physical: I get a wide variety of results for the exact same input scenario! So 1 physically completely determined input with only different cell orientations gives 3 completely different outputs. They can't be all physical, can't they? But which is the right one?

Originally Posted by kwardle

Just because they disappear when Calpha=0 does not mean they should not exist. The Calpha=0 case is in itself unphysical since it has interface governed by numerical diffusion.

Thats true, but at least if I set this to zero and do simulation with a variety of cell sizes i get the same results if I do the same inputs and only vary cell aspect ratios.

Originally Posted by kwardle

Also, I would avoid using high aspect ratio cells (>~2) with interFoam or its variants.
Just a few thoughts.
-Kent

I also recognized that worse then 2 yields really bad results, but for complex geometries this is hard to do, sometimes up to 3.5 is just necessary.
Thanks for your feedback.

Edit: Because one picture tells you more than thousand words, in the attachment there are the results of 3 simulations with completely same input data besides cell size and orientation. The first one has cells with aspect ratio 2 because they are two times long in x direction, the second one has 2 times longer cells in y direction and the last one has the cubes mentioned above. This is a cut diagonally through the domain, time is 0.3s for all cases.

[akidess]

Your mesh might be way too coarse to get a mesh-independent solution. Only in the limit of a very fine, mesh independent solution you would expect changes to the mesh to not affect your solution.

Furthermore, did you adhere to the capillary time step criterion by Brackbill, or even better the more strict critertion by Galusinski & Vigeneaux (2008, also see Deshpande 2012)?

[Ralinus]

Thanks, I ordered the papers and will work through them and report here if they empower me to make progress in solving this problem. This will take some time.

[Ralinus]

I read through the papers and calculated the time step criteria for my above problem.
Though Latex is not supported here, i attached screenshots of my calculation results.
In fact, my simulations did not meet the criterion by brackbill and are currently rerunning with smaller timesteps. I will post the result video here as soon they finished.

I also am no longer sure if the walls around aren't too close and influence the scenario in a bad way, I will try open atmosphere boundaries at all 4 sides afterwards.

Edit:
I now have results adhering to this criteria and it turns out that they are still inconsistent and furthermore the whole problem is time-step-size-dependent, that means the solution changes its shape depending on the size of timesteps even if all timestep sizes used adhere to the brackbill criteria.

How is it possible to use variable time step sizes in a consistent manner if the solution changes depending own the time step size?
Result Video can be found here:
http://pizzard.de/forumQuestion/Aspe...bill2Video.mkv

[Ralinus]

I not further simulated the scenario with open atmosphere boundaries at all sides. and the results still exhibit the same kind of time dependence , time step dependence and oscillations. (again adhering to Brackbill and Galusinski)

Im really stuck here. What am I doing wrong, that i can't get consistent results in such a well-known and often researched scenario like a blob of liquid falling on the ground?

[akidess]

It is quite possible your mesh or time step are still too coarse. You need to continue refining until you reach asymptotic convergence. This of course assumes your setup is otherwise ok. I would focus on perfectly cubic cells with cAlpha=1, and TVD schemes such as vanLeer except for the interfaceCompression term in a first step to see if you can get a convergent solution by this way. You already made a good step by simplifying your problem and using solid walls on the boundaries. Check your pressure field for any odd results and check the residuals too see if they are reducing.

Good luck