
[Sponsors] 
All Mach number implicit solver with KurganovTadmore scheme  pisoCentralFoam 

LinkBack  Thread Tools  Search this Thread  Display Modes 
October 8, 2015, 06:39 

#61 
Senior Member
Matvey Kraposhin
Join Date: Mar 2009
Location: Moscow, Russian Federation
Posts: 355
Rep Power: 20 
For pressure equation, it will be better to:
1) Move evaluation of phi to the branch of last non orthogonal corrector: Code:
if (pimple.finalNonOrthogonalIter()) { phiPos = pEqn_pos.flux(); phiNeg = pEqn_neg.flux(); phi = phiPos + phiNeg; } Code:
aphiv_pos = phiPos / (p_pos * psi_pos); aphiv_neg = phiNeg / (p_neg * psi_neg); amaxSf = max(mag(aphiv_pos), mag(aphiv_neg)); surfaceScalarField amaxSfbyDelta ( mesh.surfaceInterpolation::deltaCoeffs()*amaxSf ); surfaceScalarField Maf ( mag(phi) / (psi_pos*p_pos*a_pos + psi_neg*p_neg*a_neg) / (cSf_pos*a_pos + cSf_neg*a_neg) ); Info << "max/min Maf: " << max(Maf).value() << "/" << min(Maf).value() << endl; kappa = min ( Maf / (amaxSfbyDelta/mesh.magSf() * runTime.deltaT()), scalar(1.0) ); forAll(kappa.boundaryField(), iPatch) { fvsPatchField<scalar>& kappapf = kappa.boundaryField()[iPatch]; if (isA<coupledFvsPatchField<scalar> > (kappapf)) { forAll (kappapf, iFace) { kappapf[iFace] = 0.0; } } } Info << "max / min kappa: " << max(kappa).value() << "/" << min(kappa).value() << endl; phiPos = phiPos + (1.0  kappa) * phiNeg; phiNeg = kappa * phiNeg; https://github.com/unicfdlab/reactingCentralFoam Also, you will find simple test 'mixing channel' in that git, and you can try to run your solver in this test to see check energy conservation Last edited by mkraposhin; October 8, 2015 at 06:42. Reason: forgot link for multicomponent version 

October 8, 2015, 09:30 

#62 
Senior Member
Join Date: Oct 2013
Posts: 397
Rep Power: 17 
Thank you for your extensive comments, much appreciated! I will need some time to understand your posts (need to finish another issue first), but I will get back to you as soon as possible.
One thing I can say right now is that the other species don't influence any properties right now (missing material data ), I only use them as tracers. This might change in the future though. 

October 19, 2015, 07:50 

#63 
Senior Member
Join Date: Oct 2013
Posts: 397
Rep Power: 17 
Can you explain the meaning of the fvc::laplacian(alphahEff*T,Cp) term?
My heat capacity depends very much on pressure and temperature, so I guess this might be relevant. 

October 20, 2015, 05:56 

#64 
Senior Member
Matvey Kraposhin
Join Date: Mar 2009
Location: Moscow, Russian Federation
Posts: 355
Rep Power: 20 
Hi, as you know,
heat flux q =  lambda * grad(T), also, dh = Cp*dT, and if you will construct diffusive flux of enthalpy, you will get alpha * grad(h) = (lambda / Cp) (1 / Cp) grad(T) + alpha * T * grad(Cp) Thus, if you want to express diffusive heat flux in enthalpy, you need next expression: q =  alpha*grad(h) + alpha * T * grad(Cp)
__________________
MDPI Fluids (Q2) special issue for OSS software: https://www.mdpi.com/journal/fluids/..._modelling_OSS GitHub: https://github.com/unicfdlab Linkedin: https://linkedin.com/in/matveykraposhin413869163 RG: https://www.researchgate.net/profile/Matvey_Kraposhin 

October 21, 2015, 10:11 

#65 
Senior Member
Join Date: Oct 2013
Posts: 397
Rep Power: 17 
Ok I get what you mean, thanks for that.
Just as a clarification, there was a small error in your formula, it should be alpha * grad(h) = (lambda / Cp) * Cp * grad(T) + alpha * T * grad(Cp) instead. Does this term disappear when the inner energy is used instead? Because I don't see it in the sonicFoam solver. I guess we would just have Cv instead and the heat capacities would not cancel out? Is this an error in sonicFoam or is it neglected because its value is commonly low? Also, is there a reason you didn't include this term in the default solver but only in the reacting one? Last edited by chriss85; October 21, 2015 at 11:15. 

