Modifying the energy equation to remove mass without changing temperature

Benben · August 14, 2018, 12:58

Hello,
I am currently modifying reactingFoam, on version 6.0 of OpenFoam.
I am trying to implement an oxydation reaction :
H2O(g) + Metal(s) -----> H2(g) + MetalOxyde(s)
which for the gas phase is simply :
H2O(g) -----> H2(g)

To do so, i declared a "volScalarField" reaction rate R, which is my reaction rate in the volume in mol/m^3/s.

Code:

volScalarField R
(
    IOobject
    (
        "R",
        runTime.timeName(),
        mesh,
        IOobject::NO_READ,//Dont read at the begining
        IOobject::AUTO_WRITE//Save every time step
    ),
    mesh,
    dimensionedScalar("R", dimensionSet(0,-3,-1,0,1,0,0), scalar(0.0))
);

I initialized it to positive values before the main loop of reactingFoam.
I also defined the molar masses of H2O and H2 as "dimensionedScalar" variables, in kg/mol :

Code:

dimensionedScalar M_H2("M_H2", dimensionSet(1,0,0,0,-1,0,0), composition.W(i_H2) * 1.E-3);
dimensionedScalar M_H2O("M_H2O", dimensionSet(1,0,0,0,-1,0,0), composition.W(i_H2O) * 1.E-3);

I modified my equation in "YEqn.H" to have formation of H2 and destruction of H2O :

Code:

//Defining the source term of mass of specie i, in kg/mol/s.
volScalarField SourceTerm = R * dimensionedScalar("Molar_mass",dimensionSet(1,0,0,0,-1,0,0),0.);//Copy the volScalarField "R", change its unit to kg/m^3/s, and set it to 0.
if (Yi.name()=="H2")
{
    SourceTerm = R * M_H2;
}
else if (Yi.name()=="H2O")
{
    SourceTerm = R * -M_H2O;
}

fvScalarMatrix YiEqn
(
    fvm::ddt(rho, Yi)
  + mvConvection->fvmDiv(phi, Yi)
  - fvm::laplacian(turbulence->muEff(), Yi)
 ==
    fvOptions(rho, Yi)
  + SourceTerm
);

And then modified my "RhoEqn.H" to remove the mass of oxygen atoms that go to solid phase (as H2O is heavier than H2):

Code:

//Defining the sink term of mass due to the chemical reaction. (H2O(g) + Zn(s) ---> H2(g) + ZnO(s) => mass of gas is reduced)
volScalarField SinkTerm_reaction = R * (M_H2 - M_H2O);

fvScalarMatrix rhoEqn
(
    fvm::ddt(rho)
  + fvc::div(phi)
  ==
    fvOptions(rho)
  + SinkTerm_reaction
);

But then when I tested it, there was an unwanted modification of the temperature inside of my test case. (I set up a test case between 6 walls, with no velocity at all).
I realized that it was because i didn't remove the enthalpy of the oxygen atoms from my system, in the energy equation.
I dont really know how i should proceed.
Anyone has any clue of how i should do this ?

Benben · August 14, 2018, 13:09

I tried to do this :

created a new "volScalarField" variable he_R which is a energy source due to the oxydation reaction in J/m^3/s :

Code:

//Enthalpie variation due to the reaction (in J/m^3/s) [J = kg.m^2/s^2 => J/m^3/s = kg.m^-1.s^-3]
volScalarField he_R
(
    IOobject
    (
        "he_R",
        runTime.timeName(),
        mesh,
        IOobject::NO_READ,//Dont read at the begining
        IOobject::NO_WRITE//Dont save
    ),
    mesh,
    dimensionedScalar("he_R", dimensionSet(1,-1,-3,0,0,0,0), scalar(0.0))
);

I calculated it using the total enthalpy of specie H2 and H2O :

Code:

he_R[Cell_ID] = R[Cell_ID] * (
                              composition.Ha(i_H2,p[Cell_ID],Temp)*M_H2
                              - composition.Ha(i_H2O,p[Cell_ID],T [Cell_ID])*M_H2O
                              )

And added this source term in my energy equation "EEqn.H" :

Code:

fvScalarMatrix EEqn
(
    fvm::ddt(rho, he) + mvConvection->fvmDiv(phi, he)
  + fvc::ddt(rho, K) + fvc::div(phi, K)
  + (
        he.name() == "e"
      ? fvc::div
        (
            fvc::absolute(phi/fvc::interpolate(rho), U),
            p,
            "div(phiv,p)"
        )
      : -dpdt
    )
  - fvm::laplacian(turbulence->alphaEff(), he)
 ==
    Qdot
  + fvOptions(rho, he)
  + he_R
);

But the temperature is still changing during the simulation

