Simulation Development Roadmap MP  HE

1. Simulation Development Roadmap MPHE • High Explosives Applications • Binder and Grain Interactions • Kel-F, Estane • Equation of State • HMX, TATB • Reaction Dynamics and Molecular Processes • Vibrational Analysis of Rapid Post-Shock Events

2. Milestones, HE • Reaction Mechanism for RDX/HMX • TATB/RDX/HMX Molecular Structure • MD of inert and reactive explosives

3. Modeling HE Detonation Muller, Shepherd, Chakraborty, Dasgupta, Goddard Provide a detailed reaction mechanism: 72 Species 445 Reactions Ultimately reduce (Eckett) Understand detonation on a reaction by reaction basis Understand how additional species (binder, aging, etc.) alters detonation properties. HMX

4. Starting Points GRI nitromethane mechanism 49 species 300 reactions Melius nitromethane mechanism 27 species 130 reactions Yetter RDX mechanism 48 species 240 reactions Guirguis nitromethane mechanism & data Chakraborty/Lin reaction modeling

5. Additional Species

6. Thermochemical Data • B3LPY/6-31G** calculation • Compute Cp, H, S from 300-6000 K • Fit to 14-term NASA form

7. HMX/RDX Combustion

8. Validation: Nitromethane Shock Data

9. Chemkin with More Accurate EOS • Ideal gas law poor approximation • Underestimates volume • Overestimates density, reaction rates, factor of 15 (?) • CV uses Ideal Gas EOS • Can’t describe high P,T • Obtained CK-Real Gas • BKW, vdW, Peng-Robinson • BKW EOS:

10. Extension • New decomposition scheme for RDX, HMX based on ab initio DFT calculation • Incorporation of detailed reaction kinetics for early stage C, H, N, O reactions • CH2N reations with NOx, N2O, OH etc. • Decomposition of CH2NNO2, CH2NNO etc.

11. Unimolecular decomposition of RDX

12. Unimolecular decomposition of RDX

13. Unimolecular decomposition of RDX

14. Unimolecular decomposition of HMX

15. Unimolecular decomposition of HMX

16. TATB (1,3,5-triamino-2,4,6-trinitrobenzene) Molecular Structure QM Calculations DFT/B3LYP • Nitro-benzene - completely planar structure - due to resonance of the NO2 groups with the benzene ring • Aniline - non-planar NH2 group due to lack of resonance • Geometry is in good agreement with MW data

17. RDX Molecular Structure

18. Nitro-benzene PES from DFT; torsion fit for force field • Rigid rotation PES for C-NO2 bond • Fit by 2- and 4-fold torsion terms V2 = 3.81 V4 = -0.31

19. Aniline PES from DFT; torsion fit for force field • Rigid rotation PES for C-NH2 bond • Inversion barrier to planarity at N: 1.3 kcal (0), 2.7 kcal (90) • PES curve corrected for inversion and then fit by 1-, 3-, 4- and 5-fold torsion terms V1 = 4.86 V3 = 1.93 V4 = -0.88 V5 = 0.17

20. Ab Initio PES for Inter-molecular H-bond • Planar aniline::nitro-benzene scan • BSSE corrected E

21. TATB Work in Progress • Fit force field H-bond potential to map ab initio PES • Use C-N exocyclic torsions from fits to ab initio PES • Adapt remaining force constants from benzene and nitro-methane work • Calculate quantum corrected P,V,E surfaces (EOS) from NVT MD • Combine with Kel-F to study binder-grain interaction

22. MD and Force Field Development for HMX • Level 0 - Generic Force Field (Dreiding) - Quantum Corrected EOS from NVT MD using • Cerius2 using Ewald method for non-bonded interactions • Level 1 - Vibrationally accurate force field Developed for DMN, HMX & RDX • DFT (B3LYP) calculations on Dimethyl-nitramine (DMNA) • FFOPT parameterization of FF to fit QM frequencies and normal modes

23. Pressure Dependence of Crystallographic Cell Parameters for b-HMX (Theory using Dreiding, open symbols) with Experiments from Cady & Ollinger (filled symbols) B C A

24. Elastic Constants and Bulk Modulus of b-HMX Experimental Data (Joe Zaug, LANL), Theory (Dreiding)

25. Quantum Corrected EOS from NVT MD • Classical MD gives incorrect EOS because all modes are equally populated • Quantum correction from Velocity Auto Correlation gives proper EOS

