Strange solid particle behavior at centerline of CYLINDRICAL flow

[rmcdermo]

This is a case I'm working that is looking at using sodium bicarbonate for fire suppression. So, this is a bit new in terms of using solid particles for suppression where we would normally use the droplet model. Nevertheless, there seems to be some behavior at the centerline that does not look right. At about 1.69 s in the simulation, you see the particles on the centerline look like they start falling due to gravity and create a backflow that leads to unstable behavior. The simulations keeps going, but something is not right.

There also seems to be problems whenever the particles get close to the centerline downstream of the burner. In this particular case, it is hard to tell if the explosions are the result of problems upstream. But other variants of this case where the upstream looks ok produce spurious issues downstream, again, whenever the particle gets close to the centerline.

I tried taking CYLINDRICAL out of play, but the OPEN bc does not act correctly for 2D calculations (this is another issue that has never really been solved).

Anyway, at this point, I'm fishing for some suggestions. I think this is a flow we should be able to count on working without too much fuss.

test.fds.txt

[rmcdermo]

Here is an image from a bit later in the sim where things are getting worse.

[mcgratta]

Yes, I see the problem. Working it.

[mcgratta]

I don't think this is a CYLINDRICAL problem per se. Consider this simple 2d case, no reactions, MIRROR boundaries. Very simple. The particles fall back against the flow, but I don't have a sense of whether in this case with the particles so small we should expect this.

&HEAD CHID='test2'/

&MESH XB=0.00,0.060,0.0,0.001,0.06,0.15, IJK=30,1,45 /

&TIME T_END=10.0 /

&SURF ID='inlet',
      VEL=-0.1,
      TAU_V=0.
      PART_ID='KSA'
      PARTICLE_MASS_FLUX = 0.005 /

&VENT SURF_ID='OPEN', DB='ZMAX'/
&VENT SURF_ID='inlet', DB='ZMIN'/
&VENT SURF_ID='MIRROR', DB='XMAX'/
&VENT SURF_ID='MIRROR', DB='XMIN'/

&MATL ID               = 'Sodium Bicarbonate'
      FYI              = 'NaHCO3'
      DENSITY          = 2200.
      CONDUCTIVITY     = 0.20
      SPECIFIC_HEAT    = 1.4 /

&SURF ID='KSA'
      HEAT_TRANSFER_COEFFICIENT=50
      MATL_ID = 'Sodium Bicarbonate'
      THICKNESS=2e-06
      GEOMETRY='SPHERICAL' /

&PART ID='KSA'
      SURF_ID='KSA'
      COLOR='BLUE'
      SAMPLING_FACTOR=1
      DRAG_LAW ='SPHERE'
      MINIMUM_DIAMETER=1
      QUANTITIES='PARTICLE TEMPERATURE','PARTICLE VELOCITY','PARTICLE DIAMETER','PARTICLE MASS','PARTICLE WEIGHTING FACTOR' /

&SLCF PBY=0.0005,QUANTITY='TEMPERATURE', CELL_CENTERED=T /
&SLCF PBY=0.0005,QUANTITY='VELOCITY', VECTOR=.TRUE. /

&TAIL /

[mcgratta]

Also, the Re number for these particles is about 0.02, thus the drag coefficient is 24/Re or about 1200. I don't know if the drag law is appropriate for such low Re, and then there is the fact that the particles are clustered. Each particle in the simulation represents about 1300 actual particles.

[drjfloyd]

24/Re is Stokes flow drag, it should be appropriate.

[mcgratta]

Right, when the sphere is isolated. In our case, we have a cloud of spheres.

[drjfloyd]

With 0.002 m grid spacing the ratio of cell volume to the volume of 1,300 particles is ~200,000. That gives about 50 diameters between particles on average. I don't think we would expect much particle-particle interaction.

[mcgratta]

Right, so another question is the 0.02 kg/m2/s mass flux of particles. This is about 1/6 the mass flux of air that is injected with the particles (1.2 kg/m3 x 0.1 m/s = 0.12 kg/m2/s). I'm wondering if we're violating some assumption about particles taking up no volume or something similar?

[drjfloyd]

The particles are 2000 times more dense than air. 1/6 mass flux will be ~1/12,000 volume flux.

[rmcdermo]

For the test2 case, if I set NPPC=10, the instability is much delayed, but eventually does happen after 5 seconds. With NPPC=20, the instability is pushed out to about 9 s.

I think that with NPPC=1, you easily end up with cells that do not have particles, and this would not happen in reality. The variation in f_b then creates a velocity fluctuation (one cell with a downward force next to a cell without one creates vorticity) that leads to spurious turbulence in the numerical solution.