What am i doing wrong ?
Is it really the total enthalpy "composition.Ha(...)" that i have to considerate ? (there is also Hs the sensible enthalpy, and Hc the chemical enthalpy, which i dont realy know what they are)
Should i use the partial pressure of species to calculate the enthalpy instead of total pressure of the gas ?

atulkjoy · August 16, 2018, 03:18

Quote:

Originally Posted by Benben

I tried to do this :

created a new "volScalarField" variable he_R which is a energy source due to the oxydation reaction in J/m^3/s :

Code:

//Enthalpie variation due to the reaction (in J/m^3/s) [J = kg.m^2/s^2 => J/m^3/s = kg.m^-1.s^-3]
volScalarField he_R
(
    IOobject
    (
        "he_R",
        runTime.timeName(),
        mesh,
        IOobject::NO_READ,//Dont read at the begining
        IOobject::NO_WRITE//Dont save
    ),
    mesh,
    dimensionedScalar("he_R", dimensionSet(1,-1,-3,0,0,0,0), scalar(0.0))
);

I calculated it using the total enthalpy of specie H2 and H2O :

Code:

he_R[Cell_ID] = R[Cell_ID] * (
                              composition.Ha(i_H2,p[Cell_ID],Temp)*M_H2
                              - composition.Ha(i_H2O,p[Cell_ID],T [Cell_ID])*M_H2O
                              )

And added this source term in my energy equation "EEqn.H" :

Code:

fvScalarMatrix EEqn
(
    fvm::ddt(rho, he) + mvConvection->fvmDiv(phi, he)
  + fvc::ddt(rho, K) + fvc::div(phi, K)
  + (
        he.name() == "e"
      ? fvc::div
        (
            fvc::absolute(phi/fvc::interpolate(rho), U),
            p,
            "div(phiv,p)"
        )
      : -dpdt
    )
  - fvm::laplacian(turbulence->alphaEff(), he)
 ==
    Qdot
  + fvOptions(rho, he)
  + he_R
);

But the temperature is still changing during the simulation

What am i doing wrong ?
Is it really the total enthalpy "composition.Ha(...)" that i have to considerate ? (there is also Hs the sensible enthalpy, and Hc the chemical enthalpy, which i dont realy know what they are)
Should i use the partial pressure of species to calculate the enthalpy instead of total pressure of the gas ?

Are you updating fields like

