Numerical instability in simulating solid particle ejection

[ahsadeghi]

Hello,
Following the below discussions on FDS google group
https://groups.google.com/g/fds-smv/c/iYC4bYUagkk
https://groups.google.com/g/fds-smv/c/iT7dmOoJuCQ/m/CxZIStqyAAAJ

I still have problems with particle ejection. Thanks to Kevin and dr_jfloyd for the comments, but the solutions did not work for my case. Because I cannot change the particle properties, particle diameter, mass flow rate, and vent area, as they are from experiments and I need them for validating the simulations. I have tried different ways such as decreasing DRAG_COEFFICIENT, defining VEL_PART, changing MAX_PRESSURE_ITERATIONS, including SUSPEND_PRESSURE_ITERATIONS=T, but they didn't help. I wonder if you would please let me know if you have any solutions. Thank you for your consideration. Please find attached the input file. It is noteworthy that thermal conductivity is high to have the lumped capacitance assumption.

Model_small.txt

Regards,
Hosein

[drjfloyd]

What version of FDS are you using?
What is it exactly you are trying to do with the particles? What behavior are you expecting to see and what behavior are you seeing?

[ahsadeghi]

The version is FDS-6.8.0-1238-g64f51df
I am trying to simulate particle ejection from a 18650-type battery. I expect the HRR, particle mass flow, the total particle mass of 4 g and the particle size shown in the figure below.

the majority of these data are from the following paper:
https://www.sciencedirect.com/science/article/pii/S0378775323007334

Regards,
Hosein

[drjfloyd]

How can you know you are not getting that mass flux and size distribution? Your input file has no outputs measuring anything about particles. Your input file has no outputs specified at all.

&MATL ID='Particle_mat',
DENSITY=1.0,
CONDUCTIVITY=1e10,
SPECIFIC_HEAT=0.8,
EMISSIVITY= 0.9/

The DENSITY is not realistic. This is the density of air. A solid particle is going to have a density much larger than this. The paper gives 2000 kg/m3, If your goal is to replicate the paper, you should be using the values in the paper unless you have a good argument as to why those values are not correct.

You don't need such a high CONDUCTIVITY. Non-physical values of parameters can lead to solver issues. FDS sets the wall node size using the equation in 8.3.8 of the guide. You only need CONDUCTIVITY ~0.1 W/m/K for FDS to have just one wall node.

I would suggest you start with a much simpler input file. Just inject monodisperse particles from a surface and figure out how to measure the mass of particles being injected. Then maybe add a fire to that surface as well. Then maybe add a CNF distribution. Then maybe add time depdendent ramps for the the fire and particle flux.

&VENT ID='Burner', SURF_ID='Burner', XB=0.032,0.042,0.03,0.04,0.066,0.066/

You have a 1 cm x 1 cm surface for your particle flux. While this may be roughly the size of an 18650 cell, unless the top of the cell blew off during the experiment, the particles and gas flow are coming from a smaller area. Using a larger area will reduce the injection velocity.

You are using a CNF_RAMP. CNF stands for cummulative number fraction. The plot in the paper you are basing the ramp off of is a mass fraction plot. The two are not the same. Say you have 10 um particles and 100 um particles with a 50 % mass fraction in each. A 10 um particle has 0.1 % of the mass of a 100 um particle. In number fraction this means you are not 50 % 10 um and 50 % 100 um but 99.9 % 10 um and 0.1 % 100 um.

Why are you setting N_STRATA=1?

&REAC ID='Reaction1',
FUEL='VENTED GASES'
FORMULA='C1H4'/

Why 100 % methane your fuel (METHANE is also a predefined species and you should use that rather than the formula)? The paper states: "the composition of the gas mixture generated by the NCM LIB with 50%SOC is 20.28% H2, 38.65% CO, 30.16% CO2, 4.66% CH4, 5.9% C2H4 and 0.35% C2H6"

[ahsadeghi]

Thank you for your comment. I have done a simpler case before this case, please see the following input file, this case has a bigger surface area and the boundaries are not open to measure the total mass of the particles (I know it can affect the HRR but here I was verifying the particle part)

particle_test.txt

sorry the correct density is 2000 kg/m3, I was testing with 1 kg/m3 to see if the problem is with the weight of the particles.
I created the CNF_RAMP ramp in the input file based on the figure, the data in the figure were not used in the input file directly. Then I designed the ramp for the mass flow rate to outcome a total mass of 4 g. There is a device in the newly attached input file to measure it.

I had N_STRATA=1 because of having CNF_RAMP, similar to the example in section 15.4.4. However, removing it does not help.

I had CH4 for simplicity, I included the original mixture but there is still the numerical instability. Please see the input file below.

Model_small.txt

I have been working on this problem for a couple of weeks but only decreasing the mass flow rate of the particles to 0.0000133 (in Model_small.txt input file) helped, which is not desired. This amount seems like not having particles. Otherwise, changing other parameters only changed the instability time around 0.1-0.2 seconds.

Regards,
Hosein

[drjfloyd]

