[mkraposhin]

For pressure equation, it will be better to:

1) Move evaluation of phi to the branch of last non orthogonal corrector:

                        if (pimple.finalNonOrthogonalIter())
                        {
                                phiPos = pEqn_pos.flux();
                                phiNeg = pEqn_neg.flux();
                                phi = phiPos + phiNeg;
                        }

2) Recalculate kappa with new volume fluxes:

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;

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

[chriss85]

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.

[chriss85]

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.

[mkraposhin]

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)

[chriss85]

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?

[mkraposhin]

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.

[chriss85]

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?

[mkraposhin]

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

[chriss85]

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:

 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:

 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.

[mkraposhin]

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

[chriss85]

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.

[chriss85]

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?

[mkraposhin]

Hi,

Yes, with threshold values that you specified in fvSolution, solver will not work

    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

PIMPLE
{
    nCorrectors                 3;
    nNonOrthogonalCorrectors    1;
    momentumPredictor           true; //false adds stability
    nonOrthogonalityThreshold   45;
    skewnessThreshold           0.5;
}

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;
}

[chriss85]

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?

[chriss85]

I'm still getting too low pressures on a mesh which is only slightly nonorthogonal and skewed:

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:

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.

[mkraposhin]

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

[mkraposhin]

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?

PIMPLE
{
    nCorrectors                 2;
    nNonOrthogonalCorrectors    1;
    momentumPredictor           false;
    skewnessThreshold            0.1;
    nonOrthogonalityThreshold    15;
}

[chriss85]

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.

[mkraposhin]

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.

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

[chriss85]

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

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?

solve
    (
        fvm::ddt(rho) + fvc::div(phi)
    );