Porous Media - zeta

[Er9y714]

Following the porous media section in User Guide 15.4.8, I made a table. Pressure drop for different velocities are calculated by following Eq. 15.6

delta: 0.2 m
viscosity of air: 1.8 E-5 kg/(ms)
density: 1.2 kg/m3
K: 1.5098E-08 m2
Y: 0.1

velocity	Pressure
0.10	26
0.20	56
0.30	89
0.40	127
0.50	168
0.60	213
0.70	263
0.80	316
0.90	373
1.00	434
1.10	499
1.20	567
1.50	797
2.00	1,258
3.00	2,473
4.00	4,079
5.00	6,075
6.00	8,462
7.00	11,240
8.00	14,408
9.00	17,967
10.00	21,916

Then to find the resulting loss coefficient in each calculation, I derived the loss coefficient, zeta, by:

zeta =

I populated the same table with zeta:

velocity	Pressure	zeta
0.10	26	4300
0.20	56	2313
0.30	89	1650
0.40	127	1319
0.50	168	1120
0.60	213	988
0.70	263	893
0.80	316	822
0.90	373	767
1.00	434	723
1.10	499	687
1.20	567	657
1.50	797	590
2.00	1,258	524
3.00	2,473	458
4.00	4,079	425
5.00	6,075	405
6.00	8,462	392
7.00	11,240	382
8.00	14,408	375
9.00	17,967	370
10.00	21,916	365

According to the plot, zeta converges after 4 m/s. It is very high at small velocities. Is there a way to utilize this method while keeping the zeta constant for all velocities? In that way it can be used as a ventilation element with a loss coefficient.

Thanks,

[rmcdermo]

Sorry, but I do not understand what you are trying to do. If your zeta equation worked, you would not need the linear term in Eq. 15.6. What is the point of trying to equate these?

[Er9y714]

I could not find sphere law in user guide.
Sorry for going in more detail, but I think I need to explain what I try more in detail:

In HVAC, we can use ducts and give a loss coefficient. For any velocity we can simply calculate the pressure drop by dP=0.5x
rho x zeta x v^2.
I'd like to model the same behaviour with stationary particles.

Steps:
For example for a louvre with zeta of 50, dP at 2 m/s is: dP = 0.5x1.2x50x2^2 = 120 Pa
To find a K value which will give the same pressure drop in eq.15.6:

dP=120 Pa
delta: 0.2 m
viscosity of air: 1.8 E-5 kg/(ms)
density: 1.2 kg/m3
Y: 0.1
v=2 m/s

I get K=7.55E-07

Here begins my issue;
if 0.25 m/s is applied, I would expect a dP=0.5x1.2x50x0.25^2=1.87Pa
however, Eq. 15.6 gives 2.92Pa.
So in that method of particles, I can't get the identical pressure drop compared to HVAC.

[rmcdermo]

I think you just want to set DRAG_COEFFICIENT=0.5 for the particle class. Please try and see if that sorts you out.

[rmcdermo]

Sorry, I think the 0.5 you had in your expression above is just part of your formula, so you just need DRAG_COEFFICIENT=1.

[Er9y714]

Thank you for pointing out Drag coefficient to me. However it does not achieve what I wanted to do with Eq.15.6.

However, I have found a possible solution in verification guide, 13.1.4 Stationary Spherical Particles in a Duct (sphere_drag).

On hand calculations;
Once I tune the Cd, for varying velocity the zeta remains the same and pressure drop is matches with simple eq (dP=0.5 rho zeta u^2).

No.particle	4000			No.particle	4000
Cd	159.5			Cd	159.5
pi	3.1415			pi	3.1415
r	0.005			r	0.005
u	2			u	0.25
A	1			A	1
rho	1.2			rho	1.2
dp	120.26	Pa		dp	1.88	Pa

I will continue investigating this application. Thank you for your help!

[Er9y714]

I tried higher drag coefficients in this example: sphere_drag_1.fds

However, I do not see a higher pressure drop for DRAG_COEFFICIENT higher than 20 in that example. It seems to hit a limit.
Therefore, for 2 m/s velocity a maximum of 15 Pa is achieved. This corresponds to zeta = 6.25.

Did I reach an implicit limit in this method?

&HEAD CHID='sphere_drag_1', TITLE='Planar array of spherical particles in a duct' /

&MESH IJK=40,20,20,XB=0,2,0,1,0,1/
&MESH IJK=40,20,20,XB=0,2,0,1,1.1,2.1/
&MESH IJK=40,20,20,XB=0,2,0,1,2.2,3.2/
&MESH IJK=40,20,20,XB=0,2,0,1,3.3,4.3/
&MESH IJK=40,20,20,XB=0,2,0,1,4.4,5.4/

