2.1 A mesoscale description for compressible two-phase flows

To arrive at a mesoscale description for compressible particle-laden flows, the Navier–Stokes equations are split into microscale processes that take place on the scale of a particle and below, and meso- to macroscale processes that take place on a scale much larger than the particle size. Anderson & Jackson (Reference Anderson and Jackson1967) provide such a basis for this approach by applying a local volume filtering operator to the incompressible Navier–Stokes equations, thereby replacing the point variables (fluid velocity, pressure, etc.) by filtered fields. Applying a similar approach to the viscous compressible Navier–Stokes equations, the volume-filtered conservation equations can be expressed as (see appendix A for details)

(2.1) $$\begin{eqnarray}\displaystyle & \displaystyle \frac{\unicode[STIX]{x2202}\unicode[STIX]{x1D6FC}\unicode[STIX]{x1D70C}}{\unicode[STIX]{x2202}t}+\unicode[STIX]{x1D735}\boldsymbol{\cdot }(\unicode[STIX]{x1D6FC}\unicode[STIX]{x1D70C}\boldsymbol{u})=0, & \displaystyle\end{eqnarray}$$

(2.2) $$\begin{eqnarray}\displaystyle & \displaystyle \frac{\unicode[STIX]{x2202}\unicode[STIX]{x1D6FC}\unicode[STIX]{x1D70C}\boldsymbol{u}}{\unicode[STIX]{x2202}t}+\unicode[STIX]{x1D735}\boldsymbol{\cdot }(\unicode[STIX]{x1D6FC}\{\unicode[STIX]{x1D70C}\boldsymbol{u}\otimes \boldsymbol{u}+p\mathbb{I}-\unicode[STIX]{x1D749}\})=(p\mathbb{I}-\unicode[STIX]{x1D749})\boldsymbol{\cdot }\unicode[STIX]{x1D735}\unicode[STIX]{x1D6FC}+\boldsymbol{F}, & \displaystyle\end{eqnarray}$$

and

(2.3) $$\begin{eqnarray}\displaystyle \frac{\unicode[STIX]{x2202}\unicode[STIX]{x1D6FC}\unicode[STIX]{x1D70C}E}{\unicode[STIX]{x2202}t}+\unicode[STIX]{x1D735}\boldsymbol{\cdot }\unicode[STIX]{x1D6FC}(\{\unicode[STIX]{x1D70C}E+p\}\boldsymbol{u}-\boldsymbol{u}\boldsymbol{\cdot }\unicode[STIX]{x1D749})+\unicode[STIX]{x1D6FC}\unicode[STIX]{x1D735}\boldsymbol{\cdot }\boldsymbol{q}=(\unicode[STIX]{x1D749}-p\mathbb{I})\boldsymbol{ : }\unicode[STIX]{x1D735}(\unicode[STIX]{x1D6FC}_{p}\boldsymbol{u}_{p})+\boldsymbol{u}_{p}\boldsymbol{\cdot }\boldsymbol{F}, & & \displaystyle\end{eqnarray}$$

where $\unicode[STIX]{x1D6FC}$ is the fluid-phase volume fraction, $\unicode[STIX]{x1D70C}$ is the fluid density, $\boldsymbol{u}$ the fluid velocity, $\boldsymbol{u}_{p}$ is the particle-phase velocity (in an Eulerian frame of reference) and $E$ the total energy. Interphase heat transfer (based on the flow configuration in § 3) had negligible effect on the results and is neglected here. Flow variables in equations (2.1)–(2.3) have been non-dimensionalized by ambient density $\unicode[STIX]{x1D70C}_{\infty }^{\star }$ , speed of sound $c_{\infty }^{\star }$ , a characteristic length scale $L^{\star }$ (based on vorticity thickness defined in equation (3.3)) and heat capacity at constant pressure $C_{p}^{\star }$ . Dimensional quantities are denoted by a superscript $\star$ , and the subscript $\infty$ indicates reference quantities (taken to be air). The source term $\boldsymbol{F}$ appearing in equations (2.2) and (2.3) accounts for momentum coupling between the particle and gas phases, and its form will be discussed in § 3.2. The non-dimensional viscous stress tensor is given by

(2.4) $$\begin{eqnarray}\displaystyle \unicode[STIX]{x1D749}=\frac{\unicode[STIX]{x1D707}}{Re_{c}}(\unicode[STIX]{x1D735}\boldsymbol{u}+\unicode[STIX]{x1D735}\boldsymbol{u}^{\mathsf{T}})+\frac{\unicode[STIX]{x1D706}}{Re_{c}}\unicode[STIX]{x1D735}\boldsymbol{\cdot }\boldsymbol{u} & & \displaystyle\end{eqnarray}$$

and the heat flux $\boldsymbol{q}$ is

