Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Way to freeze flame spread similar to freezing time for ramp input? #13234

Open
mcgratta opened this issue Jul 26, 2024 Discussed in #13222 · 61 comments
Open

Way to freeze flame spread similar to freezing time for ramp input? #13234

mcgratta opened this issue Jul 26, 2024 Discussed in #13222 · 61 comments
Assignees

Comments

@mcgratta
Copy link
Contributor

Discussed in #13222

Originally posted by bwklein July 25, 2024
I'm looking for a way to freeze the 'spread' of a fire growth that is based on the SPREAD_RATE parameter.
Similar to how we can freeze the time value input to a ramp for HRRPUA, is there a way to do this with the input to the spreading function. Essentially, stopping the spread of the fire, but keeping the cells active that are within the spread radius.

@mcgratta mcgratta self-assigned this Jul 26, 2024
@mcgratta
Copy link
Contributor Author

Bryan -- we're going to implement a way to control spreading fires with more features than just freezing it. But if you look at the test case called Fires/circular_burner.fds, does this do what you want?

@mcgratta
Copy link
Contributor Author

mcgratta commented Aug 6, 2024

I looked into this issue some more. Anything can be done with spread rates, but it would be more costly than what we do now because if we make SPREAD_RATE a time-dependent function, then we have to evaluate all the WALL, CFACE and PARTICLE cells each time step and change the ignition time if the fire front has swept by. Is there a strong use case for this? If you just want to stop the spread, you can just limit the VENT on which the SPREAD_RATE is specified.

@drjfloyd
Copy link
Contributor

drjfloyd commented Aug 6, 2024

I think Bryan was looking for a solution to do this when you don't know a priori how far it will spread.

@obscureed
Copy link

Sorry to hijack the thread, but I have some features I would like to see in fires with SPREAD_RATE:
• I have tried to combine SPREAD_RATE to define growth with RAMP_Q to define a sudden shutdown of the fire later. This does not work; a slow spread currently causes a slow shutdown. The memory used to set up the spread interferes with the connection of each facet to the RAMP_Q times. I would like some way to have a spreading fire that later shuts down cleanly.
• SPREAD_RATE will achieve"t-squared" growth (where the total heat release rate is proportional to time^2), but only for a circular fire. It would be convenient to create t-squared growth with non-circular shapes, using a growing, connected collection of facets.
• The initial part of a spread-rate growth curve is usually zero heat output, unless the centre of the fire bed is exactly a facet centroid. To some extent, this is just a consequence of a coarsely resolved growth, but it would be nice to have an option for some heat output from the start.
• If the centre of the fire bed is somewhere symmetric in the grid, then the growth curve goes up in jumps of 4 or 8 facets, all at the same distance from the centre. It would be nice to have a more smoothly resolved growth curve.

You may have thought of all these already. Thanks for reading!
Ed

@mcgratta
Copy link
Contributor Author

Let me push back a bit here. Currently, the user specifies a constant SPREAD_RATE originating from a single point on a VENT. The fire spreading stops when it reaches the boundary of the VENT. One can even make the VENT itself circular, as in the example case Fires/circular_burner.fds. This feature is cheap and easy to implement because during the simulation set up, all the relevant surface cells are assigned an appropriate ignition time and a ramp of heat release rate to invoke at the time of ignition.

However, all of the suggestions put forth above will cost more because they require that surface cells be scanned at each time step to determine ignition and various other behaviors. Not a lot of cost, but more than before for a feature that may or may not be widely used. In my opinion, the idea of the "t-squared" fire has been extended far beyond its original intent because of regulatory requirements and general practice. But nevertheless, we implemented what I have described above as a convenience for end users who need to conform to these requirements. But now, much like the layer height concept, this concept is beginning to go beyond what is called for by basic practice or regulatory requirements. I worry that if we add all these new variations, we risk confusing real fire spread behavior with these ad hoc mathematical shapes. Sure, we can make the fire spread like an oval rather than a circle, but where does that end? I would prefer that if fire spread is to be modeled in a more realistic way that there be some rationale for it. For example, we use a more sophisticated spread model for wildland applications, the so-called level set method, but this is based on the actual physics of spread over complex terrain and different fuel types. We are very cautious about making the fire spread in a way that cannot be justified other than as some hypothetical construct.

@obscureed
Copy link

You are right that I am driven to t-squared growth by regulatory definitions. I do not have a problem with that: we need to distinguish between "slow" and "fast" growth situations, and t-squared growth factors are a shared language that defines growth. (It also has some basis in reality, for a constant spread rate over a circular fire bed.) It does not seem unreasonable to me to want to combine t-squared growth with a later rampdown (for example because the fire is extinguished). Currently I do this by using SPREAD_RATE to define the growth and by creating an OBST on a timer over the fire bed to kill it suddenly, effectively like dropping a drain cover onto a fire bed. This does the job.

One route to satisfy the features I mentioned would be a connection between RAMP_Q and facet-by-facet timing. (Currently there is a connection between t-squared growth rate and facet-by-facet timing, or between an arbitrary RAMP_Q and variable Heat Release Rate Per Unit Area. The current implementation of RAMP_Q by using variable HRRPUA is less realistic than anything I have suggested.)
• The way I envisage it, the RAMP_Q time variations of heat release rate could be discretised to changes in quanta matching the number of facets. The facets of the fire bed would be put into a fixed order -- for example, sorted by distance from some coordinate, with a tie-breaker based on angle. If the ramp function requires N of the facets to be active, then the first N on the sorted list are active. If necessary, all this can be pre-calculated as a time history for each facet. The number of on/off flips might be greater than currently for some ramps, but there is no scanning of distances or angles in mid simulation.
• If I want to, I can then program in t-squared growth, plateau and ramp-down into a RAMP_Q. With a bit of tweaking, I can define the RAMP_Q to approximate t-squared growth using at least one facet from the start.
• I cannot comment on how widely this would be used, but I know I would use it a lot.

@mcgratta
Copy link
Contributor Author

