firemodels / fds

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

Suggestion: Global vent/leak parameter for pressure zones #10968

Closed FireResearch-BK closed 1 year ago

FireResearch-BK commented 2 years ago

This is an enhancement suggestion stemming from issues with pressure zones / enclosed compartments.

Situation: We get many CAD models from clients, which are not built manifold/solid. These are imported into Pyrosim and ultimately exported as FDS input files. This may generate lots of thin-walled, enclosed compartments, resulting in a lot of pressure zones (in some cases more than 700). These enclosed compartments/zones cause all kinds of pressure-related issues during simulation runs when close to the fire, causing frequent crashes due to numerical instability.

Attempted workarounds: After months of reading the user guide, technical reference guide, user group and bug tracker posts, as well as countless hours of experimentation (not to mention the frustration at seeing a simulation crash after two weeks of computation), we are aware of the following workarounds to deal with pressure zones and enclosed compartments:

Suggestion: Implement an optional global flag in FDS that automatically deals with excessive pressure in zones or enclosed compartments. Here's some ideas:

Other suggestions to deal with this issue on a global basis in future FDS releases are of course also welcome.

P.S A simple test case with thin-walled enclosed compartments and an exaggerated fire near them that can be used for experimentation is attached below. It crashes at around 118 seconds in FDS 6.7.9 for me and even disabling pressure zones causes density warnings (which are the first indication of an impending instability and crash) near the end of the simulation.

Test_Case.txt

emanuelegissi commented 2 years ago

I confirm the issue. The suggested solutions would prevent many headaches. Thank you. Emanuele Gissi

Il gio 29 set 2022, 11:24 IFAB-BK @.***> ha scritto:

This is an enhancement suggestion stemming from issues with pressure zones / enclosed compartments.

Situation: We get many CAD models from clients, which are not built manifold/solid. These are imported into Pyrosim and ultimately exported as FDS input files. This may generate lots of thin-walled, enclosed compartments, resulting in a lot of pressure zones (in some cases more than 700). These enclosed compartments/zones cause all kinds of pressure-related issues during simulation runs when close to the fire, causing frequent crashes due to numerical instability.

Attempted workarounds: After months of reading the user guide, technical reference guide, user group and bug tracker posts, as well as countless hours of experimentation (not to mention the frustration at seeing a simulation crash after two weeks of computation), we are aware of the following workarounds to deal with pressure zones and enclosed compartments:

  • Rebuild geometry in simpler fashion purely with obstacles -> Unfeasible because of the necessary extra time.
  • Manually close resulting pressure zones by obstruction -> Zones may be identified by FDS and viewed in Smokeview, but there is no way to dump a human-readable list of zones (with cell locations) from FDS, which makes this approach a time intensive trial & error guessing game.
  • Manually assign vents and leakage paths to zones to simulate venting or bursting -> Unfeasible due potentially very high to number of zones and the aforementioned lack of human-readable data related to zone locations.
  • Force thickening of geometry when converting into blocks (in Pyrosim) -> Sometimes helps, sometimes doesn't and won't prevent zones from being created in the first place as well.
  • Disable pressure zones altogether or use an alternate pressure solver -> Won't prevent numerical instability due to excessive pressure in enclosed compartments as well.
  • Material burn-away -> The most logical solution to open enclosed compartments, but very hard to implement, increases computation time by an unfeasible amount, is not required for most of our simulations, and won't prevent crashes from these pressure issues either.

Suggestion: Implement an optional global flag in FDS that automatically deals with excessive pressure in zones or enclosed compartments. Here's some ideas:

  • "&PRES ZONE_LEAKAGE=0.00001/": Assign a global leakage rate to the ambient zone (zone 0) to all zones.
  • "&PRES ZONE_EQUALIZE=T/": Keep zone pressure equal to the pressure of any surrounding ambient cell at all times.
  • "&PRES ZONE_AUTOVENT=(P_THRESH,dP_VENT)/": Fake the pressure-related cracking or bursting of zones. P_THRESH is a pressure threshold and dP_VENT is the subsequent pressure loss rate per time step (or second?) once pressure has exceeded P_THRESH. The enclosed compartment would then vent until equal with ambient (zone 0) pressure and remain there for the rest of the simulation.

Other suggestions to deal with this issue on a global basis in future FDS releases are of course also welcome.

P.S A simple test case with thin-walled enclosed compartments and an exaggerated fire near them that can be used for experimentation is attached below. It crashes at around 118 seconds in FDS 6.7.9 for me and even disabling pressure zones causes density warnings (which are the first indication of an impending instability and crash) near the end of the simulation.

Test_Case.txt https://github.com/firemodels/fds/files/9673141/Test_Case.txt

