search

firemodels / fds

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

FDS6 transport scheme and extra species #846

Closed gforney closed 9 years ago

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

Application Version: FDS5.4.1
SVN Revision Number: 4708
Compile Date: 10 Sep 2009
Operating System: Linux

This issue is related to issue 864.

If you run the attached cup burner case for one second of real time, and 
look at the helium slice, you will see helium in the flame zone where 
there should be none. 

I tested the various FDS6 options individually, and found that with 
FLUX_LIMITER of -1 (the FDS5 default) all works fine, but changing that to 
4 produces the funny result.

The attached case is 2D cylindrical, but I observed the same effect with 
2D and 3D cartesian runs.

Original issue reported on code.google.com by jukka.vaari@vtt.fi on 2009-09-17 09:26:35

gforney commented 9 years ago 
Thanks -- Randy will take a look at this.

Original issue reported on code.google.com by mcgratta on 2009-09-17 11:50:22

gforney commented 9 years ago 
The FLUX_LIMITER is a red herring.  The case also works with DNS=F or with DNS=T but
N_REACTIONS=0 (I commented out the call to COMBUSTION).  The problem is somewhere in
the new combustion routines.  I have no idea why things still work with FL=-1.

What specifically goes wrong is that the diffusion flux for FX=(1,J,K) becomes
nonzero (note that FX(0,...) is the boundary value) and this is because DYDX becomes
nonzero, and this is because YY(1,J,K,N=HE) becomes finite after passing through the
combustion routine.  Oddly, this problem is somehow localized to the first off-wall
gradient.

This is as much as I have discerned at this point.  I am copying Jason on this in
case he has any ideas.

Original issue reported on code.google.com by randy.mcdermott on 2009-09-17 18:19:08

gforney commented 9 years ago 
Does this problem occur with combustion2 = .FALSE.?  
and the problem is that Y(:,:,:,HE) is 0 before combustion2 and then nonzero after?

COMBUSTION2 hasn't been officially released yet, in fact the version in the
repository looks nothing at all like my current working copy (unfortunately that
version overhauls other global array structures so that is definitely not ready for
use).  My guess is there is a pointer / index error somewhere.  I'll spend a few
moments to see if I notice something obvious.

Original issue reported on code.google.com by drjfloyd on 2009-09-17 18:41:48

gforney commented 9 years ago 
Yes, still a problem with COMBUSTION2=F.

Original issue reported on code.google.com by randy.mcdermott on 2009-09-17 18:49:12

gforney commented 9 years ago 
hmm....well that is interesting.  I'll dig into it.

Original issue reported on code.google.com by drjfloyd on 2009-09-17 18:52:01

gforney commented 9 years ago 
I set TWFIN to 1 and reran.  At around 0.27 s, at a short distance above the burner,
helium is being created, but debug is showing that this not occurring within the
combustion routine.  The combustion routine is being passed a YY with He already
present. Still looking.

Original issue reported on code.google.com by drjfloyd on 2009-09-17 19:15:20

gforney commented 9 years ago 
If it's any help.  I cut the resolution in half and still have the problem.  Also, I
set the DYDX to zero for I==1 in divg Line 180 and this also solves the problem.  It
is weird.

Original issue reported on code.google.com by randy.mcdermott on 2009-09-17 19:50:39

gforney commented 9 years ago 
I've noticed the following.  For the initial period Helium above the burner looks like:
x=0 x=1 x=2 x=3 x=4
0.0000000E+00  0.5347502E-23  0.2357265E-22  0.1751564E-22  0.1004066E-21
0.0000000E+00  0.2503321E-21  0.9426514E-22  0.7191823E-21  0.1198281E-21
0.0000000E+00  0.3066339E-21  0.4900950E-20  0.0000000E+00  0.1095634E-19

If you set the slice file to very low values this appears to just be numerical
diffusion of the much larger concentrations at the INLET.

Just before things go off the rails Helium starts to look like this:
x=0 x=1 x=2 x=3 x=4
0.0000000E+00  0.1296900E-10  0.0000000E+00  0.1174806E-10  0.0000000E+00
0.0000000E+00  0.0000000E+00  0.4141022E-10  0.0000000E+00  0.2183409E-10
0.0000000E+00  0.7702138E-10  0.0000000E+00  0.5458276E-10  0.0000000E+00
0.0000000E+00  0.0000000E+00  0.8165819E-10  0.0000000E+00  0.4037966E-10

