Numerical instability - HVAC Louvered Vent

[gforney]

FDS Version: 6.0.1
SVN Revision Number: 17534
Compile Date:Tue, 26 Nov 2013
Operating System: OSX 10.9.3

Hello,

I am trying to model a naturally-ventilated, underground, train station. Vertical shafts
are used to connect the station to the surface. 
The top of the shafts are capped and have louvers on their sides.
I am modelling the inside of the shafts as well as the outside space in order to perform
a sensitivity analysis on the effects of wind.
I have used HVAC in order to model the losses and change of direction caused by the
louvers.

I have attached a simplified version of my input file, with a single shaft.
As dr_jfloyd discovered, the problem seems to be solved once the &PRES line is removed,
which I cannot explain.

Kind Regards,
Evangelos 

Original issue reported on code.google.com by evantzanidis on 2014-06-05 13:11:14

[gforney]

There is no attachment.

Original issue reported on code.google.com by mcgratta on 2014-06-05 13:31:39

[gforney]

Oh sorry,
Here it is

Original issue reported on code.google.com by evantzanidis on 2014-06-05 13:34:16

---

• _Attachment: VENT_CAP.fds_

[gforney]

Attached is a simplified version of the case from the group.
The original case had 4 louvers and two meshes, this case has one mesh and two louvers

It dies with numerical instability if the PRES line is present.  If I remove the second
louver it will run.  If I remove the PRES line it will run.  

Also, the first louver shows upward flow when UVW is negative for W.  The second louver
shows downward flow (as did the louvers on the other two faces removed for simplifying).

Original issue reported on code.google.com by drjfloyd on 2014-06-05 13:35:52

---

• _Attachment: VENT_CAP.fds_

[gforney]

(No text was entered with this change)

Original issue reported on code.google.com by drjfloyd on 2014-06-05 13:41:08

[gforney]

The attached file is a small shaft with a vent at the bottom introducing mass.  Two
louvers are on the side of  the shaft venting to ambient.  I have turned off noise.
 The case should remain symmetric.  If I comment out the PRES line, the case runs fine
with slight asymmetries (E-6 differences in velocity) developing.  If I leave the pres
line in, it dies.  Large asymmetries develop in the velocity field with pressure iterations
turned on and those asymmetries play havoc with the HVAC model.  Adding CHECK_VN does
not solve the problem.  Any thoughts?

Any thoughts on what to check?

Original issue reported on code.google.com by drjfloyd on 2014-06-09 15:29:16

• Labels added: Application-FDS

---

• _Attachment: VENT_CAP.fds_

[gforney]

I'll have a look.

Original issue reported on code.google.com by randy.mcdermott on 2014-06-09 15:45:59

[gforney]

Jason,

Can this case be done without HVAC?  That is, can you specify VENTS with MASS_FLOW
or whatever and have the same problem?  Just wondering if we can take the HVAC code
out of the equation for now.

R

Original issue reported on code.google.com by randy.mcdermott on 2014-06-09 15:47:57

[gforney]

I can mock up an equivalent case with MASS_FLOW, but my guess is that will work as any
asymmetric in the velocity field won't change the vent flow.

Original issue reported on code.google.com by drjfloyd on 2014-06-09 15:49:25

[gforney]

It runs fine with mass flux specified rather than using HVAC.

Original issue reported on code.google.com by drjfloyd on 2014-06-09 15:55:47

---

• _Attachment: VENT_CAP_MF.fds_

[gforney]

I've tried older versions, no angle on vent, smoother vel ramp.  No luck.

One thing to note is that, despite the PRES line, you never actually get a solution
to the pressure equation that satisfies the velocity tolerance.  So that's obviously
not good.  I need to understand your UVW bc a bit better in order to help.

Original issue reported on code.google.com by randy.mcdermott on 2014-06-09 16:43:27

[gforney]

UVW isn't the issue.  Get rid of it and the case will still die.