— Reply to this email directly, view it on GitHub https://github.com/firemodels/fds/issues/10968, or unsubscribe https://github.com/notifications/unsubscribe-auth/ACXMVY7ZNAT4765SD27HVILWAVN3VANCNFSM6AAAAAAQYSF74M . You are receiving this because you are subscribed to this thread.Message ID: @.***>

mcgratta commented 2 years ago

I'll take a look at this. I have had similar thoughts. The problem is usually that some solutions violate both the physics and/or numerics with unexpected consequences. The idea of default leakage might be a good solution because it does have physical and numerical validity.

drjfloyd commented 2 years ago

A leakage area default should be tied to the 2/3 power of the zone volume. Something like 1E-5*VOL**TWH. Otherwise a default area that works well for a larger zone might not work well for a smaller zone and vice versa.

mcgratta commented 2 years ago

In this test case, you can add

&PRES SOLVER='ULMAT' /

and the case will run. Why? The default pressure solver in FDS does not strictly enforce a zero normal velocity at internal obstructions. Thus, in your case, a small amount of gas penetrates the obstruction above the fire, eventually causing the density to rise to a fictitiously high value and FDS shuts down. The ULMAT solver enforces the no-flux solid boundary condition. So why don't we use it all the time -- cost. It is more expensive than the default solver. Depending on the case and mesh size, it can be much more. For a trivial test case like this, it is affordable. For more explanation, search the User's Guide for the word ULMAT.

All this being said, we need to do a few things

  1. The default leakage is still probably a good idea for cases where the ULMAT solver might not work or be too expensive.
  2. Improve guidance on when to switch to an alternative pressure solver. It is not obvious which one to use for a given scenario. The default solver is always the cheapest (so far), but there might be other considerations, like actually finishing the calculation.
drjfloyd commented 2 years ago

Item 1 will take some restructuring of HVAC. Right now we read ZONE just before HVAC so HVAC knows all the ZONE leakage as it goes to read in HVAC inputs and set up the HVAC networks. We'd have to push defining zone leakage until after we find all the pressure zones and then reallocate all the HVAC arrays to handle things that weren't defined during READ_HVAC.

mpachera commented 2 years ago

I had the same issue which I reported to the user's group and at that time the solution I found was to automatically detect all the zones and their phase through a 3D slice file and then to create a small script what was filling automatically all the unwanted zones which weren't already solid. Such scripting is not complicated, however it required to run the simulation for a fraction of time and to (automatically) process the different results. The other drawback is that the users still has to decide the zones to keep. Do you think this might be a viable solution? Matteo

FireResearch-BK commented 2 years ago

Thank you for all the replies so far! I'm really glad that the issue is being looked at because it's been racking our brains for months now.

 

Related question: Until there's an auto-leakage function in FDS, I'm going to try to create a script for generating "&ZONE" entries from FDS2ASCII output data, but I'm slightly confused why the output file contains zones with intermediate indices. For the test case above, for example, the geometry (and FDS) basically defines four pressure zones (indexed as 0, 0.5, 1.0 and 1.5), but the output file from FDS2ASCII contains zone indices that lie in between (e.g. 0.125, 0.25, 0.375, 0.75). If I was to define zones from that output file, I'd assign those intermediate values to the closest "regular" pressure zone, but what am I supposed to do with the indices like 0.25 and 0.75 that are exactly in between two major indices? Assign them to both? Disregard them?

 

In this test case, you can add

&PRES SOLVER='ULMAT' /

and the case will run. Why? The default pressure solver in FDS does not strictly enforce a zero normal velocity at internal obstructions. Thus, in your case, a small amount of gas penetrates the obstruction above the fire, eventually causing the density to rise to a fictitiously high value and FDS shuts down. The ULMAT solver enforces the no-flux solid boundary condition. So why don't we use it all the time -- cost. It is more expensive than the default solver. Depending on the case and mesh size, it can be much more. For a trivial test case like this, it is affordable. For more explanation, search the User's Guide for the word ULMAT.

All this being said, we need to do a few things

1. The default leakage is still probably a good idea for cases where the `ULMAT` solver might not work or be too expensive.

2. Improve guidance on when to switch to an alternative pressure solver. It is not obvious which one to use for a given scenario. The default solver is always the cheapest (so far), but there might be other considerations, like actually finishing the calculation.

We've tried ULMAT for at least one problematic case before, which, however, couldn't prevent a crash. But then again, the materials in that case were set to combustible, so I assume that the associated hot gas release within a burning thin walled, enclosed compartment caused its own excessive pressure problems.

I can try to run another problematic case (without combustible materials) with ULMAT and see if it gets the simulation any farther.

 

