[nakul]

Hi,

I have modifed the rhoCentralFOAM solver for viscous supersonic flows to incorporate combustion and mixing. This solver is capable of simulating supersonic combustion. The combustion code has been taken from reactingFOAM. I have named it "rhoCentralSpeciesFOAM".

I have done a test case simulation with this and it works. But there are still few issues with the solver. The solver is not very stable. It is still quite sensitive to test case setup. Although the solver works but the results are not very accurate as the chemical reactions become more complex.

I request the OF community to have look at this solver. Any suggestions regarding any modifications are highly appreciated. Steps needed to integerate it with your OF have been specified in the following post.

-Nakul

[holger_marschall]

Hi Nakul,

I encourage you to found a working group and a sub-project on this topic on the Extend-Project portal (www.extend-project.de). There you'll find all infrastructure to organize the development. E.g.

• group discussion board (+group mailing)
• group file sharing
• project management (task management, messaging)

best regards,
Holger

[ranjansm]

Hi Nakul

I was wondering if you have developed the supersonic combustion solver any further, or put it in the openfoam-extend or anything like that ?

Thanks
Ranjan

[ranjansm]

Originally Posted by nakul

Hi,

I have modifed the rhoCentralFOAM solver for viscous supersonic flows to incorporate combustion and mixing. This solver is capable of simulating supersonic combustion. The combustion code has been taken from reactingFOAM. I have named it "rhoCentralSpeciesFOAM".

I have done a test case simulation with this and it works. But there are still few issues with the solver.

-Nakul

Nakul,

Do you happen to have the test case that you ran for this solver at all ?
Thanks
Ranjan

[nakul]

Hi,

No I didn't find time to work on it any further but I am interested in its further development.

I have run a test case using this solver. You may tell me what you expect to do with the solver and I would see if I could help with that.

-Nakul

[ranjansm]

Hi Nakul

I may have some dedicated time to work with you and get this rhoCentralSpeciesFoam in order.

But if you can give me the test case that you used, it will help me greatly.

I compiled your solver with OF-1.7.1 and ran a case, which I have, but I am sure my BCs are a big problem,

Quickly go temperature out of bounds etc.

We can work together, and as you have been suggested already, that it may be a good idea to put it in a open repository so multiple users can suggest changes etc.

Let me know!!
Thanks
Ranjan

[nakul]

Hi Ranjan

Please tell me your mail-id and I will mail you the test case by Saturday. We can discuss its setup after that. But I would advice that if you have dedicated time to work on supersonic combustion, then create a test case on your own.

Besides I don't have all the validation data for my test case so it would be a good idea if you prepare a case with large amount of validation data readily available. We can then discuss how to solve that problem. Also you may tell me what exactly are you planning to study in supersonic combustion because there are quite a few problems that one can solve in this area.

-Nakul

[qinfei]

Hi Nakul and ranjansm

I interest in simulation of combustion in high speed flow, I think rhoCentralSpeciesFoam will be very helpful to me. Could you give me the test case that you used, it will help me greatly. Then I will post the results of the simulation, and may develop it further. Thanks in advance.

[ranjansm]

Hi qinfei,
I am working iwth Nakul and one more engineer here.
I think we have ironed out most of the issues. There are some good validation cases available too.
As soon as things look good enough to getting out, I will let you know
Thanks!
Ranjan

Originally Posted by qinfei

Hi Nakul and ranjansm

I interest in simulation of combustion in high speed flow, I think rhoCentralSpeciesFoam will be very helpful to me. Could you give me the test case that you used, it will help me greatly. Then I will post the results of the simulation, and may develop it further. Thanks in advance.

[qinfei]

Hi ranjansm,

Thank your reply, I am waiting for your results, good luck!

[ChrisA]

Hi Nakul, Ranjan,

Your work sounds like the exact thing I'm working on over here, implementation of combustion/reactions from reactionFoam in our modified rhocentralfoam solver. Where are you guys at with the solver? Would it be possible to get one of your test cases?

[hg2lf]

When you calculate the sound speed c,you do like this:
c=sqrt(thermo.Cp()/(hermo.Cp()-R)*rPsi),R is the commom specific heat : 8.314j/mol.K,You should do like this:
c=sqrt(thermo.Cp()/(thermo.Cp()-rPsi/T)*rPsi).
Best regards.

[hg2lf]

when I modified the calculation of sound speed c, I use rhoCentralSpeciesFoam to simulate a cold flow in a laval nozzle,but the temperature field seems to be wrong. H =hs+1/2*U^2 does not keep constant. Somebody who research more deeply into the solver?

[tatu]

Originally Posted by hg2lf

when I modified the calculation of sound speed c, I use rhoCentralSpeciesFoam to simulate a cold flow in a laval nozzle,but the temperature field seems to be wrong. H =hs+1/2*U^2 does not keep constant. Somebody who research more deeply into the solver?

