2017-10-22 00:11 BST

View Issue Details Jump to Notes ]
IDProjectCategoryView StatusLast Update
0002687OpenFOAM[All Projects] Bugpublic2017-10-10 10:14
Reportermboesi 
Assigned To 
PrioritynormalSeveritymajorReproducibilityalways
StatusnewResolutionopen 
PlatformGNU/LinuxOSUbuntuOS Version14.04
Product Versiondev 
Target VersionFixed in Version 
Summary0002687: Error in pyrolysisChemsitryModel: Ys0 calculation will lead to negativ mass for consumed species
DescriptionIn the pyrolysisChemistryModel Ys0 represents the total mass of a solid species in the domain. The total mass is updated each time step when solving the pyrolysis chemistry:
if (updateC0)
{
    Ys0_[si][celli] -= sr*omegai;
}
Since the unit of omega and sr are kg/s and kg_product/kg_educt, this update procedure does ignore the time step size. Therefore, the calculated solid mass in the system will be wrong. The educt species mass eventually will go below zero. In case of produced species, the mass in the domain will be to high.

From my point of view, this issue could be easily solved by multiplying the RHS by the time step size, e.g.:
    Ys0_[si][celli] -= sr*omegai*dt[celli];

best regards,
Markus


Tagspyrolysis, solidChemsitryModel
Attached Files

-Relationships
+Relationships

-Notes

~0008702

henry (manager)

I agree that there are issues with the update of Ys0 but this is not just with the missing dt the bigger question is why is it updated explicitly at that point rather than being included in the ODE integration?

Also why is Ys0 update is both solve and calculate? Doesn't this mean it is updated at twice per time step? Could you check this?

I changed the code as you suggested:

    Ys0_[si][celli] -= sr*omegai*dt[celli];

but found that it does not compile. Could you provide a patch which includes all the changes needed to compile this change?

~0008703

henry (manager)

There is an additional problem with the current incrementing approach, even if it is corrected, in that it must be called exactly once per time-step. If the pyrolysis is included in a PIMPLE loop and called several time per timestep Ys0 would be incorrect and the error would accumulate.

~0008704

mboesi (reporter)

One can calculate the source term by using either the method calculate() or solve(). If you use solve() the source term is calculated in a coupled way.
So Ys0 is only updated once per chemistry update/correction.

I agree that there is another issue when using PIMPLE. This could be solved by just correcting the chemistry at the first or last PIMPLE loop.

The reason why it does not compile is that there is no variable dt defined. I will provide a patch for the pyrolysis class hopefully until tomorrow.

Concerning the PIMPLE loop issue I think that solving it only once is sufficient, since the pyrolysis model is used for solid regions only. But it would be good the get some other opinions on this.

~0008705

henry (manager)

> One can calculate the source term by using either the method calculate() or solve().

Calling both is not disallowed and depending on the particular way in which the chemistry is integrated and coupled with the other sources it might be useful. Having the integration of Ys0 implemented in such an ad hoc, inconsistent and restrictive manner is very poor programming. Ideally pyrolysisChemistryModel should be rewritten to avoid this nonsense and bring it up to the quality of the rest of OpenFOAM.

> Concerning the PIMPLE loop issue I think that solving it only once is sufficient

Surely that is case dependent. If the restriction of 1 call per time-step is maintained in pyrolysisChemistryModel then this should be checked for and a FatalError generated if it called more than once as this would introduce a cumulative error in the solution.

~0008706

henry (manager)

I think the best way forward would be to include Ys0 in the "c" array as an additional specie so that it is integrated consistently with the rest of the composition by the ODE solver.

~0008707

mboesi (reporter)

I thought about the implementation and think a good aproach would be the calculation according to the phase fraction and phase density. This would also allow the pyrolysis model to be used for Eulerian phases. The current implementation does not support any changes in the Ys0 due to transport phenomena. The new Ys0 should be calculated after calculating the source terms.

Including the Ys0 in the "c" array would double its size, since the Ys0 is per specie.

> Calling both is not disallowed and depending on the particular way in which the chemistry is integrated and coupled with the other sources it might be useful.

I am a heavy user of the chemistry classes and I think both methods should be retained. Both methods are also included in the gas phase chemistry, for a system of slow reactions one can use calculate() to be more efficient. Fast reactions have to be solved by sovle() since they require some sub time stepping. The method for source term calculation is specified in the combustion model (Switch integrateReactionRate_).