I had the same issue which I reported to the user's group and at that time the solution I found was to automatically detect all the zones and their phase through a 3D slice file and then to create a small script what was filling automatically all the unwanted zones which weren't already solid. Such scripting is not complicated, however it required to run the simulation for a fraction of time and to (automatically) process the different results. The other drawback is that the users still has to decide the zones to keep. Do you think this might be a viable solution? Matteo

Are you referring to sl3d_read.py? We've tried to modify that script to generate the zone definitions with leak areas for the FDS input file, but even for simple cases didn't get very far because due to (admittedly) insufficient understanding of the output data format (section 25.9 of the User Guide didn't help that much). Since then, we found out that FDS2ASCII also does the trick, albeit less elegantly with a few more manual steps in between. We just need to come up with some script to process that data.

drjfloyd commented 2 years ago

By default, SLCF are interpolated to cell corners. So 0.125 means one of the eight grid cells for that corner is zone 1 and the other 7 are zone 0. Try adding CELL_CENTERED=T which turns off the interpolation.

mpachera commented 2 years ago

The output from the sl3d_read is a 4D matrix (x,y,z,time), you can forget about the last axis and you end up with a 3D mapping of the zone ids and phases. The FDS2ASCII simply takes the matrix and reshape it so it's plotted on a csv file, but the structure of the file behind is the same. From there if you have the mesh size and resolution you can reconstruct the blocks you need to fill the holes. You'll need some trials with the indexes because it works differently for cases CELL_CENTERED=T or CELL_CENTERED=F.

FireResearch-BK commented 2 years ago

So with the data obtained from FDS2ASCII, I've tried to determine coordinates for zone definitions which worked out as I could use my calculated values for placing valid zones in Pyrosim. As the FDS User Guide states that each manually defined pressure zone also needs a leakage path, I've created the necessary air leak surfaces and applied them to a single of the six walls of each compartment. FDS, however, complains that the surface can not be applied to a thin obstruction when trying to launch the modified test case. Unless I did something fundamentally wrong (see the attached, modified test case), I assume this is simply a limitation within FDS 6.7.9?

Test_Case.txt

mcgratta commented 2 years ago

Yes, leakage cannot be directly applied to a thin obstruction. A hack would be to leak through a thickened portion of the boundary. My idea for "automatic" leakage would probably involve a single Lagrangian particle that injects or removes gas to retain ambient pressure and density.

drjfloyd commented 2 years ago

Rather than make particle which would just wind up changing the D_SOURCE and M_DOT_G_PPP arrays we could just directly adjust those arrays. Or make a D_SOURCE_ZONE and M_DOT_G_PP_ZONE if we wanted to track this seperately.

mcgratta commented 2 years ago

Yes, that is a good idea. We'd probably want to track this separately.

FireResearch-BK commented 2 years ago

Update on ULMAT: Started running the problematic case, but stopped it after a few days because the time required to compute a single time step increased by a minimum factor of 10, making this no option at all.

Since finding out about #10982, I'm pretty sure that the two issues are related. The surface properties in the problematic case were using "BACKING='INSULATED'", which, combined with thin-walled enclosed compartments and heat exposure on the outside, throw the pressure solver in the resulting pressure zone off by so much that even &PRES VELOCITY_TOLERANCE=0.00005,MAX_PRESSURE_ITERATIONS=50/ can only delay the inevitable simulation instability.

I'm running the problematic case with 'BACKING='EXPOSED' on all surfaces at the moment as I hope that the zone pressure issue can be worked around if heat can be freely transported through obstacles.

mcgratta commented 2 years ago

There are some thin-walled obstructions for which only an exact solver like ULMAT will work because the penetration error must be exactly zero. If ULMAT is too expensive, the thin walls should be thickened or the hollow space should be filled. At the moment, we have no alternative, cheaper pressure solver.

FireResearch-BK commented 2 years ago

Thickening is a quick and global fix, but only works if the result won't block any critical flow paths.

Is VELOCITY_TOLERANCE=0.0 a feasible, computationally cheaper alternative to ULMAT or does its value need to be greater than zero, negating any error accumulation?

mcgratta commented 2 years ago

It would take infinitely long to drive the error down to zero. If thinness is important and you cannot fill the void and the volume of the ZONE is relatively small, there is no alternative at the moment than the ULMAT solver.

mcgratta commented 1 year ago

I made a slight change to the way thin obstructions are treated. You can try you case again with these builds.

FireResearch-BK commented 1 year ago

Running a case that is prone to crash with 6.7.9 with 6.7.9-933 at the moment; will report once it crashed or finished.

FireResearch-BK commented 1 year ago

While the case has progressed well beyond the previous point of crashing, the "computing time per timestep" factor is considerably poorer with 6.7.9-933 than with 6.7.9, up to a factor of 2 (i.e. cases take twice as long as before). The number of pressure solver iterations for each time step has also increased.

Does this mean that FDS' nightly by default tries harder to minimize velocity errors near thin obstacles by running more solver iterations?

I'd like to compare the nightly's FFT performance with ULMAT, but 6 x ~160000 cells appear to be too much for FDS, causing a forrtl related crash before the first time step.

mcgratta commented 1 year ago

I will test the two versions.

mcgratta commented 1 year ago

Can you post the input file of the case you are timing.

FireResearch-BK commented 1 year ago

I sadly can't post it here as it contains client data. Is there a way I can send it to you confidentially?

mcgratta commented 1 year ago

Anything you send me is public because I work under the Freedom of Information Act in the U.S. Can you scrub the case and send a simpler version?

FireResearch-BK commented 1 year ago

Sigh. Fine. See the attachment. This is a case that crashes after ~40 seconds with the regular 6.7.9 release due to density errors and runs very poorly with 6.7.9-933 because the FFT solver iterates itself to death.

Case.txt

mcgratta commented 1 year ago

I'm running the cases. Might take a while to complete. One suggestion is to ease up on the VELOCITY_TOLERANCE. The default value is half the cell size, which in your case is about 0.02 m/s. You could try 0.01 or 0.005. The value you have chosen is so tight that it makes the pressure solver hit 100 iterations each time step.

mcgratta commented 1 year ago

I noticed that we changed a numerical parameter since our last official release. If you add

SUSPEND_PRESSURE_ITERATIONS=T

on the PRES line, FDS will stop iterating the pressure solver if progress is too slow. In 6.7.9-933, by default, FDS is going to iterate for as many times as you specify with MAX_PRESSURE_ITERATIONS.

All this being said, I am running your test case in a variety of ways. They still haven't reached the point of failure.

FireResearch-BK commented 1 year ago

One suggestion is to ease up on the VELOCITY_TOLERANCE. The default value is half the cell size, which in your case is about 0.02 m/s. You could try 0.01 or 0.005. The value you have chosen is so tight that it makes the pressure solver hit 100 iterations each time step.

The tight values are a force of habit as it was the only way to stabilize a simulation in 6.9.7 and before.

 

I noticed that we changed a numerical parameter since our last official release. If you add

SUSPEND_PRESSURE_ITERATIONS=T

on the PRES line, FDS will stop iterating the pressure solver if progress is too slow.

Testing this at the moment. It does seem to help somewhat.

In 6.7.9-933, by default, FDS is going to iterate for as many times as you specify with MAX_PRESSURE_ITERATIONS.

I'm not seeing that in any .out file of any case I ran with 6.7.9-933 so far (running without "SUSPEND_PRESSURE_ITERATIONS"). The "Pressure Iterations" counter for each time step mostly remains below MAX_PRESSURE_ITERATIONS.

All this being said, I am running your test case in a variety of ways. They still haven't reached the point of failure.

I'm running my case in 6.7.9-933, too using your suggestions for velocity tolerance and a maximum of 50 pressure iterations, but I'm not near the point of previous failure yet. With the less tight values, things proceed at a more normal rate though.

mcgratta commented 1 year ago

I started the case 48 h ago and it has run 60 s. I am not setting a VELOCITY_TOLERANCE. I am using the default.

The code commit that is referenced above changes the way that we handle the mass flux at thin obstructions that border pressure zones. I believe that this change is preventing the numerical instability, at least so far. I will continue running the case. I am also running the case with the new code and a tight tolerance, and as you observed, it is about twice as slow.

FireResearch-BK commented 1 year ago

I started the case 48 h ago and it has run 60 s. I am not setting a VELOCITY_TOLERANCE. I am using the default.

My case has been running for 30 hours now (50 max iterations, suspend iterations enabled and 0.005 velocity tolerance) and is at ~52 seconds, which is past the former crash point. So far, so good.

mcgratta commented 1 year ago

My computer is a bit old.

mcgratta commented 1 year ago

The case for which I am using the default velocity tolerance finished successfully.

mcgratta commented 1 year ago

The case with the very tight tolerance has run out 56 s.

FireResearch-BK commented 1 year ago

All 4 variations of the case (1483 pressure zones each) I threw at 6.7.9-933 using &PRES MAX_PRESSURE_ITERATIONS=50, VELOCITY_TOLERANCE=0.005, SUSPEND_PRESSURE_ITERATIONS=T/ completed successfully.

Another case using a similar setup (1080 pressure zones) and the same &PRES line also completed successfully.

Temperature plots at the devices also do not look too different from the runs with 6.7.9.

So good work on the change!