(2.5) $$\begin{eqnarray}\displaystyle \boldsymbol{q}=-\frac{\unicode[STIX]{x1D707}}{Re_{c}Pr}\unicode[STIX]{x1D735}T, & & \displaystyle\end{eqnarray}$$

where $Pr\equiv C_{p}^{\star }\unicode[STIX]{x1D707}^{\star }/k^{\star }=0.7$ is the Prandtl number, with $\unicode[STIX]{x1D707}^{\star }$ and $k^{\star }$ the dynamic viscosity and thermal conductivity, respectively. The Reynolds number used in the formulation is defined as $Re_{c}=Re/M$ , where $Re=\unicode[STIX]{x1D70C}_{\infty }^{\star }\unicode[STIX]{x0394}U^{\star }L^{\star }/\unicode[STIX]{x1D707}_{\infty }^{\star }$ is the flow Reynolds number with $\unicode[STIX]{x0394}U^{\star }$ a characteristic velocity defined later, and $M=\unicode[STIX]{x0394}U^{\star }/c_{\infty }^{\star }$ is the Mach number. The non-dimensional viscosity is modelled as a power law $\unicode[STIX]{x1D707}=[(\unicode[STIX]{x1D6FE}-1)T]^{n}$ , with $n=0.666$ as a model for air and $\unicode[STIX]{x1D6FE}=1.4$ is the ratio of specific heats. The second coefficient of viscosity is given by $\unicode[STIX]{x1D706}=\unicode[STIX]{x1D707}_{B}-\frac{2}{3}\unicode[STIX]{x1D707}$ , where the bulk viscosity $\unicode[STIX]{x1D707}_{B}=0.6\unicode[STIX]{x1D707}$ is chosen as a model for the bulk viscosity of air. Assuming an ideal gas, the thermodynamic pressure and temperature depend on

(2.6) $$\begin{eqnarray}\displaystyle p=(\unicode[STIX]{x1D6FE}-1)(\unicode[STIX]{x1D70C}E-{\textstyle \frac{1}{2}}\unicode[STIX]{x1D70C}\boldsymbol{u}\boldsymbol{\cdot }\boldsymbol{u}) & & \displaystyle\end{eqnarray}$$

and

(2.7) $$\begin{eqnarray}\displaystyle T=\frac{\unicode[STIX]{x1D6FE}p}{\unicode[STIX]{x1D70C}(\unicode[STIX]{x1D6FE}-1)}. & & \displaystyle\end{eqnarray}$$

Since pressure is a key observable in the simulations, its transport equation is useful for examining the competition of the mechanisms generating it. Differentiating equation (2.6) with respect to time and substituting equations (2.1)–(2.3) yields the pressure evolution equation in the presence of a disperse phase,

(2.8) $$\begin{eqnarray}\displaystyle \frac{\text{D}p}{\text{D}t}=-\unicode[STIX]{x1D6FE}p\unicode[STIX]{x1D735}\boldsymbol{\cdot }\boldsymbol{u}+{\mathcal{D}}-\unicode[STIX]{x1D6FE}p\frac{\text{D}\ln \unicode[STIX]{x1D6FC}}{\text{D}t}+\frac{\unicode[STIX]{x1D6FE}-1}{\unicode[STIX]{x1D6FC}}(\boldsymbol{u}_{p}-\boldsymbol{u})\boldsymbol{\cdot }\boldsymbol{F}, & & \displaystyle\end{eqnarray}$$

where $\text{D}/\text{D}t=\unicode[STIX]{x2202}/\unicode[STIX]{x2202}t+\boldsymbol{u}\boldsymbol{\cdot }\unicode[STIX]{x1D735}$ is the material derivative operator and terms involving molecular transport effects are combined into

(2.9) $$\begin{eqnarray}\displaystyle {\mathcal{D}}=\frac{\unicode[STIX]{x1D6FE}-1}{\unicode[STIX]{x1D6FC}}(\unicode[STIX]{x1D749}\boldsymbol{ : }\unicode[STIX]{x1D735}\unicode[STIX]{x1D6FC}\boldsymbol{u}+\unicode[STIX]{x1D749}\boldsymbol{ : }\unicode[STIX]{x1D735}\unicode[STIX]{x1D6FC}_{p}\boldsymbol{u}_{p})+\unicode[STIX]{x1D6FE}\unicode[STIX]{x1D735}\boldsymbol{\cdot }\frac{\unicode[STIX]{x1D707}}{Re_{c}\,Pr}\unicode[STIX]{x1D735}(p/\unicode[STIX]{x1D70C}). & & \displaystyle\end{eqnarray}$$