UVW is just a way of specifying VEL_T without having to go through the contortions
of figuring out which velocity component to give to VEL_T(1) and which to VEL_T(2).
 UVW is just the outward point normal from the HVAC vent.  Other than that the code
logic is basically the same as for VEL_T.

Original issue reported on code.google.com by drjfloyd on 2014-06-09 16:47:33

[gforney]

OK.  Then point me to the code where you specify the normal component of velocity for
a wall cell.

Original issue reported on code.google.com by randy.mcdermott on 2014-06-09 16:49:47

[gforney]

HVAC_BC in wall.f90; it is implemented as a specified mass flux boundary.

Original issue reported on code.google.com by drjfloyd on 2014-06-09 16:55:39

[gforney]

This may not be related to the pressure iterations.  Made a modified case that dies
without iterations.

Original issue reported on code.google.com by drjfloyd on 2014-06-09 17:33:34

[gforney]

I tried to make a case with length = 0 and loss = 0,0 and it would not allow flow through
the ducts.  Wouldn't you just have in=out with no losses?  HVAC does not seem to like
this situation.

Original issue reported on code.google.com by randy.mcdermott on 2014-06-09 20:01:13

[gforney]

The concept of s momentum length is integral to the assumptions in the
model.  A zero length duct is not a duct.

Original issue reported on code.google.com by drjfloyd on 2014-06-09 20:10:57

[gforney]

If I set length=1, all behaves fine.  See attached case.

Looking at your solution procedure in the Tech Guide, it is hard for me to tell whether
this scheme is fully implicit or not.  Meaning, do you have a time step restriction
for stability of your node solve?  If you do, then it is related to L_j, which is in
the denominator.  I'm sort of grasping here, because I am not that familiar with your
scheme.  But I'd say it would be good if you could handle the case with L=0 and K=0.
 This would mean you are pretty robust.

Original issue reported on code.google.com by randy.mcdermott on 2014-06-09 20:19:19

---

• Attachment: save1.fds
• Attachment: Screen Shot 2014-06-09 at 4.16.02 PM.png
  

[gforney]

I may have misspoke about the losses.  I guess you always have at least a K=1.

But what I am saying is: how should your scheme handle a duct of L=0.001 vs. L=1000.
 They are treated the same?

Original issue reported on code.google.com by randy.mcdermott on 2014-06-09 20:25:19

[gforney]

Scheme is iterative semi implicit.  Velocity is a matrix per the equation
in the tech guide,  mass and energy update explicitly with new velocity.

Length affects acceleration.  Less length more acceleration.  Hmm... Maybe
there needs to be some thought on this if length of 1 works.  The &PRES was
just a red herring I guess.

Original issue reported on code.google.com by drjfloyd on 2014-06-09 20:52:12

[gforney]

If you take the matrix you generate and plug it into Matlab, it will tell you the eigenvalues.
 Do eig(A).  If any of your eigenvalues are positive, you will have stability problems.
 You could also look at the condition number, cond(A).  Large condition numbers are
bad for stability.  I suspect your condition number depends on L_j.

To me, Eq. (10.4) looks explicit with a CFL given by dp*dt/(rho*L).  If this were a
purely explicit scheme, you would need to make sure this number is less than 0.5. 
Do you subcycle in time or just iterate?  If L is small, you might need to take small
subtimesteps.

Original issue reported on code.google.com by randy.mcdermott on 2014-06-09 21:04:38

[gforney]

The underlying problem appears to have been the fact that the HVAC and FDS pressure
solutions are indirectly coupled.  Numerical error resulted in one face of the shaft
having a slightly different velocity at the HVAC vent than its opposite.  This fed
back into the HVAC solution which fed back into the FDS pressure solution.  The result
was an undamped feedback between them.  I added some relaxation to the HVAC boundary
condition for the perturbation pressure. 

