Analysis and comparisons of various models in cold spray simulations : towards high fidelity simulations
AAnalysis and comparisons of various models in cold spray simulations: towards high fidelity simulations
Louis-Vincent Bouthier ∗ , Elie Hachem † August 19, 2020
Abstract
Cold spray technology is a quickly growing manufacturing technology which impacts lots of industries. Despite manyyears of studies about the comprehension of the phenomena and the improvements of the performance of the system, ensur-ing high fidelity simulations remains a challenge. We propose in this work a detailed high fidelity modeling and simulationsgiving more insight of the phenomena appearing in cold spray such as turbulence, oblique shocks, bow shocks, fluctuations,particles motion and particles impacts. It is mainly based on a richer model known as the Detached Eddy Simulation (DES)model. Moreover, we present several analysis of various existing models for both validations and comparisons purposes.Finally, this high fidelity framework will allow us to deal with a new configuration showing an improved performanceassessed with the previous models.
Keywords: Cold spray,
Fluent , RANS, IDDES, Turbulence, High-fidelity, CFD
Notations ‖⋅‖
Norm of a vector ⋅ Mean of a quantity in RANS 𝐶 Correction factor taking into account compressible effects (-) 𝐶 𝐷 Drag coefficient (-) 𝑐 Speed of sound ( m.s −1 ) 𝑐 𝑝 Mass thermal capacity of solid particles (
J.kg −1 .K −1 ) 𝐷 Nozzle diameter ( m ) 𝑑 𝑝 Solid particle diameter ( m ) dd𝑡 Particle derivative for time- and space-dependent quantities ( s −1 ) 𝑑 𝑤 Distance to the nearest wall ( m ) 𝑒 Internal energy per unit of mass of the gas (
J.kg −1 ) 𝜀 Turbulent dissipation rate ( m .s −3 ) 𝜂 Ratio of the length of the stagnation chamber to the length of the convergent (-) 𝛾 Ratio of the specific heats (-) ℎ max Maximum edge length of a cell ( m ) 𝑘 Turbulent kinetic energy ( m .s −2 ) ∗ Ecole Polytechnique, Route de Palaiseau, 91120 Palaiseau, France,email: [email protected] † MINES ParisTech, PSL - Research University, CEMEF - Centre for material forming, CNRS UMR 7635, CS 10207, rue Claude Daunesse, 06904 Sophia-Antipolis Cedex, France a r X i v : . [ phy s i c s . f l u - dyn ] A ug Volume viscosity relative to expansion ( kg.m −1 .s −1 ) 𝐿 Total nozzle length ( m ) 𝜆 Thermal conductivity of gas (
W.m −1 .K −1 ) 𝜆 𝑝 Thermal conductivity of solid particles (
W.m −1 .K −1 ) 𝑀 Mach number (-) 𝜇 Dynamic gas viscosity ( kg.m −1 .s −1 ) 𝜇 𝑡 Turbulent gas viscosity ( kg.m −1 .s −1 ) Nu Nusselt number (-) 𝜔 Specific turbulent dissipation rate ( s −1 ) Ω Magnitude of the vorticity tensor ( s −1 ) 𝑝 Gas pressure ( Pa ) 𝑝 𝑐 Critical gas pressure ( Pa ) Pr Prandtl number (-) 𝑟 Ideal gas specific constant (
J.kg −1 .K −1 )𝑅 𝑖𝑗 Reynolds tensor ( kg.m −1 .s −2 ) Re Reynolds number (-) 𝜌 Density of gas ( kg.m −3 ) 𝜌 𝑐 Critical density of gas ( kg.m −3 ) 𝜌 𝑝 Density of solid particles ( kg.m −3 ) 𝑆 Magnitude of the strain rate tensor ( s −1 ) 𝑡 Time variable ( s ) 𝑇 Gas temperature ( K ) 𝑇 𝑐 Critical gas temperature ( K ) 𝑇 𝑝 Particle temperature ( K ) 𝒖 Gas velocity ( m.s −1 ) 𝒖 𝒑 Solid particle velocity ( m.s −1 ) 𝜐 Acentric factor (-) 𝜉 𝑖 Ratio of the inlet diameter to the throat diameter (-) 𝜉 𝑜 Ratio of the outlet diameter to the throat diameter (-) 𝜁 Ratio of the length of the divergent to the total length of the nozzle (-)2igure 1: Diagram of a convergent-divergent channel De Laval (cf (Ref [1]))
Cold spray is a manufacturing process that began at the end of the 20 th century. It consists of manufacturing objects witha wide range of materials, from plastics to metals and more complex materials. Its fields of application are vast and includestart-ups and fab-labs for prototyping, aeronautics, construction or the automotive industry. The stakes of this technologyare multiple because it allows the design of complex geometry parts relatively quickly and often at a lower cost. Nevertheless,limits appear with regard to the cost of machines and materials, the surface state of the parts after the operation or theinternal state of the built structure in terms of residual stresses, porosity, solidity and control of reproducibility.In this work, we provide a rich overview on the cold spray technology in terms of experiments, numerical simulationsand challenges to tackle, presented by the following sub sections. While it is focusing mainly on high-pressure cold spraybut the review remains rather general.We recall first that the operating principle of dynamic cold spray is relatively simple and has been studied since the 1980sin (Ref [2]). The system requires the use of a convergent-divergent channel (also called De Laval) as shown in Fig 1 in orderto reach supersonic velocities. The main gas, which can be air, nitrogen or even helium, is introduced at the inlet of thechannel with a high pressure 𝑝 ≈ 50 bar and a high temperature 𝑇 ≈ 800 ◦ C to prime the nozzle and obtain supersonic flow.The inlet pressure and temperature conditions are highly dependent on the nozzle geometry, the substrate to be coated, theparticles to be sprayed and the desired coating condition.The previously introduced gas will be used to accelerate solid particles with an average size 𝑑 𝑝 ≈ 30 µ m made of a metalalloy, ceramic or polymer (see (Ref [3])). Despite the variability in application conditions, most authors agree that there is acritical particle impact velocity beyond which adhesion to the substrate can occur (see (Ref [4, 5, 6, 7, 8, 9])). The value of thiscritical velocity under experimental conditions is still unknown and the previous authors each mention an empirical modelto estimate it as a function of the temperature at impact and the characteristics of the material. Moreover, it is essential torecall that the adhesion of particles to the substrate is achieved in the solid state, in particular thanks to the strong plasticdeformation of the particles due to their impact velocity. In (Ref [10]), for example, the authors simulated the impact andadhesion of a particle by taking into account the data from their previous study using a finite element calculation.In summary, the goal of the whole system is to accelerate solid particles using a supersonic gas flow to make them adhereto a substrate due to their high impact velocity. The challenge is then to understand what are the relevant parameters toconsider in order to optimize this process according to the specifications. We will establish a work similar to the oneperformed in (Ref [9]) to assess the current knowledge. As mentioned above, a necessary issue is the establishment of supersonic flow within the nozzle. In an approximate way, in(Ref [9, 10, 11]), the authors presented the relationships in the De Laval nozzle in the case of a perfect fluid flow allowing tohave a first order of magnitude of the parameters involved. Knowing that a certain velocity is desired for the solid particles,it thus fixes a gas velocity inside the nozzle. By the preceding relations, this velocity is simply given, as a first approximation,by the ratio of the diameter of the nozzle outlet to the diameter of the throat. Then, all the thermodynamic quantities arededuced from this speed and thus make it possible to know in an approximate way the behavior of the gas inside the nozzle.This method, although effective, neglects phenomena such as turbulence illustrated in (Ref [12]) or shocks. This is whythe majority of the authors used in parallel numerical simulation tools to account more finely for the phenomena involved.3owever, it is worth mentioning that a large part of the authors modeled gases in supersonic flow by the law of ideal gas.The relevance of this law in many fields is no longer to be proven, however, in the context of cold spray, it can vacillate. Inthe context of high pressures and high temperatures, so-called real state laws would undoubtedly be better able to predictthe phenomena precisely, knowing also that many calculations are made numerically. Examples of such laws are Van derWaals’ state law (see (Ref [13])), the state law presented in (Ref [14, 15]) and its version modified in (Ref [16]) or the statelaw presented in (Ref [15, 18]). In fact, the numerical modeling of turbulence plays a major role in all previous works dealing with the cold spray simulation.While some approached the subject by laminar flow as in (Ref [10]) or chose to directly simulate the Navier-Stokes equations(i.e. DNS) as in (Ref [12]) , many models were used such as the standard 𝑘 − 𝜀 model in the Reynolds averaged equations (orRANS) with (Ref [3, 19]), the realizable 𝑘 − 𝜀 model with (Ref [6]), the Reynolds Stress model (see (Ref [9])) or others notablyin (Ref [20]). These models have the advantage over a DNS to reduce considerably the computation time: during a DNS,the number of nodes of the mesh is proportional to Re and the simulation is necessarily unsteady and three-dimensionalwhereas a turbulence model with RANS can be driven in two dimensions and in steady state, which limits considerably thecomputation time. The pieces of advice in (Ref [9]) on comparing these turbulence models for each application are thereforerelevant.However, the models are not limited to turbulence modeling alone. The viscosity of the gas, for example, can be consid-ered constant or temperature dependent according to the law proposed in (Ref [21]) and used in (Ref [22]) or Maxwell’s law(see (Ref [23, p.25])).An additional remark must be made concerning Stokes’s hypothesis. Indeed, this hypothesis imposing the nullity of thevolume viscosity, notably used in (Ref [12]), is relevant in the context of incompressible flows, i.e. with a very small Machnumber in front of the unit and low density monoatomic gas where it is even a property (see (Ref [15, 23, 24])). However,in (Ref [25, 26, 27]), the authors presented room temperature measurements of 𝜅/𝜇 for polyatomic gases such as nitrogen,used in cold spray. All these measurements show that this ratio is often of the order of the unit or even much higher in thecase of dihydrogen and invite to question the Stokes hypothesis in the case of gases at high temperatures and high pressuresused in cold spray.Also, thermomechanical coupling is unavoidable in compressible flow which adds more models to take into account allrelevant phenomena. In (Ref [20]), the author cited for example the 𝑘 𝜃 − 𝜀 𝜃 model, the 𝜀 𝑡 model of turbulent dissipation orthe 𝑅 𝑖𝑗 − 𝜀 models with the Reynolds tensor. The complexity of the flow along the nozzle is then a challenge but one of the objectives of the cold spray technology isto accelerate solid particles to adhere onto the substrate. Therefore, the gas flow in the nozzle is not single-phase and thesupply of solid particles to be accelerated is not insignificant. The coupling between the particles and the gas can be takeninto account or not. In (Ref [28]), the authors considered that the particles have no influence on the flow of the fluid due tothe very low mass proportion of particles in the flow. On the contrary, in (Ref [29]), the authors set up a fully coupled modelbetween compressible flow and particle displacement.
This same displacement, calculated by numerical simulations, has variations according to the numerical scheme used. Par-ticle tracking can be carried out according to the Lagrangian description (see (Ref [10])) or according to the Eulerian de-scription (see (Ref [11])). This displacement is also subject to an equation involving the compressible drag force of the flow.From this force comes a large number of models seeking to evaluate the drag coefficient. These models take into account, forexample, the fact that the particles under consideration may have irregular shapes, that the Mach number has a dominantinfluence especially at the location of the shock waves, or that the Reynolds number is also relevant in the study. For anoverview, in (Ref [8, 9, 10, 11, 28, 30, 31, 32]), each author presented a number of models for the drag coefficient of a particle.In (Ref [29]), the authors also proposed to consider the effect of shock waves by adding an additional force to a particlerelated to the local expansion of the fluid. This law is notably available in
Ansys’ Fluent software (see (Ref [17])) Even if the authors considered their simulation as DNS using
Fluent , the
Ansys developers do not encourage the customers to use their softwares forDNS (see (Ref [17])). .3.2 Particle thermodynamics Second, since we are considering a cold spray process, thermal effects within the particle are paramount. Almost all authorsagreed to consider the temperature within a particle as uniform. This fact is related to the very low value of the Biot numberwhich compares the exchange coefficient with the flow and the internal conductivity of the particle. However, despite thisgreat simplification, the exchange coefficient between the flow and a particle must be evaluated to account for the heatexchanges between the two bodies. In contrast to drag coefficient models, few exchange coefficient models are used by theauthors and the main one is the Ranz-Marshall model (see (Ref [33])) involving the Nusselt number, the Prandtl number andthe Reynolds number.
Furthermore, the interaction between the nozzle walls and the particles is still unknown. In (Ref [8]), the authors illustratedthe effects of nozzle clogging which are strongly dependent on boundary conditions, especially pressure. Together withthe authors of (Ref [12]), they suggested cooling the nozzle walls to limit clogging. In (Ref [5]), the authors demonstratedexperimentally that the material covering the nozzle walls has a significant impact on the performance of cold gas spraying,all else being equal. They interpreted this effect through variations in the thermal conductivity of the wall material. Inaddition, many authors assumed elastic rebound in three-dimensional simulations, but, in (Ref [31]), the authors proposedto use elastic-plastic impact restitution coefficients.
We present in Tab 1 a new and an interesting summary on all the numerical details used until now in the literature. Wecan notice some similarities in the previous works. We have added our approach in the last row to clearly highlights thenovelty proposed in this work. Another comment must be made about the drag coefficients. Indeed, except for the Morsiand Alexander model with compressibility correction, the Crowe model and the Henderson model (see (Ref [9])), none ofthe previous drag coefficients takes into account compressible effects which could be not relevant in the shocks for example.
After all the phenomena being presented, we recall that many researchers and engineers aimed at optimizing the nozzles andapplication conditions to obtain the best performing coating. For example, in (Ref [34]), the authors focused on the impactof the spray angle, in (Ref [7]), the authors studied the effect of the pressure with which the particles were injected intothe flow and, in (Ref [32]), the authors wanted to build a nozzle of reduced size in order to obtain refined coating grooves.Empiricism and parametric studies dominate the papers. In (Ref [32]), the authors mentioned a totally empirical rule ofthumb for constructing a cold spray nozzle with respect to Mach number and the ratio of divergent length to nozzle outletdiameter. In order to give an overview of the multiplicity of models, the Tab 2 gathers parameters characterizing the coldspray nozzles. A remark can be made about the geometries mentioned in the Tab 2: even if 𝜉 𝑖 , 𝜁 and 𝐿 seem quite variable, 𝐷 and 𝜉 𝑜 are relatively stable and suggest either a tacit agreement of the authors on the experimental conditions of applicationof cold spray, or a lack of investigation of different geometries. Also, the dominant shape of the nozzles is a linear one. Thereal added value for the application of cold spray still seems to be unknown. The main objective of this paper is to reach a new step in numerical simulations comparing its fidelity with previous models.This kind of comparison has never been done before in the field which adds originality to this work. Even if the wholeproblem remains complex, it is now possible to advance further into the comprehension of the problem. The paper presentsfirst a commonly used model to validate the approach before increasing the fidelity with a new high fidelity model. This gainin fidelity will also allow to present a new type of configuration with an increased performance. Hence, the paper shows,in Section 2, the models with the hypothesis, the geometries, the boundary conditions and the CFD conditions. Then, inSection 3, the results of the two previous models will be discussed. Finally, in Section 4, a conclusion will be drawn accordingto which has been presented.
As mentioned in the previous section, we propose two approaches. The first is related to the common modeling of coldspray, referred as the RANS model. Detailed analysis and comparisons with (Ref [35]) will be presented for validations.5 o u r c e T u r b u l e n c e m o d e l T y p e o f s i m u l a t i o n G e o m e t r y P a r t i c l e - fl o w i n t e r a c t i o n T h e r m a l fl o w b o u n d a r y c o n d i t i o n G a s m o d e l V i s c o s i t y m o d e l O r d e r o f t h ee l e m e n t s S o l v e r P a r t i c l e t r a c k i n g s c h e m e D r a g c o e ffi c i e n t ( R e f [ ] ) 𝑘 − 𝜀 - D -- I -- D L - ( R e f [ ] ) 𝑘 − 𝜀 S t D A --- D L S N ( R e f [ ] ) 𝑘 − 𝜀 S t D - A I -- P -- ( R e f [ ] ) 𝑘 − 𝜀 / R S M S t D - A I C - P L - ( R e f [ ] ) 𝑘 − 𝜀 S t D - I C Q D L - ( R e f [ ] ) 𝑘 − 𝜀 T a n d S t D A I - F D L C l ( R e f [ ] ) L a m i n a r T D --- C - P L S N / W Y / N S / P M ( R e f [ ] ) L a m i n a r T D - I C - P E S N / G / D F / B / K H ( R e f [ ] ) L a m i n a r T D -- I S u - D -- ( R e f [ ] ) 𝑘 − 𝜀 / R S M S t D A I S u S e / Q P L - ( R e f [ ] ) 𝑘 − 𝜀 S t D - I - S e - L M ( R e f [ ] ) R S M - D - I S u S e : Q P L C r + G ( R e f [ ] ) 𝑘 − 𝜔 S t D ----- L B i ( R e f [ ] ) 𝑘 − 𝜀 S t D - I -- P L H ( R e f [ ] ) 𝑘 − 𝜀 S t D A I - S e P L - ( R e f [ ] ) 𝑘 − 𝜀 S t D A I - s e D L M C ( R e f [ ] ) 𝑘 − 𝜀 - D - I --- E C r ( R e f [ ] ) 𝑘 − 𝜀 - D -------- ( R e f [ ] ) 𝑘 − 𝜀 S t D - A I --- L - P r e s e n t w o r k D E S T D T e m p e r a t u r e fi x e d o u t s i d e L a w i n ( R e f [ ] ) S u T h i r d o r d e r D L M C + G + A M T a b l e : L i s t o f t h e d i ff e r e n t n u e m r i c a l p a r a m e t e r s u s e d b y t h e a u t h o r s . T h e " - " s i g n m e a n s t h e i n f o r m a t i o n w a s n o t a v a i l a b l e i n t h e a r t i c l e . “ R S M ” m e a n s R e y n o l d s S t r e ss M o d e l , “ S t ” m e a n ss t e a d y , “ T ” m e a n s t r a n s i e n t , “ ” m e a n s t w o - w a y c o u p l e d , “ ” m e a n s o n e - w a y c o u p l e d , “ A ” m e a n s a d i a b a t i c , “ I ” m e a n s i d e a l g a s l a w , “ C ” m e a n s c o n s t a n t , “ S u ” m e a n s S u t h e r l a n d ’ s l a w ( R e f [ ] ) , “ F ” m e a n s fi r s t o r d e r e l e m e n t s , “ S e ” m e a n ss e c o n d o r d e r e l e m e n t s , “ Q ” m e a n s Q U I C K ( a d i s c r e t i z a t i o n s c h e m e d e s c r i b e d i n ( R e f [ ] )) , “ P ” m e a n s p r e ss u r e - b a s e d s o l v e r , “ D ” m e a n s d e n s i t y - b a s e d s o l v e r , “ L ” m e a n s L a g r a n g i a n , “ E ” m e a n s E u l e r i a n , “ D E S ” m e a n s D e t a c h e d E dd y S i m u l a t i o n , “ S N ” m e a n s S c h i ll e r - N a u m a nn m o d e l ( s ee R e f [ ] ) , “ G ” m e a n s G i d a s p o w m o d e l ( s ee ( R e f [ ] )) , “ D F ” m e a n s D i F e l i c e m o d e l ( s ee ( R e f [ ] )) , “ B ” m e a n s B ee s t r a m o d e l ( s ee ( R e f [ ] )) , “ K H ” m e a n s K o c h - H i ll m o d e l ( s ee ( R e f [ ] )) , “ W Y ” m e a n s W e n - Y u m o d e l ( s ee ( R e f [ ] )) , “ N S ” m e a n s N o n - S p h e r e m o d e l ( s ee ( R e f [ ] )) , “ P M ” m e a n s D u P l e ss i s a n d M a s l i y a h m o d e l ( s ee ( R e f [ ] )) , “ C r ” m e a n s C r o w e m o d e l ( s ee ( R e f [ ] )) , “ C l ” m e a n s C l i f t e t a l . m o d e l ( s ee ( R e f [ ] )) , “ M ” m e a n s M o r s i a n d A l e x a n d e r m o d e l ( s ee ( R e f [ ] )) , “ M C ” m e a n s M o r s i a n d A l e x a n d e r m o d e l w i t h c o m p r e ss i b i l i t y c o rr e c t i o n ( s ee ( R e f [ ] )) , “ B i ” m e a n s B i r d e t a l . m o d e l ( s ee ( R e f [ ] )) , “ G ” m e a n s g r a d i e n t f o r c e , “ H ” m e a n s H e n d e r s o n m o d e l ( s ee ( R e f [ ] )) , “ A D ” m e a n s a dd e d m a ss . 𝐿 (mm) 𝐷 (mm) 𝜉 𝑜 𝜉 𝑖 𝜁 Shape(Ref [3])
Linear(Ref [4]) - Linear(Ref [5]) - - Piecewise constant(Ref [6])
224 2.7 2.22 6.74 0.759
Linear(Ref [7])
Linear(Ref [8])
Linear(Ref [11])
140 2.7 2.37 ≈ 1.5 0.929
Linear(Ref [12])
137 2.55 1.88 3.14 0.956
Piecewise constant(Ref [28])
50 2
Linear(Ref [29]) - Linear - Linear(Ref [31])
120 2 2.0 5.0 0.933
Linear(Ref [32]) - - Linear(Ref [34])
70 2.0 2.0 5.0 0.571
Linear(Ref [35])
210 2 3.0 11 0.857
Linear
115 2.7 3.0 6.67 0.565
Linear
87 4 1.15 5.0 0.805
Linear
35 2.2 1.73 9.09 0.571
Linear(Ref [37])
LinearTable 2: List of the different nozzle models proposed by the authors. The "-" sign means the information was not availablein the article. When 𝜁 is not known, the length 𝐿 given is the length of the divergent. The "linear" shape means thatthe diameter changes linearly between the inlet and the neck and between the neck and the outlet. The form "piecewiseconstant" means that the diameter evolves by sections over which it is constant.7hen, we will provide a more complex model for the high fidelity simulations known as the Improved Delayed DetachedEddy Simulation, and referred here as the IDDES model. Assuming the gas as an ideal gas, neglecting any dissipation, i.e. 𝜆 = 0 , 𝜇 = 0 , 𝜅 = 0 , and assuming an isentropic flow with aone dimension hypothesis, according to (Ref [10, 11, 15, 43]), we can write = 𝑇 𝑇 = ( 𝑝 𝑝 ) 𝛾−1𝛾 = ( 𝜌 𝜌 ) 𝛾−1 , 𝑀 = ‖𝒖‖√𝛾 𝑟𝑇 (1)which are called the Barrï¿œ de Saint-Venant’s relations. The refers to the values imposed at the beginning of the nozzlein the chamber of stagnation. Hence, all the quantities are depending on the Mach number 𝑀 which is determined by thefollowing relation, ( 𝐷𝐷 𝑡 ) = 1𝑀 ( 2𝛾 + 1 (1 + 𝛾 − 12 𝑀 )) 𝛾+1𝛾−1 (2)where 𝐷 𝑡 is the diameter of the throat. Therefore, 𝑀 < 1 before the throat and
𝑀 > 1 after the throat. The Eq (1) and the Eq(2) will be used afterwards to compare and to discuss the predictions of both numerical models.
We consider an axisymmetric modeling of the nozzle neglecting gravity. The gas considered is nitrogen and the Stokeshypothesis is applied, i.e. 𝜅 = 0 . The dynamic viscosity 𝜇 evolves according to (Ref [21]) to agree with other papers (see Eq(3)), 𝜇 = 𝜇 ref ( 𝑇𝑇 ref ) 𝑇 ref + 𝑇 𝑆 𝑇 + 𝑇 𝑆 (3)where 𝜇 ref = 1.663 × 10 −5 kg.m −1 .s −1 , 𝑇 ref = 273.11 K and 𝑇 𝑆 = 106.67 K the Sutherland temperature of the model. Froma purely numerical point of view, the resolution method used is a density-based method which, according to (Ref [35]), issupposed to approach compressible phenomena such as shocks better than the pressure-based resolution method. On eachwall, the no-slip boundary condition is applied. The particles are tracked according to a Lagrangian description because oftheir low volume fraction in the flow. Their motion is tracked with a stochastic walk to account for turbulence effects. Thetemperature of the particles is also governed by Eq (4), d𝑇 𝑝 d𝑡 = 6𝜆𝜌 𝑝 𝑐 𝑝 𝑑 (𝑇 − 𝑇 𝑝 ) Nu (4)with Nu = 2.0 + 0.6Re Pr , Re 𝑝 = 𝜌𝑑 𝑝 ‖‖𝒖 − 𝒖 𝒑 ‖‖𝜇 (5)according to (Ref [33]). The size of the particles is uniform with 𝑑 𝑝 = 20 µ m . These particles are made of copper. We consider the stationary RANS (see (Ref [15, 17, p.12-9])) to model the fluid flow given by Eq (6) to Eq (8), d𝜌d𝑡 + 𝜌 𝜕𝑢 𝑖 𝜕𝑥 𝑖 = 0 (6) 𝜌 d𝑢 𝑖 d𝑡 = − 𝜕𝑝𝜕𝑥 𝑖 + 𝜕𝜕𝑥 𝑗 (𝜇 ( 𝜕𝑢 𝑖 𝜕𝑥 𝑗 + 𝜕𝑢 𝑗 𝜕𝑥 𝑖 − 23 𝜕𝑢 𝑘 𝜕𝑥 𝑘 𝛿 𝑖𝑗 ) + 𝑅 𝑖𝑗 ) , (7) 𝜌 d𝑒d𝑡 = −𝑢 𝑖 𝜕𝑝𝜕𝑥 𝑖 + 𝜕𝜕𝑥 𝑖 (𝜆 eff 𝜕𝑇𝜕𝑥 𝑖 ) + 𝜕𝜕𝑥 𝑖 ((𝜇 + 𝜇 𝑡 ) ( 𝜕𝑢 𝑖 𝜕𝑥 𝑗 + 𝜕𝑢 𝑗 𝜕𝑥 𝑖 − 23 𝜕𝑢 𝑘 𝜕𝑥 𝑘 𝛿 𝑖𝑗 ) 𝑢 𝑗 ) , (8)8here 𝑅 𝑖𝑗 = 𝜇 𝑡 ( 𝜕𝑢 𝑖 𝜕𝑥 𝑗 + 𝜕𝑢 𝑗 𝜕𝑥 𝑖 − 23 𝜕𝑢 𝑙 𝜕𝑥 𝑙 𝛿 𝑖𝑗 ) − 23 𝜌𝑘 𝜕𝑢 𝑙 𝜕𝑥 𝑙 𝛿 𝑖𝑗 (9)according to Boussinesq’s hypothesis (see (Ref [15, p.12-9])) and 𝜆 eff = 𝜆 + 𝑟𝛾 𝜇 𝑡 Pr 𝑡 (𝛾 − 1) . (10)The turbulence model used is the realizable 𝑘 − 𝜀 model which can be achieved by taking into account the effects of com-pressibility and modeling the gas as an ideal gas with Eq (11), 𝑝 = 𝜌𝑟𝑇 . (11)The turbulence model is presented in Eq (12) to Eq (19), 𝜇 𝑡 = 𝐶 𝜇 𝜌 𝑘 𝜀 , (12) 𝜌 d𝑘d𝑡 = 𝜕𝜕𝑥 𝑖 ((𝜇 + 𝜇 𝑡 𝜎 𝑘 ) 𝜕𝑘𝜕𝑥 𝑖 ) + 𝜇 𝑡 ̃𝑆 − 𝜌𝜀 − 2𝜌𝜀𝑘𝛾 𝑟𝑇 , (13) 𝜌 d𝜀d𝑡 = 𝜕𝜕𝑥 𝑖 ((𝜇 + 𝜇 𝑡 𝜎 𝜀 ) 𝜕𝜀𝜕𝑥 𝑖 ) + 𝜌𝐶 ̃𝑆𝜀 − 𝜌𝐶 𝜀 𝑘 + √𝜇𝜀/𝜌 , (14) Pr 𝑡 = 0.85, 𝐶 = max (0.43, 𝜏𝜏 + 5 ) , 𝜏 = ̃𝑆 𝑘𝜀 , (15) 𝐶 𝜇 = (𝐴 + 𝐴 𝑠 𝑘𝑈 ∗ 𝜀 ) −1 , 𝑈 ∗ = √𝑆 + Ω , (16) 𝐴 = 4.04, 𝐴 𝑠 = √6 cos ( 13 arccos (√6𝑊 )) , (17) 𝑊 = 𝑆 𝑖𝑗 𝑆 𝑗𝑘 𝑆 𝑘𝑖 𝑆 , ̃𝑆 = √2𝑆, 𝑆 𝑖𝑗 = 12 (𝜕 𝑗 𝑢 𝑖 + 𝜕 𝑖 𝑢 𝑗 ) , (18) 𝐶 = 1.9, 𝜎 𝑘 = 1.0, 𝜎 𝜀 = 1.2. (19)Heat transfers between the gas and the nozzle walls are neglected as well as heat transfers with the substrate. Also, thequantities sought are calculated with second-order elements. The particles do not have any influence on the flow and theirmotion is governed by Eq (20), d𝒖 𝒑 d𝑡 = 3𝜌4𝜌 𝑝 𝑑 𝑝 ‖‖𝒖 − 𝒖 𝒑 ‖‖ (𝒖 − 𝒖 𝒑 ) 𝐶 𝐷 (20)where 𝐶 𝐷 = (𝑎 + 𝑎 Re 𝑝 + 𝑎 Re ) 𝐶 −1 (21)according to (Ref [41, 42]) taking into account the compressible effects with 𝐶 = 1 + √ 𝜋𝛾2 𝑀 𝑝 Re 𝑝 (2.514 + 0.8 exp (−0.55√ 2𝜋𝛾 Re 𝑝 𝑀 𝑝 )) , (22) 𝑀 𝑝 = ‖‖𝒖 − 𝒖 𝒑 ‖‖𝑐 , 𝑐 = √𝛾 ( 𝜕𝑝𝜕𝜌 ) 𝑇 = √𝛾 𝑟𝑇 (23)with Eq (11), 𝛾 = 1 − 𝑇𝜌 ( 𝜕𝑒𝜕𝑇 ) −1𝜌 ( 𝜕𝑝𝜕𝑇 ) 𝜌 ( 𝜕𝜌𝜕𝑇 ) 𝑝 = 75 (24)with Eq (11) for nitrogen and (𝑎 𝑖 ) 𝑖∈ (cid:74) (cid:75) given by Tab 3. The particles do not exchange any heat with the nozzle walls or thesubstrate. 9 e 𝑝 𝑎 𝑎 𝑎 [0; 0.1] [0.1; 1] [1; 10] [10; 100] [100; 1, 000] [1, 000; 5, 000] [5, 000; 10, 000] [10, 000; +∞[ (𝑎 𝑖 ) 𝑖∈ (cid:74) (cid:75) according to Ref [42] and Ref [17] In this section, we present the details about the new proposed model which aims to increase the fidelity of the simulationwith more representative assumptions and governing equations.Therefore, we propose the implementation of the Improved Delayed Detached Eddy Simulation (or IDDES) to model thefluid flow. The IDDES equations are given by Eq (25) to Eq (44)(see (Ref [17, 44, 45])), 𝜌 d𝑘d𝑡 = 𝜕𝜕𝑥 𝑖 ((𝜇 + 𝜎 ′𝑘 𝜇 𝑡 ) 𝜕𝑘𝜕𝑥 𝑖 ) + 𝑃 𝑘 − 𝜌𝑘 𝑙 IDDES (25) 𝜌 d𝜔d𝑡 = 𝜕𝜕𝑥 𝑖 ((𝜇 + 𝜎 𝜔 𝜇 𝑡 ) 𝜕𝜔𝜕𝑥 𝑖 ) + (1 − 𝐹 ) 2𝜌𝜎 𝜔2 𝜕 𝑖 𝑘𝜕 𝑖 𝜔𝜔 + 𝛼𝜌𝜇 𝑡 𝑃 𝑘 − 𝛽𝜌𝜔 (26)where 𝜇 𝑡 = 𝜌𝐴 𝑘max (𝐴 𝜔, 𝐹 𝑆) , 𝐹 = tanh (arg ) (27) arg = min (max ( √𝑘𝐶 𝜇 𝜔𝑑 𝑤 , 500𝜇𝑑 𝜔𝜌 ) , 4𝜌𝜎 𝜔2 𝑘𝐶𝐷 𝑘𝜔 𝑑 ) (28) 𝐶𝐷 𝑘𝜔 = max ( 2𝜌𝜎 𝜔2 𝜕 𝑖 𝑘𝜕 𝑖 𝜔𝜔 , 10 −10 ) , (29) 𝛼 = 59 𝐹 + 0.44 (1 − 𝐹 ) , 𝛽 = 0.075𝐹 + 0.0828 (1 − 𝐹 ) , (30) 𝜎 ′𝑘 = 0.85𝐹 + 1 − 𝐹 , 𝜎 𝜔 = 0.5𝐹 + 𝜎 𝜔2 (1 − 𝐹 ) , (31) 𝜎 𝜔2 = 0.856, 𝐹 = tanh (max ( 2√𝑘𝐶 𝜇 𝜔𝑑 𝑤 , 500𝜇𝑑 𝜔𝜌 ) ) , (32) 𝑃 𝑘 = min (𝜇 𝑡 𝑆 , 10𝜌𝑘𝜔) , (33) 𝑙 IDDES = 𝑓 𝑑 (1 + 𝑓 𝑒 ) 𝑙 RANS + (1 − 𝑓 𝑑 ) 𝑙 LES , (34) 𝑙 LES = 𝐶
DES min (𝐶 𝑤 max (𝑑 𝑤 , ℎ max ) , ℎ max ) , (35) 𝑙 RANS = √𝑘𝐶 𝜇 𝜔 , 𝐶 DES = 𝐶
DES1 𝐹 + 𝐶 DES2 (1 − 𝐹 ) , (36) 𝑓 𝑑 = max (1 − 𝑓 𝑑𝑡 , 𝑓 𝑏 ) , 𝑓 𝑑𝑡 = 1 − tanh ((𝐶 𝑑𝑡1 𝑟 𝑑𝑡 ) 𝐶 𝑑𝑡2 ) , (37) 𝑟 𝑑𝑡 = 𝜇 𝑡 𝜘 𝜌𝑑 𝑤 √(𝑆 + Ω ) /2 , (38)10 𝑏 = min (2e −9𝛼 ′2 , 1.0) , 𝛼 ′ = 14 − 𝑑 𝑤 ℎ max , (39) 𝑓 𝑒 = 𝑓 𝑒2 max (𝑓 𝑒1 − 1.0, 0.0) , 𝑓 𝑒1 = {2e −11.09𝛼 ′2 𝛼 ′ ≥ 02e −9.0𝛼 ′2 𝛼 ′ < 0 , (40) 𝑓 𝑒2 = 1.0 − tanh (max ((𝐶 𝑟 𝑑𝑡 ) , (𝐶 𝑟 𝑑𝑙 ) )) , (41) 𝑟 𝑑𝑙 = 𝜇𝜘 𝜌𝑑 𝑤 √(𝑆 + Ω ) /2 , 𝐶 𝑤 = 0.15, 𝐶 𝑑𝑡1 = 20, (42) 𝐶 𝑑𝑡2 = 3, 𝐶 𝑙 = 5.0, 𝐶 𝑡 = 1.87, 𝐶 𝜇 = 0.09, 𝜘 = 0.41, (43) 𝐴 = 0.31, 𝐶 DES1 = 0.78, 𝐶
DES2 = 0.61. (44)This model of turbulence is designed to take advantage of both Large Eddy Simulation (or LES) and RANS simulation. TheEq (34) and the Eq (25) show length scales which are coming from either LES or RANS simulation. In the detached areas,the model calculates completely the flow until the smallest local cell size and simulates the lower scales. Otherwise, nearthe walls, to properly get the boundary layer and its effects, the RANS model takes the dominance. The whole IDDES modelpurpose is to manage the transition between the detached areas and the boundary layer. We recall that the steady statesolution of the previously described RANS model can be used as an initialization for the IDDES model. The gas is modeledwith the Redlich-Kwong-Aungier’s law (see (Ref [16, 17]) ) given by Eq (45) to Eq (48) 𝑝 = 𝜌𝑟𝑇1 + 𝜌 (𝑐 − 𝑏) − 𝛼 𝜌 𝑐 ) −𝑛 (45)where 𝑐 = 𝑟𝑇 𝑐 (𝑝 𝑐 + 𝛼 𝜌 𝑐 𝑏 ) −1 + 𝑏 − 1𝜌 𝑐 , (46) 𝛼 = 0.42747𝑟 𝑇 𝑝 𝑐 , 𝑏 = 0.08664𝑟𝑇 𝑐 𝑝 𝑐 , (47) 𝑛 = 0.4986 + 1.1735𝜐 + 0.4754𝜐 . (48)The wall of the nozzle and the substrate are modeled as thick steel walls with an external temperature of
300 K . Also,the quantities sought are calculated with third-order elements. The particles have an influence on the flow — a two-waycoupled system — and their motion is governed by Eq (49), (1 + 𝐶 vm 𝜌𝜌 𝑝 ) d𝒖 𝒑 d𝑡 = 3𝜌4𝜌 𝑝 𝑑 𝑝 ‖‖𝒖 − 𝒖 𝒑 ‖‖ (𝒖 − 𝒖 𝒑 ) 𝐶 𝐷 + (1 + 𝐶 vm ) 𝜌𝜌 𝑝 𝒖 𝒑 div (𝒖) (49)where 𝐶 vm = 0.5 is the virtual mass factor assuming spherical particles (see (Ref [46])) and 𝐶 𝐷 is defined by Eq (21) with 𝑐 given by Eq (23) with the new state law Eq (45). The geometries used in this paper are given on Fig 2. The dimensions of the various configurations are given by Tab 4. 𝜂 iskept between the four configurations. The outlet boundaries are set as far as possible from the nozzle exit to set the pressureto the atmospheric one and the temperature to the ambient one. The configurations A to C correspond to three configurationsproposed in Ref [35]. However, the last case D corresponds to a new configuration keeping the same characteristics as Abut changing the shape of the channel in the convergent. The boundary conditions, sketched on Fig 3, are summarized in Tab 4. It is interesting to point out that the exact position ofthe injector when placed in the stagnation does not matter because the major acceleration occurs in the divergent part. Forconfigurations A, B and D, the particles enter into the nozzle at the inlet where the pressure and the temperature are fixed.For configuration C, the particles have their own injector and start at the inlet of the injector where the pressure and thetemperature are fixed. The mass flow rate of the particles is fixed for each configuration at −1 to agree with commonvalues found in the literature. 11 a) Configuration A(b) Configuration B(c) Configuration C(d) Configuration D. The convex part in the convergent is visible on the left. Figure 2: Geometries of the four configurations. The scale is .Figure 3: Sketch of the boundary conditions12onfiguration A B C D
𝐿(mm) 210 115 35 210𝐷 (mm) 2 2.7 2.2 2𝜁 0.857 0.565 0.571 0.857𝜉 𝑖 𝑜
11 6.67 9.09 11𝜂 0.3 0.3 0.3 0.3
Injector diameter (mm)
Injector position (mm) S S +0.8
SInjector pressure (bar)
Standoff distance (mm)
40 100 70 40
Inlet pressure (bar)
30 30 20 30
Outlet pressure (bar)
Inlet temperature (K)
300 300 300 300
Outlet temperature (K)
300 300 300 300
Shape Linear Linear Linear Convex in the convergentand linear in the divergentTable 4: Dimensions and boundary conditions of the configurations. “S” means in the stagnation chamber and +𝑥 meansafter the throat of 𝑥mm . The configuration A to C corresponds to the configurations 1, 2 and 4 presented by Ref [35] The mesh, depicted on Fig 4, is 1,000 times wider in the buffer area to gain in computation time and refined along the outsidethe nozzle to be more accurate around the shocks and the substrate (the factor compared to the inside of the nozzle is around4 times smaller). Attempts were made with a refined mesh in the inside of the nozzle but there is not a significant differencein terms of the performance of the nozzle. Moreover, the mesh is refined near the nozzle, which is not very relevant forthe RANS model but useful for the IDDES model. Therefore, there are between 20,000 and 35,000 elements according to theconfiguration.
The results to compare with those of (Ref [35]) are presented in this subsection. Fig 5 presents the obtained velocity mag-nitude of configuration A to C. As expected, we note a strong agreement between the results and we show that the RANSsimulations have succeeded to well capture the shocks in the flow, the global acceleration of the flow and its shape. Thevelocity magnitude increases progressively from the throat to the exit until the first shock.A close comparison between both configurations A shows approximately 1% of error. Nevertheless, even if the RANSmodel has tried to mimic the model of (Ref [35]), there are several approximations which are important to take into accountbecause they create differences with real experiments. Looking at the model purely, the gas law is a first step to analyzebecause if the Van der Waals’ gas law is used (see (Ref [13])), comparing the density given by the ideal gas law 𝜌 ig and thedensity given by the Van der Waals’ law 𝜌 VdW , we get approximately |||𝜌
VdW − 𝜌 ig |||𝜌 VdW ≈ 𝑝𝑇 𝑐 𝑐 𝑇 ||||| 27𝑇 𝑐
8𝑇 − (1 − 𝑝𝑇 𝑐 𝑐 𝑇 ) −1 ||||| = 1.7% (50)with 𝑝 = 30 bar , 𝑇 = 300 K and for nitrogen, 𝑝 𝑐 = 34 bar and 𝑇 𝑐 = 126 K . Therefore, looking for the lowest error between the13 a) Global mesh(b) Zoom in the convergent (c) Zoom in the divergent (d) Zoom in the exit Figure 4: Screenshots of the mesh used (a) Configuration A(b) Configuration B(c) Configuration C
Figure 5: Velocity magnitude of the flow between the exit of the nozzle and the substrate for the configuration A to C. Thescale is . 14 a) RANS model(b) (Ref [35]) results
Figure 6: Velocity magnitude of the flow from the exit to the substrate for configuration A with RANS model comparingwith the results in (Ref [35])simulations is not relevant because even a calculation with orders of magnitude gives slight differences. It is also particularlyrelevant to compare the flow without adding particles because those ones do not have influence on the flow.
In this sub section, we focus on plotting the turbulent kinetic energy for configuration A. It is in good agreement with (Ref[35]). The shape are rather similar between both results. Precisely, there is a peak of turbulent kinetic energy just afterthe throat where the flow becomes supersonic and where shocks may appear. This increase of turbulent kinetic energy istherefore an increase of turbulent viscosity which enhance dissipation from turbulent effects. It is also a sign there is a lotof mixing around the throat and confirms the fact that the particle beam would completely fill the nozzle divergent duringthe particles fly. To analyze furthermore the behavior of these configurations, the Fig 7 compares configurations A to C interms of turbulent kinetic energy 𝑘 .The U shape in the convergent is clearly visible on Fig 7. Near the throat, it starts to increase because of the supersonicflow beginning. Then, in the diverging part, it is approximately constant, demonstrating that there is no additional phenom-ena to dissipate like shocks or something else. For each configuration, after the diverging part, there is a brutal increase ofturbulent kinetic energy. Using the Fig 5, these brutal increases occur on the location of the first oblique shock. Continuingalong the axis of the nozzle, the fluctuations of the turbulent kinetic energy occur exactly where there is a shock. Finally,near the substrate, there is a strong increase of the turbulent kinetic energy caused by the bow shock at the impact. Consid-ering the values of 𝑘 , this bow shock has a important effect on the dissipation and on turbulent effect. It has therefore a nonnegligible influence on the particles impact velocity because they are slowed down by this discontinuity. It confirms thenthe conclusion drawn in the literature that lighter particles will be easily slowed down in comparison to heavier particleswith bigger inertia. However, these heavier ones gain less speed during their acceleration thanks to the same reason. Anoptimum must be found to reach the highest impact velocity with a given configuration. Comparing configurations A to C,it happens configuration B has more fluctuations along the flow which can be problematic for the nozzle in use and all themachinery nearby. We would then prefer to use configuration C which has less steep variations of turbulent kinetic energyto avoid breaking the material in use. Pursuing the comparison of the results, Fig 8 presents the velocity magnitude of the particles during the fly outside the nozzlefor configuration A . The shapes of all the particle tracks is conserved between the RANS model results and those in (Ref[35]). In terms of values, the error reaches a maximum among all the configurations of 4.5% which is pretty good knowingall the reasons mentioned above about the model, the mesh, the approximations, etc. The distribution of the velocities acrossthe particle beam looks rather similar between the two types of results. The behaviour is similar between the configurations . . . . . . . − k ( m / s ) Configuration AConfiguration BConfiguration C
Figure 7: 𝑘 as a function of the dimensionless position along the nozzle on the axis of symmetry from RANS model resultson configuratons A to C. The dimensionless position corresponds to the position divided by 𝐿 . Dimensionless position 0corresponds to the nozzle throat (a) RANS model(b) (Ref [35]) results Figure 8: Particles track colored according to their velocity magnitude between the exit of the nozzle and the substrate forconfiguration A. The scale is . 16 a) RANS model(b) (Ref [35]) results
Figure 9: Particles track colored according to their velocity magnitude between the exit of the nozzle and the substrate forconfiguration B. The scale is (a) RANS model(b) (Ref [35]) results
Figure 10: Particles track colored according to their velocity magnitude between the exit of the nozzle and the substrate forconfiguration C. The scale is I m p ac t v e l o c it y ( m / s ) Lupoi’s trend lineRANS model trend lineLupoi’s CFDRANS model CFD
Figure 11: Particle impact velocity as a function of the distance from the axis of symmetry for configuration A with RANSmodel
Finally, it is possible to compare the particle impact velocities against their distance from the axis of symmetry. All thesesresults are presented in Fig 11. Again, the results look very similar between the configurations. The maximum velocityin each case of the RANS model is close to the maximum velocity of the results in (Ref [35]) with an error close to 1%.Moreover, the distribution of particles are identical in terms of proportion; that is to say the configurations A and B havea denser proportion of particles near the center for both results and the configuration C has a more uniform distributionof particles in both results. The width of the band in configuration A is precisely respected in both results. Besides, forconfiguration B, the width seems a little narrower whereas configuration C has a wider width. Configuration C does notseem to be problematic because the authors of (Ref [35]) carried out experiment with each configuration and showed thatthe width of the band for configuration C is experimentally wider than its predictions. The RANS model results are closerto the experimental results than those in (Ref [35]). About configuration B, the RANS model results show -wide bandinstead of the -wide one in (Ref [35]). Hence, the RANS model disagree with (Ref [35]) on configuration B for particleimpact velocity against their distribution.
In this section, we present all the new results obtained using the proposed IDDES model. For comparisons, the previousRANS model will be used in most of the cases.
We recall that the IDDES simulation is transient. Therefore, Fig 12 shows the comparison of the velocity magnitude betweenthe RANS model and the IDDES model at different time steps for configuration A . Globally, thanks to the increased accuracyfor the space discretization and the improved turbulent model, the results in this case seem less diffusive and capture preciselythe flow.The slow radial decrease of velocity around the substrate is now a clear difference, like a jump. It is also visible that asteady state of the flow could be approximately reached. We studied a detailed comparison of all the previous models with the plots of 𝑀 , 𝑝/𝑝 and 𝑇 /𝑇 .The general behavior of the plots are rather similar. Even if the theoretical model needs stronger assumptions, it givesrather relevant orders magnitude for the different variables under study. Furthermore, on each configuration and eachvariable, the range of evolution is wider for the theoretical model presented in 2.1.1 because we assumed no dissipation byturbulence, shocks or thermal conduction. Thus, the flow is allowed to gain more velocity with a higher 𝑀 lowering thendrastically 𝑝 and 𝑇 . Between the RANS model and the IDDES model, the latter tends to lead to more dissipation with a lower The behaviour is similar between the configurations a ) C o n fi g u r a t i o n A a t 𝑡 = m s ( b ) C o n fi g u r a t i o n A a t 𝑡 = m s ( c ) C o n fi g u r a t i o n A a t 𝑡 = m s ( d ) C o n fi g u r a t i o n B a t 𝑡 = m s ( e ) C o n fi g u r a t i o n B a t 𝑡 = m s ( f ) C o n fi g u r a t i o n B a t 𝑡 = m s ( g ) C o n fi g u r a t i o n C a t 𝑡 = m s ( h ) C o n fi g u r a t i o n C a t 𝑡 = m s ( i ) C o n fi g u r a t i o n C a t 𝑡 = m s F i g u r e : V e l o c i t y m a g n i t u d e f o r c o n fi g u r a t i o n A t o C a t d i ff e r e n tt i m e s t e p s f r o m t h ee x i t o f t h e n o zz l e t o t h e s u b s t r a t e w i t h t h e R A N S m o d e l a n d t h e I DD E S m o d e l . T h e b e h a v i o r i n s i d e t h e n o zz l e i s r a t h e r s i m i l a r b e t w ee n t h e m o d e l s . T h e I DD E S m o d e l i s i n t h e u pp e r h a l f p a r t o f e a c h p i c t u r e a n d t h e R A N S m o d e l i s i n t h e l o w e r h a l f p a r t o f e a c h p i c t u r e . 𝑀 , 𝑝 and 𝑇 . It can be surprising with what has been said previously but the turbulent effects are morecontrolled by the IDDES model and allow a more accurate simulation. Compared to the RANS model, the IDDES model allows to catch the phenomena occurring near the substrate like the bowshock. This shock, visible in Fig 13 in terms of pressure , creates a huge deceleration of the flow which impacts directlythe particles coming onto the substrate. Because the IDDES model is less dissipative, the bow shock near the substrate isstrong and has a more important effect on the small particles because of their lower inertia. The shape of the bow shockdepends on the topology of the whole flow as presented by the Fig 13, but the global phenomenon is pretty similar betweenconfiguration A to C. One of the steps to compare both models is to compare 𝑘 on the axis of symmetry as shown in Fig 14. Talking only aboutthe differences, 𝑘 obtained with the IDDES is clearly lower than 𝑘 got by the RANS model. The lines follow approximatelythe same path that is to say there are constant in the diverging part and increase suddenly where there is a shock. Furthermore, we propose in Fig 15 to compare the pressure at the intersection of the axis of symmetry and the substratebetween the two models. Generally, the pressure obtained with the IDDES model is higher than the one resulting from theRANS model except for configuration C. It must be related to the less dissipative turbulent effects modeled by the IDDESmodel compared to the RANS model. Also, the fluctuations due to the turbulence is clearly visible of the solid lines. Thesefluctuations oscillate around an average value and there is a change in the evolution of the pressure from 𝑡 = 20 ms forconfiguration A and 𝑡 = 12 ms for configuration B. This change is related to the arrival of huge amount of particles. In fact,the first particles start to fly at 𝑡 = 0 ms and must travel along the nozzle before reaching the substrate. Therefore, becausethere is a two-way coupling between the flow and the particles, the effect of the particles on the turbulence appears whenthere are enough particles colliding with the substrate.
Fig 16 shows the particles coloured according to their velocity magnitude. It is visible that the particles have not reached yetthe substrate at 𝑡 = 10 ms for configuration A and B which was confirmed by Fig 15 with the fluctuations of pressure on thesubstrate. The configuration C has already particles onto the substrate at 𝑡 = 10 ms as the fluctuations in Fig 15 have shown.The particles are quite organized in the convergent but tend to collide and to deviate from the throat. From this point, theturbulent effects are strong and produce a lot of mixing in the particles, as expected. The distribution of the particles alongthe radial position is then strongly dependent on the path of each particle which will be discussed afterwards.
Finally, to compare the two models, we can present the Fig 17 showing the particle impact velocity as a function of thedistance from the axis of symmetry according to the IDDES model and the RANS model. On configuration A, the impacts ofthe particles are really similar in terms of distribution of velocity and radial position. Thanks to the transient simulation, theIDDES model has succeeded to capture accurately the distribution of particles. On configuration B, which slightly differsfrom the results of (Ref [35]), the impacts of particles are rather different. One of the reasons can be the strong presence ofthe bow shock near the substrate as illustrated by Fig 13. Inside the bow shock, the pressure and the density of the fluidare much higher and the bow shock is deeper near the center. Hence, the particles close to the center spend more time ina fluid which decelerates them and have a lower impact velocity than those far from the center with a thinner bow shock.Compared to the RANS model which does not show any bow shock for configuration B, the distribution of velocity is lessuniform but the width of the band is kept thanks to the similar flow topology between both models. On configuration C,this discussion is close to the one for configuration B. The presence of the bow shock in the IDDES model creates a gradientof impact velocity from the center because the bow shock gets deeper when approaching the center. The width of the trackseems higher because of the strong turbulent effects in the flow. Nevertheless, the IDDES simulation has succeeded to catchthe width of the track given by the experimental results of (Ref [35]) which gives more credit again to the IDDES modelagainst the RANS model. the same behavior exists for density with the same relative orders of magnitude a ) C o n fi g u r a t i o n A a t 𝑡 = m s ( b ) C o n fi g u r a t i o n A a t 𝑡 = m s ( c ) C o n fi g u r a t i o n A a t 𝑡 = m s ( d ) C o n fi g u r a t i o n B a t 𝑡 = m s ( e ) C o n fi g u r a t i o n B a t 𝑡 = m s ( f ) C o n fi g u r a t i o n B a t 𝑡 = m s ( g ) C o n fi g u r a t i o n C a t 𝑡 = m s ( h ) C o n fi g u r a t i o n C a t 𝑡 = m s ( i ) C o n fi g u r a t i o n C a t 𝑡 = m s F i g u r e : P r e ss u r e f o r c o n fi g u r a t i o n A t o C w i t h t h e I DD E S m o d e l a t d i ff e r e n tt i m e s t e p s n e a r t h e s u b s t r a t e . T h e s u b s t r a t e i s a t b o tt o m o f e a c h i m a g e a n d t h e a x i a l p o s i t i o n c o rr e s p o n d s t o t h e d i s t a n c e f r o m t h e s u b s t r a t e . . . . . . . . . − − k ( m / s ) Configuration AConfiguration BConfiguration C
Figure 14: 𝑘 as a function of the dimensionless position along the nozzle on the axis of symmetry from IDDES modelon configuratons A to C. The dimensionless position corresponds to the position divided by 𝐿 . Dimensionless position 0corresponds to the nozzle throat . . . . P r e ss u r e ( M P a ) Configuration AConfiguration B Configuration C
Figure 15: Pressure as a function of time at the intersection of the axis of symmetry and the substrate. The solid curvescorrespond to the IDDES model and the dotted lines correspond to the RANS value constant along time.22 a ) C o n fi g u r a t i o n A a t 𝑡 = m s ( b ) C o n fi g u r a t i o n A a t 𝑡 = m s ( c ) C o n fi g u r a t i o n A a t 𝑡 = m s ( d ) C o n fi g u r a t i o n B a t 𝑡 = m s ( e ) C o n fi g u r a t i o n B a t 𝑡 = m s ( f ) C o n fi g u r a t i o n B a t 𝑡 = m s ( g ) C o n fi g u r a t i o n C a t 𝑡 = m s ( h ) C o n fi g u r a t i o n C a t 𝑡 = m s ( i ) C o n fi g u r a t i o n C a t 𝑡 = m s F i g u r e : P a r t i c l e v e l o c i t y m a g n i t u d e f o r c o n fi g u r a t i o n A t o C w i t h I DD E S m o d e l a t d i ff e r e n tt i m e s t e p s . I m p ac t v e l o c it y ( m / s ) IDDES modelRANS model
Figure 17: Particle impact velocity as a function of the distance from the axis of symmetry for configuration A with IDDESmodel and RANS modelFigure 18: Particle tracks colored according to their velocity magnitude between the exit of the nozzle and the substrate forconfiguration D with RANS model. The scale is .As a partial conclusion, the IDDES model is confirmed to be a more accurate simulation than the RANS model whichgives more insight of the phenomena appearing in cold spray such as turbulence, oblique shocks, bow shocks, fluctuations,particles motion and particles impacts. All these results can be used now to assess the performances of the last new proposedconfiguration.
While relying on the obtained results with the configurations A to C modeled first with the RANS model and then with theIDDES model, it is now possible to propose and to study the new configuration D. The topology of the flow is rather similarin terms of velocity magnitude, series of shocks, bow shock near the substrate, etc. Nevertheless, there is a strong differencefor the particles and their impact on the substrate which are going to be discussed in the following.
Fig 18 shows the particle tracks coloured according to the velocity magnitude of the particles with RANS model. Comparingthis results with Fig 8, the velocity magnitude of the particles seems rather similar between both pictures. However, the widthof the particle jet is much narrower. Thanks to the convex convergent, the particles did not suffer of the strong turbulenceoccurring at the throat. The convergent has succeeded to maintain the jet as narrow as possible. Also, because the particlesremain in the center of the flow, they take advantage of the maximum flow velocity to accelerate. Then, considering onlythe RANS model, configuration D seems to be a good improvement to carry on accurate cold spraying on the substrate.Fig 19 shows the particles coloured according to their velocity magnitude at different time steps obtained using the IDDESmodel. Again, comparing Fig 8 and Fig 19, the behavior remains the same but the width of the tracks looks narrower forconfiguration D. The orders of magnitude for the velocity are pretty similar and do not change between the configurationsor the models. 24 a) 𝑡 = 10 ms (b) 𝑡 = 20 ms (c) 𝑡 = 30 ms Figure 19: Particle velocity magnitude for configuration D with IDDES model at different time steps.25 . . . . . . . . I m p ac t v e l o c it y ( m / s ) Trend line Configuration DTrend line Configuration ARANS model Configuration DRANS model Configuration A
Figure 20: Particle impact velocity as a function of the distance from the axis of symmetry for configuration A and D withRANS model
To carry on an accurate comparison of the configurations A and D, it is now interesting to plot the evolution of the particleimpact velocity versus the distance from the axis of symmetry. These data are given in Fig 20 for RANS model. It is clearthat the particles coming from configuration D have a higher impact velocity than those of configuration A with RANSmodel with a much narrower width of the track on the substrate. The comments given previously are verified here thusthe configuration D gives improved results for the impact velocity magnitude and the width of the track according to RANSmodel. However, the particles are not distributed as usual: the particles are more present near from the center thanin the center. The flow creates a hole of particles in the center which can be problematic under certain conditions. It isimportant to point out that these conclusions are valid in the context of these simulations. Hence, using a real powderwith different sizes of particles or adding some fluctuations in the model will change drastically the results. That is why acomparison with the IDDES model helps to understand.Therefore, in Fig 21, the evolution of the particle impact velocity versus the distance from the axis of symmetry is givenwith the RANS model and the IDDES model. The particle impact velocity is a little bit lower with the IDDES model andreaches the same values as in (Ref [35]). Nevertheless, the width of the track and the distribution of particles is kept withthe RANS model. There are more particles in the center than in the RANS model which avoid the hole presented before.Some particles appear outside the track due to turbulent effects but the main track is much narrower than the one given byconfiguration A. With the IDDES model, because a cold spray system is used with a certain motion to create a coating, thesmall number of particles outside the thin track will not appear significantly and the user will get a narrower track withsimilar performance in terms of solidity, porosity or resistance.Both models conclude that configuration D has an improved performance against configuration A in terms of width oftrack. Because the particle impact velocities are kept, the adhesion of the particles remains the same but with a narrowertrack that is to say a more accurate coating system.
In this paper, a new high fidelity modeling for the cold spray process is presented. A review on several existing methods isalso proposed showing the capability of the new model by giving more insight of the phenomena appearing in cold spraysuch as turbulence, oblique shocks, bow shocks, fluctuations, particles motion and particles impacts. Several test cases andcomparisons were presented and discussed. The improvements of this proposed model are clearly visible on several fieldsof the problem from the topology of the flow, the fluctuations in time to the two-way coupled influence of the particles.These results underline the potential of this approach, and shall be pursued for new configurations. The extension of thecurrent work to three-dimensional can also be considered, using parallel computing and mesh adaptation. Also, exploitinga Large Eddy Simulation (LES) may lead to even better results.26 . . . . . . I m p ac t v e l o c it y ( m / s ) IDDES modelRANS model
Figure 21: Particle impact velocity as a function of the distance from the axis of symmetry for configuration D for IDDESmodel and RANS model
References [1] Wikipedia. Cold spraying — Wikipedia, the free encyclopedia, 2020. URL https://en.wikipedia.org/wiki/Cold_spraying . [Online; available on june-25-2020][2] A. Papyrin, V. Kosarev, S. Klinkov, A. Alkhimov, and V. Fomin.
Cold spray technology . Elsevier Science, The BoulevardLangford Lane. Kidlington Oxford 0X5 1CB UKRadarwcg 29 PO Box 21 L 10CO AE Amsterdam The Netherlands,première ed., 2007. ISBN 978-0-08045155-8[3] A. S. Alhulaifi, G. A. Buck, and W. J. Abegast. Numerical and Experimental Investigation of Cold Spray Gas DynamicEffects for Polymer Coating.
J. Therm. Spray Technol. , 2012. (5), pp. 852–862[4] M. Faizan-Ur-Rab, S.-H. Zahiri, S.-H. Masood, T.-D. Phan, M. Jahedi, and R. Nagarajah. Application of a holistic 3Dmodel to estimate state of cold spray titanium particles. Mater. Des. , 2015. , pp. 1227–1241. doi:10.1016/j.matdes.2015.10.075[5] D. MacDonald, S. Leblanc-Robert, R. Fernández, A. Farjam, and B. Jodoin. Effect of Nozzle Material on Down-stream Lateral Injection Cold Spray Performance. J. Therm. Spray Technol. , 2016. (6), pp. 1149–1157. doi:10.1007/s11666-016-0426-4[6] C. Singhal, Q. Murtaza, and Parvej. Simulation of Critical Velocity of Cold Spray Process with Different TurbulenceModels. Mater. Today , 5. pp. 17371–17379[7] S. Yin, Q. Liu, H. Liao, and X. Wang. Effect of injection pressure on particle acceleration, dispersion and deposition incold spray.
Comput. Mater. Sci. , 2014. , pp. 7–15. doi:10.1016/j.commatsci.2014.03.055[8] P. Liebersbach, A. Foelsche, V. K. Champagne, M. Siopis, A. Nardi, and D. P. Schmidt. CFD Simulations of Feeder TubePressure Oscillations and Prediction of Clogging in Cold Spray Nozzles. J. Therm. Spray Technol. , 2020. , pp. 400–412.doi:10.1007/s11666-020-00992-0[9] S. Yin, M. Meyer, W. Li, H. Liao, and R. Lupoi. Gas Flow, Particle Acceleration, and Heat Transfer in Cold Spray: Areview. J. Therm. Spray Technol. , 2016. (5), pp. 874–896. doi:10.1007/s11666-016-0406-8[10] K. H. Leitz, M. O’Sullivan, A. Plankensteiner, H. Kestler, and L. S. Sigl. OpenFOAM Modeling of Particle Heating andAcceleration in Cold Spraying. J. Therm. Spray Technol. , 2017[11] K.-H. Leitz, M. O’Sullivan, A. Plankensteiner, T. Lichtenegger, S. Pirker, H. Kestler, and L.-S. Sigl. CFDEM modellingof particle heating and acceleration in cold spraying.
Int. J. Refract. Met. Hard Mater , 2018. , pp. 192–198. doi:10.1016/j.ijrmhm.2018.02.003 2712] R. N. Raoelison, L. L. Koithara, S. Costil, and C. Langlade. Turbulences of the supersonic gas flow during cold sprayingand their negative effects: A DNS CFD analysis coupled with experimental observation and laser impulse high-speedshadowgraphs of the particles in-flight flow. Int. J. Heat Mass Transfer , 2020. , pp. 118894–118912. doi:10.1016/j.ijheatmasstransfer.2019.118894[13] D. Berthelot. Sur la notion des états correspondants et sur divers points correspondants remarquables (About thenotion of corresponding states and various noticeable corresponding points).
Journal de Physique Théorique et Ap-pliquée , 1903. (1), pp. 186–202. doi:10.1051/jphystap:019030020018600. URL hal.archives-ouvertes.fr/jpa-00240740/document . (in French)[14] O. Redlich and J. N. S. Kwong. On the Thermodynamics of Solutions. V. An Equation of State. Fugacities of GaseousSolutions. the Symposium on Thermodynamics and Molecular Structure of Solutions , vol. 44. Division of Physical andInorganic Chemistry, the 114th Meeting of the American Chemical Society, Portland, Oregon, pp. 233–244. doi:10.1021/cr60137a013[15] R. W. Johnson. Handbook of fluid dynamics . Electrical engineering handbook series. CRC Press, CRC PressTaylor andFrancis Group6000 Broken Sound Parkway NW, Suite 300Boca Raton, FL 33487-2742, 2nd ed., 2016. ISBN 978-1-4398-4957-6. URL https://books.google.fr/books?id=JBTlucgGdegC [16] R. H. Aungier. A Fast, Accurate Real Gas Equation of State for Fluid Dynamic Analysis Applications.
J. Fluids Eng. ,1995. (2), pp. 277–281. doi:10.1115/1.2817141[17] ANSYS Fluent User Guide. Site internet, 2020. URL [18] M. Benedict, G. B. Webb, and L. C. Rubin. An Empirical Equation for Thermodynamic Properties of Light Hydrocarbonsand Their Mixtures I. Methane, Ethane, Propane and 𝑛 -Butane. J. Chem. Phys. , 1940. (4), pp. 334–345. doi:10.1063/1.1750658[19] M. Faizan-Ur-Rab, S. Zahiri, S. H. Masood, M. Jahedi, and R. Nagarajah. Development of 3D Multicomponent Modelfor Cold Spray Process Using Nitrogen and Air. Coatings , 2015. , pp. 688–708. doi:10.3390/coatings5040688[20] A. Gerschenfeld. Mécanique des fluides pour l’énergie (Fluid mechanics for the energy field), 2020. (in French)[21] W. Sutherland. LII. The viscosity of gases and molecular force. Philos. Mag. , 1893. (223), pp. 507–531. doi:10.1080/14786449308620508. Series 5[22] V. Varadaraajan and P. Mohanty. Design and optimization of rectangular cold spray nozzle: Radial injection angle,expansion ratio and traverse speed. Surf. Coat. Technol. , 2017. , pp. 246–254. doi:10.1016/j.surfcoat.2017.03.005[23] R. B. Bird, W. E. Stewart, and E. N. Lightfoot.
Transport Phenomena . John Wiley and Sons, Chemical EngineeringDepartment University of Wisconsin-Madison, deuxième ed., 2002. ISBN 0-471-41077-2[24] Wikipedia. Hypothèse de Stokes (Stokes’ Hypothesis) — Wikipedia, the free encyclopedia, 2020. URL https://fr.wikipedia.org/wiki/Hypothese_de_Stokes . [Online; available on april-29-2020] (in French)[25] G. Billet, V. Giovangigli, and G. de Gassowski. Impact of Volume Viscosity on a Shock/Hydrogen Bubble Interac-tion.
Revue Interne de l’École Polytechnique , 2007. . URL [26] J. Shang, T. Wu, H. Wang, C. Yang, C. Ye, R. Hu, J. Tao, and X. He. Measurement of temperature-dependent bulkviscosities of nitrogen, oxygen and air from spontaneous Rayleigh-Brillouin scattering.
IEEE Access , 2019. , pp. 136439–136451. doi:10.1109/ACCESS.2019.2942219. URL [27] Z. Gu and W. Ubachs. Temperature-dependent bulk viscosity of nitrogen gas determined from spontaneous Rayleigh-Brillouin scattering. Opt. Lett. , 2013. (7), pp. 1110–1112[28] W.-Y. Li and C.-J. Li. Optimal Design of a Novel Cold Spray Gun Nozzle at a Limited Space. J. Therm. Spray Technol. ,2004. (3), pp. 391–396 2829] B. Samareh, O. Stier, V. Lüthen, and A. Dolatabadi. Assessment of CFD Modeling via Flow Visualization in Cold SprayProcess. J. Therm. Spray Technol. , 2009. (5-6), pp. 934–943. doi:10.1007/s11666-009-9363-9[30] M. Faizan-Ur-Rab, S. H. Zahiri, S.-H. Masood, M. Jahedi, and R. Nagarajah. 3D CFD Multicomponent Model for coldspray additive manufacturing of titanium particles. CFD Model. Simul. Mater. Process. 2016, Proc. Symp. pp. 213–220[31] O. C. Ozdemir and C. A. Widener. Influence of Powder Injection Parameters in High-Pressure Cold Spray.
J. Therm.Spray Technol. , 2017. , pp. 1411–1422. doi:10.1007/s11666-017-0606-x[32] A. Sova, I. Smurov, M. Doubenskaia, and P. Petrovskiy. Deposition of aluminum powder by cold spray micronozzle. TheInternational Journal of Advanced Manufacturing Technology , 2017. , pp. 3745–3752. doi:10.1007/s00170-017-1443-2[33] W. E. Ranz and W. R. Marshall. Evaporation from drops, Part I. Chem. Eng. Prog. , 1952. (3), pp. 141–146. URL http://dns2.asia.edu.tw/~ysho/YSHO-English/1000CE/PDF/CheEngPro48,141.pdf [34] S. Yin, X.-F. Wang, W.-Y. Li, and B.-P. Xu. Numerical Study on the Effect of Substrate Angle on Particle Impact Velocityand Normal Velocity Component in Cold Gas Dynamic Spraying Based on CFD. J. Therm. Spray Technol. , 2010. (6),pp. 1155–1162. doi:10.1007/s11666-010-9510-3[35] R. Lupoi and W. O’Neill. Powder stream characteristics in cold spray nozzles. Surf. Coat. Technol. , 2011. , pp.1069–1076. doi:10.1016/j.surfcoat.2011.07.061[36] B. Samareh and A. Dolatabadi. Dense Particulate Flow in a Cold Gas Dynamic Spray System.
J. Fluids Eng. , 2008. .doi:10.1115/1.2957914. URL http://citeseerx.ist.psu.edu/viewdoc/download?doi=10.1.1.822.8162&rep=rep1&type=pdf [37] L. Cui, A. G. Gerber, and G. C. Saha. Cold Gas Dynamic Spray Technology: The Simulation of Aerodynamics of Flow.
Key Eng. Mater. , 2019.
J. Therm. Spray Technol. , 2012. (3-4), pp. 541–549. doi:10.1007/s11666-011-9707-0[39] J. P. Du Plessis and J. H. Masliyah. Mathematical Modelling of Flow Through Consolidated Isotropic Porous Media. Transp. Porous Media , 1998. , pp. 145–161[40] C. T. Crowe. Drag Coefficient of Particles in a Rocket Nozzle. techreport 5, AIAA Journal, 1967. doi:10.2514/3.4119[41] R. Clift, J. R. Grace, and M. E. Weber. Bubbles, Drops, and Particles . Academic Press, Academic Press, Inc. (London) LTD.24/28 Oval Road London NW1, 1978. ISBN 0-12-176950-X. URL https://pdfs.semanticscholar.org/0f9f/6ab31a685e621053a6ebe89a9e00fc29dc96.pdf?_ga=2.145792564.1879786883.1588854142-707512567.1588854142 [42] S. A. Morsi and A. J. Alexander. An Investigation of Particle Trajectories in Two-Phase Flow Systems.
J. Fluid Mech. ,1972. (2), pp. 193–208. doi:10.1017/S0022112072001806[43] L. Jacquin. Mécanique des fluides (Fluid mechanics) . École Polytechnique, 2018. (in French)[44] M. L. Shur, P. R. Spalart, M. K. Strelets, and A. K. Travin. A hybird RANS-LES approach with delayed-DES andwall-modelled LES capabilities.
Int. J. Heat Fluid Flow , 2008. (29), pp. 1638–1649. doi:10.1016/j.ijheatfluidflow.2008.07.001. URL [45] M. S. Griskevich, A. V. Garbaruk, J. Schütze, and F. R. Menter. Development of DDES and IDDES Formula-tions for the 𝑘 − 𝜔 Shear Stress Transport Model.
Flow, Turbul. Combust. , 2011. (3), pp. 431–449. doi:10.1007/s10494-011-9378-4. URL https://cfd.spbstu.ru/agarbaruk/doc/2012_Gritskevich-et-al._Development-of-DDES-and-IDDES-Formulations-for-the-k-w-Shear-Stress-Transport-Model.pdf [46] P. D. Tam. Quelques rappels sur la notion de masse ajoutée en mécanique des fluides (Some reminders about the no-tion of added mass in fluid mechanics). Report, CEN Saclay BP n ◦ https://inis.iaea.org/collection/NCLCollectionStore/_Public/08/340/8340824.pdfhttps://inis.iaea.org/collection/NCLCollectionStore/_Public/08/340/8340824.pdf