Simulation of solid propellant microstructures by combining the collective rearrangement method with the discrete element method

The polydisperse particulate components in solid propellant are incompact and randomly packed, which determines the microstructural features of the propellants. A packing method, combining the discrete element method (DEM) and collective rearrangement method, was applied to model propellant microstructures. The validity of this method was investigated by comparing the calculated and experimental properties of the monodisperse, bidisperse, and polydisperse random close packed sphere systems. The propellant models were generated using a stepwise approach, and their homogeneity, local randomness, and long-range pattern were analyzed. A statistical study of aluminum (Al) particle distribution was also conducted. The results indicated that this packing method can effectively determine the microscopic characteristics of random close packed monodisperse spheres. The maximum packing fraction of bidisperse and polydisperse spheres had similar trends to those reported in experimental studies and using other packing algorithms. In addition, this method was capable of generating non-compacted propellant structures with uniformly distributed polydisperse particles. The radial distribution functions (RDFs) for Al-Al particles provided information about the Al distribution, but this was mainly related to the size and content of the large particle components.The polydisperse particulate components in solid propellant are incompact and randomly packed, which determines the microstructural features of the propellants. A packing method, combining the discrete element method (DEM) and collective rearrangement method, was applied to model propellant microstructures. The validity of this method was investigated by comparing the calculated and experimental properties of the monodisperse, bidisperse, and polydisperse random close packed sphere systems. The propellant models were generated using a stepwise approach, and their homogeneity, local randomness, and long-range pattern were analyzed. A statistical study of aluminum (Al) particle distribution was also conducted. The results indicated that this packing method can effectively determine the microscopic characteristics of random close packed monodisperse spheres. The maximum packing fraction of bidisperse and polydisperse spheres had similar trends to those reported in experimental studies and using other packing...


I. INTRODUCTION
Solid propellants are heterogeneous composite materials filled with particulate components and polymer binders. Many propellant properties are closely related to their microstructures, such as aluminum (Al) agglomeration, 1,2 mechanical properties, 3 and burning rate. 4 Therefore, obtaining detailed structural information for solid propellants is critical for a comprehensive understanding of their performance. However, it is a complicated and challenging task to reconstruct propellant heterogeneity, even using tomography. 5,6 Due to the rapid advances in computational power, numerical simulation is assumed to be a promising alternative. In recent decades, various packing methods have been proposed and are capable of generating 2D 7,8 and 3D propellant packing models. [9][10][11] Based on these models, several Al agglomeration, propellant combustion, and propellant mechanical models have been established. 3,[11][12][13][14] The current 3D models are typically generated by a particle expansion method, such as those proposed by Lubachevsky-Stillinger. 15 However, the packing fraction of these models can be affected by the growth rate of particle size. It has been reported that higher growth rates reduce the packing capability, while slower growth rates enable higher packing fractions. 9,10 Practically, the polymer binder, plasticizer, and other continuous components in a composite solid propellant typically account for between 10-20% of the propellant. Therefore, the particulate components are incompact and the packing fraction is clearly determined by the propellant formulation. We considered that a packing method using the actual sizes and contents of particles according to the propellant formulation to construct a propellant packing model would be more appropriate. The packing method reported in Ref. 16, which can be classified in the collective rearrangement category, [17][18][19] seems to be appropriate. The spherical particles are first generated ARTICLE scitation.org/journal/adv according to a given number and size distribution in a domain with many overlaps. Then, a relaxation iteration is applied to reduce or eliminate the overlaps. The packing process is completed once the mean overlap value is small enough. Using this method, any packing fraction (lower than the dense packing limit) of particle packing models can be generated by the direct input of given parameters. However, although there are several approaches to eliminate overlaps by using the collective rearrangement method, they are typically based only on geometric calculations and are challenging for non-professional programmers. The discrete element method (DEM) is a dynamic method for dealing with the problems of collisions and motion between granular materials. 20 The usefulness of the DEM in particle packing has been demonstrated in many studies. [21][22][23][24] More importantly, there are several general purpose discrete-element modeling software packages (PFC, 25,26 EDEM), with dynamic algorithms embedded that can be adopted for particle packing. Most of the DEM-based packing methods in the literature focus on dense packing, such as gravitational deposition, 22 centripetal packing, 27 and compression. 28 They are not suitable for modeling incompact propellant structures.
Therefore, a packing method combining the DEM and collective rearrangement method was applied to simulate the propellant microstructures. The particle generation and calculations were implemented using the PFC3D 5.0 software, with precise control of the packing fraction. The method was validated by studying monodisperse, bidisperse, and polydisperse particle packs. Three series of propellants were modeled according to literature formulations, and their microstructures were analyzed. To avoid boundary effects, 24,29 all packs were generated with periodic boundaries.

