[mkraposhin]

Originally Posted by chriss85

Right now I'm using a simple mesh with orthogonal cells only. I have used surface-snapped meshes before but I used to get slower calculation speeds because the variance in cell sizes leads to higher maximum Courant numbers.

It will be interested to test solver on 3D meshes with tetrahderal elements and high non-orthogonality

[chriss85]

Which schemes do you suggest? VanLeer? I tried a few combinations in the shockTube and noticed that there are some requirements on the schemes to get meaningful results. Can you elaboarate on this topic?

[mkraposhin]

I found that usage of vector modifications of limiters (vanLeerV, vanAlbadaV and so on) can lead sometimes to instabilities or oscillations. Also, i found that vanAlbada limiter is more stable, than vanLeer.

I'm planning to test scheme on tetrhahderal meshes next week or later

[chriss85]

I've noticed a different problem now where I see an increase in temperature over the whole case
Right now I believe this is a problem in the first initial step. I will see an increase in temperature immediately, afterwards it remains the same. If I resume the calculation after stopping it, temperature increases again.

[mkraposhin]

Please, check that your are using enthalpy in termophysical properties and that you are specifiying Cp (not Cv) for energy equation

[chriss85]

Yes my error was there, I just figured this out myself as well. My code was still using inner energy instead of enthalpy in thermo.calculate() function...I will have to check this to make sure it works for both cases, but now the results are fine atleast

[mkraposhin]

Originally Posted by chriss85

Yes my error was there, I just figured this out myself as well. My code was still using inner energy instead of enthalpy in thermo.calculate() function...I will have to check this to make sure it works for both cases, but now the results are fine atleast

My advice is to use only enthalpy as primary variable for energy conservation, because it is much easier to account for expansion/contraction work - dp/dt instead of div(p*U) for internal energy.

where dp/dt - is the partial derivative of p by time

[chriss85]

Can you explain why you call field.oldTime() in your solver? I don't see this in e.g. rhoCentralFoam or sonicFoam.

[mkraposhin]

Can you explain why you call field.oldTime() in your solver? I don't see this in e.g. rhoCentralFoam or sonicFoam.

I do it, because i want to be sure, that old time values of these fields are stored in a proper way. Sometime, they are not stored and current values are used for explicit temporal derivatives, which leads to wrong results

[chriss85]

I have not seen this effect before but I don't know exactly how the mechanism works for caching the old field values.

[mkraposhin]

Originally Posted by chriss85

I have not seen this effect before but I don't know exactly how the mechanism works for caching the old field values.

.oldTime() method is executed when you are calling fvm::ddt operator or GeometricField::= () operator. But sometimes, when you update your fields manually, old time values are not cached, and current (new time step) and old time step values becomes equal. I found that this can happen for fvm::ddt(psi,p) .

That's why i forced storing of values at old time step before any changes to fields were made

[chriss85]

Oh, this is good to know, I need to figure out if this affects other equations in my code. I think I only have time derivatives in the flow and I don't use relaxation at the moment though.

[mkraposhin]

Dear FOAMers, i'm happy to announce next update of hybrid Kurganov-Tadmore / PISO Scheme. We tried this scheme for several industrial applications and found that it gives reasonably good results.

Examples are stored in git
https://github.com/unicfdlab/realLifeExamples
You tube video can be viewed here
Convergin-Diverging Nozzle http://www.youtube.com/watch?v=QgtsqaMp6dw
Compressor Simulation https://youtu.be/Egf8vtIGJL8

Also, we successfully applied this scheme and MULES for reacting flows, solver (reactingCentralFoam) can be downloaded from
https://github.com/unicfdlab/reactingCentralFoam

To run examples, you need additional boundary conditions (for example for cases whith subsonic inlet and partially supersonic outlet), this B.C. can be downloaded here
https://github.com/unicfdlab/libcompressibleTools

[chriss85]

Thanks for keeping up the good work

Have you seen improvements in meshes with bad quality? I think I will have to update my code to include the latest changes, judging by the commits in the repo.

I have experienced pressure and temperature dropping below ambient conditions sometimes as a result of a wrong calculation, but I am not sure if this is caused by your solver, as I have already seen something similar with sonicFoam before. Restarting the calculation before this point with a smaller time step size usually fixes this.

[chriss85]

Here are some screenshots of the minimum temperature and pressure curves I was talking about in the last post. In my thermophysical models I clamp the values so they can't become negative.

This is actually also happening on an orthogonal mesh:

    Overall domain bounding box (0 0 0) (0.0501 0.0351 0.0025)
    Mesh (non-empty, non-wedge) directions (1 1 1)
    Mesh (non-empty) directions (1 1 1)
    Boundary openness (4.844674461e-19 -1.299500967e-16 -1.35506872e-15) OK.
    Max cell openness = 1.918099971e-16 OK.
    Max aspect ratio = 4.014832162 OK.
    Minimum face area = 1.64695608e-08. Maximum face area = 2.3004639e-07.  Face area magnitudes OK.
    Min volume = 4.1173902e-12. Max volume = 5.75115975e-11.  Total volume = 3.659847041e-06.  Cell volumes OK.
    Mesh non-orthogonality Max: 0 average: 0
    Non-orthogonality check OK.
    Face pyramids OK.
    Max skewness = 0.0006070333015 OK.
    Coupled point location match (average 0) OK