If you make 'RAMP_Q' a step function whose duration is dictated by the fuel load (mass per area), then you have a ring of fire that expands outwards until it reaches the boundary of the VENT and goes out. If you make the duration longer than the time to reach the boundary, you will get t-squared growth, steady burning, and then a ramp down of some sort. Why is this not adequate to the purpose?

@drjfloyd
Copy link
Contributor

Wtih DEVC and CTRL functions you should be able to script out any desired rate of growth over any arbitrary set of VENT or OBST including decay. For example, for oval spread you could define rings of VENTs each with a different SURF_ID and RAMP_Q.

@obscureed
Copy link

@mcgratta : I am not completely sure I understand what you're suggesting. If you're suggesting a combination of SPREAD_RATE=... to define growth and a RAMP_Q to define extinguishing, then the reason it is not adequate is that it does not work. The events defined in the RAMP_Q are smeared out by the growth timescale defined in the SPREAD_RATE. (My simplistic way of thinking about it is that FDS has reset the timers of all the facets of the fire bed, so they respond late to the RAMP_Q.) Maybe this is a bug -- it certainly surprised me the first time I saw it.

Here is an example -- I'll put the code at the end of this post. This is intended to give "fast" t-squared growth (as defined in UK regulations), which should reach 1000 kW in 146 seconds. Then a plateau, then sudden extinguishing, ramping down from 100% at t=300s to 0% at t=310s. This is what happens instead:
graphHRR

There are several aspects I would like to change about this HRR curve:
• The sudden extinguishing has been smeared out badly.
• There is zero HRR until after t=20s.
• The HRR growth goes up in leaps of 4 or 8 facets, which is the same as 4% or 8% of the eventual fire size.
— If the discretisation of the HRR curve used individual pixels, each 1% of the eventual fire size in this case, then the first facet might turn on after a few seconds, which would be an improvement. I would prefer at least one facet to be active from the start. I can (and I do) achieve this by adjusting the position of the fire bed so that the centroid exactly matches a facet centroid, but this is fiddly.
• The t-squared growth that I would like is initially achieved, but it slows down when the spread approaches the corners of the square fire bed. To reach the full HRR takes 168s instead of 146s. This is certainly not a bug -- spread rate is doing what was specified -- but I would like to have a way to specify t-squared, facet-by-facet growth for non-circular fire beds.
• The overshoot on the first leap at t=21s is a bit strange, but that's just numerical stuff and not a priority.

@drjfloyd -- Yes, I could define rings of VENTS to define ovals. I could even define a different RAMP_Q for every facet. However, I still think it is the kind of feature that the software could take care of productively.

Here's the example code. It runs in a few minutes:

&HEAD CHID= 'firegrowth_fast', TITLE= 'firegrowth_fast'  /

! A fire (1.0m*1.0m square, 1.0MW) in an open box.
! This combines SPREAD_RATE corresponding to "fast" growth (as in BS7974-1:2019, Table 2)
!           and RAMP_Q      intended to extinguish the fire suddenly (over 10s) at t= 300s [but SPREAD_RATE interferes]

&TIME T_END= 500., RESTRICT_TIME_STEP= .FALSE./
&DUMP NFRAMES= 250 /
&MATL ID= 'CONCRETE',  FYI= 'Typical Concrete', CONDUCTIVITY= 1.4, SPECIFIC_HEAT= 0.84, DENSITY= 2100.0 /
&SURF ID= 'WALL', RGB= 200,210,220, MATL_ID= 'CONCRETE',  THICKNESS= 0.25, DEFAULT= .TRUE.   /

&MESH XB =   -2.000,   2.000,  -2.000,   2.000,   0.000,   4.000, IJK= 20,20,20, MPI_PROCESS= 0, ID= 'mesh1' /

&REAC ID= 'FUEL', C= 1.0, H= 2.0, O= 0.5, SOOT_YIELD= 0.100, CO_YIELD= 0.040, HEAT_OF_COMBUSTION= 20000.000 /
&SPEC ID= 'SOOT' /

&SURF ID= 'fire',  HRRPUA= 250.000 COLOR= 'RED',  RAMP_Q= 'fire_ramp',  /
&VENT XB =   -1.000,   1.000,  -1.000,   1.000,   0.000,   0.000, SURF_ID= 'fire', SPREAD_RATE= 0.00773578 /

&RAMP ID= 'fire_ramp', T=    0.000, F=    0.000, /
&RAMP ID= 'fire_ramp', T=    1.000, F=    1.000, /
&RAMP ID= 'fire_ramp', T=  300.000, F=    1.000, /
&RAMP ID= 'fire_ramp', T=  310.000, F=    0.000, /

&VENT XB =   -2.000,  -2.000,  -2.000,   2.000,   0.000,   4.000, SURF_ID= 'OPEN' /
&VENT XB =    2.000,   2.000,  -2.000,   2.000,   0.000,   4.000, SURF_ID= 'OPEN' /
&VENT XB =   -2.000,   2.000,  -2.000,  -2.000,   0.000,   4.000, SURF_ID= 'OPEN' /
&VENT XB =   -2.000,   2.000,   2.000,   2.000,   0.000,   4.000, SURF_ID= 'OPEN' /
&VENT XB =   -2.000,   2.000,  -2.000,   2.000,   4.000,   4.000, SURF_ID= 'OPEN' /
&SLCF PBX=    0.000, QUANTITY= 'TEMPERATURE', ID= 'x0' /
&TAIL /

@mcgratta
Copy link
Contributor Author

What you are doing is what I suggested, even though it is not exactly what you want. To do what you want would require making this feature far more complicated and prone to error and CPU bloat. I am uncomfortable going to all this effort to replicate a fictitious and unnatural heat release rate curve. If you want to mimic exactly any HRR curve, do not use the spread feature and just ramp up the fire over a given area to exactly match the curve.

@obscureed
Copy link

"Not exactly what you want" is an understatement -- I specified the fire to die out by t=310s, but it is still at full power at t=320s and it does not die completely until after t=470s. I could argue with your other points as well, but for now I'll be quiet.

@drjfloyd
Copy link
Contributor

