CFD Online Logo CFD Online URL
www.cfd-online.com
[Sponsors]
Home > Forums > 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 Tree5Likes

Reply
 
LinkBack Thread Tools Display Modes
Old   August 14, 2015, 09:15
Default All Mach number implicit solver with Kurganov-Tadmore scheme - pisoCentralFoam
  #1
Senior Member
 
mkraposhin's Avatar
 
Matvey Kraposhin
Join Date: Mar 2009
Location: Moscow, Russian Federation
Posts: 321
Rep Power: 10
mkraposhin is on a distinguished road
Dear OpenFOAM users!

We finished development of hybrid Kurganov-Tadmore/PISO semi-implicit solver for all-mach number flows

You can download it with tutorials at github https://github.com/mkraposhin/pisoCentralFoam

git link https://github.com/mkraposhin/pisoCentralFoam.git

Solver description located in git archive and wiki https://github.com/mkraposhin/pisoCentralFoam/wiki

Model was tested in Mach numbers from 0.01 to 6 for cases:
  1. Laminar flow in channel
  2. Laminar flow around cylinder
  3. Turbulent flow around cylinder
  4. Shock Tube test (Sod's problem)
  5. forward supersonic test
  6. backward step supersonic case
  7. Converging-Diverging nozzle
  8. Oblique shock

Model was implemented for OpenFOAM of next versions:
- OpenFOAM 2.3.0
- OpenFOAM 2.4.0
- OpenFOAM-dev (thanks to Hassan Kassem)
- foam-extend 3.1

Youtube video


forwardStep
ssss likes this.

Last edited by mkraposhin; August 29, 2015 at 14:54.
mkraposhin is offline   Reply With Quote

Old   August 14, 2015, 10:44
Default
  #2
Senior Member
 
anonymous
Join Date: Aug 2014
Posts: 191
Rep Power: 3
ssss is on a distinguished road
Thank you very much for this solver, and all the documentation and validation developed. This is a very interesting work.

I have some questions about the solver. Do you think that using the PIMPLE algorithm could help the accuracy of the solver or if it could help to use a higher timeStep? Furthermore, is this solver suitable for small unstructured meshes?
ssss is offline   Reply With Quote

Old   August 14, 2015, 10:48
Default
  #3
Senior Member
 
mkraposhin's Avatar
 
Matvey Kraposhin
Join Date: Mar 2009
Location: Moscow, Russian Federation
Posts: 321
Rep Power: 10
mkraposhin is on a distinguished road
This solver can't be used PIMPLE method.

This solver can be used with any grid. But, to work with non-orthogonal grids, nNonOrthogonalCorrectors must be set to 1 or higher
mkraposhin is offline   Reply With Quote

Old   August 17, 2015, 02:50
Default
  #4
Member
 
Join Date: Jun 2012
Posts: 68
Rep Power: 5
maHein is on a distinguished road
Very interesting work!

Have you tested the dissipation at very low Mach numbers? For instance, calculate the drag coefficient of a cylinder at various Mach numbers at constant Reynolds numbers and compare it to experimental values?

As far as I understood, your solver is not bound to the acoustic Courant number but rather the flow Courant number, am I correct?

Kind regards,

Martin
maHein is offline   Reply With Quote

Old   August 17, 2015, 03:42
Default
  #5
Senior Member
 
mkraposhin's Avatar
 
Matvey Kraposhin
Join Date: Mar 2009
Location: Moscow, Russian Federation
Posts: 321
Rep Power: 10
mkraposhin is on a distinguished road
Hi, thank you for your opinion

Quote:
Originally Posted by maHein View Post

Have you tested the dissipation at very low Mach numbers? For instance, calculate the drag coefficient of a cylinder at various Mach numbers at constant Reynolds numbers and compare it to experimental values?
Well i only found, that Cd of cylinder changes with Mach number. But i didn't compared results with experimental data. If you shall provide me with reference (pdf file), i shall try

Quote:
Originally Posted by maHein View Post
As far as I understood, your solver is not bound to the acoustic Courant number but rather the flow Courant number, am I correct?
Yes
mkraposhin is offline   Reply With Quote

Old   August 17, 2015, 07:11
Default
  #6
Senior Member
 
Join Date: Oct 2013
Posts: 331
Rep Power: 6
chriss85 is on a distinguished road
This is very interesting work! I haven't looked into it very much yet, but I plan do so.

A few questions:

1) Is it possible to modify this solver for use in a multiphysical problem? I have a couple of mass, impulse and energy source terms that I need to include, which in turn also depend on the flow (or on properties derived from it). As this solver is using a PISO approach, we cannot use a consistent coupling by using multiple iterations per timestep and solving the coupled equations in a segregated manner, so the coupling will use the lagged fields from the previous time step. Do you think this can still be a sensible approach?