Something I have been thinking we will need for embers anyway is to add a random walk based on the local fluid eddy viscosity. Particles should experience a subgrid force due to unresolved eddies. This randomizes the position and minimizes the chances of these "channels" forming. This would then minimize the required NPPC, which is a net speed up. The random walk is very cheap.

[mcgratta]

Hmm -- seemed like there were enough particles. Are you setting SAMPLING_FACTOR=1 on the PART line?

[rmcdermo]

I keep SAMPLING_FACTOR=NPPC, so as not to slow SMV.

[mcgratta]

Also, another thing to do is set DT_INSERT to a value less than the default 0.01 s on the SURF line. If set low enough, you get particles every time step. So a mix of NPPC and DT_INSERT should help populate the domain reasonably well.

[rmcdermo]

With DT_INSERT=0.005 and NPPC=10 (SAMPLING_FACTOR=100), things look reasonable.

[mcgratta]

The number of particles per grid cell will be approximately

DZNPPC/(W0DT_INSERT) = 0.002 m x 10 / (0.1 m/s x 0.005 s) = 40

Your original case had 2 particles per cell, if my math is right.

[rmcdermo]

Kevin, We are not out of the woods. Take a look at this case (just unzip and run SMV, look at particles and velocity vectors). It is the original cylindrical test, but now with a few more particles. I have also added an OBST for the pedestal of the cup, so it takes the centerline issue out of play. As the particles get very close to the wall, their upward drag is reduced and then eventually settle back toward the inlet. Once they are falling back through the inlet, things go haywire. Not sure if there is something we need to do in terms of removing these particles before they cause this havoc. Thoughts?

test_save.zip

[mcgratta]

Set SAMPLING_FACTOR=1 and you'll see a distinct problem near the wall, something that appears to be an interpolation thing. I'll dig in moire.

[mcgratta]

I understand what is happening at the wall. When a particle is within one half a grid cell to the vertical wall, it takes the gas phase value of w instead of interpolating the gas and ghost cell value of w. Thus, we don't see particles "sticking" to the wall. That being said, there seems to be a build up of particles within the "boundary layer" which eventually causes trouble. Not sure why.

[rmcdermo]

Things are looking pretty good so far using the UGLMAT pressure solver Marcos just pushed up (I am assuming the Poisson errors might be a diagnostic mistake due to the unstructured grid).

[rmcdermo]

Doh! Spoke too soon.

[mcgratta]

I don't think this has to do with the pressure solver. It's something about the momentum transfer from particles to gas, which is tied in with the ODE solver for particles.

[rmcdermo]

My hunch was that any velocity penetration error could be pushing the particles too close to the boundary, or something along those lines.

[mcgratta]

Add these lines to your case

&SLCF PBY=0.0005,QUANTITY='H', CELL_CENTERED=T /
&SLCF PBY=0.0005,QUANTITY='F_X' /
&SLCF PBY=0.0005,QUANTITY='F_Y' /
&SLCF PBY=0.0005,QUANTITY='F_Z' /

At the boundary, these acceleration terms fluctuate due to the fluctuation of H in the solid. I have to check when the particle acceleration (ACCEL_X, etc) is added. This might be a bad feedback loop.

[mcgratta]

I checked and made sure that FVX and FVZ at solid boundaries do not have contributions from ACCEL_X, ACCEL_Z. Still, things seem awfully noisy near these boundaries.

[mcgratta]

If you load and view alot of particles, it appears that the instability originates at the separation point of a counter-clockwise vortex that forms in the corner below the burner. I suspect that the del dot F term in the Poisson equation gets very large or small due to the clustering of particles there.

[rmcdermo]

Something really seems off with the particle transport. While investigating another issue, I set up this very simple case. It should just blow neutrally buoyant particles (DENSITY=1.2 kg/m3) up and out of the domain. But the case fails almost immediately, just a few time steps, with an instability.

In this case, I would expect everything to be isothermal, but the divergence goes off the rails. Not sure if that is the chicken or the egg.

&HEAD CHID='test'/

&MESH IJK=10,10,10, XB=0,1,0,1,0,1/

&TIME T_END=10/

&MISC PARTICLE_CFL=T, PARTICLE_CFL_MAX=0.5/

&SURF ID='inlet'
      COLOR='RED'
      VEL=-0.1
      PART_ID='spheres'
      PARTICLE_MASS_FLUX=0.02
      NPPC=1
      DT_INSERT=0.01/

