Simulation of GAP/HTPB phase behaviors in plasticizers and its application in composite solid propellant

2.1 Theory of DPD

A mesoscale simulation program, dissipative particle dynamic (DPD), provides a dynamic algorithm that incorporates hydrodynamics for studying coarse-grained systems over long length and time scales (19), (20), (21). This technique is frequently used to simulate the phase behavior of polymer blends. For the coarse-grained particle (DPD bead) i and j, the interactions f_i between them can be mainly divided into five terms (22), (23): conservative force Fijc, dissipative FijD and random force FijR, particle i is also confined by bond potential fis and angle potential fiA. All the potential or force exerted on the particles can be expressed as

The terms in the bracket describing the bond or non-bond interaction of the bead i with its neighbor j are given by

where r_ij is the distance between beads i and j; e_ij is the unit vector (r_i−r_j)/r_ij; ω(r) is a function that determines the radial dependence of the force; a_ij is the repulsion amplitude that reflects the chemical characteristics of the interacting particles; γ_ij and σ_ij are the friction coefficient and the amplitude of the noise, respectively; Δt is the time step; and ξ_ij is a random number with zero mean and unit variance. It is noteworthy that a_ij has relation with Flory-Huggins parameter χ_ij, a_ij=a_ii+3.27χ_ij(ρ=3), where a_ii has the value of 25k_BT. Total bond and angle potential energy can be expressed as

where C and k are the coefficient describing the strength of the bond and angle potential, respectively; r⁰ and θ⁰ are the equilibrium bond length and angle, respectively; and r and θ are the bond length and angle, respectively.

2.2 Theory of MD

Molecular dynamic (MD) is generally used to simulate the motions of atoms or molecules. It is also adopted in this article to investigate the thermal diffusion in the blends of polymers and plasticizers. In MD simulations, Newton’s equation of motion (24) is used to calculate the force of atom i:

where F_i is the force, m_i is the mass and a_i is the acceleration of atom i. The force on atom i can be computed directly from the derivative of the potential energy V with respect to the coordinate r_i:

where t is the simulation time. A classical algorithm named as Verlet velocity algorithm (25) can solve the equation of motion. It overcomes the out-of-synchrony shortcoming of the Verlet leapfrog method. The Verlet velocity algorithm is as follows:

v(t+Δt)=v(t)+Δt[a(t)+a(t+Δt)]/2,[11]

where r(t), v(t) and a(t) are respectively the position, velocity and acceleration at time t. This algorithm was adopted in our MD simulations.

For a simulation system at equilibrium, atoms tend to diffuse away from their original location. The mean square displacement (MSD) of these atoms with respect to their original position is obtained as the second moment of their distribution at t>0 and is related to the diffusion coefficient as follows:

where r_i denotes the position vector of ith atom, and angular brackets denote an ensemble average. Then the limiting slope of the mean square displacement as a function of time can be used to evaluate the self-diffusion coefficient of an atom undergoing random Brownian motion in three dimensions.

2.3 Simulations method

2.3.1 DPD method

In our research, two kinds of polymers GAP (X̅_n=37) and HTPB (X̅_n=56), which adopted as the propellant binders, were utilized to simulate their phase behaviors in the presence of DOS or A3/DOS plasticizers. By means of the DPD code in Material Studio 7.0 software (26), we first presented the molecular structures and coarse-grained models of the polymers (Figure 1) and plasticizers (Figure 2). HTPB chains contain different configurations of trans-2,3-butadiene (l), cis-2,3-butadiene (m) and vinyl-1,2-butadiene (n). All of these groups can be represented by H₁, H₂ and H₃ beads (l+m+n=55), respectively. In the case of GAP, the ether and the side azidomethyl groups are the significant difference from HTPB. Soft polyether segments (-OCHCH₂-) make the chains flexible, and therefore we simplified them with G₁ coarse grain. However, the rigid azidomethyl groups (-CH₂N₃) in the side chain decreases the flexibility and weight percentage of the main chains (27); thus, the G₂ bead was used to characterize it. The coarse-grained models of GAP and HTPB are linear and branched, respectively.

Figure 1:

Coarse-grained modeling of GAP and HTPB in DPD simulations.

Figure 2:

Coarse-grained modeling of DOS and A3 in DPD simulations.