26. Cv from NVT MD

27. Molecular structure for DMNA • DFT calculation gives non-planar structure as a minimum - planar structure has -117cm-1 NC2 inversion mode • Experimental MW data indicates planar DMNA, possibly due to fluxional averaging • Barrier to planarity is 0.51 kcal/mol from DFT • Inversion mode at 225cm-1 for non-planar structure

28. Optimization of FF to reproduce vibrational frequencies

29. Comparison of Structural Parameters for DMNA

30. Impulse dynamics – the redistribution of internal vibrational energy • 2 DMNA molecules in crystallographic orientation • All H atoms in molecule 1 imparted a 50A/ps impulse towards O atoms of neighboring DMNA molecule • Redistribution of energy monitored in various modes • N-N bond gains maximum energy of 25kcal/mol and has periodicity of 80fs

31. Polymer binders - Overview • Many polymers are used as binders in a variety of High Explosive systems. In this regard, studies of materials & impact properties of these materials are important. • Of particular interest is Kel F-800 (3M) and Estane (BF Goodrich). • Estane is a polyurethane and work has just started on this material. Explosive Crystal Polymeric Binder

32. 0.157 0.066 -0.361 Cl F H F 0.800 C C C C -0.318 0.786 0.204 F F H F 3 -0.391 1 -0.348 n Total charge = 0.00 Polymer binders - Overview Kel-F 800 • Kel F-800 is a random copolymer of chlorotrifluoroethylene and vinylidene fluoride monomer units in a 3:1 ratio. • The presence of the vinylidene fluoride disrupts the crystallinity of the chlorotrifluroethylene to form an essentially amorphous polymer. • Although amorphous, the polymer is very dense due to the presence of the Cl and F atoms. • Quantum Mechanical (Hartrees Fock, 6-31G** basis set) calculations were done using Jaguar to obtain Mulliken and ESP charges.

33. EoS from Atomistic simulations for Kel-F • Strategy - use pure polymer systems • Force field choice and optimization. • The choice of Force Field is critical for these systems • Charges • In order to get good charges, Quantum Mechanical calculations were done on representative oligomers and electrostatic potential derived charges were used in the polymer. This is particularly important for halogenated polymers. • Packing • To best represent the amorphous polymer, the packing of the chains have to be optimized. We have used in-house developed algorithms for this purpose and tested by comparing with experimental Cohesive Energy Density for a variety of polymers. • Simulations • The Molecular Dynamics simulations was run using in-house developed MD and commercial software to obtain data for Equation of State.

34. Kel F-800 - Building • The polymer chains are built under periodic boundary conditions. • The regular “Amorphous Builder” in Cerius2 builds using RIS, but does give proper density for Kel-F, most likely because the packing is not optimal. • The MSC has developed a module (in Cerius2) to calculate optimize packing for calculating Cohesive Energy Density (CED) of polymers. • This builds amorphous cells starting with half the expected density. • The system is minimized and NVT dynamics run for a short period. • The dimensions of the unit cell are slightly compressed and the minimization and dynamics steps repeated. • The unit cell is repeatedly compressed, minimized and NVT dynamics run until experimental density is obtained. • The net result is a much better representation of the packing in amorphous polymers.

35. Kel F-800 - Packing 24 monomers - 2 chains 24 monomers - 16 chains • Using 2 chains causes alignment within the unit cell giving regions of crystallinity. • Using more chains in the unit cell overcomes this problem. • Two open questions are how many chains to use and how long must each chain be? 20 Å 50 Å

36. Kel F-800 - CED Validation • Initial choice of force field and validation studies were conducted on pure PCTFE, where experimental data is available. • For further validation, a whole range of polymers have been studied and their CED values compared to experiment. • For the similar PCTFE, the calculated CED is well within the experimental range, as is the case for PVC and PMMA. • This lends confidence that the FF for Kel-F is reasonable and the calculated CED value should compare well with experimental value for Kel-F. Validation Solubility parameters 12 11 10 9 Solubility parameter/ (cal/cc)1/2 8 7 6 5 PIB PS PE PVC PVEA PMMA PCTFE Kel-F Polymer type Experimental value/range Calculated value

37. Kel F-800 Chain number dependence of Solubility Parameter 9.3 9.2 9.1 9 Solubility (cal/cc)1/2 8.9 8.8 polydisperse 8.7 8.6 4 8 12 16 20 24 Number of chains Kel F-800 - dependence of CED on number of chains • CED value is independent of number of chains in unit cell provided there are sufficient molecular dynamics steps for equilibration.

