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

LinkBack  Thread Tools  Search this Thread  Display Modes 
August 25, 2015, 09:07 

#21 
Senior Member
Join Date: Oct 2013
Posts: 394
Rep Power: 15 
Ok, then I need to check if the definition of psi is correct in my solver. I'm using tabulated material models for everything. Right now, I have precalculated speed of sound to avoid calculating the derivative in each step, which is a bit expensive to calculate as it requires more table lookups. This solution is correct for my case, but if it's possible to keep the old formulation (with the same performance) it would be more elegant.


August 26, 2015, 04:12 

#22 
Senior Member
Join Date: Oct 2013
Posts: 394
Rep Power: 15 
I just checked sonicLiquidFoam (was a bit short on time yesterday). For EOS they use rho=rho0+psi(pp0) and I noticed that an additional term div(phi) with phi=(rho00/psi)*phid=(rho0/psip0)*phid was included in the pressure equation, which is likely used to account for the offsets in this formulation. In my case I have an equation of state based on tabular data. For low pressures and temperatures it approximates the ideal EOS. Up until now I have simply implemented psi as psi(p,T) = p / (rho(p,T) + SMALL) with rho(p,T) being an interpolating data table. Is this sufficient or is a correction to the pEqn needed here? I have not fully understood the calculation of the pressure yet unfortunately, so I'm a bit clueless here.


August 26, 2015, 05:10 

#23  
Senior Member
Matvey Kraposhin
Join Date: Mar 2009
Location: Moscow, Russian Federation
Posts: 343
Rep Power: 18 
Quote:
Yes, because you must implement psi as a coefficient between rho and p, rho = psi*p. In such assumption, terms d (rho) / dt `=` d(psi*p) / dt and d( rho U) / dx = d ( psi*U*p) / dx can be approximated with pisoCentralFoam scheme without changing pressure equation. No, because in your formulation psi is function of pressure and this adds additional nonlinearity to system. For example in my solver psi = 1 / (R*T) and it depends only on T. 

August 26, 2015, 05:31 

#24 
Senior Member
Join Date: Oct 2013
Posts: 394
Rep Power: 15 
Oh I just noticed that I made an error in my post, of course I implemented psi as rho(p, T) / (p + SMALL) and not the inverse.
I'm not sure I understood your post correctly. I agree that this kind of EOS introduces nonlinearity. The derivatives you talked about should still be valid with this EOS, as they don't derive in regards to p or rho. Is this nonlinearity in determining the pressure already being solved by using multiple corrector iterations, or are additional terms necessary like in sonicLiquidFoam? I guess I really need to spend some time trying to understand the derivation of the pressure, unfortunately I'm a one man show here 

August 26, 2015, 06:38 

#25  
Senior Member
Matvey Kraposhin
Join Date: Mar 2009
Location: Moscow, Russian Federation
Posts: 343
Rep Power: 18 
Quote:


August 26, 2015, 10:08 

#26 
Senior Member
Join Date: Oct 2013
Posts: 394
Rep Power: 15 
Yes that's clear. But do you think I can get a correct solution with the current implementation? And would I need more corrector iterations compared to a linear EOS?
Today I have done some reading on the wiki and now I believe I understand the incompressible PISO algorithm and the derivation of the pressure equation, but I have not transferred the method to the compressible case yet to compare it with the source code. 

August 27, 2015, 04:57 

#27  
Senior Member
Matvey Kraposhin
Join Date: Mar 2009
Location: Moscow, Russian Federation
Posts: 343
Rep Power: 18 
Quote:
to get equation for pressure, do next things: 1) put EoS in continuity equations: d(rho)/dt + d(rho U )/dx = d(psi*p)/dt + d(psi*p*U)/dx = d(psi*p)/dt + d(psi*U*p)/dx = S_rho 2) relate change in velocity with pressure gradient and other terms: A*U = H  d(p)/dx 3) put expression for velocity into continuity equation, formulated for pressure d(psi*p)/dt + d( (psi*H/A) p )/dx  d ( (psi*p/A) * d(p)/dx )/dx = S_rho Fluxes in the last equation multiplied by pressure are mass fluxes, satisfying continuity 

