firemodels / fds

Fire Dynamics Simulator
https://pages.nist.gov/fds-smv/
Other
660 stars 622 forks source link

char oxidation consistency in vegetation #11326

Closed ericvmueller closed 1 year ago

ericvmueller commented 1 year ago

This is a modified version of the char_oxidation verification cases. A single particle of CHAR is exposed to a heat flux.

In the 'lumped' case, the char oxidation follows the shortcut of consuming AIR and producing PRODUCTS (and ASH). In the 'tracked' case, CHAR is converted to CARBON DIOXIDE and ASH, while consuming OXYGEN, with both gas species tracked directly.

The mass consumption and generation all seem consistent, but the temperature of the particle during the peak of the reaction is very different. I've tried using the 'tracked' approach but mimicking the 'lumped' approach (so some WATER VAPOR is produced instead of all CARBON DIOXIDE), but the result is still very different to the 'lumped' approach.

image

char_oxidation_lumped.fds.txt char_oxidation_tracked.fds.txt

rmcdermo commented 1 year ago

Are the specific heats supposed to be different for the CHAR MATL?

ericvmueller commented 1 year ago

No sorry, that was something I was messing around with. The uploaded files should be fixed now, but the issue is the same

ericvmueller commented 1 year ago

My guess right now is it comes down to the enthalpy calculation here. In the lumped approach the values are larger because of the larger NU_GAS, and the particle may be more sensitive to Q_DOT_S_PPP since the energy is applied directly in the 1D heat transfer while Q_DOT_G_PPP is added to the divergence. So even though energy is conserved the response of the particle is different.

https://github.com/firemodels/fds/blob/7ca9172e26a98618483e1bbc1be878191cdda2f1/Source/wall.f90#L2971-L2975

drjfloyd commented 1 year ago

Eric, I think you are correct. The lumped case is not indentical to the tracked case. In the lumped case you have PRODUCTS and AIR and in the tracked case you have CARBON DIOXIDE and OXYGEN. PRODUCTS is not the same as CARBON DIOXIDE and AIR is not the same as OXYGEN.

mcgratta commented 1 year ago

Right, but the question is --- should we expect such a large temperature difference?

drjfloyd commented 1 year ago

I'll have to run and look at the cases.

mcgratta commented 1 year ago

The lines you've cited were added in the past few months. I'm trying to understand why, in this particular case, the energy associated with the consumption of AIR and production of PRODUCTS should be assigned to the solid and gas heating terms, respectively.

drjfloyd commented 1 year ago

We produce gas at the solid temperature in the solid phase but in the gas phase we add the gas at the current temperature. So if the solid and gas temperature are not the same we need to account for that as an energy source term. We do the same thing for droplet evaporation (the DELTA_H_G line 3182 in part).

The same is true for gas consumption. We pull the gas into the solid at the gas temperature but once it become solid it is at the solid temperature. We again need an energy source term, this time in the solid, to account for this.

mcgratta commented 1 year ago

Let's try this

Q_DOT_G_PPP = Q_DOT_G_PPP + ML%ADJUST_BURN_RATE(NS,J)*ML%NU_GAS(NS,J)*RHO_DOT*(H_S-H_S_B)

instead of the IF statement that you cite above. That is, let's put all the changes in sensible enthalpy into the gas.

ericvmueller commented 1 year ago

Okay, I will compare the two cases again with that change

drjfloyd commented 1 year ago

It is not correct to put it all in the gas. Say your solid phase reaction is Fe(s) + O2(g) -> FeO2(s).

The specific heat of CO2 is slightly less than the specific heat of AIR. The gas temperature is about 1000 K away from the solid temperature at the peak of the reaction (see the attached plots shown temperature difference and reaction rate). With 1000 K delta that small difference is going to add up.

image

ericvmueller commented 1 year ago

The differences are certainly reduced with that change, but I also appreciate what Jason is saying. Maybe we need to rethink this lumped shortcut for vegetation cases.

image

drjfloyd commented 1 year ago

With the lumped shortcut you are making CO2 and H2O as products and not just CO2.

We now have the ability to specifiy multiple simple chemistry reactions. If you want ultimately a case where you have both a reaction for virgin foliage and a reaction for char oxidation, you could define two simple chemistry reactions, one with the CHO for the virigin pyrolyzate and the second just being C. Then for the char oxidation you would specify the products for that second reaction rather than the products for the first.

drjfloyd commented 1 year ago

The underlying issue is we have inputs that say all of air is consumed to make all of products, but in reality, only part of air (the O2) is consumed to make part of products (the CO2). The rest is in reality staying in the gas phase; however, the wall model doesn't know that. I changed the REAC to C->CO2. Then went into wall and just hard coded that 23 % of the air (the O2) applies in the Q_DOT_S_PPP term and that 29 % of the products (the CO2 from the char) applies and reran. When I do that the tracked and lumped become very close.