A checkerboard pattern of positive and zero valued helium. I think the problem lies
in whatever is causing this.

Original issue reported on code.google.com by drjfloyd on 2009-09-17 20:54:42

gforney commented 9 years ago 
Randy, Jason -- does this case work with the release version, 5.4.1?

Original issue reported on code.google.com by mcgratta on 2009-09-17 21:28:26

gforney commented 9 years ago 
As far as I can tell, it appears that the release version with default options works
and anything other than flux_limiter=-1 results in a failure.

Original issue reported on code.google.com by drjfloyd on 2009-09-17 21:51:15

gforney commented 9 years ago 
True for this exact case (first one that Jukka posted).  And this would point to FL
as the culprit.  BUT, as I said above, with DNS=F, FL=2 works fine.  Further, with
DNS=T, FL=2, and the call the COMBUSTION commented out, all works fine.  Also, take
out diffusion in the transport scheme and all works fine (hence it is not the fault
of numerical diffusion or dispersion because this is still present if physical
diffusion is not).

Original issue reported on code.google.com by randy.mcdermott on 2009-09-18 12:13:40

gforney commented 9 years ago 
Try running the attached for a couple of seconds. No fire but a lot of activity on 
the He front.

Original issue reported on code.google.com by jukka.vaari@vtt.fi on 2009-09-18 13:08:41

gforney commented 9 years ago 
Jukka,

I don't understand why, but at half the resolution this case looks fine to me.  At
the resolution you posted I get a belch of He from the bottom, but no problems at the
top of the burner like before.

R

Original issue reported on code.google.com by randy.mcdermott on 2009-09-18 13:38:34

gforney commented 9 years ago 
Yes, I've had problems reproducing the problems with 2mm grid. This is why I keep 
posting these 1 mm cases. I'll try to do more testing on 2 mm.

Original issue reported on code.google.com by jukka.vaari@vtt.fi on 2009-09-18 13:42:55

gforney commented 9 years ago 
Withe fine grid, I still get that DNS=F works fine.  We are back to the diffusion
term.  ''In theory'', the FLUX_LIMITER scheme does not care if we are DNS or LES.

Original issue reported on code.google.com by randy.mcdermott on 2009-09-18 13:47:20

gforney commented 9 years ago 
I put some debug into the combustion routine to monitor He in the region where it is
being created.  I didn't notice that the He values going into the combustion routine
being changed by the routine, so I don't think the problem lies in combustion but
rather something to do with Q or species changes in combustion is leading to the
problem somewhere else.

Original issue reported on code.google.com by drjfloyd on 2009-09-18 13:50:27

gforney commented 9 years ago 
Dollars to donuts, this has something to do with DIVERGENCE_PART_1 (where the
diffusion term is initialized and FX, etc. are computed) is called AFTER SCALARF
(which needs FX to have the current diffusion term) but WITHIN the CHANGE_TIME_STEP
loop.  This is why the resolution matters and DNS matters, because this changes how
many times we go through this loop.

I swear I have had this same problem in the past...

Original issue reported on code.google.com by randy.mcdermott on 2009-09-18 14:07:22

gforney commented 9 years ago 
You may be right. Try the attached test case and see how the total mass of He 
evolves inside the box (the devc file). If you don't specify DT, the mass will 
increase over 10-fold and do so rapidly in the beginning when diffusion is fast and

time step changes a lot. But if you set DT small enough, He mass is conserved.

Original issue reported on code.google.com by jukka.vaari@vtt.fi on 2009-09-18 14:15:20

gforney commented 9 years ago 
Yep.  This is it.  When I run with 

&TIME T_END=2., DT=0.001, LOCK_TIME_STEP=.TRUE. /

with your original 1 mm case, everything looks ok.

Original issue reported on code.google.com by randy.mcdermott on 2009-09-18 14:19:59

gforney commented 9 years ago 
Well, things look DIFFERENT anyway.

Original issue reported on code.google.com by randy.mcdermott on 2009-09-18 14:23:34

gforney commented 9 years ago 
Jukka,

In your diffusion test, you can also use

