CFD Online Discussion Forums (https://www.cfd-online.com/Forums/)
-   OpenFOAM (https://www.cfd-online.com/Forums/openfoam/)
-   -   Small time step in interFoam (https://www.cfd-online.com/Forums/openfoam/88429-small-time-step-interfoam.html)

 Andrea_85 May 17, 2011 11:28

Small time step in interFoam

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

Code:

``` 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 May 19, 2011 05:35

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

 Andrea_85 May 19, 2011 08:05

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 May 19, 2011 09:17

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 May 19, 2011 09:28

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 May 19, 2011 09:36

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 May 19, 2011 11:28

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

Code:

```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 May 19, 2011 12:27

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

Bye

 Andrea_85 May 20, 2011 03:47

2 Attachment(s)
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 May 20, 2011 04:42

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 May 20, 2011 05:18

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 May 20, 2011 08:13

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 May 20, 2011 09:29

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.

Regards.

 Andrea_85 May 23, 2011 10:17

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 May 23, 2011 11:13

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 May 24, 2011 04:19

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 May 24, 2011 05:25

Quote:
 Originally Posted by Andrea_85 (Post 308877) [...] 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 July 11, 2011 08:19

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

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 November 13, 2012 05:12