You specified a RAMP that stops 310 s after the wall cell ignites. The T in RAMP_Q is time since ignition not overall simulation time. My guess is the last wall cell is igniting at 160 s (160 + 310 = 470).

What you are asking for is to have multiple control function inputs governing just the t2 growth vent type plus a bunch of logic to deal with a list of facets and times. This is not as trivial as change as you might think. We implemented the detailed control function logic in part so that user's can implement more complex needs without us having to develop new code, make verificaiton test cases for new code, document new code, and maintain new code in perpetuity.

@rmcdermo
Copy link
Contributor

rmcdermo commented Aug 21, 2024

Discussed this a bit and we think we can make a fairly simple change to how ramps are processed that would move us in the right direction. The key issue is decoupling the "T" in the RAMP from the "TSI" or Time Since Ignition in the ramp that multiplies the mass flux at a wall cell. One thought is to have a logical on the RAMP line that tells FDS whether T is meant to be TIME or TSI. Another option would be to use different parameters (T for TIME and TSI for Time Since Ignition---this is my preference but might not be very backward compatible).

If I make this hack in the code with your simple ramp,

&RAMP ID='fire_ramp', T=  0.0, F=1.0/
&RAMP ID='fire_ramp', T=300.0, F=1.0/
&RAMP ID='fire_ramp', T=310.0, F=0.0/

then I get this for the result (with twice the resolution; getting rid of the initial spike is also resolution dependent, but as you said we can tackle that later).
Figure_1

In general, the ramp function can be set to handle different spread rates and vent areas as follows. Suppose we have a vent with area A_VENT and a spread rate S. The circular area of the spread is A_SPREAD = pi*(S*T)^2. If F is the ramp function we would use for a t^2 fire with A_VENT, then the new ramp function is simply,

F_NEW(T) = F(T) * A_VENT/MIN(A_VENT,A_SPREAD[T])

Until A_SPREAD>A_VENT, this is just a constant that depends on the spread rate (if S is used for F[T] the constant is just 1) and beyond that you get back the original F(T).

I'm not sure if this solves all problems, but it seems to me it would get us moving toward what you are asking for. Let us know if this is worth implementing.

@drjfloyd
Copy link
Contributor

drjfloyd commented Aug 21, 2024 via email

@obscureed
Copy link

A few comments:

• I do not understand the usefulness of a new parameter SPREAD_DISTANCE, or new limits on the spread area. If the user wants a circular fire bed, smaller than the XB patch, this can already be specified in &VENT RADIUS=....

• "The T in RAMP_Q is time since ignition not overall simulation time."
— Yes, that is how RAMP_Q currently behaves when combined with SPREAD_RATE, but it is not what I expected from the User Guide, and not what I wanted.
— On the other hand, I can believe that it is exactly what some users might want -- for example, if the fuel at each facet has a finite life and burns with a certain history regardless of its neighbors.
— So now there is a bit of untangling to do. I cannot see a completely clean solution: either you break backward compatibility (where previous input files that used T in RAMP_Q to mean TSI will behave differently) or you break consistency (where T in a ramp usually means simulation time, except when applied to a RAMP_Q of a fire-bed VENT with SPREAD_RATE). Either way is going to be a bit ugly, and the User Guide will have to be updated.

• I accept @drjfloyd 's point that this is combining two controls that affect fire size -- the facet-by-facet spread, and the ramp -- and this is potentially complicated.
— I do not see how the existing logic controls currently help here.
— I have not completely worked out what @rmcdermo 's suggestion does, but it certainly looks useful. Basically, I think this is assuming a fixed rule to combine the two controls, and the rule is multiplication? Or the spread rate defines which facets are active, and the ramp (as a function of simulation time) defines the HRRPUA that applies to active facets?

• One point I want to make is that all this is helping to reach more realistic fire models without massive user input.
— For a fire size that starts small and grows slowly to large, SPREAD_RATE just looks right. The simple alternative is to have the large fire bed struggling along with a small HRRPUA spread uniformly across it while the fire size is small -- and this does not resemble a real fire at all well. (The complicated alternative is for the user to build a patchwork of fire beds, each with different ramps and fire sizes -- this is really fiddly to set up. To check when the combined patchwork reaches a certain size, you would need a spreadsheet.)
— None of the current discussion helps the 4- and 8-facet jumps in HRR. I still think there is potential for facet-by-facet response to an arbitrary RAMP_Q, without needing geometric scanning in mid-simulation. This is a separate feature from SPREAD_RATE, so perhaps it can wait, though it would be a good fit for grow-plateau-extinguish timelines.

• One extra feature that would be convenient would be a freeze on spread rate, in the same way that a NO_UPDATE_DEVC_ID can be applied to a ramp. I think that's a feature request for another day, but I mention it now in case it needs to be factored in.

@rmcdermo
Copy link
Contributor

@drjfloyd It depends on what you mean by "correct". For the period when the A_S < A_V, the fire grows with t^2 simply by following F=1, with T in the ramp being TIME. If you want each newly activated cell to have its own ramp based on TSI, then my approach will not work. But I think that using TIME instead of TSI gives the tight control over the shutdown that Ed is after. So, it's a bit of a compromise: easy coding and controls shutdown, but it does not solve all problems.

@drjfloyd
Copy link
Contributor

@obscureed you can define mutiple vents each with their own set of ramps and if you wanted those ramps to experience decay based on some event in FDS, that can be also be done using control functions. I recognize that this is more work for the user, but it can be done and gives much greater flexibility.

@rmcdermo At the instant of spread to the next cells in the current approach, the effective diameter of the fire does not change. We get a smooth change in the effective physical size of the fire which with the heat release rate drives entrainment. With your proposed approach we will a see jump up in the size of the fire diameter and a jump down in the burning intensity when each new ring ignites. This is going to cause some kind of step change in the entrainement. I see this change as potentially taking us further from the concept of what a t2 fire is (which I recognize is in itself just a convenient way of looking at real fires on solid fuels). Is this change in behavior going to be a net good or net bad for the typical user of SPREAD_RATE?

@rmcdermo
Copy link
Contributor