38. Kel F-800 - effect of poly-dispersity Shortest chain • Dispersity: • 2  12 mer • 5  18 mer • 10  24 mer • 5  30 mer • 2  36 mer } 24 chains in total 3504 atoms An ethylene group represents a monomer unit Molecular weight distribution 30000 Mn = 2471.74 Mw =2640.37 25000 (exp = 29360) (exp = 75660) Longest chain 20000 15000 Weight of polymer Volume = 5.169x10-3 m3.kg-1 10000 Density = 2.0 g.cm-3 (3M experimental value) 5000 0 12 18 24 30 36 Length of chain

39. Kel F-800 3504 atom case PVE curves at T = 0K Volume 200000 70 Energy 60 150000 P = 0Gpa => E min 50 100000 40 Energy Pressure (GPa) (cal/mol) 50000 30 20 0 10 -50000 0 -100000 -10 40791 26379.3 30131.1 33825.1 37395.4 44019.1 46621.6 48423.6 50907.6 6.82E+05 1.01E+06 Volume (A^3) Kel F-800 - FF Validation Cold Compression • Dreiding -EXP6 force field is the force field of choice - works in both Cerius2 & MPSim. • Energy is at a minimum with no external pressure. Pressure vs Volume Cold Compression 70 60 50 40 Pressure /Gpa 30 20 10 0 0.5 0.6 0.7 0.8 0.9 1 V/Vo

40. Kel F-800 Pressure-Volume Isotherms Pressure vs Volume Equation of State - 3504 atoms (100ps, data averaged over the last 30ps) 80 100K 200K 70 300K 400K 60 500K 600K 50 Pressure /GPa 40 30 20 10 0 0.5 0.6 0.7 0.8 0.9 1 V/Vo

41. Kel F-800 Velocity Auto Correlation • Compute quantum thermodynamics from classical MD. • Uses the velocity autocorrelation function to get the vibrational Density of State. • From this we get the thermodynamic properties using the harmonic approximation.

42. Kel F-800 - Density of States Simulated IR spectrum of single chain Density of States from VAC for 6 chain unit cell case

43. Kel F-800 Energy-Volume curve Kel F-800 Energy vs Volume 3000 100K 2900 200K 2800 300K 2700 400K 2600 Energy /kJ/kg 500K 2500 600K 2400 2300 0.7 0.75 0.8 0.85 0.9 0.95 1 V/Vo

44. Kel F-800 EoS • Solids can be modeled by the Mie-Gruneisen equation of State. • The definition of the Gruneisen parameter G is given by the thermodynamic relationship • G is a function of Volume and Temperature. • A “cold” approximation to G makes it only a function of volume

45. Kel F-800 Kel F-800 Gruneisen Parameter 4000 100K 200K 3500 300K 3000 400K 2500 500K [P(v)-P(c)/E(v)-E(c)] /(Gpa.kg/kJ) 2000 600K 1500 1000 500 0 1900 2100 2300 2500 2700 1/V (kg/m^3)

46. Gruneisen Parameter 100K 1400 200K 1200 300K 1000 [P(v)-P(c)/E(v)-E(c)] /(Gpa.kg/kJ) 400K 800 500K 600 600K 400 Kel F-800 200 Gruneisen G vs Temperature 0 2.65 1900 2100 2300 2500 2.6 1/V (kg/m^3) 2.55 2.5 2.45 Gruneiscen parameter G 2.4 2.35 2.3 2.25 2.2 0 100 200 300 400 500 600 700 Temperature /K Kel F-800 Gruneiscen Slope of this line gives G (dimensionless)

47. Kel F-800 Heat Capacity, Cv 1.1 1 0.9 Cv /kJ/(K kg) 0.8 0.7 0.6 0.5 200 250 300 350 400 450 500 550 Temperature / K Kel F-800 Heat Capacity Exp Tm = 378K Exp Tg = 303K

48. Large scale MD simulations of HE • High explosives: HMX,TATB grains in a polymeric binder (kelF, estane) • Model HE System • periodic in all 3 directions • x- and y-direction extent (10 nm) • z-direction extent (1-10 mm) • Allow axial compressive loading using new steady state algorithms (Cagin, Goddard) polymer grains polymer 10 nm 1-10 mm