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 10, 2019, 23:10
Question How to add Source term (2) for PYROLYSIS - reactingOneDim
  #1
Senior Member
 
Kumaresh
Join Date: Oct 2016
Posts: 138
Rep Power: 4
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, 14:57
Default
  #2
Member
 
Atul Kumar
Join Date: Dec 2015
Location: National Centre for Combustion Research and Development
Posts: 44
Rep Power: 5
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, 22:30
Default
  #3
Senior Member
 
Kumaresh
Join Date: Oct 2016
Posts: 138
Rep Power: 4
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, 02:57
Question SOURCE TERM 2 - Compilation ERROR
  #4
Senior Member
 
Kumaresh
Join Date: Oct 2016
Posts: 138
Rep Power: 4
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, 3 views)
File Type: png ATTACHMENT_2_Convection.png (153.7 KB, 3 views)

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

Old   September 17, 2019, 09:10
Default Modification (1)
  #5
Senior Member
 
Kumaresh
Join Date: Oct 2016
Posts: 138
Rep Power: 4
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   Yesterday, 08:24
Default Modification (2)
  #6
Senior Member
 
Kumaresh
Join Date: Oct 2016
Posts: 138
Rep Power: 4
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   Yesterday, 09:49
Default
  #7
Senior Member
 
Kumaresh
Join Date: Oct 2016
Posts: 138
Rep Power: 4
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

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


All times are GMT -4. The time now is 04:06.