Compared to the single-phase pressure transport equation (e.g. Pantano & Sarkar Reference Pantano and Sarkar2002), the last two terms on the right-hand side of equation (2.8) represent new contributions due to the presence of particles. The second-to-last term on the right-hand side involves convection aligned with volume fraction gradients as well as a mechanism associated with the so-called $pDV$ work term (Lhuillier, Theofanous & Liou Reference Lhuillier, Theofanous and Liou2010). The last term on the right-hand side represents work due to drag which is only active when the slip velocity magnitude between the phases is $|\boldsymbol{u}_{p}-\boldsymbol{u}|>0$ .

2.2 Averaged pressure intensity

For examining mechanisms affecting the pressure intensity in turbulence, we use Reynolds decomposition with spatial averaging operator $\overline{(\cdot )}$ , such that any quantity $A$ can be decomposed into $A=\overline{A}+A^{\prime }$ as is often used to examine Reynolds stress or turbulent kinetic energy transport.

Using this decomposition, multiplying equation (2.8) by $p$ and rearranging yields the following pressure intensity transport equation

(2.10) $$\begin{eqnarray}\displaystyle & & \displaystyle \frac{\unicode[STIX]{x2202}\overline{p^{\prime }p^{\prime }}}{\unicode[STIX]{x2202}t}+\overline{\boldsymbol{u}}\boldsymbol{\cdot }\unicode[STIX]{x1D735}\overline{(p^{\prime }p^{\prime })}=-2\unicode[STIX]{x1D6FE}\overline{p}~\overline{p^{\prime }\unicode[STIX]{x1D735}\boldsymbol{\cdot }\boldsymbol{u}^{\prime }}-(2\unicode[STIX]{x1D6FE}-1)\overline{p^{\prime }p^{\prime }\unicode[STIX]{x1D735}\boldsymbol{\cdot }\boldsymbol{u}^{\prime }}\nonumber\\ \displaystyle & & \displaystyle \quad -\,2\overline{\boldsymbol{u}^{\prime }p^{\prime }}\boldsymbol{\cdot }\unicode[STIX]{x1D735}\overline{p}-\unicode[STIX]{x1D735}\boldsymbol{\cdot }(\overline{p^{\prime }p^{\prime }\boldsymbol{u}^{\prime }})-2\unicode[STIX]{x1D6FE}\overline{p^{\prime }p^{\prime }}\unicode[STIX]{x1D735}\boldsymbol{\cdot }\overline{\boldsymbol{u}}+\underbrace{2\overline{p^{\prime }{\mathcal{D}}^{\prime }}}_{\hspace{-72.0pt}\text{viscous and heat conduction}}\nonumber\\ \displaystyle & & \displaystyle \quad \times \,\underbrace{-2\unicode[STIX]{x1D6FE}\left(\overline{p^{\prime }p^{\prime }\frac{\text{D}\ln (\unicode[STIX]{x1D6FC})}{\text{D}t}}+\overline{p}~\overline{p^{\prime }\frac{\text{D}\ln (\unicode[STIX]{x1D6FC})}{\text{D}t}}\right)}_{\text{volume}\;\text{displacement}}\nonumber\\ \displaystyle & & \displaystyle +\,\underbrace{2(\unicode[STIX]{x1D6FE}-1)\left(\overline{\frac{p^{\prime }\boldsymbol{F}}{\unicode[STIX]{x1D6FC}}}\boldsymbol{\cdot }(\overline{\boldsymbol{u}_{p}}-\overline{\boldsymbol{u}})+\overline{\frac{p^{\prime }\boldsymbol{F}}{\unicode[STIX]{x1D6FC}}\boldsymbol{\cdot }(\boldsymbol{u}_{p}^{\prime }-\boldsymbol{u}^{\prime })}\right)}_{\text{drag}\;\text{coupling}}.\end{eqnarray}$$

The first two lines of equation (2.10) have the same form as for single-phase flows, and they have been analysed for several turbulence configurations (e.g. Sarkar Reference Sarkar1992). As shown by Sarkar (Reference Sarkar1992), $\overline{p^{\prime }\unicode[STIX]{x1D735}\boldsymbol{\cdot }\boldsymbol{u}^{\prime }}$ acts to transfer energy between the kinetic and internal energy of turbulence. The term $\unicode[STIX]{x1D735}\boldsymbol{\cdot }\overline{(p^{\prime }p^{\prime }\boldsymbol{u})}$ has a turbulent-flux-like form for pressure intensity, $-2\unicode[STIX]{x1D6FE}\overline{p^{\prime }p^{\prime }}\unicode[STIX]{x1D735}\boldsymbol{\cdot }\overline{\boldsymbol{u}}$ has a production-like form for pressure intensity, although depending on mean flow expansion rather than mean shear, and $2\overline{p^{\prime }{\mathcal{D}}^{\prime }}$ combines all of the heat conduction and viscous mechanisms based on (2.9). The terms involving $\text{D}/\text{D}t(\ln \unicode[STIX]{x1D6FC})$ show how pressure intensity is modified from the material transport of volume fraction (depending on its correlation to the local pressure) and the last term describes the effects due to drag. The analysis of equation (2.10) in subsequent sections focuses on groups of terms: (i) volume displacement; (ii) drag coupling; (iii) viscous and heat conduction; and (iv) remaining terms without $\overline{(p^{\prime }p^{\prime })}_{t}$ that comprise what we call the inviscid hydrodynamics. The following sections introduce the simulations used in evaluating the terms in (2.10); their relative importance on modifying pressure fluctuations will be presented in § 4.2.