A3 is the blend of BDNPF/BDNPA with the mole ratio of 1:1.

The coarse graining method of plasticizers (DOS and A3) follows the similar principle that distinguishes the characteristic groups and general groups. Ester groups in the DOS molecule have relatively stronger polarity and, thus, were simplified as bead D₂. Two kinds of nonpolar alkyl groups in the middle or end site of DOS molecules were represented as D₃ and D₁ beads respectively. A3 is a eutectic mixture of BDNPF and BDNPA. The characteristic groups in A3 are energetic (-NO₂), ether bond (-O-CH₂-O-) and (-O-CHCH₃-O-). Because the coarse grain beads of plasticizers share a similar volume with GAP or HTPB unit beads, we modeled the coarse grains of the nitro and methyl end groups in BDNPF/A as B₂. Lastly, the groups contain ether bond on the middle site in BDNPF and BDNPA were respectively represented as B₂ and B₃.

After the coarse grain models were completely built, the interaction (or repulsion) parameter a_ij, which contribute mainly to the conservative force Fijc in equation (2), should be defined. Groot and Warren (28) associate a_ij with Flory-Huggins parameters χ_ij through the linear equation a_ij(ρ=3)=25+3.50χ_ij. Then we constructed the molecular structures from these “units” mentioned above and computed the χ_ij pairwise interactions of these beads (H₁, H₂, H₃, G₁, G₂, D₁, D₂, D₃, B₁, B₂ and B₃) using the BLEND code in Material Studio 7.0. Ultimately, a triangular matrix was obtained and presented in Table 1.

Table 1:

Calculated repulsion parameter matrix for different beads.

	aij
	G1	G2	H1	H2	H3	B1	B2	B3	D1	D2	D3
G1	25.0
G2	26.1	25.0
H1	32.8	31.4	25.0				Symmetry
H2	33.6	32.2	25.1	25.0
H3	33.8	32.3	25.0	25.2	25.0
B1	35.3	31.0	28.4	29.4	28.4	25.0
B2	25.4	27.5	34.5	35.3	35.1	37.3	25.0
B3	25.2	26.7	32.0	32.7	32.6	34.5	25.1	25.0
D1	36.4	34.4	25.9	25.9	25.7	26.3	38.7	35.4	25.0
D2	34.1	40.4	56.9	58.5	57.9	56.2	33.5	34.6	62.1	25.0
D3	37.7	35.7	26.9	27.1	26.7	27.1	40.0	36.6	25.5	63.2	25.0

Next, in order to explore the phase aggregation or separation behavior of propellant binders in different plasticizers, two periodic cells with a volume of 30×30×30 ų were constructed to pack the binders and different plasticizers for DPD dynamics. For comparison with the experiments, the mass fraction of these binders and plasticizers in our simulations was in conformity with that used in experiments. One cell was filled with GAP, HTPB and DOS with the mass ratio of 0.5:0.5:1.2. The mass ratio of the system containing mixed plasticizers was GAP/HTPB/A3/DOS=0.5:0.5:0.6:0.6. It is worth noting that two types of compounds, BDNPF and BDNPA, account for the same mass fraction in the blend. A time step Δt of 0.05 was used; the r_c and grid spacing were respectively determined as 4.64 Å and 1.0. One hundred and five DPD steps were exerted to relax the coarse grain configurations. At the end of simulations, no fluctuation in pressure or temperature occurred, which means that the systems reached equilibrium. Then the analysis can be performed.

2.3.2 MD method

Thermal diffusion, which is related to the free volume of molecules, also has an influence on the phase behaviors of blends. Thermal diffusion could accelerate the process of phase separation; that is, higher diffusion rate of molecules makes blend get to phase equilibrium faster, or vice versa. We therefore also performed all atom molecular dynamics (MD) simulations (29) to study the diffusion of HTPB and GAP in the blend of GAP/HTPB/DOS and GAP/HTPB/A3/DOS using Forcite code in Material Studio 7.0 (26). To be specific, first, GAP, HTPB, A3 and DOS molecules were separately constructed and relaxed using COMPASS forcefield. Then two amorphous cells (50×50×50 ų) with periodic boundary conditions were built, one was filled with GAP, HTPB and DOS, another one was filled with GAP, HTPB, A3 and DOS. The densities of these cells containing GAP/HTPB/DOS and GAP/HTPB/A3/DOS were 0.99 g·cm⁻³ and 1.13 g·cm⁻³ (practical case). The percentage content of each molecule was the same as the DPD beads. Modeling procedures for MD simulations are given in Figure 3.

