[mkraposhin]

Originally Posted by chriss85

Just a short question on this topic. Why are you not separating the flux in positive/negative parts in the continuity equation?

solve
    (
        fvm::ddt(rho) + fvc::div(phi)
    );

Because phi = phiPos + phiNeg
This means, that for continuity equation, formulated with explicit scheme, both of two versions are equals mathematically and numerically:
div(phi) = div(phiPos + phiNeg) = div(phiPos) + div(phiNeg)

[chriss85]

Thanks, that makes sense. I need to read up on the Kurganov-Tadmor schemes when I find some time.

[levvis]

Hi Matvey

I'm Lewis, a student from the UK studying Aerospace Engineering.

For my Final Year Project, I have been allocated "Development of an OpenFOAM code for Aerospace Engineering".

I have run some tutorial cases and got to grips with OpenFOAM and ParaView, using Ubuntu and the terminal window.

As your code is a newer hybrid solver, with scope for expansion, I am very interested in using your published code as a foundation. Therefore I will adapt it to my own need, the beauty of OpenFOAM.

I see you've already run cases using a cylinder, which is ideal for me, as my plan with your code is to run it with a cylinder within the geometry. And progress this shape into an aerofoil, for greater use in the Aviation Industry.


What are your opinions on use of your code to adapt and create this application? As I would like to be able to run the simulations at a sub-sonic velocity, and then supersonic.

Many thanks in advance

[mkraposhin]

Hi, Lewis

I'm not sure that i undestand your question correctly.

Because you are talking about

Development of an OpenFOAM code for Aerospace Engineering

,
or

 I am very interested in using your published code as a foundation

but you didn't mentioned which physical case you want to consider. So i can't understand for which case you want to extend pisoCentralFoam solver.

As for simulation of external flow over bodies, i can tell you that this solver was test for simple and complex geomtries for Mach humbers from 0.01 to 6 by me and my colleagues.

So, if you need to use it only for gas dynamics, then you don't need to change it. In this case, we need to create a separate thread on forum in "Solving" section for questions about setup of numerical case.

And if you are talking about extending this solver, we can proceed here, but you must supply more information on your physical problem

[levvis]

Hi Matvey

My question is, do you think there is an opportunity for me to use this for airflow over an aerofoil, at transonic speeds? Will it converge?

I suppose you have already answered that question, when you said you and your colleagues used it with simple and complex geometries.

I'd like to use this solver, just have to edit some of the code to show I know what I'm doing, and why I did it.


Thanks

Lewis

[mkraposhin]

Originally Posted by levvis

Hi Matvey

My question is, do you think there is an opportunity for me to use this for airflow over an aerofoil, at transonic speeds? Will it converge?

I suppose you have already answered that question, when you said you and your colleagues used it with simple and complex geometries.

I'd like to use this solver, just have to edit some of the code to show I know what I'm doing, and why I did it.


Thanks

Lewis

'

Hi, i think, yes, it will converge. We have plans to make it LTS, but not for this year.
Also, we make publication YSC 2015 procedia, so you can download pdf from here

PisoCentralFoam at researchgate

or from elsevier

doi:10.1016/j.procs.2015.11.007

Also, we made application of this scheme for multicomponent reacting flow

It will be nice if you will help us during your work to investigate range of application of this scheme (solver)

[mkraposhin]

Dear colleagues!

We made update of pisoCentralFoam solver for OpenFOAM 3.0.0. In new version LTS time scheme is included and pisoCentralFoam can be used for steady state simulations

[Sethu]

Dear Matvey Kraposhin,

Thank you very much for your pisoCentralFoam LTS version. Im wondering whether I can use your solver for application like steady state simulation of NASA rotor 37 is an isolated transonic axial compressor rotor where i need to account for rotational effect either in SRF or MRF (for rotor+stator). Please give me your comments on this.