II. PACKING METHODS
The packing processes were mainly implemented using the "ball distribute" command in PFC 3D 5.0. 25 The particles were represented by balls. Specifically, the ball size distribution, volume fraction of each type of ball and the model porosity (packing fraction) are pre-set. Then, balls are placed randomly in the model domain according to the pre-set parameters. There were a large number of overlapping particles at this stage ( Fig. 1(a)) because there was no check for overlaps when a sphere was placed. Repulsion and collision between spheres were attained by iteratively using dynamic calculations based on the DEM. This procedure ceased when the relative deviation of the average contact distance was less than 0.05 between each of the 5,000 steps (see section II C). Finally, rearranged spheres which met the target criterion were generated ( Fig. 1(b)). The dynamic calculations in the DEM were cyclic processes between each step, including the iterative application of a forcedisplacement law at the contacts and Newton's second law to particle motion.
A. Force-displacement law Only particle-particle interactions were considered because the model domain was periodic. The magnitude of the overlap and the relative displacement at the contact point were related to the contact force. The simplest linear contact model was adopted without regard to shear force and dashpot components. The contact plane was centered within the interaction volume of the two pieces (Fig. 2). The contact location Xc and normal contact force direction nc were determined as follows: where, R and X are the sphere radius and sphere center, respectively; δn is the contact distance, which is the minimal signed distance separating the sphere surfaces such that the gap is negative if the surfaces overlap; and d is the distance between the center of two spheres. Then, the normal linear contact force at the contact position is: where, kn is the normal stiffness coefficient. During a time-step Δt, the relative displacement increment at the contact point is given by Δδn=δn −(δn) 0 , and Fn is updated using Eqs. (1)-(5).

B. Law of motion
According to the resultant force and Newton's law of motion, each particle position could be determined at the end of each time step Δt. The translational motions of particles in the DEM were solved via the central difference method.
A damping-force term was added to the equation of motion to remove additional kinetic energy. The damped motion equation can be written as: where, F d is the damping force, which is proportional to the Fn, and X is the acceleration.
where, α is a damping constant, andẊ is the particle velocity. The acceleration at time t can be described using the 1/2 step velocity, Putting Eq. (6) into Eq. (8) we obtained: Finally, the particle position at t+Δt can be written as: The positions of the particles and contact forces acting upon them were updated at the end of each step. In this study, the time steps were chosen automatically by PFC3D using a time-stepping algorithm. 25 Because only a structural model was required in this work, a detailed physical representation was not necessary. In fact, the value range of kn and α that can achieve this purpose is very wide. The values of kn is related to the speed of particles while α is related with the energy dissipation efficiency. A proper value of α should be between 0 and 1 (α > 0). Therefore, kn =1000 and α =0.9 were selected because experiences show that they can eliminate most of overlaps at a relatively fast speed. In this work, for particles with density ρp = 100, radius r =1, kn = 1000 and α=0.9, the timestep increases from 0.043s/step to 0.4s/step after 5000 steps, indicating the particles gradually slowed and model tend to be static.

C. Stopping criteria
The determination of when to stop a packing process is an open issue. The average contact distance (d, total contact distance divided by total contact number) decreases as the calculation process proceeds and becomes vanishingly small when particles are not closely packed. In this study, d was checked every 5,000 steps, and the relative deviation (δ) of d was defined as follows: δ = d − dp dp (11) where, dp is the average contact distance of the previous check. It can be inferred that the calculation can be stopped when δ is small enough. Therefore, a critical value of δ (δc) was selected as the stopping criteria (i.e., calculations stopped when δ < δc). The δc value was determined by investigating d when calculations were stopped under different values of δc, and the effects of model size and δc on d were compared (Table I). The sphere radius r = 1, and particle density ρp = 100 were chosen to generate packing models, and all parameters were dimensionless. Table I indicates that d values were almost the same when δc decreased from 0.1 to 0.001, and the effects of model size and packing density on d were insignificant. However, a larger number of calculations were required when δc decreases, especially for crowded packing. Therefore, a δc value of 0.05 was selected in this study in consideration of the computation efficiency and accuracy.