3.1 Flow configuration

The simulations are designed to provide a model flow to represent phenomenology of full-scale droplet-injected jet applications. To this end, we focus on what can be considered the near-nozzle exit region of a high-speed shear layer. Assuming the turbulent mixing layers from these jet exits are thin relative to their curvature, the turbulence would develop similarly as a planar spatially developing shear layer. Using a coordinate transformation based on the mean streamwise flow, we consider the temporally developing frame of reference, as shown in figure 1. Focusing on the thin shear layer turbulence enables the direct simulation of a broader range of turbulence scales than could be represented in a full jet simulation; the largest scales there being the size of the jet diameter rather than the shear layer thickness. Thus, the present configuration provides a Reynolds-number realistic representation of a section of a high-Reynolds-number jet, and it enables probing of the sound generation mechanisms of high-Reynolds-number turbulence. We mention that the acoustic far field from the current mixing layers will not include any geometric propagation associated with that from a round jet nor from the closing of shear layers at the end of the potential core. For the very near-field analysis considered here (i.e. within $20\,\unicode[STIX]{x1D6FF}_{m}$ where the momentum thickness $\unicode[STIX]{x1D6FF}_{m}$ is given in equation (3.2)), results are expected to compare well with those near jets. One such comparison was made between temporally developing shear layer direct numerical simulation and round jet large-eddy simulation (LES) which showed similar pressure skewness $S_{k}$ with Mach number (Buchta & Freund Reference Buchta and Freund2017). Although such a configuration precludes a one-to-one comparison to far-field jet sound, it provides a detailed description of underlying particle–turbulence and acoustic interactions and motivates analysis for more complex configurations.

Figure 1.

Shear layer configuration: (a) iso-view of the computational domain, (b) zoomed-in view of the $x$ – $y$ plane and (c) detail of (b). The vorticity magnitude is coloured by $(\unicode[STIX]{x1D735}\times \boldsymbol{u})<0.5\unicode[STIX]{x0394}U/\unicode[STIX]{x1D6FF}_{m}$ increasing from blue to red with particle position in black. Grey scales of dilatation are from $|\unicode[STIX]{x1D735}\boldsymbol{\cdot }\boldsymbol{u}|<0.01\unicode[STIX]{x0394}U/\unicode[STIX]{x1D6FF}_{m}$ decreasing from black to white.

In practice, water droplets are introduced into the jet exit shear layers by water-jet injection ports positioned azimuthally near the nozzle exit (Krothapalli et al. Reference Krothapalli, Venkatakrishnan, Lourenco, Greska and Elavarasan2003; Greska Reference Greska2005), which would be challenging to consider in the planar temporally developing configuration here. However, experimental measurements (Krothapalli et al. Reference Krothapalli, Venkatakrishnan, Lourenco, Greska and Elavarasan2003) indicate that the injected liquid breaks up quickly in the high-speed shear layer and forms a cloud of small droplets (mainly $d_{p}<5~\unicode[STIX]{x03BC}\text{m}$ ) at the end of the injection region where turbulence attenuation is observed. Therefore, by seeding particles directly inside developed turbulence, focus is placed on turbulence–particle coupling mechanisms, multiphase development in a compressible regime and near-field sound changes under such conditions. Future work may be warranted to assess how the current turbulence–particle–acoustic coupling generalizes to other injection–turbulence interactions and for evaporation in heated jets.

The gas phase is initialized with velocity

(3.1) $$\begin{eqnarray}\displaystyle [u,v,w](\boldsymbol{x})=\left[\frac{\unicode[STIX]{x0394}U}{2}\tanh \left(\frac{y}{2\unicode[STIX]{x1D6FF}_{m}^{o}}\right)+u^{\prime }(\boldsymbol{x}),v^{\prime }(\boldsymbol{x}),w^{\prime }(\boldsymbol{x})\right], & & \displaystyle\end{eqnarray}$$