&PART ID='spheres', SURF_ID='part surf',
      QUANTITIES='PARTICLE VELOCITY','PARTICLE DIAMETER','PARTICLE WEIGHTING FACTOR'
      SAMPLING_FACTOR=1/
&SURF ID='part surf', RADIUS=2.e-6, GEOMETRY='SPHERICAL', MATL_ID='stuff' /
&MATL ID='stuff',
      DENSITY       = 1.2
      CONDUCTIVITY  = 0.20
      SPECIFIC_HEAT = 1.4/

&VENT DB='ZMIN', SURF_ID='inlet'/
&VENT DB='ZMAX', SURF_ID='OPEN'/

&SLCF PBY=0.5, QUANTITY='VELOCITY', VECTOR=T/
&SLCF PBY=0.5, QUANTITY='TEMPERATURE', CELL_CENTERED=T/
&SLCF PBY=0.5, QUANTITY='DIVERGENCE', CELL_CENTERED=T/

&TAIL/

[drjfloyd]

We don't have a buoyancy term for particles. The particles start falling down and it looks like the divergence issue is happening when a lot of particles collect on the air inlet. Is this just an issue with drag?

[rmcdermo]

I understand about the buoyancy---I was just trying to make a particle that would more or less act like a tracer. I agree, something is odd with particles falling a backward at the inlet. The input parameters here are pretty close to that of particle_flux.fds in the verification suite, but there we just push the particles out of a vertical surface and they fall to the ground.

This is interesting... I change the input file as follows:

&SURF ID='inlet'
      COLOR='RED'
      VEL=-5
      PART_ID='spheres'
      PARTICLE_MASS_FLUX=0.02
      RAMP_PART='ramp pmf'
      NPPC=1
      DT_INSERT=0.01/

&RAMP ID='ramp pmf', T=0, F=0/
&RAMP ID='ramp pmf', T=5, F=0/
&RAMP ID='ramp pmf', T=10, F=1/

Now if you run this case, all is fine for the first 5 sec and Smokeview shows particles being injected, albeit without much randomness. At 5 sec, when the ramp kicks in, the particles look more random, then the instability occurs.

[mcgratta]

Yes, I think this is all about drag. One way to look at it is this -- add solid particles to a single grid cell and blow on it. At first, you expect to see a bit of a flow diversion. Add more particles, and you see more. Eventually, the flow goes unstable. We're just adding momentum without accounting for wake reduction etc.

[rmcdermo]

But at a minimum, something is wrong with the RAMP. There should be no particle mass flux for the first 5 sec.

[drjfloyd]

With a particle density of 1.2 kg/m3, a radius of 2 microns, and a flux of 0.02 kg/m2/s you are injecting a projected area of particles of 15,000 m2 each second in a domain that is 1 m2.

The ramp may just be setting the particle weight to 0.

[rmcdermo]

Yes, that is exactly what is happening. You can see it in Smokeview. Is this the behavior we want? I would think that if PWT<TWO_EPSILON_EB, we should not inject the particle.

[drjfloyd]

We should not inject particles if the flux is zero.

[mcgratta]

I'll fix that.

[mcgratta]

I added T_INSERT to the SURF line parameters to enable one to insert particles starting at a particular time. This same functionality exists for particles introduced within a volume on the INIT line. Otherwise, if a ramp is used to control the particle mass per unit time per unit area (PARTICLE_MASS_FLUX), and the ramp is zero, the code still works and the particles are massless, essentially. It's a bit tricky to remove the particles entirely if the mass ramp is zero, since the particles have already been initialized by the time the ramp is evaluated. I'd rather just leave things as is, and use T_INSERT to delay insertion.

Returning to the issue at hand, these examples fail because the particles create too much drag on the air because we are not currently using any kind of drag reduction algorithm. There is one in the code that was designed for droplets, but we haven't had much experience with it. I'm going to see if we can use a calculated packing ratio to determine if too many particles are being inserted.

[mcgratta]

Run this case with a variety of DT from 0.005 to 0.05. I understand why there is a time dependence, but this is really severe and it's complicating this case alot

&HEAD CHID='part2'/
&MESH IJK=10,10,10, XB=0,1,0,1,0,1/
&TIME T_END=10, DT=0.05 /
&SURF ID='inlet', COLOR='RED', VEL=-0.1 /
&INIT XB=0.4,0.6,0.4,0.6,0.6,0.7, N_PARTICLES_PER_CELL=1, T_INSERT=3., PART_ID='spheres' /
&PART ID='spheres', SURF_ID='part surf', DRAG_LAW='SPHERE'
      QUANTITIES='PARTICLE VELOCITY','PARTICLE DIAMETER','PARTICLE WEIGHTING FACTOR'
      SAMPLING_FACTOR=1/
