A blended drag coecient model is constructed using a series of empirical relations based on Reynolds number, Mach number, and Knudsen number. When validated against experiments, the drag coecient model produces matching values with a standard deviation error of 2.84% and a maximum error of 11.87%. The model is used in a Lagrangian code which is coupled to a hypersonic aerothermodynamic CFD code, and the particle velocity and trajectory are validated against experimental results. The comparative results agree well and show that the new blended drag coecient model is capable of predicting the particle motion accurately over a range of Reynolds number, Mach number, and Knudsen number. Keywords: Lagrangian Particle Trajectory, Particle-laden ows, Drag coecient, Hypersonics, Aerothermodynamics Nomenclature 1 Symbols 2 a Speed of sound, m/s A Area of cross-section, m2 C Compressibility factor CD Drag coecient d Diameter, m E Internal energy, J/kg F Force, N F Force vector, N FcFcFc Convective ux matrix FdFdFd Di usive ux matrix Corresponding author Email addresses: raghava.davuluri@uky.edu (Raghava S. C. Davuluri), sean.bailey@uky.edu (Sean C. C. Bailey), kaveh.tagavi@uky.edu (Kaveh A. Tagavi), alexandre.martin@uky.edu (Alexandre Martin) Preprint submitted to International Journal of Heat and Fluid Flow September 24, 2020 k Slip coecient Kn Knudsen number L Characteristic length, m M Mach number _m Rate of change of mass due to chemical reactions, kg/s m Mass, kg pforce Power required to overcome force, W _q Heat rate, W Q Conservative variables vector r Radius, m Re Reynolds number S Molecular speed ratio (= Mp 2) _S Source term vector for KATS T Temperature, K (u;v;w ) Velocity components in axial, radial and z-directions, m/s U State vector V Velocity vector, m/s W Source vector Yi Mass fraction, kg/kg (x;y;z ) Position components in axial, radial and z-directions, m Ratio of speci c heats  Error 0Correction term for k  Molecular mean free path, m  Density, kg/m3  Rarefaction parameter Subscripts 3 conv Convection exp Experimental f Surrounding uid num Numerical p Particle r Relative 2 rad Radiation rxn Chemical reactions u Uncertainty w Wall 1 Free-stream 1. Introduction 4 Particle-laden ows are an essential class of multiphase ows in which small, discrete particles are im- 5 mersed within a carrier uid. Several notable examples can be found in nature or technology (e.g., rocket 6 nozzle ows [1, 2]; resuspension phenomena [3]; dusty-laden environmental ows [4, 5]; sprays [6]; granular 7 ows [7]; powder ows [8]; and slurry ows [9]). Seeding particles are also intentionally introduced into 8 uids to enable measurement of the uid velocity (e.g., Particle Image Velocimetry - PIV [10, 11]; Particle 9 Tracking Velocimetry - PTV [12]; and Laser Doppler Anemometry - LDA [13, 14]). 10 The primary approaches used to model particle-laden ows are Eulerian-Eulerian, Eulerian-Lagrangian, 11 and fully Eulerian with uid-structure interaction (FSI) methods. In the Eulerian-Eulerian method, the 12 carrier uid, as well as discrete particles, is treated as continuous and interpenetrating continua. Both the 13 particle phase and the carrier uid phase are solved using the Navier-Stokes equations. The set of equations 14 for the two-phase ow is solved in the Eulerian frame of reference, and the method is primarily used when the 15 volume fraction of the particles is high. In this methodology, various bins of particles, each with a speci c 16 size, are assumed, and additional relationships are used to account for particle-wall and particle-particle 17 interactions. 18 A fully Eulerian with FSI method is used in the cases where discrete particles undergo large structural 19 deformation and phase change in carrier uid ow. There are various approaches to model the particle 20 ows using this methodology [15, 16, 17, 18]. Based on the approach, carrier uid employs Eulerian or 21 arbitrary Lagrangian-Eulerian formulations, whereas Lagrangian formulation is used for discrete particles. 22 This method can be used for particles that vary in size. 23 In the case of the Eulerian-Lagrangian method, the carrier uid is treated as a continuum, and governing 24 Navier-Stokes equations are solved in the Eulerian frame of reference. The interactions of the discrete 25 particles with the carrier uid are then modeled in the Lagrangian frame of reference, a posteriori using 26 information from the solution for the carrier uid. The method can be employed for particles that vary in 27 size when interactions between the two phases are essential to consider. 28 The particles interact with the carrier uid through mass, momentum, and energy. The mass interaction 29 between the particles and the uid takes place in the form of chemical reactions. The chemical products 30 3 from the reactions change the composition of the uid and size (and shape) of the particles. On the other 31 hand, the momentum interactions take place in the form of force exerted by the uid and reaction force from 32 the particles. The uid dictates the motion of the particle through the exertion of forces. The various types 33 of forces exerted on the particle by the uid are discussed in Section 2.3 brie y. The reaction force from 34 the discrete particles may change the velocity of the uid. Lastly, the work done by the exerted forces and 35 the heat transfer between the uid and particle constitutes the energy interactions. The energy interactions 36 results in a change in the temperature of discrete particles and the carrier uid. 37 Accurately modeling particle behavior can be essential in the study of the thermal protection system 38 (TPS) behavior during atmospheric entry. The brous ablative materials that are used as a TPS in many 39 space vehicles counter the high heat rates and have been observed to expel particles in a process referred to 40 as \spallation" [19, 20, 21, 22, 23, 24]. Experiments have found that this phenomenon is a non-negligible 41 mechanism of surface degradation [25], and particles shed from the surface could introduce additional chem- 42 ical species outside of the surface boundary layer [26]. Modeling the particle trajectories in these types of 43 ows, however, is complicated by the potential for individual particles to experience ow regimes ranging 44 from subsonic to hypersonic, within gases which can be rare ed or continuum. 45 Here, we present and validate a Lagrangian particle trajectory prediction and interaction model im- 46 plemented within a computational uid dynamics framework following the Euler-Lagrangian method. To 47 facilitate the use of this model in the study of the spallation of TPS materials [27], the model was designed 48 to predict particle trajectories for the full range of subsonic to hypersonic ow regimes and a wide range of 49 Knudsen numbers. The objective of this study is to demonstrate that the numerical model can predict the 50 particle motion. The objective is achieved in two steps. First, a blended drag coecient model is constructed 51 using a series of empirical relations based on the ow and particle conditions. Validation is performed by 52 comparing the results of the model with experiments, ensuring that the model is valid over various regimes. 53 Then, the model is used in a Lagrangian particle code and coupled to a hypersonic aerothermodynamic CFD 54 code so that particle properties (velocity and trajectory) can be validated against experimental results. 55 2. Numerical Approach 56 2.1. Particle motion 57 The particle model simulates the dynamics of a particle by employing a Lagrangian formulation [28] and includes the chemical interaction of the particle with the ow eld [27]. The model assumes the particle to be spherical in shape. In order to build the model, the governing equations are cast in the form of @U @t=W (1) 4 where Uis the state vector, and Wis the source term vector. The elements of the vectors are then represented as U=0 BBBBBBBBB@mp mpup mpvp mpwp mpEp1 CCCCCCCCCA; W=0 BBBBBBBBB@_mchemistry Fpx Fpy Fpz _qconv+pforce_qrad+ _qrxn1 CCCCCCCCCA(2) wherempis the mass of the particle, with ( up;vp;wp) its Cartesian velocity components, and Epits total 58 energy. The particle mass can be altered through the chemical reaction source term _ mchemistry ; the velocity 59 through the Cartesian components of total force acting on the particle Fpx;Fpy;Fpz ; and the total energy 60 through the convective ( _ qconv), radiative ( _ qrad), and reaction ( _ qrxn) heat rates, as well as through the power 61 required to overcome force acting on the particle ( pforce). The chemistry model accounts for oxidation, 62 nitridation, and sublimation of carbon particle. Though the chemistry model was previously used on carbon 63 particles, similar models can be employed for particles of other compositions. The discretization of the system 64 of equations in Eq. 2 is performed using the backward Euler method. The mass conservation equation is 65 uncoupled from momentum and energy conservation equations, and the block Gauss-Seidel method is used 66 to solve the two sets of equations. The set of momentum and energy equations are solved using Newton's 67 method. This approach has been veri ed using the method of manufactured solution [29], which con rmed 68 the order of accuracy of the discretization and correctness of the numerical code. 69 2.2. Flow eld 70 The ow elds used in this study are determined using converged solutions of the aerothermodynamic computational uid dynamics (CFD) code KATS. KATS is a laminar aerothermodynamics Navier-Stokes solver that uses a nite-volume approach with various modules [30, 31, 32, 33] for simulating complex phenomena in ablation-related problems, including spallation studies [27, 28, 26, 24]. The module KATS- CFD is used for computing compressible viscous ows including thermo-chemical nonequilibrium ows by solving the corresponding governing equations in the form @Q @t+r(FcFcFcFdFdFd) =_S; (3) where Qis a vector of conservative variables, FcFcFcandFdFdFdare convective and di usive ux matrices, and _S 71 is a source term vector. For solution, the system of equations represented in Eq. 3 is discretized rst-order 72 in time and second-order in space. The solver uses the PETSc library [34, 35, 36] to solve the linear system 73 of equations, with ParMETIS [37] used for domain decomposition, and MPI [38] for message passing. More 74 information regarding the solver is given in Ref. [39]. The coupling of the Lagrangian particle code and the 75 5 hypersonic aerothermodynamic CFD code was developed in previous works, and used to study the e ects of 76 the ow eld on the particles, and vice-versa [26, 27, 28]. 77 2.3. Drag coecient 78 To obtain Fp= (Fpx;Fpy;Fpz), we note that the rate of change of momentum of a particle can be described by the Basset-Boussinesq-Oseen (BBO) equation [40] which de nes Fp=Fpdrag+Fppressure +Fpadd+Fpbassett +Fext; (4) where Fpdragis the steady-state drag force acting on the particle that assumes no change of relative velocity 79 between the particle and carrier uid, and a uniform pressure eld. Fpdragis equal to the sum of pressure 80 and shear forces acting on the particle surface. Fppressure is the force due to the local pressure gradient that 81 accelerates the particle { also referred to as the buoyancy force { and is equivalent to the weight of the 82 uid displaced that would otherwise occupy the volume of the particle. This force plays a signi cant role 83 in slurry ows and bubbly ows [40]. Fpaddis an unsteady force term called the added mass force . When a 84 particle accelerates, it displaces uid around it. Therefore, the force acting on the particle depends on an 85 e ective mass, which is equivalent to the sum of particle mass and added mass that comes from the change 86 in inertia in the uid. For spheres, the added mass is equal to half of the mass of the uid displaced by 87 the particle [41]. Fpbassett is also an unsteady force term and is called Basset force or Basset history force. 88 While Fpaddaccounts for pressure e ects on the particle due to acceleration, Fpbassett accounts for viscous 89 e ects [40]. When the particle accelerates through the uid, there is a change in its relative velocity with 90 respect to time. The change in relative velocity causes the boundary layers around the particle to change. 91 However, it takes some time for the boundary layers to respond to the change, and Fpbassett accounts for this 92 temporal delay. Lastly, Fextdenotes external forces such as gravity. 93 The drag force Fpdragis a function of uid and particle velocities, and is included in Eq. 4 for all conditions. 94 The added mass force Fpaddis present when a relative acceleration is detected between the phases. This 95 term becomes more signi cant when there is an increase in volume fraction and particle size [42]. The Basset 96 forceFpbassett only plays an important role in unsteady ows. This force increases when the volume fraction, 97 particle size, uid density, or uid viscosity increases [42, 43]. The magnitudes of FpaddandFpbassett are 98 directly dependent on ratio of uid to particle density [44]. 99 Several studies were performed [45, 46, 47, 48], which determined that for smaller particles (order of size: 0.1-100 m) and high ow velocities, all but the drag force were insigni cant. The drag force acting on the particle can be expressed as Fpdrag=1 2CDfApjVrjVr (5) so that the direction of Fpdragis inverse to that of Vr, the relative velocity of the particle to the surrounding 100 uid, and is de ned as ( VpVf), where Vprepresents the particle velocity vector and Vfis the velocity 101 6 vector of the surrounding uid. The drag coecient, as used in the above equation, includes all the physical 102 mechanisms that are responsible for the interphase force between the particles and the uid, including 103 pressure, viscous, and wave drags. The drag coecient is also a function of the shape of the body. However, 104 the model and the following discussion pertains only to spherical particles. The drag coecient depends on 105 various uid ow conditions, which can be summed up its dependence on three dimensionless parameters: 106 Reynolds number (Re), Mach number (M), and Knudsen number (Kn). 107 Reynolds number is de ned as the ratio of inertial to viscous forces. In the case of spheres, it is de ned as: Re =Inertial forces Viscous forces=fjVrjdp f(6) wherefis the density of the uid owing over the sphere, dpis the diameter of the sphere, and fis the 108 dynamic viscosity of the uid. For low Reynolds number (Re <1) where viscous forces are dominant and 109 inertial forces are negligible, CDis inversely proportional to Re and is within what is referred to as Stoke's 110 regime [49, 44]. In this regime { in the context of continuum { the ow remains attached to the particle 111 with no wake formation. As the Reynolds number increases (1 <Re<1000), the inertial forces become 112 e ective with ow separation, wake formation, and vortex shredding taking place. This regime is referred 113 to as the intermediate regime , and the inverse relationship between CDand Re continues until the former 114 reaches a constant value. In this regime, the sum of form drag and viscous drag contribute to CD. With 115 the further increase in Re, the inertial forces become dominant over viscous forces. As Re increases beyond 116 1000, the wake region becomes turbulent, whereas the boundary layer region in front remains laminar. This 117 ow regime is called Newton's regime , whereCDis only due to form drag and remains relatively constant 118 until the critical Re is reached. At the critical Re ( 3105), both the boundary layer and wake region 119 become turbulent. The separation point is further moved rearward, which causes a sharp decline in form 120 drag and, thus, a sudden decrease in CD. Following the decline, CDincreases with the increase in Re. The 121 dependency of CDon Re, as described above, indicate the inertial/viscous e ects . The surface roughness, 122 which also a ects the drag coecient, is not considered in this work, and the particles are assumed to be 123 smooth. 124 Mach number is de ned as the ratio between the relative velocity of the body and speed of sound and is de ned as: M =jVrj af(7) whereafis the speed of sound of the uid through which the body travels. The Mach number dependency 125 denotes the compressibility e ects on the drag coecient. It should be noted that the behavior of CD 126 with respect to Re, as explained above, is for incompressible and continuum ow. With compressibility 127 included,CDdoes not follow the same trend as described above. However, simpli cations assumed in 128 7 describing the behavior of CDare often justi ed. For low Mach numbers (M <0:5),CDis relatively 129 independent of M [40]. With the increase in M, CDincreases until the ow reaches the sonic condition 130 (M = 1) where there is a signi cant jump caused by the formation of the shock wave in front of the body. 131 The presence of shock wave produces a sudden rise in the forebody and afterbody pressure distribution that 132 a ects the drag force. However, this change in the CDvalue is entirely dependent on Re. The transition 133 is smooth for higher Reynolds number ows (Re >105), non-smooth for moderate Reynolds number ows 134 (200<Re<104), and accompanied by signi cant changes to a small change in M for low Reynolds number 135 ows (Re<200) [50, 51, 52, 53, 54]. With an increase in Mach number above 1, there is a rapid increase 136 inCDuntil it reaches a point where the increase is relatively small. In these high Mach number regimes, 137 CDhas a stronger dependence on M, followed by Re. Also, for low Reynolds number ows, CDis seen 138 to decrease with an increase in Mach number without reaching a maximum value near M = 1 [40]. This 139 behavior is caused by the rare ed condition, which will be discussed below. 140 The nature of Fpdragis highly dependent on the size of the particle relative to the intermolecular spacing of the carrier uid. Hence, to characterize this ratio, we employ the Knudsen number, de ned as the ratio of the molecular mean free path, , and characteristic length, L[55]. The Knudsen number is de ned as Kn = L=8 >< >:M Rep 2;if Re1; Mp 2Re;if Re>1:(8) The characteristic length, L, corresponds to the length scale through which the carrier uid interacts with 141 the particle. The characteristic length is approximately equal to the boundary layer thickness ( ) for the 142 ows where Re is higher than one and is equivalent to the particle diameter( dp) for highly rare ed ows or 143 for the ows where Re is less than or equal to 1 [56]. 144 It should be noted that some studies used particle diameter and some used particle radius as characteristic 145 length in de ning Kn [55]. However, the study conducted by Macrossan [57] indicated that Eq. 8 is a more 146 valid de nition of Kn. 147 As Kn increases, the molecular interactions at the surface of the particle can change with di erent regimes 148 of interaction de ned based on the Knudsen number: continuum ow (Kn 0:01); slip ow (0 :01<Kn 149 0:1); transition ow (0 :1<Kn10); and free-molecular ow (Kn 10). In the continuum and slip 150 ow regimes, Fpdragis caused by both pressure and viscous forces. However, the compressibility of the 151 ow modi es the pressure distribution into forebody and afterbody. For subsonic, incompressible ows, CD152 loosely follows the \standard curve" [51] where its only dependence is on Re. However, as M increases, the 153 compressibility e ects become increasingly important and must be taken into account. 154 Within transition and free molecular ow regimes, the continuum concepts of pressure and viscosity lose 155 meaning. Fpdragis better described as the di erence between the momentum imparted by the incoming 156 molecules and momentum imparted by molecules that recoil from the body [58]. Within compressible ow 157 8 behind a normal shock wave [59, 60, 61], it has been observed that CDdepends on the wall temperature of the 158 particle and Re. As the Kn increases from transition to free-molecular regimes, Fpdragbecomes increasingly 159 dependent on these parameters. 160 At low Mach numbers (M <<1),CDdecreases as Kn increases. In the slip and transition regimes, the 161 CDvalue increases as Kn increases and M decreases (for large Mach numbers). In the free-molecular regime, 162 CDis more dependent on the molecular speed ratio and the surface temperature of the particle since the 163 drag no longer depends on the viscosity of the uid [40, 54]. The molecular speed ratio is de ned as the 164 ratio of relative velocity of the body and \most probable" thermal velocity of the gas molecules [62]. The 165 values ofCDincreases with molecular speed ratio until it reaches a value of unity (high transonic ows), and 166 then decreases, approaching 2 as molecular speed ratio approaches in nity [63]. Similarly, CDincreases as 167 the ratio of the surface temperature of the particle to surrounding uid temperature increases [54, 61]. The 168 above-described behavior is usually referred as rarefaction e ects on the drag coecient. 169 In the incompressible and continuum regime, the behavior of CDis dominated by inertial/viscous e ects. 170 The compressible e ects dictate the CDbehavior for ows with higher Reynolds number (Re >4560), 171 whereas rarefaction e ects dictate the behavior for ows with low Reynolds number (Re <4560) [50, 54, 64]. 172 It can be inferred from Eq. 8 that the Kn regimes can be due to lower Re and/or higher M. Therefore, 173 rarefaction e ects occur along with inertial/viscous or compressible e ects, and vice-versa. The regions 174 where such a combination of e ects occur are the compressible regime with low to moderate Re, and the 175 incompressible regime with low Re. 176 Numerous empirical models for CDhave been proposed which cover the range of Re, M, and Kn, and these models are summarized in Table 1. The formulations of these empirical models correspond to the drag coecient of spherical particles. In this table, (Re r;Mr;Sr) and (Re1;M1;S1) indicate the Reynolds number, Mach number, and molecular speed ratio based on relative velocity and free-stream conditions, respectively. The ratio of wall temperature (in this case, the temperature of the particle, assuming it is isothermal) to the free-stream temperature of the gas is represented by Tw T1 . More speci cally, these are de ned as Rer=fjVrjdp f;Mr=jVrj af;and S r= M rr 2(9) Re1=1jV1jdp 1;M1=jV1j a1;and S1= M1r 2(10) where,,a, are the density, dynamic viscosity, speed of sound, ratio of speci c heats of the carrier uid, 177 anddpthe diameter (or size) of the particle. The subscript \ f" indicates the properties of uid surrounding 178 the particle, whereas \ 1" represents the properties of the uid in the free-stream region. 179 Stokes developed an expression for CDthrough the solution of the incompressible Navier-Stokes equation 180 for the condition where Re <<1. Oseen [65] included some convective e ects and produced an expression 181 9 with extended validity (Re <5). The Reynolds number range was further extended by Schiller and Nau- 182 man [66, 67] to Re <200 by empirically using experimental data with a similar expression also developed by 183 Walsh [68]. These formulations, being based on Stokes initial solution of the incompressible Navier-Stokes 184 equation, implicitly assume continuum ow at low M. 185 Cunningham [69] assumed a constant slip velocity of the uid, due to rarefaction, on the surface of the sphere, which is equal to (1 k)Vrwherekis the slip coecient. Cunningham derived the relation between the Kn and kwhich is expressed as k= 1 +9 2Kn1 (11) As his model is a modi cation of the Stokes solution, it is valid for Re <<1 and M<0:3, but only valid for 186 Kn0:1. 187 Carlson and Hoglund [71] developed a drag model by adding correction terms for continuum and free- 188 molecular regimes. The term for rarefaction e ect was added from Millikan's model [73] with the constants 189 determined from the experiment [74]. The inertial e ect term was the same as Schiller and Nauman's 190 model [66], and the compressibility term was taken from Hoerner's model [62]. The model works well for 191 M2 and 0:1Re100. 192 Henderson's [72] CDmodel was constructed using two sets of equations, one for M <1 and one for 193 M>1:75, with linear interpolation, applied to obtain CDin the range 1M1:75. For M<1, the Oseen 194 model [65] was used for the continuum regime, Millikan's model [73] for the slip and transition regimes, and 195 an exponential correction term added to account for compressibility e ects. For the free molecular regime, 196 the Langevin [75] and Epstein [76] models were used. Similarly, for M >1:75, an empirical relationship as a 197 function of Re 1, M1, and S1was developed using the experimental data of Refs. [61, 74]. This model works 198 well for continuum, slip, transition, and free-molecular ow for M <6 and Re up to the laminar-turbulence 199 transition 3105 . It should be noted that all the empirical models listed in Table 1 are functions of 200 Rer, Sr, and M r(based on the relative velocity of the particle). However, the Henderson CDmodel uses 201 Rer;Sr;Mrfor ows with M <1, but uses Re 1;S1;M1for ows with M >1:75. The model also includes 202 the e ect of the particle's temperature on the drag coecient through the variable ( Tw=T1). 203 Tedeschi [70] developed an empirical relationship by using the same methodology as Cunningham [69] but for an extended Re range, and applied it on the Schiller and Nauman [66] formulation. Furthermore, to include compressibility e ects, parameter Cwas added. Rarefaction e ects were introduced through a rarefaction parameter and a slip coecient k. The empirical relationships for the parameters were 10 Table 1: Previously proposed drag-coecient formulations for particles Model Formulation Stokes CD=24 Rer(12) Oseen [65] CD=24 Rer 1 +3 16Rer (13) Schiller and Nauman [66] CD=24 Rer 1 + 0:15Re0:687 r (14) Cunningham [69] CD=24 Rer 1 +9 2Kn1 (15) Tedeschi [70] CD=24 Rerkh 1 + 0:15 (kRer)0:687i (Kn)C (16) Carlson and Hoglund [71] CD=24 Rer 1 + 0:15Re0:687 rh 1 + exp 0:427 M4:63 r3 M0:88 ri 1 +Mr Rerh 3:82 + 1:28 exp 1:25Rer Mri (17) Henderson [72]CD=8 >>>>>>>>>>>>>>>>>>>>>< >>>>>>>>>>>>>>>>>>>>>:242 4Rer+ Sr8 < :4:33 +3:651:53 Tw T1 1 + 0:353 Tw T1exp 0:247Rer Sr9 = ;3 51 + exp 0:5MrpRer" 4:5 + 0:38 0:03Re r+ 0:48pRer 1 + 0:03Re r+ 0:48pRer+ 0:1M2 r+ 0:2M8 r# + 0:6Srh 1exp Mr Reri ;if M1 CD(M=1) +4 3(M11) CD(M=1 :75)CD(M=1) ;if 1<M1:75 0:9 + 0:34 M2 1 + 1:86 M1 Re11 2 2 +2 S2 1+1:058 S1 Tw T11 21 S4 1 1 + 1:86 M1 Re11 2;if M>1:75 (18) 11 developed using the experimental results of Refs. [61, 51, 77, 50]. The parameters are formulated as C(Re;M) = 1 +Re2 r Re2 r+ 100exp0:225 M2:5 r (19) (Kn) = 1:177 + 0:1770:851Kn1:161 0:851Kn1:16+ 1(20) b1k1:687+b2k1 = 0 (21) where the slip coecient kis derived by solving the above non-linear equation (Eq. 21). The coecients in the equation are de ned as b1=9 4(0:15)L rpKn 02rp LSrp Kn0:687 andb2= 1 +9 4L rpKn 0(22) whereLis the characteristic length used in de ning Kn, rpis the radius of the particle, S ris the molecular speed ratio based on the relative velocity of the particle, and 0is the correction term applied to the slip coecient with S0= (1k) Sr. The ratio between the characteristic length, L, and radius of the particle, rp, can be expressed as L rp=8 >< >:2;if Re r1 2=pRer;if Re r>1(23) and the expression for 0is given as 0=3 8p S0 1 + S02 erf S0 +exp S02 4: (24) The Tedeschi drag coecient model is intended to be valid in the continuum, slip, and transition regimes. 204 The accuracy of the Tedeschi CDmodel has been demonstrated successfully for ows with M 1 and 205 Re200 for all Kn. 206 As noted, the drag coecient for each model discussed above is given in Table 1, with each having a 207 distinct range of valid Re, Kn, and M. The objective here is to use these models to prepare a blended CD 208 model that can be implemented numerically by choosing an appropriate model for a given regime. The 209 results of this amalgamated model are then compared to experimental results for validation. The speci c 210 studies used for comparison to model predictions are presented in Table 2. 211 3. Blended Drag Coecient Model Construction 212 The collection of experimental data [78, 79, 80, 58, 77, 61] presented in Table 2, covers a range of values 213 0:12M6:4, 3Re50;000, and Kn2, that is, from subsonic to the hypersonic regime, from low 214 to high Reynolds number, and from continuum to transitional ow regime. All the experiments shown in 215 Table 2 were conducted on spherical particles. It should be noted that Aroesty's [58] experimental data for 216 12 Table 2: Drag-coecient measurements used for comparison to model results. Experiment Flow conditions Uncertainty error 2:1<M1<2:8 Kane [78] 15 <Re1<800 5% 0:1<Kn1<1 3:8<M1<4:3 Wegener and Ashkenas [79] 50 <Re1<1000 10% 0:006<Kn1<0:106 M1= 2 Sreekanth [80] 3 <Re1<60 4% 0:1<Kn1<1 M1= 2, 4, 6 Aroesty [58] 10 <Re1<10;000 DB - 0:001<Kn1<0:3 7%, 4%, 5% MM - 0:1<Kn1<1:2 7%, 6%, 5% 0:17<M1<0:99 Lawrence [77] 185 <Re1<11;600 2% 0:002<Kn1<0:1 0:12<M1<6:39 Bailey and Hiatt [61, 51] 15 <Re1<50;300 2% 0:001<Kn1<2 13 the sphere drag coecient is classi ed into two sets: the Drag Balance (DB) and the Moving Model (MM). 217 For the DB data, the drag was measured on mounted spheres in a wind tunnel. The MM data corresponds 218 to experiments where drag was measured on small spheres that were freely falling through a wind tunnel jet. 219 The MM data also consists of drag coecient values for di erent Tw=T1values. 220 Sreekanth [80], and Wegener and Ashkenas [79] provided the data in the form of a triplet (M ;Re;Kn) 221 where Kn was calculated from the measurements made. Other sets of experimental data [78, 58, 77, 61] were 222 in the form of (M ;Re) doublets, and Kn was calculated using Eq. 8. In the case of the Sreekanth [80] data, 223 Kn was calculated using the particle diameter as characteristic length for Re >1. This resulted in signi cant 224 discrepancies between CDand its theoretical free-molecular ow value at a low Kn. Similar disagreements, 225 forCDas a function of Kn, were also found for the Wegener and Ashkenas [79] data as they also used the 226 same form of characteristic length to calculate Kn. Since the formulation was chosen incorrectly for both 227 sets of data, Kn was recalculated using Eq. 8. 228 To produce a blended model for particle CDthat can be implemented numerically, the experimental data given in Table 2 was compared to di erent empirical models, as given in Table 1, at the respective ow conditions. The best performing empirical model is identi ed by comparing model values with experimental data using the relative error r= CDexpCDnum CDexp (25) whereCDexpandCDnumare the experimental and numerical values of CD, respectively. The values of rwere 229 then examined for di erent Mach number regimes, allowing selection of the best performing empirical model 230 along with the corresponding range of Re and Kn, where it produces the minimum value of r. Note that 231 some experimental data points produced values of r>50% and thus were not used for model selection. It 232 should be noted that these excluded experimental data points constitute only about 3% of the total data. 233 To minimize discontinuities in the CDpredicted by the blended model, smoothing functions are applied at the interval values of Kn and Re between the empirical models used. For example, Model 1 and Model 2 are used for ows with Re or Kn less than and greater than N, respectively, where Nis a numerical value. The smoothing function takes the form: CD= (1f)CD1+fCD2 (26) where f(N0) =log(N0)log(NN) log(N+ N)log(NN)(27) Here,CD1andCD2represent the drag coecients from Model 1 and Model 2, N0represents any value 234 in the interval [ NN;N + N]. Based on the proximity of other interfaces,  Nis chosen for each 235 interface. The smoothing function, as shown in Eq. 26, converts the discontinuities to a smooth curve. In 236 14 some instances, the discontinuities are still present. For example, there is a sudden discontinuous change in 237 CDwhen Kn crosses 1 for Eq. 16. In such cases, the CDcalculated using Eq. 26 is multiplied with a factor 238 to bridge the discontinuity. The values used for factors are determined through numerical study and occur 239 mostly for ows with M ranging between 1.7 and 5. The components of the blended CDmodel are provided 240 for the di erent Mach number ranges below. For most of the Mach number regimes, it was found that Eq. 18 241 works well for most of the cases. For the rest of the regions, Eq. 18 was replaced with other empirical models 242 whose numerical values gave a better agreement with the experimental ones. 243 The ow conditions given in Table 2 and discussed in the following sections are those of the free-stream. 244 These, however, are equivalent to relative ow conditions. The relative velocities in Aroesty's MM ex- 245 periments [58] are equivalent to particle velocities, and the relative velocities in the rest of the experi- 246 ments [80, 79, 58, 61, 77, 78], are equivalent to the free-stream velocity of the uid. 247 3.0.1. Subsonic ow (M0:8) 248 Table 3 lists the blended model components used for subsonic ow (M 0:8). In this (still compressible) 249 regime, Eq. 18 was found to best predict CDover the widest Re range in the continuum regime, and at 250 high Re in the slip regime. For Reynolds number below 2000, Eq. 16 was found to perform best for ows 251 with 0:01<Kn<0:04, and Eq. 14 for ows with 0 :04<Kn<0:1. No experimental data was found for 252 the transition regime in subsonic ow. Because of its exclusive dependency on Kn, Eq. 16 was used in the 253 transition regime. 254 Table 3: Drag coecient model components for subsonic ow Re Kn Eq. # 8Re Kn <0:01 18 Re<2000 0:01<Kn<0:04 16 Re<2000 0:04<Kn<0:1 14 Re>2000 0:01<Kn<0:1 18 The blended model prediction for CDat M<0:630 are compared to the experimental data in Fig. 1. 255 Figures 1(a) and 1(b) present the behavior of the model's CDas a function of Re and Kn, respectively, for 256 di erent Mach numbers. Figure 1(b) indicate continuum and slip regimes. As shown, the model matches 257 well with the experimental data points for all Re and Kn available within the stated uncertainty. The most 258 signi cant di erence occurs at Re = 4000 and M = 0 :204, where an error of 12% occurs. 259 15 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 1.1 1.2 1.3 100 1000 10000 100000CD ReModel Bailey & Hiatt Lawrence M=0.141 M=0.204 M=0.346M=0.423 M=0.594 M=0.630(a)CDvs. Re 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 1.1 1.2 1.3 0.001 0.01 0.1CD KnModel Bailey & Hiatt Lawrence M=0.141 M=0.204 M=0.346 M=0.423M=0.594 M=0.630 (b)CDvs. Kn Figure 1: Comparison between blended drag coecient model and experimental data for subsonic ow (M 0:8). Di erent colors indicate di erent Mach numbers. 16 3.0.2. Transonic ow (0:8<M1:2) 260 The transonic ow regime is divided into two parts: ow with Mach number below 1 and ow with Mach 261 number above 1. The components of the blended CDmodel for both regimes are presented in Table 4. 262 For transonic ow with M <1, Eq. 18 was found to be a suitable model to predict the experimental 263 results for both continuum and slip regimes, whereas Eq. 16 performed best for the transition regime. For 264 transonic ow with M 1, Eq. 17 performed best for ows at higher Re values in slip regime (Kn >0:03), 265 whereas Eq. 16 performed best for ows in transition regime for Re <200. Apart from that, CDfor the rest 266 of the regions is predicted adequately with Eq. 16 and Eq. 18 for M <1 and M>1 ows, respectively. 267 Table 4: Drag coecient model components for transonic ow M Re Kn Eq. # 0:8<M<18Re Kn <0:1 18 Kn>0:1 16 8Re Kn <0:01 18 Re<1000 0:01<Kn<0:1 18 1M<1:2 Re>1000 0:01<Kn<0:03 18 Re>1000 0:03<Kn<0:1 17 Re<200 Kn >0:1 16 Re>200 Kn >0:1 18 In Figs. 2 and 3, the blended model predictions for CDare compared to the experimental data from 268 Bailey and Hiatt [61] and Lawrence [77]. Figure 2 shows the results for Mach numbers spanning the range of 269 0:826M0:973 and Fig. 3 shows the results for Mach numbers spanning the range of 1 :013M1:164. 270 Figure 2(b) and 3(b) show slip and transition regimes. For transonic ow with M <1 (Fig. 2), most 271 of the experimental data points match well with the blended model, within the stated uncertainty of the 272 experiment. The maximum error among the points in this range was 13%, which corresponds to the data 273 from Lawrence [77] at M = 0 :896. For transonic ow with M >1, only the Bailey and Hiatt [61] data is 274 available for comparison. As shown in Fig. 3, most of the data agree with the blended model, with the 275 maximum error of 13% occurring for M = 1 :164 and Re = 95. 276 3.0.3. Supersonic ow (1:2<M3:0) 277 In order to increase the capability of the blended model for best prediction, the supersonic regime 278 is divided, and M = 1 :7 is arbitrarily selected for the division. The components of the blended model 279 are provided in Table 5. For 1 :2<M1:7, Eq. 18 was found to provide the best agreement with the 280 experiments within the continuum, slip, and transition regimes when Kn <0:08 and Eq. 17 for Kn >0:08. 281 17 0.5 0.6 0.7 0.8 0.9 1 1.1 1.2 1.3 1.4 1.5 100 1000 10000CD ReModel Bailey & Hiatt Lawrence M=0.826 M=0.867 M=0.896 M=0.918M=0.951M=0.973(a)CDvs. Re 0.5 1 1.5 2 2.5 3 3.5 4 4.5 0.01 0.1 1CD KnModel Bailey & Hiatt Lawrence M=0.826M=0.867M=0.896M=0.918 M=0.951 M=0.973 (b)CDvs. Kn Figure 2: Comparison between blended drag coecient model and experimental data for transonic ow (0 :8<M1:0). Di erent colors indicate the di erent values of Mach number. 18 0.6 0.8 1 1.2 1.4 1.6 1.8 2 2.2 2.4 2.6 2.8 10 100 1000 10000CD ReModel Bailey & Hiatt M=1.013M=1.028M=1.164(a)CDvs. Re 0.4 0.8 1.2 1.6 2 2.4 2.8 3.2 3.6 0.01 0.1 1CD KnModel Bailey & Hiatt M=1.013 M=1.028 M=1.164 (b)CDvs. Kn Figure 3: Comparison between blended drag coecient model and experimental data for transonic ow (1 :0<M1:2). Di erent colors indicate the di erent values of Mach number. 19 For 1:7<M3:0, Eq. 16 is used in transition regime for ows where Re <20 and Kn <1, and for ows 282 where 20<Re<100 and 0:5<Kn<0:6. The other regions ranging from continuum to transition regimes 283 perform well with Eq. 18. 284 Table 5: Drag coecient model components for supersonic ow M Re Kn Eq. # 1:2<M1:78Re Kn <0:08 18 Kn>0:08 17 8Re Kn <0:1 18 Re<20 Kn <1 16 Re>20 Kn >1 18 1:7<M3 20<Re<100 0:1<Kn<0:5 18 20<Re<100 0:5<Kn<0:6 16 Re>20 Kn >0:6 18 Re>100 0:1<Kn<0:6 18 Figures 4 and 5 compare the blended model predictions for CDwithin the supersonic regime to the 285 experimental data from Bailey and Hiatt [61], Aroesty [58], Sreekanth [80], and Kane [78]. The model 286 prediction for Mach numbers within the range 1 :211<M<1:698 is shown in Fig. 4 and for 1 :727<M< 287 2:991 in Fig. 5. It should be noted that the Reynolds number dependence shows non-monotonic behavior 288 within the range of 100 <Re<1000. However, the model predicts the experimental data well, and the 289 maximum error is 10% for a point from Bailey and Hiatt [61] data and for Mach number 1.474 at Re = 100. 290 At higher Mach numbers, as shown in Fig. 5, this non-monotonic behavior becomes less prevalent, and the 291 maximum error is 10% corresponding to Aroesty [58] at M = 2 :174. Figure 4(b) and 5(b) indicate slip and 292 transition regimes. The same non-monotonic behavior is also observed for Knudsen number dependence of 293 CDas shown in Fig. 4(b) 294 Two observations can made from Fig. 4. First, the CDbehavior is non-monotonic. Second, with an 295 increase in M, CDdecreases at low Re and increases at higher Re. The drag from fore-body pressure, after- 296 body pressure, and viscous components sum up to produce the total drag. The contribution of viscous drag 297 decreases with an increase in Re and M and is minimal for ows above Re 100 [62, 81]. In the supersonic 298 regime, the fore-body pressure drag decreases with an increase in Reynolds number until Re 300. After 299 that, it remains relatively constant [81]. At high Reynolds number (Re >104), the fore-body pressure drag 300 is a linear function of 1 =M2and increases with an increase in M [82]. On the other hand, the after-body 301 pressure drag increases with an increase in M until it reaches 1 [62]. After that, it shows a non-monotonic 302 20 behavior at low Reynolds number and increases with Re and decreases with M [58, 83, 84]. The rise in after- 303 body pressure drag is signi cantly smaller for ows above Re = 2 105[85]. TheCDreaches a maximum at 304 a Mach number between 1.5 and 2, after which fore-body drag becomes the main contributor [50, 62]. Thus, 305 at low Reynolds number, the viscous drag plays a signi cant role and decreases with M due to an increase 306 in rarefaction e ects [55]. On the other hand, at high Reynolds number, viscous drag is relatively small, and 307 pressure drags play a signi cant role, which increases with an increase in M. 308 It should be noted that a point of intersection of all the CDcurves can be observed in Fig. 4. This 309 behavior is due to Eq. 18, which has three sets of equations. For ows with 1 <M<1:75, Eq. 18 linearly 310 interpolates the values of CDusing the other two sets of equations. It is this interpolation equation that 311 causes the intersection of the CDcurves. It is hypothesized that Henderson [72] introduced the interpolation 312 equation to capture the transition of CDvariation with M from subsonic to supersonic ows. 313 3.0.4. Hypersonic ow (M>3:0) 314 The hypersonic ow regime is divided into two parts: one where the Mach number is below 5 and another 315 where the Mach number is above 5. For 3 <M5, Eq. 18 performs well for ows in regimes ranging from 316 continuum to transition regimes, except for a region in transition regime for Re >200 and 0:33<Kn<0:45 317 where Eq. 17 performs well. On the other hand, for M >5 ows, Eq. 18 is used for continuum, slip, and 318 transition regimes for Kn <0:45, and Eq. 17 is used for the ows with Kn >0:45. 319 Table 6: Drag coecient model components for hypersonic ow M Re Kn Eq. # 8Re Kn <0:1 18 Re<200 Kn >0:1 18 3<M<5 Re>200 0:1<Kn>0:33 18 Re>200 0:33<Kn<0:45 17 Re>200 Kn >0:45 18 M>58Re Kn <0:45 18 8Re Kn >0:45 17 Figure 6 presents the model for Mach number below 5 and Fig. 7 presents the model for Mach number 320 above 5 along with experimental data points from Bailey and Hiatt [61], Aroesty [58], and Wegener [79]. 321 The Mach numbers considered for Fig. 6 are within the range of 3 :083M4:957 , and those for Fig. 7 322 are 5:193M5:989 . Figure 6(b) and 7(b) both indicate transition regime. 323 Figure 6 shows that most of the points match well with the model, and the maximum error is 10%, for the 324 points reported by Bailey [61], for a Mach number of 4.957. However, for the latter part of the hypersonic 325 21 0.6 0.8 1 1.2 1.4 1.6 1.8 2 2.2 2.4 10 100 1000 10000 100000CD ReModel Bailey & Hiatt Aroesty - DB M=1.211 M=1.250M=1.398M=1.474M=1.573M=1.629M=1.656M=1.698(a)CDvs Re 0.6 0.8 1 1.2 1.4 1.6 1.8 2 2.2 2.4 0.01 0.1 1CD KnModel Bailey & Hiatt Aroesty - DB M=1.211M=1.250 M=1.398M=1.474M=1.573M=1.629M=1.656M=1.698 (b)CDvs Kn Figure 4: Comparison between blended drag coecient model and experimental data for supersonic ow (1 :2<M1:7). Di erent colors indicate di erent Mach numbers. 22 0.5 1 1.5 2 2.5 3 3.5 1 10 100 1000 10000CD ReModel Aroesty - DB Aroesty - MM Sreekanth Kane Bailey & HiattM=1.727 M=1.931 M=2.020 M=2.174M=2.310M=2.649 M=2.720 M=2.991(a)CDvs Re 0.5 1 1.5 2 2.5 3 3.5 0.01 0.1 1 10CD KnModel Aroesty - DB Aroesty - MM Sreekanth Kane Bailey & Hiatt M=1.727M=1.931M=2.020 M=2.174 M=2.310 M=2.649M=2.720M=2.991 (b)CDvs Kn Figure 5: Comparison between blended drag coecient model and experimental data for supersonic ow (1 :7<M3:0). Di erent colors indicate di erent Mach numbers. 23 0.8 1 1.2 1.4 1.6 1.8 2 2.2 10 100 1000 10000CD ReModel Wegener Aroesty - MM Aroesty - DB Bailey & Hiatt M=3.083M=3.230 M=3.381 M=3.659M=3.864 M=4.170 M=4.957(a)CDvs Re 0.8 1 1.2 1.4 1.6 1.8 2 2.2 2.4 2.6 0.1 1 10CD KnModel Wegener Aroesty - MM Aroesty - DB Bailey & HiattM=3.083 M=3.230 M=3.381 M=3.659 M=3.864M=4.170 M=4.957 (b)CDvs Kn Figure 6: Comparison between blended drag coecient model and experimental data for hypersonic ow (3 :0<M5:0). Di erent colors indicate di erent Mach numbers. 24 0.95 1 1.05 1.1 1.15 1.2 1.25 1.3 1.35 100 1000 10000CD ReModel Bailey & Hiatt Aroesty - MM Aroesty - DBM=5.193 M=5.616 M=5.898(a)CDvs Re 0.95 1 1.05 1.1 1.15 1.2 1.25 1.3 1.35 0.1 1CD KnModel Bailey & Hiatt Aroesty - MM Aroesty - DB M=5.193M=5.616 M=5.898 (b)CDvs Kn Figure 7: Comparison between blended drag coecient model and experimental data for hypersonic ow (M >5:0). Di erent colors indicate di erent Mach numbers. 25 regime, it was observed that the error was consistently higher for all the data points. The maximum error 326 for this case is 9%, corresponding to Aroesty [58] MM data, for M = 5 :898. 327 4. Blended model uncertainty 328 To quantify the uncertainty of the blended model, we use the relative errors between the model and each 329 set of experimental data, including the uncertainty errors. The error is de ned as 330 =8 >< >:0;ifr<u ru;ifru:(28) The errorrepresents the deviation in the percentage of the model for every set of experimental data. 331 Both a standard percentage deviation (SD) and maximum percentage deviation in errors are calculated. A 332 mean value for all the values, from every set of experimental data, is calculated. Then, an average of 333 squared di erences of all from the calculated mean is calculated for each set. The square root of this new 334 average is then used as the SD of the respective data set. The maximum of all the 's yields the maximum 335 percentage deviation. 336 The maximum error and SD for each set of experimental data calculated using Eq. 28 is given in Table 7. 337 Moreover, all the data from each set of experiments are merged into one nal set, and the maximum error 338 and SD were calculated for it. It is inferred from the calculations that the SD of the model is 2.84%. 339 The calculated SD value indicates that the of all the experimental data is distributed close to the mean 340 (=2.48%). A low SD value demonstrates the reliability of the model when predicting the value of the drag 341 coecient. It is also noted that the maximum percentage deviation is 11.87%. Henderson [72] did a similar 342 study by comparing his and Carlson & Hoglund's model [71] to several experiments [61, 73, 86]. Henderson 343 concluded that the maximum percentage deviation for Carlson & Hoglund model is 117% and for Henderson 344 model 16%. On the other hand, the blended drag model presented here has a lower maximum percentage 345 deviation of 11 :87% and is validated by comparing with a greater number of experiments that covers a wider 346 range of M, Re, and Kn numbers. 347 The best way to quantify the error in the drag coecient model is to compare the computed particle 348 motion with the experimental data. The following sections assess the validity of the model by comparing 349 the velocity and trajectories of particles to experimental data [87, 70, 88]. 350 5. Validation of Particle Trajectory Model 351 To predict the trajectory of the particle, Re, Kn, and M are used to nd CDat that instant in time using 352 the blended model described in Section 3. The net resulting force on the particle can then be found through 353 26 Table 7: Statistics of errors between the model and experiment Experiment SD Maximum error Kane [78] 2.70% 9.10% Wegener and Ashkenas [79] 0.72% 4.79% Sreekanth [80] 2.50% 8.83% Aroesty-DB [58] 2.56% 10.38% Aroesty-MM [58] 3.37% 11.10% Lawrence [77] 3.12% 10.72% Bailey and Hiatt [61, 51] 2.56% 11.87% All 2.84% 11.87% Eq. 5 and Eq. 4, which is then inserted into Wof Eq. 1 to determine the momentum change of the particle. 354 The system of equations represented by Eq. 1 is solved that calculates mass, velocity, and temperature of 355 the particle for each time step. The velocity of the particle is equal to the rate of change of position. The 356 discretized version is inserted in Eq. 1 and also solved to obtain the position of the particle. 357 To validate the particle trajectory model, we compare the predicted velocity (Section 5.1) and the tra- 358 jectory (Section 5.2). The experiments chosen use supersonic ows over wedge samples. 359 Flows that have large velocity gradients provide an opportunity to measure the dynamics of the particle. 360 It is particularly true for ows containing an aerodynamic shock. Upstream of the shock, one can assume 361 the particle velocity matches that of the uid and Vr= 0. However, when crossing the shock, there is an 362 abrupt change in ow velocity. However, the inertia of the particle prevents it from responding instantly. 363 In the downstream region, the velocity of the particle is initially the same as the upstream velocity. At this 364 point, Vr6= 0 and is referred to as the particle slip error, and the distance the particle travels before jVrjis 365 less than 1% of the surrounding uid velocity ( jVfj) is de ned as the relaxation distance. 366 A compressible ow, with a Mach number greater than 1, over an inclined surface such as a wedge, 367 produces an oblique shock and an expansion fan. When a particle passes through this region, it encounters a 368 range of velocity gradients, which is a suitable domain not only for measuring the dynamics but also to test 369 a particle trajectory model. As the particle travels through the relaxation distance, it experiences signi cant 370 changes in its dynamics, and it is seen that these changes magnify the errors of the model. 371 The tracer particles employed in the experiments used for validation were non-reactive. Therefore, the 372 chemical interactions of the particles are not taken into consideration. 373 27 5.1. Velocity 374 The experiments performed by Thomas et al. [87, 45] and Tedeschi et al. [70, 47, 46] were used to validate 375 the model velocity prediction. 376 Thomas et al. [87] used seeding particles of a known diameter in a supersonic ow over a wedge sample and 377 measured their average velocity using laser doppler anemometry (LDA). The experiments were conducted 378 for three di erent wedge samples with half-angles of 6.3, 12, and 14. Di erent types of seeding particles 379 were also used: atomized olive oil, incense smoke, Titanium oxide (TiO 2), and \Blanc Fixe, Micro" (re ned 380 BaSO 4). The experiment using \Blanc Fixe, Micro" as seeding particles was only performed for M = 1 :95 381 in air, over a wedge with a half-angle of 14, and is considered here for validation. The test conditions for 382 the experiment are given in Table 8, and the test geometry is shown in Fig. 8(a). Table 8 lists the values of 383 temperature ( T1), density ( 1), mass fractions ( Y), velocity ( V1), and Mach number (M 1) of air in the 384 free-stream (or upstream) and the wall (or surface) temperature ( Tw) of the wedge. 385 Tedeschi et al. [70] used LDA to measure seeding particle velocities in a supersonic ow. However, in 386 this case, the ow was over the shock generator. The shock generator consisted of a plate that is inclined 387 at an angle of 8in the free stream direction, as illustrated in Fig. 8(b). The seeding particles used for the 388 experiments were incense smoke and latex spheres. The uncertainty for experiments was given as 0.1% in 389 thex-direction and 4% in the y-direction. For the validation, the experimental results obtained using latex 390 spheres as seeding particles for a ow of M = 2 :3 and total pressure of 0 :5105in the air were considered. 391 The test conditions for the experiment are given in Table 8. 392 h14◦xy O (a) Thomas et al. Flowθ= 8◦O lxy Shock generator (b) Tedeschi et al. Figure 8: Wedge sample used by Thomas et al. [87] (where h= 0:07 m) and wind-tunnel setup for experiment used by Tedeschi et al. [70] (where l= 0:167 m) where \O" is the origin. To validate the particle trajectory model, the converged steady-state CFD solutions were computed using 393 KATS for test conditions of both the experiments, as given in Table 8. The ow eld data was extracted and 394 used in the particle trajectory code to determine Vpvelocity pro les. For the simulation, the average size of 395 the seeding particles in the respective experiments and the free-stream conditions of the ow were taken as 396 initial values. Single-particle simulation is performed, and the particle's velocity variation with respect to 397 28 Table 8: Test conditions of LDA experiments Experiment T1, K1, kg/m3YN2YO2V1, m/s M1Tw, K Thomas et al. [87] 167.0 0.2804 0.7671 0.2329 505.92 1.95 295.2 Tedeschi et al. [70] 147.0 0.0956 0.7671 0.2329 559.20 2.30 303.0 its position in the x-direction is compared with the experimental data. 398 The solution of the ow eld for the experiments of Thomas et al. [87], obtained using KATS, is presented 399 in Fig. 9(a). The solution of the ow eld for the shock generator experiment of Tedeschi et al. [70] is shown 400 in Fig. 9(b). It is to be noted that both ow eld solutions presented in Fig. 9 are grid-independent. The 401 velocity measurements for the Tedeschi et al. [70] experiment were performed across the region of oblique 402 shock wave induced by the inclined shock generator. Since the ow over the inclined surface is equivalent 403 to the ow over a wedge with 8half-angle, the latter was used to simulate the ow eld. Fig. 9(a) shows 404 the ow velocity across the wedge where the velocity decreases from free-stream after the oblique shock 405 and increases to free-stream (and more) as it goes through the expansion fan. Although the same general 406 behavior is observed in Fig. 9(b), the ow velocity does not decrease as much in the region between oblique 407 shock and expansion fan due to the large size of the wedge. 408 (a) Thomas et al. [87] (b) Tedeschi et al. [70] Figure 9: Velocity contours of the simulated ow eld solution corresponding to test conditions of the experiments. 29 Thomas et al. [87] used \Blanc Fixe, Micro" particles with an average size of 0.7 m and a density of 409 4400 kg/m3, and the same values were used for simulation. These particles are of the same magnitude when 410 compared to those present in rocket nozzle ows (0 :110m) [89, 90] as well as those in Martian dust 411 (0:11m) [4, 5] . Figure 10 compares the results of the simulation to the experiment for the velocity of 412 the particle. 413 In Fig. 10(a), the experimental results indicate that the velocity of the particle decreases and increases 414 near the shock location. However, it should be noted that this uctuation is not physical, and was attributed 415 to the optical distortion of the LDA beams caused by the shock [87]. Since the vertical velocity at this location 416 is low, the measurements do not see any such uctuations in the vertical direction, as shown in Fig. 10(b). 417 For both vertical and horizontal velocity components, the numerical model produces excellent agreement 418 with the experimental results. 419 Tedeschi et al. [70] used latex particles with an average size of 0.523 m and a density of 1050 kg/m3, 420 and the simulation uses the same values. The results for the longitudinal velocity and transversal velocity 421 are presented in Fig. 11. It is observed that simulated results show good agreement in predicting the 422 Tedeschi et al. [70] results. Figure 11(a) shows that for the longitudinal velocity, the simulated trajectory 423 agrees well with the experimental data. Similarly, the transversal velocity of both simulation and experiment 424 match each other well, as seen in Fig. 11(b). 425 The location of velocity measurements for each experiment was at a transversal distance from the center 426 axis of the wind tunnel ranging from upstream to downstream (after the expansion fan) of the shock. 427 Therefore, these experimentally measured velocities are a function of speci c locations in space. On the 428 other hand, the particle model follows a Lagrangian frame reference, and the simulated trajectory does 429 not necessarily pass through all the measurement locations. However, it was observed that the maximum 430 di erence between the trajectory and measurement points was of the order of magnitude 103m for both 431 Thomas et al. [87] and Tedeschi et al. [70] experiments. Figure 12 shows the location of measurement points 432 and simulated trajectories based on Thomas et al. [87] and Tedeschi et al. [70] experimental conditions. The 433 variation of simulated ow eld velocity ( VfV1), for both Thomas et al. [87] and Tedeschi et al. [70], are 434 also provided in Fig. 12, and illustrates the abrupt change in the value after the shock. Since the simulated 435 points were very close to the measurement points, it can be inferred that the velocities predicted by the 436 model closely resemble the experimental ones. 437 It should be noted that, irrespective of the free-stream ow regime, the particle's Mach number M p (= M r) starts from 0 and increases to subsonic, then to transonic, and so on when it crosses the shock or ejects from a surface. Figure 13 illustrates the simulated Re p, Mp, and Kn pbased on Thomas et al. [87] and Tedeschi et al. [70] experimental conditions. It can be seen that as the particle travels and crosses the shock, M pincreases to 0 :45 and 0:2 for Thomas et al. [87] and Tedeschi et al. [70], respectively. These points 30 380 400 420 440 460 480 500 520 540 -0.004-0.002 0 0.002 0.004 0.006 0.008 0.01up, m/s x, mUncertainty Range Experimental Data Numerical Data(a) Longitudinal velocity -20 0 20 40 60 80 100 120 -0.004-0.002 0 0.002 0.004 0.006 0.008 0.01vp, m/s x, mUncertainty Range Experimental Data Numerical Data (b) Transversal velocity Figure 10: Comparison of longitudinal and transversal components of particle's velocity between simulated and the experiments of Thomas et al. [87]. The uncertainty range was calculated by adding the stated experimental uncertainty to the average value of the velocity at every given location. 31 480 490 500 510 520 530 540 550 560 570 580 0 0.002 0.004 0.006 0.008 0.01 0.012 0.014up, m/s x, mUncertainty Range Experimental Data Numerical Data(a) Longitudinal velocity 0 10 20 30 40 50 60 70 0 0.002 0.004 0.006 0.008 0.01 0.012 0.014vp, m/s x, mUncertainty Range Experimental Data Numerical Data (b) Transversal velocity Figure 11: Comparison of longitudinal and transversal velocity components of particle's velocity between simulated and exper- iments of Tedeschi et al. [70]. The uncertainty range is calculated by adding the error of 4% to the measured mean velocity values. 32 -0.0025-0.002-0.0015-0.001-0.0005 0 -0.01 -0.005 0 0.005 0.01 0.015-100-80-60-40-20 0∆y, m ∆Vf, m/s x, mMeasured points ∆y - Thomas ∆Vf - Thomas ∆y - Tedeschi ∆Vf - TedeschiFigure 12: Comparison of simulated trajectories based on Thomas [87] and Tedeschi [70] experiments with respect to location of measurement points. of maximum M pdenote the downstream regions after the shock, as can be seen by Vfpro les in Fig. 12. The M ppro les then decrease before the particle adapts to the ow eld velocity. Similarly, Kn pchanges from a free-stream value and decreases while the particle crosses the shock. As the particle travels through an expansion fan, M pand Kn pagain increase to di erent values and follow the same behavior after that. It is interesting to note that even though M pand Re pare zero in the free-stream region, Kn phas a constant value which can be calculated by inserting Eq. 6 and 7 into Eq. 8 to obtain Knp= lim jVrj!0Mp Repr 2=1 1dpa1r 2(29) where1,a1,1represent free-stream density, speed of sound, and dynamic viscosity of the uid, whereas 438 dpis the particle size. As the particle travels through the computational domain of the ow eld as shown 439 in Fig. 9, for both the experiments, the simulation shows that the particle is in transition regime, and Kn p 440 ranges between 0 :2 and 0:8 for Thomas et al. [87] and between 1 and 4 :5 for Tedeschi et al. [70]. 441 The results, as shown in Fig. 10 and 11, suggest that the trajectory model captures the velocity changes 442 accurately { despite the particle being in non-continuum regime { and successfully reproduces the particle 443 velocities with M pin the subsonic regime. 444 33 0 0.5 1 1.5 2 2.5 3 -0.005 0 0.005 0.01 0.015 0 0.05 0.1 0.15 0.2 0.25 0.3 0.35 0.4 0.45Rep, Knp Mp x, mRep - Thomas Mp - Thomas Knp - Thomas Rep - Tedeschi Mp - Tedeschi Knp - TedeschiFigure 13: Variation of Reynolds number, Mach number, and Knudsen number of particles travelling through ow elds based on Thomas [87] and Tedeschi [70]. 5.2. Trajectories 445 The comparison of the measured particle velocity with simulated results is not sucient to evaluate the 446 accuracy of the trajectory model. It is also necessary to extend the validation to particle trajectories. The 447 experimentally determined particle trajectories can show an accumulation of error due to uncertainty in 448 velocity measurements. The comparison with the simulated trajectory and its closeness to the experimental 449 one determines the capability of the trajectory model in predicting the particle path and its properties. 450 In order to validate the particle trajectories, the velocity measurements conducted by Ross et al. [88, 91] 451 were used. The measurements were obtained using Particle Image Velocimetry (PIV) and examined the slip 452 velocity of ow around a wedge. The experiments were performed in a high-speed blowdown wind tunnel, 453 and a converging-diverging axisymmetric nozzle was used to produce uniform free-stream of velocity M = 2. 454 Four di erent wedge samples were used in the experiment, shown in Fig. 14. They had the same base 455 thickness (h), but with di erent half-angles (5, 10, 15, and 22:5). It can be seen that the size of the 456 wedge samples decreases with the increase in the half-angle of the wedge. For these experiments, aluminum 457 oxide (Al 2O3) powder was used as the seeding material. 458 Ross et al. [88] calculated the ow pathlines by integrating the velocity data from the experiments. An 459 imaginary particle position in the upstream region was initially considered, with the velocity interpolated 460 34 hO 5◦O 10◦O15◦O22.5◦xyFigure 14: Con guration of 5, 10, 15, and 22 :5half angle wedge airfoils with a constant base thickness h(= 0.25 in) where \O" is the origin. Table 9: Test conditions of PIV experiments CaseT1, K1, kg/m3YN2YO2V1, m/s M1Tw, K 1 160.0 2.197 0.7671 0.2329 508.076 2.0 288.71 2 295.5 1.210 0.7671 0.2329 690.338 2.0 533.15 from the data. At a  tlater, the velocity and position were calculated from the initial position, and the 461 calculation procedure repeated until the particle crossed the data region. The experiments were conducted 462 for the test conditions given in Table 9. While keeping the free-stream Mach number as 2, the two tests were 463 conducted by changing the total temperatures ( Tw), which in turn changed the free-stream temperatures 464 (T1) and velocities ( V1). It can be noticed from Table 9 that Case-1 represents a low-temperature ow, 465 whereas Case-2 corresponds to a heated ow. Though the free-stream Mach number is the same, the shock 466 strength is di erent. The velocity and temperature gradients are high for Case-2 when compared to Case-1. 467 Thus, the particles traveling in Case-1 have a higher Re and lower Kn and, therefore, a di erent CD. 468 The errors in the experiments, as given by Ross et al. [88], were typically 3%. Therefore, the uncertainty range along the particle trajectory is calculated as x=rx 0:5 (0:03)2andy=ry 0:5 (0:03)2(30) where (x;y) is the position of the experimental points, and 0.5 is the magni cation factor used on the images 469 to reduce the signal-to-noise ratio. 470 The steady-state solutions of the ow eld over di erent wedge samples were simulated using KATS-CFD 471 based on the test conditions in Table 9. The results of the grid-independent solutions are shown in Fig. 15. 472 It can be seen in Fig. 15 that as the half-angle of the wedge increases from 5to 22:5, there is an increase 473 in the velocity intervals across the sample. The maximum de ection angle for a Mach 2 ow is 22 :9, and 474 therefore, the shock appears to be almost detached for the 22 :5half-angle wedge sample. 475 Both the simulation and the PIV experiments were conducted using an average particle size of 0 :8m, 476 and a density of 4000 kg/m3. Figure 16 presents the simulated and experimental trajectories for di erent 477 35 (a) 5half angle wedge (b) 10half angle wedge (c) 15half angle wedge (d) 22 :5half angle wedge Figure 15: Velocity contours of Mach 2 ow, for two test conditions, over 5, 10, 15, and 22 :5half angle wedges 36 wedge samples. Figure 16(a) shows the trajectories for test conditions corresponding to Case-1, whereas 478 Fig. 16(b) corresponds to Case-2. 479 It can be observed that the simulated trajectories pass through most of the PIV data points and are 480 inside the uncertainty range for 5, 10, and 15half-angle wedge samples. However, there seems to be a 481 disagreement for the 22 :5half-angle wedge sample in both cases. Maxwell and Seasholtz [92] conducted 482 studies that claimed that the relaxation distance for ows with normal Mach numbers, ahead of the shock, 483 ranging from 1.15 to 1.43, was relatively constant. Based on the measured data of other samples, it was 484 noted that the size of the 22 :5half-angle wedge sample is of the same order as the relaxation length of the 485 particles. This resulted in particles reaching the re-accelerating region of the expansion fan before they could 486 adapt to the velocity change in the uid, thus producing a relaxation error. These particle relaxation errors 487 caused the overlapping of velocity vectors, giving inconsistent measurements for velocity after the shock. 488 Due to skewed velocity measurements for the 22 :5wedge sample, the pathline calculated did not accurately 489 depict the path traversed by the tracer particle. 490 Other than the skewed data for the 22 :5wedge sample, the particle-tracking code was able to match 491 the trajectories for other samples accurately, for both cases. Therefore, it can be implied that the simulated 492 trajectory for a 22 :5wedge sample resembles the one followed by the tracer particle. 493 Similar to Fig. 13, Fig. 17 show the range of Re p, Mp, and Kn pfor the test-cases used to validate the 494 particle trajectories. These plots show that the particles travel from slip regime to transition regime with 495 Kn ranging from 0.02 to 0.14 for Case-1 and from 0.06 to 0.22 for Case-2, respectively. The behaviors of Re p496 and M pfor di erent wedge con gurations are similar to the ones observed in Fig. 13. 497 6. Conclusions 498 To predict particle trajectories in complex hypersonic ows, drag coecients needs to be accurately 499 modeled over a wide range of Mach number, Reynolds number, and Knudsen number. Since no single 500 empirical model is able to reproduce the drag coecient behavior over these diverse regimes, a blended model 501 was constructed. The new model combines the best performing empirical model for di erent combinations 502 of regime. Experimental data that envelops a wide range of Re, M, and Kn was selected for comparison, 503 and the model agrees within a standard deviation and a maximum error of 2.84% and 11.87%, respectively. 504 The blended model was then inserted into a Lagrangian particle trajectory model, and validations were 505 performed for particle velocity and trajectory. The comparisons for ten di erent experimental test-cases 506 were in excellent agreement, with results well within the uncertainties of the experiments. These results 507 provide con dence that a Lagrangian particle trajectory model, equipped with the blended drag coecient 508 model, can accurately predict particle dynamics. Moreover, the approach is numerically ecient, allowing 509 it to be used to predict the particle dynamics for particle-laden ows, at least for cases where inter-particle 510 37 0.0000.0010.0020.0030.0040.0050.0060.007-0.0025-0.0020-0.0015-0.0010-0.00050.0000y, m x, m PIV Data PIV Uncertainty Range Numerical Data5o10o15o22.5o(a) Case - 1 0.0000.0010.0020.0030.0040.0050.0060.007-0.0025-0.0020-0.0015-0.0010-0.00050.0000y, m x, m PIV Data PIV Uncertainty Range Numerical Data5o10o15o22.5o (b) Case - 2 Figure 16: Comparison of simulated trajectories with the ones retrieved from PIV data 38 0 1 2 3 4 5 6 -0.001 0 0.001 0.002 0.003 0.004 0.005 0.006 0.007 0 0.02 0.04 0.06 0.08 0.1 0.12 0.14 0.16Rep Mp, Knp x, mRep - Case 1 Mp - Case 1 Knp - Case 1 Rep - Case 2 Mp - Case 2 Knp - Case 2(a) 5half angle wedge 0 2 4 6 8 10 12 14 -0.001 0 0.001 0.002 0.003 0.004 0.005 0.006 0.007 0 0.05 0.1 0.15 0.2 0.25 0.3Rep Mp, Knp x, mRep - Case 1 Mp - Case 1 Knp - Case 1 Rep - Case 2 Mp - Case 2 Knp - Case 2 (b) 10half angle wedge 0 5 10 15 20 25 -0.001 0 0.001 0.002 0.003 0.004 0.005 0.006 0.007 0 0.05 0.1 0.15 0.2 0.25 0.3 0.35 0.4 0.45Rep Mp, Knp x, mRep - Case 1 Mp - Case 1 Knp - Case 1 Rep - Case 2 Mp - Case 2 Knp - Case 2 (c) 15half angle wedge 0 5 10 15 20 25 30 35 40 45 -0.001 0 0.001 0.002 0.003 0.004 0.005 0.006 0.007 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7Rep Mp, Knp x, mRep - Case 1 Mp - Case 1 Knp - Case 1 Rep - Case 2 Mp - Case 2 Knp - Case 2 (d) 22 :5half angle wedge Figure 17: Variation of Reynolds number, Mach number, and Knudsen number of particles travelling through ow elds for di erent half-angle wedge con gurations and for two test conditions 39 interactions and structural deformations can be neglected. 511 The new blended model allows to predict the interactions of the particles with the ow eld in situations 512 where those capabilities are needed, such as the dust around an atmospheric entry capsule, the spalled 513 particles of an ablating heat shield, or the metal oxide particulate inside a solid rocket plume. 514 Acknowledgements 515 Financial support for this work was provided by NASA Kentucky EPSCoR Award NNX10AV39A and 516 NASA Award NNX13AN04A. The authors acknowledge G. Tedeschi for providing access to his thesis. The 517 rst author would also like to thank A. D. Ahmad, U. D uzel, and S. Schmitt for their assistance, as well as 518 the librarians in charge of the InterLibrary Loan (ILL) at the University of Kentucky for going out of their 519 way to provide the necessary material. 520