CFD Online Logo CFD Online URL
www.cfd-online.com
[Sponsors]
Home > Forums > OpenFOAM

Small time step in interFoam

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

Like Tree4Likes

Reply
 
LinkBack Thread Tools Display Modes
Old   May 17, 2011, 11:28
Default Small time step in interFoam
  #1
Senior Member
 
Andrea Ferrari
Join Date: Dec 2010
Posts: 275
Rep Power: 7
Andrea_85 is on a distinguished road
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
Andrea_85 is offline   Reply With Quote

Old   May 19, 2011, 05:35
Default
  #2
Member
 
Sabin Ceuca
Join Date: Mar 2010
Location: Munich
Posts: 42
Rep Power: 7
sabin.ceuca is on a distinguished road
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,
sabin.ceuca is offline   Reply With Quote

Old   May 19, 2011, 08:05
Default
  #3
Senior Member
 
Andrea Ferrari
Join Date: Dec 2010
Posts: 275
Rep Power: 7
Andrea_85 is on a distinguished road
Hi Sabin,
thanks for reply. I tried to implement MULES::implicit following this interFoam behavior in micro-dimensions.
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
Andrea_85 is offline   Reply With Quote

Old   May 19, 2011, 09:17
Default
  #4
Member
 
Sabin Ceuca
Join Date: Mar 2010
Location: Munich
Posts: 42
Rep Power: 7
sabin.ceuca is on a distinguished road
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,
sabin.ceuca is offline   Reply With Quote

Old   May 19, 2011, 09:28
Default
  #5
Senior Member
 
Andrea Ferrari
Join Date: Dec 2010
Posts: 275
Rep Power: 7
Andrea_85 is on a distinguished road
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
Andrea_85 is offline   Reply With Quote

Old   May 19, 2011, 09:36
Default
  #6
Member
 
Sabin Ceuca
Join Date: Mar 2010
Location: Munich
Posts: 42
Rep Power: 7
sabin.ceuca is on a distinguished road
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
sabin.ceuca is offline   Reply With Quote

Old   May 19, 2011, 11:28
Default
  #7
Senior Member
 
Andrea Ferrari
Join Date: Dec 2010
Posts: 275
Rep Power: 7
Andrea_85 is on a distinguished road
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
Andrea_85 is offline   Reply With Quote

Old   May 19, 2011, 12:27
Default
  #8
Member
 
Sabin Ceuca
Join Date: Mar 2010
Location: Munich
Posts: 42
Rep Power: 7
sabin.ceuca is on a distinguished road
Hi,
try to switch from fvSolution > PISO > momentumPredictor yes; to no.

Bye
sabin.ceuca is offline   Reply With Quote

Old   May 20, 2011, 03:47
Default
  #9
Senior Member
 
Andrea Ferrari
Join Date: Dec 2010
Posts: 275
Rep Power: 7
Andrea_85 is on a distinguished road
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
Attached Images
File Type: jpg alpha.jpg (38.3 KB, 138 views)
File Type: jpg U.jpg (40.9 KB, 133 views)
Andrea_85 is offline   Reply With Quote

Old   May 20, 2011, 04:42
Default
  #10
Member
 
Duong A. Hoang
Join Date: Apr 2009
Location: Delft, Netherlands
Posts: 92
Rep Power: 8
duongquaphim is on a distinguished road
Send a message via Yahoo to 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
duongquaphim is offline   Reply With Quote

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

Old   May 20, 2011, 08:13
Default
  #12
Senior Member
 
Pei-Ying Hsieh
Join Date: Mar 2009
Posts: 271
Rep Power: 9
phsieh2005 is on a distinguished road
Hi,

There is a thread on parasitic currents:

parasitic currents

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

Pei
phsieh2005 is offline   Reply With Quote

Old   May 20, 2011, 09:29
Default
  #13
Senior Member
 
santiagomarquezd's Avatar
 