Figure 3:

Modeling for the amorphous cells of GAP/HTPB/DOS and GAP/HTPB/A3/DOS in MD simulations.

After the models were built completely, energy minimization was performed to relax the molecular structures by Smart algorithm, the maximal iteration step was 5000. Subsequently, annealing simulations including five cycles between 300 and 500 K for the relaxed structures were carried to overcome energy barriers, and the obtained configuration with the lowest energy was used to execute a 40 ps NVT dynamic process. During the process, the initial velocities of atoms followed Maxwell distributions and were calculated by velocity Verlet algorithm, Coulombic and van der Waals interactions were calculated by Atom based method (30) and Nosé thermostat (31) was used to control the temperature. Once the systems were at equilibrium, the final 12 frames of the results were extracted to calculate the mean square displacement (MSD) curves of HTPB and GAP in two cells.

2.4 Materials and experiments

GAP prepolymer has a hydroxyl group content of 0.55 mmol·g⁻¹ and an average molecular weight of 3700 g·mol⁻¹. HTPB prepolymer has a hydroxyl group content of 0.75 mmol·g⁻¹ and an average molecular weight of 3000 g·mol⁻¹. Desmodur N100 (curing agent I) has an average molecular weight of 728 g·mol⁻¹ and an average functionality of 3.87. The geminal dinitro energetic plasticizer A3 (BDNPF/A) has an acid value of 3.93×10⁻⁵ mmol·g⁻¹ and 0.04% wt. H₂O content, and ester plasticizer DOS has a purity of 99.6%. All the reagents above were purchased from Liming Research Institute of Chemical Industry in China. Toluene diisocyanate (TDI, curing agent II) with a molecular weight of 174 g·mol⁻¹ and functionality of 2, curing catalysts dibutyltin dilaurate (DBTDL) and triphenyl bismuth (TPB), DOS and catalysts were purchased from the Institute of Chemistry Chinese Academy of Sciences in Shanghai China. Cyclotrimethylenetrinitramine (RDX), ammonium perchlorate (AP) with average particle size of 180 μm and 130 μm, respectively, were produced from Xing’an Chemical Industry Co., Ltd., Shanxi. Aluminum (Al) with average particle size of 13.7±3 μm was purchased from Ansteel Group Corp., Liaoning. Binder and plasticizers were dried in vacuum oven for 8 h at 60–80°C before use. A mixture of DBTDL (dissolved in DOS, 0.5% wt.) and TPB (dissolved in DOS, 0.5% wt.) with ratio of 1:2 was used as the curing catalyst. Details of the propellant formulations in our experiments are listed in Table 2.

All the experimental composite propellant mixtures were prepared in 400 g batches in a vertical mixer. The liquid ingredients (except curing agents and catalysts) were charged and stirred for 30 min followed by vacuum mixing for another 30 min to remove entrapped air or moisture. Subsequently, the Al power, RDX and AP were respectively added into the mixture with a time interval of 20 min, followed by vacuum mixing. N100 and TDI were added simultaneously with further mixing of 20 min. Lastly, the catalysts DBTDL and TPB were dropped slowly to the slurry, and then a 10-min mixing was required. A small amount of slurry (40 g) was taken for rheological measurements, which were carried out under 60°C. The remaining slurry was cast into a mold and cured at 60°C for 6 days to obtain a propellant grain.

3.1 Phase behavior of binders/plasticizers

The miscibility or the aggregation morphology of the liquid components has a dramatic effect on the fluidity of the propellant slurry and the molding quality of the propellant grains, whereas the liquid microarea among solid particles determines the momentum transfer efficiency in the two-phase flow process and the distribution of components (27). For instance, if phase separation occurs in the liquid phase of propellant slurry, the fluidity of the slurry may be insufficient for casting (32), and also, molding defects may appear after curing. Therefore, we explored the phase behaviors of the GAP/HTPB binder in different plasticizers using DPD technique and selected the optimal plasticizer for this binder. First, the morphology of GAP (red beads) and HTPB (blue beads) bulks with energy evolution are illustrated in Figure 4.

