
[Sponsors] 
May 17, 2011, 11:28 
Small time step in interFoam

#1 
Senior Member
Andrea Ferrari
Join Date: Dec 2010
Posts: 296
Rep Power: 8 
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 1e7 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 

May 19, 2011, 05:35 

#2 
Member
Sabin Ceuca
Join Date: Mar 2010
Location: Munich
Posts: 42
Rep Power: 8 
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, 

May 19, 2011, 08:05 

#3 
Senior Member
Andrea Ferrari
Join Date: Dec 2010
Posts: 296
Rep Power: 8 
Hi Sabin,
thanks for reply. I tried to implement MULES::implicit following this interFoam behavior in microdimensions. 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=1e4 m/s and the velocities that i find at the interface are from 0.05 up 0.8 m/s!!...this is completely nonphysical!! 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 

May 19, 2011, 09:17 

#4 
Member
Sabin Ceuca
Join Date: Mar 2010
Location: Munich
Posts: 42
Rep Power: 8 
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, 

May 19, 2011, 09:28 

#5 
Senior Member
Andrea Ferrari
Join Date: Dec 2010
Posts: 296
Rep Power: 8 
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 

May 19, 2011, 09:36 

#6 
Member
Sabin Ceuca
Join Date: Mar 2010
Location: Munich
Posts: 42
Rep Power: 8 
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 

May 19, 2011, 11:28 

#7 
Senior Member
Andrea Ferrari
Join Date: Dec 2010
Posts: 296
Rep Power: 8 
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.78795e29 Max(alpha1) = 1 MULES: Solving for alpha1 Liquid phase volume fraction = 0.00786006 Min(alpha1) = 6.78795e29 Max(alpha1) = 1 Using alpha1 smooth and evaluate K DILUPBiCG: Solving for Ux, Initial residual = 0.586111, Final residual = 5.56451e07, No Iterations 13 DILUPBiCG: Solving for Uy, Initial residual = 0.281552, Final residual = 8.30934e07, 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.03802e05, No Iterations 116 time step continuity errors : sum local = 1.93059e07, global = 1.37252e08, cumulative = 5.59463e05 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.15263e05, No Iterations 161 time step continuity errors : sum local = 7.36691e08, global = 9.963e09, cumulative = 5.59363e05 DICPCG: Solving for p_rgh, Initial residual = 0.000386962, Final residual = 1.70171e05, No Iterations 15 DICPCG: Solving for p_rgh, Initial residual = 2.10838e05, Final residual = 9.12595e08, No Iterations 246 time step continuity errors : sum local = 5.816e10, global = 3.18267e12, cumulative = 5.59363e05 ExecutionTime = 234.19 s ClockTime = 234 s Courant Number mean: 0.000148511 max: 0.174501 Interface Courant Number mean: 6.46745e05 max: 0.174501 deltaT = 1.51997e07 Time = 0.000245267 MULES: Solving for alpha1 Liquid phase volume fraction = 0.00786006 Min(alpha1) = 6.78794e29 Max(alpha1) = 1 MULES: Solving for alpha1 Liquid phase volume fraction = 0.00786007 Min(alpha1) = 6.78793e29 Max(alpha1) = 1 Using alpha1 smooth and evaluate K DILUPBiCG: Solving for Ux, Initial residual = 0.741821, Final residual = 5.45149e07, No Iterations 14 DILUPBiCG: Solving for Uy, Initial residual = 0.313413, Final residual = 6.99609e07, 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.01163e05, No Iterations 30 time step continuity errors : sum local = 6.38824e07, global = 5.72751e08, cumulative = 5.59936e05 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.44036e05, No Iterations 225 time step continuity errors : sum local = 1.9181e07, global = 3.94815e09, cumulative = 5.59896e05 DICPCG: Solving for p_rgh, Initial residual = 0.0011261, Final residual = 5.27881e05, No Iterations 8 DICPCG: Solving for p_rgh, Initial residual = 5.96673e05, Final residual = 9.5319e08, No Iterations 259 time step continuity errors : sum local = 7.3896e10, global = 8.87931e13, cumulative = 5.59896e05 ExecutionTime = 234.39 s ClockTime = 235 s Regards andrea 

May 19, 2011, 12:27 

#8 
Member
Sabin Ceuca
Join Date: Mar 2010
Location: Munich
Posts: 42
Rep Power: 8 
Hi,
try to switch from fvSolution > PISO > momentumPredictor yes; to no. Bye 

May 20, 2011, 03:47 

#9 
Senior Member
Andrea Ferrari
Join Date: Dec 2010
Posts: 296
Rep Power: 8 
Hi,
that is not the solution. Again the time step is of the order of 1e7s. 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 

May 20, 2011, 04:42 

#10 
Member

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 

May 20, 2011, 05:18 

#11 
Senior Member
Andrea Ferrari
Join Date: Dec 2010
Posts: 296
Rep Power: 8 
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 1e9/1e10 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 

May 20, 2011, 08:13 

#12 
Senior Member
PeiYing Hsieh
Join Date: Mar 2009
Posts: 311
Rep Power: 10 
Hi,
There is a thread on parasitic currents: parasitic currents This worked for some. Maybe it will work for you too. Pei 

May 20, 2011, 09:29 

#13 
Senior Member
Santiago Marquez Damian
Join Date: Aug 2009
Location: Santa Fe, Santa Fe, Argentina
Posts: 430
Rep Power: 15 
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. Research Scientist Research Center for Computational Methods (CIMEC)  CONICET/UNL Tel: 543424511594 Int. 7032 Colectora Ruta Nac. 168 / Paraje El Pozo (3000) Santa Fe  Argentina. http://www.cimec.org.ar 

May 23, 2011, 10:17 

#14 
Senior Member
Andrea Ferrari
Join Date: Dec 2010
Posts: 296
Rep Power: 8 
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 1e6 s) and also the velocities are smaller, but the interface is too diffuse (1520 cells instead of 23). 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 

May 23, 2011, 11:13 

#15 
Senior Member
Santiago Marquez Damian
Join Date: Aug 2009
Location: Santa Fe, Santa Fe, Argentina
Posts: 430
Rep Power: 15 
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 nonlinear 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. Research Scientist Research Center for Computational Methods (CIMEC)  CONICET/UNL Tel: 543424511594 Int. 7032 Colectora Ruta Nac. 168 / Paraje El Pozo (3000) Santa Fe  Argentina. http://www.cimec.org.ar 

May 24, 2011, 04:19 

#16 
Senior Member
Andrea Ferrari
Join Date: Dec 2010
Posts: 296
Rep Power: 8 
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 

May 24, 2011, 05:25 

#17  
Senior Member
Anton Kidess
Join Date: May 2009
Location: Delft, Netherlands
Posts: 1,121
Rep Power: 20 
Quote:
 Anton 

July 11, 2011, 08:19 

#18 
Senior Member
Andrea Ferrari
Join Date: Dec 2010
Posts: 296
Rep Power: 8 
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), (i1,j), (i,j+1), (i,j1) for an orthogonal mesh...but please correct me if i'm wrong) 2) The second question concerns this paper: "A balancedforce 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 

November 13, 2012, 04:22 

#19 
New Member
Farhad N.
Join Date: Apr 2010
Posts: 6
Rep Power: 8 
Andrea !
I wonder if you finally made to run your model, I am dealing the sam problem how have you solved your problem ? 

November 13, 2012, 05:12 

#20 
Senior Member
Andrea Ferrari
Join Date: Dec 2010
Posts: 296
Rep Power: 8 
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 

Thread Tools  
Display Modes  


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 