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

How to add Source term (2) for PYROLYSIS - reactingOneDim

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

Like Tree2Likes
  • 1 Post By Kummi
  • 1 Post By atulkjoy

Reply
 
LinkBack Thread Tools Search this Thread Display Modes
Old   September 11, 2019, 00:10
Question How to add Source term (2) for PYROLYSIS - reactingOneDim
  #1
Senior Member
 
Kumaresh
Join Date: Oct 2016
Posts: 208
Rep Power: 6
Kummi is on a distinguished road
Send a message via Yahoo to Kummi
Hello Foamers,
My topic of research is to build the 1D mathematical model for simulating the coal pyrolysis in OpenFOAM. Based on fireFOAM (OF211), I developed own solver for only pyrolysis. The code construction for pyrolysis was built under reactingOneDim (~/OpenFOAM/OpenFOAM-2.1.1/src/regionModels/pyrolysisModels/reactingOneDim).

By default, TEqn. in reactingOneDim.C
Quote:
fvScalarMatrix TEqn
(
fvm::ddt(rhoCp, T_) //Unsteady term
- fvm::laplacian(K_, T_) //Laplacian term
==
chemistrySh_ //Heat loss rate - Source term
// + fvc::div(phiQr)
+ fvc::div(phiGas) //Source term 1
+ r*cp*dT/dx //Source term 2
);
Quote:
ENERGY EQU: rho*Cp*(dT/dt) = del/delx[K*dT/dx ] + rho*dQ/dT + r*cp*dT/dx (W/m3)
Unsteady conduction = Diffusion term + SOURCE TERM 1 + SOURCE TERM 2
My query is all about the SOURCE TERM 2. Concerning it, my query is all about.

SOURCE TERM 1 (Pyrolysis) = heat loss due to pyrolysis solved by Arrhenius-like degradation chemistry (SOLVED this source term based on FireFOAM 1D pyrolysis model code- NO ISSUES HERE)
SOURCE TERM 2 (Evaporation) = phase change - moisture to vapor (moisture embedded in the wet coal).
I have figured out that the source terms can be included inside the pyrolysis solver under "reactingOneDimCoeffs"in case file (~/OpenFOAM/OpenFOAM-2.1.1/tutorials/combustion/fireFoam/les/oppositeBurningPanels/constant/pyrolysisZones)
Quote:
pyrolysis
{
active true;
pyrolysisModel reactingOneDim;
regionName panelRegion;
reactingOneDimCoeffs
{
filmCoupled false;
gasHSource false; //Energy source term due to pyrolysis gas (found in higher OF version)
QrHSource false; //Energy source term due in depht radiation (found in higher OF version)
radFluxName Qr;
minimumDelta 1e-8;
reactionDeltaMin 1e-8; // if ("filmDelta0" > reactionDeltaMin_)
moveMesh false;
}
infoOutput true;
}
Quote:
CONDITION for SOURCE TERM 2: T= 100 ==> alpha (moisture content) =0

Based on it, mass and heat balances are calculated at the boiling plane (surface) as,
Mass balance (r) = rho * alpha * Velocity [kg/m2.s]
Heat balance (-K dT/dx) = (r) * latent heat [W/m2]
Moisture Evaporation is calculated at the surface (Boiling Plane). The calculated mass balance (r) is multiplied with cp and dT/dx, which will be taken as SOURCE TERM 2 carrying the unit of W/m3.

REF - Clear details here in this MANUSCRIPT ~https://sci-hub.tw/10.1016/0016-2361(83)90225-9 (ATKINSON, B., & MERRICK, D. (1983). Mathematical models of the thermal decomposition of coal4. Heat transfer and temperature profiles in a coke-oven charge. Fuel, 62(5), 553561)

Have anyone come across such problems in OpenFOAM? Anyone tried coding SOURCE TERM for pyrolysis under reactingOneDim.C file ?

Correct me if I'm wrong anywhere please. I'm trying my best to add source term under reactingOneDimCoeffs for pyrolysis solver.
Please share your ideas, it will be highly helpful.
Thank you.
atulkjoy likes this.
Kummi is offline   Reply With Quote