Figure 4:

Morphology of GAP/HTPB/DOS (A) and GAP/HTPB/A3/DOS (B) blend in DPD simulations.

Beads of A3 and DOS are hidden to facilitate the observation.

As can be seen from Figure 4, the GAP and HTPB beads tend to aggregate as spheres when the DPD cells relax to equilibrium. This aggregation behavior mainly proceeds before 6.7×10⁴ steps. From the view of thermodynamics, the aggregation of the same beads is proved to be a spontaneous process. The spheroidal profile [especially in (B)] may decrease the surface tension between binders and plasticizers. However, the aggregation morphology of GAP and HTPB beads in the two figures shows distinct characteristics. Particularly, two kinds of beads have a more uniform dispersion in A3/DOS than in DOS plasticizer. For quantitative analysis, the concentration of GAP and HTPB beads in the cells is extracted and subsequently plotted in Figure 5.

Figure 5:

Density profiles of GAP and HTPB beads in GAP/HTPB/DOS (A) and GAP/HTPB/A3/DOS (B) in DPD simulations.

The density distribution is counted along the x-axis of cells.

The concentration profiles in Figure 5A and B correspond to the equilibrium condition of the systems in Figure 4A and B, separately. There are some similarities between the two plots in that both the density of GAP and HTPB beads presents a gradient distribution. Specifically, GAP beads aggregate at the distance of about 8 Å and 22 Å, whereas the HTPB beads concentration peaks appear at the distance of about 5 Å, 15 Å or 24 Å, where the content of GAP beads is small.

The reason for distinct density difference between GAP and HTPB in Figure 5A can be attributed to the pairwise miscibility between binders and plasticizers. The high repulsion parameters (a_ij) between the GAP beads (G₁ and G₂) and the HTPB units (H₁, H₂ and H₃) in Table 1 reveal weak interactions between polar groups (azido methyl and polyether units) and non-polar groups (unsaturated alkene groups). In the case of DOS, two kinds of alkyl groups have similar chemical characteristics to alkene groups. As a result, there is strong affinity between them. However, despite the strong polarity of ester group on DOS molecules and azido methyl groups on GAP chains, the mobility of the electronegative azide groups (-N₃) is limited due to its steric hindrance, thus leading to weaker interactions and a higher value of a_ij (40.4) between them. Therefore, phase separation of GAP/HTPB binders plasticized by DOS occurs distinctly.

By contrast, the bead density of two binders exhibits a more uniform distribution in Figure 5B than that in Figure 5A, which indicates that the addition of A3 can enhance the miscibility of GAP and HTPB by a certain degree. Both A3 and GAP molecules contain ether (or similar) groups. In addition, the azido methyl groups interact strongly with ether groups. This results in a good miscibility between GAP and A3; that is, GAP can be dissolved with A3 in any proportion, that is the reason why GAP is generally plasticized by A3 (33), (34). It should be emphasized that although the end of A3 (BDNPF/A) molecules contains polar nitro groups, the symmetry of -NO₂ substitution sites and methyl groups in the ends may enhance the affinity of A3 molecules and the cis- and trans-1, 4-butadiene units on HTPB chains (a_ij of 28.4 in Table 1). In addition, A3 and DOS liquid in our experiments are mutually miscible due to the compatibilizing effect of B₁ units in A3 molecules. As a result, the distribution of GAP and HTPB chains in A3/DOS is more uniform than that in DOS.

Next, we will analyze the effect of thermal diffusion on the phase behaviors of GAP/HTPB/ DOS and GAP/HTPB/A3/DOS blends. Based on the MD simulation results, MSD profiles of GAP and HTPB in the two blends are plotted in Figure 6, and the diffusion coefficient (k) of two polymers can be obtained through linear regressions.

Figure 6:

MSD profiles of GAP (A) and HTPB (B) in two blends in MD simulations.