August 28, 2015, 04:07 

#28 
Senior Member
Join Date: Oct 2013
Posts: 394
Rep Power: 15 
Thanks for the hints, I will try to follow them.
In the meantime, I have done some tests on a real world application using my tabulated data with some inconclusive results. Please see the attached image. I have calculated an expanding flow using an initial high temperature, high pressure region at the top left of the geometry and calculated only flow and radiation. When I use 3 correctors I notice a change in the density, even in the region the shockwave hasn't reached yet. Any ideas about that? It's somewhat difficult to tell which result is best. 2nd order CrankNicolson scheme didn't seem to have much of an influence (I went completely 2nd order). However, both changing the maxCo value and changing the number of correctors has an influence on the results. Computing time is somewhat limited for me so I would prefer to go with the maxCo=0.25 cases, but it seems that the time step error is still too large? Edit: Forgot to mention that the mesh is completely orthogonal. Last edited by chriss85; August 28, 2015 at 06:35. 

August 28, 2015, 06:17 

#29 
Senior Member
Join Date: Oct 2013
Posts: 394
Rep Power: 15 
Your explanation together with this old thread helped me to understand the derivation of the pressure equation in compressible flow solvers like sonicFoam. This actually made me realize that I need the density source terms in the pressure equation, which I didn't put there yet
I think this also made me understand the problems involved with the nonlinearity in the material models here. Basically, the pressure equation is not completely implicit anymore because the psi(p,T) dependence is handled explicitly, while the terms in pEqn are implicit. I still haven't completely understood the significance of this though. Does it mean that my time step size is limited by this explicit treatment? Or would multiple correctors be able to fix this nonlinearity? I still need to check at which point psi is reevaluated. 

August 28, 2015, 08:58 

#30 
Senior Member
Join Date: Oct 2013
Posts: 394
Rep Power: 15 
Ok I understood the problem some more now. Psi gets updated inside of thermo.correct(), which is called in hEqn.H in your solver, outside of the PISO loop. Inside the PISO loop, rho is calculated with rho=psi*p, using the new p, but the compressibility from the old pressure value of the previous time step. I think this is what is causing the wrong pressures in the results I showed above. This means that thermo.correct needs to be called inside of the PISO loop to update the compressibility after the pressure equation has been solved. It's still not clear to me if this is a solution that will converge properly, and I'm also not sure if I need to call it inside the nonorthogonal corrector loop (would have to update a few terms defined outside of this loop then, possibly hurting performance). thermo.correct is also comparitively expensive to call because more than the compressibility is updated in there, so it might be better to update only the compressibility somehow. This would require a different structure of the thermophysical libraries though.
I'm going to run a few calculations with this change and I'll report back if this fixes the density errors. 

August 28, 2015, 09:46 

#31 
Senior Member
Matvey Kraposhin
Join Date: Mar 2009
Location: Moscow, Russian Federation
Posts: 343
Rep Power: 18 
Hi, give me 23 days to write answer


August 31, 2015, 03:30 

#32 
Senior Member
Join Date: Oct 2013
Posts: 394
Rep Power: 15 
Sure, you don't have to hurry
I have run some calculations in which I call thermo.correct() inside the PISO corrector loop and it seems to give consistent results when I use atleast 3 correctors. Using only 2 correctors gave me some differences in results. Everything above 2 correctors looks very similar. I think this can be regarded as a general flaw in many OpenFOAM solvers, as the thermophysical APIs explicitly allow a pressure dependence in the compressibility but most solvers won't be able to process this dependency properly. Do you think I should file a bug report? 

August 31, 2015, 04:27 

