[ulugbey]

Hi,
I have a inhouse solver which is based on Openfoam 2.1.
When the energy equation is solved according to internal energy, the solver can be compiled. But When the energy equation is solved according to enthalpy, the solver cannot be compiled.
Energy equation according to internal energy is as follows:

{
// total enthalpy equation:
    
        //volScalarField k("k", thermo.Cp()*turbulence->muEff()/PrT); // thermal conductivity
            //volScalarField divrhoHPFlux = fvc::div(Riemann.rhoHFlux()+Riemann.pFlux());
            volScalarField divrhoHFlux = fvc::div(Riemann.rhoEFlux());
            //volScalarField divPFlux = fvc::div(Riemann.pFlux());
      
      solve
      (
        fvm::ddt(rhoE)     
        + divrhoHFlux                
        - fvc::laplacian(turbulence->alphaEff(), hs) // thermal Diffusion
        //+ fvc::laplacian(turbulence->alphaEff(), hs) // thermal Diffusion
        //- fvc::laplacian(k,T) // thermal Conduction
        ==
        fvc::div(sigmaDotU)
        + combustion->Sh()  // Source for enthalpy equation in reactive flows
      );
        
        hs.dimensionedInternalField() = rhoE.dimensionedInternalField() / rho.dimensionedInternalField() - 0.5*magSqr(U.dimensionedInternalField())+p.dimensionedInternalField()/rho.dimensionedInternalField(); //- turbulence->k()().dimensionedInternalField();
        Info<< " Old hs  min...max   = " << min(hs).value() << " ... " << max(hs).value() << endl;    
        Info<< " Old T   min...max   = " << min(T).value() << " ... " << max(T).value() << endl;

        if(temperatureFix)
        {
                  scalar Tlow = 280.0;
            Info << "TEMPERATURE ACTIVE" << endl;
            volScalarField dummyT(T);  //
            forAll(dummyT,i) dummyT[i] = Tlow;    // a volScalarField of the correct size and dimension with value = Tlow            
            //volScalarField lowestE(e);
            volScalarField lowestHs(hs);
            thermo.hsFromT(lowestHs,dummyT); // hs is computed to result in a temperature of Tdummy for the local composition in each cell
            //thermo.eFromT(lowestE,dummyT); // e is computed to result in a temperature of Tdummy for the local composition in each cell
            //e = max(e,lowestE);
            hs = max(hs,lowestHs);
            rhoE = rho * (hs + 0.5*magSqr(U)-p/rho);
        }

            //hs=min(hs,hsMax);
            //hs=max(hs,hsMin);
            //rhoH = rho * (hs + 0.5*magSqr(U));
        
       //  Info<< " New hs  min...max   = " << min(hs).value() << " ... " << max(hs).value() << endl;
      //   Info<< " New T   min...max   = " << min(T).value() << " ... " << max(T).value() << endl;

            hs.correctBoundaryConditions();
        //thermo.correct();
            rhoE.boundaryField() = rho.boundaryField() * ( hs.boundaryField() + 0.5*magSqr(U.boundaryField())+p.boundaryField()/rho.boundaryField());
        //rhoH = rho*(hs + 0.5*magSqr(U));
   
}

Energy equation according to enthalpy is as follows:

{
// total enthalpy equation:
    
        //volScalarField k("k", thermo.Cp()*turbulence->muEff()/PrT); // thermal conductivity
            //volScalarField divrhoHPFlux = fvc::div(Riemann.rhoHFlux()+Riemann.pFlux());
            volScalarField divrhoHFlux = fvc::div(Riemann.rhoHFlux());
            volScalarField divPFlux = fvc::div(Riemann.pFlux());
      
      solve
      (
        fvm::ddt(rhoH)     
        + divrhoHFlux    
        + divPFlux            
        - fvc::laplacian(turbulence->alphaEff(), hs) // thermal Diffusion
        //+ fvc::laplacian(turbulence->alphaEff(), hs) // thermal Diffusion
        //- fvc::laplacian(k,T) // thermal Conduction
        - fvc::div(sigmaDotU)
        ==
         fvc::ddt(p)
        + combustion->Sh()  // Source for enthalpy equation in reactive flows
      );
        
        hs.dimensionedInternalField() = rhoH.dimensionedInternalField() / rho.dimensionedInternalField() - 0.5*magSqr(U.dimensionedInternalField()); //- turbulence->k()().dimensionedInternalField();
        Info<< " Old hs  min...max   = " << min(hs).value() << " ... " << max(hs).value() << endl;    
        Info<< " Old T   min...max   = " << min(T).value() << " ... " << max(T).value() << endl;

        if(temperatureFix)
        {
                  scalar Tlow = 280.0;
            Info << "TEMPERATURE ACTIVE" << endl;
            volScalarField dummyT(T);  //
            forAll(dummyT,i) dummyT[i] = Tlow;    // a volScalarField of the correct size and dimension with value = Tlow            
            //volScalarField lowestE(e);
            volScalarField lowestHs(hs);
            thermo.hsFromT(lowestHs,dummyT); // hs is computed to result in a temperature of Tdummy for the local composition in each cell
            //thermo.eFromT(lowestE,dummyT); // e is computed to result in a temperature of Tdummy for the local composition in each cell
            //e = max(e,lowestE);
            hs = max(hs,lowestHs);
            rhoH = rho * (hs + 0.5*magSqr(U));
        }

            //hs=min(hs,hsMax);
            //hs=max(hs,hsMin);
            //rhoH = rho * (hs + 0.5*magSqr(U));
        
       //  Info<< " New hs  min...max   = " << min(hs).value() << " ... " << max(hs).value() << endl;
      //   Info<< " New T   min...max   = " << min(T).value() << " ... " << max(T).value() << endl;

            hs.correctBoundaryConditions();
        //thermo.correct();
            rhoH.boundaryField() = rho.boundaryField() * ( hs.boundaryField() + 0.5*magSqr(U.boundaryField()));
        //rhoH = rho*(hs + 0.5*magSqr(U));
   
}

What could be a problem ? How can I be able to compile the solver when the energy equation is solved according to enthalpy?


Best regards