CFD Online Logo CFD Online URL
www.cfd-online.com
[Sponsors]
Home > Forums > Software User Forums > OpenFOAM > OpenFOAM Programming & Development

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

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

Like Tree20Likes

Reply
 
LinkBack Thread Tools Search this Thread Display Modes
Old   August 31, 2015, 06:54
Default
  #41
Senior Member
 
mkraposhin's Avatar
 
Matvey Kraposhin
Join Date: Mar 2009
Location: Moscow, Russian Federation
Posts: 343
Rep Power: 18
mkraposhin is on a distinguished road
Quote:
Originally Posted by chriss85 View Post
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
mkraposhin is offline   Reply With Quote

Old   August 31, 2015, 06:59
Default
  #42
Senior Member
 
Join Date: Oct 2013
Posts: 394
Rep Power: 15
chriss85 will become famous soon enough
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?
chriss85 is offline   Reply With Quote

Old   August 31, 2015, 07:08
Default
  #43
Senior Member
 
mkraposhin's Avatar
 
Matvey Kraposhin
Join Date: Mar 2009
Location: Moscow, Russian Federation
Posts: 343
Rep Power: 18
mkraposhin is on a distinguished road
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
mkraposhin is offline   Reply With Quote

Old   September 1, 2015, 05:16
Default
  #44
Senior Member
 
Join Date: Oct 2013
Posts: 394
Rep Power: 15
chriss85 will become famous soon enough
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.
chriss85 is offline   Reply With Quote

Old   September 1, 2015, 05:25
Default
  #45
Senior Member
 
mkraposhin's Avatar
 
Matvey Kraposhin
Join Date: Mar 2009
Location: Moscow, Russian Federation
Posts: 343
Rep Power: 18
mkraposhin is on a distinguished road
Please, check that your are using enthalpy in termophysical properties and that you are specifiying Cp (not Cv) for energy equation
mkraposhin is offline   Reply With Quote

Old   September 1, 2015, 06:04
Default
  #46
Senior Member
 
Join Date: Oct 2013
Posts: 394
Rep Power: 15
chriss85 will become famous soon enough
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
chriss85 is offline   Reply With Quote

Old   September 1, 2015, 06:27
Default
  #47
Senior Member
 
mkraposhin's Avatar
 
Matvey Kraposhin
Join Date: Mar 2009
Location: Moscow, Russian Federation
Posts: 343
Rep Power: 18
mkraposhin is on a distinguished road
Quote:
Originally Posted by chriss85 View Post
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
mkraposhin is offline   Reply With Quote

Old   September 14, 2015, 09:20
Default
  #48
Senior Member
 
Join Date: Oct 2013
Posts: 394
Rep Power: 15
chriss85 will become famous soon enough
Can you explain why you call field.oldTime() in your solver? I don't see this in e.g. rhoCentralFoam or sonicFoam.
chriss85 is offline   Reply With Quote

Old   September 14, 2015, 09:29
Default
  #49
Senior Member
 
mkraposhin's Avatar
 
Matvey Kraposhin
Join Date: Mar 2009
Location: Moscow, Russian Federation
Posts: 343
Rep Power: 18
mkraposhin is on a distinguished road
Quote:
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
mkraposhin is offline   Reply With Quote

Old   September 14, 2015, 09:31
Default
  #50
Senior Member
 
Join Date: Oct 2013
Posts: 394
Rep Power: 15
chriss85 will become famous soon enough
I have not seen this effect before but I don't know exactly how the mechanism works for caching the old field values.
chriss85 is offline   Reply With Quote

Old   September 14, 2015, 09:45
Default
  #51
Senior Member
 
mkraposhin's Avatar
 
Matvey Kraposhin
Join Date: Mar 2009
Location: Moscow, Russian Federation
Posts: 343
Rep Power: 18
mkraposhin is on a distinguished road
Quote:
Originally Posted by chriss85 View Post
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
mkraposhin is offline   Reply With Quote

