INTRODUCTION
Turboshaft engines operating under environmentally degraded conditions (e.g., sand storms, soot clouds, etc.) experience severe performance losses and critical loss of operability due to the ingestion of foreign particulates into the propulsion system. The most common types of ingested particulates are volcanic ash, desert sand, dust, and salt, in the form of calcia–magnesia–alumino-silicates (CMAS) that often vary in composition depending on the geographical location. Although engine particle separator systems offer a potential solution to limit the amount of ingested contaminants, smaller particles (<80 μm) evade removal and are entrained onto the turbine stage. Extensive use of particle separation systems negatively impacts the overall system performance due to a reported high loss of inlet pressure and inlet mass flow, thus limiting their usage to aircraft regimes where the risk of contaminant ingestion is high [Reference Filippone and Bojdo1]. To date, contaminated environments remain a persistent challenge for the aviation community, not only in desert regions, but also to enable safer flights in complex megacities and urban canyon environments across the world (see Fig. 1 showing high-intensity regions of particulate matter less than 2.5 μm in diameter (PM2.5) extending across North Africa, South West Asia, and South East Asia [Reference Van Dokelarr, Martin, Brauer, Kahn, Levy, Verduzco and Villeneuve2]). Thus, the development of advanced thermal protection systems (TPS) and durable materials able to withstand all-weather environments remains a critical area of research for military and civilian aerospace sectors.
Figure 1: Global satellite-derived PM2.5 averaged over 2001–2006 showing local intensity of particulate matter less than 2.5 μm. White space indicates water or locations containing <50 measurements [Reference Van Dokelarr, Martin, Brauer, Kahn, Levy, Verduzco and Villeneuve2].
Due to the urgent need for higher power and efficiency, new technologies will continue to push the operating limits to higher temperatures to reap the benefits from Brayton cycle thermodynamics. Operating at such conditions introduces various thermophysical and thermochemical challenges that significantly affect the character of the CMAS particle-laden flow and its interaction with engineered protective materials such as thermal and environmental barrier coatings (T/EBCs). Once ingested, the particulates exposed to the harsh combusting environment undergo phase transformation to melt, with the possibility of CMAS depositing and reacting with the coatings on the combustor liner, turbine shroud, blades, and nozzle guide vanes [Reference Walock, Barnett, Ghoshal, Murugan, Swab, Pepi, Hopkins, Gazonas, Rowe, Kerner, J, D, M and A3,Reference Nieto, Walock, Ghoshal, Zhu, Gamble, Barnett, Murugan, Pepi, Rowe and Pegg4]. The ingested particulates also get into the bleed air intake from the compressor and settle mostly as unreacted but sintered sand in the cooling channels in the turbine vanes and blades. This leads to cooling channels blockage in the hot section, thereby reducing bleed air mass flow while further increasing bulk material temperatures [Reference Ghoshal, Murugan, Walock, Nieto, Barnett, Pepi, Swab, Zhu, Kerner, Rowe, Shiao, Hopkins and Gazonas5]. Thermal protection strategies used in today's turboshaft engines use T/EBCs along with bleed air film cooling of super alloy or ceramic matrix composite (CMC) components to protect the coatings from undergoing degradation [Reference Ghoshal, Murugan, Walock, Nieto, Barnett, Pepi, Swab, Zhu, Kerner, Rowe, Shiao, Hopkins and Gazonas5,Reference Murugan, Ghoshal, Walock, Barnett, Pepi and Kerner6]. Without protective T/EBC coatings, the moisture in the combusting atmosphere will rapidly degrade the underlying SiC–SiC CMC substrate material due to oxidation [Reference Larry Fehrenbacher, Kutsch, Vesnovsky, Fehrenbacher, Ghoshal, Walock, Murugan, Nieto, Mahmoud, Sridharan, Colorado, Bhalla, Singh, Gupta, Langhorn, Jitianu and Jose Manjooran7]. Although T/EBCs are effective at its current capability of 2400 °F, these strategies do not yet provide a robust solution to protect against increased CMAS attack at future engine hot section components exposed to 3000 °F (or higher temperatures) air flow, and as a result still undergo severe degradation.
To develop superior T/EBC systems that are resistant to CMAS attack at higher temperatures (>3000 °F), an integrated computational and experimental approach is essential to provide validated models and insights into the interaction physics. CMAS transport in turboshaft engines involve nonlinear physical and chemical processes that can be categorized into a sequence of interrelated physical stages comprising of (i) environmental entrainment and transport of inertial particles in turbulent flow regimes, followed by strongly coupled processes, including (ii) phase change, adherence to the surface and subsequent CMAS buildup, chemistry and phase evolutions at the surface contact, (iii) migration and infiltration into TBC porous structures, chemistry and phase change inside the T/EBC multilayers, and (iv) stress buildup and fracture mechanics of T/EBC spallation and hot tribocorrosion. These multi-physics processes occur over a wide range of spatio-temporal scales, thus requiring a well-defined canonical approach to provide the basis for understanding in more practical problems. From this description, and taking a modeling standpoint, it is evident that the turbulent transport of particulates sets the initial conditions for all subsequent stages and physical processes. Thus, the accurate modeling of hot particle transport and its deposition is critical for the overall modeling sequence of CMAS interactions with engineered materials. This work highlights recent progress with the current state-of-the-art computational approaches specifically for inertial particle transport in turbulent flows, deposition, and wettability using continuum mechanics. Although the complete spectrum of physical stages is inherently important, it is currently outside the scope of this paper and will be included in future reviews.
Annular multiphase flows, such as in turbine disks, can contain up to billions of particles that are transported and dispersed within the carrier phase of a turbulent flow field. One of the most important features of such flow is the presence of a wide spectrum of length scales and timescales associated with the smallest (particle) scale motions as well as a wide range of anisotropic turbulent flow structures. Depending on the governing flow parameters, such as particle-to-fluid length scales (l p/l k), timescale ratios (t p/t k), and particle Reynolds number (Rep), turbulence modulation effects become important due to various induced mechanisms, such as (a) enhanced dissipation due to the presence of particles, (b) the transfer of kinetic energy to the fluid from the particle, and (c) self-induced vortex shedding behind the particles, that can dramatically impact the dynamics [Reference Balachandar and Eaton8,Reference Guha9]. The level of interaction between the phases is typically determined by the fractional volume occupied by the dispersed phase (ϕv), as ϕv increases the level of interaction between particles tend to become stronger [Reference Elghobashi10]. The dispersion of particles and their evolution in unconfined homogeneous turbulent flows with one-way coupling are reasonably understood. On the other hand, flows in two-way and four-way coupling regimes still present a formidable challenge [Reference Balachandar and Eaton8–Reference Loth11], thus requiring detailed experimentation and numerical studies to improve their understanding. State-of-the-art methods are based either on resolving the smallest particle fluid motions, such as in direct numerical simulation (DNS), or using a macroscale approach which resolves larger turbulent flow structures, such as large-eddy simulation (LES), often coupled to a mesoscopic or subgrid scale (SGS) model for the dispersed phase [Reference Fox12,Reference Passalacqua, Laurent, Madadi-Kandjani, Heylmun and Fox13]. The ultimate option is to perform a fully resolved DNS, in which all scales of turbulence, and flow scales introduced by particles, are completely resolved. Although such simulations have been performed [Reference Uhlmann14], they are severely limited to a collection of up to ~103 particles. In Eulerian–Eulerian approaches, a two-fluid (TF) method is used to treat the carrier and dispersed phases as interpenetrating fluid media and the particle phase properties are given a field representation. The formulation requires additional momentum and energy equations for the particle phase, with momentum and energy exchange between the phases taken into account as source or sink terms [Reference Fox, Laurent and Massot15]. The Eulerian–Lagrangian (EL) approach is one of the most popular methods used for modeling turbulent dispersed flow computations [Reference Elghobashi16]. In this approach, the motion of the particles in fluids is based on Newton's equations for tracking along a trajectory with relevant forces acting on the particles. The main advantage of this approach is that there are no inherent limitations on particle size or distributions (e.g., mono or polydisperse mixtures) due to the direct treatment of particle motion, as demonstrated in Ref. [Reference Kozak, Dammati, Bravo, Hamlington and Poludnenko17]. Purely Lagrangian methods such as Smooth Particle Hydrodynamics (SPH) are formulated based on interparticle forcing and have been recently reported to capture two-phase flow dynamics with sufficient accuracy [Reference Ma, Zhang, Liu and Yao18,Reference Yang, Dai and Kong19]. Both Eulerian and Lagrangian approaches are often formulated under point particle approximation assumptions and are, therefore, to particles of a size smaller than the Kolmogorov scale in DNS, or the smallest resolved eddy in LES.
The evolution of the dispersed flow in a confined environment, such as those found in gas turbine engines, will naturally lead to particulate wall impingement and deposition. Deposition is governed by many physical processes, including thermophoresis, thermal diffusion, surface roughness, and electrostatic charges. These processes have been extensively reviewed in Ref. [Reference Guha9]. An understanding of particle deposition behavior in materials is vital to identifying vulnerable regions where infiltration and thermal stress are likely to occur. Another important factor is material wettability, often characterized by its contact angle (hydrophobic or hydrophilic) or moving contact line (MCL) behavior at the triple coexistence point. The challenge in modeling MCL arises from the so-called Huh–Scriven paradox which predicts a logarithmic singularity for the energy dissipation rate at the triple contact point. Multiple theories have been proposed to overcome this limitation by introducing boundary slip models, thin precursor films [Reference Snoeijer and Andreotti20]. Recent deposition studies using LES models have been developed based on the collision Weber number criteria for rebound or deposition [Reference Bravo, Xue, Murugan, Ghoshal, Walock and Flatau21]. These studies have shown the effects of particle size and Weber number on the deposition characteristics when considering a stator–rotor configuration at high temperatures. Other family of deposition models include stochastic eddy interaction models [Reference Graham and James22] and Markov chain-type turbulent random walk models [Reference Mofakham and Ahmadi23], in which the particle's trajectory is modeled as a succession of stochastic interactions with the turbulent eddies. The use of these models, however, has been limited to ensemble-averaged flows in dispersed solid particles with small fractional volumes largely due to the lack of physical model closures needed for strongly coupled flow. Recently, Bons et al. [Reference Bons, Prenter and Whitaker24] developed a “bounce stick” model. This formulation is based upon elasto-plastic deformation including the effects of adhesion and shear removal along with impact. The model is highly dependent on key material properties, such as particle yield stress σy, composite Young's modulus E c, surface free energy γ, and impulse ratio μimp, which remain difficult to measure experimentally (little data is available in the literature). In Ref. [Reference Singh and Tafti25], Singh et al. presented a novel deposition model that computes sticking probability based on both the particle local impact velocity and temperature. The model is implemented into an LES framework and the sticking efficiencies successfully demonstrated at temperatures as high as 1050 °C (1922 °F). Models capturing the sensitivity to physical and chemical composition properties (e.g., sand, dust, and ash) in complex turbulent transport still remain scarce. Simplified algebraic models implementing physio-chemical properties have been proposed [Reference Bojdo and Filippone26] based on the concept of thermal Stokes number to examine the relationship between time to equilibrate, residence time, and particulate size. This model has been used with some success to estimate deposition on practical nozzle vane engine components.
The objective of the present study is to investigate key physical aspects of CMAS particle transport and deposition at a hierarchy of scales by presenting novel physics-based models of wettability and deposition dynamics. Once demonstrated, the wettability models will inform the material design process to identify a TBC with sandphobic properties to be vetted at sand-laden gas turbine engine operating conditions. The work presented in this manuscript is the first step toward this goal. This manuscript is organized as follows: Section “Introduction” (current section) presents the motivation and a literature review on the physical processes and state-of-the-art models while highlighting its relevance to turboshaft rotorcraft engines; Section “Results and discussion” presents the technical results of material wettability (at the droplet scale) and deposition physics using LES models; Section “Conclusions and future work” summarizes key findings from our results and provides perspectives on future work; and Section “Computational methods” presents a detailed summary of the physics-based models and mathematical methods used in the study.
RESULTS AND DISCUSSION
To identify the local material vulnerability regions and degradation mechanisms due to CMAS attack, a physics-based model is needed that considers the phenomena at multiple scales. In this work, we have identified two important configurations that address this need by conducting studies of material wettability and turbulent multiphase deposition. In material wettability, a Lagrangian particle model is significantly modified to accommodate fluid-solid interparticle energy exchange terms required to model the behavior at the triple coexistence line. The configuration is a sessile droplet that is often used to study the material contact angle with ample experimental data available for comparison. In these simulations, the interparticle force models in the Lagrangian framework are validated with the droplet evolution properties (such as surface tension and contact angle). In deposition modeling, an LES approach is adopted to model the CMAS transport in turbulent flow environment and kinematic impact on turbine blades. It enables a systematic study of the impact of particle size, characterized by the Stokes number, on the deposition behavior. This includes identifying areas within the turbine section that are susceptible to particle deposition and the particle sizes that are primarily responsible for most of the damage. A demonstrated modeling capability is an important step toward the development of novel sandphobic coatings.
Material wettability analysis
This section presents numerical results obtained from a recently developed physics-based model for molten particle deposition to predict the behavior associated with contact angle dynamics. Contact angle is an important parameter to determine the wettability of material-fluid pairs. The SPH method is applied to conduct wettability analysis considering fluid-solid energy exchange mechanisms. Although other interface tracking models are available for two-phase flows, SPH provides various modeling advantages as described in the section “Computational methods”. This paper introduces interparticle force formulations to model the droplet evolution properties such as surface tension and wetting at the triple coexistence point. The computational method is presented in the section “Computational methods”.
The interaction of liquids over solid surfaces introduces an unphysical stress when using continuum hydrodynamic models at the triple point. This phenomenon is referred to as the Huh–Scriven paradox [Reference Huh and Scriven27]. Therefore, accurate contact dynamics requires a resultant nonzero force to model proper boundary slip at the triple point. Since the effect of gas particles at the triple point remains negligible for the present work, this force is obtained by determining the parameter σ s separately (see Section “Computational methods” for definition) for liquid–liquid interactions and for liquid–solid interactions. The model constants used are C ST for liquid–liquid interactions and C LS for liquid–solid interactions (see Section “Computational methods” for definition). The framework is evaluated by comparison with the sessile drop experiment conducted [Reference Du, Wang, Wen, Wang, Zhang and Yin28] at ambient conditions. These experiments were selected because the properties of the materials used are well-understood and reliable steady-state contact angle measurements are provided. In the experimental setup, water droplet of volume 4 × 10−6 L (droplet diameter D = 1.97 mm) is placed on a flat ultra-high-molecular-weight polyethylene (UHMWPE) substrate and allowed to settle at 25 °C before calculating the contact angle. The process is repeated five times and the average contact angle is measured to be 96.5 ± 2°. The surrounding medium is air at room temperature. The physical properties for water and polymer at room temperature are provided in Table 1.
TABLE 1: Properties of materials used in SPH sessile drop simulations.
CMAS, calcia–magnesia–alumina-silicate compound; SPH, smoothed particle hydrodynamics; UHMWPE, ultra-high-molecular-weight polyethylene polymer; YSZ, yttria-stabilized zirconia ceramic coating.
A series of droplet evolution studies are conducted to determine the interparticle physical forcing terms. These simulations require initialization of the droplet to an arbitrary shape and monitoring the shape evolution subject to various degrees of forcing to calculate the surface tension. To minimize energy, the initial droplet shape evolves to a spherical geometry at steady state due to the effects of surface tension. At equilibrium, the surface tension γ is calculated with the Young–Laplace equation and the pressure P inside the droplet. In this study, five tests are conducted using the forcing parameter values C ST, ranging from 0 to 100, indicating a linear relationship between C ST and γ as shown in Fig. 2(a). The value of C ST = 15.5 captures the correct surface tension behavior and evolution of the water droplet in these studies shown in the normalized pressure contours (Fig. 2(b)).
Figure 2: (a) Variation of the calculated surface tension with the model constant C ST and (b) visual representation of the shape evolution of liquid starting from arbitrary shape.
A domain comprising of dimensions 20D × 10D is selected for the wettability analysis. The simulations are initiated either with a circular or semi-circular sessile drop resting on the polymer surface. In both cases, the droplet is initialized with zero impact velocity corresponding to the experiment. A sweep of the liquid–solid strength parameter C LS is conducted to determine the effect on the droplet wettability properties. During the simulations, the droplet interface interacts with the polymer surface under gravity forces. Once the droplet reaches equilibrium, the curvature of the droplet near the triple coexistence point is extracted to determine the contact angle. Figure 3 shows the model wettability behavior determined as the contact angle subject to C LS forcing parameter variations. Note the behavior of the forcing term with contact angle is nonlinear in this case due to the interaction with the material properties. Such studies determine the forcing value of C LS = 9.3 to captures the correct wetting behavior with an accuracy of 2%. Figure 3 also shows the good agreement with experiment for the water–polymer pair.
Figure 3: (a) Sweep of the steady-state solution of sessile water droplet on a polymer surface modeled and compared with experiment by Du et al. [Reference Du, Wang, Wen, Wang, Zhang and Yin28] (square) and (b) SPH steady simulation and comparison with experiment with C ST = 15.5 and C LS = 9.3.
Furthermore, the wettability response and droplet shape at equilibrium are shown in Fig. 4 for select forcing terms. This shows a visual representation of each of the contact angles obtained through the parametric sweep. The results [Reference Mofakham and Ahmadi23,Reference Bons, Prenter and Whitaker24] show that the model is able to qualitatively capture the droplet hydrophobic to hydrophilic behavior over the reference material surface at ambient conditions.
Figure 4: Steady-state solution obtained from an SPH simulation of sessile water droplet on a polymer surface with C ST = 15.5 and C LS = 6, 10, and 14 (from left to right).
High-temperature sessile drop simulations are conducted to obtain model constants for the CMAS-T/EBC coating pair at 1200 °C. The experimental data available in the literature is very limited. Only recently [Reference Kang, Bai, Du, Yu, Bao, Wang and Ding29] researchers have investigated the wetting behavior of CMAS over Yttria-stabilized zirconia (YSZ) ceramic coatings through a sessile drop experiment. The study was conducted using argon gas as the inert medium and reported the temporal evolution of contact angle for the CMAS-YSZ. Figure 5 shows the physical sequence of selected images showing the wettability behavior in this study. The physical properties used of the CMAS-YSZ pair used in the current simulation was obtained from a literature survey and presented in Table 1. The modeled sessile drop is presented at steady state when the force equilibrium between particles is reached and the droplet no longer spreads. The most suitable comparison between simulated data and experimental data should be done at 4000 s (contact angle 60°) since the experimental droplet also reaches steady state at this time. Figure 5 illustrates the steady-state result of successful SPH simulation at high temperature for which the surface tension value and the contact angle predictions are the closest to the experiment. Further work is underway to determine the correct value of the CMAS surrogate mixture surface tension at high temperatures.
Figure 5: Steady-state solution of molten CMAS sessile droplet over YSZ coating at 1275 °C obtained from the experiment by Kang et al. [Reference Kang, Bai, Du, Yu, Bao, Wang and Ding29] (a) and from SPH simulation (b).
Turbulent multiphase deposition analysis with varying Stokes number
The analysis in this section considers the behavior of a dispersed cloud of CMAS particulates in a confined turbulent carrier fluid as well as during the initial entrainment into the boundary layer. Classical analysis of the transport of particles by the outer layer (inviscid flow) suggests that particles will be segregated by their convective Stokes number [Reference Guha9]. The Stokes number is defined as the ratio of characteristic time of particle to the characteristic time of the flow, St = (touo)/lo, where to is the relaxation time of the particle, uo is the freestream velocity, and lo is the characteristic length scale. Particles with small Stokes numbers (St = 0.2) will behave as tracer particles and will turn around the blade, while particles with large Stokes numbers (St = 20) will ballistically impinge on the blade surface.
This behavior is verified by considering particle-laden flow around a reference prototypical transonic guide vane described in Ref. [Reference Arts and De Rouvroit30]. The study is conducted by considering a von Karman Institute (VKI) turbine geometry with periodic boundary conditions to simulate a linear cascade array. Turbulence is modeled using the LES approach with coupled Lagrangian particle tracking solver and feedback accounting for drag and heat transfer. The inflow conditions are simulated using a cylindrical array to trip the flow state to turbulent transition. The nominal flow regime represents the first stator row downstream of the combustor section (Mach Number1 ≈ 0.2, Mach Number2 ≈ 0.9, and Rec, ≈ 3 × 105), where a molten sand particle in thermal equilibrium with combustion products at T∞ ≈ 2000 K and a blade temperature of Tw ≈ 1300 K is specified. Figure 6 shows the normalized mass deposition rate for three different Stokes numbers, ranging from 0.2 to 20, where the injected particles are convected axially toward the vane.
Figure 6: Normalized particle mass deposition rate on the vane surface. The gray scale here indicates surface particle count upon impact (black: 0 counts; white: 200 particle counts per unit surface area).
Deposition concentrations increase with increasing Stokes number and peak deposition occurs near the trailing edge of the pressure side, which is the natural impingement location based on the initial particle trajectories. This suggests that initial particle size distributions and trajectories play an important role in predicting the final deposition pattern. This has important implications in predicting related effects, specifically agglomeration or preferential concentration of particles due to large-scale anisotropic flow structures in the combustor.
Furthermore, the quantity of mass deposited across the blade vane surface is reported via the particle count per unit area along the blade vane surface (Fig. 7). Following the experimental convention, the particle count is reported in terms of suction side (0–86 mm) and pressure side (85–155 mm) along the blade surface. It is evident that the peak deposition occurs in the pressure side of the blade and that it decreases with Stokes number. At St = 20, peak deposition count at order of magnitude 109 is observed at 100 mm, St = 2 peak deposition decreases by an order of magnitude to 108 and shifts further upstream to 90 mm. At St = 0.2, the particles no longer affect have enough inertia to cross streamlines and entrain into the boundary layer region; hence, the particulate count is limited at 106. The results identify the vulnerable localized regions on the blade as the leading and trailing edge of the vane surface. These local deposit regions have also been observed in engine-scale experiments. Ballistic entrainment of the particulates in these regions is driven by the fluid sensitivity to the Stokes number and the combined effect of stagnation flow at the leading edge along with pressure effects from adjacent blades in the linear cascade. This analysis does not consider the effect of surface roughness and porosity of the T/EBC layer.
Figure 7: Particle count per unit area along the blade vane surface for Stokes number 0.2, 2, and 20. Convention specifies 0–86 mm as the suction side and 86–155 mm as the pressure side.
CONCLUSIONS AND FUTURE WORK
In this research work, we have conducted multiscale simulations of two important configurations useful for T/EBC coating technology development in gas-turbine engines. The configurations included (i) a sessile droplet to study material wettability phenomena at relevant conditions, followed by a study of (ii) multiphase turbulent flow and deposition on a turbine blade cascade. The wettability models presented provide the groundwork to determine the sandphobicity of candidate TBCs informing future material design. Gas turbine engines equipped with sandphobic TBCs will provide more protection against sand induced damage when operating in austere environmental conditions.
In material wettability, a Lagrangian SPH particle model was significantly extended to handle fluid–solid interparticle energy exchange terms to model the behavior at the triple coexistence line. This behavior was studied with a sessile droplet configuration and its results successfully compared with experiments. The study showed that determination of forcing constants are a critical step through fundamental studies of droplet shape evolution. The results agreed well with experiments and were able to capture a range of hydrophobic and hydrophilic behaviors by physically varying the intermolecular forcing terms. Turbulent flow particle deposition studies conducted using LES focused on VKI turbine array cascade. This studied the behavior of particle entrainment into the boundary layer. The results showed a high degree of sensitivity to the deposition characteristics with variation in particle Stokes number. Surface particle count results across the blade demonstrated that Stokes St = 20 particles lead to the highest deposits, as well as identified peak deposited regions in the leading and trailing edge regions.
In summary, much of the previous research has focused on predictions of multiphase flows of ideal Newtonian fluids and efficient algorithms are available for one- and two-way coupled unconfined flows. However, multiphase flows exhibiting non-Kolmogorov turbulence behavior in confined environments present many interesting challenges for computations due to their inherent multiscale, nonequilibrium, and nonlocal nature. These challenges are associated with describing the behavior and evolution of the droplet transport across multiple regimes, phase interface from the molecular scale (e.g., at the contact line), to the transport of a heterophase dispersed phase, to the description of larger anisotropic flow structures, and their mutual interactions. The development of a computational framework must be augmented concurrently with physio-chemical property uncertainty quantification which strongly impacts the adhesion and deposition characteristics [Reference Nieto, Walock, Ghoshal, Zhu, Gamble, Barnett, Murugan, Pepi, Rowe and Pegg4]. Direct simulations across the full spectrum of spatio-temporal scales, and encompassing all of the physics, in a unified framework is not affordable now and will remain inaccessible in the foreseeable future. Hence, the solution to describe such complex nonlinear systems is to utilize separate computational approaches where the essential information is exchanged to ultimately provide predictive insights toward the development of novel material systems.
COMPUTATIONAL METHODS
Smooth particle hydrodynamics
SPH is a mesh-free Lagrangian approach in which any field variable (scalar or vector) is represented by a set of discrete particles and evaluated by a finite kernel approximation to the smoothing function over neighboring particles. The advantages of using SPH over contemporary approaches (e.g., Volume of Fluids or Level Set methods) are in its ability to handle large free surface and deforming boundaries without mesh dependencies. This is important for wettability since it will capture the evolution of the triple-phase contact line based purely on particle force equilibrium without requiring any additional treatments to account for density difference between the phases. Since the domain does not require specific connectivity, SPH can easily simulate large deformation and track the physical evolution of a deformed drop and secondary behaviors upon any interactions [Reference Ma, Zhang, Liu and Yao18]. The schematic of the finite kernel and smoothing function radius for SPH is shown in Fig. 8.
Figure 8: Schematic view of an SPH kernel and smoothing function radius [Reference Ray, Yang, Kong, Bravo and Kweon31].
The continuum field property 
$A\lpar {{{r}}_a} \rpar$ of a particle a at position 
${ {r}}_a$ is obtained by performing convolution integral with a smoothing function 
$W\lpar {{ {r}}_a-{ {r}}_b\comma \;h} \rpar$ also known as the kernel,
$$A\lpar {{{r}}_a} \rpar = \sum\limits_{b = 1}^N {A\lpar {{{r}}_b} \rpar W\lpar {{{r}}_a-{{r}}_b\comma \;h} \rpar \displaystyle{{m_b} \over {\rho _b}}}, $$
where the summation is performed over N particles surrounding the particle i. 
$m_b\comma \;{ {r}}_b$, and ρ b are the mass, position, and density of the neighboring particle b, respectively, and h is the smoothing function. Therefore, the gradient approximation of function 
$A\lpar {{ {r}}_a} \rpar$ can be computed from Eq. (1) to obtain the governing equations in Lagrangian form.
The continuity equation and momentum governing equations are given by
$$\displaystyle{{{\rm d}\rho _a} \over {{\rm d}t}} = \sum\limits_b {m_b\lpar {{ {u}}_a-{ {u}}_b} \rpar \cdot \nabla _aW_{ab}}, $$
where ρ is the density, u is the velocity field, and 
$W_{ab} = W\lpar {{ {r}}_a-{ {r}}_b\comma \;h} \rpar$ is the smoothing function.
$$\displaystyle{{{\rm d}{ {u}}_a} \over {{\rm d}t}} = \sum\limits_b {m_b\left({\displaystyle{{P_a} \over {{\rm \rho}_a^2 }} + \displaystyle{{P_b} \over {{\rm \rho}_b^2 }}} \right)\nabla _aW_{ab} + g_a}, $$where P is the pressure and g a is the body force acting on the particle of interest. A complete description of the SPH framework is given in Refs. [Reference Ray, Yang, Kong, Bravo and Kweon31,Reference Yang and Kong32]. In the present work, an interparticle forcing term formulation is introduced to model the droplet evolution properties such as surface tension and wetting at the triple coexistence point. The interparticle force introduced is given by
$$F_{ab}^I = \displaystyle{{{\rm \sigma} _sh} \over {s_2}}\displaystyle{{r_{ab}} \over {r_{ab}}} \times \left. {\left({\matrix{ {{\lpar {s_{ab}-s_1} \rpar }^2-s_2\comma} \hfill & {0 \le s_{ab} \lt 1 + s_1/2} \hfill \cr {{\lpar {s_{ab}-s_1} \rpar }^2\comma} \hfill & {1 + s_1/2 \le s_{ab} \lt 2} \hfill \cr {0\comma} \hfill & {{\rm else}} \hfill \cr } } \right.} \right\},$$where s ab = 2r ab/kh, kh is 1.5 times the smoothing distance h, s 2 = (2 − s 1)2/2. The parameter σs here defines the strength of the interparticle force and is determined independently for the liquid–liquid interactions CST and for the liquid–solid interactions C LS [Reference Yang, Dai and Kong19]. Verification of the forcing parameters are presented in detail in the "Results" section.
LES methods
The multiphase turbulent flow simulations were conducted using the compressible, unstructured-mesh LES solver CharLES developed by Cascade Technologies Inc. CharLES is a fully explicit solver, using a third-order Runge–Kutta formulation in time and a low-dissipative finite-volume scheme in space [Reference Brès, Bose, Emory, Ham, Schmidt, Rigas and Colonius33–Reference Jain, Bravo, Kim, Murugan, Ghoshal, Ham and Flatau36]. A blending method that combines nondissipative central flux with a dissipative upwind flux is used, providing computational stability in coarse grain regions. The numerical dissipation due to the use of the upwind scheme is minimized since the upwind flux is applied only in limited regions, depending on the local quality of the mesh. The solver is second-order accurate in space. The Vreman SGS model [Reference Vreman A37] is used to represent SGS structures. For the multiphase particulate modeling, a hybrid EL approach is used with a two-way coupling that accounts for drag forces and heat transfer exchange as follows. The Lagrangian liquid particle motion is simulated using Basset–Boussinesq–Oseen (BBO) equations [Reference Crowe, Sommerfeld and Tsuju38]. It is assumed that the density of the particle is much larger than that of the gaseous fluid (O[1000]), that the particle size is small compared with the turbulence integral length scale, and that the effect of shear on the particle motion is negligible. The large density ratio means that the Basset force and the added mass term are small and, therefore, are neglected. Under these assumptions, the particle motion is governed by the Lagrangian equations below:
$$\displaystyle{{{\rm d}x_p} \over {{\rm d}t}} = u_p,$$
$$\displaystyle{{{\rm d}u_p} \over {{\rm d}t}} = D_p\lpar u-u_p\rpar, $$where u p and u are the droplet velocities and gas phase velocities at the droplet location. The drag force on the particle is modeled by a drag coefficient Cd [Reference Mofakham and Ahmadi23] for a spherical particle as
$$D_p = \displaystyle{3 \over 4}C_d\displaystyle{{{\rm \rho} _g} \over {{\rm \rho} _l}}\displaystyle{{\vert {{ u}-{ u}_p} \vert } \over {d_p}},$$where dp is the particle diameter. Deposition modeling is achieved using a simple stick model that is based on the particle impact speed and the local Weber number. Further details of the Lagrangian particle tracking methodology in Ref. [Reference Bravo, Kim, Ham and Kerner39].
A Voronoi approach is used to generate the mesh topology. This methodology generates a unique grid by taking the geometry surface description and a set of seed points as input. For the given geometry, the surface topology was generated first followed by seeding of specified points into the computational domain. The Voronoi concept was then employed to generate the complete grid by creating polygons in a manner that each polygon comprises exactly one seed and every point in that polygon is closer to its seed that any other seeds. The grid was then smoothed using Lloyd iterations [Reference Jain, Bravo, Kim, Murugan, Ghoshal, Ham and Flatau36,Reference Du, Faber and Gunzburger40] leading to anisotropic nearly hexagonal body-fitted meshes. During the smoothing operation, the generating points were moved to the centroids of their Voronoi cells, and the Voronoi mesh was recomputed (see Fig. 9).
Figure 9: Mesh smoothing using Lloyd iterations where seed points (red) move toward cell centroids (black +) with each iteration and the Voronoi diagram is recomputed [Reference Jain, Bravo, Kim, Murugan, Ghoshal, Ham and Flatau36].
An important benefit of including Lloyd iterations is that this smoothing step can be used to produce boundary alignment of the generating points near surfaces. It can also be restricted to a subset of the point cloud so that it does not disturb pre-existing point spacing.
The Voronoi-based grid generation technique is readily employed to generate mesh for complex geometries by specifying only a set of seed points. The refinement of the mesh in the desired region is also easily achievable by changing the location of seed points. The resulting mesh always maintains the orthogonality of face normal and cell displacement vectors, thus enabling high computational efficiencies for a dynamic mesh generator and a fluid solver. Coupled with a proper flux discretization method, Voronoi mesh points can be moved at arbitrary speeds independent of the fluid velocity. As such, one can have a boundary-fitted mesh component for the complex geometrical surface that interacts dynamically with the mesh component in the region of the moving fluid. The connectivity between the two components can be quickly reconstructed at each time step removing relatively slow heuristics/interpolations normally associated with unstructured grids.
Acknowledgments
This work was supported in part by resources from the DoD High-Performance Computing Modernization Program (HPCMP) FRONTIER Award to ARL with project title “Petascale High Fidelity Simulation of Atomization and Spray/Wall Interactions” . L.B., A.G., and M.M. were supported by the VTD 6.1 basic research mission program in propulsion sciences and a DoD Laboratory University Collaborative Initiative (LUCI) Fellowship. The simulations were run on the Centennial HPC System at the ARL DSRC. This work was conceptualized through ARL participation in the 2018 Center for Turbulence Research Summer Program at Stanford University.