October 21, 2015, 12:56 

#66 
Senior Member
Matvey Kraposhin
Join Date: Mar 2009
Location: Moscow, Russian Federation
Posts: 355
Rep Power: 20 
Hi,
thank you for correcting formula This term vanishes only when Cp or Cv is constant, what is the very common case when sonicFoam or pisoCentralFoam are used. But for multicomponent mixture, mixtures where Cp or Cv strongly depends on mass fractions, temperature or pressure, this term is not zero. Assuming this term to be equal to zero can lead to false flux.
__________________
MDPI Fluids (Q2) special issue for OSS software: https://www.mdpi.com/journal/fluids/..._modelling_OSS GitHub: https://github.com/unicfdlab Linkedin: https://linkedin.com/in/matveykraposhin413869163 RG: https://www.researchgate.net/profile/Matvey_Kraposhin 

October 22, 2015, 03:49 

#67 
Senior Member
Join Date: Oct 2013
Posts: 397
Rep Power: 17 
I have implemented your suggested changes now, but I am still seeing this problem where the values ocassionally drop to zero
What I haven't yet implemented is the pressure gradient limiter from your repo. Do you think this may be related? 

October 22, 2015, 04:26 

#68 
Senior Member
Matvey Kraposhin
Join Date: Mar 2009
Location: Moscow, Russian Federation
Posts: 355
Rep Power: 20 
No, I think you must start from examining what causes values in your simulations to drop down:
 Which effect leads to this behaviour  flow?, sources in equations? which sources?  Who diverges first? Pressure? Mass fraction? Temperature Of course, you must perform your checks on mesh with orthogonal, same size (aspect ratio ~ 1), nonskewed cells. I think, that maybe you made simple programming error in code and you don't see it because your code has become complex
__________________
MDPI Fluids (Q2) special issue for OSS software: https://www.mdpi.com/journal/fluids/..._modelling_OSS GitHub: https://github.com/unicfdlab Linkedin: https://linkedin.com/in/matveykraposhin413869163 RG: https://www.researchgate.net/profile/Matvey_Kraposhin 

October 22, 2015, 05:12 

#69 
Senior Member
Join Date: Oct 2013
Posts: 397
Rep Power: 17 
It's a possibility. I have just created an orthogonal mesh for the geometry I'm currently investigating and will check if this problem occurs here as well. Then I can try to simplify my model to locate the problem. This problem is somewhat hard to track down because it doesn't happen with every geometry.
Edit: I have fixed another bug that caused random crashes in the calculation of mass and energy source terms. I think this might actually be a bug in OpenFOAM itsself. When calling turbulence>kappaEff(patchI) (probably also other quantities), it is returning the boundary field of a temporary volScalarField. I think the volScalarField might get freed and this could lead to memory errors. I'm not sure if the tmp mechanism is smart enough here to keep the volScalarField in memory. Anyway this leads to random crashes at the first iteration step but not every time I try to execute it, so I guess this is some kind of memory error. The fix for this is getting the whole field and storing it while the patch values are used. So far I have not seen the problem with zero values ocurring on two different orthogonal meshes yet, but I will check if it occurs sometime later. Edit2: I'm now seeing some other crashes on orthogonal meshes. One occurs where the pressure equation cannot be solved. Before it everything seems to be fine, but the final residual becomes NaN and after that no sensible values result. Another crash seemed to have a slightly different cause. Here the enthalpy equation started with initial residual = 1 and after that the Mach number which was calculated directly before solving the pressure equation was way off already. Then the pressure equation fails to solve again. Here are some example logs. Please note that the pressure and temperature ranges and the source terms are within reasonable limits for my case. 1st case: Code:
Min/max T:305.698429 49545.39054 Min/max p:106671.6502 78716106.47 band 0 Min/max S_R:6.011111455e+14 1.287920047e+15 Min/max kappa0:1 442152.7186 band 1 Min/max S_R:1.35734141e+14 3.262267452e+14 Min/max kappa1:1 218791.9251 band 2 Min/max S_R:6.064729679e+13 1.456700554e+14 Min/max kappa2:1 65945.29143 band 3 Min/max S_R:7.389808006e+13 1.543602339e+14 Min/max kappa3:1 31630.64568 band 4 Min/max S_R:3.237639662e+13 6.344276695e+11 Min/max kappa4:1 544.8898607 band 5 Min/max S_R:3.85204465e+12 5.835938786e+11 Min/max kappa5:1 14016.68309 DILUPBiCG: Solving for h, Initial residual = 0.0003407718613, Final residual = 3.441480603e11, No Iterations 4 Calculating temperature Temperature calculation took 40 ms S_R_min = 8.901135056e+14, S_R_max = 1.531648129e+15 S_Ohm_min = 1.454568452e76, S_Ohm_max = 2.279551212e+15 mach_min = 0.003847853739, mach_max = 2.331126917, mach_avg = 0.351727519 Mach number calculation took 0.04 s. GAMG: Solving for p, Initial residual = 2.152126861e05, Final residual = nan, No Iterations 20 2nd case: Code:
Min/max T:302.2874232 46906.74148 Min/max p:102652.6941 53403746.71 band 0 Min/max S_R:5.230251072e+14 2.451615746e+15 Min/max kappa0:1 306117.5729 band 1 Min/max S_R:7.076602005e+13 3.74588935e+14 Min/max kappa1:1 161244.1424 band 2 Min/max S_R:2.450812831e+13 7.081934198e+13 Min/max kappa2:1 48453.6687 band 3 Min/max S_R:3.141509763e+13 1.111786329e+14 Min/max kappa3:1 23165.60069 band 4 Min/max S_R:1.539103448e+13 1.558906401e+11 Min/max kappa4:1 305.4959672 band 5 Min/max S_R:2.709040834e+12 1.383976203e+11 Min/max kappa5:1 7397.124793 DILUPBiCG: Solving for h, Initial residual = 1, Final residual = 1.936204819e10, No Iterations 5 Calculating temperature Temperature calculation took 370 ms S_R_min = 6.600383089e+14, S_R_max = 2.740915073e+15 S_Ohm_min = 1.333567423e85, S_Ohm_max = 2.835578292e+15 mach_min = 386576.0505, mach_max = 3.144371603e+16, mach_avg = 5.632320259e+13 Mach number calculation took 0.11 s. GAMG: Solving for p, Initial residual = 0.8781873941, Final residual = nan, No Iterations 20 Do you think we could write some external wrapper that detects such a problem and switches the linear solver automatically? This would require some parsing or solver modifications I suppose, but might greatly increase stability. Alternatively, one could implement a safety guard directly in the linear solver code that tries a different solver when convergence cannot be archieved. Last edited by chriss85; October 23, 2015 at 04:37. 

