[Sk Hossen Ali]

Dear All,

I am running a case with reactingTwoPhaseEulerFoam(rTPEF) solver to simulate the plug flow reactor condition on H2/Air combustion.

The model is a 1D constant area adiabatic duct of 5mm length consisting of an inlet at the left end and an outlet at the right end at the pressure of 1 atm filled with N2(g) initially.

Input::------------------
At the inlet, the unburnt stoichiometric H2/Air (Y_H2 = 0.02865, Y_O2= 0.23, Y_N2= 0.74135) mixture at 1000 K and 1m/sec enters.

Expected Output::------
I am expecting the gases leaving out at the outlet to be at adiabatic flame temperature and constant pressure equilibrium composition.

(Note::I have commented out all the diffusion terms from the EEqn, UEqn, and YiEqn to ensure the solvers solves for an ideal plug flow reactor model)

The results obtained in the above reactingTwoPhaseEulerFoam, is compared with the equilibrium calculation of NASA CEA and Cantera.
I have solved the same problem using reactingFoam solver under the same condition.

Results::-----------------
The reactingFoam simulation results were almost in line with NASA CEA and Cantera results which proves its solution to be accurate, so that the velocity and dynamic results of reactingFoam simulation can be used to compare the same from rTPEF simulation.

The temperature at the outlet is almost reaching the adiabatic flame temperature, T = 2665 K compared to the adiabatic flame temperature T_ad = 2696 K.
The composition at the outlet is also matching... (The details of the comparison is given in the attached file.)
The velocity at the outlet for the case of rTPEF is not matching with reactingFoam.


Findings and Warnings in rTPEF::-------------------
Although the mass, momentum and energy residuals have shown to be decreasing and very less at the final time steps results.
One of the major concerns to me is the density calculated for the reactingPhase in the rTPEF case.
I have output of rho object for both of the phases in the rTPEF simulation (note that the other phase is of no importance to us since its phase fraction throughout the domain is zero initially and maintained thereafter)

The density of individual species is being calculated correctly, while in the case where there is a mixture of species the density is being calculated incorrectly.
While I did similar things in reactingFoam it's able to calculate properly all density.

For example the density at the inlet at the stoichiometric condition is 0.254839 kg/m^3, but the density calculated in the rTPEF is 0.0710424.

I am providing the case file for both reactingFoam and rTPEF and a tabulated comparison of these with NASA CEA and Cantera for H2/Air and Ch4/Air.
https://github.com/alihossen1995/compPDF

Note::
To output rho in rTPEF and rF we need to change few things in creatFieldsRefs.H (like NoWrite to MustWrite) files of the respective folders.

Queries : -
I believe that the wrong density calculated for the reactingPhase in rTPEF is the main reason why velocity at the exit in the case of rF and rTPEF is not.
Since both rF and rTPEF are using the same thermophysicalProperties libraries what is the going wrong in calculating reactingPhase density in rTPEF?

[Sk Hossen Ali]

Dear All,

the above problem of the density calculation disparity was solved successfully by editing the code segment in

src/thermophysicalModels/reactionThermo/mixtures/multiComponentMixture/multiComponentMixture.H

template<class ThermoType>
const ThermoType& Foam::multiComponentMixture<ThermoType>::cellVolMixture
(
    const scalar p,
    const scalar T,
    const label celli
) const
{
//    scalar rhoInv = 0.0;
//    forAll(speciesData_, i)
//    {
//        rhoInv += Y_[i][celli]/speciesData_[i].rho(p, T);
//    }
//
//    mixtureVol_ =
//        Y_[0][celli]/speciesData_[0].rho(p, T)/rhoInv*speciesData_[0];
//
//    for (label n=1; n<Y_.size(); n++)
//    {
//        mixtureVol_ += 
// Y_[n][celli]/speciesData_[n].rho(p,T)/rhoInv*speciesData_[n];
//    }

    mixtureVol_ = Y_[0][celli] * speciesData_[0];

    for (label n=1; n<Y_.size(); n++)
    {
        mixtureVol_ += Y_[n][celli]*speciesData_[n];
    }
    return mixtureVol_;
}


template<class ThermoType>
const ThermoType& Foam::multiComponentMixture<ThermoType>::
patchFaceVolMixture
(
    const scalar p,
    const scalar T,
    const label patchi,
    const label facei
) const
{
//    scalar rhoInv = 0.0;
//    forAll(speciesData_, i)
//    {
//rhoInv += Y_[i].boundaryField()[patchi][facei]/speciesData_[i].rho(p,T);
//    }
//
//mixtureVol_ =Y_[0].boundaryField()[patchi][facei]/speciesData_[0].rho(p, 
//     T)/rhoInv * speciesData_[0];
//
//    for (label n=1; n<Y_.size(); n++)
//    {
//        mixtureVol_ +=
//            Y_[n].boundaryField()[patchi][facei]/speciesData_[n].rho(p,T)
//          / rhoInv*speciesData_[n];
//    }

    mixtureVol_ = Y_[0].boundaryField()[patchi][facei] * speciesData_[0];

    for (label n=1; n<Y_.size(); n++)
    {
        mixtureVol_ += Y_[n].boundaryField()[patchi][facei]*speciesData_[n];
    }

    return mixtureVol_;
}