The cell vent is 9 mm2 and you are modeling it at 100 mm2. The particles are going to be coming out at a decent velocity from the real vent. I added a VEL_PART, turned on PARTICLE_CFL, and injected more particles (NPPC=10). I think your problem was with a small number of particles being injected at near zero velocity that they were collecting on the top and eventually the resulting drag made things go unstable.

Model_small.fds.txt

[rmcdermo]

Is there any reason why the particles would not be exiting the VENT at the same velocity as the fuel gas? To me, that would mean ramping the part velocity as well. The fuel velocity is (HRRPUA/HOC)/RHO_FUEL, which should peak around 0.3 m/s = (10,000 kW/m2 / 50,000 kJ/kg) / (0.667 kg/m3). But maybe I'm missing something.

[drjfloyd]

In the real world there is a very tiny vent hole on the battery top that fails first. It is much smaller than the overall battery diameter. In the FDS input the VENT is placed over the entire battery top which means the battery gases are at a much lower velocity than if they just came from the vent hole.

Also at early times the RAMP for the flows is near zero.

[rmcdermo]

OK, but in the real case, with a much smaller orifice, still the gas and particles are at roughly the same velocity, I would think. I am just saying that we should probably keep that relative velocity in the model if we want to avoid spurious drag. Yes, accurate account of the jet momentum is lost once you smear it out to 10 times the area.

[ahsadeghi]

Thank you @rmcdermo and @drjfloyd. The simulation is still running with the changes you made @drjfloyd. Yes, you are right, I should fix the vent area.

[maria-justo-alonso]

Hello,
This is a kind of follow up question
I am beginner in FDS.I am simulating the ventilation in a battery room. I would like to set up the thermal runaway of a battery module where the cells of a module enter in TR when the neighboring wall of the cell has a temperature higher than 380C. I know that it is possible to set a timer. But is possible to set a minimum temperature for starting the fire?
Then as there are two parts of the room very differentiated, one with a combustion and the other where there is only the mixing of particles due to ventilation I guess it is better to have different mesh's sizes. Would you have any advice how to do that?
Thanks in advance! Maria

[mcgratta]

You can use IGNITION_TEMPERATURE and a HRRPUA ramp on a SURF line to define the boundary conditions of your battery pack. When the surface temperature reaches the specified ignition temperature, the object will burn with the given HRRPUA.

[maria-justo-alonso]

So would you make every cell in the battery pack as something like this:
&OBST XB= 0.05, 0.15, 1.8, 1.9, 0, 0.15, COLOR='RED', SURF_ID='fire1' /
&SURF ID='fire1',IGNITION_TEMPERATURE =380, HRRPUA=250 /
And then one mesh containing the battery pack and another for the rest of the room?

[maria-justo-alonso]

Hi again. I have managed the delayed ignition in the cell, however, I need some help with the other cells in teh pack. I have defined them as objects as it follows:
&MATL ID='STEEL',
FYI='Drysdale, Intro to Fire Dynamics - ATF NIST Multi-Floor Validation',
SPECIFIC_HEAT=0.46,
CONDUCTIVITY=45.8,
DENSITY=7850.0,
EMISSIVITY=0.95/

&SURF ID='fire1', TMP_INNER=40.0,
MATL_ID(1,1)='STEEL',
MATL_MASS_FRACTION(1,1)=1.0,
THICKNESS(1)=6.0E-4, IGNITION_TEMPERATURE =300, HRRPUA=0.380513346 /

&OBST XB= 0.4475, 0.8000, 2.1742, 2.2000, 1.3998, 1.5000, COLOR='BLUE',THICKEN=.TRUE., SURF_ID='fire1' /

However, it seems that there is no heat transfer to these objects. They did not increase the temperature when the neighboring cell is burning and even when the room air has temperature over 300 close to the burning cell, the other cells do not change the temperature at all. Am I defining the object wrong and should define it in a different manner? Thanks!

[mcgratta]

Try without TMP_INNER and make sure you have an external heat source that will obviously change the temperature. If that does not work, create a simple test case and start an issue.

[maria-justo-alonso]

I have removed it as recommended, however, the temperature in the objects is not changing. When I plot the slices and I ask to display only the contours of the objects, I have no color information within the objects outlines, the slices are empty within the object.

[drjfloyd]

An SLCF of TEMPERATURE when you have 1D heat tranfers is a gas phase only output. The volume covered by the OBST is a solid; there is no gas there. See the table of solid phase outputs in the User's Guide.

[mcgratta]

I do not suggest you mix meshes at first. Try something simple to see if the concept works. Then you might explore grid refinement. But I warn you that putting one mesh inside another is a fragile construct and generally avoid it.

[maria-justo-alonso]

Ok, I will try that then!Thanks a lot!

[ahsadeghi]

Dear @drjfloyd
Previously you mentioned here that CONDUCTIVITY ~0.1 W/m/K is enough for FDS to have just one wall node (lumped capacitance assumption), i.e., Bi=hL/k <0.1. Would you please explain how it is done in FDS? Do you have an approximation for h?

Regards,
Hosein

[drjfloyd]

8.3.8 Solid Phase Numerical Gridding Issues