Edit:
Only a few cells are affected by this behavior. Since the Mach number depends on the temperature and pressure, I will get completely wrong Mach numbers in these cells (right now on the order of 100), so this might prevent the solver from returning to a somewhat valid result after such an error ocurred? Would it make sense to enforce stricter limits on temperature/pressure and possibly Mach number?

[mkraposhin]

Hi, i'm sorry for late response!

Well, i fixed many things, and one of the change is the new procedure for update of kappa field, which blends standard and KT scheme.

Also, pressure solution procedure has been updated

You can see here example for multicomponent flow

https://github.com/unicfdlab/reactin...OpenFOAM-2.3.0

I see, that your mesh is rather good, so it should not affect solution

If you can share solver code and simple example, then i shall try to check correctness of implementation of time integration procedure

UPD.
And yes, i tried that scheme for bad meshes with non-orthogonality ~70 degr and skewness ~2.0 and it works good

[chriss85]

Ok, then I can probably rule out the solver for my problem.
I don't think I can share code right now, it's rather complex and I'm probably not allowed to do so, however I could post some excerpts if this is helpful. What exactly do you want to check in the time integration procedure? For the flow solving, I basically use your code but I store the fields in lists inside a class since my solver supports multiple regions. There are some additional source terms from electromagnetic effects which are calculated before the flow at each time step. I'm also including radiation which is calculated in hEqn and an ablation model which is also done there. I will try to isolate this effect by disabling the various submodels.

[mkraposhin]

Hi, can you post code parts, which corresponds to:

- pressure equation and density update
- blending field kappa update
- energy equation
- continuity equation
- Courant evaluation

Also, can you check that this strange behaviour is not observed when heat and mass sources are zero?

[chriss85]

Here you go. I've removed some parts related to profiling and data output for better readability. I also noticed that I didn't yet include continuity error calculation for your solver, I will include that ASAP.

I'm going to perform some calculations with disabled source terms to see if I can find the cause there.

[mkraposhin]

Hi,

i studied files that you uploaded and went to next conculsions:

Density eqn - ok
Enthalpy equation - ok, but you need addtitional terms for diffusion for muticomponent fluids and for cases where Cp highly depends on pressure or temperature

      fvm::ddt(rho,h)
    + fvm::div(phiPos,h)
    + fvm::div(phiNeg,h)
    - fvm::laplacian(turb.alphaEff(), h)
    + fvc::laplacian(alphahEff*T,Cp)
     + fvc::div(dzetaPhi) 
    + EkChange
    - S_ohm
    - S_R
    - SeMetal
    - SePolymer
    - S_sheath
    ==
      dpdt
    + fvc::div( ((-turb.devRhoReff()) & U) )

Momentum equation - OK, but you don't need to recalculate pressure gradient from previous iteration

Where dzetaPhi - is a diffusion internal energy flux

   /*
     *
     * Diffusive fluxes of energy due to mass diffusion
     *
     */
    PtrList<volScalarField> hi (Y.size());
    forAll (hi, iCmpt)
    {
	hi.set
	(
	    iCmpt,
	    new volScalarField
	    (
		IOobject
		(
		    "h" + Y[iCmpt].name(),
		    runTime.timeName(),
		    mesh,
		    IOobject::NO_READ,
		    IOobject::NO_WRITE
		),
		mesh,
		dimEnergy/dimMass
	    )
	);
	
	scalarField& hiIF = hi[iCmpt].internalField();
	const scalarField& pIF  = p.internalField();
	const scalarField& TIF  = T.internalField();
	
	forAll(hiIF, iCell)
	{
	    hiIF[iCell] = thermo.composition().Hs(iCmpt, pIF[iCell], TIF[iCell]);
	}
	
	forAll(hi[iCmpt].boundaryField(), iPatch)
	{
	    fvPatchScalarField& hip = hi[iCmpt].boundaryField()[iPatch];
	    const fvPatchScalarField& pp = p.boundaryField()[iPatch];
	    const fvPatchScalarField& Tp = T.boundaryField()[iPatch];
	    forAll(hip, iFace)
	    {
		hip[iFace] = thermo.composition().Hs(iCmpt, pp[iFace], Tp[iFace]);
	    }
	}
    }
    
    surfaceScalarField dzetaPhi
    (
	"dzetaPhi",
	fvc::flux
	(
	    diffusiveFlux[0],
	    hi[0],
	    "div(rhoi*Uri,hi)"
	)
    );
    
    for (label iCmpt = 1; iCmpt < Y.size(); iCmpt++)
    {
	dzetaPhi +=
	fvc::flux
	(
	    diffusiveFlux[iCmpt],
	    hi[iCmpt],
	    "div(rhoi*Uri,hi)"
	);
    }

    volScalarField alphahEff ("alphahEff", turbulence->alphaEff());
    volScalarField Cp ("Cp", thermo.Cp());

For pressure equations, i would suggest several changes, in next message