Santiago Marquez Damian
Join Date: Aug 2009
Location: Santa Fe, Santa Fe, Argentina
Posts: 418
Rep Power: 14
santiagomarquezd will become famous soon enough
Andrea, I've been struggling with this same problem from I started to use interFoam

Velocity spikes at interfase (interFoam)

About interFoam solver #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.
__________________
Santiago MÁRQUEZ DAMIÁN, Ph.D.
Post-doctoral Fellow
Research Center for Computational Mechanics (CIMEC) - CONICET/FICH-UNL
T.E.: 54-342-4511594 Ext. 1005
Güemes 3450 - (3000) Santa Fe
Santa Fe - Argentina
http://www.cimec.org.ar
santiagomarquezd is offline   Reply With Quote

Old   May 23, 2011, 10:17
Default
  #14
Senior Member
 
Andrea Ferrari
Join Date: Dec 2010
Posts: 275
Rep Power: 7
Andrea_85 is on a distinguished road
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
Andrea_85 is offline   Reply With Quote

Old   May 23, 2011, 11:13
Default
  #15
Senior Member
 
santiagomarquezd's Avatar
 
Santiago Marquez Damian
Join Date: Aug 2009
Location: Santa Fe, Santa Fe, Argentina
Posts: 418
Rep Power: 14
santiagomarquezd will become famous soon enough
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.
__________________
Santiago MÁRQUEZ DAMIÁN, Ph.D.
Post-doctoral Fellow
Research Center for Computational Mechanics (CIMEC) - CONICET/FICH-UNL
T.E.: 54-342-4511594 Ext. 1005
Güemes 3450 - (3000) Santa Fe
Santa Fe - Argentina
http://www.cimec.org.ar
santiagomarquezd is offline   Reply With Quote

Old   May 24, 2011, 04:19
Default
  #16
Senior Member
 
Andrea Ferrari
Join Date: Dec 2010
Posts: 275
Rep Power: 7
Andrea_85 is on a distinguished road
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
Andrea_85 is offline   Reply With Quote

Old   May 24, 2011, 05:25
Default
  #17
Senior Member
 
akidess's Avatar
 
Anton Kidess
Join Date: May 2009
Location: Delft, Netherlands
Posts: 919
Rep Power: 17
akidess will become famous soon enough
Quote:
Originally Posted by Andrea_85 View Post
[...]
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
akidess is offline   Reply With Quote

Old   July 11, 2011, 08:19
Default
  #18
Senior Member
 
Andrea Ferrari
Join Date: Dec 2010
Posts: 275
Rep Power: 7
Andrea_85 is on a distinguished road
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_85 is offline   Reply With Quote

Old   November 13, 2012, 04:22
Default
  #19
New Member
 
Farhad N.
Join Date: Apr 2010
Posts: 5
Rep Power: 7
farhhad is on a distinguished road
Andrea !

I wonder if you finally made to run your model, I am dealing the sam problem how have you solved your problem ?
farhhad is offline   Reply With Quote

Old   November 13, 2012, 05:12
Default
  #20
Senior Member
 
Andrea Ferrari
Join Date: Dec 2010
Posts: 275
Rep Power: 7
Andrea_85 is on a distinguished road
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

parasitic currents

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
Andrea_85 is offline   Reply With Quote

Reply

Thread Tools
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 On
Pingbacks are On
Refbacks are On


Similar Threads
Thread Thread Starter Forum Replies Last Post
Superlinear speedup in OpenFOAM 13 msrinath80 OpenFOAM Running, Solving & CFD 18 March 3, 2015 06:36
How to write k and epsilon before the abnormal end xiuying OpenFOAM Running, Solving & CFD 8 August 27, 2013 15:33
Differences between serial and parallel runs carsten OpenFOAM Bugs 11 September 12, 2008 11:16
AMG versus ICCG msrinath80 OpenFOAM Running, Solving & CFD 2 November 7, 2006 16:15
unsteady calcs in FLUENT Sanjay Padhiar Main CFD Forum 1 March 31, 1999 12:32


All times are GMT -4. The time now is 22:30.