&TIME T_END=10. /
&RADI RADIATION=.FALSE./
&MISC PARTICLE_CFL_MAX=0.9/
&VENT MB='XMIN',SURF_ID='VELIN'/
&VENT MB='XMAX',SURF_ID='OPEN'/
&VENT MB='YMIN',SURF_ID='SLIP'/
&VENT MB='YMAX',SURF_ID='SLIP'/
&VENT MB='ZMIN',SURF_ID='SLIP'/
&VENT MB='ZMAX',SURF_ID='SLIP'/

&SURF ID='VELIN',VEL=-2/
&SURF ID='SLIP',FREE_SLIP=.TRUE. /

&SLCF PBY=0.5,QUANTITY='VELOCITY',VECTOR=.TRUE./
&SLCF PBY=0.5,QUANTITY='PRESSURE'/

&SURF ID='PARTICLE',THICKNESS=0.005, GEOMETRY='SPHERICAL',MATL_ID='PARTICLE' /
&MATL ID='PARTICLE',DENSITY=1000.,CONDUCTIVITY=1,SPECIFIC_HEAT=1. /
&PART ID='DRAG 5',  SURF_ID='PARTICLE', DRAG_COEFFICIENT=5,  STATIC=.TRUE.  /
&PART ID='DRAG 10', SURF_ID='PARTICLE', DRAG_COEFFICIENT=10, STATIC=.TRUE.  /
&PART ID='DRAG 20', SURF_ID='PARTICLE', DRAG_COEFFICIENT=20, STATIC=.TRUE.  /
&PART ID='DRAG 50', SURF_ID='PARTICLE', DRAG_COEFFICIENT=50, STATIC=.TRUE.  /
&PART ID='DRAG 100', SURF_ID='PARTICLE', DRAG_COEFFICIENT=100, STATIC=.TRUE.  /

&INIT N_PARTICLES_PER_CELL=10, CELL_CENTERED=.TRUE., PART_ID='DRAG 5',  XB=1.01,1.02,0.0,1.0,0.0,1.0 /
&INIT N_PARTICLES_PER_CELL=10, CELL_CENTERED=.TRUE., PART_ID='DRAG 10', XB=1.01,1.02,0.0,1.0,1.1,2.1/
&INIT N_PARTICLES_PER_CELL=10, CELL_CENTERED=.TRUE., PART_ID='DRAG 20', XB=1.01,1.02,0.0,1.0,2.2,3.2/
&INIT N_PARTICLES_PER_CELL=10, CELL_CENTERED=.TRUE., PART_ID='DRAG 50', XB=1.01,1.02,0.0,1.0,3.3,4.3/
&INIT N_PARTICLES_PER_CELL=10, CELL_CENTERED=.TRUE., PART_ID='DRAG 100', XB=1.01,1.02,0.0,1.0,4.4,5.4/

&DUMP DT_DEVC=0.1 /

&DEVC XYZ=0.50,0.5,0.50, QUANTITY='PRESSURE', ID='FDS 5'  /
&DEVC XYZ=0.50,0.5,1.55, QUANTITY='PRESSURE', ID='FDS 10' /
&DEVC XYZ=0.50,0.5,2.60, QUANTITY='PRESSURE', ID='FDS 20' /
&DEVC XYZ=0.50,0.5,3.70, QUANTITY='PRESSURE', ID='FDS 50' /
&DEVC XYZ=0.50,0.5,4.80, QUANTITY='PRESSURE', ID='FDS 100' /

&TAIL /

[rmcdermo]

This is a time step problem. I changed this line in your input file and get the slice results below. I'll see if we can automate this.

&TIME T_END=10., DT=0.001 /

[Er9y714]

Great, looking forward to it!

I would expect PARTICLE_CFL_MAX=0.9 would prevent this issue, but did not.

If I set DT=0.001, FDS might increase the timestep if needed in longer runs. I suppose, this might not be sufficient as a work-around.

[rmcdermo]

The latest test bundle should work. But note I had to use PARTICLE_CFL_MAX=0.1 to get accurate results. At first I thought this was due to the velocity ramp, which defaults to TAU_V=1. But after playing with this awhile it seems that indeed you need this tight time step for the high drag case just to maintain the steady state. We will think about this more, perhaps there is something we can do to stabilize this better. But for now, at least PARTICLE_CFL_MAX affects the result. But for extremely high packing ratio, you ware going to need a tight time step to get accurate results.

https://github.com/firemodels/test_bundles/releases/tag/FDS_TEST

[Er9y714]

