[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

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
);

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)

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;
}

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), 553–561)

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]

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]

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]

Dear Atul Kumar,
I have created the structure for Source TERM 2 shown below, which is created under reactingOneDim.C for the following condition:

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]

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

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.

[Kummi]

Dear Atul Kumar,
Thank you for directing. I made some corrections.

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.

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]

Dear Atul Kumar,
Thank you for your help. I have compiled successfully with further modifications.

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]

After compiling case file, I came up with the error while running it.

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]

Resolved above error by adding appropriate linking files

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]

Chosen Gauss linear scheme basically,

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

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,

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]

Source TERM 2 is created with the following condition:

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.

METHOD 1
if (Tb[cellI] == Tsat) //Tsat = 100deg
{
alphap[cellI] = scalar(0.0);
}
evapRate +=rhob[cellI]*alphap[cellI]*Vxp[cellI];

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

[charles4allme]

Hi Kunmi,

Did you ever resolve the issue,

I am facing a similar issue with my own case.