with $\unicode[STIX]{x0394}U$ the velocity difference between the upper and lower streams and velocity perturbations ( $[u^{\prime }(\boldsymbol{x}),v^{\prime }(\boldsymbol{x}),w^{\prime }(\boldsymbol{x})]$ ) are comprised of a broadband Fourier based form localized in the shear layer. Details of the additive perturbations used for triggering turbulence are explained in full elsewhere (Capecelatro & Buchta Reference Capecelatro and Buchta2017). The initial pressure is assumed constant and density is computed from the Crocco–Busemann relationship, which is an established initialization of similar shear layer turbulence simulations (Vishnampet, Bodony & Freund Reference Vishnampet, Bodony and Freund2015; Buchta & Freund Reference Buchta and Freund2017).

For each simulation presented herein, the domain lengths in the $x$ (streamwise), $y$ (cross-stream) and $z$ (spanwise) directions are $L_{x}=768\unicode[STIX]{x1D6FF}_{m}^{o}$ , $L_{y}=1600\unicode[STIX]{x1D6FF}_{m}^{o}$ and $L_{z}=128\unicode[STIX]{x1D6FF}_{m}^{o}$ with $768\times 1601\times 192$ grid points, respectively, based on the initial momentum thickness ( $\unicode[STIX]{x1D6FF}_{m}^{o}$ ) using

(3.2) $$\begin{eqnarray}\displaystyle \unicode[STIX]{x1D6FF}_{m}(t)=\frac{1}{\unicode[STIX]{x1D70C}_{\infty }\unicode[STIX]{x0394}U^{2}}\int _{-Ly/2}^{L_{y}/2}\overline{\unicode[STIX]{x1D70C}}\left[\frac{1}{2}\unicode[STIX]{x0394}U-\widetilde{u}(y,t)\right]\left[\frac{1}{2}\unicode[STIX]{x0394}U+\widetilde{u}(y,t)\right]\,\text{d}y, & & \displaystyle\end{eqnarray}$$

evaluated at $t=0$ . Tests indicated that near-field pressure fluctuations of the single-phase (unladen) case were independent of the domain size; moreover, they are similar to those parameters used in previous single-phase shear layer direct numerical simulations (DNS) (Buchta Reference Buchta2016; Buchta & Freund Reference Buchta and Freund2017). The Reynolds number based on the velocity difference and the initial momentum thickness is set to $Re_{\unicode[STIX]{x1D6FF}_{m}^{o}}=120$ . Turbulence spectra discussed in § 4.1 support a sufficient grid for resolving turbulence based on this initial Reynolds number. To model water droplets in air, the particle density is taken to be $\unicode[STIX]{x1D70C}_{p}=1000\,\unicode[STIX]{x1D70C}_{\infty }$ . The particle diameter, $d_{p}$ , is then chosen to yield a desired initial Stokes number $St=\unicode[STIX]{x1D70F}_{p}/\unicode[STIX]{x1D70F}_{f}$ based on the particle and fluid time scales, $\unicode[STIX]{x1D70F}_{p}=\unicode[STIX]{x1D70C}_{p}d_{p}^{2}/(18\unicode[STIX]{x1D707})$ and $\unicode[STIX]{x1D70F}_{f}=\unicode[STIX]{x1D6FF}_{\unicode[STIX]{x1D714}}/\unicode[STIX]{x0394}U$ , respectively. The fluid time scale $\unicode[STIX]{x1D70F}_{f}$ is based on the vorticity thickness of the shear layer

(3.3) $$\begin{eqnarray}\displaystyle \unicode[STIX]{x1D6FF}_{\unicode[STIX]{x1D714}}(t)=\frac{\unicode[STIX]{x0394}U}{|\text{d}\tilde{u} (y,t)/\text{d}y|_{max}}, & & \displaystyle\end{eqnarray}$$

with $\tilde{u} (y,t)$ the Favre-average streamwise velocity. For the majority of the simulations, the initial Stokes number is set to $St=1.0$ . Sensitivity to Stokes number for $0.25\leqslant St\leqslant 4$ is also examined for $M=1.5$ . The Stokes number based on the Kolmogorov time scale ( $\unicode[STIX]{x1D70F}_{\unicode[STIX]{x1D702}}$ ), $St_{\unicode[STIX]{x1D702}}=\unicode[STIX]{x1D70F}_{p}/\unicode[STIX]{x1D70F}_{\unicode[STIX]{x1D702}}$ , is in the range $3.4\leqslant St_{\unicode[STIX]{x1D702}}\leqslant 4.4$ (for $St=1$ ) at the time of particle seeding, where $\unicode[STIX]{x1D70F}_{\unicode[STIX]{x1D702}}$ depends on the local dissipation rate in the shear layer. For reference, the Stokes number for the water-droplet injection in high-speed jets, based on $10\,\unicode[STIX]{x1D6FF}_{w}^{o}$ and $d_{p}=4~\unicode[STIX]{x03BC}\text{m}$ , was estimated to be 0.2 (Krothapalli et al. Reference Krothapalli, Venkatakrishnan, Lourenco, Greska and Elavarasan2003). Initially, all of the simulations developed from vorticity thickness $\unicode[STIX]{x1D6FF}_{w}^{o}\equiv \unicode[STIX]{x1D6FF}_{w}(t=0)=1$ , based on equations (3.1) and (3.3), to thickness $\unicode[STIX]{x1D6FF}_{w}(t=t^{\bullet })\approx 10$ , which was sufficient time to establish turbulence; at this time, $t=t^{\bullet }$ , the particles were positioned inside the shear layer as previously described in Capecelatro & Buchta (Reference Capecelatro and Buchta2017). Rather than add particles at $t=0$ , which can affect turbulence transition, seeding them in turbulence is similar to micro-jet injection application for high-Reynolds-number full-scale jets. The total number of particles considered in each simulation is based on the volume fraction $\unicode[STIX]{x1D6F7}_{v}$ within the shear layer, defined as