"With your proposed approach we will a see jump up in the size of the fire diameter and a jump down in the burning intensity when each new ring ignites."

I don't think we are on the same page. Taking T to be TIME in the RAMP, yes, with each new ring we get an increase in size, that is just how a t^2 fire grows. But we will not get a reduction in mass flux of fuel at the new wall cell. In this case F is 1 and the HRRPUA is a constant.

@drjfloyd
Copy link
Contributor

I think we are not on the same page. Say I am using a RAMP plus SPREAD_RATE to give a t2 growing fire. At some time t I have 100 kW and 0.1 m2 at t and at t+dt jumps to the next ring of cells giving 0.2 m2 and grows slightly (we generally have small timesteps) so the new HRR is 101 kW. My HRRPUA has gone from 1000 kW/m2 to 505 kW/m2 and the effective fire diameter has jumped by a factor of 2 over a timestep. In the old method, the original 0.1 m2 would keep the 1000 kW/m2 and the new would start at 5 kW/m2. My effective fire diameter and burning intensity will hardly change.

@rmcdermo
Copy link
Contributor

rmcdermo commented Aug 22, 2024

With the idea I'm proposing T_IGN is used for activation of a wall cell, but not for the RAMP of Q. I probably should not have complicated the idea with a discussion of rescaling F, forget about that for now. All F does is go to zero at a specific time (when Ed wants to turn off the fire). While the A_S(T) = pi * (S * T)^2 is less than the XB area of the VENT (assume a huge VENT area for now), each active cell has an HRRPUA of whatever is on the SURF line. The fire grows in size only by an increase in the number of active wall cells (which increases the AREA by t^2).

@drjfloyd
Copy link
Contributor

drjfloyd commented Aug 22, 2024 via email

@rmcdermo
Copy link
Contributor

Yes, for now at least, in my approach the step changes are resolution dependent. As I said, my method does 2 things: it's easy to code and it tightly controls the shutdown of the fire.

@mcgratta
Copy link
Contributor Author

mcgratta commented Sep 9, 2024

Jason, can you give a simple example of your approach. Does it require any code change?

@drjfloyd
Copy link
Contributor

drjfloyd commented Sep 9, 2024

By my approach do you mean using multiple vents and ramps?

Fast growing fire (1 MW in 150 s) over an 0.2 m x 0.2 m burner at 0.1 m grid resolution. This gives two rings of grid cells. The inner ring gives 40 kW and grows from 0 to 30 s, The outer ring gives 120 kW and grows from 30 s to 60 s.

&SURF ID='FIRE1',HRRPUA=1000,RAMP_Q='FIRE1'/
&SURF ID='FIRE2',HRRPUA=1000,RAMP_Q='FIRE2'/

&VENT XB=-0.1, 0.1,-0.1, 0.1, 0.0, 0.0, SURF_ID='FIRE1'/
&VENT XB=-0.2,-0.1,-0.1, 0.1, 0.0, 0.0, SURF_ID='FIRE2'/
&VENT XB= 0.1, 0.2,-0.1, 0.1, 0.0, 0.0, SURF_ID='FIRE2'/
&VENT XB=-0.2, 0.2, 0.1, 0.2, 0.0, 0.0, SURF_ID='FIRE2'/
&VENT XB=-0.2, 0.2,-0.2,-0.1, 0.0, 0.0, SURF_ID='FIRE2'/

&RAMP ID='FIRE1',T=0,F=0.00/
&RAMP ID='FIRE1',T=10,F=0.11/
&RAMP ID='FIRE1',T=20,F=0.44/
&RAMP ID='FIRE1',T=30,F=1.00/

&RAMP ID='FIRE2',T=30,F=0.00/
&RAMP ID='FIRE2',T=40,F=0.26/
&RAMP ID='FIRE2',T=50,F=0.59/
&RAMP ID='FIRE2',T=60,F=1.00/

You could add DEVC/CTRL_ID or DEVC/CTRL_DEP_ID to either or both RAMPs and/or DEVC/CTRL_ID to the VENTs to get whatever complex behavior you wanted on top of the t^2 growth.

@mcgratta
Copy link
Contributor Author

mcgratta commented Sep 9, 2024

Referring to this section of the User's Guide

image

could we use the output of a device with SPATIAL_STATISTIC='AREA' and use that in conjunction with a user-specified RAMP of the desired total HRR and apply the quotient to any burning cell?

@drjfloyd
Copy link
Contributor

drjfloyd commented Sep 9, 2024

That is a good idea.

You could integrate BURNING RATE or HRRPUA over the VENT with SPREAD_RATE using a small, non-zero QUANTITY_RANGE(1) to get the burning area. Then have a custom control with a time input that gives the total HRR over time. Then another control function to divide by the area to get HRRPUA. Then for your SURF for the VENT you could have HRRPUA=1 and for the RAMP_Q use CTRL_ID_DEP and give it the HRRPUA result.

For more complex behavior (like turning off the fire) you could have intermediate control functions between the HRR function and the division function that adjust the HRRPUA. For example if you had a CTRL function to turn the fire off that went from true with fire on to false with fire off, then you could output that CTRL using DEVC with QUANTITY=CONTROL which outputs 0 false and 1 true. Then you would add a function that multiplies the HRR by the DEVC before the division function.

@mcgratta
Copy link
Contributor Author

mcgratta commented Sep 27, 2024

I am trying to implement this in a test case. I am stuck. If I create a RAMP of the global desired HRR, how can I put the dependent variable of the ramp, say F, into a CTRL line with FUNCTION_TYPE='DIVIDE'? I have the fire area, but I need the output of the HRR ramp to somehow serve as the first INPUT_ID.

@drjfloyd
Copy link
Contributor

Have a CUSTOM CTRL that is the fire RAMP. The output of that CTRL will be the fire size at that time and can be used as an input for another control

@drjfloyd
Copy link
Contributor

Something like this plus DEVC/CTRL logic to get BURNING AREA:

&SURF ID='FIRE',HRRPUA=1,RAMP_Q='HRRPUA RAMP'/
&RAMP ID='HRRPUA RAMP',T=0,F=0,CTRL_DEP_ID='HRRPUA CTRL'/