III. RESULTS AND DISCUSSION
A. Monodisperse packs

Packing fraction
Numerous studies have focused on the random close packing (RCP) of equal sized hard spheres, but the packing fraction is still not well-defined. 30 It is widely accepted that the packing fraction varies from 0.60 to 0.68 and is method dependent. 31 Considering that d will increase significantly once the packing limit is exceeded, it is possible to measure the maximum packing fraction by investigating the inflection point of d.
The pre-set packing fraction (ρ0) of the models varied from 0.60 to 0.70, while the radius (r) of spheres and the length of the cubic domain were set to 1 and 30 (dimensionless), respectively. The number of particles produced varied from 3,867 to 4,512 depending on ρ0. The results are shown in Fig. 3. Figure 3 shows that the average contact distance became vanishingly small when ρ0 was less than 0.635. It increased sharply (by about 10 ∧ 8 times) from 1.75e-11 to 2.16e-3 when ρ0 increased from 0.635 to 0.64. The d value continuously increased with ρ0 when it was higher than 0.635. This implies that the particles became crowded and overlaps became apparent when the pre-set packing fraction was around 0.635. In other words, the RCP packing fraction of singlesized particles was about 0.635 using this method. This result was consistent with those reported by other studies. 31,32 Further proof, as depicted in Fig. 3, was provided by the relative deviation of the packing fraction ε, which is defined as follows: where, ρ0 is the prescribed packing fraction and ρ is the real packing fraction when packing is completed. The ε value was relatively constant when ρ0 < 0.64 and increased when ρ0 > 0.64. This means that the overlap volume among particles will continue to grow when ρ0 is above the packing limit, and the model becomes invalid. It should be noted that when packs were not dense, ε was lower than 0.025%. This deviation was so small that it could be ignored when compared with the whole system; thus, the following analysis did not distinguish between ρ and ρ0.

The coordination number (CN) and radial distribution function (RDF)
To validate the RCP models, certain microscopic properties are usually measured and compared with experimental results. The most commonly used properties are the CN (or average near-contact number) and RDF. 24 For comparison with other literature results, 0.64 was used as the dense packing fraction.
The CN refers to the average number of particles that are in "contact" with a certain particle. The "contact" is defined by a tolerance e. If the center distance between two spherical particles is less than (1+e) D, where D is the diameter of particles in an equal sized pack, then they should be treated as being in "contact". Obviously, the CN varies significantly with the choice of e. (ρ = 0.640) packs were compared. The CN increased as the e increased, and the denser the pack, the higher the CN became. It should be noted that unlike many other packing algorithms, e = 0 was possible in this study where overlaps were allowed. Therefore, a dramatic increase in the CN from 1.00 to 1.01 occurred when ρ = 0.626, which indicated that overlaps diminished significantly when the pack was loose. Table II presents a comparison between the results from this method and data selected from Ref. 10, and shows that good agreements were obtained. Although the CN was an average value, the effects of pack size and measurement region were negligible. The RDF g(r) represents the probability of finding a particle (or local density of particles) at a given distance from a reference center. It provides information about the local concentration of particles and is widely used in particle packing studies. The normalized RDF g(r) is defined in Eq. (13): where dN(r) is the number of particle centers within a distance between r and r+dr from the reference particle center, and φ is the number density of particles in the model domain.
In Ref. 10, the RDFs of POLIPack and Aste's experimental results are compared. In the present study, the RDFs of packs were also compared with the results mentioned above using the overprint method. For the sake of clarity, the RDF data in this study are presented in a scatter form in Figs. 5(a), (b) and their line charts are presented in Fig. 6, in which the noise data were processed for comparison. The horizontal axis was normalized with respect to the particle diameter. Figure 5(a) shows that the g(r) curve for the DEM-based pack agreed well with reported results when the packing fraction was 0.64. The curve displays two separate peaks at distances of √ 3 and 2.0, which are recognized as key short-range features of the random packing of hard particles. The peaks correspond to two types of particle connection: edge-sharing in-plane equilateral triangles and three particles along a line. When the packing fraction was 0.626, two smoothed peaks were apparent (Figs. 5(b) and 6), which was due to the slight positional uncertainty of particles when the pack was relatively loose. Compared with POLIPack, the g(r) in this work had a smoother peak at √ 3 and was more similar to experimental data, which may be a consequence of the allowance of overlaps.