Thank you! I will test it, towards the end of the week.
Do you refer "extremely high packing ratio" to N_PARTICLES_PER_CELL ? If yes, we can reduce it to 1 and increase the DRAG_COEFFICIENT or THICKNESS to get the same total value.

[Er9y714]

I tested the bundle.

• With PARTICLE_CFL_MAX=0.1, the higher limit of N_PARTICLES_PER_CELL is around 10. The higher the N_PARTICLES_PER_CELL, the higher the pressure loss error. (the input below is an example, in which I would expect 18.8 Pa, but measured 16 Pa)
• When I set the N_PARTICLES_PER_CELL to 1, pressure loss is correctly calculated between DRAG_COEFFICIENT = 1 and 10000 without adding any parameter for PARTICLE_CFL_MAX.
  • For static particle applications, it seems there is no advantage of N_PARTICLES_PER_CELL larger than 1.

&HEAD CHID='sphere_drag_1', TITLE='Planar array of spherical particles in a duct' /

&MESH IJK=40,20,20,XB=0,2,0,1,4.4,5.4/ 

&TIME T_END=1. /
&RADI RADIATION=.FALSE./ 
&MISC PARTICLE_CFL_MAX=0.1/

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

&SURF ID='VELIN',VEL=-2 TAU_V=0.01/
&SURF ID='SLIP',FREE_SLIP=.TRUE. /


&SLCF PBY=0.5,QUANTITY='VELOCITY',VECTOR=.TRUE./
&SLCF PBY=0.5,QUANTITY='PRESSURE'/

&SURF ID='PARTICLE',THICKNESS=0.005, GEOMETRY='SPHERICAL',MATL_ID='PARTICLE' /
&MATL ID='PARTICLE',DENSITY=1000.,CONDUCTIVITY=1,SPECIFIC_HEAT=1. / 
&PART ID='DRAG E', SURF_ID='PARTICLE', DRAG_COEFFICIENT=10, STATIC=.TRUE.  / 
 
&INIT N_PARTICLES_PER_CELL=25, CELL_CENTERED=.TRUE., PART_ID='DRAG E', XB=1.01,1.02,0.0,1.0,4.4,5.4/ 

&DUMP DT_DEVC=0.01 /
 
&DEVC XYZ=0.50,0.5,4.80, QUANTITY='PRESSURE', ID='FDS E'/ 

&TAIL /

[rmcdermo]

The reason the pressure gradient hits a limit as the drag force is increased is seen in this block on code in part.f90.

! Add summed particle accelerations to the momentum equation. Limit the value to plus/minus abs(u)/dt to prevent a sudden
! change in gas direction.

DO K=0,KBAR
   DO J=0,JBAR
      DO I=0,IBAR
         UODT = ABS(UU(I,J,K)*RDT)
         VODT = ABS(VV(I,J,K)*RDT)
         WODT = ABS(WW(I,J,K)*RDT)
         FVX(I,J,K) = FVX(I,J,K) + MIN(UODT,MAX(-UODT,FVX_D(I,J,K)))
         FVY(I,J,K) = FVY(I,J,K) + MIN(VODT,MAX(-VODT,FVY_D(I,J,K)))
         FVZ(I,J,K) = FVZ(I,J,K) + MIN(WODT,MAX(-WODT,FVZ_D(I,J,K)))
      ENDDO
   ENDDO
ENDDO

I played with removing this limit and the case immediately goes unstable. This is one of those things that bites us, when we hit these limits silently. I have no idea how often this happens in practice. But it is worrisome because the drag on static particles is the foundation of the "physics-based" approach to modeling wildfires with Lagrangian particles. Here we are seeing pressures off by a factor of more than 2 for high drag forces unless a tight PARTICLE_CFL_MAX is used.

My current thinking is that we need to try to apply a smoothing (or relaxation) factor to the drag for these static particles. We would like to avoid needing infinitesimal time steps to maintain an accurate steady state.

[rmcdermo]

@Er9y714 I added your additional cases to the sphere_drag_1 verification case and set PARTICLE_CFL_MAX=0.1 for this case. I think we have to leave it here for now. Thanks for bringing this up. Before digging into this I did not appreciate that we had this limit on our drag force term. I wonder if this is biting us in WUI cases with vegetation. Now that we at least have PARTICLE_CFL_MAX hooked up to normal stationary particles, we may see an effect on the WUI cases. I've nominated this for a "Boris". Thanks

[Er9y714]

I appreciate it your time and the nomination!

As I posted above, I have tested the bundle. It provides a solution for static particles for my intended application without using the PARTICLE_CFL_MAX. I look forward to seeing the next release to make use of it and maybe a complete solution to the issue.