&CTRL ID='HRRPUA CTRL',FUNCTION_TYPE='DIVIDE',INPUT_ID='FIRE SIZE','BURNING AREA'/

&CTRL ID='FIRE SIZE',FUNCTION_TYPE='CUSTOM', INPUT_ID='TIME',RAMP_ID='FIRE SIZE'/
&RAMP ID='FIRE SIZE',T=0,F=0/Fire size in kW as a function of time
&RAMP ID='FIRE SIZE',T=1,F=1/
&RAMP ID='FIRE SIZE',T=2,F=4/

&DEVC XYZ=...,QUANTITY='TIME',ID='TIME'/

@drjfloyd
Copy link
Contributor

In terms of the AREA function, it doesn't operate on pixels in smokeview. It operates on grid cells in FDS. It is a good point that the safest way to doing this would be to ensure that you could never have a division by zero which would be good to add to the example. Could change:

&CTRL ID='HRR', FUNCTION_TYPE='CUSTOM', INPUT_ID='TIMER', RAMP_ID='HRR' /

to

&CTRL ID='HRR-RAMP', FUNCTION_TYPE='CUSTOM', INPUT_ID='TIMER', RAMP_ID='HRR' /
&CTRL ID='HRR', FUNCTION_TYPE='SUM', INPUT_ID='HRR-RAMP','CONSTANT',CONSTANT=1E-10/ Prevent divide by zero error

@obscureed
Copy link

Sorry, yes, I was using "pixels" as a shorthand for "grid-cell facets" -- nothing to do with graphics or smokeview. The &VENT is given an XYZ coordinate at the corners of some grid cells, and SPREAD_RATE activates facets on the basis of their centroids, so no facets are active initially.

The confusing part for me is that the protection from division by zero is not needed here. (I no longer think that the answer can be that the error value is detected and automatically reset to zero. If that was the case, then the initial Face Area would be zero, so HRRPUA CTRL would be zero, so Face Area, as calculated from HRRPUA, would never go above zero even when spread occurred.) I do not have an explanation for that.

I agree that it is cleaner to specify only the lower limit of QUANTITY_RANGE. However, changing the lower limit to a tiny value does not remove the spike in HRRPUA at low times (around t=4s). This spike is a consequence of the inputs: the area of the fire is proportional to t^2 initially (to within facet resolution), but the HRR target is proportional to t initially, based on the first two values in the HRR ramp. So the HRRPUA, calculated from (HRR/area) actually has a (1/t) asymptote at small times if we ignore facet resolution. The first few values are saved from being even larger by whatever copes with zero facets.

@drjfloyd
Copy link
Contributor

The plot I attached was running with QUANTITY_RANGE(1)=1E-10 and FDS 6.9.1.

@drjfloyd
Copy link
Contributor

I guess we don't need the div zero protection in the input file. We have this in the source:

CF%INSTANT_VALUE=0._EB
...
IF (ABS(DV%SMOOTHED_VALUE)>TWO_EPSILON_EB) CF%INSTANT_VALUE = CF%INSTANT_VALUE / DV%SMOOTHED_VALUE

So when the area is zero, the function just gets set to zero.

@drjfloyd
Copy link
Contributor

In this input file we have the added factor that the fire plateaus at 100 s but the spread rate is 0.1 m/s on a 6 m x 6 m square vent. So the fire finishes spreading at ~30 s. If you wanted to really mock up a t^2 spreading fire, you could use a circular VENT and set the spread rate to get the desired growth time to the fire plateau. However, if you want to guarantee the HRR while having the area grow as close as it can to a spreading circle, you cannot expect a clean looking HRRPUA curve. The VENT area can only increase stepwise in integer multiples of the area of a grid cell.

If the VENT was:

&VENT SURF_ID='FIRE', XB=1.0,7.0,1.0,7.0,0.0,0.0, XYZ=4.0,4.0,0.0, SPREAD_RATE=0.03, RADIUS=3 /

then HRRPUA plot would look more like a constant HRRPUA; however, early on as the fire spreads you will have large relative jumps in area which are going to make the HRRPUA curve have spikes. There is no way to avoid this as we can only set boundary conditions on the entire grid cell. As the fire grows, each new ring of grid cells adds less relative area and the HRRPUA curve approaches a constant value.

image

@obscureed
Copy link

TLDR: To avoid the HRR values being permanently incorrect, make sure that the area of your requested fire bed exactly matches the area that is achieved as an integer number of facets.

I have found one pitfall with this method. This pitfall occurs if the area of the fire bed reached by the spreading fire does not equal the area implied in the &VENT definition of the fire bed. The consequence is that the plateau maximum fire size is permanently wrong.

This is easiest to see with a square fire bed. Here are three alternative definitions of a fire bed:

&VENT SURF_ID='FIRE', XB=2.200,4.400,2.200,4.400,0.000,0.000, XYZ=3.3,3.3,0.0, SPREAD_RATE=0.0155563 /
&VENT SURF_ID='FIRE', XB=2.101,4.499,2.101,4.499,0.000,0.000, XYZ=3.3,3.3,0.0, SPREAD_RATE=0.0155563 /
&VENT SURF_ID='FIRE', XB=2.299,4.301,2.299,4.301,0.000,0.000, XYZ=3.3,3.3,0.0, SPREAD_RATE=0.0155563 /

With cell size of 0.2m (starting at XYZ=0.,0.,0.), these three definitions will all result in the same number of facets in the fire bed (because the spreading fire includes all the facets whose centroids fit inside the definition). The achieved area of the fire bed is 4.84m^2, which is 121 pixels each 0.2m*0.2m. This area, A_achieved, exactly equals the requested area, A_requested, for the first definition above. For the second and third definitions, however, A_requested is 5.75m2 and 4.01m2 respectively.

The problem comes because FDS automatically applies a correction factor, AREA_ADJUST, when A_achieved does not equal A_requested. Suppose you define a fire bed with a constant HRRPUA in a conventional way:

