firemodels / fds

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

Cp in condensed phase #321

Closed gforney closed 9 years ago

gforney commented 9 years ago
Please complete the following lines...

Application Version: 5.1.2 Serial
SVN Revision Number:1338
Compile Date: Feb 19, 2008
Operating System: Linux

Describe details of the issue below:

Please check out the attached plot (and input file c01.data).  I was doing
calculations in which I varied physical parameters to test the sensitivity
of the calculated burning rate and ignition time.  As one set of
calculations, I varied the SPECIFIC_HEAT.  

The different curves in the plot change only the value of the
SPECIFIC_HEAT, as 0.1, 0.5, 1.0, and 5.0 kJ/kg-K (calculation was unstable
for Cp=10, which granted, is unrealistically high);  the HEAT_OF_REACTION
is 2000 kJ/kg.

My concern is this:  the overall effective heat of gasification, Hg,eff, is
about:

Hg,eff = heat_of_reaction + Cp(Ts-Tamb) or

       = 2000 + 0.1(350-21) = 2033 kJ/kg with Cp=0.1 or
       = 2000 + 5.0(350-21) = 3650       with Cp=5.0.

Hence, Hg,eff varies by almost a factor of two.  To first order, I would
expect the steady burning rate to vary inversely with Hg,eff.  That is, I
would expect the burning rate with Cp=5.0 to be about half that with
Cp=0.1.  Referring to the figure attached, one can see that the calculation
does not give this result; rather, the mdot is about the same for these
values of Cp.  This result is contrary to my seat-of-the pants feel for
what should happen.  Any thoughts? 

Original issue reported on code.google.com by greg8394 on 2008-04-21 20:44:52


gforney commented 9 years ago
I'll assign this to Simo, and cc Jason Floyd and KM.

Original issue reported on code.google.com by mcgratta on 2008-04-21 21:27:23

gforney commented 9 years ago
Looking at the plot, the delay in ignition increases with specific heat and seems to

increase appropriately.  Suggests that the specific heat is being used OK.

Just mulling here, but once the ignition temperature is reached, the mass loss rate

is driven by the incident heat flux and the heat loss into the solid.  I would look

at the wall temperature profiles.  Conductivity is 0.1 W/m^2 (yes I rounded off). 

You can't have to of a gradient in the sample to conduct away much of the 100 kW/m^2

flux as the higher the surface temp the higher the reaction rate.  So if indeed the

wall temperature gradient is conducting only a small fraciton of the incident heat

away, then the results you have would make sense.

Original issue reported on code.google.com by drjfloyd on 2008-04-21 21:56:36

gforney commented 9 years ago
Ran your case with c=0.01 and c=1.  c=0.01 the surface temperature is 385 and c=1 it

is 387.  This temperature difference given the change in material properties does 
not seem unreasonable.  When you compute k dT/dx for the surface node (using the 
_prof.csv file) you find that with c=0.01 only about 1 % of the incident flux is 
being conducted into the solid and for c=1, 0.5 %.  There is nothing to fix here. 

Original issue reported on code.google.com by drjfloyd on 2008-04-22 13:44:59

gforney commented 9 years ago
I still think this result is not heuristically satisfying.  The Cp=5 case has an
effective Hg of about twice that of the Cp=0.1 case, yet the calculation shows
approximately the same mass loss rate and total burning time.  For Cp=5, I would have
expected about twice the burning time and half the burning rate.  To first order,
this is a B-number thing: driving force is the net heat flux (external flux minus the
losses; e.g. re-radiation, conduction, etc.  I picked a high-flux case (100 kW/m2)
so
the mass loss term would dominate).  The resistance to mass loss is the effective
heat of gasification (HEAT_OF_REACTION + Cp(Ts-T0)).  So if the resistance to mass
loss is twice, the burning rate should be half, and the burning time, twice.  

Original issue reported on code.google.com by greg8394 on 2008-04-22 16:36:17

gforney commented 9 years ago
p.s. the conductivity "losses" are a small effect for this case.  For me, a good way
to think about this is to integrate mdot over the burning time.  In that case,
conductivity losses do not matter, and the question is: should the _average_ burning
rate be lower with a higher Hg,eff?  I think they should, but the calculated result
is not showing that.  I'll come up with other ways to show this if I am not being clear.

Original issue reported on code.google.com by greg8394 on 2008-04-22 17:01:41

gforney commented 9 years ago
We're still discussing this issue. We'll keep it on hold for now.

