[Andrea_85]

Hi all,
i have a question regarding the interFoam solver. Why the alpha equation is solved explicit? From alphaEqn.H we have:

 for (int aCorr=0; aCorr<nAlphaCorr; aCorr++)
00012     {
00013         surfaceScalarField phiAlpha =
00014             fvc::flux
00015             (
00016                 phi,
00017                 alpha1,
00018                 alphaScheme
00019             )
00020           + fvc::flux
00021             (
00022                 -fvc::flux(-phir, scalar(1) - alpha1, alpharScheme),
00023                 alpha1,
00024                 alpharScheme
00025             );
00026 
00027         MULES::explicitSolve(alpha1, phi, phiAlpha, 1, 0);
00028 
00029         rhoPhi = phiAlpha*(rho1 - rho2) + phi*rho2;
00030     }

I'm running simulation of water and air in a microchannel (200micron long and 50 micron wide, 2 dimensional): the time step is very small (not higher than 1e-7 s) and the simulation takes a very long time. I also found very high velocities (spurious) at the interface. I'm able to reduce these velocities by reducing the Courant number but the simulation time increase more (i'm using adjustable time step). My question is: using an implcit MULES to solve the alpha equation (for expample is used in interPhaseChangeFoam) would not be better? Normally using an implicit solver allows you to have fewer restrictions on the time step, so why in interFoam an explicit is used?

From my experience, simulations at the scale of millimeters do not give this kind of problem and the time step is well described by the CourantNumer's formula and i dont understand why decreasing the dimension should mean decrease the time step.

Thanks

andrea

[sabin.ceuca]

Ciao Andrea,
regarding the first part of your question I really don't know why MULES is solved explicitly in interfoam unlike in interPhaseChangeFoam...but you can change it to implicitly if you want and recompile the code

@ 2nd part of your question: dt=Courant*dx/U according to this small dx decreases your time step

Hope I could help you somehow,

[Andrea_85]

Hi Sabin,
thanks for reply. I tried to implement MULES::implicit following this http://www.cfd-online.com/Forums/ope...ensions-2.html.
Unfortunately the results are not better, at least for my case.
Of course if dx is small dt has to be small, but also the velocity has to be smaller. The point is that i injected the fluid from one side of my microchannel at U=1e-4 m/s and the velocities that i find at the interface are from 0.05 up 0.8 m/s!!...this is completely non-physical!! I known about the parasite currents but 4 order of magnitude higher than the inlet velocity, in my opinion, are too big. Increasing the dimension of the domain (from microns to millimiters) i also found some "not true" velocities at the interface but non more higher than 1 order of magnitude of the inlet velocity, and the results seem to be ok.
So if someone else has experience with very small channel, i'm interesting to know how to solve problem (or at least limit) of these high velocities at the interface.

Thanks

andrea

[sabin.ceuca]

Hi,
are you starting your simulations from some already converged initial values or have you just created your case? Because if you start from some "bad" initial values I would recommand you to increase your velocity at the inlet little by little, the same trick for the turbulence model (if used).
Regards,

[Andrea_85]

Hi,
Basically i start my simulation with a flat interface. On the wall i have fix contact angle, so in the first time steps the interface moves to the right position with some oscillation. Found very high velocity during this sort of relaxation i guess that is normal due to the inertial term, the problem is that the velocity remains big "until eternity" and with the correct contact angle, non only at beginning of the simulation. But of course i can try your suggestion. In that case how can i specify different velocity for different time step? How can i increase the velocity little by little?...i have fix value boundary condition at the inlet for U.

Thanks

andrea

[sabin.ceuca]

Hello,
so the easiest way is described by alberto on his home page:
http://albertopassalacqua.com/?p=69

or you could use the groovyBC:
http://openfoamwiki.net/index.php/Contrib_groovyBC

I personaly prefer for this simple case "timeVaryingUniformFixedValue" as alberto explains it

Ciao,
fammi sapere se funziona

[Andrea_85]

Ciao,
The unsteady boundary condition works fine, unfortunately even this does not solve the problem. Here part of my log file:

MULES: Solving for alpha1
Liquid phase volume fraction = 0.00786006  Min(alpha1) = -6.78795e-29  Max(alpha1) = 1
MULES: Solving for alpha1
Liquid phase volume fraction = 0.00786006  Min(alpha1) = -6.78795e-29  Max(alpha1) = 1
Using alpha1 smooth and evaluate K

DILUPBiCG:  Solving for Ux, Initial residual = 0.586111, Final residual = 5.56451e-07, No Iterations 13
DILUPBiCG:  Solving for Uy, Initial residual = 0.281552, Final residual = 8.30934e-07, No Iterations 12
DICPCG:  Solving for p_rgh, Initial residual = 0.0119703, Final residual = 0.00052313, No Iterations 7
DICPCG:  Solving for p_rgh, Initial residual = 0.000610758, Final residual = 3.03802e-05, No Iterations 116
time step continuity errors : sum local = 1.93059e-07, global = 1.37252e-08, cumulative = -5.59463e-05
DICPCG:  Solving for p_rgh, Initial residual = 0.00360174, Final residual = 0.000176314, No Iterations 18
DICPCG:  Solving for p_rgh, Initial residual = 0.000231001, Final residual = 1.15263e-05, No Iterations 161
time step continuity errors : sum local = 7.36691e-08, global = 9.963e-09, cumulative = -5.59363e-05
DICPCG:  Solving for p_rgh, Initial residual = 0.000386962, Final residual = 1.70171e-05, No Iterations 15
DICPCG:  Solving for p_rgh, Initial residual = 2.10838e-05, Final residual = 9.12595e-08, No Iterations 246
time step continuity errors : sum local = 5.816e-10, global = 3.18267e-12, cumulative = -5.59363e-05
ExecutionTime = 234.19 s  ClockTime = 234 s

Courant Number mean: 0.000148511 max: 0.174501
Interface Courant Number mean: 6.46745e-05 max: 0.174501
deltaT = 1.51997e-07
Time = 0.000245267

MULES: Solving for alpha1
Liquid phase volume fraction = 0.00786006  Min(alpha1) = -6.78794e-29  Max(alpha1) = 1
MULES: Solving for alpha1
Liquid phase volume fraction = 0.00786007  Min(alpha1) = -6.78793e-29  Max(alpha1) = 1
Using alpha1 smooth and evaluate K

DILUPBiCG:  Solving for Ux, Initial residual = 0.741821, Final residual = 5.45149e-07, No Iterations 14
DILUPBiCG:  Solving for Uy, Initial residual = 0.313413, Final residual = 6.99609e-07, No Iterations 13
DICPCG:  Solving for p_rgh, Initial residual = 0.032685, Final residual = 0.00153892, No Iterations 4
DICPCG:  Solving for p_rgh, Initial residual = 0.00160929, Final residual = 8.01163e-05, No Iterations 30
time step continuity errors : sum local = 6.38824e-07, global = -5.72751e-08, cumulative = -5.59936e-05
DICPCG:  Solving for p_rgh, Initial residual = 0.00963796, Final residual = 0.000471136, No Iterations 11
DICPCG:  Solving for p_rgh, Initial residual = 0.000542228, Final residual = 2.44036e-05, No Iterations 225
time step continuity errors : sum local = 1.9181e-07, global = 3.94815e-09, cumulative = -5.59896e-05
DICPCG:  Solving for p_rgh, Initial residual = 0.0011261, Final residual = 5.27881e-05, No Iterations 8
DICPCG:  Solving for p_rgh, Initial residual = 5.96673e-05, Final residual = 9.5319e-08, No Iterations 259
time step continuity errors : sum local = 7.3896e-10, global = -8.87931e-13, cumulative = -5.59896e-05
ExecutionTime = 234.39 s  ClockTime = 235 s

As you can see the time step is very small. Now i also used an alpha field already in the right position (with the correct contact angle between the interface and the wall), so it does not move from a flat interface to a curved interface (that maybe is the cause of these high velocities). I keep thinking about it.

Regards


andrea

[sabin.ceuca]

Hi,
try to switch from fvSolution > PISO > momentumPredictor yes; to no.

Bye

[Andrea_85]

Hi,
that is not the solution. Again the time step is of the order of 1e-7s. I leave the simulation running all the night and after 54000s of "real" time of simulation, the simulated time is 0.04s.
The velocity at the interface is again too big. See the attached picture.

Regards

andrea

[duongquaphim]

Dear Andrea,

All of those stories are about parasitic current which originally comes from the oscillation of curvature (or surface force). As my experience, reduce Courant number (till 0.1) and use smooth function for alpha1 can solve a bit.

Regards,

Duong

[Andrea_85]

Hi Duong,
I agree with you, i'm using a smooth function for alpha and my Co is 0.2. If i set Co=0.005 i'm able to reduce these velocities but the time step falls down to 1e-9/1e-10 s. Is the first time since i'm using OF that i have to work with such small domain dimension and the standard techniques to reduce spurious currents do not seem to work for this case. Problably would be interesting try to implement a new curvature calculation following other approaches (there are a lot of literature about that).

Regards

andrea

[phsieh2005]

Hi,

There is a thread on parasitic currents:

http://www.cfd-online.com/Forums/ope...-currents.html

This worked for some. Maybe it will work for you too.

Pei

[santiagomarquezd]

Andrea, I've been struggling with this same problem from I started to use interFoam

http://www.cfd-online.com/Forums/ope...interfoam.html

http://www.cfd-online.com/Forums/ope...-solver-6.html #113

I'd suggest you to run first with c_alpha=0, and div schemes like this


divSchemes
{
div(rho*phi,U) Gauss upwind;
div(phi,alpha) Gauss vanLeer;
div(phirb,alpha) Gauss interfaceCompression;
}

so, in this way you won't have influence of the nonlinear part in the advection equation and the sharpness of the free surface will be assured only by the vanLeer scheme. As you showed in #1 when c_alpha <> 0 it influences phiAlpha and then rhoPhi. Finally even when alpha is smooth, if rhoPhi is noisy it will affect in momentum equation via convective term. The other way is to retain c_alpha <> 0 and try to go up with nAlphaCorr. Increasing this parameter should improve the solution of nonlinear advection equation via a fixed point iteration in alpha, but I'm not quite sure about the improvement of this method. I did some tests in 1D phase advection and the rule in not clear for me.

Hope this could help you (and me).

Regards.

[Andrea_85]

Hi Santiago,
I'm glad to hear that I'm not the only one with this kind of problems!!

I tried to run my simulation using upwind and c_alpha=0 (and 3 nAlphaCorr). The simualtion is quite faster (delta t of the order of 1e-6 s) and also the velocities are smaller, but the interface is too diffuse (15-20 cells instead of 2-3). Using only the VanLeer scheme is not sufficient to keep the interface sharp in my opinion, but it seems that these parasite velocities come from the compression term in the alpha equation as you've suggested.

I am still pretty confused about how to solve the problem.

Regards

andrea

[santiagomarquezd]

Andrea, glad to hear your parasitic velocities are small. My suggestion was to run with c_alpha=0 and nAlphaCorr=0, since the last one is only needed in case if non-linear term is activated. The other idea is c_alpha <> 0 and nAlphaCorr > 0. In case using c_alpha > 0 is useless the only choice is to change the div(phi,alpha) discretization scheme until have the sharpness you want obtain, one I've tried is gamma (the one from Hrv Thesis).

Regards.

[Andrea_85]

Hi Santiago,
Nice to learn new things every day!

So, you said use Gamma scheme instead of Van Leer and c_alpha=0. Something like this:

divSchemes
{
div(rho*phi,U) Gauss upwind;
div(phi,alpha) Gauss Gamma01;
div(phirb,alpha) Gauss interfaceCompression;
}

PISO
{
momentumPredictor no;
nCorrectors 3;
nNonOrthogonalCorrectors 1;
nAlphaCorr 0;
nAlphaSubCycles 2;
cAlpha 0;
}


Regards

andrea

[akidess]

Originally Posted by Andrea_85

[...]
I tried to run my simulation using upwind and c_alpha=0 (and 3 nAlphaCorr). The simualtion is quite faster (delta t of the order of 1e-6 s) and also the velocities are smaller, but the interface is too diffuse (15-20 cells instead of 2-3). Using only the VanLeer scheme is not sufficient to keep the interface sharp in my opinion, but it seems that these parasite velocities come from the compression term in the alpha equation as you've suggested.
[...]

Andrea, spurious currents arise due to inaccuracies in curvature calculation. The curvature calculate is based on the gradient of the volume fraction. Steep gradients have larger discretization errors, so by turning off the interface compression you improve the curvature calculation by smoothing your alpha field. This is actually somewhat common practice (see e.g. Brackbill or Kothe on VOF), except that often a smoothed volume fraction is used for curvature calculation, and a compressed volume fraction is advected in the time stepping algorithm.

- Anton

[Andrea_85]

Hi all,
sorry for the late response, but i was testing a bit all your advice.
I tried different type of smoothing function, different schemes, different values for c_alpha (using zero is not the solution as suggested by akidess) and i also tried to add the density in the surface force calculation. For all cases i have not seen a great reduction of these spurious current (at most a factor of 1.2/1.5...).

I have also read some literature and i have 2 questions for you:

1)the current curvature calculation is implemented over a 5 points stencil or 9 points stencil? (the normal vector is calculated at the cell faces from the interpolated gradient, so i think they are 5 points...values of alpha1 in cells (i,j) (i+1,j), (i-1,j), (i,j+1), (i,j-1) for an orthogonal mesh...but please correct me if i'm wrong)

2) The second question concerns this paper:
"A balanced-force algorithm for continuous and
sharp interfacial surface tension models within a
volume tracking framework" (Francois et al. 2006, Journal of Computational Physics 213 (2006) 141–173), in which basically they say to reduce the spurious current until 10^(-18)/10^(-20) using a balance force algorithm. Does anyone of you have experience of that? would be difficult to implement in OF (hard question probably...)?


thanks

andrea

[farhhad]

Andrea !

I wonder if you finally made to run your model, I am dealing the sam problem how have you solved your problem ?

[Andrea_85]

HI Farhad,

well, actually i didn't solve the problem. Let's say that i am trying to live with it. The problem is worse if you work with small capillary numbers or very small grid cells. Smooth the color function in the curvature calculation helps a bit but it doesn't eliminate the problem. I suggest you to read this thread

http://www.cfd-online.com/Forums/ope...-currents.html

in which, based on the work done by Ali Q Raeini et al., a new version of interFoam has been proposed. At the moment i think it works just for the static case and not for dynamic. I don't know if someone has adavaced in the meantime. Maybe akidess could add his comments on that.

good luck

andrea