I tested this with a modified case where I changed one of the OPEN boundaries to be
a inflow (represent wind).  Looking at the duct mass flows on opposite faces the two
are equal until the wind condition reaches the shaft at which point the upstream face
 has decreased (and ultimately reverses) flow.  Also I upped your loss to 5 which is
more representative of a louvered vent than 1 (which is for a fully open vent).

Original issue reported on code.google.com by drjfloyd on 2014-06-17 15:05:15

---

• _Attachment: VENT_CAP.fds_

[gforney]

changing to fixed

Original issue reported on code.google.com by drjfloyd on 2014-06-17 15:21:55

[gforney]

I tested some cases that include HVAC, on FDS6.1.1. I noticed that the relaxation factor
makes a significant improvement but it still becomes unstable if you reduce the length
of the duct enough. In the last modified case that you uploaded, a numerical instability
occurs for a duct length of 0.05 meters. 

Furthermore, increasing the size of the vent, also increases the minimum duct length
that allows for a stable solution. For exampple, in a case of mine (attached - simplified)
I have placed 20m^2 vents that connect my domain to ambient conditions through HVAC
and represent louvers. The minimum  duct length that does not cause an instability
is 2 meters, which is quite large, considering that the thickness of the louvers should
be around 0.1 meters. For lower lengths it does not result in numerical instability,
but rather it continues with high velocity errors and constantly at the maximum number
of pressure iterations.

As you said, length is inversely proportional to acceleration, thus small lengths cause
large accelerations. I believe that large accelerations, along with a relatively large
time steps, cause the instability.