Original issue reported on code.google.com by mcgratta on 2008-04-22 17:55:07

gforney commented 9 years ago
The basis of my argument can be seen in the simple analysis of Tewarson/Pion
("Flammability of plastics—I. Burning intensity:  Combust.Flame, 26:85-103 1976)or
of
Rhodes and Quintiere (Burning Rate and Flame Heat Flux for PMMA in the Cone
Calorimeter.NIST GCR 95-664; Thesis; 125 p. December 1994. on the NIST BFRL
Publications page:http://www.fire.nist.gov/bfrlpubs/firedoctoc.html)

Referring to Eq. 6.2 (page 47) of Rhodes, 

mdot" L = qdot, ext" + qdot, fl" - eta sigma Tv^4

mdot" is the steady-state mass-loss rate.  This is the method commonly used to
extract L (effective heat of gasification) from cone data.  In our case, the flame
radiation qdot, fl" is 0 (no flame, just an imposed flux), the qdot, ext" is so high
(100 kW/m2) that to first approximation, the re-radiation term (around 10 kW/m2) can
be neglected.  
This gives:

mdot" L = qdot, ext"

so that for a constant qdot, ext", halving L should double mdot".

To me, this part--that reducing the effective heat of gasification by two should
increase the mdot by about two--is obvious.  Perhaps the problem is really in my
input file.  Perhaps I have some mistake there.  But I am still trying to get
concurrence on what the _expected_ behaviour is.   Comments?

Original issue reported on code.google.com by greg8394 on 2008-04-24 14:50:25

gforney commented 9 years ago
FYI, Jason Floyd identified a possible bug in the "shrinking" algorithm. He and Simo

are looking it over.

Original issue reported on code.google.com by mcgratta on 2008-04-24 14:54:44

gforney commented 9 years ago
I can’t tell if Simo and Jason’s efforts with regard to the shrinking algorithm are
related to the issue I raised, or an unrelated bug.  

I'm frustrated because no-one has acknowledged the validity of what I believe to be
a
non-physical result.  That is, the response above indicates to me that you disagree
that the result should be as I suggested.    So, I’ll try again, with a few different
arguments.

1. Physical argument: 
For _steady_ gasification of an infinitely thick 1-D solid that vaporizes
around 350 C, subjected to a radiant high flux (100 kW/m2),

a.) re-radiation is a small (< 10%) loss; 
b.) conduction is not a loss term; 
c.) a c.v. around the pyrolysis shown shows that 
    heat input ~= mdot x Hg,eff
        where mdot is the mass loss rate, and Hg,eff is the effective heat
        of gasification {   delH,reaction + Cp(Treaction -Tinitial)    }
d.) at a given heat input, if you double Hg,eff, then mdot will be cut by two.

2. Consensus:
Both Takashi and Bernhard Schartel (the "German Takashi") agree that the
argument in 1.d)above explains what should happen.

3. Another calculation:
See the attached plot.  All I did was vary the Cp from 0.1 to 5.0 kJ/kg K.  The first
figure is the FDS result, the second, the results using Stas's code.  The result in
the latter is what I expected to see in the FDS result, but I did not.

If someone would agree that what Stas's calc. shows is what we
_expect_, I would feel like my message got through.  

Beyond that... the most likely source of the error here is a mistake on my
part in setting up the problem.  If someone could please take a look at my input
files (also attached), I'd appreciate it.    (I used JQ's PMMA physical property
numbers.)

Thanks,

Greg

Original issue reported on code.google.com by greg8394 on 2008-04-30 20:36:41


gforney commented 9 years ago
We're all working as fast as we can. I think you have identified a legitimate 
problem, but there is only so much time in the day. We'll use the additional input

to help debug.

Original issue reported on code.google.com by mcgratta on 2008-04-30 20:52:41

gforney commented 9 years ago
Greg, I'm sorry we can't proceed as fast as you and us wanted on this. As Kevin said,
we are running out of hours in a day. Unfortunatelly, FDS usually must wait for the
other projects to get done because there is no such a thing for me as "FDS-maintenance".

What comes to Cp problem, there seems to be two or three separate but related problems:
1) Is the formulation of chemical source term in solid solver correct? If you look
at
the tech guide, what would you say?
2) Did we implement that correctly in case of shrinking wall?
3) In case of shrinking wall, how should we distribute the energy source (sink)
between the remaining solid and the gas phase? This is slightly a problem in FDS,
because gas phase is not explicitly present in the model, and it is hard to tell if
conserve the energy or not.