Old   September 14, 2015, 09:51
Default
  #52
Senior Member
 
Join Date: Oct 2013
Posts: 394
Rep Power: 15
chriss85 will become famous soon enough
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.
chriss85 is offline   Reply With Quote

Old   September 28, 2015, 08:06
Default
  #53
Senior Member
 
mkraposhin's Avatar
 
Matvey Kraposhin
Join Date: Mar 2009
Location: Moscow, Russian Federation
Posts: 343
Rep Power: 18
mkraposhin is on a distinguished road
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
hk318i, chriss85 and ssss like this.

Last edited by mkraposhin; September 29, 2015 at 09:45.
mkraposhin is offline   Reply With Quote

Old   October 2, 2015, 05:13
Default
  #54
Senior Member
 
Join Date: Oct 2013
Posts: 394
Rep Power: 15
chriss85 will become famous soon enough
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 is offline   Reply With Quote

Old   October 5, 2015, 04:57
Default
  #55
Senior Member
 
Join Date: Oct 2013
Posts: 394
Rep Power: 15
chriss85 will become famous soon enough
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:
Code:
    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?
Attached Images
File Type: png temperature.png (12.8 KB, 31 views)
File Type: png pressure.png (16.2 KB, 22 views)

Last edited by chriss85; October 5, 2015 at 09:28.
chriss85 is offline   Reply With Quote

Old   October 5, 2015, 12:24
Default
  #56
Senior Member
 
mkraposhin's Avatar
 
Matvey Kraposhin
Join Date: Mar 2009
Location: Moscow, Russian Federation
Posts: 343
Rep Power: 18
mkraposhin is on a distinguished road
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
mkraposhin is offline   Reply With Quote

Old   October 6, 2015, 05:26
Default
  #57
Senior Member
 
Join Date: Oct 2013
Posts: 394
Rep Power: 15
chriss85 will become famous soon enough
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.
chriss85 is offline   Reply With Quote

Old   October 6, 2015, 15:16
Default
  #58
Senior Member
 
mkraposhin's Avatar
 
Matvey Kraposhin
Join Date: Mar 2009
Location: Moscow, Russian Federation
Posts: 343
Rep Power: 18
mkraposhin is on a distinguished road
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?
mkraposhin is offline   Reply With Quote

Old   October 7, 2015, 03:43
Default
  #59
Senior Member
 
Join Date: Oct 2013
Posts: 394
Rep Power: 15
chriss85 will become famous soon enough
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.
Attached Files
File Type: h pisoCentralFoam.H (3.9 KB, 12 views)
File Type: h pisoCentral_hEqn.H (827 Bytes, 16 views)
File Type: h pisoCentral_rhoEqn.H (1.1 KB, 14 views)
File Type: h pisoCentral_UEqn.H (342 Bytes, 10 views)
File Type: h compressibleMultiRegionCourantNo.H (643 Bytes, 8 views)
chriss85 is offline   Reply With Quote

Old   October 8, 2015, 06:33
Default
  #60
Senior Member
 
mkraposhin's Avatar
 
Matvey Kraposhin
Join Date: Mar 2009
Location: Moscow, Russian Federation
Posts: 343
Rep Power: 18
mkraposhin is on a distinguished road
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
Code:
      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
Code:
   /*
     *
     * 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
mkraposhin is offline   Reply With Quote

Reply

Thread Tools Search this Thread
Search this Thread:

Advanced Search
Display Modes

Posting Rules
You may not post new threads
You may not post replies
You may not post attachments
You may not edit your posts

BB code is On
Smilies are On
[IMG] code is On
HTML code is Off
Trackbacks are Off
Pingbacks are On
Refbacks are On


Similar Threads
Thread Thread Starter Forum Replies Last Post
DPMFoam - Serious Error --particle-laden flow in simple geometric config benz25 OpenFOAM Running, Solving & CFD 27 December 19, 2017 20:47
foam-extend_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


All times are GMT -4. The time now is 13:27.