&SURF ID='part surf', THICKNESS=2.e-6, GEOMETRY='SPHERICAL', MATL_ID='stuff' /
&MATL ID='stuff', DENSITY=1.2, CONDUCTIVITY=0.20, SPECIFIC_HEAT=1.4 /
&VENT DB='ZMIN', SURF_ID='inlet'/
&VENT DB='ZMAX', SURF_ID='OPEN'/
&SLCF PBY=0.5, QUANTITY='VELOCITY', VECTOR=T/
&TAIL/

[mcgratta]

SECOND_ORDER_PARTICLE_TRANSPORT does not change things.

[rmcdermo]

Wow. It flips sign! Interesting find. That does not seem right.

[mcgratta]

I'm going to dig into this. Could be that with particles so small, the approximations and truncations of the derivation are not appropriate. Perhaps there is some expression we can use to choose the sub-time step based on the particle diameter, density and relative velocity.

[rmcdermo]

This looks to me like we have some kind of time splitting bug. We should be able to tell right away, by comparing the vertical component of the drag force and the gravity force, whether a particle should move up or down. We might have to treat these two situations differently. A truly explicit scheme should not care about the order. But we sometimes update variables and then use them again in the routine.

[mcgratta]

In this case, the particle velocity is more or less the upward flow, w_gas=0.1 m/s. The z particle position more or less is governed by

z_new = z_old + w_gas*dt_p - 0.5*g*dt_p^2

If dt_p>0.01 s, z_new is less than z_old. That is why DT=0.05 has the particles falling and DT=0.005 has them rising.

[rmcdermo]

But... particles should not have a gravity term in the position equation. The position evolves by

dx/dt = u_p

The velocity evolves by

du_p/dt = sum of forces (gravity goes here). If done this way, the force balance will always point the particle in the right direction, regardless of time stop (I think).

[mcgratta]

See the notes. This is not how it's done.

[rmcdermo]

I understand. Blame me if you want. I'm just suggesting an approach to try to fix the current issue, which was not realized at the time of the fancy derivation. If we break the 2nd order convergence of the flat fire case, but make the code more stable, that would be fine.

[rmcdermo]

Also, at the time those derivations were done, we were not doing the substepping (as I recall). It is quite possible that the substepping resolves the need for the higher order temporal scheme and the verification cases will all work reasonably well.

[mcgratta]

I implemented the old routine we used to use, where we linearize the equation for the particle velocity. The results of the flat fire and terminal velocity cases is below. Note that there is no second order option with this scheme as yet, so both the first and second order cases are the same. The terminal velocity does not "converge" because it is more or less exact by construction. The terminal velocity position is first order because the scheme is at present first order.

position_convergence.pdf terminal_velocity_convergence.pdf flat_fire_trajectory.pdf flat_fire_u.pdf flat_fire_w.pdf flat_fire_x.pdf flat_fire_z.pdf

[rmcdermo]

Here are the current flat fire results. Is there a way through the LES time step or particle sub time step to get this kind of accuracy? Just curious.

[mcgratta]

I'm looking into it.

On Thu, Oct 15, 2020 at 9:31 AM Randy McDermott notifications@github.com wrote:

Here are the current flat fire results. Is there a way through the LES time step or particle sub time step to get this kind of accuracy? Just curious.

[image: Screen Shot 2020-10-15 at 9 27 36 AM] https://user-images.githubusercontent.com/4418497/96133998-1ba5d400-0ec9-11eb-8d4a-16df262fed8d.png

— You are receiving this because you modified the open/close state. Reply to this email directly, view it on GitHub https://github.com/firemodels/fds/issues/8761#issuecomment-709326318, or unsubscribe https://github.com/notifications/unsubscribe-auth/ACWPCF2KH3OGN56SO4Q7HADSK32R5ANCNFSM4RU73GYA .

[mcgratta]

I made a slight adjustment. When we compute an approximate solution of the trajectory equation, we get something like this in the dragless case:

w(t) = w(0) - g*t

What I was doing was this

z_new = z_old + dt*w(t)

But what I should have done is

z_new = z_old + dt*(w(t)+w(0))/2

That improves things almost to the point where we were before, but with much less baggage:

position_convergence.pdf flat_fire_trajectory.pdf flat_fire_u.pdf flat_fire_w.pdf flat_fire_x.pdf flat_fire_z.pdf

[rmcdermo]

Looks good!