# All Mach number implicit solver with Kurganov-Tadmore scheme - pisoCentralFoam

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

 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; } 2) Recalculate kappa with new volume fluxes: 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& kappapf = kappa.boundaryField()[iPatch]; if (isA > (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; You can compare your solver with reactingCentralFoam - multicomponent version of KT/PISO scheme 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/matvey-kraposhin-413869163 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 multi-component 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/matvey-kraposhin-413869163 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), non-skewed 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/matvey-kraposhin-413869163 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.441480603e-11, 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.454568452e-76, 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.152126861e-05, Final residual = nan, No Iterations 20 Switching the pressure solver from GAMG to PBiCG apparently fixed this problem. 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.936204819e-10, 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.333567423e-85, 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 Here I had to switch the enthalpy solver from PBiCG to GAMG to fix this problem. 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/matvey-kraposhin-413869163 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 5e-7s on a single core. It should crash shortly before and the pressure value will be low. The non-orthogonality 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.file-upload.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; Actually, FVM implemented in OpenFOAM will not work with non-orthogonality larger then 60-70, skewness larger 3 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/matvey-kraposhin-413869163 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 (non-closed singly connected) ignition 168 228 ok (non-closed singly connected) isolation 560 627 ok (non-closed singly connected) elektrode1 840 912 ok (non-closed singly connected) outlet 41 59 ok (non-closed singly connected) POM 41399 41642 ok (non-closed singly connected) Checking geometry... Overall domain bounding box (-0.0188051 -0.0174541 -0.00260018) (0.0161754 0.0175068 0.00439982) Mesh (non-empty, non-wedge) directions (1 1 1) Mesh (non-empty) directions (1 1 1) Boundary openness (-2.937652783e-15 3.668005964e-15 -4.747549016e-17) OK. Max cell openness = 3.530366807e-16 OK. Max aspect ratio = 1.134600845 OK. Minimum face area = 4.132480545e-09. Maximum face area = 7.501945466e-08. Face area magnitudes OK. Min volume = 3.013267064e-13. Max volume = 1.928490921e-11. Total volume = 8.868828418e-07. Cell volumes OK. Mesh non-orthogonality Max: 26.7596668 average: 4.114130965 Non-orthogonality check OK. Face pyramids OK. Max skewness = 0.3333333333 OK. Coupled point location match (average 0) OK. Mesh OK. Once the flow reaches the nonorthogonal cells the pressure begins to drop in some cells to somewhat around 0.1 Bar. I'm currently using the Minmod scheme and the following settings in fvSolution: Code: PIMPLE { nCorrectors 3; nNonOrthogonalCorrectors 2; momentumPredictor true; skewnessThreshold 0.5; nonOrthogonalityThreshold 45; } It's only happening on nonorthogonal meshes as far as I can tell.

November 7, 2015, 06:08
#76
Senior Member

Matvey Kraposhin
Join Date: Mar 2009
Location: Moscow, Russian Federation
Posts: 355
Rep Power: 20
Quote:
 Originally Posted by chriss85 One more thing, what is the (dis)advantage of using the different kind of Courant numbers in your solver?
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

 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/matvey-kraposhin-413869163 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:
 Originally Posted by mkraposhin 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
So would you suggest using the faceCourant option for better determination of the Courant number?

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:
 Originally Posted by chriss85 So would you suggest using the faceCourant option for better determination of the Courant number?
This may help, but i'm sure that problem in different place. You can detect infuence of massFlux Courant number determination algorithm with the next behaviour: velocity is not changing, but Co is increasing.

Quote:
 Originally Posted by chriss85 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.
Well, i think you must test your meshes with original pisoCentralFoam. Because i tested it with bad meshes and i did not get such strange results, including with mesh that you provided to me

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

November 12, 2015, 10:24
#80
Senior Member

Join Date: Oct 2013
Posts: 397
Rep Power: 17
Quote:
 Originally Posted by mkraposhin Hi, you can look into hEqn.H for the hint on how to implement advection-diffusion for scalar in pisoCentralFoam (rhoCentralFoam). Convection for scalar Y (mass fraction, for example) can be implemented like this Code: fvm::ddt(rho,Y) + fvm::div(phiPos, Y) + fvm::div(phiNeg,Y) - fvm::laplacian(turbulence->muEff / Sc, Y)
Just a short question on this topic. Why are you not separating the flux in positive/negative parts in the continuity equation?
Code:
solve
(
fvm::ddt(rho) + fvc::div(phi)
);