&VENT SURF_ID='FIRE', XB=2.101,4.499,2.101,4.499,0.000,0.000, XYZ=3.3,3.3,0.0, SPREAD_RATE=0.0155563 /
&SURF ID='FIRE', HRRPUA=347.8 /

When FDS notices that A_achieved=4.84 does not equal A_requested=5.75, then it applies a correction factor to the requested HRRPUA -- in this case, HRRPUA_applied will be higher than HRRPUA_requested, by a factor of (5.75/4.84). The assumption is that the user cares more about the total HRR than about HRRPUA, so HRRPUA_applied is adjusted to give the correct total HRR. This causes an inaccuracy in the current technique, because we calculate 'Fire Area' from A_achieved, but we use it to calculate HRRPUA_requested. The graph below is for the larger XB, and it shows that HRR is incorrect by about 20% compared to the intended fire size of 2000kW. (This is the worst-possible XB for the fire bed, which has coarse but not horrifying resolution: 11*11 cells. It is a smaller effect for finer cell sizes.)

One way to avoid this pitfall is to make sure that A_requested equals A_achieved. For a rectangular fire bed, this is easy enough: make sure that the XB coordinates match grid coordinates (that is, facet edge coordinates). For a circular fire bed, you can adjust RADIUS so that the area (pi*RADIUS^2) exactly equals the achieved area (which is an integer multiple of the facet area) -- but there is no easy formula for the number of facets that will fit inside RADIUS, so this takes some trial and error.

graphHRR

@obscureed
Copy link

Here is my latest version of the setup, using CTRL arithmetic to define t^2 growth. The graphs at the end show that HRR is lower than the intended 2000 kW by about 3%. The HRRPUA judders around significantly, as @drjfloyd predicted.

&HEAD CHID='spread2', TITLE='Sample fire spread case' /

&MESH IJK=20,20,20, XB=0.0,4.0,0.0,4.0,0.0,4.0, MULT_ID='mesh' /
&MULT ID='mesh', DX=4.0, DY=4.0, DZ=4.0, I_UPPER=1, J_UPPER=1, K_UPPER=0 /

&TIME T_END=200. /

&DUMP DT_DEVC=0.01, DT_HRR=0.01 /

&REAC FUEL='N-HEPTANE', SOOT_YIELD=0.02 /

&VENT MB='XMIN', SURF_ID='OPEN' /
&VENT MB='XMAX', SURF_ID='OPEN' /
&VENT MB='YMIN', SURF_ID='OPEN' /
&VENT MB='YMAX', SURF_ID='OPEN' /
&VENT MB='ZMAX', SURF_ID='OPEN' /

!! New fire: final size Q=2000 kW, final radius r=1.3 m, "ultrafast" growth, alpha=0.1876 kW/s^2
!! Hence final HRRPUA H=Q/(pi*r^2)=376.698 kW/m^2
!! Hence spread rate v=sqrt(alpha/(pi*H))=0.0125906 m/s
!! Unnecessary to specify XYZ at the centre of the square face, but it reminds us it is at a facet centroid.
!! RADIUS has not been adjusted; therefore its area does not exactly match the area of the facets that it contains
!! -- therefore total HRR will be slightly wrong at all times
!! -- A_requested = pi*1.3^2 = 5.3093 m^2
!! -- A_achieved =             5.4800 m^2 (from 137 facets, each 0.2*0.2) -- no simple formula to predict this
!! -- so HRR will be lower than the intended values by a factor (5.3093/5.4800)
!! -- eg peak 1938 kW instead of 2000 kW, 3% error
&VENT SURF_ID='fire', XB=2.0,4.6,2.0,4.6,0.0,0.0, XYZ=3.3,3.3,0.0, RADIUS=1.3, SPREAD_RATE=0.0125906 /

&SURF ID='fire', HRRPUA=1., RAMP_Q='ramp_hrrpua' /
&RAMP ID='ramp_hrrpua', T=0, F=0, CTRL_ID_DEP='ctrl_hrrpua' /

&DEVC ID='TIMER', XYZ=1,1,1, QUANTITY='TIME' /

&CTRL ID='ctrl_hrr_growth', FUNCTION_TYPE='MULTIPLY', INPUT_ID='CONSTANT','TIMER','TIMER', CONSTANT=0.1876 /
&DEVC ID='devc_hrr_growth', QUANTITY='CONTROL VALUE', CTRL_ID='ctrl_hrr_growth', XYZ=1,1,1 /

&CTRL ID='ctrl_hrr_limit', FUNCTION_TYPE='CUSTOM', INPUT_ID='TIMER', RAMP_ID='ramp_hrr_limit' /
&RAMP ID='ramp_hrr_limit', T=  0., F= 2000.0 /
&RAMP ID='ramp_hrr_limit', T=150., F= 2000.0 /
&RAMP ID='ramp_hrr_limit', T=180., F=    0.0 /

&CTRL ID='ctrl_hrr', FUNCTION_TYPE='MIN', INPUT_ID='ctrl_hrr_growth','ctrl_hrr_limit' /
&DEVC ID='devc_hrr', QUANTITY='CONTROL VALUE', CTRL_ID='ctrl_hrr', XYZ=1,1,1 /

&DEVC ID='FireArea', QUANTITY='HRRPUA', QUANTITY_RANGE(1)=1e-10,
   XB=0,8,0,8,0,1, SURF_ID='fire', SPATIAL_STATISTIC='SURFACE AREA' /

&CTRL ID='ctrl_hrrpua', FUNCTION_TYPE='DIVIDE', INPUT_ID='ctrl_hrr','FireArea' /
&DEVC ID='devc_hrrpua', QUANTITY='CONTROL VALUE', CTRL_ID='ctrl_hrrpua', XYZ=1,1,1 /

&BNDF QUANTITY='BURNING RATE', CELL_CENTERED=T /
&DEVC XB=1.0,7.0,1.0,7.0,0.0,0.0, QUANTITY='HRRPUA', SPATIAL_STATISTIC='MAX', ID='hrrpua_max' /

&TAIL /

h__HRR

g__hrrpua

@drjfloyd
Copy link
Contributor

drjfloyd commented Oct 1, 2024

