The presence of spalled particles might a ect the ow eld and the material, thereby in uencing the aero- dynamic heating rates of the thermal protection system. In order to study the impact of particles on the ow eld, a two-way coupling is performed between a Lagrangian particle trajectory model and a hyper- sonic aerothermodynamics ow solver. Time-accurate solutions are computed for argon and air ow elds. Single-particle and multiple-particle simulations are performed and results are studied. The studies in the argon environment indicated that the particles start sublimating soon after they cross the shock and keep sublimating during the remainder of their time in the ow. In the air environment, the particles start re- leasing carbon vapor soon after their ejection, and the magnitudes of di erent carbon products vary based on the particle's location relative to the shock. Similar behavior is observed for multiple-particle simulations but with an increase in the amount of released carbon vapor that di uses and convects over a larger area. A single-particle simulation is also run by adding additional gas phase reactions, and it is found that the production rate of certain carbon products (C 1and CN) increases. A parametric study is conducted based on parameters that a ect the motion of the particles. The results of the comprehensive study show that the carbon vapor released by spalled particles tends to change the composition of the ow eld, particularly the upstream region of the shock, which a ects the heat ux incident on the test sample. Keywords: Spallation, Ablation, Thermal Protection System, Time Accurate Solutions, Lagrangian Particle Trajectory, Particle-laden ows, Aerothermodynamics Nomenclature Symbols A Face area, m2 E Speci c total energy, J/kg Corresponding author Email addresses: raghava.davuluri@uky.edu (Raghava S. C. Davuluri), zhanghb28@mail.sysu.edu.cn (Huaibao Zhang), kaveh.tagavi@uky.edu (Kaveh A. Tagavi), alexandre.martin@uky.edu (Alexandre Martin) Preprint submitted to engrxiv October 3, 2022 Eve Speci c vibrational-electron-electronic energy, J/kg FD Drag force, N FcFcFc Convective ux matrix FdFdFd Di usive ux matrix h Enthalpy, J/mol Ji Mass di usive ux, kg/(m2{s) M Mach number _mC=i Mass production rate due to surface reactions, kg/s mp Mass of the particle, kg Mw Molar mass, kg/mol n Normal face vector P Primitive variables vector p Static pressure, N =m2 pdrag Drag power, J/s Q Conservative variables vector q Heat ux vector _q Heat transfer rate, J/s R Residual vector S Source term vector in KATS t Time, s T Temperature, K (u;v;w ) Velocity components in axial, radial and z-directions, m/s U State vector V Volume of the mesh cell, m3 W Source vector (x;y;z ) Position components in axial, radial and z-directions, m  Density, kg/m3 Viscous tensor _! Mass production rate in KATS, kg/s Subscripts conv Convection f Flow eld g gas phase 2 p Particle rad Radiation rxn Chemical reaction s solid phase tr translational-rotational energy mode ve vibrational-electronic energy mode 1. Introduction Ablative thermal protection system materials shield entry vehicles from high heat loads through the 5 combination of heat and mass transfer phenomena usually referred to as \ablation." As the surface of the material undergoes chemical transformation, mostly through oxidation, it has been observed that solid particles are ejected from the surface into the ow. This process { typically included in the ablation umbrella term { is de ned as \spallation." It is believed that the thermal, mechanical shear and inner pressure stresses cause the ejection of the particles. The spalled particles are likely produced by disconnected bers or chunks 10 of material, although they could also be formed from soot, a by-product of the pyrolysis of the resin that constitutes heat shield material. The presence of spalled particles around the material is believed to a ect the ow eld, and thus the surface heating rates since the near-surface chemistry is modi ed. Moreover, the particles' ejections take place before the material fully oxidizes, thereby increasing the heat rates at the surface and accelerating the 15 material recession. The particles travel through a high-temperature zone and re-radiate energy back to the surface as they heat up, thus increasing the heat ux onto the capsule. Through these phenomena, spallation a ects the ablative material directly. Additionally, the particles also tend to react chemically and physically with the species in the ow along their paths, thus changing the chemical composition of the ow eld, and hence altering the radiative heat ux. 20 Discrepancies were found when experimentally measured heating rates and temperature pro les of Pioneer- Venus [1] and Galileo Probe [2, 3, 4] heat shields were compared with the theoretical predictions. The value of the parameters was over-predicted near the stagnation point and under-predicted in the downstream region. Although the experimental results were re-evaluated by varying the e ect of turbulence, no reasonable expla- nation was found for the discrepancy. The results suggested that additional mechanisms caused an increase 25 in radiation or turbulence in the downstream region. Spallation was believed to be one of those mechanisms. Furthermore, when laser attenuation and emission spectroscopy measurements were performed by Raiche and Driver [5], the results were unexpected as it was found that with an increase in heating rates, the optical attenuation increased for a PICA sample. This unexpected phenomenon was speculated to be due to the 3 scattering, absorptive and re ective nature of spalled particles. Also, the optical emissions in the inviscid 30 region corresponded to black body radiation at about 3800 K and were likely due to the presence of spalled particles. Spectroscopic measurements by Yoshinaka [6] on ablating models in air ow eld demonstrated the presence of CN emission spectra in the upstream region of the surface. Similar spectroscopic experiments by Kihara et al. [7] in the region ahead of the shock detected the presence of CN emission spectra in the nitrogen and air environments and the presence of C emission spectra in the argon environment. The most 35 likely reason for the presence of CN and C in the upstream region of the shock is due to spalled particles, as they are the only carbon source that can reach a distance that far from the sample. Based on above experimental observations, numerical models [8, 9, 7, 10] were developed to compute the dynamics of spalled particles. These models used the extracted ow eld data to study the spallation behavior. Pace et al. [11] developed an Eulerian particle model of uniform mass, and loosely coupled it to 40 Computational Fluid Dynamics (CFD) code DPLR (Data Parallel Line Relaxation) [12], to evaluate the impact of spalled particles on radiative heating. Following these studies, a Lagrangian particle trajectory model [13] was recently developed to compute the properties of the spalled particles. The code takes into account the chemical reactivity of the particles with the ow eld species. The code was one-way coupled in previous work to the Fluid Dynamics (FD) 45 module of Kentucky Aerothermodynamic and Thermal-response Solver (KATS) [14] to evaluate the e ect of ow eld on the spalled particles. In order to assess the e ect of spalled particles on the ow eld, a two-way coupling between the spallation model and KATS-FD model is performed and is presented here. The coupling is performed for a single particle and for multiple particles in argon and air environments. Additionally, a parametric study is conducted for single particle simulation in air environment to summarize 50 the e ect based on the parameters that govern the particle's motion. It should be noted that the trajectories presented in the present work might not represent the ones from spallation experiments [15, 16, 17]. However, the main objective of the present work is to assess the impact of particle dynamics on ow eld and not necessarily to replicate the experimental trajectories. It is for the same reason only gas-particle interactions and no other ablative mechanisms are considered in this work. 55 2. Methodology 2.1. Numerical models 2.1.1. KATS { FD The thermo-chemical non-equilibrium ow eld in the continuum regime is computed using KATS- FD [18], a laminar hypersonic aerothermodynamic modi ed Navier-Stokes CFD solver. The governing equation of the model is of the form: @Q @t+r(FFFF dFdFd) =Sf (1) 4 where Qis a vector of conservative variables, FFFandFdFdFdare convective and di usive ux matrices, respectively, andSfis the source term vector. The weak form of the governing equation (Eq. (1)) is obtained by integrating it over a nite volume Vfor an arbitrary mesh cell and is given by: Z V@Q @tdV+I A(FFFF dFdFd)ndA=Z VSfdV (2) where the Gauss theorem has been applied to the ux term. Assuming uniform physical properties within the control volume, the left-hand side term of the equation is integrated with respect to time by employing the rst-order backward Euler method, and integration of the ux integral is performed by adding the uxes across each surface. The nal form of the equation is V t@Q @Pn @R @Pn P=Rn: (3) where Pis a vector of primitive variables,  tis the time step size, and Ris the residual vector which is expressed as: RX j2cell(FdFdFdFFF)njAj+VSf (4) where nandAare the face normal and face area respectively. The vectors of conservative variables, primitive variables, and source terms are of the forms: Q=0 BBBBBBBBBBBBBBBBBB@1 ... ngs gu gv gw E Eve1 CCCCCCCCCCCCCCCCCCA; P=0 BBBBBBBBBBBBBBBBBB@1 ... ngs u v w T Tve1 CCCCCCCCCCCCCCCCCCA; Sf=0 BBBBBBBBBBBBBBBBBB@_!f 1 ... _!f ngs 0 0 0 0 _!f ve1 CCCCCCCCCCCCCCCCCCA(5) whereiis the density of species i, (u;v;w ) are the components of bulk velocity, and EandEveare the total energy and vibrational-electron-electronic energy per unit mass characterized by temperatures Tand Tve, respectively. _ !f iis the mass production rate of species iand _!f veis the vibrational energy transfer rate between two di erent energy modes. Subscripts 1 to ngsrepresent the species number, while subscript g 5 represents the gas mixture. The convective and di usive ux matrices in Eq. (1) are given as: FFF=0 BBBBBBBBBBBBBBB@1u  1v  1w ......... ngsu  ngsv  ngsw gu2+p  gvu  gwu guv  gv2+p  gwv guw  gvw  gw2+p (E+p)u(E+p)v(E+p)w1 CCCCCCCCCCCCCCCA;FdFdFd=0 BBBBBBBBBBBBBBB@J1 ... Jngs xx xy xz yx yy yz zx zy zz uqPngs i=1(Jihi)1 CCCCCCCCCCCCCCCA(6) wherepis the static pressure, is the viscous tensor, Jiis the di usive ux of species i, and qis the heat ux vector. The code employs second-order spatial discretization and rst-order temporal integration. 60 2.1.2. Lagrangian particle trajectory model The model simulates the dynamics of a spalled particle by employing a Lagrangian formulation [19]. The model also includes the chemical interaction of the particle with the ow eld. The governing equation of the model is @U @t=W (7) where Uis the state vector and Wis the source term vector. The elements of the vectors are represented as: U=0 BBBBBBBBB@mp mpup mpvp mpwp mpEp1 CCCCCCCCCA; W=0 BBBBBBBBB@_mC FDx FDy FDz _qconv+pdrag_qrad+ _qrxn1 CCCCCCCCCA(8) wherempis the mass of the particle, ( up;vp;wp) are velocity components of the particle, Epis the speci c total energy of the particle, _ mCis the mass source term, FDx;FDy;FDz are components of drag force acting on the particle, and _ qconv, _qrad, _qrxn,pdrag are convective, radiative, reaction heat rates, and power drag. The reactions at the particle surface considered in this model are oxidation, nitridation, and sublimation. 65 The reactions are irreversible in nature and are given in Table 1. The mass recession rates due to oxidation and nitridation are computed by employing parameters from Driver's chemistry model [20, 21], whereas the Knudsen-Langmuir equation [22] is used to calculate the surface recession rate due to sublimation. Therefore, the nal form of the mass source term can be expressed as: _mC=_mC=CO_mC=CN_mC=C1;C2;C3 (9) 6 Table 1: Surface chemistry model considered for Spallation model Type Reaction Oxidation C(s) + O ! CO Oxidation 2 C(s) + O 2! 2 CO Nitridation C(s) + N ! CN Sublimation C(s) ! C1 Sublimation 2 C(s) ! C2 Sublimation 3 C(s) ! C3 where _mC=idenotes the mass recession rates due to the reaction producing species i. 2.2. Two-way coupling method The two-way coupling scheme is performed by adding a source term vector Spto Eq. (1), which accounts for the particle's contributions that a ect the uid. Therefore, the modi ed form of the governing equation of KATS-FD can be expressed as @Q @t+r(FFFF dFdFd) =Sf+Sp (10) The elements of Spare represented as Sp=0 BBBBBBBBBBBBBBBBBB@_!p 1 ... _!p ngs _!p mom;x _!p mom;y _!p mom;z _!p tr 01 CCCCCCCCCCCCCCCCCCA(11) where _!p iis the mass rate of species produced or consumed due to surface reactions, _ !p mom;i is the rate of change of momentum, and _ !p tris the rate of change of total energy of the particle. 70 2.2.1. Solution technique Initially, a steady state solution is computed for the ow eld environment using KATS-FD. The ow eld parameters are then used by the spallation model to calculate the dynamics of the particle. The two-way coupling is achieved by inserting the elements of source term vector Winto the source term vector Spof the CFD code. The coupling is performed according to the following: 75 7 Mass coupling The products released from the surface reactions of the particle tend to a ect the composition of the ow eld environment. The mass recession rate due to the spalled particle is inserted into the CFD code as: _!p i=_mC=i V(i= CO;CN;C1;C2;C3) (12) whereVis the volume of the mesh cell in which the spalled particle is present [23, 24]. The mass coupling also accounts for the change in concentrations of oxygen (atomic and molecular) and atomic nitrogen responsible for the particle surface reactions, which are given as: _!p O=MwO MwC_mC=CO V(13) _!p N=MwN MwC_mC=CN V(14) _!p O2=0:5MwO2 MwC_mC=CO V(15) whereMwiis the molar mass of species i. Momentum coupling The motion of the spalled particle is a ected by the drag force acting on it. Hence, momentum coupling deals with reaction of drag force on the ow eld. The elements of Spcorresponding to the momentum terms determined from the spallation model are: _!p mom;i =FDi V(16) whereirefers to the x-,y-, andz-directions. Energy coupling The energy rate terms of the spalled particle from the vector Wwhich a ect the ow eld are _ qconv, pdrag, and _qrxn. The convective heat rate ( _ qconv) accounts for the heat transfer rate between the ow gases and the spalled particle. The heat energy released or absorbed due to particle reactions per unit time ( _ qrxn) a ects the temperature of the ow eld. Also, the work done by the particle to overcome the drag force ( pdrag) a ects the kinetic energy of the ow eld. Hence, the energy coupling is performed by inserting the energy rate terms from WintoSpas: _!p tr=_qconv+pdrag+ _qrxn V(17) It is assumed that the radiative heat rate from the particle impacts the ablative material, and is less e ective on the ow eld. Hence, _ qradterm is not considered in the energy coupling. 8 Starting from the initially converged CFD solution, KATS is used in a transient mode, and the spallation 80 source terms are added to the cell centers as given by Eqs. (12) { (17) along the path of the particle. To maintain a time-accurate solution, the Courant-Friedrichs-Lewy (CFL) number in the CFD code is kept under 1. 2.2.2. Cell location algorithm The source vector Spis computed at the cell centers of the mesh, whereas the spallation model provides 85 results based on nodal properties. A new algorithm was developed to locate the center of the mesh cell in which the spalled particle is present. Using the algorithm, the source terms calculated by the spallation code are added to the center of the cell in which the particle is present. Figure 1 illustrates the trajectory of a 30m particle ejected normally from the surface, 5 mm from the center axis, at an initial velocity of 90 m/s in the air ow eld. The cell centers located by the new algorithm, denoted by red squares, are shown at 90 three di erent locations with respect to the position of the particle. Figure 1: Cell-centers located with respect to the trajectory of the particle 2.2.3. Veri cation To verify that an accurate coupling procedure was employed, a constant spallation source term was added to a zero velocity ow eld. An integration was performed over the cell volumes to evaluate the mass 9 deposited by the spalled particle in the ow eld. It was veri ed that the mass deposited in the ow eld was 95 equal to the total mass added from the source term le (constant source term time step sizetotal time steps). Since the mass added by the spallation code was equivalent to the mass deposited in the ow eld, the applied coupling method was considered veri ed. Figure 2 shows a time-accurate solution of zero velocity ow eld when a 14 m particle is ejected with a velocity of 370 m/s normally, 13 mm from the center axis. The particle deposits a source term of magnitude 1018kg/s throughout its travel. The particle trajectory 100 takes 2697 iterations of time step size 1 107s to complete. Figure 2 corresponds to the solution after 600 time steps. The values of numerically calculated mass (mass integrated over cell volumes), theoretical mass, and their relative error for various iterations are given in Table 2. Figure 2: Carbon vapor deposited in a zero-velocity ow eld after 600 time-steps The mass deposited by the particle is transported through the ow eld by advection and di usion. The chemical vapor produced is advected in the direction of the bulk velocity of the uid, and is di used in 105 the direction of the density gradient. In this particular case, since it is a zero velocity ow eld, the mass deposited is not advected and is only propagated by di usion along the particle trajectory. 10 Table 2: Veri cation results Iterations Numerical Mass (kg) Theoretical Mass (kg) Absolute error (kg) 0 1.83E-36 0.00E+00 1.83E-36 100 1.00E-23 1.00E-23 2.28E-26 200 2.00E-23 2.00E-23 8.70E-28 300 3.00E-23 3.00E-23 8.22E-27 400 4.00E-23 4.00E-23 2.14E-26 500 5.00E-23 5.00E-23 3.69E-26 600 6.00E-23 6.00E-23 3.75E-26 700 7.00E-23 7.00E-23 4.38E-26 3. Results and discussion The time-accurate two-way coupling solutions were computed for Mach 5 high enthalpy Argon and Air environments. The boundary conditions used to compute the steady state solution for the argon environment 110 are given in Table 3. These conditions correspond to the ow of speci c total enthalpy of 8 :06 MJ/kg and result in a cold-wall heat ux of 246 W/cm2and pressure of 23 :343 kPa at the stagnation point. It should be noted that the ionization reactions of argon are not included for the simulation, resulting in a ow eld solution with higher temperatures and greater shock stand-o distance. Table 3: Uniform ow conditions of Mach 5 argon ow f(kg/m3)U1(m/s)Ttr(K)Tw(K) 3:451033860:0 1470 500 The ow conditions considered to compute the steady state solution for the air environment are presented 115 in Table 4. These conditions yield a ow of speci c total enthalpy 5 :6 MJ/kg that would result in a cold wall heat ux and pressure of 74 W/cm2and 7:155 kPa, respectively, at the stagnation point. Table 4: Uniform ow conditions and mole fractions of species of Mach 5 air ow f(kg/m3)U1(m/s)Ttr(K)Tve(K)Tw(K)YArYN2YO2YNOYOYN 1:491033162.7 896 896 600 0.047 0.686 0.051 0.042 0.000 0.173 Simulations for a single particle as well as for multiple particles were performed. For the case of the single particle simulation, the coupling was performed in three stages: only mass, mass and energy coupling, and 11 total (mass, momentum, and energy) coupling. This systematic procedure helps track changes in ow eld 120 properties. For the multiple particles simulation, only the total coupling was performed. 3.1. Single-particle simulation 3.1.1. Argon ow eld The mass decrease of the spalled particle while traveling in the argon environment only occurs through sublimation. The ejection parameters and the total physical time taken by the particle to travel through the 125 computational domain are given in Table 5. The physical time mentioned hereafter denotes the time taken by the particle to travel from its ejection point until it exits the computational domain. Table 5: Ejection parameters and physical time for single particle simulation in argon environment Size ( m) Velocity (m/s) Position (mm) Angle () Physical time (ms) 14 370 (0,13,0) 0 0.2697 To ensure time accuracy, the solution was computed with the time step size of 4 109s, which cor- responds to a maximum CFL of 0.824. The results for mass coupling of the particle given in Table 5 are illustrated in Figs. 3, 4, 5, and 6. These gures consist of a temperature pro le and density pro les (C 1, C2, 130 and C 3) at di erent travel times. The temperature pro le is used to indicate the position of the particle with respect to shock location, whereas the density pro les are used to determine the mass of the vapor deposited by the spalled particle. The animation of this simulation is provided in the Supplementary Material. Sublimation is only a function of the temperature of the particle. Due to the very low thermal conductivity of argon, the thermal interaction between the particle and ow eld is negligible. At 0 :07 ms into the particle's 135 travel, its temperature remains relatively low, resulting in negligible recession. This can be observed in Fig. 3. The temperature of the spalled particle increases while it crosses the shock. At about 0 :15 ms of travel time, the particle starts sublimating, as seen in Fig. 4. Owing to a high temperature, the released vapor is di used rapidly into the ow eld. Additionally, the bulk velocity of the uid guides the di used vapor along the downstream region of the shock. 140 At 0:20 ms of travel time, the particle is in the upstream region of the shock and is about to re-enter (Fig. 5). With an increase in the particle's surface temperature, it sublimates more vapor, thus increasing the concentration of the sublimed vapor. This di used carbon vapor increases in magnitude when the particle is in the downstream region, imme- diately after re-entering the post-shock region. Figure 6 shows the particle at 0 :26 ms travel time, with the 145 additional concentration of carbon species around it. It can be seen from the above simulation that the order of magnitude of di used vapor decreases from C3species to C 1species. It can be concluded that when the particle sublimates, the concentration of C 3 12 (a) Temperature (b) C 1Density (c) C 2Density (d) C 3Density Figure 3: Trajectory of the particle in an argon ow eld at 0.07 ms with heterogeneous carbon reactions. (An animation of this gure is provided in the Supplementary Material.) (a) Temperature (b) C 1Density (c) C 2Density (d) C 3Density Figure 4: Trajectory of the particle in an argon ow eld at 0.15 ms with heterogeneous carbon reactions. (An animation of this gure is provided in the Supplementary Material.) vapor dominates other species. The simulation results for the mass and energy coupling and the total coupling are virtually identical to 150 13 (a) Temperature (b) C 1density (c) C 2density (d) C 3density Figure 5: Trajectory of the particle in an argon ow eld at 0.20 ms with heterogeneous carbon reactions. (An animation of this gure is provided in the Supplementary Material.) (a) Temperature (b) C 1density (c) C 2density (d) C 3density Figure 6: Trajectory of the particle in an argon ow eld at 0.26 ms with heterogeneous carbon reactions. (An animation of this gure is provided in the Supplementary Material.) the ones obtained for mass coupling. The temperature varies by 0 :01 K and no change is detected for the velocities when spallation source terms are coupled. 14 The ratio of the source term vector ( SpV) and residual vector ( R), as de ned in Eqs. (11) and (4), respectively, is called the coupling term vector [25]. The absolute magnitudes of these elements determine the impact of the two-way coupling. A low absolute magnitude of these elements indicates that the corresponding 155 coupling term does not a ect the uid. The presence of carbon species in the ow eld is only due to the particle's surface reactions. The mass uxes and the production rates of these species from the uid side are negligibly small, thereby making the mass coupling terms high in magnitude. Hence, it can be seen from the above simulation that the carbon species produced by the particle a ects the ow eld by changing its composition. On the other hand, the momentum and energy uxes of the uid are higher when compared 160 to the drag force and energy rate of the particle, making the momentum and energy coupling irrelevant. Similar behavior is also observed in the cases presented later in this paper. However, since the spallation phenomenon involves the ejection of thousands of particles whose combined source term is comparable to the ow eld uxes, a noticeable change to the temperature and dynamics of the gas mixture might become more apparent. 165 3.1.2. Air ow eld The mass removal of a spalled particle in an air environment occurs as a result of oxidation, nitridation and sublimation. The ejection parameters and physical time of the particle considered for this test case are shown in Table 6. Table 6: Ejection parameters and physical time for single particle simulation in air environment Size ( m) Velocity (m/s) Position (mm) Angle () Physical time (ms) 20 120 (0,10,0) 0 0.4580 To achieve time accuracy, the solution was computed with a time step size of 5 1010s, which corresponds 170 to a maximum CFL of 0.625. Figures 7, 8, 9, and 10 represent the simulation results of the mass coupling for the particle given in Table 6. These gures consist of temperature pro les and density pro les of CO, CN, C 1, C2, and C 3at di erent particle travel times. As opposed to the argon ow eld, which only accounted for sublimation, the air ow eld also includes oxidation and nitridation which depend on the concentrations of reaction species. 175 The animation of this simulation is provided in the Supplementary Material. Due to a high concentration of atomic oxygen and atomic nitrogen near the surface of the sample, the particle reacts and releases vapor soon after its ejection. However, the particle sublimates insigni cantly due to its relatively low temperature. This behavior is shown in Fig. 7 at 0 :10 ms of its travel time. As the particle crosses the shock, at 0 :20 ms of its travel, it starts to sublimate and release vapor of 180 varying concentrations. The CO vapor's concentration decreases in the upstream region as atomic oxygen 15 (a) Temperature (b) CO density (c) CN density (d) C 1density (e) C 2density (f) C 3density Figure 7: Trajectory of the particle in a reacting air ow eld at 0.10 ms with heterogeneous carbon reactions. (An animation of this gure is provided in the Supplementary Material.) is nearly absent, and oxidation occurs only due to molecular oxygen. Also, the nitridation is minimal in the post-shock region due to nearly zero concentration of atomic nitrogen. Therefore, only a minuscule amount of CN vapor is produced as seen in Fig. 8. 16 (a) Temperature (b) CO density (c) CN density (d) C 1density (e) C 2density (f) C 3density Figure 8: Trajectory of the particle in a reacting air ow eld at 0.20 ms with heterogeneous carbon reactions. (An animation of this gure is provided in the Supplementary Material.) At 0:30 ms, the particle is still in the upstream region. As its temperature increases, the concentration 185 of sublimate vapor increases. However, the concentrations of CO and CN vapor remain almost constant in this region. This behavior can be observed in Fig. 9. 17 (a) Temperature (b) CO density (c) CN density (d) C 1density (e) C 2density (f) C 3density Figure 9: Trajectory of the particle in a reacting air ow eld at 0.30 ms with heterogeneous carbon reactions. (An animation of this gure is provided in the Supplementary Material.) Figure 10 shows the particle after 0 :44 ms of travel time. At that time, the particle is in the downstream region, and the amount of vapor released by oxidation, nitridation, and sublimation reactions increases. From these gures (Figs. 7-10), it can be seen that the dominating species is CO, which is directly 18 (a) Temperature (b) CO density (c) CN density (d) C 1density (e) C 2density (f) C 3density Figure 10: Trajectory of the particle in a reacting air ow eld at 0.44 ms with heterogeneous carbon reactions. (An animation of this gure is provided in the Supplementary Material.) released as a result of oxidation. However, the spectroscopic experiments [7] observed high CN spectra in the upstream region. Since nitridation does not generate the amount of CN necessary to produce the observed quantities, other chemical processes must be present. It is hypothesized that the presence of CN is 19 due to the transformation of CO through a gas-phase exchange reaction [26, 27]: CO + N*)CN + O (18) In order to evaluate the conversion of CO to CN, this single equation is added to the ow eld chemistry 190 model of the CFD code, and the exact same simulation is performed again. Figure 11 compares the solutions of CN density pro le when the exchange reaction given by Eq. 18 is considered or not. (a) CN density at 0.20 ms (without exchange reaction) (b) CN density at 0.20 ms (with exchange reaction) (c) CN density at 0.30 ms (without exchange reaction) (d) CN density at 0.30 ms (with exchange reaction) Figure 11: CN density pro les without and with the exchange reaction at 0.20 and 0.30 ms Figures 11 ( a) and 11 (b) represent the simulation of the particle at 0 :20 ms of its travel time while Figs. 11 (c) and 11 (d) depict the simulation at 0 :30 ms. Figures 11 ( a) and (c) represent the CN density pro le without adding the exchange reaction, which shows that the concentration of CN vapor released and 195 subsequently di used is very small. Figures 11 ( b) and (d) illustrate the CN density pro le when the exchange reaction is included. It is noted that there is a considerable amount of CN in the downstream region and very little in the upstream region. The CO vapor released by the particle within the shock is converted to CN. However, due to a very small concentration of N in the post-shock region, the exchange reaction is insigni cant, and a very low concentration of CN is observed. It is, however, postulated that if more reactions 200 are considered, more CN would be observed in the upstream region. Therefore, a list of dissociation and exchange gas-phase reactions [28] is added to the chemistry model of CFD, and listed in Table 7. It should be noted that the rst three of the dissociation reactions and the rst two of the exchange reactions were 20 already included in the earlier simulation (Figs. 7-10). The reactions involving carbon products are the ones that are added for this simulation. The simulation with same ejection parameters as given in Table 6 are 205 used to understand the e ects of additional chemistry, and the results are shown in Figures 12, 13, 14, and 15. The animation of this simulation is provided in the Supplementary Material. Table 7: Dissociation and exchange reactions considered for single particle simulation in air environment Type Reactions Dissociation Reactions8 >>>>>>>>>>>>>>>< >>>>>>>>>>>>>>>:N2+ M*)2 N + M O2+ M*)2 O + M NO + M*)N + O + M C2+ M*)C1+ C 1+ M CN + M*)C1+ N + M CO + M*)O + C 1+ M Exchange Reactions8 >>>>>>>>>>>>>>>>>>>>>>>>>>>>>>< >>>>>>>>>>>>>>>>>>>>>>>>>>>>>>:N2+ O*)NO + N NO + O*)O2+ N CO + N*)CN + O C1+ N 2*)CN + N C1+ NO*)CN + O C1+ NO*)CO + N C1+ O 2*)O + CO C2+ C 2*)C3+ C 1 C2+ N 2*)CN + CN CN + C 1*)C2+ N At 0:10 ms of the particle's travel time, the added reactions increase the concentration of C 1as seen in Fig. 12 (d) in spite of the temperature being low. Since half of the dissociation reactions produce C 1, it is believed that these may be responsible for an increase in C 1concentration. Figure 12 also shows the other 210 products' density pro les as a result of particle's oxidation, nitridation, and sublimation (C 2and C 3) which 21 closely resemble the earlier simulation shown in Fig. 7. (a) Temperature (b) CO density (c) CN density (d) C 1density (e) C 2density (f) C 3density Figure 12: Trajectory of the particle in a reacting air ow eld at 0.10 ms with heterogeneous and homogeneous carbon reactions. (An animation of this gure is provided in the Supplementary Material.) When the particle crosses the shock at 0 :20 ms, a change in the C 1and CN concentrations is observed in 22 both upstream and downstream locations. However, the change in CN density pro le is similar to what was observed in Fig. 11 (b). While comparing with Fig. 8, the only signi cant dissociation reactions are the ones 215 in which C 2and CO are converted to C 1vapor. The sublimate vapor further undergoes exchange reactions to produce more of CN and C 1which can be seen in Fig. 13. The other concentrations remain unchanged with minimal e ect from added chemical reactions. At 0:30 ms of particle's travel time, there is a further increase in the C 1and CN concentration in both upstream and downstream of the shock. It should be noted that when compared to the CN concentration 220 as seen in Fig. 11 (d), the added reactions produce more CN. It can be inferred that the presence of CN in the upstream region of the shock, as observed from the spectroscopic experimental data [7], is indeed due to the additional gaseous phase reactions between carbon products and air. The oxidation and sublimation of the particle undergo similar changes when compared with Fig. 9. The above-mentioned behavior can be observed in Fig. 14. 225 At 0:44 ms of its travel as shown in Fig. 15, the particle enters the downstream region with higher CO, C2, and C 3concentrations, and lower CN and C 1concentrations when compared to Fig. 10. It is believed that the exchange reactions involving C 1and CN play a signi cant role in this region. These results indicate that along the particle's path, their chemistry in conjunction with the gas-phase reactions increases the concentration of carbon products which could not be explained by other ablative 230 phenomena. Similar to the argon environment, when the mass and energy coupling and the total coupling results are compared to only mass coupling simulation results, they appear virtually identical: the temperature is di erent by nearly 0 :01 K when spallation source terms are applied. From these single-particle simulation results, it can be seen that the particle releases vapor in both 235 upstream and downstream regions. This reinforces the hypothesis that the presence of spalled particles ahead of the shock is the reason for the optical emissions of the carbon species observed in that region. The conclusions in this section are based on the simulation of one particular trajectory of one particle. The trajectory of the particle is a function of its size and ejection parameters. In order to assess the e ect of the particle on the ow eld, a parametric study on the parameters that a ect the particle's motion is 240 conducted. 3.1.3. Parametric study The parametric study is performed by varying each parameter while keeping others constant. The objective of this study is to summarize the e ect of the particle's behavior on the ow eld by including all parameters that a ect its motion. The study is conducted in a Mach 5 air environment where the total 245 density of the carbon products released by the particle is considered. A baseline case for simulation is used and other simulations where each parameter is changed are compared to the baseline case. The results 23 (a) Temperature (b) CO density (c) CN density (d) C 1density (e) C 2density (f) C 3density Figure 13: Trajectory of the particle in a reacting air ow eld at 0.20 ms with heterogeneous and homogeneous carbon reactions. (An animation of this gure is provided in the Supplementary Material.) are compared at the same timestamps since the simulation is transient in nature. The initial parameters and total physical time for simulation considered as the baseline case for the parametric study are given in 24 (a) Temperature (b) CO density (c) CN density (d) C 1density (e) C 2density (f) C 3density Figure 14: Trajectory of the particle in a reacting air ow eld at 0.30 ms with heterogeneous and homogeneous carbon reactions. (An animation of this gure is provided in the Supplementary Material.) Table 8 and the simulation results are shown in Fig. 16. All the simulations considered for the parametric 250 study were computed with a time step size of 8 1010s, which corresponds to the maximum CFL of 1.00. 25 (a) Temperature (b) CO density (c) CN density (d) C 1density (e) C 2density (f) C 3density Figure 15: Trajectory of the particle in a reacting air ow eld at 0.44 ms with heterogeneous and homogeneous carbon reactions. (An animation of this gure is provided in the Supplementary Material.) The simulation of the baseline case shows the particle's path at 0 :10 ms, as it travels towards the shock and releases carbon vapor soon after the ejection. As the particle crosses the shock at 0 :20 ms, there is a 26 Table 8: Ejection parameters and physical time of baseline case for parametric study Size ( m) Velocity (m/s) Position (mm) Angle () Physical time (ms) 24 100 (0,10,0) 0 0.5651 decrease in C density pro le. The concentration of C is increased again as the particle enters the shock area as observed at 0 :44 ms. The complete trajectory of the particle is presented over the temperature contour 255 of the ow eld in Fig. 16 (a). The chemical behavior of the particle is similar to the one discussed in Section 3.1.2. Variation of particle size. A time-accurate solution is computed by varying only particle size while other parameters are kept constant. The ejection parameters and total physical time for the simulation are given in Table 9. 260 Table 9: Ejection parameters and physical time for parametric study with variation of size Size ( m) Velocity (m/s) Position (mm) Angle () Physical time (ms) 20 100 (0,10,0) 0 0.4624 The result of the above simulation is shown in Fig. 17. It should be noted that due to the decrease in size and thus reduced inertia, the particle interacts with the shock but does not cross it. Compared to Fig. 16, there is an overall increase in the released carbon vapor by the particle. The particle travels through temperature gradients in the downstream region of the shock which causes an increase in its temperature and more surface recession. However, the particle's total physical time to pass through the computational 265 domain is less when compared to the baseline case. Therefore, comparing at 0 :44 ms is not possible. Instead, the simulation at 0 :32 ms is presented in Fig. 17 (e), before it exits the computational domain. Variation of ejection velocity. The ejection velocity is changed and a time-accurate solution is computed. The ejection parameters and total physical time for the simulation are given in Table 10 and the simulation results in Fig. 18. 270 Table 10: Ejection parameters and physical time for parametric study with variation of ejection velocity Size ( m) Velocity (m/s) Position (mm) Angle () Physical time (ms) 24 75 (0,10,0) 0 0.3075 The results show a similar behavior to the study conducted by changing the size of the particle. With low ejection velocity, the particle travels through the downstream region of the shock. The path through 27 (a) Temperature (b) C density at 0 :10 ms (c) C density at 0 :20 ms (d) C density at 0 :30 ms (e) C density at 0 :44 ms Figure 16: Trajectory of the particle in an air ow eld for baseline case high-temperature regions allows the particle to react more and release more carbon vapor which can be observed in Fig. 18. The particle's behavior is not compared at 0 :44 ms as it exits the computational domain well before that time. 275 28 (a) Temperature (b) C density at 0 :10 ms (c) C density at 0 :20 ms (d) C density at 0 :30 ms (e) C density at 0 :32 ms Figure 17: Trajectory of the particle in an air ow eld for 20 m particle (other ejection parameters same as baseline case). Variation of ejection position. The ejection parameters and total physical time for the time accurate simu- lation where the ejection position is varied are presented in Table 11. The results of the simulation are illustrated in Fig. 19. The trajectory of the particle crosses the shock 29 (a) Temperature (b) C density at 0 :10 ms (c) C density at 0 :20 ms (d) C density at 0 :30 ms Figure 18: Trajectory of the particle in an air ow eld for ejection velocity of 75 m/s (other ejection parameters same as baseline case) Table 11: Ejection parameters and physical time for parametric study with variation of ejection position Size ( m) Velocity (m/s) Position (mm) Angle () Physical time (ms) 24 100 (0,5,0) 0 0.6137 for this case. In comparison with Fig. 16, it can be seen that the particle releases a large amount of carbon with di used vapor also traveling towards the center-line direction. The temperature of the uid and 280 concentrations of atomic oxygen and nitrogen increase as the position moves towards the centerline. The particle's chemical reactivity increases when it is ejected from the positions that are closer to the centerline. Also, with the bulk velocity of the uid being low at this position, the released carbon vapor is di used along the downstream region of the shock and around the sample. Variation of ejection angle. For this case, three di erent time-accurate solutions with di erent ejection 285 angles are considered and their details are given in Table 12. The simulation results for ejection angle 30, 60, and 90are presented in Fig. 20, 21, and 22, respectively. The particle with an ejection angle of 30passes through the shock and travels through the upstream region until it exits the computational domain. Based on the trajectory and in comparison to Fig. 16, it can be seen that the density pro le of C released by the particle increases as it travels towards the shock and goes 290 30 (a) Temperature (b) C density at 0 :10 ms (c) C density at 0 :20 ms (d) C density at 0 :30 ms (e) C density at 0 :44 ms Figure 19: Trajectory of the particle in an air ow eld for ejection position at (0 ;5;0) (other ejection parameters same as baseline case) on to decrease thereafter. This behavior of the particle can be seen in Fig. 20. The path of the particle ejected at an angle of 60is similar to that of a particle ejected at an angle of 30as seen in Fig. 21. However, the particle passes through lower temperature zones when compared to that 31 Table 12: Ejection parameters and physical time for parametric study with variation of ejection position Size ( m) Velocity (m/s) Position (mm) Angle () Physical time (ms) 24 100 (0,10,0) 30 0.4748 24 100 (0,10,0) 60 0.3303 24 100 (0,10,0) 90 0.2199 of the 30angle trajectory. The density pro le of C increases as the particle moves towards the shock and decreases as it travels through the upstream region. The comparison to 0 :44 ms travel time is not performed 295 as the particle exits the computational domain before that time. In the case of the ejection angle 90as shown in Fig. 22, the particle travels through the downstream region of the shock and exits the computational domain in a very short time. With the low bulk velocity of the uid near the sample, the released carbon vapor from the particle is di used along the downstream region and around the sample. The concentration of the released vapor keeps increasing until it exits the 300 domain. The parametric study conducted indicates that the closer the ejection position of the particle to the center-line, the greater the chemical reactivity of the particle and hence, larger amount of vapor released. With the decrease in mass and ejection velocity, the particle tends to travel through the downstream region and does not cross the shock, thereby, releasing more vapor. In addition, the particle's reactivity increases 305 and the particle releases more vapor with decrease in ejection angle from 90to 0. 3.2. Multiple particles simulation In this section, the same process described earlier is repeated, but more than one particle is ejected in the ow eld at the same time. 3.2.1. Argon ow eld 310 The ejection parameters and physical time of the particles considered for multiple particle simulations are given in Table 13. Table 13: Ejection parameters and physical time for multiple particles simulation in argon environment Size ( m) Velocity (m/s) Position (mm) Angle () Physical time (ms) 12 380 (0,10,0) 0 0.2370 14 370 (0,13,0) 0 0.2697 16 360 (0,07,0) 0 0.3323 32 (a) Temperature (b) C density at 0 :10 ms (c) C density at 0 :20 ms (d) C density at 0 :30 ms (e) C density at 0 :44 ms Figure 20: Trajectory of the particle in an air ow eld for 30ejection angle (other ejection parameters same as baseline case) Figures 23, 24, 25, and 26 illustrate the multiple particles simulation in an argon ow eld at di erent travel times. The simulation was performed at a time step size of 4 109s, corresponding to a maximum CFL of 0.824, which ensures time accuracy. The animation of this simulation is provided in the Supplementary 315 33 (a) Temperature (b) C density at 0 :10 ms (c) C density at 0 :20 ms (d) C density at 0 :30 ms Figure 21: Trajectory of the particle in an air ow eld for 60ejection angle (other ejection parameters same as baseline case) (a) Temperature (b) C density at 0 :10 ms (c) C density at 0 :20 ms Figure 22: Trajectory of the particle in an air ow eld for 90ejection angle (other ejection parameters same as baseline case) 34 Material. The behavior of the particles is similar to the results observed for the single particle simulation. At 70 s from ejection, the particles approach the shock. With their temperatures being relatively low, they show negligible signs of recession, as seen in Fig. 23. (a) Temperature (b) C 1density (c) C 2density (d) C 3density Figure 23: Trajectories of the particles in an argon ow eld at 0.07 ms with heterogeneous carbon reactions. (An animation of this gure is provided in the Supplementary Material.) The particles tend to sublimate in the post-shock region (0 :15 ms of travel time), releasing di erent vapor 320 concentrations. The vapor then convects and di uses in a much larger region downstream, as seen in Fig. 24. At a travel time of 0 :23 ms, two of the particles re-enter the shock region as shown in Fig. 25. The concentrations of the carbon vapor released increase signi cantly and propagate to a much larger area. As soon as these two particles complete their travel in the computational domain, the di used vapor area decreases and the contribution only comes from the third particle, as shown in Fig. 26. 325 It is noticed that even with multiple spalled particles in the ow eld, the temperature changes by a maximum of only 0 :1 K. The other parameters remain the same throughout the simulation. 3.2.2. Air ow eld The ejection parameters for the particles simulated in this section are given in Table 14. The solution is computed at a time step size of 8 1010s, which corresponds to a maximum CFL of 1 :000 to maintain 330 time accuracy. The particles of size 20 m, 25 m, and 30 m penetrate through the shock, whereas the 15 m does 35 (a) Temperature (b) C 1density (c) C 2density (d) C 3density Figure 24: Trajectories of the particles in an argon ow eld at 0.15 ms with heterogeneous carbon reactions. (An animation of this gure is provided in the Supplementary Material.) (a) Temperature (b) C 1density (c) C 2density (d) C 3density Figure 25: Trajectories of the particles in an argon ow eld at 0.23 ms with heterogeneous carbon reactions. (An animation of this gure is provided in the Supplementary Material.) not. The behaviors of these particles are shown in Figs. 27, 28, 29, and 30 at the respective times. The animation of this simulation is provided in the Supplementary Material. 36 (a) Temperature (b) C 1density (c) C 2density (d) C 3density Figure 26: Trajectories of the particles in an argon ow eld at 0.33 ms with heterogeneous carbon reactions. (An animation of this gure is provided in the Supplementary Material.) Table 14: Ejection parameters and physical time for multiple particles simulation in air environment Size ( m) Velocity (m/s) Position (mm) Angle () Physical time (ms) 15 70 (0,05,0) 0 0.1866 20 120 (0,10,0) 0 0.4580 25 100 (0,13,0) 0 0.5590 30 85 (0,08,0) 0 0.5744 At 0:09 ms, the particles start producing large concentrations of CO and CN vapors. However, the 335 particles hardly sublimate as their temperatures are relatively low. The simulation at this time is illustrated in Fig. 27. The released vapor is seen to di use over a larger area and also observed to di use towards the center-line direction. Figure 28 shows the simulation at 0 :19 ms. At this time, the 20 m particle passes through the shock, the 25 m particle nearly comes close to the shock, and the 30 m particle is still traveling towards the 340 shock. The 15 m particle, however, has already changed direction while in the high temperature region and will never reach the shock. The contribution to the CO vapor density decreases for particles that cross the shock and remains almost the same for the particles that are in the downstream region. Similarly, the particles that cross the shock generate very little CN vapor while the other particles contribute moderately. The 15 m particle starts sublimating, and since it is located in a low-velocity zone, the di used vapor tends 345 37 (a) Temperature (b) CO density (c) CN density (d) C 1density (e) C 2density (f) C 3density Figure 27: Trajectories of the particles in a reacting air ow eld at 0.09 ms with heterogeneous carbon reactions. (An animation of this gure is provided in the Supplementary Material.) to move towards the centerline. The 20 m and 25 m particles also sublimate, contributing a very small amount of C 1, C2, and C 3vapors. At 0:35 ms, 20 m, 25 m, and 30 m particles are in the post-shock region, whereas the 15 m particle 38 (a) Temperature (b) CO density (c) CN density (d) C 1density (e) C 2density (f) C 3density Figure 28: Trajectories of the particles in a reacting air ow eld at 0.19 ms with heterogeneous carbon reactions. (An animation of this gure is provided in the Supplementary Material.) completes its travel, moving out of the computational domain. The CO vapor concentration decreases, and CN vapor becomes negligible. This behavior can be seen in Fig. 29. It is also observed that a small 350 concentration of CO and CN remains close to the surface of the ablator. 39 (a) Temperature (b) CO density (c) CN density (d) C 1density (e) C 2density (f) C 3density Figure 29: Trajectories of the particles in a reacting air ow eld at 0.35 ms with heterogeneous carbon reactions. (An animation of this gure is provided in the Supplementary Material.) At 0:55 ms time, 20 m particle exits the domain, and the other two particles (25 m and 30 m) enter the downstream region. The CO, CN, C 1, C2, and C 3concentrations seem to increase in magnitude, as illustrated in Fig. 30. However, a small amount of concentration of CO and CN remains close to the ablator. 40 (a) Temperature (b) CO density (c) CN density (d) C 1density (e) C 2density (f) C 3density Figure 30: Trajectories of the particles in a reacting air ow eld at 0.55 ms with heterogeneous carbon reactions. (An animation of this gure is provided in the Supplementary Material.) Apart from the e ect on the composition of the ow eld, the multiple particles simulation has a very 355 small e ect on the temperature and velocity elds. It is observed that the translational temperature changes by 0:3 K and vibrational temperature changes by 0 :05 K. Also, the velocity eld in the y-direction shows 41 regions of change in velocity in the post-shock region. Figure 31 illustrates these changes in the upstream region at 0 :35 ms. The positions of the particles are shown in Fig.29 ( a). It is noticed that there is a di erence of 0 :03 m/s of velocity in the y-direction. 360 Figure 31: Flow eld Velocity in y-direction for the air simulation at 0.35 ms The results from the multiple particle simulations illustrate that the products released from the particle convect and di use over a larger area. Although the number of particles in these simulations were three and four, the actual number in the spallation phenomena will be in the thousands. Therefore, the concentration of carbon products in total would be signi cant. The presence of this carbon vapor tends to increase the radiative heat ux directed toward the sample. Since these particles release vapor as they travel, the shift in 365 the radiative heat ux is transient in nature. This also explains why the heat rates have been under-predicted in the downstream region in past studies. Therefore, evaluating the radiative heat ux contribution from these spalled particles is an essential exercise that will be performed in future work. 42 4. Conclusion A two-way coupling was achieved using a Lagrangian particle trajectory code and a CFD code. The 370 results obtained through this coupling allowed assessment of the spalled particle e ect on the ow eld. Initially, the coupling was performed for a single particle spalled in Mach 5 argon and air environments. The results in the argon ow showed that the particle starts releasing the sublimate vapor once it crosses the shock, and its concentration keeps increasing as it moves through the computational domain. It was also observed that the magnitude of concentration of sublimate vapor increases from C 1species to C 3species 375 at any point of the particle's travel. In an air environment, the particle reacts and releases the vapor soon after the ejection. It was observed that the particle releases CO vapor throughout its entire travel whereas CN vapor is released only in small quantities ahead of the shock. The sublimation of the particle follows the same trend as seen in the argon ow eld. The magnitude of concentration of CO vapor released was higher than that of the other vapors. Additional gas-phase reactions were added to the CFD code and the same 380 simulation was run. The results including added chemistry showed an increase in CN and C 1vapor along with CO in the air environment. However, the coupling results show a very small change in temperature of about 0:01 K and no di erence in the ow eld velocity. To summarize the impact of the particle on the ow eld, a parametric study was conducted. The study concluded that the closer the ejection position to the centerline and the lower the size and the ejection velocity, the more mass that is shed through chemical 385 reactions. Similar behavior was also observed with the ejection angle of the particle decreasing from 90to 0. The coupled simulation was also performed for multiple particles spalled in argon and air environments. The results for an argon environment showed a pattern similar to that of the single-particle simulation. The simulation for three particles showed that the vapor released was convected and di used over a larger 390 area and resulted in the change of temperature of 0 :1 K. Similarly, for an air environment, the simulation was performed using four particles. The CO and CN vapor released di used over a larger area, and a small concentration of these species remained close to the surface of the ablator. The particles also tend to sublimate in the post-shock region, and the concentration of the sublimate species increased as the particles moved through the computational domain. The particles that did not reach the shock sublimated near the 395 surface. These coupling results indicate a likely relationship between the presence of spalled particles and optical emissions, observed ahead of the shock. Also, the di usive uxes encompassed over larger areas provide a plausible explanation for the under-prediction of heat rates and temperature pro les in the downstream region. 400 43 5. Acknowledgements Financial support for this work was provided by NASA Kentucky EPSCoR Award NNX10AV39A, and NASA award NNX13AN04A. The authors would like to thank L. P. Askins (University of Kentucky), P. Ghosh (VIT-AP), and D. A. Saunders (AMA, Inc.) for reviewing the manuscript.