2) I suppose that impulse source terms need to be added directly into UEqn (that is, before the momentum predictor)?

3) Which ranges for the flow Courant number are sensible to use with this solver? I would like to use an adaptive time step size using the flow Courant number. Do you see any issues here?

4) Is this solver suitable for running transient simulations where parts of the solution can have zero velocity? I have a very wide range of flow conditions (plasma simulation), high velocities (on the order of 1e3-1e4m/s, mostly subsonic due to higher speed of sound in plasma) and strong energy source terms leading to expansion. However, flow can be zero during some parts of the solution.

Last edited by chriss85; August 17, 2015 at 09:37.
chriss85 is offline   Reply With Quote

Old   August 17, 2015, 16:49
Default
  #7
Senior Member
 
mkraposhin's Avatar
 
Matvey Kraposhin
Join Date: Mar 2009
Location: Moscow, Russian Federation
Posts: 321
Rep Power: 10
mkraposhin is on a distinguished road
Hi, chriss85!

Let me answer your questions step by step

1) Yes, you can add sources to rho, U and h. I'm not sure that PIMPLE approach is appicable for such scheme. So, lagged approach is prefered. You can use fvOptions API to insert sources as it done in other OpenFOAM solvers

2) Yes, it is better to insert momentum source before momentum predictor step.

3) For subsonic flow Co must be less then 1, for transonic or supersonic Co < 0.5

4) This conditions are similar to those one, for which this solver was designed
mkraposhin is offline   Reply With Quote

Old   August 19, 2015, 03:58
Default
  #8
Senior Member
 
Join Date: Oct 2013
Posts: 331
Rep Power: 6
chriss85 is on a distinguished road
Thanks for your answers, much appreciated. One more thing, how does one best implement a scalar transport equation (with a source term) in such a solver? I have read that such an equation should preferably use the central upwind schemes here: rhoCentralFoam transport equation .
Should the equations be implemented like the equation for rho itself, or is the default scalar transport equation sufficient?
chriss85 is offline   Reply With Quote

Old   August 20, 2015, 07:22
Default
  #9
Senior Member
 
Hassan Kassem
Join Date: May 2010
Location: UK
Posts: 149
Rep Power: 7
hk318i is on a distinguished road
Great work!
I will try it for inviscid transonic flow instead of rhoCentralFoam. Although rhoCentralFoam gives good results, it is really hard to maintain a stable solution.

Did you published any results using this new solver? Could you please recommend any papers describe this algorithm, if available?
hk318i is offline   Reply With Quote

Old   August 20, 2015, 12:47
Default
  #10
Senior Member
 
mkraposhin's Avatar
 
Matvey Kraposhin
Join Date: Mar 2009
Location: Moscow, Russian Federation
Posts: 321
Rep Power: 10
mkraposhin is on a distinguished road
Quote:
Originally Posted by hk318i View Post
Great work!
I will try it for inviscid transonic flow instead of rhoCentralFoam. Although rhoCentralFoam gives good results, it is really hard to maintain a stable solution.

Did you published any results using this new solver? Could you please recommend any papers describe this algorithm, if available?
Hi, thank you for your opinion.

If you want to see more examples, then you must look into the Tutorials-2.3.0 folder, because it contains more examples, then Tutorials-2.4.0

I sent paper with description of solver to PROCEDIA journal and i hope they will accept it until november.

Anyway, i put draft version of description of alghorithm in git repository. You can find it here
https://github.com/mkraposhin/pisoCe...-CFDONLINE.pdf

Last edited by mkraposhin; August 20, 2015 at 12:50. Reason: grammar
mkraposhin is offline   Reply With Quote

Old   August 20, 2015, 12:59
Default
  #11
Senior Member
 
Hassan Kassem
Join Date: May 2010
Location: UK
Posts: 149
Rep Power: 7
hk318i is on a distinguished road
Thanks for sharing your solver. Please once your paper get accepted inform us so we could cite it.
I am thinking of implementing dynamicMesh in your solver but I am still studying it. Do you have an suggestions? Did you think about including a dynamic mesh is this algorithm before?
hk318i is offline   Reply With Quote

Old   August 20, 2015, 13:07
Default
  #12