forAll(T.ref(),cellID)
{

he_R[Cell_ID] = R[Cell_ID] * (
composition.Ha(i_H2,p[Cell_ID],Temp)*M_H2
- composition.Ha(i_H2O,p[Cell_ID],T

}
before ?? or after ??? you van use mole/mass fraction of species to check that

Benben · August 16, 2018, 03:19

I do this at the begining of the pimple loop, before "UEqn.H".

Edit : my bad, temperature is still changing

Benben · August 16, 2018, 09:03

I changed the total enthalpy for the sensible enthalpy, as that is what i selected in my case setting (constant/thermophysicalProperties.energy = sensibleEnthalpy).
I also added the kinetic energy that is removed by the reaction, even if there is no velocity.

Code:

he_R[Cell_ID] = R[Cell_ID] *
(
    composition.Hs(i_H2,p[Cell_ID],Temp)*M_H2
    - composition.Hs(i_H2O,p[Cell_ID],T [Cell_ID])*M_H2O
    + K[Cell_ID]*(M_H2 - M_H2O)
)

My temperature is still decreasing.
An interesting fact is that even thought my mass is supposed to decrease, and my temperature decrease, the pressure increase.

Benben · August 17, 2018, 05:43

So, actually, the reason why the pressure increase is probably due to the fact that reactingFoam uses a pressure equation at the end of the loop to correct the pressure (poisson pressure equation). I didnt modify this equation to take into account the mass variation of the systeme. I will look into it once i managed to have my temperature variation fixed.

I am a bit stuck here, as according to this :
How exactlly is enthalpy calculated
The hentalpy variation to take into account is indeed the sensible enthalpy, but my code still doesnt work.

I checked that the mass is indeed removedof the system, and it is (rho goes from 0.393963 Initially to 0.394123 in my volume).

I instrumented the code to display internal enthalpy in a cell, before and fater the first time step, as well as the theoretical enthalpy variation as i calculate it, and it clearly isnt the same thing :

Code:

debug (YEqn_RDMP.H) : Hs before = 529701
debug (YEqn_RDMP.H) : T before = 600
debug (EEqn_RDMP.H) : Delta Hs theoretical = 16.7804
debug (EEqn_RDMP.H) : Hs after = 529812
debug (EEqn_RDMP.H) : T after = 599.991

I even displayed every calculated term of the energy equation, and everything is of the order of 10^(-20) or less, except for my energy source term.

Code:

debug (EEqn_RDMP.H) : fvc::ddt(rho, K) = -1.34174e-24
debug (EEqn_RDMP.H) : fvc::div(phi, K) = 0
debug (EEqn_RDMP.H) : -dpdt = -0
debug (EEqn_RDMP.H) : Qdot = 0
debug (EEqn_RDMP.H) : he_R = 16.7804

So how the hell does my Hs change so much ??

Benben · August 17, 2018, 07:01

According to this :
https://cfd.direct/openfoam/energy-equation/#x1-60003
The energy equation in term of enthalpy that is derived for openfoam injects the mass conservation at some points in it.
So, as my mass conservation equation is different, this equation is wrong.

I will try to obtain the equivalent equation with a mass sink.

Z.Q. Niu · July 27, 2019, 02:10

Hi benoit,
Have you find the correct way to code right energy source? Thanks!

Best Regards!

August 14, 2018, 12:58	Modifying the energy equation to remove mass without changing temperature	#1
Benben Member benoit favier Join Date: Jun 2017 Posts: 64 Rep Power: 9	Hello, I am currently modifying reactingFoam, on version 6.0 of OpenFoam. I am trying to implement an oxydation reaction : H2O(g) + Metal(s) -----> H2(g) + MetalOxyde(s) which for the gas phase is simply : H2O(g) -----> H2(g) To do so, i declared a "volScalarField" reaction rate R, which is my reaction rate in the volume in mol/m^3/s. Code: volScalarField R ( IOobject ( "R", runTime.timeName(), mesh, IOobject::NO_READ,//Dont read at the begining IOobject::AUTO_WRITE//Save every time step ), mesh, dimensionedScalar("R", dimensionSet(0,-3,-1,0,1,0,0), scalar(0.0)) ); I initialized it to positive values before the main loop of reactingFoam. I also defined the molar masses of H2O and H2 as "dimensionedScalar" variables, in kg/mol : Code: dimensionedScalar M_H2("M_H2", dimensionSet(1,0,0,0,-1,0,0), composition.W(i_H2) * 1.E-3); dimensionedScalar M_H2O("M_H2O", dimensionSet(1,0,0,0,-1,0,0), composition.W(i_H2O) * 1.E-3); I modified my equation in "YEqn.H" to have formation of H2 and destruction of H2O : Code: //Defining the source term of mass of specie i, in kg/mol/s. volScalarField SourceTerm = R * dimensionedScalar("Molar_mass",dimensionSet(1,0,0,0,-1,0,0),0.);//Copy the volScalarField "R", change its unit to kg/m^3/s, and set it to 0. if (Yi.name()=="H2") { SourceTerm = R * M_H2; } else if (Yi.name()=="H2O") { SourceTerm = R * -M_H2O; } fvScalarMatrix YiEqn ( fvm::ddt(rho, Yi) + mvConvection->fvmDiv(phi, Yi) - fvm::laplacian(turbulence->muEff(), Yi) == fvOptions(rho, Yi) + SourceTerm ); And then modified my "RhoEqn.H" to remove the mass of oxygen atoms that go to solid phase (as H2O is heavier than H2): Code: //Defining the sink term of mass due to the chemical reaction. (H2O(g) + Zn(s) ---> H2(g) + ZnO(s) => mass of gas is reduced) volScalarField SinkTerm_reaction = R * (M_H2 - M_H2O); fvScalarMatrix rhoEqn ( fvm::ddt(rho) + fvc::div(phi) == fvOptions(rho) + SinkTerm_reaction ); But then when I tested it, there was an unwanted modification of the temperature inside of my test case. (I set up a test case between 6 walls, with no velocity at all). I realized that it was because i didn't remove the enthalpy of the oxygen atoms from my system, in the energy equation. I dont really know how i should proceed. Anyone has any clue of how i should do this ? atulkjoy likes this.

August 14, 2018, 13:09		#2
Benben Member benoit favier Join Date: Jun 2017 Posts: 64 Rep Power: 9	I tried to do this : created a new "volScalarField" variable he_R which is a energy source due to the oxydation reaction in J/m^3/s : Code: //Enthalpie variation due to the reaction (in J/m^3/s) [J = kg.m^2/s^2 => J/m^3/s = kg.m^-1.s^-3] volScalarField he_R ( IOobject ( "he_R", runTime.timeName(), mesh, IOobject::NO_READ,//Dont read at the begining IOobject::NO_WRITE//Dont save ), mesh, dimensionedScalar("he_R", dimensionSet(1,-1,-3,0,0,0,0), scalar(0.0)) ); I calculated it using the total enthalpy of specie H2 and H2O : Code: he_R[Cell_ID] = R[Cell_ID] * ( composition.Ha(i_H2,p[Cell_ID],Temp)M_H2 - composition.Ha(i_H2O,p[Cell_ID],T [Cell_ID])M_H2O ) And added this source term in my energy equation "EEqn.H" : Code: fvScalarMatrix EEqn ( fvm::ddt(rho, he) + mvConvection->fvmDiv(phi, he) + fvc::ddt(rho, K) + fvc::div(phi, K) + ( he.name() == "e" ? fvc::div ( fvc::absolute(phi/fvc::interpolate(rho), U), p, "div(phiv,p)" ) : -dpdt ) - fvm::laplacian(turbulence->alphaEff(), he) == Qdot + fvOptions(rho, he) + he_R ); But the temperature is still changing during the simulation What am i doing wrong ? Is it really the total enthalpy "composition.Ha(...)" that i have to considerate ? (there is also Hs the sensible enthalpy, and Hc the chemical enthalpy, which i dont realy know what they are) Should i use the partial pressure of species to calculate the enthalpy instead of total pressure of the gas ? atulkjoy likes this.

August 16, 2018, 09:03		#5
Benben Member benoit favier Join Date: Jun 2017 Posts: 64 Rep Power: 9	I changed the total enthalpy for the sensible enthalpy, as that is what i selected in my case setting (constant/thermophysicalProperties.energy = sensibleEnthalpy). I also added the kinetic energy that is removed by the reaction, even if there is no velocity. Code: he_R[Cell_ID] = R[Cell_ID] * ( composition.Hs(i_H2,p[Cell_ID],Temp)M_H2 - composition.Hs(i_H2O,p[Cell_ID],T [Cell_ID])M_H2O + K[Cell_ID]*(M_H2 - M_H2O) ) My temperature is still decreasing. An interesting fact is that even thought my mass is supposed to decrease, and my temperature decrease, the pressure increase.

August 17, 2018, 05:43		#6
Benben Member benoit favier Join Date: Jun 2017 Posts: 64 Rep Power: 9	So, actually, the reason why the pressure increase is probably due to the fact that reactingFoam uses a pressure equation at the end of the loop to correct the pressure (poisson pressure equation). I didnt modify this equation to take into account the mass variation of the systeme. I will look into it once i managed to have my temperature variation fixed. I am a bit stuck here, as according to this : How exactlly is enthalpy calculated The hentalpy variation to take into account is indeed the sensible enthalpy, but my code still doesnt work. I checked that the mass is indeed removedof the system, and it is (rho goes from 0.393963 Initially to 0.394123 in my volume). I instrumented the code to display internal enthalpy in a cell, before and fater the first time step, as well as the theoretical enthalpy variation as i calculate it, and it clearly isnt the same thing : Code: debug (YEqn_RDMP.H) : Hs before = 529701 debug (YEqn_RDMP.H) : T before = 600 debug (EEqn_RDMP.H) : Delta Hs theoretical = 16.7804 debug (EEqn_RDMP.H) : Hs after = 529812 debug (EEqn_RDMP.H) : T after = 599.991 I even displayed every calculated term of the energy equation, and everything is of the order of 10^(-20) or less, except for my energy source term. Code: debug (EEqn_RDMP.H) : fvc::ddt(rho, K) = -1.34174e-24 debug (EEqn_RDMP.H) : fvc::div(phi, K) = 0 debug (EEqn_RDMP.H) : -dpdt = -0 debug (EEqn_RDMP.H) : Qdot = 0 debug (EEqn_RDMP.H) : he_R = 16.7804 So how the hell does my Hs change so much ??

Similar Threads
Thread	Thread Starter	Forum	Replies	Last Post
whats the cause of error?	immortality	OpenFOAM Running, Solving & CFD	13	March 24, 2021 08:15
Inconsistencies in reading .dat file during run time in new injection model	Scram_1	OpenFOAM	0	March 23, 2018 23:29
Ansys CFX problem: unexpected very high temperatures in premix laminar combustion	faizan_habib7	CFX	4	February 1, 2016 18:00
Difficulty In Setting Boundary Conditions	Moinul Haque	CFX	4	November 25, 2014 18:30
Two-Phase Buoyant Flow Issue	Miguel Baritto	CFX	4	August 31, 2006 13:02

August 16, 2018, 03:19		#4
Benben Member benoit favier Join Date: Jun 2017 Posts: 64 Rep Power: 9	I do this at the begining of the pimple loop, before "UEqn.H". Edit : my bad, temperature is still changing

August 17, 2018, 07:01		#7
Benben Member benoit favier Join Date: Jun 2017 Posts: 64 Rep Power: 9	According to this : https://cfd.direct/openfoam/energy-equation/#x1-60003 The energy equation in term of enthalpy that is derived for openfoam injects the mass conservation at some points in it. So, as my mass conservation equation is different, this equation is wrong. I will try to obtain the equivalent equation with a mass sink.

July 27, 2019, 02:10		#8
Z.Q. Niu Member Niu Join Date: Apr 2014 Posts: 55 Rep Power: 12	Hi benoit, Have you find the correct way to code right energy source? Thanks! Best Regards!