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++) 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 |
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, |
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 |
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, |
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 |
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 ;) |
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 Regards andrea |
Hi,
try to switch from fvSolution > PISO > momentumPredictor yes; to no. Bye |
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 |
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 |
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 |
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 |
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. |
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, 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. |
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 |
Quote:
- Anton |
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 ? |
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 |
interfoam performance
Hi to all,
I have a similar problem whit the performance of inteFoam. My problem concern a circular pipe filled with liquid and closed to both extremity. At initial time the right wall of the pipe is opened and gas allowing gas to go inside and liquid outside. My aim is to validate OpenFOAM for this kind of problem so I compared its performance with Fluent. The comparison seems very good but there is a strange difference. In my simulations I use a fixed Courant number and adaptive time step. The difference between is that Fluent give me a deltaT around 10e-04 while OpenFOAM give me a deltaT around 10e-6.So the simulation run very slow with OpenFOAM...Why this different behavior?There i a way to speed simulation? I obtain a little accelerate through decreasing tolerance but I would obtain i larger time step if possibile and kwon why there is this difference between two code. thanks to all |
If you are using the solver settings from the tutorials, you may have the maxCo and maxAlphaCo both set to 0.5. I typically will push these to 2.0 and 1.0 (at least) and then set nAlphaSubCycles to 5 in fvSolutions. Depending on how you have configured the BCs what is likely happening is that you are getting high velocities in the gas phase (probably inflow on a pressure outlet) that are driving down the dt. Setting maxCo to a higher value than maxAlphaCo will let the time step be controlled more by the velocities near the interface.
|
Thanks for your reply
I increase my Courant number to 0.8 in order to obtained a faster solution but with this Co I obtained a fluctuating solution. In fact, for this kind of problem the Courant number must not be more than 0.5, other experience says that courant number must not be more than 0.35 so I set its to 0.25. Regarding the velocity fields it seems ok. With the calculated velocity the courant number should be about 10e-04 but is 10e-06 |
This precisely why you use sub-time stepping on alpha. So if you set maxAlphaCo to 0.5 and nAlphaSubCycles is 3 then you have a Co per alpha step of 0.5/3 = 0.167. This is why I say use 1.0/5 = 0.2. You could push it to 1.5/5 = 0.3 and still probably be OK.
|
I set nAlphaSubCycles=2.
Two questions about you reply. 1)There is a way to verify if the real Courant numeber is Co/nAlphaSubCycles? Beacause in the output file the max Courant number remain fixed to 0.25 and not to 0,1225 2)I know that nAlphaSubCycles are the subcycle inside alpha equation. Through subcycle we have a more stable solution without incresing time step. Why nAlphaSubCycles should have to change Courant number? |
So when you say Courant number are you meaning maxCo or maxAlphaCo? interFoam should print out both values.
Here, "Courant Number" is limited by the maxCo setting in controlDict and "Interface Courant Number" is limited by maxAlphaCo. In the VOF-style algorithm, only the Co at the interface is the problem. The Co has to be small enough that it takes ~4 time 'steps' for the interface to travel through a cell. So it comes from the local mesh size and local velocity--not sure if you are taking that into account when calculating by hand. Since the solution of the volume fraction transport will be the limiting thing on time step size AND since all the computational time is spent in the pressure solve with the alpha solve being relatively cheap, you solve the alpha eqn multiple times with smaller time-steps for every one pressure solve. This is nAlphaSubCycles. Anyway, if you have nAlphaSubCycles set to 2 and you put the maxAlphaCo to 0.8 then you are having and effective max Co on each alpha solve of 0.4 which is unstable as you have seen. In any case, the value printed at each timestep is the OVERALL Co so to get the Co on a given sub-time step of alpha you divide this by nAlphaSubCycles. It changes Co because the time step size on each sub-cycle step is dt/nAlphaSubCycles. If you really want to check what the Co field looks like at any point in time, you would have to tweak the solver or add in a functionObject. One other thing to keep in mind is that things get crazy (interface oscillations and parasitic currents) if cAlpha is greater than 1.0 (it may be 2.0 in some of the old tutorial cases). Hope this helps. -Kent |
when I talking about Courant number I mean maxCo.
So you suggest to set maxCo=0.9 (explicit solver) and maxAlphaCo=2 with nAlphaSubCycles=10? thanks you very much for your replies |
No, I suggest:
maxCo=2.0 maxAlphaCo=1.0 nAlphaSubCycles=5 |
I understand why I must set maxAlphaCo=1.0 and nAlphaSubCycles=5.
last question: interFoam is an explicit solver so maxCo have to set <1,why you suggest maxCo=2.0? the nAlphaSubCycles acts also in p and U equations?I thinked no |
Hi all,
Could anyone please explain why we use MULES in alpha equation and not solve method? Thanks, Hrushi |
The volume fraction equation has special requirements with regard to boundedness, conservation and diffusivity. Thus we use a special algorithm to solve it.
|
Hi Anton,
Thank you for a quick reply. I have couple of follow up questions: 1. Do you know any reference that discusses these special requirements? I have read Rusche's thesis but did not find anything for this. 2. Will these requirements hold true when we change the range of fraction from 0:1 to -1:1? Thanks Hrushi |
Rusche's thesis does contain a brief section on how to deal with the issues - 3.2.6. Any paper dealing with the VoF method or how to improve it likely has information on the subject. In short:
1. Conservation - As an example if you have a bubble, you do not want the bubble volume to grow or shrink. It may significantly change the outcome of your simulation. 2. Boundedness - since all material properties are interpolated between the two phases, unboundedness can very quickly lead to divergence (e.g. due to negative density) 3. Diffusivity - Numerical diffusivity will eventually get rid of your fluid interface if not dealt with properly. Of course that's not acceptable for immiscible fluids. |
Well, solve with the proper advection TVD schemes should also do the work. In fact TVD schemes were designed to do exactly that. The problem is that, in practice, this is only true in 1D not in multi-dimensions, in this case FCT based schemes behave better. I discuss the topic in my PhD. thesis, that is provided in my user page at openfoamwiki.net: http://openfoamwiki.net/index.php/User:Santiagomarquezd
Regards. |
Hi Anton and Santiago,
Thank you for your quick reply. Really appreciate it. I am trying to solve a modified alpha equation with interfoam solver. I have no compression velocity and a chemical potential term in my alpha equation. I could not get any results with solve method, which is why I posted the question. Yesterday, I tried solving the alpha equation using mules. for (int aCorr=0; aCorr<nAlphaCorr; aCorr++) { surfaceScalarField phiAlpha (fvc::flux(phi,alpha1,alphaScheme)); /* solve ( fvm::ddt(alpha1)+ (fvc::div(phiAlpha)) - fvc::div(diffflx) );*/ MULES::implicitSolve(geometricOneField(), alpha1, phi, phiAlpha, volScalarField("Sp",mobility*fvc::laplacian(G)), zeroField(), 1, 0); rhoPhi = phiAlpha*(rho1 - rho2) + phi*rho2; } But my output says that I have maximum of alpha1 more than 1 at the start, which gradually reaches one but still stays more than 1 at equilibrium. Due to this, my results are really weird, and I do not know how to proceed. Any ideas? Thanks, Hrushi |
Hi Andrea :-
Now, as we have the fifth version of OpenFOAM, do you think we overcame that problem? I am trying to simulate the water flow in a 20-micron capillary tube but admirably, I shocked when running the capilaryRise example on OF5 and have the interface oscillations ( and it is just a simple case of 20 mm channel 2D flow!). When I run mine 20-micron case with very fine mesh, I have overly increased or decreased alpha values. what are your recommendations? my best! |
All times are GMT -4. The time now is 07:22. |