Something isn't right in the wall boundary condition. The control logic is use the desired t^2 fire and divide by the current area to get the HRRPUA needed to achieve that fire size. Whatever a 1.3 m radius winds up being on that grid is going to wind up being the area used in the control logic so the fire size should be right. If you look at the BNDF data you will see that for some reason the burning rate is not the same over the four meshes. It should be. Something is amiss.

spread2_0518

@drjfloyd
Copy link
Contributor

drjfloyd commented Oct 1, 2024

For some reason the B1%AREA_ADJUST is not the same for the four meshes and it should be. This is resulting in the lower burning rate on meshes 2,3, and 4.

@drjfloyd
Copy link
Contributor

drjfloyd commented Oct 1, 2024

@mcgratta in this case what is happening is the circular vent is split over four meshes with most of it in mesh 1. VT%INPUT_AREA and VT%FDS_AREA are being computed mesh-by-mesh. In mesh 4 (upper right in image) the circle just barely clips the mesh so the ratio of the clipped area to that of a single grid cell area is low. In mesh 1, the circle is mostly in the mesh and the clipped area vs grid are is close so the ratio is close to 1. If RADIUS isn't used on the VENT line, then this works as both the INPUT_AREA and FDS_AREA wind up the same. Instead of what we do now in READ_VENT, should we just compute INPUT_AREA using the XB or RADIUS, on each mesh compute FDS_AREA, gather all the FDS_AREA to get the full VENT FDS area, and then push back out to all meshes?

@mcgratta
Copy link
Contributor Author

mcgratta commented Oct 1, 2024

Is this logic needed only for a spreading fire, or is this a problem with vents split across meshes in general?

@drjfloyd
Copy link
Contributor

drjfloyd commented Oct 1, 2024

Looks like this would occur any time you had a vent split unequally over multiple meshes when RADIUS was being used. The AREA_ADJUST factor is not correct in that case which gets applied to any mass flux boundary condition.

@drjfloyd
Copy link
Contributor

drjfloyd commented Oct 1, 2024

&HEAD CHID='vent_error', TITLE='Sample fire spread case' /

&MESH IJK=20,20,20, XB=0.0,4.0,0.0,4.0,0.0,4.0, MULT_ID='mesh' /
&MULT ID='mesh', DX=4.0, DY=4.0, DZ=4.0, I_UPPER=1, J_UPPER=1, K_UPPER=0 /

&TIME T_END=10. /

&REAC FUEL='N-HEPTANE', SOOT_YIELD=0.02 /

&VENT MB='XMIN', SURF_ID='OPEN' /
&VENT MB='XMAX', SURF_ID='OPEN' /
&VENT MB='YMIN', SURF_ID='OPEN' /
&VENT MB='YMAX', SURF_ID='OPEN' /
&VENT MB='ZMAX', SURF_ID='OPEN' /

&VENT SURF_ID='fire', XB=2.0,4.6,2.0,4.6,0.0,0.0, XYZ=3.3,3.3,0.0,RADIUS=1.3/

&SURF ID='fire', HRRPUA=1000,COLOR='RED' /

&BNDF QUANTITY='BURNING RATE', CELL_CENTERED=T /

&TAIL
vent_error_0643

@mcgratta
Copy link
Contributor Author

mcgratta commented Oct 1, 2024

Should we just add an IF RADIUS>0 after this block in read.f90:

               ! Vent area

               IF (ABS(XB_USER(1)-XB_USER(2))<=SPACING(XB_USER(2))) &
                  VT%UNDIVIDED_INPUT_AREA = (XB_USER(4)-XB_USER(3))*(XB_USER(6)-XB_USER(5))
               IF (ABS(XB_USER(3)-XB_USER(4))<=SPACING(XB_USER(4))) &
                  VT%UNDIVIDED_INPUT_AREA = (XB_USER(2)-XB_USER(1))*(XB_USER(6)-XB_USER(5))
               IF (ABS(XB_USER(5)-XB_USER(6))<=SPACING(XB_USER(6))) &
                  VT%UNDIVIDED_INPUT_AREA = (XB_USER(2)-XB_USER(1))*(XB_USER(4)-XB_USER(3))

               IF (RADIUS>0._EB) VT%UNDIVIDED_INPUT_AREA = PI*RADIUS**2

               IF (ABS(VT%X2-VT%X1)<=SPACING(VT%X2) ) VT%INPUT_AREA = (VT%Y2-VT%Y1)*(VT%Z2-VT%Z1)
               IF (ABS(VT%Y2-VT%Y1)<=SPACING(VT%Y2) ) VT%INPUT_AREA = (VT%X2-VT%X1)*(VT%Z2-VT%Z1)
               IF (ABS(VT%Z2-VT%Z1)<=SPACING(VT%Z2) ) VT%INPUT_AREA = (VT%X2-VT%X1)*(VT%Y2-VT%Y1)

               ! Check the SURF_ID against the list of SURF's

@mcgratta
Copy link
Contributor Author

mcgratta commented Oct 1, 2024

But it appears we're doing something like that

      SELECT CASE(ABS(VT%IOR))
         CASE(1)
            DO K=K1+1,K2
               DO J=J1+1,J2
                  IF ( VT%RADIUS>0._EB) THEN
                     VT%INPUT_AREA = VT%INPUT_AREA + CIRCLE_CELL_INTERSECTION_AREA(VT%Y0,VT%Z0,VT%RADIUS,Y(J-1),Y(J),Z(K-1),Z(K))
                     IF (((YC(J)-VT%Y0)**2+(ZC(K)-VT%Z0)**2)>VT%RADIUS**2) CYCLE
                  ENDIF
                  VT%FDS_AREA = VT%FDS_AREA + DY(J)*DZ(K)
               ENDDO
            ENDDO

@drjfloyd
Copy link
Contributor

drjfloyd commented Oct 1, 2024