image

I'll give some thought about how we might deal with that either fully within FDS or possibly with a little more input needed from users (i.e., a second list of species for MATL giving the component of a lumped species that is really participating).

rmcdermo commented 1 year ago

Nice. Seems we should be able to make this work without explicitly tracking CO2, yes?

Similar to how evaluate the enthalpy of the stoichiometric pocket for the extinction model, here we can just grab the reactant and product gases from the lumped species for Q_DOT_S_PPP.

drjfloyd commented 1 year ago

Yeah that would be the idea.

drjfloyd commented 1 year ago

We might be able to do this automatically. We could sum NU(SM)*SM%MASS_FRACTION for all positve NUs and then again for all negative NUs and look to see what is consumed and what is made.

The simple, but not rigorous, approach would be to do my hack comcept where we use the computed fractions of mass that particpate and the lumped enthalpies. This would have some error in global conservation in that the primtive species specific heat is likely to be different than the lumped species specific heat.

The more complex approach would be to use the primitive species enthalpies rather than the lumped enthalpies. This would require storing primitive species enthalpies (not a huge memory requirement and we already do this for liquid specific heat for use in evaporation), creating a set of NU arrays for primitive species, and in wall loop over primitive rather than tracked species when computing the Q_DOT terms. The M_DOT can stay as tracked.

rmcdermo commented 1 year ago

Couple of questions:

  1. What is the sensitivity of the particle temperature to the uncertainty in HEAT_OF_REACTION?
  2. How would an example MATL look like that used the CHAR + O2 => CO2 reaction with lumped species only in gas phase? Would it be something like
    
    &SPEC ID='OXYGEN', LUMPED_COMPONENT_ONLY=T/
    &SPEC ID='CARBON DIOXIDE', LUMPED_COMPONENT_ONLY=T/

&MATL ID='CHAR' ... NU_SPEC = 2.61,-1.65 SPEC_ID = 'CARBON DIOXIDE','OXYGEN' HEAT_OF_REACTION = -12000. /


The only restriction I can see is that 'CARBON DIOXIDE' and 'OXYGEN' should each only be part of one lumped species.
mcgratta commented 1 year ago

Suppose we make a slight modification of the syntax

NU_SPEC   = 2.65,-1.65
SPEC_ID   = 'PRODUCTS%CARBON DIOXIDE','AIR%OXYGEN'

Could you work with this information? I think the parsing of the character strings would be fairly easy.

drjfloyd commented 1 year ago

@rmcdermo 1. is a good question. I can run a couple variants of the case and see.

For 2. I think we can do this without more user inputs. We know how much gas wind up in the solid (if the net gas nu is less than the net solid nu we know some gas is in the solid). We know from the lumped species what primitive species is consumed and what primitive species is made. If this doesn't work Kevin's idea is a pretty good one.

rmcdermo commented 1 year ago

Might be easier to have a corresponding SMIX_ID list, for example...

NU_SPEC   = 2.65,-1.65
SPEC_ID   = 'CARBON DIOXIDE','OXYGEN'
SMIX_ID   = 'PRODUCTS','AIR'
rmcdermo commented 1 year ago

You don't need more user inputs unless one of the SPEC in the reaction is part of more than one LUMPED SPEC. Then you have to figure out which one to push out of the solid.

ericvmueller commented 1 year ago

To Randy's first point... I made a quick hack as Jason described to scale the enthalpy change by the sub-component mass fraction and tested with HR=12 MJ/kg and HR=18 MJ/kg. I also tested a case with the enthalpy change entirely commented out (Kevin was curious about this).

For this particular test case, the effect of HR uncertainty (which I dont think is unreasonable for the wide range we see suggested in wildfire literature) is on the same order as the effect of uncertainty enthalpy.

image

drjfloyd commented 1 year ago

With the current inputs the user is giving tracked (lumped) species. With NU and SM%MASS_FRACTION. I know the mass of each primitive species being consumed and the mass of each primitive species being produced. I can sum the two over all species and any primitive species that are net consumed go to the solid and any that are net produced go to the gas. I don't think we will need new inputs.

mcgratta commented 1 year ago

Keep in mind that if this is done in the PYROLYSIS subroutine, it could add up.

drjfloyd commented 1 year ago

This can be done in read and we make a second NU array for primitive species. The tracked species array gets use for M_DOT_G and the primitive for Q_DOT_G and Q_DOT_S

rmcdermo commented 1 year ago