Original issue reported on code.google.com by shostikk on 2008-05-01 11:45:10

gforney commented 9 years ago
Simo:

Thanks for your response.  My frustration was not with the rate of progress, but
rather, that I didn't know if anyone believed me about what should happen.  The issue
had been closed, and was only re-opened after I went down to talk to Kevin--but even
then, I didn't feel like he was convinced yet.  Hence, my further brow-beating.  I'll
look at the chemical source term in the technical guide and see if I can offer any
thoughts.

Greg

Original issue reported on code.google.com by greg8394 on 2008-05-01 13:09:14

gforney commented 9 years ago
1-c) CV around pyrolysis zone:
If we neglect re-radiation/convection and conduction into the solid, then the net 
heat flux is the incident heat flux minus the enthalpy flux of the gas leaving the

pyrolysis zone.  From discussions with Simo, I think the crux of this lies in the 
proper manner to deal with the enthalpy flux.

Original issue reported on code.google.com by drjfloyd on 2008-05-01 14:07:50

gforney commented 9 years ago
Greg,

To expand on Simo's comment 12, he had sent this to me in an email earlier on in 
this issue.

Assume we have a constant flow, adiabatic reactor where we feed in a solid at one 
end, add just enough energy to pyrolyze it so that we have gas at the pyrolysis 
temperature leaving the other end.  Being adiabatic the energy input (Q + dm/dt 
c_solid (T_amb - T_amb)) must equal the energy ouput (dm/dt * [heat of vaporization

+ c_gas (T_pyrolysis - T_amb)]).  

Therefore:
Q = dm/dt *[heat of vaporization + c_gas (T_pyrolysis - T_amb)]
Note c_solid is not seen in the equation at all.  The question is, if this is not 
correct, what is the correct form for Eq 3.92 in the Tech Guide.

Original issue reported on code.google.com by drjfloyd on 2008-05-01 17:12:43

gforney commented 9 years ago
I think that the problem you have defined has the same solution that I was talking
about above:

energy input = dm/dt x [  c_solid x (T_pyrolysis - T_amb)  +  (heat of vaporization)
] 

so it is the specific heat of the solid, not the gas.  I am not sure I understand the
first parenthetical statement above: (Q + dm/dt c_solid (T_amb - T_amb)) 
What is the term after the + sign supposed to represent?

As for Eqn. 3.92 in the Tech ref. guide, I don' t see how the last part (the integral
of the Cp term) comes into it at all; i.e., I think it should be dropped.  If the
reaction at a specific location is occurring at a given temperature (which I think
it
is), then the energy liberated (or absorbed) is just reaction rate x heat_of_reaction
, that is, the first term in the brackets in 3.92.   Comments?

Original issue reported on code.google.com by greg8394 on 2008-05-01 18:44:20

gforney commented 9 years ago
dm/dt : mass loss rate

In the reaction you take material A with c_a and with some heat of reaction convert

it to material B with c_b.  At the end of the day we need to account for the 
difference in energy content as c changes from c_a to c_b.  If that accounting does

not go into 3.92, where does it go?

Original issue reported on code.google.com by drjfloyd on 2008-05-01 19:47:09

gforney commented 9 years ago
d Enthalpy = Cp dT

so if temperature is constant (which is the case here during the reaction), there is
no sensible enthalpy change associated with any difference in Cp for the different
substances a and b.  

Take water as an example.  Cp_water is 1 cal/g-K, and Cp_steam is about 0.5 .  If we
vaporize water at around 100 C, the the heat of vaporization is about 540 cal/gm, and
that number does not depend upon the Cp of the two phases.  If were were talking
about total enthalpy, and we changed the temperatures (start and finish) away from
100 C, then we'd have to think about Cp(T1-T1), but I am considering only isothermal
examples.  Sound OK?

Original issue reported on code.google.com by greg8394 on 2008-05-01 21:13:46

gforney commented 9 years ago
In your inputs you have PMMA at 0.1 to 5 kJ/kg-K and Z_1 (the fuel gasses being 
produced) at 2.2 kJ/kg-K with the total energy requried to melt and vaporize being

unchanged.  

Original issue reported on code.google.com by drjfloyd on 2008-05-02 12:11:11

gforney commented 9 years ago
I'm not sure I follow the nomenclature your using here.  In my input files, I am using:

