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]

you are right! The exhaust _2 was wrong and now it makes the flow as I want:)
REgarding the ignition, I understand you point that the ignition happens when teh wall temperatrue is over 300 for the surfaces called fire. However, for the surface called burner, I understand it should be automatically ( with 10 seconds delay ) as teh surface is defined as follows: &SURF ID='Burner',
COLOR='RED',
PART_ID='Particle',
NPPC = 10,
VEL_PART=-10,
PARTICLE_MASS_FLUX=1.33
RAMP_PART= 'Particle_ramp_mass_flux'
HRRPUA= 7656.25,
RAMP_Q='Burner_RAMP_Q', THICKNESS(1)=0.3/. Therefore, for me is strange that &OBST XB=0.1, 0.116, 1.8000, 2.1520, 0.0250, 0.1250, COLOR='GRAY', THICKEN=.TRUE., SURF_ID='fire1' /
is getting to an air temperatrue nearby around 600 , while the burner is not yet starting. I also understand that you are recommending to change &DEVC ID='temp_obst_2_1', XYZ=0.108, 1.9760, 0.31, QUANTITY='TEMPERATURE', SETPOINT=300 / to wall temperature?

[drjfloyd]

In your input file, the 'Burner' surface is only being applied to one OBST. All of your other OBST have the 'fire1' surface which requires a surface temperature of 300 C for a wall cell to start burning.

This isn't doing anything in your simulation other than recording the gas temperature at the XYZ in the devc.csv output file. Yes you have set a SETPOINT, but no other input in your fds file uses the output of this device to do anything. Your 'fire1' surface uses the wall temperature for each cell in the OBST to decide if that wall cell should ignite.

&DEVC ID='temp_obst_2_1', XYZ=0.108, 1.9760, 0.31, QUANTITY='TEMPERATURE', SETPOINT=300 /

I suggest you make a very simple input file so that you can test your object igniton.

[maria-justo-alonso]

thanks a lot for your patience!
I didnt use the &DEVC ID='temp_obst_1_1', XYZ=0.108, 2.1520, 0.1250, QUANTITY='WALL TEMPERATURE', IOR=2, SETPOINT=300 / because I was getting this error messsage ERROR(427): DEVC temp_obst_1_1 requires repositioning.

ERROR: FDS was improperly set-up - FDS stopped (CHID: sample1c2e)
And despite changing the coordinates many times, I did not mange to get the positioning, so I moved to the air temperature as an incorrect proxi.
However, I dont understand why if during the 10 first seconds when my single cell 'Burner' is not burning, I get a rise in the air temperature near the obst_1_1 although the heat should only start when teh burner starts to burn after 10 seconds according to the Burner_RAMP_Q. Also after the 10 seconds I expect the 'Burner to start giving away heat ( and thus warming the air around) but I see no change in it.

[drjfloyd]

Please carefully review your inputs:

&RAMP ID='Burner_RAMP_Q', T=5.0, F=0.0/
&RAMP ID='Burner_RAMP_Q', T=6.814195, F=0.010001/

This is not a 10 second delay. T=5 is a 5 s delay

&MESH XB=0, 3.67, 0, 4.69, 0, 2.4, IJK= 60, 80, 60 /,MULT_ID='M'/

3.67/60 = 0.061 m x resolution
4.68/80 = 0.058 m y resolution

&OBST XB=0.4360, 0.4520, 1.8000, 2.1520, 0.2050, 0.3050, COLOR='RED', THICKEN=.TRUE., SURF_ID='Burner'/

The closet grid boundaries to this OBST are 0.4282 m and 0.4893 m

Look at the OBST after this one:

&OBST XB=0.4600, 0.4760, 1.8000, 2.1520, 0.2050, 0.3050, COLOR='YELLOW', THICKEN=.TRUE., SURF_ID='fire1' /

It's closet grid boundaries are 0.4282 m and 0.4893 m. The same as the previous OBST. What does the User's Guide say happens when there are overlapping OBST?

Again I highly suggset you start with a much simpler input file where your opportunities for errors are less. Test one thing at a time. Once you are confident that you know how your inputs work, then add complexity.

[maria-justo-alonso]

Now I understand!! When I change my mesh, suddenly it works as it should!
thanks for all these great advices! I will rebuild the case stepwise as you advice!

[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