p.s. In my case I use multiple vents per shaft, in order to allow for different velocities
(i.e. air may be supplied from one side of the shaft, while smoke is being exhausted
from the other side. The problem persists, even if you keep a single vent and delete
the rest.

Original issue reported on code.google.com by evantzanidis on 2014-07-15 15:16:06

---

• _Attachment: base_station.fds_

[gforney]

How sensitive are your results to the duct length that is used? If you are just looking
to get some form of time-averaged of impact of wind on the flow out the louvers, then
maybe duct length won't have much impact other than increasing the time before you
reach steady-state.  

Original issue reported on code.google.com by drjfloyd on 2014-07-15 22:33:54

[gforney]

Its a parametric study investigating the sensitivity of natural smoke ventilation to
both design parameters and external conditions. I was hoping to compare the times that
it takes for the station to become untenable but that seems difficult since increasing
the length of the ducts introduces a damping to acceleration through the HVAC system.
I might run a few simulations to check the impact of duct length to the rate of smoke
filling.  

Original issue reported on code.google.com by evantzanidis on 2014-07-16 09:42:18

[gforney]

I'm looking at your geometry in Smokeview.  The LENGTH parameter is used to determine
inertia to acceleration within a duct.   In your case, the flow out the louver is being
driven by the flow in the shaft.  In some sense the inertia you need to overcome to
change flow at the louver is something on the order of the shaft width by the louvers
which is around 10 m.  In your case, having to use a duct length around 2 m for stability
is not a bad physical picture of the flow conditions.  I think about what I might do
to improve stability.  I'm reluctant to just increase the relaxation as we do want
to capture short, real, pressure transients and too much relaxation could smooth those
out.

Original issue reported on code.google.com by drjfloyd on 2014-07-17 14:52:17

[gforney]

I wonder if an alternate meshing strategy would help you.  Your case is maxing out pressure
iterations.  Your shafts are overlapping multiple meshes and the way you have divided
up around the fire is giving multiple mesh overlaps at the ends of the center region.

The attached image shows your current mesh layout in the top half.  The bottom half
is an alternate layout that should make the mesh to mesh coupling easier on FDS.  It
might help.

Also looking at the Smokeview velocity vectors near the fire, I think your fire mesh
might benefit from being a little larger.   

Original issue reported on code.google.com by drjfloyd on 2014-07-17 17:20:48

---

• Attachment: Clipboard01.png
  

[gforney]

I agree that overdoing it with the relaxation factor will introduce a damping, thus
making changes in velocity sluggish. The HVAC calculation procedure seems very robust,
I guess it is the LES solver that struggles to cope with the BCs that are imposed to
it.

I applied the changes that you suggested and the pressure iterations have reduced (I
havent completed the simulation yet). Thank you for the feedback!
Do you think that placing multiple vents one after the other along the same wall could
cause an instability?

Original issue reported on code.google.com by evantzanidis on 2014-07-17 20:16:45

[gforney]

I run the case with the meshing strategy that you suggested and the problems persist.
I also tried assigning BAROCLINIC=.FALSE., which only improved the situation slighly.
Furthermore, I tried increasing the length of the ducts, first to 4 and then to 8 meters,
which did not help either.

The behaviour is always the same, the pressure iterations max out after about 400 seconds,
at which time the pressure (vertical slicefile) starts varying excesively in space
and time.

I have attached my last input file.

Original issue reported on code.google.com by evantzanidis on 2014-07-18 11:37:39

---

• Attachment: VCoarse.fds

[gforney]

I tried increasing the mesh resolution. The refined mesh has double the number of cells
in each direction. This caused the pressure iterations to reduce but they still max
out at some time-steps. Furthermore, the pressure continues to vary erratically.

I have attached the input file and 2 screenshots of pressure (474 ans 475.2 seconds)

Original issue reported on code.google.com by evantzanidis on 2014-07-18 17:38:31

---

• Attachment: pressure1.png
  
• Attachment: pressure2.png
  
• Attachment: Coarse.fds

[gforney]

Thought of a possible alternative approach.   Rather than using louvered HVAC vents
you could try using the screen drag model.  You won't get the louvered effect for the
outflow but you should be able to do something reasonable for the pressure drop.  
What you would need to do is estimate what the flow will be out your louvers during
the fire.  The using the HVAC momentum equation you can do a hand calc to determine
the pressure drop you expect.  Then take the screen_drag example case in Examples\Sprinkler_and_Sprays
and modify it to have the same grid size and dimensions as one of your louvers.  Use
the velocity you estimated and play with the screen drag model inputs to get the correct
pressure drop.  This won't be perfect, but it should hopefully avoid the issue you
are currently having with the HVAC routines

Original issue reported on code.google.com by drjfloyd on 2014-07-24 12:03:57

[gforney]

This is a very good idea. Do you know what kind of boundary this translates to? I couldnt
understand from the tech guide.

I found something that might be useful. The pressure iterations never max-out if I
replace all HVAC-Vents with open vents. The domain might be split up into too many
meshes but it performs decently when these boundaries are removed

Original issue reported on code.google.com by evantzanidis on 2014-07-24 13:21:48

[gforney]

An open vent isn't going to account for the pressure loss through the louver.  You will
overestimate the flows by probably quite a bit.

The screen drag model is a non-isotropic porous media model where drag is imposed on
the flow only in the direction normal to the screen.  You have a louver and not a wire
mesh screen, but it is similar in the sense that you have a plane where a large pressure
drop is imposed on flow normal to the plane.

The reason I am suggestion the screen drag is that it is fully coupled with the FDS
pressure solver, so you should avoid the bad feedback problems that you are currently
getting trying to do this with the HVAC solver.  There may be more to be done with
HVAC development in the longer term to get this working with HVAC.

Original issue reported on code.google.com by drjfloyd on 2014-07-24 13:31:50

[gforney]

I understand, I will try it out

Original issue reported on code.google.com by evantzanidis on 2014-07-24 13:34:41

[gforney]

ruchinkumar99,

The issue tracker is for reporting bugs.  Questions on using FDS should be submitted
to the discussion forum (https://groups.google.com/d/forum/fds-smv)

Also, when using the issue tracker please do not add a second topic to an existing
issue.  

Original issue reported on code.google.com by drjfloyd on 2015-04-03 11:53:43