October 29, 2015, 07:36 

#70 
Senior Member
Matvey Kraposhin
Join Date: Mar 2009
Location: Moscow, Russian Federation
Posts: 355
Rep Power: 20 
Hi, don't think that linear algebra solver is responsible for such strange behaviour.
I think that the reason lies in a bad matrix (and solution). You can try to turn off momentum predictor after momentum equation and also you can try to increase number of PISO correctors. From my point of view, it is necessary to separate several parts of program: 1) KT/PISO hybrid scheme 2) New physical properties, that you implemented in your solver 3) External sources I can try to test part #1 if you will propose case for it
__________________
MDPI Fluids (Q2) special issue for OSS software: https://www.mdpi.com/journal/fluids/..._modelling_OSS GitHub: https://github.com/unicfdlab Linkedin: https://linkedin.com/in/matveykraposhin413869163 RG: https://www.researchgate.net/profile/Matvey_Kraposhin 

October 29, 2015, 07:39 

#71 
Senior Member
Join Date: Oct 2013
Posts: 397
Rep Power: 17 
I will try to find something, but (un)fortunately the problem isn't ocurring all the time. Right now I'm doing calculations on orthogonal meshes and I haven't seen this problem in the last ~10 calculations.
I think it might also be possible that the bug I mentioned in my last post concerning the user of thermo.kappa(patchI) might be responsible. If the solver is accessing a deallocated object here and it doesn't crash, then it might receive invalid input which could result in wrong values being introduced in the calculation. This might also explain why I haven't seen it lately. I will let you know if/when it occurs again and try to separate out my additions to the solver. 

October 30, 2015, 09:24 

#72 
Senior Member
Join Date: Oct 2013
Posts: 397
Rep Power: 17 
I have managed to find a reproducible case where in the original solver the pressure drops to about 0.1 Bar in a single cell on a nonorthogonal mesh and the solver crashes afterwards. I believe that the problem is similar to the ones I saw before with my modified version on other meshes. Note that this hasn't occured on orthogonal meshes as far as I remember, so it might be a numerical problem with some of the schemes. I will generate an orthogonal mesh based on the shape of this one to check.
Please run this case with the allRun.sh script until 5e7s on a single core. It should crash shortly before and the pressure value will be low. The nonorthogonality and skewness thresholds present in the latest version of your solver were not used yet because I haven't implemented them in my code so far. The case uses values which are above the values present in the mesh to avoid these features with the latest version. I couldn't post it as attachment because of the mesh size, sorry for the hoster. http://www.fileupload.net/download...crash.zip.html I just tried a short test of the threshold features and the problem ocurred a bit later. Do you think these can fix such a problem? What effects on result quality can be expected? Can you give some general guidelines on required mesh quality for your solver? 

