[henry]

There are many things you can play with in interFoam to speed it up, you could run with a higher Courant number but you may need to increase the number of gamma sub-cycles or reduce the compression effect either by reducing the cGamma or by using upwind on the compression term. You might also want to make the solver tolerances less tight but be careful that the continuity error doesn't introduce unboundedness in gamma. Depending on you case you might also find using the AMG solver on the pressure and tuning the parameters helps.

[sergei]

I have tried already a higher Courant number (above 0.35 I am considering not wise to go), and I have used compression coefficients of 1.5 and 1.0 (may be I have to try lowers ones). I have not varied much differencing schemas (I am using Gamma schema at the moment). This is because of analogy with Star, where I found that MARS has a superior performance for VOF, and I thought that Gamma and MARS are rather similar.
I will try also upwinding.
Thanks again for your suggestions, those are very helpful.

Sergei

[henry]

The performance of the code and of the interface capturing are quite sensitive to the choice of scheme for both the div terms in the gamma equation. While Gamma01 with a coefficient of 1 is a good choice for the normal convection term it is less clear if this is also appropriate for the compression term because NVD/TVD analysis is not valid in this case.

I am developing a new type of scheme for the compression term which does not rely on NVD/TVD to aid boundedness and will release a preliminary version with 1.2.

[sergei]

Hi Henry,
I have been testing interFoam with different discretisation schemes, following your advices.

1. Basically, I am currently using Gamma scheme for divergence (on gamma) with a coefficient between 0.8 and 1 (probably, for div(phi,gamma) I will use 0.5). For velocity, I am using linearUpwind. This seems to speed up the computations (in my case), especially during a number of first iterations.

2. According to your advice, I have also tried AMG for the pressure. Indeed, varying parameters speeds up the computations by some 10-15 %. I have always been thinking that AMG is more effective on large grids (I have about 100k 3D test-case, which is no t that large).

3. Results seem to be similar to STAR-CD, however, computation time is still very different (factor of about 4-5 slower). The only principal difference is that in STAR I have used an option where Courant number is restricted for the interface only. (Actually, I suspect that for interface it is kept under, let say, 0.3, and for the bulk it is probably restricted by 1.).
I have tried some variants of biasing, but not with visible success. When the Courant number is not restricted in the bulk, but only at the interface, then after a while solution does not converge (it reaches 5001 iterations, which is probably limit you built into the code, and then - diverges). This is basically, what you have been warning me about.
I tried also some gradient based biasing, but probaly made some mistakes.
Now I am thinking of something simple like maximum Courant number of 0.3 for "heavy" fluid, and 0.9, for instance, for air.
Can you advise whether the following construction is correct (I would like to multiply your definition inside the max() function in CourantNo.H by)

(1.0 - Const * neg(fvc::interpolate(gamma)))

where Const is a bias constant.
The above expression is supposed to be 1, if gamma>0, and a lower value when gamma = (or <) 0.

Can I use fvc::snGrad(gamma) for calculating gradient of the vof-scalar (if it is not zero, then there is the interface)?

Thanks in advance (it is rahter long message so)
Sergei

[henry]

1) I would recommend using Gamma with a coefficient of 1 on both terms of the gamma equation partly because in this way the scheme is TVD conformant but also because with this value it produces near-symmetric phase-fraction distributions. I have better schemes in the 1.2 release but Gamma should be finr for your purpose.

For velocity I wouldn't recommend linearUpwind, try Gamma2V 1.0. Also it would be interesting to see what the bahaviour is like if you try upwind.

What is the Reynolds number of your case?

3) is the difference between the STAR-CD and interFoam speed due to the cost per time-step or the number of time-steps required, i.e. the Courant number?

As far as case optimisation goes I can't offer much assitance unless I play with the case and get a feel for the problem. I would be happy to do this for you if you with to contact me to do so.

[sergei]

Hi Henry,
In my case the Reynolds number is not higher than 300.
It is basically multiphase flow in a microfluidic device. In a simple test it is a t-like junction of channels, which has one inlet (with a given inlet velocity for the "heavy" phase), and two pressure boundaries with different pressure values.
This means that there is simultaneously air flow from one pressure boundary to another, and water inflow from the inlet.
At the beginning typical water velocity is around 1 m/s, whereas the air velocity is some 30-35 m/s.
If I am computing time steps for typical cells, which is based on the water velocity and taking the Courant number of 0.3, it is about 1.5*10^-5 sec.
In STAR-CD I am using time step of 10^-5 sec. Typical computation lasts for 8 hours on 4 processors.
Now about the same case in FOAM. After some few iterations, the time step goes down to about 2-3*10^-7 sec.
CPU time required per time step is shorter in FOAM, but obviously not by a factor of 30.
On one CPU case is running for longer than 5 days. This would mean that FOAM solver is more than 5 times slower.
Actually, it seems to me clear that this is primerily (I hope) because of maximum global Courant number restricted by 0.3.