~0008708

henry (manager)

> Including the Ys0 in the "c" array would double its size, since the Ys0 is per specie.

Correct. Currently Ys0 is integrated using a 1 step Euler explict approach, why? Why not integrate Ys0 consistenly with the rest of the species? Why should it be treated fundamentally differently?

> I think both methods should be retained.

I wasn't suggesting removing either but currently if both are called for pyrolysisChemistryModel or either called more than once the results will be wrong. If the current approach for the integration of Ys0 is maintained then calling both of the functions or either more than once must be disallowed somehow.

~0008709

mboesi (reporter)

> Why not integrate Ys0 consistently with the rest of the species? Why should it be treated fundamentally differently?

The pyrolysisChemistryModel class calculates the chemical source terms (RR(Yi)) for the species conservation equations. The species conservation equation is solved for Y. Ys0 = Y * rho, thus I think it would be redundant to include Ys0 in the "c" array, since it could be easily reconstructed by Ys0 = Y * rho. Assuming a solid region, e.g. wood, the mass conservation would be ddt(rho,Yi) = RRs(Yi). So including the Ys0 to "c" would mean you solve the same thing twice. Once in when the species conservation is solved and once when pyrolysis source term is corrected.

> If the current approach for the integration of Ys0 is maintained then calling both of the functions or either more than once must be disallowed somehow.

I totally agree with that. If the Ys0 would be calculated by Ys0 = Y * rho, it should work with the current implementation. An inconsistency would exist, since the Y from the previous time step would be used to update the Ys0. In order to fix this, The update of Ys0 has to be done after solving the species conservation.

~0008710

henry (manager)

As Ys0 = Y*rho which can be rewritten in terms of c why is it needed at all? surely the reaction rate expressions can be written in terms of c only.

Alternatively if Ys0 is treated as Y*rho this can be calculated within the reaction rate function as required so there would be no need to store and update Ys0 or am I missing something?

~0008711

mboesi (reporter)

The pyrolysis calculates the reaction rates based on total mass (Ys0). I looked it up, Ys0 is used to calculate the mass dependent part of the reaction rate: kf *= pow(c[i]/Ys0[i], exp)*Ys0.
Since c[i] = rho * Y[i] * mesh.V(), Ys0 should be equal to c. From my point of view also the mass dependent part seems to be in valid, because if c = Ys0 the above equation would simplify to kf *= pow(1, exp)*Ys0. Therefore, I think Ys0 could be omitted for the calculations: kf *= pow(c1, exp) should be good. Additionally, there would be no need to store it.

For solid reactions it would be nice to have an information about the absolute mass in the system, this might be one reason to keep it.

PS: I had the calculation of Ys0 wrong, the correct formula is: Ys0 = Y * rho * V

~0008712

henry (manager)

I agree with your assessement, Ys0 should be factored out of the expressions, there should be no need for it.

Also the local renaming of cell volumes to delta makes no sense, why not call it V?

There is also a lot of inefficiency in the implementation, particularly with the local allocation of the c array which should be cached as work space. Another inefficiency is

    scalarField delta(this->mesh().V());

which should be

    const scalarField& delta(this->mesh().V());

because there is no need to copy the cell volumes.

~0008713

henry (manager)

Actually it would be better if delta (V) were also factored out of the implementation altogether; it is not needed and adds complexity, confusion and inefficiency.

~0008714

mboesi (reporter)

I thought about that too, but typical reaction rates for solids are mass based (kg/kg^exp/s). I am not sure if the results are the same when using partial densities (omitting V), since the reaction order might be different from 1.

I mean it makes totally sense to do this calculations mass based, due to the nature of solid reaction rates.

~0008715

henry (manager)

Given that V does not change the results would be the same with or without it included in the formulation. The resulting reaction rates are divided by V at the end of the solution so why have V included before then? It looks like clutter to me.

Why should solid reaction rated by mass based but not liquid or gas reaction rates? I don't see an advantage in this inconsistency.

~0008716

mboesi (reporter)

> Why should solid reaction rated by mass based but not liquid or gas reaction rates?