(3.4) $$\begin{eqnarray}\displaystyle \unicode[STIX]{x1D6F7}_{v}=\frac{N_{p}\unicode[STIX]{x03C0}d_{p}^{3}}{6L_{x}L_{z}\unicode[STIX]{x1D6FF}_{\unicode[STIX]{x1D714}}^{o}}. & & \displaystyle\end{eqnarray}$$

A summary of relevant parameters used in each case is provided in table 1. Across all of the laden simulations considered, the total number of particles varies in the range $3\times 10^{6}\lesssim N_{p}\lesssim 250\times 10^{6}$ depending on mass loading and Stokes number. The corresponding mass loading within the shear layer

(3.5) $$\begin{eqnarray}\displaystyle \unicode[STIX]{x1D6F7}_{m}=\frac{\unicode[STIX]{x1D70C}_{p}}{\unicode[STIX]{x1D70C}}\frac{\unicode[STIX]{x1D6F7}_{v}}{1-\unicode[STIX]{x1D6F7}_{v}}, & & \displaystyle\end{eqnarray}$$

ranges from $\unicode[STIX]{x1D6F7}_{m}=0.1$ to $10$ . Most of the statistics reported throughout § 4 will concentrate on $\unicode[STIX]{x1D6F7}_{m}=0$ , $1$ and $10$ for clarity. However, results for intermediate mass loadings will be included to highlight intermediate behaviour especially regarding sound intensity and TKE reduction.

Table 1.

Simulation configurations and corresponding initial fluid–particle parameters. The superscript ‘ $\bullet$ ’ corresponds to the particle introduction time when the shear layer growth reaches $\unicode[STIX]{x1D6FF}_{m}/\unicode[STIX]{x1D6FF}_{m}^{o}=10$ . The Kolmogorov length scale at this time ( $\unicode[STIX]{x1D702}^{\bullet }$ ) is computed from midplane ( $y=0$ ) turbulence data.

3.2 Particle-phase description

In this work we consider monodisperse, spherical, rigid particles with diameters smaller than the Kolmogorov length scale. The displacement of an individual particle  $i$ is calculated according to Newton’s second law of motion by

(3.6) $$\begin{eqnarray}\displaystyle \frac{\text{d}\boldsymbol{x}_{p}^{(i)}}{\text{d}t}=\boldsymbol{v}_{p}^{(i)}, & & \displaystyle\end{eqnarray}$$

and

(3.7) $$\begin{eqnarray}\displaystyle \frac{\text{d}\boldsymbol{u}_{p}^{(i)}}{\text{d}t}=\frac{\boldsymbol{f}_{drag}^{(i)}}{m_{p}}-\frac{1}{\unicode[STIX]{x1D70C}_{p}}\unicode[STIX]{x1D735}p[\boldsymbol{x}_{p}^{(i)}]+\frac{1}{\unicode[STIX]{x1D70C}_{p}}\unicode[STIX]{x1D735}\boldsymbol{\cdot }\unicode[STIX]{x1D749}[\boldsymbol{x}_{p}^{(i)}], & & \displaystyle\end{eqnarray}$$

where $\boldsymbol{x}_{p}^{(i)}$ and $\boldsymbol{v}_{p}^{(i)}$ are the instantaneous position and velocity of the $i$ th particle, respectively, and $m_{p}$ is the particle mass. Here, the fluid-phase velocity, pressure gradient and viscous stress tensor are taken at the centre position of particle $i$ . The particle equations are non-dimensionalized using the same reference quantities used in equations (2.1)–(2.3). The drag force is given by