#33  
Senior Member
Matvey Kraposhin
Join Date: Mar 2009
Location: Moscow, Russian Federation
Posts: 343
Rep Power: 18 
Quote:
D(rho)/Dt = S_rho D(rhoU)/Dt = S_U  dp/dx D(rho*e)/Dt = S_e  d(p*U)/dx p = rho*R*T Implicit solvers of OpenFOAM are uses PISO method to solve this equations sequentially, for example: Code:
D(rho)/Dt = S_rho > new density D(rhoU)/Dt = S_U  dp/dx > new velocity prediction D(rho*e)/Dt = S_e  d(p*U)/dx > new energy and temperature rho = p / (R*T) and D(rho)/Dt = S_rho > new pressure and new velocity correction If you want to couple energy with pressure in PISO loop, then, maybe you have to couple density and other equations (such as mass fractions) with new fluxes. From my experience, it is better to couple only pressure and velocity with inner PISO iterations, while other variables can be coupled with outer iterations (PIMPLE). In other words, all variables, that affects compressibility (psi) must be computed once before PISO loop. 

August 31, 2015, 05:00 

#34  
Senior Member
Join Date: Oct 2013
Posts: 394
Rep Power: 15 
Quote:
Your solver (and other PISO based ones) currently do the following: Code:
continuity equation > new density velocity predictor > new velocity energy equation > new energy thermo.correct() > new temperature, compressibility, ... PISO loop { pressure equation > new pressure and velocity rho = psi*p > new density } Code:
continuity equation > new density velocity predictor > new velocity energy equation > new energy thermo.correct() > new temperature, compressibility, ... PISO loop { pressure equation > new pressure and velocity thermo.correct() > new temperature, compressibility, ... //calculating compressibility would suffice, but needs a change in the implementation. rho = psi*p > new density } 

August 31, 2015, 05:16 

#35  
Senior Member
Matvey Kraposhin
Join Date: Mar 2009
Location: Moscow, Russian Federation
Posts: 343
Rep Power: 18 
Quote:


August 31, 2015, 05:27 

#36 
Senior Member
Join Date: Oct 2013
Posts: 394
Rep Power: 15 
That's basically what I would like to do, however this would require modifications in the thermophysical library, or a manual implementation of the compressibility model in the solver. To keep things simple I call thermo.correct() right now. This has some overhead because data table lookups are somewhat expensive, but with 3 corrector iterations per time step it's less than a second per time step.
Do you think I should file a bug report for other solvers regarding this issue? The implementation of the thermophysical library API in which psi has p and T parameters hints that the solvers should support such a dependency, while in reality they don't. 

August 31, 2015, 05:32 

#37 
Senior Member
anonymous
Join Date: Aug 2014
Posts: 205
Rep Power: 10 
Well indeed if your PSI depends on p it should be better to correct the psi term after updating the pressure in each corrector. The problems are:
1) It can reduce the convergence rate if there are huge variations of psi with p 2) It increases the computational cost Either way you can maintain 2 nCorrectors in the PISO reducing the timeStep 

August 31, 2015, 06:09 

#38 
Senior Member
Join Date: Oct 2013
Posts: 394
Rep Power: 15 
Fortunately my psi only depends slightly on p, so there is not much of a problem with convergence. I would rather stay at large time steps to get shorter calculation times. Using maxCo=0.25 with 3 correctors gives me results which only show a small temporal error. Depending on the case I have Mach numbers ranging between 0 and 3. Using this solver I get a speedup between 24x compared to a sonicFoambased solver. I have also run some shockTube comparisons between this solver and sonicFoam and I get somewhat better results with this solver, so overall it seems to work better


August 31, 2015, 06:23 

#39  
Senior Member
Matvey Kraposhin
Join Date: Mar 2009
Location: Moscow, Russian Federation
Posts: 343
Rep Power: 18 
Quote:


August 31, 2015, 06:37 

#40 
Senior Member
Join Date: Oct 2013
Posts: 394
Rep Power: 15 
Right now I'm using a simple mesh with orthogonal cells only. I have used surfacesnapped meshes before but I used to get slower calculation speeds because the variance in cell sizes leads to higher maximum Courant numbers.
The mesh I'm using right now has about 400k cells. The speedup I'm observing is mostly because this solver avoids the additional iterations per time step. I'm also calculating radiation using the P1 model, which is taking 5060% of the calculation time per time step due to bad convergence. Avoiding the repeated calculation of the P1 equations helps the performance significantly here. 

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 