Old   September 15, 2019, 15:57
Default
  #2
Member
 
Atul Kumar
Join Date: Dec 2015
Location: National Centre for Combustion Research and Development
Posts: 48
Rep Power: 6
atulkjoy is on a distinguished road
fvScalarMatrix TEqn
(
fvm::ddt(rhoCp, T_) //Unsteady term
- fvm::laplacian(K_, T_) //Laplacian term
==
chemistrySh_ //Heat loss rate - Source term moss loss rate*LatentHet
// + fvc::div(phiQr) // Radiation diffusion inside solid
+ fvc::div(phiGas) //Source term 1 // Cpg*(Ts-Tg) // sp. heating due to gas temp
+ r*cp*dT/dx //Source term 2 // CpsTs // sp capicity of solid
);

in setting of control/pyrolysisZoneDict term 1 fvc::div(phiGas) and term 2 fvc::div(phiQr) are optional.


fvc::div(phiQr)+fvc::div(phiGas)+chemistrySh_ = m'''Hp, where Hp is total heat of pyrolysis or (CpsTs - L + Cpg*(Tg-Ts))
Kummi likes this.
atulkjoy is offline   Reply With Quote

Old   September 15, 2019, 23:30
Default
  #3
Senior Member
 
Kumaresh
Join Date: Oct 2016
Posts: 208
Rep Power: 6
Kummi is on a distinguished road
Send a message via Yahoo to Kummi
Dear Atul Kumar,
I have no words enough to Thank you. ^^
Thank you for your time and response. It helps me a lot to proceed further.
Kummi is offline   Reply With Quote

Old   September 17, 2019, 03:57
Question SOURCE TERM 2 - Compilation ERROR
  #4
Senior Member
 
Kumaresh
Join Date: Oct 2016
Posts: 208
Rep Power: 6
Kummi is on a distinguished road
Send a message via Yahoo to Kummi
Dear Atul Kumar,
I have created the structure for Source TERM 2 shown below, which is created under reactingOneDim.C for the following condition:
Quote:
CONDITION for SOURCE TERM 2: T= 100 ==> alpha (moisture content) =0

Based on it, mass and heat balances are calculated at the boiling plane (surface) as,
Mass balance (r) = rho * alpha * Velocity [kg/m2.s]
Heat balance (-K dT/dx) = (r) * latent heat [W/m2]
Quote:
reactingOneDim.C - STRUCTURE

void reactingOneDim::updateMoistureEvaporation()
{
moisEva_ == dimensionedScalar("zero", moisEva_.dimensions(), 0.0);
scalar Tsat = 373; //AT Tsat=100deg, moisture evaporates
dimensionedScalar Enthalp("Enthalp",dimEnergy/dimMass,2.257e6); // Enthalp for coal decomposition [J/kg]

forAll(intCoupledPatchIDs_, i) // Propagate evaporated moisture through 1-D regions
{
const label patchI = intCoupledPatchIDs_[i]; // Calculation over specified PATCH --> patchI (labelID) creating patchID
const scalarField& moisEvap = moisEva_.boundaryField()[patchI]; //boundaryField()[patchI] - "values on the patch surface" with boundary patch variable called "moisEvap" - script to access boundary field datas
scalarField& alphap = alpha_.boundaryField()[patchI]; //moisture fraction declaration
scalarField& Vxp = Vx_.boundaryField()[patchI]; //velocity declaration
forAll(moisEvap, faceI) // "addressed moisEvap as boundary patch - implemented for multiple faces" with faceID as faceI
{
const labelList& cells = boundaryFaceCells_[faceI]; //boundaryFaceCells_[faceI] - EVERY "SINGLE CELL" IS ACCESSED IN "ALL FACEs"
scalar evapRate = 0.0; //mass flux
scalar heatFlux = 0.0; //heat flux
// scalar alpha = 0.1; // 10% moisture content
// scalar Vxp
forAll(cells, k)
{
const label cellI = cells[k]; //patchID = cellI EVERY "SINGLE CELL" IS ACCESSED --> to calculate
//forAll(mesh.cells(), celli)
if (T[cellI] = Tsat) //Tsat.value() - line 173 - LOOP starts
{
alphap[cellI] = scalar(0.0); //Not sure about this alpha declaration
}

//volScalarField VX = V.component(vector::X) ~ volScalarField VX(0.0);
evapRate =rho_*alphap[cellI]*Vxp; //kg/m2.s - line 178 - calculation of mass flux
heatFlux +=evapRate*Enthalp.value(); //W/m2 - heat flux calculation
moisEva_[cellI] +=evapRate;
}}}}
After compiling, I got the following error
Quote:
reactingOneDim/reactingOneDim.C: In member function ‘void Foam::regionModels:yrolysisModels::reactingOneDi m::updateMoistureEvaporation()’:
reactingOneDim/reactingOneDim.C:173:17: error: invalid types ‘<unresolved overloaded function type>[const label {aka const int}]’ for array subscript
reactingOneDim/reactingOneDim.C:178:40: error: cannot convert ‘Foam::tmp<Foam::Field<double> >’ to ‘Foam::scalar {aka double}’ in assignment
reactingOneDim/reactingOneDim.dep:566: recipe for target 'Make/linux64Gcc47DPOpt/reactingOneDim.o' failed
make: *** [Make/linux64Gcc47DPOpt/reactingOneDim.o] Error 1
The error is all due to the line 178 and 173.
Line 173 - Temperature loop starts
Line 178 - Calculation of mass flux

I am not sure how the temperature loop can be applied inside such kind of above structure ?
Kindly check the attachments 1 and 2 for understanding about the phenomenon and share your ideas please.
Attached Images
File Type: jpg ATTACHMENT_1_Boiling.jpg (110.9 KB, 7 views)
File Type: png ATTACHMENT_2_Convection.png (153.7 KB, 8 views)

Last edited by Kummi; September 17, 2019 at 06:11.
Kummi is offline   Reply With Quote

Old   September 17, 2019, 10:10
Default Modification (1)
  #5
Senior Member
 
Kumaresh
Join Date: Oct 2016
Posts: 208
Rep Power: 6
Kummi is on a distinguished road
Send a message via Yahoo to Kummi
Dear Atul Kumar,
Thank you for directing. I made some corrections.

Quote:
forAll(cells, k)
{
const label cellI = cells[k];
//if (T[cellI] = Tsat)
//if (T.internalField[cellI] = Tsat) //T is a volScalarField
if (T.regionMesh[cellI] = Tsat) // line 173 - LOOP starts
{
alphap[cellI] = scalar(0.0);
}
//else alphap[cellI] = XX

//evapRate =rho_*alphap[cellI]*Vxp[faceI]; //kg/m2.s
evapRate =alphap[cellI]*Vxp[faceI]; // line 178 - calculation of mass flux

heatFlux +=evapRate*Enthalp.value(); //W/m2
moisEva_.boundaryField()[patchI][faceI] += evapRate; //boundaryField()[patchI][faceI]
}}}}
>>In line 173, neither of internalField[cellI] nor regionMesh[cellI] helps to resolve it.
*PS: if-loop is possible inside this structure? not quite sure about alpha[cellI]=0, which indicates the above mentioned condition for moisture content (embedded in coal) is zero when T = 100deg

>>rho_ = volScalarField
alphap = scalarField
Vxp = scalarField,

When adding rho_ in line 178, error pops up. Because rho_ is volScalarField, whereas alphap and Vxp are scalarField. Eliminating rho_ in line 178, has no error, however not realistic. I am checking as how to convert the volScalarField rho_ into scalarField function.
Quote:
ERROR: //reactingOneDim/reactingOneDim.C:173:11: error: ‘((Foam::regionModels:yrolysisModels::reactingOn eDim*)this)->Foam::regionModels:yrolysisModels::reactingOneD im::T’ does not have class type
Kindly check it in your free time please.

Thank you
Kummi is offline   Reply With Quote

Old   September 19, 2019, 09:24
Default Modification (2)
  #6
Senior Member
 
Kumaresh
Join Date: Oct 2016
Posts: 208
Rep Power: 6
Kummi is on a distinguished road
Send a message via Yahoo to Kummi
Dear Atul Kumar,
Thank you for your help. I have compiled successfully with further modifications.
Quote:
forAll(intCoupledPatchIDs_, i)
{
const label patchI = intCoupledPatchIDs_[i];
const scalarField& moisEvap = moisEva_.boundaryField()[patchI];
const scalarField& Tb = T_.boundaryField()[patchI]; // (1) - changed
const scalarField& rhob = rho_.boundaryField()[patchI]; // (2) - changed
scalarField& alphap = alpha_.boundaryField()[patchI]; //moisture fraction declaration
scalarField& Vxp = Vx_.boundaryField()[patchI]; //velocity declaration
forAll(moisEvap, faceI)
{
const labelList& cells = boundaryFaceCells_[faceI];
scalar evapRate = 0.0; //mass flux
scalar heatFlux = 0.0; //heat flux
forAll(cells, k)
{
const label cellI = cells[k];
if (Tb[cellI] == Tsat) // (3) - changed
{
alphap[cellI] = scalar(0.0);
}
evapRate +=rhob[cellI]*alphap[cellI]*Vxp[cellI]; // (4) - changed
heatFlux +=evapRate*Enthalp.value(); //W/m2
}
moisEva_.boundaryField()[patchI][faceI] += evapRate; //boundaryField()[patchI][faceI]
}
}
>>Corrections made: rho_ and T_ = volScalarField are converted into scalarField at boundary using declaration by boundaryField()[patchI].


However, have a long way to go


Thank you
Kummi is offline   Reply With Quote

Old   September 19, 2019, 10:49
Default
  #7
Senior Member
 
Kumaresh
Join Date: Oct 2016
Posts: 208
Rep Power: 6
Kummi is on a distinguished road
Send a message via Yahoo to Kummi
After compiling case file, I came up with the error while running it.
Quote:
kummiGasFoam: symbol lookup error: /home/kumaresh/OpenFOAM/OpenFOAM-2.1.1/cokeovenGasFOAM/lib/libpyrolysisModels.so: undefined symbol: _ZNK4Foam12regionModels15pyrolysisModels14reacting OneDim7moisEvaEv
Seems some linking problem.


Can anyone tell me what exactly this error is all about?
Kummi is offline   Reply With Quote

Old   September 21, 2019, 00:24
Default
  #8
Senior Member
 
Kumaresh
Join Date: Oct 2016
Posts: 208
Rep Power: 6
Kummi is on a distinguished road
Send a message via Yahoo to Kummi
Resolved above error by adding appropriate linking files
Quote:
virtual const surfaceScalarField& moisEva() const;
const surfaceScalarField& noPyrolysis::moisEva() const
{
FatalErrorIn("const volScalarField& noPyrolysis::moisEva() const")
<< "moisEva field not available for " << type() << abort(FatalError);
return surfaceScalarField::null();
}
noPyrolysis.C and noPyrolysis.H files.
Kummi is offline   Reply With Quote

Old   September 21, 2019, 01:21
Default Compiled successfully
  #9
Senior Member
 
Kumaresh
Join Date: Oct 2016
Posts: 208
Rep Power: 6
Kummi is on a distinguished road
Send a message via Yahoo to Kummi
Chosen Gauss linear scheme basically,
Quote:
div((interpolate(moisEva)*Cpgas),T) Gauss linear;
//Gauss linear, Gauss upwind, Gauss limitedLinear 1
Compiled successfully, but when the simulation gets started there comes error
Quote:
Evolving pyrolysis in region: oven
diagonal: Solving for rho, Initial residual = 0, Final residual = 0, No Iterations 0
#0 Foam::error:rintStack(Foam::Ostream&) at ??:?
#1 Foam::sigFpe::sigHandler(int) at ??:?
#2 in "/lib/x86_64-linux-gnu/libc.so.6"
#3 Foam::regionModels:yrolysisModels::reactingOneDi m::updateMoistureEvaporation() at ??:?
#4 Foam::regionModels:yrolysisModels::reactingOneDi m::evolveRegion() at ??:?
#5 Foam::regionModels:yrolysisModels:yrolysisMode lCollection::evolve() at ??:?
#6
at ??:?
#7 __libc_start_main in "/lib/x86_64-linux-gnu/libc.so.6"
#8
at ??:?
Guess, there is some problem in internal coding under updateMoistureEvaporation structure. Should look into it.

I guess the problem is due to alpha (moisture content). Because in my case, alpha has no expression (field operation), alpha is only solved based on the condition that, when T = 100deg, alpha =0 such that H2O(l) = H2O (g) @ T=100deg.Based on the condition, evaporation rate is calculated.

As seen in the code,
Quote:
if (Tb[cellI] == Tsat) { alphap[cellI] = scalar(0.0); }
evapRate +=rhob[cellI]*alphap[cellI]*Vxp[cellI];
In this equation, if the alphap[cellI] is equated with 0 (zero) , then the "evapRate" will also become zero right? If then, how to solve this condition?

If anyone have some ideas about it. Kindly do share. It will be highly helpful.
Thank you
Kummi is offline   Reply With Quote

Old   September 21, 2019, 04:00
Default
  #10
Senior Member
 
Kumaresh
Join Date: Oct 2016
Posts: 208
Rep Power: 6
Kummi is on a distinguished road
Send a message via Yahoo to Kummi
Source TERM 2 is created with the following condition:
Quote:
CONDITION T= 100 ==> alpha (moisture content) =0.
Based on it, mass and heat balances are calculated at the boiling plane (surface) as,
Mass balance (r) = rho * alpha * Velocity [kg/m2.s]
Heat balance (-K dT/dx) = (r) * latent heat [W/m2]
This condition is implemented in OpenFOAM by adopting two different Methods separately.
Quote:
METHOD 1
if (Tb[cellI] == Tsat) //Tsat = 100deg
{
alphap[cellI] = scalar(0.0);
}
evapRate +=rhob[cellI]*alphap[cellI]*Vxp[cellI];
Quote:
METHOD 2
if (alphap[cellI] >= scalar(0.1)) //10% of moisture
{
Tb[cellI] == Tsat; //Tsat = 100deg
}
if (alphap[cellI] < scalar(0.1)) //10% of moisture
{
alphap[cellI] = scalar(0.0);
}
evapRate +=rhob[cellI]*alphap[cellI]*Vxp[cellI];
Both structures are compiled successfully, however while running the case file - it leads to error as shown above. I am not sure how the temperature loop can be applied for such condition ?
Attachments 1 and 2 gives clear overview about the phenomenon and condition.

Can someone tell me ~ to make this condition work, which FOAM I should check ?

Thank you
Attached Images
File Type: jpg ATTACHMENT_1_Boiling.jpg (110.9 KB, 3 views)
File Type: png ATTACHMENT_2_Convection.png (153.7 KB, 2 views)
Kummi is offline   Reply With Quote

Reply

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
[swak4Foam] Installation Problem with OF 6 version Aurel OpenFOAM Community Contributions 9 November 30, 2018 11:20
Add different source term in diffusion equation at each time step Lewis Liang OpenFOAM Programming & Development 1 June 7, 2018 11:10
[swak4Foam] funkyDoCalc with OF2.3 massflow NiFl OpenFOAM Community Contributions 11 November 1, 2016 07:43
[Other] How to use finite area method in official OpenFOAM 2.2.0? Detian Liu OpenFOAM Meshing & Mesh Conversion 4 November 3, 2015 04:04
DxFoam reader update hjasak OpenFOAM Post-Processing 69 April 24, 2008 02:24


All times are GMT -4. The time now is 13:16.