I had modified rhoSimpleFoam as rhoSRFSimpleFoam for a single rotor case. But, when i executed the solver im getting error as follows,
[Maximum number of iterations exceeded

From function thermo<Thermo, Type>::T(scalar f, scalar T0, scalar (thermo<Thermo, Type>::*F)(const scalar) const, scalar (thermo<Thermo, Type>::*dFdT)(const scalar) const, scalar (thermo<Thermo, Type>::*limit)(const scalar) const) const
in file /home/nrcs-cae/OpenFOAM/OpenFOAM-2.3.x/src/thermophysicalModels/specie/lnInclude/thermoI.H at line 76.

FOAM aborting

#0 Foam::error:rintStack(Foam::Ostream&) at ??:?
#1 Foam::error::abort() at ??:?
#2 Foam::species::thermo<Foam::hConstThermo<Foam:er fectGas<Foam::specie> >, Foam::sensibleInternalEnergy>::TEs(double, double, double) const at ??:?
#3 Foam::hePsiThermo<Foam:siThermo, Foam:ureMixture<Foam::constTransport<Foam::speci es::thermo<Foam::hConstThermo<Foam:erfectGas<Foa m::specie> >, Foam::sensibleInternalEnergy> > > >::calculate() at ??:?
#4 Foam::hePsiThermo<Foam:siThermo, Foam:ureMixture<Foam::constTransport<Foam::speci es::thermo<Foam::hConstThermo<Foam:erfectGas<Foa m::specie> >, Foam::sensibleInternalEnergy> > > >::correct() at ??:?
#5
at ??:?
#6 __libc_start_main in "/lib/i386-linux-gnu/libc.so.6"
#7
at ??:?
Aborted (core dumped)


Any advice on selecting solvers for my steady state (SRF,MRF) turbo machinery compressible cases?


Regards,

Sethu

[mkraposhin]

Hi, Sethu

I added MRF capabilities to pisoCentralFoam, you can download new solver and tutorials from github.

As for rhoSRFSimpleFoam, you must create new thread, because this thread dedicated to pisoCentralFoam

[Bryan S]

Dear Dr. Kraposhin,
I have been using reactingCentralFoam to compute a nonreacting supersonic flow with foreign gas injection in a boundary layer. The mixtures are binary, but the foreign gas often has a molecular weight that is quite different from that of the free stream gas (nitrogen). I have a question about how species diffusion is handled. It seems to me that a unity Schmidt number assumption is made using the mixture-averaged viscosity, such that all species effectively have the same diffusivity. Am I correct in this understanding?

I also see that the equations invoke Fick's Law for diffusion, except I don't see a correction velocity as is typically done when mass fractions are implemented as opposed to mole fractions in Fick's Law. If I'm correct, do you have a sense of what the errors are associated with this both in the species transport and energy equations? Thanks.

[levvis]

Hi Matvey

I'm struggling to understand this error, all I have done is simply created a new geometry (a double-wedge aerofoil) in GMSH and meshed it. And then substituted this in for the cylinder geometry.

I've changed a few properties too, like rho = 1.225 and velocity to simulate conditions at Sea Level. And upon running the command I get a printstack error:

Mean and max Acoustic Courant Numbers = 0.0111359 0.102596
Courant Number mean: 0.0046786 max: 35.2339
Time = 0.00248

diagonal: Solving for rho, Initial residual = 0, Final residual = 0, No Iterations 0
#0 Foam::error:rintStack(Foam::Ostream&) at ??:?
#1 Foam::sigFpe::sigHandler(int) at ??:?
#2 ? in "/lib/x86_64-linux-gnu/libc.so.6"
#3 Foam::GaussSeidelSmoother::smooth(Foam::word const&, Foam::Field<double>&, Foam::lduMatrix const&, Foam::Field<double> const&, Foam::FieldField<Foam::Field, double> const&, Foam::UPtrList<Foam::lduInterfaceField const> const&, unsigned char, int) at ??:?
#4 Foam::GaussSeidelSmoother::smooth(Foam::Field<doub le>&, Foam::Field<double> const&, unsigned char, int) const at ??:?
#5 Foam::smoothSolver::solve(Foam::Field<double>&, Foam::Field<double> const&, unsigned char) const at ??:?
#6 Foam::fvMatrix<double>::solveSegregated(Foam::dict ionary const&) at ??:?
#7 Foam::fvMatrix<double>::solve(Foam::dictionary const&) at ??:?
#8 ? at ??:?
#9 ? at ??:?
#10 __libc_start_main in "/lib/x86_64-linux-gnu/libc.so.6"
#11 ? at ??:?
Floating point exception (core dumped)

i haven't changed anything else, so haven't even delved into the physics of the solver.

Can you, or anyone else, please advise?

Many thanks

[chriss85]

I suggest trying another linear equation solver. Matvey (or someone else from his group) has ported the BiCGStab solver which is more stable and crashes less often. Check out the gits, it's in one of them.

[levvis]

Originally Posted by chriss85

I suggest trying another linear equation solver. Matvey (or someone else from his group) has ported the BiCGStab solver which is more stable and crashes less often. Check out the gits, it's in one of them.

I shall try this too, but just trying to problem solve the pisoCentralFoam.
It has been edited ever so slightly and is now running without crashing, but the flow is wrong.
It's an MH45 aerofoil at an angle of attack of 10 degrees.

But when running it, the lower velocity just seems to grow out from the aerofoil.
Has anyone seen this before? Can anyone advise? Am I missing files or something?

Picture attached

EDIT: Problem solving files, the time step is very small and thus the end time is small, this is probably the reason for the simulation? Because the simulation is correct in respective to such a small time step?
Am I right in thinking so? I've changed it and have it running but obviously it's a long compute time, even with a higher time step.

Just thinking as the difference between the scale is very small, and looking at solvers like simpleFoam aerofoil case, if converted to a small time step, it shows a similar picture.

[jimteb]

Hi Matvey,

I saw from an earlier post that pisoCentralFoam is not bound to the acoustic Courant number, but rather the flow Courant number.

Does this mean that this solver is unsuitable for cases where you are interested in acoustic waves?

Thanks,

James

[mkraposhin]

Hi,

We tested this solver for several simple cases (monopole and diploe). It was found, that it can be used to reproduce acoustic waves if you will preserve next conditions:

1) Acoustics courant < 1 (Or better < 0.1)
2) Mesh resolution is > 20 cells per resolved wave