As can be seen obviously from Figure 6, the diffusion of GAP and HTPB in GAP/HTPB/DOS is faster than in GAP/HTPB/A3/DOS, resulting in a higher k value in the former blend. It indicates that the required time for phase separation between GAP and HTPB in GAP/HTPB/DOS is shorter than that in GAP/HTPB/A3/DOS, and two reasons account for this phenomenon. First internal friction and viscosity of A3 molecules are higher than that of DOS; therefore, polymer diffusion in the blend becomes difficult upon addition of A3. Second, as we discussed above, A3 improves the compatibility of the blend, and the surface tension between GAP and HTPB tends to decrease. The conclusion drawn here is consistent with the result of the DPD simulation. Additionally, compared with the blend plasticized by DOS, the blend plasticized by A3/DOS shows potential applications in solid propellant. GAP/HTPB-based propellants using A3/DOS as plasticizers will therefore be prepared in the next section.

3.2 Rheological study of propellant slurry

Miscibility of liquid phase has a significant influence on the fluidity (or rheology) of propellant slurry; moreover, better fluidity of slurry in the mixing and casting (preparation) process is critical to the molding quality of propellant grains. It is therefore necessary to study the rheological behavior of propellant slurry. In Section 2.1, the advantage of A3/DOS blend in plasticizing GAP/HTPB binders has been corroborated through DPD and MD techniques. In this section, we prepared the propellant slurry plasticized by DOS and A3/DOS for confirming the simulation results and selecting the optimal plasticizer for the GAP/HTPB based propellant. According to formulations listed in Table 2, two propellant slurries were mixed under the same conditions as described in Section 1.3, then slurry samples 1 and 2 were obtained. For convenience of investigation and comparison, the pictures of two slurries are shown in Figure 7.

Figure 7:

Flow behavior of propellant slurry plasticized by DOS (A) and A3/DOS (B).

As can be seen in Figure 7, the solid particles in the slurry plasticized by DOS cluster together. The mixture has poor fluidity and leveling in the PTFE mold until it is under forced vibration. This mixture exhibits a low viscosity but inconsecutive state; that is, solid particles, especially AP and RDX, cannot efficiently adhere to the binders so that the slurry presents a lumpy appearance. This phenomenon can be attributed to the poor physical compatibility between the pairwise units on GAP and DOS. The G₂/D₂ and G₁/D₃ sets exhibit higher repulsion parameters a_ij in Table 1, which leads to the obvious phase separation between GAP and HTPB (Figure 4A) and thus inconsecutive state (Figure 7A). As far as the application is concerned, it is difficult to guarantee a better molding quality of the propellant plasticized by DOS after curing.

By contrast, the fluidity and the leveling property of slurry plasticized by A3/DOS are far superior than that of the former. This can be mainly attributed to the better compatibility between polymers plasticized by A3/DOS as shown in Figure 4B. Thus, GAP chains get more disperse in plasticizers A3/DOS than that in DOS, further leading to a continuous state of slurry in Figure 7B. Therefore, compatibility between the liquid components is demonstrated to have a significant effect on the rheology or processing properties of propellant. Based on the results, the slurry containing A3/DOS was used for a further study in the next section.

3.3 Influence of p_o/p_l on the rheology

In Section 2.2, A3/DOS was determined as the promising plasticizers of GAP/HTPB-based propellant; therefore, high-energy GAP/HTPB-based propellant can be prepared for some further studies. However, before preparation, a quantitative rheological study of this propellant should be implemented to provide rheological parameters for the mixing and casting (or preparation) process. Because the molding quality of propellant depends on its rheological properties, integral propellant grain can be obtained with better rheological properties, but poor rheological properties may lead to some defects in the propellant. Therefore, appropriate rheological parameters are generally required for preparation processes.

As reported in some literature (35), (36), (37), (38), viscosity and flow curves are the most essential graphs to understand the rheological behaviors of a liquid. Slurries with different plasticizing ratio (p_o/p_l) were investigated and presented in Figure 8A and B.

Figure 8:

Viscosity curves (A) and flow curves (B) of slurry with varied p_o/p_l.