October 30, 2015, 10:46 

#73 
Senior Member
Matvey Kraposhin
Join Date: Mar 2009
Location: Moscow, Russian Federation
Posts: 355
Rep Power: 20 
Hi,
Yes, with threshold values that you specified in fvSolution, solver will not work Code:
nonOrthogonalityThreshold 80; skewnessThreshold 5; Also, i can give you advice on schemes selection  if your solution diverges, switch to vanAlbada, or if it still diverges, then swith to Minmod scheme. Also, you must set nNonOrthogonalCorrectors to be larger then 0  1 or 2 is enough Code:
PIMPLE { nCorrectors 3; nNonOrthogonalCorrectors 1; momentumPredictor true; //false adds stability nonOrthogonalityThreshold 45; skewnessThreshold 0.5; } Code:
divSchemes { default none; div((devRhoReff&U)) Gauss linear; div((muEff*dev2(T(grad(U))))) Gauss linear 1; //momentum equation div(phiNeg,U) Gauss Minmod; div(phiPos,U) Gauss Minmod; //energy equation div(phiNeg,h) Gauss Minmod; div(phiPos,h) Gauss Minmod; div(phiNeg,Ek) Gauss Minmod; div(phiPos,Ek) Gauss Minmod; //continuity equation div(phid_neg,p) Gauss Minmod; div(phid_pos,p) Gauss Minmod; } laplacianSchemes { default Gauss linear corrected; } interpolationSchemes { default none; interpolate(rho) linear; interpolate((rho*U)) linear; reconstruct(psi) Minmod; reconstruct(p) Minmod; reconstruct(U) Minmod; reconstruct(Dp) Minmod; }
__________________
MDPI Fluids (Q2) special issue for OSS software: https://www.mdpi.com/journal/fluids/..._modelling_OSS GitHub: https://github.com/unicfdlab Linkedin: https://linkedin.com/in/matveykraposhin413869163 RG: https://www.researchgate.net/profile/Matvey_Kraposhin 

November 2, 2015, 04:00 

#74 
Senior Member
Join Date: Oct 2013
Posts: 397
Rep Power: 17 
Thanks for your support! In this case, the skewness threshold was needed for archieving a stable calculation. Do you suggest trying to stay at the edge of stability for better precision, or is it safe to user more stable thresholds?
One more thing, what is the (dis)advantage of using the different kind of Courant numbers in your solver? 

November 5, 2015, 10:29 

#75 
Senior Member
Join Date: Oct 2013
Posts: 397
Rep Power: 17 
I'm still getting too low pressures on a mesh which is only slightly nonorthogonal and skewed:
Code:
Mesh stats points: 195935 faces: 531165 internal faces: 487332 cells: 168047 faces per cell: 6.060786566 boundary patches: 6 point zones: 0 face zones: 0 cell zones: 0 Overall number of cells of each type: hexahedra: 165118 prisms: 0 wedges: 0 pyramids: 0 tet wedges: 0 tetrahedra: 0 polyhedra: 2929 Breakdown of polyhedra by number of faces: faces number of cells 6 261 9 2082 12 453 15 115 18 18 Checking topology... Boundary definition OK. Cell to face addressing OK. Point usage OK. Upper triangular ordering OK. Face vertices OK. Number of regions: 1 (OK). Checking patch topology for multiply connected surfaces... Patch Faces Points Surface topology elektrode2 825 896 ok (nonclosed singly connected) ignition 168 228 ok (nonclosed singly connected) isolation 560 627 ok (nonclosed singly connected) elektrode1 840 912 ok (nonclosed singly connected) outlet 41 59 ok (nonclosed singly connected) POM 41399 41642 ok (nonclosed singly connected) Checking geometry... Overall domain bounding box (0.0188051 0.0174541 0.00260018) (0.0161754 0.0175068 0.00439982) Mesh (nonempty, nonwedge) directions (1 1 1) Mesh (nonempty) directions (1 1 1) Boundary openness (2.937652783e15 3.668005964e15 4.747549016e17) OK. Max cell openness = 3.530366807e16 OK. Max aspect ratio = 1.134600845 OK. Minimum face area = 4.132480545e09. Maximum face area = 7.501945466e08. Face area magnitudes OK. Min volume = 3.013267064e13. Max volume = 1.928490921e11. Total volume = 8.868828418e07. Cell volumes OK. Mesh nonorthogonality Max: 26.7596668 average: 4.114130965 Nonorthogonality check OK. Face pyramids OK. Max skewness = 0.3333333333 OK. Coupled point location match (average 0) OK. Mesh OK. Code:
PIMPLE { nCorrectors 3; nNonOrthogonalCorrectors 2; momentumPredictor true; skewnessThreshold 0.5; nonOrthogonalityThreshold 45; } 

November 7, 2015, 06:08 

#76 
Senior Member
Matvey Kraposhin
Join Date: Mar 2009
Location: Moscow, Russian Federation
Posts: 355
Rep Power: 20 
Hi, i found that in official OpenFOAM releases Courant number estimated with different comparing to extend version. Sometimes approach, implemented in official version leads to diverging solution, so i implemented both of them
__________________
MDPI Fluids (Q2) special issue for OSS software: https://www.mdpi.com/journal/fluids/..._modelling_OSS GitHub: https://github.com/unicfdlab Linkedin: https://linkedin.com/in/matveykraposhin413869163 RG: https://www.researchgate.net/profile/Matvey_Kraposhin 

November 7, 2015, 06:11 

#77 
Senior Member
Matvey Kraposhin
Join Date: Mar 2009
Location: Moscow, Russian Federation
Posts: 355
Rep Power: 20 
This is strange.
Are you getting pressures that are lower then 0.1Bar in your own solver or in pisoCentralFoam? Are sources active? Is your geometry is like what you posted above? Can you try next seetings? Code:
PIMPLE { nCorrectors 2; nNonOrthogonalCorrectors 1; momentumPredictor false; skewnessThreshold 0.1; nonOrthogonalityThreshold 15; }
__________________
MDPI Fluids (Q2) special issue for OSS software: https://www.mdpi.com/journal/fluids/..._modelling_OSS GitHub: https://github.com/unicfdlab Linkedin: https://linkedin.com/in/matveykraposhin413869163 RG: https://www.researchgate.net/profile/Matvey_Kraposhin 

November 9, 2015, 07:07 

#78  
Senior Member
Join Date: Oct 2013
Posts: 397
Rep Power: 17 
Quote:
Regarding my last comment with the wrong results, I have not tested the same case on the original pisoCentralFoam solver yet, but I believe the issue to be the same thing. I was under the impression that this should not happen with slightly nonorthogonal meshes. I can also tell you that my custom source terms don't produce this behaviour when I use sonicFoam as a flow solver, atleast I never experienced this before switching to your solver. If needed, I can run one of these meshes with sonicFoam and my source terms to verify. I still have the code in place to switch flow solvers, might have to update one or two things to get it working with the most recent code though. 

November 9, 2015, 14:35 

#79  
Senior Member
Matvey Kraposhin
Join Date: Mar 2009
Location: Moscow, Russian Federation
Posts: 355
Rep Power: 20 
Quote:
Quote:
Also, you must check that solver works good with perfect gas EOS. The problem may be in EOS implementation and/or thransport and thermophysical properties
__________________
MDPI Fluids (Q2) special issue for OSS software: https://www.mdpi.com/journal/fluids/..._modelling_OSS GitHub: https://github.com/unicfdlab Linkedin: https://linkedin.com/in/matveykraposhin413869163 RG: https://www.researchgate.net/profile/Matvey_Kraposhin 

November 12, 2015, 10:24 

#80  
Senior Member
Join Date: Oct 2013
Posts: 397
Rep Power: 17 
Quote:
Code:
solve ( fvm::ddt(rho) + fvc::div(phi) ); 

Thread Tools  Search this Thread 
Display Modes  


Similar Threads  
Thread  Thread Starter  Forum  Replies  Last Post 
DPMFoam  Serious Error particleladen flow in simple geometric config  benz25  OpenFOAM Running, Solving & CFD  27  December 19, 2017 20:47 
foamextend_3.1 decompose and pyfoam warning  shipman  OpenFOAM  3  July 24, 2014 08:14 
Solver is finishing with huge Mach number  Fonzie  CFX  1  March 12, 2007 14:15 
High Mach number solver error  Luke  CFX  3  January 31, 2007 22:26 
TVD scheme at low Mach number  Axel Rohde  Main CFD Forum  5  August 6, 1999 02:01 