[mg.mithun]

Hi mkraposhin,
Have you tried this solver for cavitating flows?

[AlvarezEdo]

Hi Mkraposhin,

I'm modeling a shock tube with multiple sections of different gases. I've been an OpenFOAM front-end user for a while, but I have never done anything complicated like mixing two solvers. I'm using your solver to model the shock with turbulent boundary layers along the shock tube, and is working very well, but now I need to add the different immiscible stages of gases.

How should I go about doing so? Is there a specific solver (compressibleMultiphaseInterFoam, for example) you would recommend mixing into your solver?



Thanks for your help, and great job in developing pisoCentralFoam!

[strakey]

Matvey - I have been testing your reactingPimpleCentralFoam solver (OF4 version) for the purpose of modeling Rotating Detonation Engines. Doing some simple non-reacting simulations of hydrogen injected into air at subsonic Ma# (0.7ish) shows tremendous numerical diffusion compared to rhoReactingFoam. I'm using the standardMachToAcCourantRatio Kappa function and maxkappa is pretty much 1.0 everywhere and minKappa is around 0.5. If I force Kappa to 0 (in the code), the jets look very much like rhoReactingFoam, so it appears to be the central fluxes causing the diffusion. I was wondering why use Ma/CFL as the kappa function as opposed to just Ma?



Thanks for posting your code, it looks like a significant contribution to the community!


Pete

[mkraposhin]

Hi, strakey!

Numerical diffusion in hybrid KNP/PIMPLE solvers arises due to next approximation techniques:

1. TVD approximation of convective fluxes
2. KNP approximation (which is actually very similar to HLL)
3. Coupling between MULES approximation for mass fractions and energy equation
4. Euler approximation for temporal derivatives
5. Spatial mesh

Therefore, you can reduce numerical diffusion by tuning OpenFOAM settings in next way:

1. Try to use different limited schemes: upwind => Minmod(V) => vanLeer(V)
2. KNP (or HLL) adds excessive diffusion to contact discontinuities, but it makes solution less diffusive near shocks, therefore, it is not avoidable
3. Very important point is MULES-KNP coupling -- transport of enthalpy/energy in multicomponent flow must synchronized with transport of mass fractions. To do that, I tried to use high order coefficients (lambdaCoeffs) from MULES scheme in energy equation. However, this approach has a significant drawback: lambdaCoeffs is equal to 0 in the region where cells are occupied only by inert specie (Ydefault). Therefore, where Ydefault is equal to 1, there energy transport is upwinded. I tried to resolve this issue for OpenFOAM+ v2112: https://github.com/unicfdlab/hybridC...gitef-dev-2112 .
4. Also, you should try backward for ddt terms instead of Euler. This option should decrease numerical diffusion significantly.
5. Finally, use uniform Cartesian (quadrangle or cubic) meshes in the region of interest