&DUMP MASS_FILE=.TRUE./

Interestingly, the total mass is conserved with the larger time step.

R

Original issue reported on code.google.com by randy.mcdermott on 2009-09-18 14:34:29

gforney commented 9 years ago 
Indeed. I got 'WARNING! RHO_AVG <=0 in BAROCLINIC_CORRECTION' followed by SIGSEGV.

Original issue reported on code.google.com by jukka.vaari@vtt.fi on 2009-09-18 14:34:46

gforney commented 9 years ago 
Ha! I think I have it.

Okay, Kevin is going to get sick of hearing me complain about this, but I am starting
to believe that this whole thing is yet another old fashioned Von Neumann stability
problem.  It has all the symptoms: DNS, DT, DX...  I feel like House (American T.V.
show)...

Notice that in the stability check after velo VN is computed using the SC and MU. 
But with the new GET_DIFFUSIVITY routines for DNS, SC*MU/RHO and D are not
necessarily the same, especially for He (I think SC=0.2 or so, yet FDS default is
SC=0.5).  One of the symptoms of this problem, both with your cup burner and the
diffusion problem is that by starting the case out with a sufficiently small DT I can
get things to work.  Alternatively, I had to set SC=0.1 (no need to specify DT in
this case) and things look fine.  I know, SC should be 0.2 you say.  Well,
alternatively, set SC=0.2 and reduce VN_MAX to 0.9.  The stability limit for an
explicit diffusion scheme is unforgiving and VN_MAX=1 is on the hairy edge.  I would
argue that 0.9 should be our new default actually.  With DYNSMAG, etc. we need to
treat the "laminar" regions of the flow correctly, or we are asking for trouble. 
But, I digress... 

So, try setting SC=0.2 and VN_MAX=0.9 and let me know if things get better for you.

p.s. It looks like Jason was right (numerical problems) and I owe everyone donuts...
Still don't know why things work with FL=-1.

Original issue reported on code.google.com by randy.mcdermott on 2009-09-18 16:19:57

gforney commented 9 years ago 
Is part of the solution that in DNS we need to do what we do in divg and rather than
using Sc or 1/Pr, compute it based on the species?

Original issue reported on code.google.com by drjfloyd on 2009-09-18 16:53:46

gforney commented 9 years ago 
I think that would be the "correct" way to do it. The stability check needs to be
consistent with what actually gets used in the transport equations. One thing you
could do is set SC = MIN(RHO*D/MU) within divg for DNS, but I still think we need to
give ourselves a little safety factor on VN_MAX.

Original issue reported on code.google.com by randy.mcdermott on 2009-09-18 17:43:39

gforney commented 9 years ago 
I mean SC = MIN(MU/(RHO*D))...

Original issue reported on code.google.com by randy.mcdermott on 2009-09-18 17:53:47

gforney commented 9 years ago 
Since the problem is a diffusive term and if we want to use what is in divg, should
the VN portion of check_stability just be moved into divg as a subroutine?  Then we
can use the RHO_D and MU arrays computed in divg and just store VN as global variable
so that CHECK_STABILITY can still access it if needed.

Original issue reported on code.google.com by drjfloyd on 2009-09-18 18:24:55

gforney commented 9 years ago 
Ran the no_fire case with SC=0.1 and VN_MAX=0.9.  All looked well up until the point
where the first eddy of helium reaches the axis of symmetry.  Two images attached
just before and just after this point in time.

Original issue reported on code.google.com by drjfloyd on 2009-09-18 19:14:51

gforney commented 9 years ago 
Jason,

Let me guess, you are running the release version of 5.4.1? In that version,
PROJECTION=T does not work with CYLINDRICAL.  Set PROJECTION=F and see if things work
better.  Or use latest SVN.

R

Original issue reported on code.google.com by randy.mcdermott on 2009-09-18 20:10:40

gforney commented 9 years ago 
This thing really does not want to die.  I ran out to 3 sec with Sc=.2 and VN_MAX=.9
and got an instability.  But with Sc=.1 and VN_MAX=.9, it worked.

All this said, I made a lot of clean up changes to the FLUX_LIMITER scheme as I was
debugging the case, so I no longer exactly have a repository version.  I want to run
some other tests before I commit my changes.