(3.8) $$\begin{eqnarray}\displaystyle \frac{\boldsymbol{f}_{drag}^{(i)}}{m_{p}}=\frac{\unicode[STIX]{x1D6FC}}{\unicode[STIX]{x1D70F}_{p}}(\boldsymbol{u}[\boldsymbol{x}_{p}^{(i)}]-\boldsymbol{v}_{p}^{(i)})F(\unicode[STIX]{x1D6FC},Re_{p}), & & \displaystyle\end{eqnarray}$$

where $F$ is the dimensionless drag force coefficient of Tenneti, Garg & Subramaniam (Reference Tenneti, Garg and Subramaniam2011), which has a nonlinear dependence on the particle Reynolds number $Re_{p}=\unicode[STIX]{x1D70C}d_{p}Re_{c}\Vert \boldsymbol{u}-\boldsymbol{v}_{p}\Vert /\unicode[STIX]{x1D707}$ .

This drag coefficient reduces to the classic Schiller & Naumann (Reference Schiller and Naumann1933) correlation in the limit of small particle concentration, and reduces to Stokes drag as $Re_{p}$ approaches $0$ . For all simulations presented in this work, the particle Mach number and Knudsen number were found to be small (i.e. ${<}0.05$ for particle Mach number), and thus compressibility effects on the drag force (e.g. Clift, Grace & Weber Reference Clift, Grace and Weber2005) can be neglected. Other models for particle motion in compressible, inhomogeneous flows have been analysed (e.g. Parmar, Haselbacher & Balachandar Reference Parmar, Haselbacher and Balachandar2012), which include important physical mechanisms for which compressible hydrodynamics couple to particle scales. Based on the current shear layer configuration, the primary energetic pressure wavelengths, for example, are $O(10^{3})$ greater than the particle diameter, suggesting a scale separation (see figures 8 and 14). From high-speed shear layer turbulence, intense sound waves steepen (e.g. Lighthill (Reference Lighthill, Batchelor and Davies1956)) and form weak shocks that have shorter wavelength. This effect shrinks the scale separation, and the particle-phase description would need modified. However, in the current simulations, particles reside in the shear layer turbulence, $|y|\lesssim 3\unicode[STIX]{x1D6FF}_{m}$ (see figure 9) and such weak shock wave and particle interactions are absent.

The Basset force and added mass effects are neglected due to the large density ratios considered, which is appropriate when the particle diameter is significantly smaller than the acoustic wavelengths (Cleckler, Elghobashi & Liu Reference Cleckler, Elghobashi and Liu2012). It should also be noted that the current numerical approach precludes exact particle–acoustic interaction (e.g. acoustic scattering by solid bodies) since boundary conditions (e.g. no-slip condition) are not applied at the particle surface. Omitting these mechanisms are not anticipated to have large effect on the current results. Again, this is supported, to a degree, by the scale separation between particle diameter and hydrodynamic scales. For heated jets or shear layers, which are not considered here, other mechanisms like droplet evaporation and heat transfer between phases may contribute to changes in turbulence and sound radiation. However, the current configuration considers equal temperature between particles and the ambient gas phase; thus, these effects are also anticipated to be small and are therefore neglected.

Finally, all the reported results neglect particle collision effects. In the small volume fraction limit $\unicode[STIX]{x1D6FC}_{p}\ll 1$ , this approximation is reasonable, and it has been applied in similar shear layer turbulence configurations (Dai et al. Reference Dai, Jin, Luo and Fan2018). However, for the largest mass loading considered here, $\unicode[STIX]{x1D6F7}_{m}=10$ with $\unicode[STIX]{x1D6F7}_{v}=O(10^{-2})$ , collisions are abundant. Despite this, numerical tests, including particle collisions, indicate a negligible statistical effect on the turbulence and near-field pressure radiation for coefficient of restitution in the range $0.2\leqslant \text{e}\leqslant 0.8$ ; thus, our conclusions appear insensitive to particle collisions. Different injection scenarios, higher Stokes numbers or a mixture of these with particle collisions may reveal a sound radiation sensitivity.

3.3 Numerical implementation

Spatial derivatives in equations (2.1)–(2.5) are approximated by high-order finite-difference operators that satisfy the summation-by-parts (SBP) property (Strand Reference Strand1994). An explicit, sixth-order, centred finite difference is used in the domain interior, and third-order, one-sided finite differences are applied at the boundary. The SBP scheme is combined with the simultaneous-approximation-term (SAT) approach (Svärd, Carpenter & Nordström Reference Svärd, Carpenter and Nordström2007; Vishnampet et al. Reference Vishnampet, Bodony and Freund2015) at the lateral domain boundaries which provides provable stability. The outflow SAT at $\pm y=L_{y}/2$ employs a characteristic boundary condition that weakly enforces a stationary target solution $\boldsymbol{Q}_{target}(y)=\boldsymbol{Q}(y,t=0)$ , where $\boldsymbol{Q}=[\unicode[STIX]{x1D6FC}\unicode[STIX]{x1D70C},\unicode[STIX]{x1D6FC}\unicode[STIX]{x1D70C}\boldsymbol{u},\unicode[STIX]{x1D6FC}\unicode[STIX]{x1D70C}E]^{\mathsf{T}}$ is the vector of fluid-phase conserved variables. To assist in the asymptotic convergence of the SAT condition and further avoid spurious reflections into the domain, an additional sponge zone (Freund Reference Freund1997) of width $w=50\unicode[STIX]{x1D6FF}_{m}^{o}$ is added to the flow equations (denoted ${\mathcal{N}}$ for brevity) according to