I have picked up one of the typical cases we have been dealing with using STAR-CD, which have also been validated expiremtally. In this way, I can see how well FOAM is performing.
Further, I am planning to look also at conjugate heat transfer and FSI (these are main areas of interest at the moment).

Thanks for your help so far. I do understand that you have quite many messages per day, and will try not to bother you much with this. If I get some progress, I will report it.

Thanks
Sergei

[henry]

Look carfully at where the velocity is causing the maximum Courant number and why. I guess there is some level of numerical instability causing spikes in the velocity field which need to be delt with. Are these spikes at boundaries? near the interface? near corners? Are they reduced or removed by using upwind? Also how good is the mesh? Try running checkMesh on it if you are unsure about the quality.

[hsieh]

Hi, Sergei,

I am wondering if you have made any progress?

I am working on a project involves a small 2D channel (3 mm wide channel converges to 1 mm wide channel with 45 degree convergences. This 1 mm wide channel is 2 mm long. Then, it diverges at 45 degree to 3 mm wide). Surface tension force is strong with 70 contact angle. Results of water flowing into the convergent section seemed OK, but, I got strange results when water/air interface flowing passed the divergent section.

Could you post the schemes you used? Thanks!

Pei

[sergei]

Hi Pei,
Basically, I had with interFoam reasonably goed resualts - my problem was (and is) that the solver is very "slow". I have tried to modify it wit relaxed Courant number for bulk flow (<1 in bulk, <0.3 at interface), and get speed up by a factor of 1.2-2 for my cases, but it becomes somemtimes unstable (especially when bubbles are formed).

Anyway, back to your question. You problem seems to be rather standard. Can you give me some more information (boundary conditions: pressure, inflow velocity, etc.), so I can try it during weekend.
At the moment (without knowing the exact boundary conditions) it looks like you have problems with mesh.

Sergei

[hsieh]

Hi, Sergei,

Thanks for the reply!

I have run the same problem with finer meshes, different div schemes - all lead to the same characteristics. Now, I am wondering whethere this is due to strong surface tension effect (smallest length scale is 1 mm with about 10 cells, that is 0.1 mm for delta X).

The liquids I am using are air and water (1000 density ratio). I am wondering how much the solution will be affected if I reduce the density raio to 100?

If you don't mind, can I email you my case file and the mp2 movies of the results? It will be highly appreciated if you can simply take a look at the schemes I used.

Pei

[ali]

I have also experienced difficulties (decreasing timestep by orders) in dealing with high density ratio cases 800:1, but when I use 40:1 I have no problem. surface tension is not high, however, maybe high shear in my problem (liq jet in air), makes it difficult to solve. The schemes I'm using:

gradSchemes
{
default Gauss linear;
grad(U) Gauss linear;
grad(gamma) Gauss linear;
}

divSchemes
{
div(rho*phi,U) Gauss upwind;
div(phi,gamma) Gauss Gamma01 0.2;
div(phirb,gamma) Gauss Gamma01 1;
}

laplacianSchemes
{
default Gauss linear corrected;
}