Senior Member
 
mkraposhin's Avatar
 
Matvey Kraposhin
Join Date: Mar 2009
Location: Moscow, Russian Federation
Posts: 321
Rep Power: 10
mkraposhin is on a distinguished road
Quote:
Originally Posted by chriss85 View Post
Thanks for your answers, much appreciated. One more thing, how does one best implement a scalar transport equation (with a source term) in such a solver? I have read that such an equation should preferably use the central upwind schemes here: rhoCentralFoam transport equation .
Should the equations be implemented like the equation for rho itself, or is the default scalar transport equation sufficient?
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

Code:
fvm::ddt(rho,Y)
+
fvm::div(phiPos, Y)
+
fvm::div(phiNeg,Y)
-
fvm::laplacian(turbulence->muEff / Sc, Y)
Where
  • Y - mass fraction
  • rho - mixture density
  • phiPos - positive direction mass flux
  • phiNeg- negative direction mass flux
  • Sc - Schmidt number field

But this discreete equation must be obtained with the same limiter as used for energy equation, otherwise temperature of the mixture will not conserve

A better way may be in adaptation of MULES technique for this KT scheme. And i'm working on it

Last edited by mkraposhin; August 21, 2015 at 02:14.
mkraposhin is offline   Reply With Quote

Old   August 20, 2015, 13:15
Default
  #13
Senior Member
 
mkraposhin's Avatar
 
Matvey Kraposhin
Join Date: Mar 2009
Location: Moscow, Russian Federation
Posts: 321
Rep Power: 10
mkraposhin is on a distinguished road
Quote:
Originally Posted by hk318i View Post
Thanks for sharing your solver. Please once your paper get accepted inform us so we could cite it.
I am thinking of implementing dynamicMesh in your solver but I am still studying it. Do you have an suggestions? Did you think about including a dynamic mesh is this algorithm before?
I can include dynamic mesh for OF 2.3.0 and 2.4.0. This will take some time - about 1-2 weeks, depending on my busyness. And i need a good examples for testing.

First example can be forwardStep with inverted motion
mkraposhin is offline   Reply With Quote

Old   August 20, 2015, 13:27
Default
  #14
Senior Member
 
Hassan Kassem
Join Date: May 2010
Location: UK
Posts: 149
Rep Power: 7
hk318i is on a distinguished road
Quote:
Originally Posted by mkraposhin View Post
I can include dynamic mesh for OF 2.3.0 and 2.4.0. This will take some time - about 1-2 weeks, depending on my busyness. And i need a good examples for testing.

First example can be forwardStep with inverted motion
That will be great, there are many tutorials for pimpleDyMFoam. The simplest one is movingCone. The movingCone is available for pimpleDyMFoam(OF2.3-2.4 also), sonicDyMFoam and rhoCentralDyMFoam in OpenFOAM-dev
hk318i is offline   Reply With Quote

Old   August 24, 2015, 13:37
Default
  #15
Senior Member
 
mkraposhin's Avatar
 
Matvey Kraposhin
Join Date: Mar 2009
Location: Moscow, Russian Federation
Posts: 321
Rep Power: 10
mkraposhin is on a distinguished road
Hello, Hassan

I made dynamic mesh version of pisoCentralFoam - pisoCentralDyMFoam and it is available for downloading in git.

Can you update OpenFOAM-dev version?
mkraposhin is offline   Reply With Quote

Old   August 24, 2015, 17:58
Default
  #16
Senior Member
 
Hassan Kassem
Join Date: May 2010
Location: UK
Posts: 149
Rep Power: 7
hk318i is on a distinguished road
Quote:
Originally Posted by mkraposhin View Post
Hello, Hassan

I made dynamic mesh version of pisoCentralFoam - pisoCentralDyMFoam and it is available for downloading in git.

Can you update OpenFOAM-dev version?
Thank you very much. I updated the OpenFOAM-dev version and the tutorials as well.
I run a quick case for movingCone case for OpenFOAM-2.4 vs OpenFOAM-dev and I got the same results.

Regarding the turbulence modelling, Does the solver use the compressible models for the whole range?
hk318i is offline   Reply With Quote

Old   August 25, 2015, 07:00
Default
  #17
Senior Member
 
Hassan Kassem
Join Date: May 2010
Location: UK
Posts: 149
Rep Power: 7
hk318i is on a distinguished road
I tested the movingCone case for compressible flow using OpenFOAM-dev tutorial for rhoCentralDyMFoam. The results have an excellent agreement.
Attached Images
File Type: jpg Mach.jpg (28.3 KB, 62 views)
hk318i is offline   Reply With Quote