We compute INPUT_AREA in two locations. The block you show and again in at Line 12217 in read.f90. That second location is where we get the intersection of the circle with the grid to get INPUT_AREA and sum up the grid cells touched by the circle to get FDS_AREA. Then in init.f90 we divide the two to get the area adjust. I think for all vents we would want the case that if the vent is being split over meshes to get the input area as specified on &VENT and have FDS_AREA be the area over all meshes. This might also impact cases without RADIUS where you have a VENT that crosses meshes of different resolutions or without being grid algined. I'll try to make a test case to see if that is also a problem.

@drjfloyd
Copy link
Contributor

drjfloyd commented Oct 1, 2024

It seems to work OK without RADIUS and when the mesh changes resolution.

@mcgratta
Copy link
Contributor Author

mcgratta commented Oct 2, 2024

Jason, in your example above, where the HRRPUA is not uniform over the four meshes, it appears that the desired HRR is achieved. So I'm not in favor of trying to "fix" this. Reopen if you think there is still some issue to deal with.

@mcgratta mcgratta closed this as completed Oct 2, 2024
@drjfloyd drjfloyd reopened this Oct 2, 2024
@drjfloyd
Copy link
Contributor

drjfloyd commented Oct 2, 2024

This means we have a dependence in the fire symmetrey based on meshing which I don't think is a desireable behavior, and the HRR was not correct in the case attached by @obscureed. It was low because of this issue.

@drjfloyd
Copy link
Contributor

Have had a chance to get back to this. A few things:

For a circular VENT an approximate method was being used to get the area of the unsnapped circular VENT that interesects both the VENT XB and the MESH XB. For mass flux the error was was small but for volume flow this was not being done correctly for the case where the VENT XB does not completely cover the circle. One fix for this is to just do the analytic function for the intersection. There are six general cases for the overlap and with a little geometry and trig you can get a solution for each case.

A second issue is that each mesh is computing the adjustment factor independent of the other meshes. As seen in my simple test case above, if the meshes do not symetrically divide the circle you will get different burning rates on different meshes even if the same grid resolution is being used. We get the correct total HRR, but there is going to be some degree of mesh dependence in the plume behavior which I think is undesirable. The solution for this is an MPI exchange of the per mesh presnapped area of circular VENTs. This takes one MPI call after reading the inputs. With that we can ensure a circular vent uses the adjustment factor for each grid cell in the vent.

The third issue and the real cause for why the CTRL case was not getting the correct HRR is how a VENT with a mass flux bc is treated. FDS adjusts the mass flux based on the VENT XB vs. how the VENT actually snaps to the grid. So if you give a 1000 kW/m2 HRRPUA to 1 m2 XB VENT, but the VENT snaps to 1.2 m2, FDS reduces the HRRPUA to 1000 x 1/1.2. This way you get the correct HRRPUA. In this case as you observed the difference in VENT area was the same as the difference in the FDS ouptut. The CTRL logic computes the burning rate for the area of the grid cells that are burning, but when that burning rate is applied it then gets adjusted by the vent snapped vs unsnapped area which in this case decreases the HRRPUA a little. We are correcting a calculation that no longer needs correction. One solution for doing this in general is you could run without SPREAD_RATE on the VENT and look at what the vent area is for the FireArea device then adjust the RAMP formula to increase/decrease the HRRPUA accordingly. Another solution is we could add a flag to VENT to disable area adjustment. I think the flag is better. It is a lot easier and it might be useful to have this type of control anyways for other cases where people might use CTRL functions with actual FDS areas to set boundary conditions.

With putting in all these changes, when I run the case with the proposed AREA_ADJUST=F on the VENT I get the correct HRR (plotted here as burning rate * HoC) and the BNDF shows uniform burning. The spikes in the HRR appear to be an order of operations thing. The fire spreads to new wall cells using the prior burning rate giving a step increase during the timestep. Then at the end of the timestep the DEVC and CTRL see the new cells and adjust the burning rate for the new area bringing it back to the curve for the next timestep. In the actual realized HRR curve those spikes get lost as they are smoothed out via the mixing process.

image
spread2_0729

@obscureed
Copy link

This is good stuff, @drjfloyd -- thanks. When the area correction factor is applied, it should be the same for all meshes, so the MPI call is needed. I agree with you that the AREA_ADJUST flag could be useful on other occasions (for example, when a user wants to specify HRRPUA and actually get that HRRPUA).

One question, though: aren't you introducing an error by using analytical formulas for circle/VENT overlap, rather than summing the area of the snapped-on facets? (I'm thinking about a fire bed where RADIUS is specified and area adjustment is active.) You calculate the area adjustment factor using A_requested and A_analytical, but this factor is later applied to A_snapped, so the HRR delivered to the model will be incorrect by a factor (A_analytical/A_snapped). It sounds like the analytical method is an improvement on the previous approximation, especially for edge cases with strange overlap or non-overlap, but it still has this error. There was an example earlier in the thread, with RADIUS=1.3 on a mesh resolution of 0.2, where this error was 3%.

@drjfloyd
Copy link
Contributor

FDS is still computing both the true area and the faceted area. The AREA_ADJUST flag fixes the HRR error for the CTRL function case earlier in the thread.

@mcgratta
Copy link
Contributor Author

Whatever changes need to be made to the code, thoroughly test them. I have spent years debugging these kinds of "adjust" factors. They have many unexpected side effects.

@drjfloyd
Copy link
Contributor

That is why I haven't pushed anything yet.

I've tested the analytic cirlce area calc with test cases that hit each of the 15 total variants of the 6 general cases (1 variant for no overlap, 1 variant for circle > rectangle, 1 variant for rectangle > circle, 4 ways for each of the 1,2, or 3 corners inside [which quadrant the excess sits in]).

The adding a flag to bypass adjustment for mass flux by itself should not impact default behavior.

It is the last change of ensuring uniform flux when a circular vent is split over meshes that could cause problems. I need to test this more and run it through firebot once enough space opens up on our cluster.

drjfloyd added a commit to drjfloyd/fds that referenced this issue Dec 20, 2024
drjfloyd added a commit that referenced this issue Dec 20, 2024
FDS Source: Fix double adjust issue in Issue #13234 by MPI exchange of VENT areas
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
Labels
None yet
Projects
None yet
Development

No branches or pull requests

5 participants