[mkraposhin]

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]

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?

[mkraposhin]

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

[maHein]

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

[mkraposhin]

Hi, thank you for your opinion

Originally Posted by maHein

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

Originally Posted by maHein

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

[chriss85]

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.

[mkraposhin]

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

[chriss85]

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: http://www.cfd-online.com/Forums/ope...-equation.html .
Should the equations be implemented like the equation for rho itself, or is the default scalar transport equation sufficient?

[hk318i]

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?

[mkraposhin]

Originally Posted by hk318i

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

[hk318i]

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?

[mkraposhin]

Originally Posted by chriss85

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: http://www.cfd-online.com/Forums/ope...-equation.html .
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

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

[mkraposhin]

Originally Posted by hk318i

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

[hk318i]

Originally Posted by mkraposhin

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

[mkraposhin]

Hello, Hassan

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

Can you update OpenFOAM-dev version?

[hk318i]

Originally Posted by mkraposhin

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]

I tested the movingCone case for compressible flow using OpenFOAM-dev tutorial for rhoCentralDyMFoam. The results have an excellent agreement.

[mkraposhin]

Originally Posted by hk318i

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

[chriss85]

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.

[mkraposhin]

Originally Posted by chriss85

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