Nakul's phiHp is equal to the one with internal energy, i.e. it is

    surfaceScalarField phiEp
        (
            aphiv_pos*(rho_pos*(hs_pos + 0.5*magSqr(U_pos)) + p_pos)
          + aphiv_neg*(rho_neg*(hs_neg + 0.5*magSqr(U_neg)) + p_neg)
          + aSf*p_pos - aSf*p_neg
        );

Shouldn't this equation look different if we replace the pressure term with internal energy (rho*e = rho*hs - p). Could removing the first two pressure terms work? In other words:

    surfaceScalarField phiHp
        (
            aphiv_pos*(rho_pos*(hs_pos + 0.5*magSqr(U_pos)) )
          + aphiv_neg*(rho_neg*(hs_neg + 0.5*magSqr(U_neg)) )
          + aSf*p_pos - aSf*p_neg
        );

This seems to bring my underestimated temperature field closer to the value it should be in my case (according to reactingFoam and experiments).

There's also something fishy going on with the diffusive energy corrector equation. I get very promising temperature values if I comment out the last two terms:

fvm::ddt(rho, hs) - fvc::ddt(rho, hs)
 - fvm::laplacian(turbulence->alphaEff(), hs) // alphaEff = alpha + alphat 
//+ fvc::laplacian(turbulence->alpha(), hs)  
//- fvc::laplacian(k, T)   // turbulent + laminar energy diffusion via T

I guess the last two terms are there to drive energy and temperature consistent, i.e. hs --> CpT + p/rho + magSqr(U). I still don't fully understand the two last terms, e.g. why k is laminar + turbulent, but alpha is laminar?



Tatu

[ChrisA]

It's worth noting that the thermal conductivity (k) here is calculated using a prandlt number equal to 1 (or constant, depending on if they fixed the bug where it was hard coded equal to 1) this could be a source of discrepancy, especially in the boundary layer. Anyway, sutherland is able to calculate k (kappa in sutherland), giving it proper variability with respect to temperature. If you're using a different viscosity model, it should have an equation to calculate k. There just need to be functions implemented to calculate it for the domain and pass it up to the solver (you can copy the way Cp is written for the relevant libraries).

Also, with regards to your question, I don't believe k has a turbulent value as that term governs the energy transfer due to conduction, which shouldn't change as a result of turbulent flow (directly anyway).

The alpha bits I'm a bit more foggy on, my focus case is laminar so I haven't checked out the turbulence much.

It's been awhile since I've checked this stuff in the libraries, feel free to correct anything you think might be wrong.

[hg2lf]

Hi,Tatu,

Thanks for your reply. But don't you think the phiHp has nothing to do with p. I think the phiHp should be like this:

surfaceScalarField phiHp
(
aphiv_pos*(rho_pos*(hs_pos + 0.5*magSqr(U_pos)) )
+ aphiv_neg*(rho_neg*(hs_neg + 0.5*magSqr(U_neg)) )
);
so ,what do you think of this?

[tatu]

Originally Posted by hg2lf

Hi,Tatu,

Thanks for your reply. But don't you think the phiHp has nothing to do with p. I think the phiHp should be like this:

surfaceScalarField phiHp
(
aphiv_pos*(rho_pos*(hs_pos + 0.5*magSqr(U_pos)) )
+ aphiv_neg*(rho_neg*(hs_neg + 0.5*magSqr(U_neg)) )
);
so ,what do you think of this?

I'm not completely sure what is the meaning of the term aSf*p_pos - aSf*p_neg (anyone?). It could be there to capture "propagation of the stress" through cells due to compressibility effects. Anyway, I think it should be left in the equation. Can anyone confirm this?

Tatu

[hg2lf]

When you compare the equation in the solver and the governing equation,you could find that the phiHp should not inclue p. Maybe the term aSf*p_pos - aSf*p_neg has something to do with the scheme of rho central upwind ?

[tatu]

Originally Posted by hg2lf

When you compare the equation in the solver and the governing equation,you could find that the phiHp should not inclue p. Maybe the term aSf*p_pos - aSf*p_neg has something to do with the scheme of rho central upwind ?

Yep, that's what I meant.

Tatu

[sharonyue]

`aSf*p_pos - aSf*p_neg` is not related to the scheme. In the following code, p is simply advected by 0.5*(phi+p+ + phi-p-) instead of 0.5*(phi(KT scheme)p+ + phi(KT scheme)p-). Written in this way, they offset with each.

Oh you can think it in this way, if you neglect `aSf*p_pos - aSf*p_neg`, KT scheme would be used to discretize \nabla\cdot(Up), but it should not.

See Eqn 27 in http://dyfluid.com/rhoCentralFoam.html

```
surfaceScalarField phiEp
(
aphiv_pos*(rho_pos*(hs_pos + 0.5*magSqr(U_pos)) + p_pos)
+ aphiv_neg*(rho_neg*(hs_neg + 0.5*magSqr(U_neg)) + p_neg)
+ aSf*p_pos - aSf*p_neg
);
```