HEAT_OF_REACTION = 2000 kJ/kg (This is from JQ and Rhodes, p47, 2800 kJ/kg - the
sensible part which is about 800 kJ/kg)

SPECIFIC_HEAT = 0.1 to 5.0 kJ/kg-K (A typical value is 2.2 kJ/kg-K, also from JQ and
Rhodes)

(I did not specify ...(1) in either of these, but I think it should defaul to 1.

What is your Z_1 ?

Original issue reported on code.google.com by greg8394 on 2008-05-02 13:16:54

gforney commented 9 years ago
At the end of the day energy must be conserved.

If I have a material with a c=5, a heat of reaction of 2, and a pyrloysis 
temperature of 350 C, then at the momment of pyrolysis we have total enthalpy of 
(350 - 20)*5 [enthalphy of solid] - 2 [heat of reaction] = 1648.  This implies the

gas should have a c of 4.99.  

If I have a material with a c=0.1, a heat of reaction of 2, and a pyrloysis 
temperature of 350 C, then at the momment of pyrolysis we have total enthalpy of 
(350 - 20)*0.1 [enthalphy of solid] - 2 [heat of reaction] = 31.  This implies the

gas should have a c of 0.094.

In the mixture fraction, the presumed c of the fuel is ~2.2 at 350 C.

Any inconsistency in this is being seen in the energy deposity/removed from the 
solid phase.  So for your c=5 case, the gas enthalpy << enthalpy of the solid at 
pyrolysis and the extra energy stays in the solid and is available for further 
pyrolysis.  I think this is what is causing what you are seeing.  The question is 
then, what changes need to be made to 3.92 or how we define inputs?

Original issue reported on code.google.com by drjfloyd on 2008-05-02 13:38:46

gforney commented 9 years ago
Jason:

I'll copy your comment here so I can comment more easily on each part. 

>>At the end of the day energy must be conserved.
>>
Maybe yes, maybe no.  It depends upon where you draw your control volume.  For
example, if the burning polymer is within your c.v., then the total enthalpy (which
includes the chemical energy in the polymer) will be conserved as the polymer
vaporizes.  On the other hand, if only the surface is on the edge of your c.v., then
mass addition will be like an energy addition, and H will not be conserved--it is
being added to the system.

>>If I have a material with a c=5, a heat of reaction of 2, and a pyrloysis 
>>temperature of 350 C, then at the momment of pyrolysis we have total enthalpy of

>>(350 - 20)*5 [enthalphy of solid] - 2 [heat of reaction] = 1648.  This implies the

>>gas should have a c of 4.99.  

I think this part above is incorrect; the gas Cp does not come into play until you
heat the gas up above Tpyrolysis.  Cp is just a variable which describes the increase
in H as T goes up (P constant); Cp =dH/dT|P,const.  So if the gas is at Tpyrolysis,
Cp gas does not matter yet.  It's easiest to cast all this in terms of the enthalpy
of each constitutent; Cp is just the convenient parameter that we can use to
calculate the variation in H(T,Xi) as T goes up.

Also, "total enthalpy" is always relative: you need to define a reference temperature
(this varies according to the thermodynamic data).  

>>If I have a material with a c=0.1, a heat of reaction of 2, and a pyrloysis 
>>temperature of 350 C, then at the momment of pyrolysis we have total enthalpy of

>>(350 - 20)*0.1 [enthalphy of solid] - 2 [heat of reaction] = 31.  This implies the

>>gas should have a c of 0.094.

This is incorrect for the same reason given above.  

>>In the mixture fraction, the presumed c of the fuel is ~2.2 at 350 C.

>>Any inconsistency in this is being seen in the energy deposity/removed from the 
>>solid phase.  So for your c=5 case, the gas enthalpy << enthalpy of the solid at

>>pyrolysis and the extra energy stays in the solid and is available for further 
>>pyrolysis.  I think this is what is causing what you are seeing.  The question is

>>then, what changes need to be made to 3.92 or how we define inputs?

I don't follow this part.  Maybe we have to agree on the above pre-requisites first.

Original issue reported on code.google.com by greg8394 on 2008-05-02 14:40:12

gforney commented 9 years ago
No matter what, any alorigthm we implement needs to conserve energy.  In the case of

pyrolysis any energy lost/gained by the solid phase must be balanced by the opposite

in the gas phase.  You cannot isolate one from the other in the software, they are

coupled.  

How can the gas Cp not come into play?  When you pyrolyze a solid you produce gas at

the pyrolysis temperature.  Are you suggesting that the gas has no enthalpy at the

momment of pyrolysis? No, it has enthalphy based on the pyrolysis temperature and 
the reference temperature (yes we know enthalpy is relative, look at Eq 3.92 and you

will see the reference temperature).

One must consider the total system.  If we take a solid and turn it into a gas, all

the energy must be accounted for in the underlying algorithm.  If I have 1 kg of 
fuel gas at 350 C in FDS, that has an enthalpy (assuming 20 C ref T) of ~726 kJ/kg.

At 350 C, the PMMA as you have defined it varies between 33 kJ/kg and 1650 kJ/kg. 

You have defined a heat of reaction of 2000 kJ/kg.  If the heat of reaction is added

to the solid enthalpy (as in your H2O example), note that in no case do the energy

balances agree.  We conserve energy, so energy must be made available from 
somewhere, and in the alogrithm that somewhere is the solid.  So again the question

to you is, what should the algorithm be?  The solid phase equations as we use them

are in the manual, the numerical method is in the manual, if these are incorrect how

should they be corrected?  

Original issue reported on code.google.com by drjfloyd on 2008-05-02 14:57:04

gforney commented 9 years ago
In the new release of FDS 5.1.5, I added a parameter called ADJUST_HOR to the MISC 
line. We cannot change the code permanently until we upgrade to 5.2, as this 
parameter will change the functionality. But for the moment, if you want to prevent

FDS from adjusting the HEAT_OF_REACTION with differences in specific heat between 
reactants and products, set ADJUST_HOR=.FALSE.

Original issue reported on code.google.com by mcgratta on 2008-05-08 22:08:43

gforney commented 9 years ago
I get results similar to what Greg got with Stas/FAA's model when cp is varied--see

01.pdf. 

I think of the heat of reaction (or heat of volatilization or vaporization, 
depending on the terminology) as the difference in total enthalpy (chemical plus 
sensible enthalpy, referenced to the same temperature datum) between the condensed-
phase at Tp and the volatiles at Tp. If you accept this definition, then the heat of

reaction/volatilization/vaporization already includes cp differences between the gas

and the solid, so it is not necessary to include a cp term in the HOR. 

Original issue reported on code.google.com by chris.lautenberger on 2008-05-09 22:17:14


gforney commented 9 years ago
Greg,

Since FDS 5.2, the default operation of the condenced phase model is to use simply
heat_of_reaction, not to add any c_p integrals to it.

Could you confirm, that you now get results that are consistent with intuition and
other existing models? Thanks.

Simo

Original issue reported on code.google.com by shostikk on 2008-06-12 11:44:53

gforney commented 9 years ago
Simo:

The attached file shows the result if I use the new keyword I use the keyword that
Kevin told me about: ADJUST_HOR=.FALSE..  The top plots shows the mass loss rate as
a
function of time with varying Cp (0.1, 0.5, 1, and 5): FDS on the left, Stass/FAA on
the right;  the bottom plot shows the two together.  The results with FDS make
physical sense now, and are much closer to Stas’s results.

Greg  

Greg  

Original issue reported on code.google.com by greg8394 on 2008-06-12 14:12:43


gforney commented 9 years ago
(No text was entered with this change)

Original issue reported on code.google.com by shostikk on 2008-06-13 11:01:23

gforney commented 9 years ago
In SFPE.Handbook.of.Fire.Protection.Engineering.-.3Ed.-.2002, p.2-233, I found two 
equations about heat conduction within solid phase which involve the enthalpy of 
volatiles. 
the equations and my comments are shown in the attached file.

Original issue reported on code.google.com by youngewong on 2008-08-18 10:45:42


gforney commented 9 years ago
Consider a gas-phase reaction -- the heat of reaction is defined as the enthalpy 
difference between the gaseous products and the gaseous reactants under ISOTHERMAL

conditions, right? The same applies here, except the reactant is solid and the 
products are gaseous (and sometimes also solid, i.e. charring reactions). 

Those equations in HOR.doc are mass conservation and energy conservation with a 
source term -- what do propose the source term should be?

Original issue reported on code.google.com by chris.lautenberger on 2008-08-19 14:36:39

gforney commented 9 years ago
I think the source term is heat of pyrolyzation.

Original issue reported on code.google.com by youngewong on 2008-08-20 00:21:42