OK, this makes an assumption about how the primitive species will be apportioned. Suppose the primitive product is 'WATER VAPOR'. You will have a 'PRODUCTS' with 'WATER VAPOR' and you may also be explicitly tracking 'WATER VAPOR' as its own lumped species. In that case, you'd probably just want to push 'WATER VAPOR' out of the solid, not 'PRODUCTS'. You would not want to push some 'PRODUCTS' and some 'WATER VAPOR'.

drjfloyd commented 1 year ago

The only primitive species that get considered are those listed in the SPEC_ID on the MATL line.

rmcdermo commented 1 year ago

I understand. And each of those could be part of multiple lumped species.

The NU and the SPEC on MATL define the reaction, correct?

drjfloyd commented 1 year ago

&REAC FUEL='WOOD GAS',C=3,H=5,O=2/ &REAC FUEL='CHAR',C=1/ &MATL ID='WOOD',A=1,E=1,HEAT_OF_REACTION=1,SPEC_ID='WOOD GAS',NU_SPEC=0.8,MATL_ID='CHAR',NU_MATL=0.2/ &MATL ID='CHAR', A=1,E=1,HEAT_OF_REACTION=1,SPEC_ID='AIR','PRODUCTS 2',NU_SPEC=-11.5334,12.5334/

The first MATL is easy, we are just making species. For the second, do NU_SPEC(AIR)SM(AIR)%MASS_FRACTION(:) and get the mass of all the species consumed. Do NU_SPEC(PRODUCTS 2)SM(PRODUCTS 2)%MASS_FRACTION(:) and get the mass of all the species consumed. Add these two together and you get that for 1 kg of CHAR you consume 2.664 kg of O2 and make 3.664 kg of CO2. Store these values in a tracked species NU array and use that NU array plus the primitive species enthalpy in wall. We'd need to store the primitive species enthalpy like we do for liquid enthalpy. Not a big deal as we already compute the primitive enthalpy in order to compute the lumped enthalpy. Just need to store it.

rmcdermo commented 1 year ago

OK, but you are assuming that O2 is the reactant in AIR that reacts with the CHAR. I don't see this as general.

drjfloyd commented 1 year ago

I am not assuming that. That is what the user has defined via the inputs.

rmcdermo commented 1 year ago

Where?

drjfloyd commented 1 year ago

The User has specified that the wall takes in 11 kg of AIR per kg of CHAR and that it spits out 12 kg of PRODUCTS 2 per kg of CHAR. AIR and PRODUCTS 2 have a speficic set of primitive species that they are composed of. The only species that gets consumed is OXYGEN.

rmcdermo commented 1 year ago

OK, sorry for being dense. Seems like we probably should have been doing this all along. Let's first write things up in the tech guide and discuss next week.

ericvmueller commented 1 year ago

One question: does this method imply that two additional lumped species be tracked for the second REAC? So we lose some computational advantage of the original lumped approach for char?

drjfloyd commented 1 year ago

Yes it means two additional species. I ran the case with 1 and 2 REAC lines each three times. The 2 REAC case took 15 % more time due to spending about 50 % more time in DIVG and MASS, about 30 % more time in FIRE, and no significant change in the other major areas of the code.

With one REAC, your char oxidiation isn't making just CO2, it is making CO2 and H2O (the products made when your REAC fuel has C and H) and possibly CO and SOOT if you have yields for those.

ericvmueller commented 1 year ago

Understood - I think the justification for this (19.1.1 of the users guide) is that the discrepancy is minor enough to justify the benefits in terms of simplicity and speed. Particularly for larger WUI cases. The discussion about enthalpy here suggests this is currently not the case. However, if we only track net enthalpy but still create all PRODUCTS from the original REAC (essentially ignore nitrogen in the enthalpy change), can we still get away with this simplification?

I think the change you are suggesting will allow us to test that though, so we are free to explore the relative importance in different cases.

drjfloyd commented 1 year ago

Made the changes for tracking products; still need to run them through firebot. I found a couple other things that needed addressing (dealing with NU_PART and cases where the yields do not sum to 1) so wouldn't surprise me if something broke.

I ran the lumped and tracked cases from the opening post to this issue where the lumped REAC was set to C=1 (changed the lumped input to set HUMIDTY and Y_X_INFTY to be the same as the tracked). The surface temperature plot is now showing nearly identical results.

image

drjfloyd commented 1 year ago

I also ran the lumped case with the original REAC (labeled 1 REAC in the plot). It is more different from the tracked than the case where the REAC is C=1, but the differences are still very minor.

image

rmcdermo commented 1 year ago

Definitely tighter. Thanks, Jason.

ericvmueller commented 1 year ago

Looks good to me. I think having a write-up of the enthalpy terms will also allow us to be more direct about the uncertainty introduced by the lumped approach, and then the user can decide what is acceptable. Thanks Jason!