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

LinkBack  Thread Tools  Search this Thread  Display Modes 
August 31, 2015, 06:54 

#41 
Senior Member
Matvey Kraposhin
Join Date: Mar 2009
Location: Moscow, Russian Federation
Posts: 355
Rep Power: 20 
It will be interested to test solver on 3D meshes with tetrahderal elements and high nonorthogonality


August 31, 2015, 06:59 

#42 
Senior Member
Join Date: Oct 2013
Posts: 397
Rep Power: 17 
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?


August 31, 2015, 07:08 

#43 
Senior Member
Matvey Kraposhin
Join Date: Mar 2009
Location: Moscow, Russian Federation
Posts: 355
Rep Power: 20 
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 

September 1, 2015, 05:16 

#44 
Senior Member
Join Date: Oct 2013
Posts: 397
Rep Power: 17 
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. 

September 1, 2015, 05:25 

#45 
Senior Member
Matvey Kraposhin
Join Date: Mar 2009
Location: Moscow, Russian Federation
Posts: 355
Rep Power: 20 
Please, check that your are using enthalpy in termophysical properties and that you are specifiying Cp (not Cv) for energy equation


September 1, 2015, 06:04 

#46 
Senior Member
Join Date: Oct 2013
Posts: 397
Rep Power: 17 
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


September 1, 2015, 06:27 

#47  
Senior Member
Matvey Kraposhin
Join Date: Mar 2009
Location: Moscow, Russian Federation
Posts: 355
Rep Power: 20 
Quote:
where dp/dt  is the partial derivative of p by time 

September 14, 2015, 09:20 

#48 
Senior Member
Join Date: Oct 2013
Posts: 397
Rep Power: 17 
Can you explain why you call field.oldTime() in your solver? I don't see this in e.g. rhoCentralFoam or sonicFoam.


September 14, 2015, 09:29 

#49  
Senior Member
Matvey Kraposhin
Join Date: Mar 2009
Location: Moscow, Russian Federation
Posts: 355
Rep Power: 20 
Quote:


September 14, 2015, 09:31 

#50 
Senior Member
Join Date: Oct 2013
Posts: 397
Rep Power: 17 
I have not seen this effect before but I don't know exactly how the mechanism works for caching the old field values.


September 14, 2015, 09:45 

#51  
Senior Member
Matvey Kraposhin
Join Date: Mar 2009
Location: Moscow, Russian Federation
Posts: 355
Rep Power: 20 
Quote:
That's why i forced storing of values at old time step before any changes to fields were made 

September 14, 2015, 09:51 

#52 
Senior Member
Join Date: Oct 2013
Posts: 397
Rep Power: 17 
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.


September 28, 2015, 08:06 

#53 
Senior Member
Matvey Kraposhin
Join Date: Mar 2009
Location: Moscow, Russian Federation
Posts: 355
Rep Power: 20 
Dear FOAMers, i'm happy to announce next update of hybrid KurganovTadmore / 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 ConverginDiverging 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 Last edited by mkraposhin; September 29, 2015 at 09:45. 

October 2, 2015, 05:13 

#54 
Senior Member
Join Date: Oct 2013
Posts: 397
Rep Power: 17 
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. 

October 5, 2015, 04:57 

#55 
Senior Member
Join Date: Oct 2013
Posts: 397
Rep Power: 17 
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 (nonempty, nonwedge) directions (1 1 1) Mesh (nonempty) directions (1 1 1) Boundary openness (4.844674461e19 1.299500967e16 1.35506872e15) OK. Max cell openness = 1.918099971e16 OK. Max aspect ratio = 4.014832162 OK. Minimum face area = 1.64695608e08. Maximum face area = 2.3004639e07. Face area magnitudes OK. Min volume = 4.1173902e12. Max volume = 5.75115975e11. Total volume = 3.659847041e06. Cell volumes OK. Mesh nonorthogonality Max: 0 average: 0 Nonorthogonality check OK. Face pyramids OK. Max skewness = 0.0006070333015 OK. Coupled point location match (average 0) OK 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? Last edited by chriss85; October 5, 2015 at 09:28. 

October 5, 2015, 12:24 

#56 
Senior Member
Matvey Kraposhin
Join Date: Mar 2009
Location: Moscow, Russian Federation
Posts: 355
Rep Power: 20 
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...OpenFOAM2.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 nonorthogonality ~70 degr and skewness ~2.0 and it works good 

October 6, 2015, 05:26 

#57 
Senior Member
Join Date: Oct 2013
Posts: 397
Rep Power: 17 
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. 

October 6, 2015, 15:16 

#58 
Senior Member
Matvey Kraposhin
Join Date: Mar 2009
Location: Moscow, Russian Federation
Posts: 355
Rep Power: 20 
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? 

October 7, 2015, 03:43 

#59 
Senior Member
Join Date: Oct 2013
Posts: 397
Rep Power: 17 
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. 

October 8, 2015, 06:33 

#60 
Senior Member
Matvey Kraposhin
Join Date: Mar 2009
Location: Moscow, Russian Federation
Posts: 355
Rep Power: 20 
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) ) 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()); 

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 