[Kummi]

My code is contructed on a region for MOISTURE EVAPORAION as shown below:

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

forAll(intCoupledPatchIDs_, i)
{
const label patchI = intCoupledPatchIDs_[i];
const scalarField& moisEvap = moisEva_.boundaryField()[patchI];
const scalarField& Tb = T_.boundaryField()[patchI];
const scalarField& rhob = rho_.boundaryField()[patchI];
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;
scalar heatFlux = 0.0;
//************************************************** ******************//
forAll(cells, k)
{
const label cellI = cells[k];
if (Tb[cellI] == Tsat) // defining at T= 100deg
{
alphap[cellI] = scalar(0.0);
}
evapRate +=rhob[cellI]*alphap[cellI]*Vxp[cellI]; //mass flux
heatFlux +=evapRate*Enthalp.value(); //W/m2 - //heat flux
}
moisEva_.boundaryField()[patchI][faceI] += evapRate;
//************************************************** ******************//
}
}
}

This above code is written in simple manner and compiled without errors, but I face errors while running the case file. Inside the commented lines, the main code structure for moisture evaporation is constructed, which seems to be not satisfied the flow physics.
Kindly look into the attached image please. Here the Equs. 6-8 define the position of plane. With the position of plane, Equs. 9 and 10 are calculated.

(1) Equ. 9 makes me think that r = rho*w*(dxb/dt) , where dxb/dt is the movement of plane, which I have not introduced in above code inside commented lines

(2) Moreover the moving plane (xb) has two sides, high (-) and low (-) temp sides, and depending upon it calculations are made, which I have not introduced in above code inside commented lines

I feel the above equations are simple ~ but I can't able to figure out the OpenFOAM coding technique as how to approach and solve.

Correct me if im wrong anywhere please.

Kindly share your ideas and ask me anything if you are not clear

Thank you

[gschaider]

Originally Posted by Kummi

My code is contructed on a region for MOISTURE EVAPORAION as shown below: This above code is written in simple manner and compiled without errors, but I face errors while running the case file. Inside the commented lines, the main code structure for moisture evaporation is constructed, which seems to be not satisfied the flow physics.
Kindly look into the attached image please. Here the Equs. 6-8 define the position of plane. With the position of plane, Equs. 9 and 10 are calculated.

(1) Equ. 9 makes me think that r = rho*w*(dxb/dt) , where dxb/dt is the movement of plane, which I have not introduced in above code inside commented lines

(2) Moreover the moving plane (xb) has two sides, high (-) and low (-) temp sides, and depending upon it calculations are made, which I have not introduced in above code inside commented lines

I feel the above equations are simple ~ but I can't able to figure out the OpenFOAM coding technique as how to approach and solve.

Correct me if im wrong anywhere please.

Kindly share your ideas and ask me anything if you are not clear

Thank you

I'm sorry. I don't do that kind of support here


Just one remark: at the first glance I've seen that you're checking for the Tsat with the == operator. This is almost never a good idea in numerical codes as it is unlikely that the temperature is EXACTLY that value when being something solved for. Use a check with a tolerance.

[Kummi]

Thank you for your response Bernhard.

Tsat == operator - This is almost never a good idea in numerical codes as it is unlikely that the temperature is EXACTLY that value when being something solved for. Use a check with a tolerance.

I hope you are mentioning if loop statement if (Tb[cellI] == Tsat) mentioned above in my code.

If you don't mind, can you please elaborate this stament "Use a check with a tolerance"
Thank you once again ^^

[gschaider]

Originally Posted by Kummi

Thank you for your response Bernhard.
I hope you are mentioning if loop statement if (Tb[cellI] == Tsat) mentioned above in my code.

If you don't mind, can you please elaborate this stament "Use a check with a tolerance"
Thank you once again ^^

Means that you should check for the absolute value of the difference of T in the cell and Tsat to be smaller than a threshold


|T-Tsat| < eps

[Kummi]

Thank you Bernhard. I have tried implementing your idea by checking with a tolerance.

The structure below is a part of my main code. It compiles well and when running case file, it leads to error at 1968.6s

