CFD Online Logo CFD Online URL
www.cfd-online.com
[Sponsors]
Home > Forums > Software User Forums > OpenFOAM > OpenFOAM Programming & Development

Modifying the energy equation to remove mass without changing temperature

Register Blogs Members List Search Today's Posts Mark Forums Read

Like Tree2Likes
  • 1 Post By Benben
  • 1 Post By Benben

Reply
 
LinkBack Thread Tools Search this Thread Display Modes
Old   August 14, 2018, 12:58
Smile Modifying the energy equation to remove mass without changing temperature
  #1
Member
 
benoit favier
Join Date: Jun 2017
Posts: 64
Rep Power: 4
Benben is on a distinguished road
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.
Benben is offline   Reply With Quote

Old   August 14, 2018, 13:09
Default
  #2
Member
 
benoit favier
Join Date: Jun 2017
Posts: 64
Rep Power: 4
Benben is on a distinguished road
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.
Benben is offline   Reply With Quote

Old   August 16, 2018, 03:18
Default
  #3
Member
 
Atul Kumar
Join Date: Dec 2015
Location: National Centre for Combustion Research and Development
Posts: 46
Rep Power: 5
atulkjoy is on a distinguished road
Quote:
Originally Posted by Benben View Post
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
atulkjoy is offline   Reply With Quote

Old   August 16, 2018, 03:19
Default
  #4
Member
 
benoit favier
Join Date: Jun 2017
Posts: 64
Rep Power: 4
Benben is on a distinguished road
I do this at the begining of the pimple loop, before "UEqn.H".

Edit : my bad, temperature is still changing
Benben is offline   Reply With Quote

Old   August 16, 2018, 09:03
Default
  #5
Member
 
benoit favier
Join Date: Jun 2017
Posts: 64
Rep Power: 4
Benben is on a distinguished road
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 is offline   Reply With Quote

Old   August 17, 2018, 05:43
Default
  #6
Member
 
benoit favier
Join Date: Jun 2017
Posts: 64
Rep Power: 4
Benben is on a distinguished road
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 is offline   Reply With Quote

Old   August 17, 2018, 07:01
Default
  #7
Member
 
benoit favier
Join Date: Jun 2017
Posts: 64
Rep Power: 4
Benben is on a distinguished road
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.
Benben is offline   Reply With Quote

Old   July 27, 2019, 02:10
Default
  #8
Member
 
Niu
Join Date: Apr 2014
Posts: 55
Rep Power: 7
Z.Q. Niu is on a distinguished road
Hi benoit,
Have you find the correct way to code right energy source? Thanks!

Best Regards!
Z.Q. Niu is offline   Reply With Quote

Reply

Tags
energy equation, reactingfoam

Thread Tools Search this Thread
Search this Thread:

Advanced Search
Display Modes

Posting Rules
You may not post new threads
You may not post replies
You may not post attachments
You may not edit your posts

BB code is On
Smilies are On
[IMG] code is On
HTML code is Off
Trackbacks are Off
Pingbacks are On
Refbacks are On


Similar Threads
Thread Thread Starter Forum Replies Last Post
whats the cause of error? immortality OpenFOAM Running, Solving & CFD 12 July 11, 2019 18:33
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


All times are GMT -4. The time now is 20:31.