Regarding your question about blending function kappa. I've tried several different ways - Ma/CCo, 1/CCo and others and found that Ma/CCo serves better for correct simulation. Considerations are simple: when Ma < 1, but CCo is also < 1 (actually it must be equal to Ma), then KNP scheme must be recovered, because it works well for subsonic flows with acoustic waves propagation. However, when M < 1, while CCo > 1, the standard PIMPLE method must be used. Finally, when M > 1, then regardless of CCo, the KNP method should work.

It would be cool if you find the way to improve this blending and share it with the community.

And just two final notes:
1) I strongly recommend to use latest versions . Since I don't have too much to support all OpenFOAM branches, I fix bugs only in the latest one.
2) I've published another solver for highspeed multicomponent flows (reactingLagrangianQGDFoam: https://github.com/unicfdlab/QGDsolv...m/gasDiffusion), which can be very good for freestream simulations with shocks and shears.

Look at the next comparison of solvers for overexpanded jet.

rhoCentralFoam:
RCF-centerline_ds.png

pimpleCentralFoam (but with Euler ddt scheme):
PCF-centerline_ds.png

QGDFoam
QGD-centerline_ds.png


Lines with different colores on these figures correspond to different mesh resolutions.

rhoCentralFoam implements KNP scheme, which is from the Gold fund of numerical algorithms. pimpleCentralFoam empirically was found to be more robust for complex 3D geometries. And QGDFoam is superior in resolving complex flows. but in relatively simple geometries.

[strakey]

Matvey - Thanks for the reply. I think I understand this better now. I downloaded the OF6 version and compiled it for OF7 with some minor changes. I also downloaded the Yeqn.H and heqn.H from the 2112 branch and this also compiled fine under OF7. The MULES modifications you made here did have a big impact on diffusion of temperature. I also have it outputting the hLambdaCoeff as a volScalarField with an fvc::average() procedure and the coefficients are near 1 now as expected. I still, however see too much diffusion in the species, much more than energy. From the Yeqn.H code, it looks like all species equations are being modified by the same set of lambdaCoeffs, which are 0 regardless of what I set my inertSpecies to. I'm not sure how a slicedSurfaceScalarField works though. I'm going to insert some more output fields to see what exactly is being computed for lambdaCoeffs for each of the species (4 in my case) but it seems as though lambdaCoeffs should be different for each species field?



Update: It looks like the lambdaCoeffs (allFacesLambda) being used to modify the transport equations in Yeqn.H is that for the last species in the loop calling the MULES limiter,

forAll(Y, iCmpt)
    {
        volScalarField& Yi = Y[iCmpt];
        if ( Yi.name() != inertSpecie )
        {
            Info << "Limiting for: " << Yi.name() << endl;
            surfaceScalarField& rhoPhiYCorr = phiY[iCmpt];
        
            mulesWithDiffusionImplicitLimiter
            (
                rho,
                Yi,
                phi_own,
                phi_nei,
                allFacesLambda,
                rhoPhiYCorr,
                diffusiveFlux[iCmpt],
                mDCf[iCmpt],
                SuSp[iCmpt]
            );
        
        if (Yi.name() == "CH4")  //pete 2-28-2022
        {
        lambdavolH2O=fvc::average(lambdaCoeffs);
        }
        }
    }

it just keeps getting overwritten.



For my case, the lambdaCoeff for the last species is 0 everywhere, so 0 is being used for all of the species fields everywhere. I think hLambdaCoeff is just left at its initial value (1.0).



I also tried backward, vanLeer and other things to minimize diffusion but these all have a small effect.


Lastly, I don't know if you have seen this work (https://www.appliedccm.com/wp-conten...M-2019-DWS.pdf), but this seems similar to what you are doing here with your hybrid solver.


Thanks
Pete