scalar Tsat = 373;
scalarField& Tb = T_.boundaryField()[patchI];
scalarField& rhob = rho_.boundaryField()[patchI]; // density
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
forAll(cells, k)
{
const label cellI = cells[k];
// if (Tb[cellI] - Tsat <= scalar(1.0))
if (Tb[cellI] - Tsat == scalar(1.0)) // if (Tb[cellI] == Tsat)
{
alphap[cellI] = scalar(0.0);
}
else
{
alphap[cellI] = alphap[cellI];
}
evapRate +=rhob[cellI]*alphap[cellI]*Vxp[cellI];
moisEva_.boundaryField()[patchI][faceI] += evapRate;

Time = 1968.6
Evolving pyrolysis in region: furnace
diagonal: Solving for rho, Initial residual = 0, Final residual = 0, No Iterations 0
diagonal: Solving for Ycoal, Initial residual = 0, Final residual = 0, No Iterations 0
DICPCG: Solving for T, Initial residual = 2.17328e-05, Final residual = 3.9941e-21, No Iterations 1
pyrolysis min/max(T) = 402.532, 1373

Evolving pyrolysis in region: oven
diagonal: Solving for rho, Initial residual = 0, Final residual = 0, No Iterations 0
diagonal: Solving for Ycoal, 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::PCG::solve(Foam::Field<double>&, Foam::Field<double> const&, unsigned char) const at ??:?
#4 Foam::fvMatrix<double>::solve(Foam::dictionary const&) at ??:?
#5 Foam::fvMatrix<double>::solve() at ??:?
#6 Foam::regionModels:yrolysisModels::reactingOneDi m::solveEnergy() at ??:?
#7 Foam::regionModels:yrolysisModels::reactingOneDi m::evolveRegion() at ??:?
#8 Foam::regionModels:yrolysisModels:yrolysisMode lCollection::evolve() at ??:?
#9
at ??:?
#10 __libc_start_main in "/lib/x86_64-linux-gnu/libc.so.6"
#11
at ??:?
Floating point exception (core dumped)

In this code,
(1) if (Tb[cellI] == Tsat) --> ERROR

(2) I have tried if (Tb[cellI] - Tsat == scalar(1.0)), the same error continues


My instinct says, there is a error in if-else loop - not sure though.

Kindly share your ideas please.

[Kummi]

And followed by above code, my energy source code structure is:

void reactingOneDim::solveEnergy()
{
if (debug)
{
Info<< "reactingOneDim::solveEnergy()" << endl;
}
const volScalarField rhoCp(rho_*solidThermo_.Cp());
const surfaceScalarField phiQr(fvc::interpolate(Qr_)*nMagSf());
const surfaceScalarField phiGas(fvc::interpolate(phiHsGas_));
scalar Tsat1 = 373; //AT Tsat=100deg, moisture evaporates
scalar Tres = 773; //AT Tres (Resolidifcation temp) =500deg, coal decomposition
dimensionedScalar Cpgas("Cpgas", dimEnergy/dimMass/dimTemperature, 2000);
fvScalarMatrix TEqn
(
fvm::ddt(rhoCp, T_) // + fvm::div(phi, T)
- fvm::laplacian(K_, T_)
==
chemistrySh_
+ fvc::div(phiQr)
+ fvc::div(phiGas)
);

////////////////////////////// NEW //////////////////////////////
if (moistureHSource_)
{
const surfaceScalarField phimoisEva(fvc::interpolate(moisEva_));
forAll(T_, cellI) //forAll(cells, k) ~ forAll(mesh.cells(), celli)
{
if ((T_[cellI] > Tsat1) && (T_[cellI] < Tres)) //(100deg < T_[cellI] < 500deg)
{
TEqn += fvc::div(phimoisEva*Cpgas, T_);
}
else if ((T_[cellI] <= Tsat1) || (T_[cellI] >= Tres)) //(100deg >= T_[cellI] >= 500deg)
{
TEqn += scalar(0.0);
}
}
}
////////////////////////////// NEW //////////////////////////////
TEqn.relax();
TEqn.solve();
Info<< "pyrolysis min/max(T) = " << min(T_).value() << ", "
<< max(T_).value() << endl;
}

The commented "NEW" in above structure is the newly introduced code.

Did my above error is due to the false construction of NEW structure? If any ideas, kindly do share please.
Thank you for your time ^^