It can be seen from Figure 8A that in the shear rate range of 0 s⁻¹ to 30 s⁻¹, the viscosity (η) values on these curves decrease sharply, which resulted from the pseudoplastic and shear-thinning characteristics. However, the variation trend of these curves is different from that of the common pseudoplastic fluid (39), (40). There are step changes on these curves. Initially, the viscosities of these slurries decrease sharply with the increment of the shear rate because the reorientation of polymer chains occurs and the physical entanglements weaken under the shear force, then the η decreases sharply. In this stage, there is a stable momentum transfer in the solid-liquid two-phase flow. With the increasing shearing speed, η of the slurry remarkably decreased with the shear rate ranging from 13 s⁻¹ to 17 s⁻¹, which indicates that distortional flow occurs in slurry. This is mainly due to the asynchrony of liquid flow and solid flow. In addition, the blend cannot be completely compatible (but better than GAP/HTPB/DOS blend). This may result in the discontinuity of the slurry at higher rate. Figure 8B shows the corresponding shear stress with the increase in shear rate. Initially, for each curve, the shear stress increases to its peak because the entanglement and the interactions between molecular chains are weakened by the rotor. Then the stress decreases sharply with the distortional flow.

Moreover, the p_o/p_l has a significant influence on η. A dramatic decrease in η can be seen with the p_o/p_l increasing from 0.8 to 1.1, indicating that the increasing A3 and DOS content can dilute the GAP and HTPB chains. The entanglement and interactions between chains are weakened, but the variation amplitude of η tends to decrease when the p_o/p_l is greater than 1.1 because the intermolecular distance of polymers is long enough.

However, the viscosity (η) itself is not the only influencing factor for the rheology of the propellant slurry. Some other rheological parameters should be considered before casting, such as the non-Newtonian index (n) and yield stress (τ_y). For a comprehensive analysis of the slurry, we also calculated or captured them from the experimental data as presented in Figure 9, noting that the apparent viscosity is extracted from Figure 8A with a shear rate of 1 s⁻¹.

Figure 9:

Rheological parameters of propellant slurry with the variation of p_o/p_l.

(A) Non-Newtonian index (n), (B) viscosity (η) and (C) yield stress (τ_y).

Figure 9 shows a stronger dependence of rheological parameters on the plasticization. With the p_o/p_l increasing from 0.8 to 1.3, the n ascends from 0.31 to 0.45, indicating that the slurry becomes more similar to the Newtonian fluid, which has less dependence on the shear rate. It is noteworthy that the non-Newtonian characteristics are presented by polymers in this blend. A greater p_o/p_l Newtonian fluid is desired in the casting process of the propellant. For the case of η, the variation trend has been discussed above. Here, we are concerned about the value of η, which is an essential parameter for casting, with the p_o/p_l increasing above 1.1, the η ranges from 157 Pa·s to 167 Pa·s, compared to the η of HTPB propellant slurry reported in literature (40), the η of GAP/HTPB propellant slurry at 60°C is small enough for mixing and casting, see Table 3, and is also beneficial to the dispersal and wetting of solid particles. A p_o/p_l range of 1.1 to 1.3 seems to be more appropriate. The last parameter τ_y is obtained from the Casson rheological model (41), (42), revealing the leveling property of the slurry; a smaller value of τ_y causes the slurry to flow under the static pressure itself, i.e. the blend tends to flow spontaneously. As can be seen from Figure 9C, the τ_y decreases dramatically as the p_o/p_l varies from 0.8 to 1.3. τ_y fluctuates with a p_o/p_l range of 1.1 to 1.3 and reaches the minimum 55.23 Pa with the p_o/p_l of 1.2, meaning a better leveling property can be obtained here.

3.4 Molding quality of GAP/HTPB propellant

Based on the rheological study above, a conclusion can be drawn that a p_o/p_l of 1.2 facilitates the mixing and casting process of the propellant slurry. Then we prepared propellant samples using this parameter. After the preparation and curing process, the propellant was sliced for observation as shown in Figure 10.

Figure 10:

Samples of GAP/HTPB propellant after curing.

It can be observed from Figure 10 that the cured samples exhibit a compact characteristic. No defects exist inside the slices, indicating a better molding quality of the propellant plasticized by A3/DOS. This demonstrates that the A3/DOS plasticizer is suitable for the GAP/HTPB propellant. There are many reasons that contribute to this. First, compared to GAP/HTPB/DOS, the GAP/HTPB/A3/DOS blend shows a better compatibility. Next, the blend has a weakened diffusion due to A3, and the required time for phase separation becomes longer. Lastly, a p_o/p_l of 1.2 facilitates the mixing and casting process.