Original issue reported on code.google.com by randy.mcdermott on 2009-09-18 20:42:47

gforney commented 9 years ago 
Jukka,

Please have a look at Issue 883.  I made some mods to FDS6 that I would like you to
verify for your cup burner case.  The case I have on my pc (attached) ran fine.

Thanks.
Randy

Original issue reported on code.google.com by randy.mcdermott on 2009-10-27 21:09:34

gforney commented 9 years ago 
Hi,

I tested SVN5002. With 2mm grid the case ran clean. Both of by test runs with 1 mm

stopped with a numerical instability around 1.5 seconds real time. There appeared to

be nothing wrong in the slices. However this is the time when helium starts reaching

the flame in appreciable amounts.

Jukka

Original issue reported on code.google.com by jukka.vaari@vtt.fi on 2009-10-28 16:10:32

gforney commented 9 years ago 
Thanks Jukka.  Can you please send me your exact input file?

Original issue reported on code.google.com by randy.mcdermott on 2009-10-28 16:15:19

gforney commented 9 years ago 
Here it comes.

Original issue reported on code.google.com by jukka.vaari@vtt.fi on 2009-10-28 16:18:43

gforney commented 9 years ago 
Jukka,

Try with SC=0.2 on MISC.  My case is up to 4.6 sec and still going.

R

Original issue reported on code.google.com by randy.mcdermott on 2009-10-28 18:37:02

gforney commented 9 years ago 
I got a clean run up to 10 s. This did not represent the extinguishing concentration

though, and I'd like to verify that as well, hopefully I can report back tomorrow.
I 
did notice that the HRRPUA slice looked a bit odd with SC=0.2 close to the cup rim

which is the all-important region.

Original issue reported on code.google.com by jukka.vaari@vtt.fi on 2009-10-29 16:02:10

gforney commented 9 years ago 
Should this case be run with DNS=T?  Note that with LES and only setting SC=.2 you
still have PR=.5 and this is likely inconsistent.

My feeling is that we have fixed the stability problem, or at least know what causes
it.  If we guarantee stability with FDS6 options, like setting VN_MAX=.5, then I get
complaints about how long thing take to run.  So, my current feeling is that it is
better to use a less restrictive default which works most of the time and then set
specific options as needed for special cases.

Again, note that with DNS=T you still need to set either SC or PR=.2 for the
stability check until we make that consistent.  You may also have to play with VN_MAX
and VN_MIN.

Original issue reported on code.google.com by randy.mcdermott on 2009-10-29 16:14:46

gforney commented 9 years ago 
Whoops, right. I didn't notice since my original post used DNS=T. So let me run a 
couple of tests with DNS=T and report back, OK? We need to get this issue closed.

Original issue reported on code.google.com by jukka.vaari@vtt.fi on 2009-10-29 16:25:44

gforney commented 9 years ago 
Randy,

Are there specific conditions that appear to lead to the instability with the default
options.  Is it because of the large difference in molecular weight or viscosity? 
If
there is any general rules of thumb as to when one might need to change VN, SC, or
PR, we should add them to the User's Guide when it is updated for FDS 6.

Original issue reported on code.google.com by drjfloyd on 2009-10-29 17:46:12

gforney commented 9 years ago 
As we discussed above, the key issue is the consistency between the diffusivity used
in the transport scheme and the one used in the stability check.  In the long run I
hope the user will not need any special modification.  For now though when D is
computed for helium or hydrogen, for example, then we use SC*D_AIR in the stability
check we have problems.  At the moment I am just taking a stab that .2 is a good SC
for helium.

Original issue reported on code.google.com by randy.mcdermott on 2009-10-29 18:09:03

gforney commented 9 years ago 
Randy,

With DNS=T I was able to run my test cases up to extinguishment and beyond. The He

extinguishing concentration has increased since SVN4792, but besides the default SC,

I used CFL_VELOCITY_NORM=1 and FLUX_LIMITER=-1 back then. 

This issue was about stability, and basically I think you've got it clinched, so I

suggest we call this one 'verified'.

Jukka

Original issue reported on code.google.com by jukka.vaari@vtt.fi on 2009-10-30 11:06:33

gforney commented 9 years ago 
Thanks Jukka!

Original issue reported on code.google.com by randy.mcdermott on 2009-10-30 12:25:17