The point is that concentrations can be more or less easily measured in gas and liquids. For determining reaction rates in gases or liquids you take concentrations. But typically solid reaction rates are determined by mass loss profiles e.g. in thermal gravimetric analysis (TGA). Therefore, effective rate expressions have the unit 1/s/kg_solid. However, the pre-exponential factors A have the unit 1/s/kg_solid^(order-1), since the total kinetic rate is calculated by k = A*exp(-Ta/T)*c^order. c is the mass of the reacting solid in kg.

If we take partial densities things would be a bit different: kdensity = A*exp(-Ta/T)*c^order/V^order. Dividing kdensity by k/V gives V^order/V which is only the same if the order is equal to one. From my experience the order is usually between -1 up to 2.5 or so.

This gives three possibilities:
-) Calculate every thing based on mass for solid chemistry
-) Use partial densities and correct the source term by V^order. This has to be done for every individual reaction and cell.
-) Introduce a scaling of the pre-exponential factor for A. This scaling factor depends on V and order, meaning you need a different A for every cell.

~0008717

henry (manager)

Does this mean that the reaction rate coefficients are directly dependent on the cell volumes? If so how are they specified but if not why is it not possible to remove the cell volumes from the formulation?

~0008718

mboesi (reporter)

From my point of view the volume is necessary as a consequence of the measurement methods and the provided kinetic parameters for solid reactions.

I am not 100% sure but I think that it is not possible to remove the volume from the calculations (except the ones mentioned above). I will check that and come back to you.

~0008719

henry (manager)

Does this mean that the kinetic parameters are directly dependent on the size of the mesh cells and if the mesh is refined or the cell volumes changed in any way the kinetic parameters have to be changed correspondingly?

~0008720

mboesi (reporter)

No the kinetic parameters are not cell size dependent. But the reaction rate is indirectly depending on the cell size, since cell size, density and mass fraction define the mass of a specie in the cell.

> -) Introduce a scaling of the pre-exponential factor for A. This scaling factor depends on V and order, meaning you need a different A for every cell.

If you were referring to this, the answer is yes, pre-exponential factors directly depend on cell size and have to be changed correspondingly.

~0008721

henry (manager)

If you have a large mesh with lots of different cell sizes how do you specify all the parameters for each cell? If you apply automatic refinement how would the parameters be updated to compensate for the change in cell volumes?

I do not think it is practical to have the solid reaction parameters to be directly dependent on the cell volumes; if this is absolutely necessary I am not sure that CFD of solid reactions is realistic.

~0008722

mboesi (reporter)

If the absolute reaction rate (kabs = k(T)*k(c)) is calculated with c = Y*rho*V, there is nothing to update or change. It should work without any problem.

I think I understand your doubts about solid reactions in CFD now. So let's assume a solid having mass m, density rho, occupying volume V is reacting. V is divided into N volume cells (having volume Vc). If we assume a constant temperature, k(T) is the same in each cell, however k(c) depends on Vc. The absolute reaction rate in a cell would be kabsCell = k(T)*(Y*rho*Vc)^order. To get the cumulative absolute reaction rate all kabsCell have to be summed up giving kabs = k(T)*(Y*rho*V)^order.

If this solid consists of one huge cell only, the absolute reaction rate would also be kabs = k(T)*(Y*rho*V)^order.

Therefore, cell refinement or similar things do not effect the cumulative rate. Yes, the rates per cell are affected, but not the cumulative rate. The fundamental thing of solid reactions is to determine the mass loss of the whole solid, which is cell size independent.

~0008723

henry (manager)

> If the absolute reaction rate (kabs = k(T)*k(c)) is calculated with c = Y*rho*V, there is nothing to update or change. It should work without any problem.

It would run but the results would depend directly on the mesh so if you wanted mesh independent results you would need to change the parameters every time the mesh is changed. If the mesh is refined and unrefined automatically then I don't see how this would be possible.

From what you are saying only the total reaction rate of the solid is valid and the distribution of reaction rate has no physical meaning as it is directly dependent on the mesh, in which case the simulation of the solid should only ever be performed on one large cell.

~0008727

mboesi (reporter)

> It would run but the results would depend directly on the mesh so if you wanted mesh independent results you would need to change the parameters every time the mesh is changed.

I think the is no way to get mesh independent results for the absolute reaction rate of solids, since they are usually based on mass loss measurements.