B. Bidisperse packs
Modeling of non-uniform particle packing is necessary because solid propellants are filled with blended particle components. In experiments conducted by McGeary, 33 the binary packing of spherical particles was simulated using two packing methods similar to those used in Ref. 10. The calculated packing density was then compared with the experimental results to validate these methods. We attempted to verify our packing method in the same way. Two series of binary packs with different fine particle volume fractions were generated, with coarse-to-fine particle ratios of 6.5:1 and 3.4:1, respectively. The length of all packs was constant and was set to be 100; thus, the particle number in each pack varied from about 20,000 to 100,000. The maximum packing fraction for each condition was obtained using the same strategy as in section III A 1, with a minimum interval of 0.005 for the pre-set packing fraction. The results are shown in Figs. 7(a), (b).  Bidisperse packing is able to attain a higher packing fraction than monodisperse packing because fine particles can fill the spaces between coarse components. The maximum packing fraction a binary packing can achieve is related to the coarse-to-fine particle ratio and relative content of the particles used. As shown in Fig. 7(a), when the coarse-to-fine particle ratio was 6.5:1 and the fine particle fraction was 25%, the maximum packing limit was 0.799. When compared with other simulation methods, a good accordance with experimental results was observed when the fine particle content varied from 25-30%. However, our method produced an overestimation in other cases, especially when the fine particle fraction was less than 25%. This maybe because greater weight is given to the contacts of small particles when calculating average contact distance. Therefore, although the average contact distance is vanishingly small, the overlaps of large particles may still be high, resulting in an overestimation of the packing fraction. For packs with a coarse-to-fine particle ratio of 3.4:1, all packing algorithms overestimated the packing limit compared with the experimental results ( Fig. 7(b)). The results for our packs were closer to Rocpack (a = 0.2), which used a slow expansion rate for particle size, enabling a higher packing fraction. 29

C. Polydisperse packs
In practical applications powder blends are polydispersed. Ammonium perchlorate (AP) and Al powder are the most common fillers in metallized propellant. The experimental and model results for the RCP fraction of AP/Al binary particles were compared. The Al and AP powders were first mixed in different ratios and then placed in a cylindrical glass vessel. The filling process was accompanied with vibration and shaking until the packing volume of the mixed powders could not be less than 20 ml. The powders inside the vessel were then weighed and the packing fraction was calculated using Eq. 14 where, m is the total weight of powder; f Al is the weight fraction of Al; and the ρ Al and ρAP values were set as 2.7 and 1.95 g/cm 3 . It should be noted that although the particle density was not accurate it did not affect the comparison of experiments with models. The packing models were generated according to the size distribution of the powders, which were measured by a Malvern Mastersizer and the values obtained are presented in Fig. 8. In addition, because the particle sizes were not uniform, the model size had to be large enough to reduce the deviation between the prescribed particle volume fraction and the results generated by the model. However, large models will dramatically reduce the calculation efficiency. Therefore, we set the model side length to 2.0-2.6 mm for Al/AP binary packing, 1.0 mm for Al packing, and 10 mm for AP packing. This resulted ARTICLE scitation.org/journal/adv FIG. 7. Comparison of the maximum packing fraction between different bidisperse packing studies. (a) coarse-to-fine particle ratio = 6.5:1; (b) coarse-to-fine particle ratio = 3.4:1.

FIG. 9.
Comparison of the maximum packing fraction of aluminum/ammonium perchlorate (Al/AP) powder between calculated and experimental results.
in the deviation of the particle volume fraction in each pack being less than 1%, and the total particle number being less than 1,000,000. The maximum packing fraction for each pack was obtained by determining the total volume fraction of particles when the average contact distance increased drastically. The comparison between experimental and model results under different Al/AP ratios is shown in Fig. 9. Figure 9 shows that the model results followed the same trend as the experimental results under different AP/Al mass ratios. The calculated maximum packing density occurred when the AP/Al was about 3/1, which was consistent with experiments. This indicated that the packing density of this model was able to imitate the powder filling in solid propellants. However, the model results were always higher than the actual situation. This may be due to the non-spherical shape of the Al and AP powders and the slight compression between particles in the models.