interpolationSchemes
{
default linear;

[hsieh]

Hi, Ali,

Thanks for the post.

I have tried out the same schemes you posted using OpenFOAM-1.1. Under OpenFOAM-1.2, I used

div(rho*phi,U) Gauss limitedLinearV;
div(phi,gamma) Gauss limitedLinear01 1;
div(phirb,gamma) Gauss gammaCompression 1;

I noticed very high vecocities near the interface. I am wondering whethere this is due to spurious currents.

Is there anything in OpenFOAM VOF that attempted to reduce spurious currents? If not, there is a paper published by Y. Renardy and M. Renardy at Virginia Tech that claims that their approach eliminates spurious currents. I am wondering if it is difficult to implement it into OpenFOAM.

Pei

[ali]

Hi Pei,

I guess for your case (surface tension dominated flow), it may be good to use filter (kernel) to smooth the volume fraction a couple of times in each time step, some improvements are observed. Some kernels are reported by Doug Kothe and Co (LANL) you may want to check and implement. They are simplest things one can do to reduce spurious currents. I believe Renardy's approach is developed for geomtery based VOF methods and many other tricks to reduce these currents are often developed for a specific case or for very simple grids. The main problem is accurate determination of interface normal and curvature where Renardy's method helps. However, it may be feasible to implement in OF.

[sergei]

Hi Pei,
Please send your case to my home address at:
s.shulepov@chello.nl
I'l try to look at it.

Basically, with other (commercial) software, which is quite similar to FOAM, we have carried out in a past very similar (to yours) simulations, and compared those with analytical solutions.
We were actually interested in interface velocity of capillary pressure driven flows, pressure distribution in both (heavy and light) fluids, etc. Results were OK.
If you are worrying about "artificial" velocities at the interface (cells) due to local curvature variations (fluctuations) - it is a different story.

Before implementing something yourself, you can try to increase then compression factor, use Gamma01 value of 1 for both cases (see advises by Henry above). I am not happy with upwind for velocities (but again, this is my experience with STAR-CD solver, and not with FOAM), albeit people of FOAM are recommending this - they know their code better.
If you want to implement your own kernel (as advised by Ali) there is a relatively simple method described in papers by Yabe and Xiao (see for instance "Description of complex and sharp interfaces...", Computers Math Applic, v.29,p15-25).
Of course references to papers by Kothe are very good, but it will cost you a bit more time, I think.
And again, if you want, I can try to modify your model and compare it with other solvers as well.

Sergei

[henry]

I am not happy with upwind for anything and do not generally recommend it except for testing and in extreme circumstances. Generally for interFoam I would recommend limitedLinearV 1 for U, limitedLinear01 1 for the normal gamma transport term and gammaCompression 1 for the compression term. If the mesh is good, the Reynolds number is low etc. etc. you should be able to use linear (central differencing) for U and if that works OK I would certainly recommend it.

[ali]

Hi Henry,
I know you're too busy, but could you please give a very brief hint about what gammaCompression scheme basically is and how did you derive it? thx.

Hi Pei,
The problem of OF with high density ratio is very much reduced in OF1.2. when I run the exactly same case I mentioned above with OF1.2, timestep remained constant and all went well for a density ratio of 350:1.

It's intersting to know what kind of magic Henry did to resolven the problem . I guess it has to do with pressure-velocity coupling and also twoPhaseMixture (maybe).

If you still have problems with OF1.2, you may want to implement the kernel, but before that make sure you play with parameters in fvSolution too, you can decrease the convergence criteria and see if it improves the results or not.

[hsieh]

Hi, Ali, Sergei, and Henry,

Thanks for the replies and suggestions.

Sergei, I have sent you the video files and case files to your home email address.

I have used both OpenFOAM-1.1 and OpenFOAM-1.2 and tried out all the schemes (div) suggested. All lead to similar flow characterisitcs. I wish I can post these video files so that you can take a look at them. I think the results look strange. There is a possibility that this could be real, but, I do not know for sure. I have to setup some experiments to test this.

Thanks again for your feedback.

Pei

[sergei]

Hi,
I have tested for a while interFam solver of version 1.2 on different meshes and different problem setups, and have to say that VOF solver performs very well. Even on rather coarse meshes results are suprisingly (at least to me) good.
It is still much slower than say STAR, but behavior of interface is much more "physical".
Henry, it is really good code, thanks.

Further, I have a question. I am actually also ineteresting in VOF with moving meshes. I have seen some discussions on interFoam and dynamic mesh, and even tried to run some examples but without success till now (because of lack of knowledge in C++, I think).
I would like to start from a simple case: squeezing flow, i.e., 2D channel, where upper wall moves down, compressing light and heavy fluid. Wall velocity may in general be a function of, let say, time. In this way, I can look at mass, momentum, etc. balance, and there are also analytical solutions for different geometries and setups (2d, axi-symmetric).

Can somebody gibve me a suggestion on how to dynamically scale the mesh during computations in a simple way (if one exests)?

Pei, may be you have already some draft version of tutorials on moving meshes? I will really appreciate it.

Again, at the moment I am interested in a simple scaling, and not in topological changes of the domain.

Thanks in advance
Sergei

PS. If this case would work (I have tried to run examples found on discussion forum, but have problems in v.1.2), and there is an interest, I can submit a simple tutorial, as a contribution into documentation.

[hsieh]

Hi, Sergei,

I wrote a draft about 2 months ago (un-complete). I had to stop because I changed to a new job. I can email you the cases Dr, Jassak sent me and the steps to run them.

But, I do not understand why you want to scale the mesh. In the situation you mentioned above, moving the wall down (or up), you only need to add mesh layers or removing mesh layers during run time.

Pei

[mattijs]

Hi Sergei,

mesh movement is done through the movingFvMesh class in the movingFvMesh library ($FOAM_SRC/movingFvMesh). The currently selected mesh gets its 'move()' function called whenever the mesh has to be moved. Have a look through the ones in the movingFvMesh/lnInclude directory to find the one that behaves the most like you want it to and copy/modify it.