(3.9) $$\begin{eqnarray}\displaystyle {\mathcal{N}}(\boldsymbol{Q})=-\unicode[STIX]{x1D70E}(y)[\boldsymbol{Q}-\boldsymbol{Q}_{target}(y)]. & & \displaystyle\end{eqnarray}$$

The sponge strength increases quadratically toward the domain extent by

(3.10) $$\begin{eqnarray}\displaystyle \unicode[STIX]{x1D70E}(y)=\left(\frac{|y|-(L_{y}/2-w)}{w}\right)^{2}. & & \displaystyle\end{eqnarray}$$

Second derivative approximations apply first-order derivative operators consecutively, which necessitates using artificial dissipation to damp the highest wavenumber energetics supported by the grid. To this end, high-order accurate SBP dissipation operators are used to provide artificial dissipation which are based on the third derivative (Mattsson, Svärd & Nordström Reference Mattsson, Svärd and Nordström2004; Vishnampet Reference Vishnampet2015). The fluid and particle equations, (2.1)–(2.3) and (3.6)–(3.7), are advanced in time with a constant time step using a standard fourth-order explicit Runge–Kutta scheme. The acoustic Courant–Friedrichs–Lewy (CFL) number was monitored throughout the simulation and remained below $\text{CFL}<0.5$ .

Fluid quantities appearing in equation (3.8) are interpolated to the location of each particle via trilinear interpolation. The interphase exchange terms appearing in equations (2.1) and (2.3) are computed by projecting the Lagrangian data onto the computational grid according to

(3.11) $$\begin{eqnarray}\displaystyle & \displaystyle \unicode[STIX]{x1D6FC}=1-\unicode[STIX]{x1D6FC}_{p}=1-\mathop{\sum }_{i=1}^{N_{p}}{\mathcal{G}}(|\boldsymbol{x}-\boldsymbol{x}_{p}^{(i)}|)V_{p}, & \displaystyle\end{eqnarray}$$

(3.12) $$\begin{eqnarray}\displaystyle & \displaystyle \boldsymbol{F}=-\mathop{\sum }_{i=1}^{N_{p}}{\mathcal{G}}(|\boldsymbol{x}-\boldsymbol{x}_{p}^{(i)}|)\boldsymbol{f}_{drag}^{(i)}, & \displaystyle\end{eqnarray}$$

and

(3.13) $$\begin{eqnarray}\displaystyle \boldsymbol{u}_{p}\boldsymbol{\cdot }\boldsymbol{F}=-\mathop{\sum }_{i=1}^{N_{p}}{\mathcal{G}}(|\boldsymbol{x}-\boldsymbol{x}_{p}^{(i)}|)\boldsymbol{v}_{p}^{(i)}\boldsymbol{\cdot }\boldsymbol{f}_{drag}^{(i)}, & & \displaystyle\end{eqnarray}$$

where ${\mathcal{G}}$ is a filter kernel, $N_{p}$ is the total number of particles and $V_{p}$ is the particle volume. In this work, ${\mathcal{G}}$ is taken to be Gaussian with a characteristic size $\unicode[STIX]{x1D6FF}_{f}=10d_{p}$ , defined as the full width at half the height of the kernel. For computational efficiency, the filtering procedure is solved in two steps (Capecelatro & Desjardins Reference Capecelatro and Desjardins2013). First, the particle data are transferred to the nearest neighbouring cells via trilinear extrapolation. The data are then diffused such that the final width of the filtering kernel is independent of the mesh size. To ensure unconditional stability and reduce computational cost, the diffusion process is solved implicitly during each Runge–Kutta sub-iteration by utilizing the approximate factorization scheme of Briley & McDonald (Reference Briley and McDonald1977). When the particle diameter is sufficiently smaller than the grid spacing ( $d_{p}<\unicode[STIX]{x0394}x/10$ ), the diffusion process is not considered and the projection method reverts to trilinear extrapolation. Details on the interphase exchange process can be found in (Capecelatro & Desjardins Reference Capecelatro and Desjardins2013).