View Issue Details

IDProjectCategoryView StatusLast Update
0002687OpenFOAMBugpublic2017-12-07 22:41
Reportermboesi Assigned Tohenry  
PrioritynormalSeveritymajorReproducibilityalways
Status closedResolutionsuspended 
PlatformGNU/LinuxOSUbuntuOS Version14.04
Product Versiondev 
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

Activities

henry

2017-09-05 18:11

manager   ~0008702

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?

henry

2017-09-05 19:35

manager   ~0008703

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.

mboesi

2017-09-06 08:58

reporter   ~0008704

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.

henry

2017-09-06 09:19

manager   ~0008705

> 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.

henry

2017-09-06 11:10

manager   ~0008706

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.

mboesi

2017-09-06 11:38

reporter   ~0008707

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_).

henry

2017-09-06 12:22

manager   ~0008708

> 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.

mboesi

2017-09-06 13:51

reporter   ~0008709

> 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.

henry

2017-09-06 14:02

manager   ~0008710

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?

mboesi

2017-09-06 14:38

reporter   ~0008711

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

henry

2017-09-06 15:03

manager   ~0008712

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.

henry

2017-09-06 15:11

manager   ~0008713

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.

mboesi

2017-09-06 15:34

reporter   ~0008714

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.

henry

2017-09-06 16:17

manager   ~0008715

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.

mboesi

2017-09-06 17:04

reporter   ~0008716

> 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.

henry

2017-09-06 17:58

manager   ~0008717

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?

mboesi

2017-09-06 19:39

reporter   ~0008718

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.

henry

2017-09-06 20:17

manager   ~0008719

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?

mboesi

2017-09-06 20:47

reporter   ~0008720

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.

henry

2017-09-06 20:53

manager   ~0008721

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.

mboesi

2017-09-06 21:53

reporter   ~0008722

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.

henry

2017-09-06 22:16

manager   ~0008723

> 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.

mboesi

2017-09-07 10:16

reporter   ~0008727

> 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.

mboesi

2017-09-07 17:55

reporter   ~0008728

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.

henry

2017-09-08 15:16

manager   ~0008731

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

mboesi

2017-09-11 13:55

reporter   ~0008748

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

mboesi

2017-09-12 20:58

reporter   ~0008769

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.

henry

2017-09-13 10:55

manager   ~0008770

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?

henry

2017-10-06 10:17

manager   ~0008820

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

mboesi

2017-10-10 10:07

reporter   ~0008848

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.

henry

2017-10-10 10:14

manager   ~0008849

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

henry

2017-12-01 11:06

manager   ~0009115

Shall I close this as not reproducible?

mboesi

2017-12-07 07:51

reporter   ~0009137

Hi, close it as not reproducible for now. I will be back when I have solved this issue.

henry

2017-12-07 22:41

manager   ~0009139

Pending update from reporter

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
2017-12-01 11:06 henry Note Added: 0009115
2017-12-07 07:51 mboesi Note Added: 0009137
2017-12-07 22:41 henry Assigned To => henry
2017-12-07 22:41 henry Status new => closed
2017-12-07 22:41 henry Resolution open => suspended
2017-12-07 22:41 henry Note Added: 0009139