Spatial structure of particle components
Previous studies have demonstrated the effectiveness of this method in modeling random close packed sphere systems. In this study, the method was applied for modeling solid propellant microstructures and the model propellant "A" data from Ref. 13 were adopted. On a volume basis, the propellant consisted of 52.41% 200 μm coarse AP (cAP), 11.75% 30 μm Al, and 9.04% 10 μm fine AP (fAP).
Because particles are dispersed randomly and evenly in practical propellants, it was necessary to investigate the distribution of particles in the model packs. Interestingly, the direct use of dynamic calculations to rearrange all particles in the model resulted in an   FIG. 11. Variation in the area fraction of fine ammonium perchlorate (fAP) along the z-axis (direct approach). inhomogeneous structure. Figure 10 presents the slices in the X-Y plane and through the center of the pack of propellant "A". It can be seen that small particles were locally concentrated, with many voids in the packs, especially for fAP. The variations in the area fraction of fAP cross sections obtained by the X-Y plane continuously passing through the cube are depicted in Fig. 11. The results indicated that there were more fAP particles at the model edge. This was because the mass of large particles was much larger than that of small particles, which caused small particles to be pushed toward the edge of the model by large ones due to particle repulsion forces. In addition, the propellant was loosely packed, which means that overlaps could be eliminated even when the particles were not evenly distributed. Therefore, small particles were concentrated around the model edge. To solve this problem, the packing method was modified into a stepwise form. Large particles (cAP) were generated first and 50,000 steps of dynamic calculation were applied to reduce most of the overlaps and dissipate energy. Then, small particles (fAP and Al) were generated and rearranged according to the previous stopping criteria. The model results are presented in Figs. 12, 13, and 14. It can be seen that the fAP particles were evenly distributed in the pack section, and the variation in the area packing fraction of fAP displayed no obvious tendency for the particles to concentrate locally. Therefore, this stepwise approach was more suitable for the incompact packing of polydisperse particles.
Three different random seeds, which controlled the randomness of the initial data, were used to generate three parallel packs. The variations in the area fraction of the AP and Al cross sections obtained by the X-Y plane continuously passing through each cube are plotted in Figs. 15(a), (b).
Although the same information was used to investigate the three parallel packs separately, the results were quite different from each other. Figure 15 shows that the area fraction of AP and Al particles varied around their volume fractions (0.6145 and 0.1175, respectively) although there was no general tendency observed in the variation. These results indicated that the propellant packing models generated by this method were locally random.
The partial RDF for Al-Al particles (g Al-Al (r)) calculated from the propellant packing model showed the statistical long-range pattern of particle distribution. Figure 16 shows the averaged g Al-Al (r) values for three parallel packs. Although a different packing method was adopted, the curve in Fig. 16 was in good agreement with the literature results. 13 The first peak was linked to the local g(r) maximum, globally around 60 μm, which was characterized as the most frequent distance between neighboring Al particles. The minimum interval of the function was around 150-160 μm. This indicates that cAP particles generally existed at this distance from the reference center, which resulted in a decrease in the density of local Al particles.

