[rkc.cfd]

Hello,

I am trying to modify rhoReactingBuoyantFoam to simulate mixing of two gases with buoyancy effects. The gases are at constant temperature and density fluctuations arise due to mass fractions. I have taken out the reaction terms and also the EEeqn.H from
rhoReactingBuoyantFoam.


Then I calculate the density based on the species mass fraction at the location and feed the mixture density into thermo.rho().

/*
 * Calculate Average rho
 */
    volScalarField Yt(0.0*Y[0]);
    volScalarField Winv(0.0*Y[0]);
    for (label n=0; n<Y.size(); n++)
    {
        Yt += Y[n];
        Winv += Y[n] / thermo.composition().W(n);
    }
    //- Mean molecular weight [kg/kmol]
    volScalarField Wmix = 1.0/Winv;

    //- Specific gas constant [J/(kg K)] for this mixture
    volScalarField Rmix = specie::RR/Wmix;

    forAll(Y,i){
             for (label n=0; n<Y.size(); n++) {
                 rhoMix[i]=0;
             }
            }

    //Calculate rhoMix
        forAll(Y,i){
          for (label n=0; n<Y.size(); n++) {
                   rhoMix[i]=p[i]/Rmix[i]*T[i];
          }
         }
thermo.rho()=rhoMix;

I have placed the above equation just before the p equation in the main solver.
The solution diverges with the following error:

Program received signal SIGFPE, Arithmetic exception.
0x00007ffff24dbab9 in Foam::divide (res=..., f1=..., f2=...) at fields/Fields/scalarField/scalarField.C:114
114    BINARY_OPERATOR(scalar, scalar, scalar, /, divide)
(gdb) bt
#0  0x00007ffff24dbab9 in Foam::divide (res=..., f1=..., f2=...) at fields/Fields/scalarField/scalarField.C:114
#1  0x000000000048f178 in Foam::divide<Foam::fvsPatchField, Foam::surfaceMesh> (res=..., gf1=..., gf2=...) at /home/krishna/OpenFOAM/OpenFOAM-2.3.x/src/OpenFOAM/lnInclude/GeometricScalarField.C:114
#2  0x000000000047bfba in Foam::operator/<Foam::fvsPatchField, Foam::surfaceMesh> (tgf1=..., gf2=...) at /home/krishna/OpenFOAM/OpenFOAM-2.3.x/src/OpenFOAM/lnInclude/GeometricScalarField.C:114
#3  0x000000000046a740 in main (argc=1, argv=0x7fffffffd128) at pEqn.H:72

Is there a better way to address the problem of updating the rho. The case runs fine with non-modified rhoReactionBuoyantFoam.


Cheers