I would not say that the reaction rate has no physical meaning. When modelling a solid by volume cells, the perstective is changed from mass based to volume based. Therefore, purely mass based kinetic data is not volume mesh independent.


I found one publication employing a two-fluid model to model the combustion of wood chips. Equations 33 - 43 explain the calculation of the solid kinetic rates in their model.
http://www.tandfonline.com/doi/full/10.1080/13647830.2011.610903?src=recsys


From my point of view, the implementation using the total mass (rho*Y*V) is perfectly fine for solid reactions.

~0008728

mboesi (reporter)

I have attached a patch fixing the issues with Ys0_ and the reaction order calculation kf *= pow(c,exp) now.
The updateC argument was also removed from the omega(c,T,p) method.

~0008731

henry (manager)

I see that Ys0_ is still calculated but not used for anything, can it be removed altogether? It would save a lot of storage.

~0008748

mboesi (reporter)

It could be removed. I am going to rework this.

~0008769

mboesi (reporter)

I have removed the Ys0 calculation and fields.

Additionally, I have changed the rate calculation, since I have messed up the units of gas-solid reactions and pyrolysis reactions. The units of pyrolysis reactions is 1/s or 1/s/kg^order. After reading some literature on pyrolysis kinetics and modelling, I think the best way to calculate the reaction rate would be k = kf * c^order * c[lhs[0].si]. Pyrolysis reactions should have only one lhs solid species.

~0008770

henry (manager)

Can you provide a validation case which demonstrates that the updated implementation is correct and that the results are better than with the current implementation?

~0008820

henry (manager)

Have you managed to create a case which demonstrates the updated implementation is correct? If not should I close this report as "not reproducible"?

~0008848

mboesi (reporter)

Hi, I am still working on it. I am quite busy at the moment but I try to provide the test case and the implementation by end of Oktober.

~0008849

henry (manager)

That's fine, I will leave this report open until you are ready.
+Notes

-Issue History
Date Modified Username Field Change
2017-09-05 15:43 mboesi New Issue
2017-09-05 15:43 mboesi Tag Attached: pyrolysis
2017-09-05 15:43 mboesi Tag Attached: solidChemsitryModel
2017-09-05 18:11 henry Note Added: 0008702
2017-09-05 19:35 henry Note Added: 0008703
2017-09-06 08:58 mboesi Note Added: 0008704
2017-09-06 09:19 henry Note Added: 0008705
2017-09-06 11:10 henry Note Added: 0008706
2017-09-06 11:38 mboesi Note Added: 0008707
2017-09-06 12:22 henry Note Added: 0008708
2017-09-06 13:51 mboesi Note Added: 0008709
2017-09-06 14:02 henry Note Added: 0008710
2017-09-06 14:38 mboesi Note Added: 0008711
2017-09-06 15:03 henry Note Added: 0008712
2017-09-06 15:11 henry Note Added: 0008713
2017-09-06 15:34 mboesi Note Added: 0008714
2017-09-06 16:17 henry Note Added: 0008715
2017-09-06 17:04 mboesi Note Added: 0008716
2017-09-06 17:58 henry Note Added: 0008717
2017-09-06 19:39 mboesi Note Added: 0008718
2017-09-06 20:17 henry Note Added: 0008719
2017-09-06 20:47 mboesi Note Added: 0008720
2017-09-06 20:53 henry Note Added: 0008721
2017-09-06 21:53 mboesi Note Added: 0008722
2017-09-06 22:16 henry Note Added: 0008723
2017-09-07 10:16 mboesi Note Added: 0008727
2017-09-07 17:55 mboesi File Added: solidChemistryModel.tar.gz
2017-09-07 17:55 mboesi Note Added: 0008728
2017-09-08 15:16 henry Note Added: 0008731
2017-09-11 13:55 mboesi Note Added: 0008748
2017-09-12 20:58 mboesi File Added: solidChemistryModel_Rate_fix.tar.gz
2017-09-12 20:58 mboesi Note Added: 0008769
2017-09-13 10:55 henry Note Added: 0008770
2017-10-06 10:17 henry Note Added: 0008820
2017-10-10 10:07 mboesi Note Added: 0008848
2017-10-10 10:14 henry Note Added: 0008849
+Issue History