Old   August 25, 2015, 08:16
Default
  #18
Senior Member
 
mkraposhin's Avatar
 
Matvey Kraposhin
Join Date: Mar 2009
Location: Moscow, Russian Federation
Posts: 321
Rep Power: 10
mkraposhin is on a distinguished road
Quote:
Originally Posted by hk318i View Post

Regarding the turbulence modelling, Does the solver use the compressible models for the whole range?
Solver uses compressible turbulence models for the whole range as it done in rhoPimpleFoam and other compressible PISO/SIMPLE/PIMPLE OpenFOAM solvers

My next step is to port solver to openfoam-extend version 3.0 and 3.1

Also, i uploaded small video in youtube with results of of forwardStep test case
mkraposhin is offline   Reply With Quote

Old   August 25, 2015, 08:26
Default
  #19
Senior Member
 
Join Date: Oct 2013
Posts: 331
Rep Power: 6
chriss85 is on a distinguished road
I have implemented your solver in my application and I'm seeing nice speedups compared to the sonicFoam solver. Results are a bit different, but atleast in the simple shockTube case your solver shows better results, so I'm more inclined to trust your solver

I'm still somewhat wondering about the definition of the speed of sound:
In OpenFOAM (namely rhoCentralFoam and your solver) the following formula is used:
c = sqrt(Cp/Cv*psi) = sqrt(Cp/Cv*p/rho)
According to http://citeseerx.ist.psu.edu/viewdoc...0.1.1.215.7150 , this is equal to the frozen speed of sound, meaning the speed of sound of a fluid that doesn't change its chemical composition when pressure/temperature change. In my case I'm assuming a fluid in local thermodynamic equilibrium, meaning that chemical composition is a function of pressure and temperature. According to the paper above, speed of sound should then be defined as
C = sqrt(Cp/Cv/(drho/dp)) (with constant T in the derivative)

Do you agree with this modification for my case?
Depending on the kind of flow, the involved reaction rates, the occuring time-, temperature- und pressure- scales, both of these extremes may apply, or if one is unlucky, something in between. I'm not sure how the cases which lie in between may be treated properly.
chriss85 is offline   Reply With Quote

Old   August 25, 2015, 08:47
Default
  #20
Senior Member
 
mkraposhin's Avatar
 
Matvey Kraposhin
Join Date: Mar 2009
Location: Moscow, Russian Federation
Posts: 321
Rep Power: 10
mkraposhin is on a distinguished road
Quote:
Originally Posted by chriss85 View Post
I have implemented your solver in my application and I'm seeing nice speedups compared to the sonicFoam solver. Results are a bit different, but atleast in the simple shockTube case your solver shows better results, so I'm more inclined to trust your solver

I'm still somewhat wondering about the definition of the speed of sound:
In OpenFOAM (namely rhoCentralFoam and your solver) the following formula is used:
c = sqrt(Cp/Cv*psi) = sqrt(Cp/Cv*p/rho)
According to http://citeseerx.ist.psu.edu/viewdoc...0.1.1.215.7150 , this is equal to the frozen speed of sound, meaning the speed of sound of a fluid that doesn't change its chemical composition when pressure/temperature change. In my case I'm assuming a fluid in local thermodynamic equilibrium, meaning that chemical composition is a function of pressure and temperature. According to the paper above, speed of sound should then be defined as
C = sqrt(Cp/Cv/(drho/dp)) (with constant T in the derivative)

Do you agree with this modification for my case?
Depending on the kind of flow, the involved reaction rates, the occuring time-, temperature- und pressure- scales, both of these extremes may apply, or if one is unlucky, something in between. I'm not sure how the cases which lie in between may be treated properly.
C = sqrt(Cp/Cv/(drho/dp)) becomes sqrt(Cp/Cv/psi) = sqrt(Cp/Cv*p/rho) in case of perfect gas EOS. Where psi = (drho / dp) = 1/(R*T) and rho = psi*p
This formulation is true for all compressible implicit OF solvers. If you will change EOS, for example - rho = A*(p/p0)^k, then you have to change equation for pressure. for example - see sonicLiquidFoam
mkraposhin is offline   Reply With Quote

Reply

Thread Tools
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 On
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 23 February 2, 2016 15:33
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 15:15
High Mach number solver error Luke CFX 3 January 31, 2007 23: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 06:17.