Statistical analysis of the Al concentrating region
The partial RDF for Al-Al particles mentioned above provides information about the spatial distribution of Al particles. According to previous studies 7,11,13 the minimum position of the function when g Al-Al (r) <1 can be regarded as the border of the statistical "pocket". The "pocket" is the interval region formed by coarse oxidizer particles. The Al particles inside the "pocket" tend to coalesce into agglomerates (see "pocket model" 2 ). In addition, the low point of g(r) is closely related to the size and content of the coarse particles in the propellant. Therefore, the g Al-Al (r) values of the two series of practical propellants were investigated and the effect of the cAP and fAP content and size on the Al distribution were considered. Example A. Hydroxyl-terminated polybutadiene (HTPB) propellants with differences in the size of cAP particles. The first example used propellants 1, 3, and 4 in Ref. 34. Three propellant packs with differences in their cAP size (250, 350, and 450 μm) were generated. The size distribution of Al and AP particles were consistent with Fig. 2 in Ref. 35. The pack length was 3,000 μm, the AP, Al, and HTPB binder contents were converted into volume fractions, and their densities were set to 1.95, 2.70, and 0.92 g/cm 3 , respectively. About 320,000 particles were generated. Five different regions per pack were analyzed, and their RDF for Al-Al particles were averaged and are plotted in Fig. 17. Figure 17 shows that with the increase in cAP size, the radius interval moved to a higher position when g Al-Al (r) <1. This indicates that the Al concentrating area was enlarged when the cAP size increased. The experiments conducted in Ref. 34 showed that the average Al agglomerate size also increased with the increase in cAP size. This indicated that the RDF results were in accordance with the calculations based on the "pocket" model, i.e., the "pocket" volume was enlarged with an increase in the coarse oxidizer particle size, leading to the formation of large agglomerates. Example B. Polybutadiene acrylic acid acrylonitrile (PBAN) propellants with differences in the size and content of fAP particles.
The second example was based on Ref. 36. Propellant packs with differences in fAP sizes (49 and 82.5 μm) and cAP/fAP contents (70/30, 80/20, and 90/10) were generated. The length of each pack was 2000 μm. The cAP were assumed to follow a normal distribution over the range of 355-425 μm, with a mean of 390 μm. The fAP were assumed to follow a normal distribution over the range of 75-90 μm, with a mean of 82.5 μm and over the range of 45-53 μm, with a mean of 49 μm. The Al particles were assumed to be same size, with a diameter of 30 μm. The density of the PBAN binder was 0.936 g/cm 3 . The average RDF (g Al-Al (r)) of each pack is presented in Figs. 18(a), (b). Ten runs per pack were analyzed. According to the "pocket" model, fine oxidizers with a sufficient size may have a package effect on Al powders, which will decrease the effective "pocket" volume. The experimental study in Ref. 36 found that the average size of an Al agglomerate decreased with an increase in the fAP content, while the inverse trend was apparent when the fAP was 82.5 μm. Figures 18(a) and (b) show that with the increase in the relative content of fAP, the radius range of each pack moved in a higher direction when g Al-Al (r) <1. This made sense because a reduction in cAP content would lead to an enlargement in the spacing between cAPs; thus, the "pocket" range constructed by the cAP would increase. Considering that an 82.5 μm fAP particle was 2.8 times larger than an Al particle, it can better pack Al particles than fAP with a diameter of 49 μm. However, there was little noticeable effect of fAP size on RDFs, and the decrease in agglomerate size when the fAP content increased was inconsistent with the RDF-determined Al concentrated region. Therefore, the RDF curves only indicate the effect of cAP content on the Al particle distribution in this example. The decrease in Al agglomerate size may be due to the existence of fAP pockets (which cannot be identified by g Al-Al (r)), and the easier ignition of Al due to the local fuel excess when the fAP content increases. 36 Based on this, we infer that RDFs are capable of roughly presenting the concentrated region of Al particles, but are not necessarily able to accurately measure the average "pocket" size.

IV. CONCLUSIONS
A packing method combining the DEM and collective rearrangement method was used to simulate RCP and incompact propellant microstructures. The model results agreed well with both experimental and other simulation methods.
(1) For monodisperse packs, the relative deviation of the packing fraction was less than 0.025% compared with the pre-set value. Their maximum random packing fraction was about 0.64, which agreed with the results of previous studies. The coordination number (CN) of packs was close to the experimental results and those obtained using other packing methods. The RDF curve displayed two peaks at distance of √ 3 and 2.0, which are key short-range features of the random packing of equal sized hard spheres.
(2) For bidisperse and polydisperse packs, the dense packing fraction of different compositions had a similar trend to those observed experimentally and when using other packing algorithms. (3) The propellant packs were generated using a stepwise approach, which was suitable for packing loose and polydisperse particles. The particulate components in the model propellant packs were locally random and long-range patterns could also be observed. (4) The partial RDF for Al-Al particles based on realistic propellant formulations could provide information on the size of the Al-rich region surrounded by cAP particles, while the size and content of fAP particles had no obvious effect on RDFs.