2. Simulation methodology

In the present investigation, fluid phase is considered to be incompressible and described by the continuity and Navier–Stokes equation as

(2.1)$$\begin{gather} \frac{\partial \boldsymbol{u}_i}{\partial \boldsymbol{x}_i}=0, \end{gather}$$

(2.2) $$\begin{gather}\frac{\partial \boldsymbol{u}_i}{\partial t} + \frac{\partial\boldsymbol{u}_i \boldsymbol{u}_j}{\partial \boldsymbol{x}_{j}} ={-} \frac{1}{\rho_f} \frac{\partial p}{\partial \boldsymbol{x}_i} + {\nu} \frac{\partial^2 \boldsymbol{u}_i}{\partial \boldsymbol{x}_{j} \partial \boldsymbol{x}_{j}} + \frac{\boldsymbol{f}_{i}({\boldsymbol x},t)}{\rho_f}, \end{gather}$$

where $\boldsymbol{u}_i$ is the velocity, $p$ is the pressure, $\rho _f$ is the density of the fluid and $\nu$ is the kinematic viscosity. The feedback force density due to the solid phase is considered using the $f({\boldsymbol x},t)$ term. The pseudo-spectral method has been used to solve the Navier–Stokes equation. A second-order Adams–Bashforth scheme for the nonlinear term and Crank–Nicholson time discretization has been used for the linear terms. The lift and drag forces are included in the feedback force term, which can be expressed as

(2.3) \begin{equation} \boldsymbol{f}_{i}({\boldsymbol x},t) ={-} \sum_I (\boldsymbol{F}_{i,I}^D + \boldsymbol{F}_{i,I}^L) \delta ({\boldsymbol x}-{\boldsymbol x}_I), \end{equation}

where ${\boldsymbol x}$ is the fluid node, ${\boldsymbol x}_I$ is the position of the $I$th particle, $\boldsymbol {F}_{i,I}^L$ and $\boldsymbol {F}_{i,I}^D$ are the lift and drag forces on the particle $I$. Additionally, $\delta ({\boldsymbol x}-{\boldsymbol x}_I)$ is the Dirac delta function in three dimensions. The fluid phase momentum equation is not corrected including the fluid volume fraction term as the particle volume fraction is much lower compared with the continuum phase volume fraction ($\phi _v^{air}\approx 1$) (Dritselis & Vlachos Reference Dritselis and Vlachos2008, Reference Dritselis and Vlachos2011; Richter Reference Richter2015; Vreman Reference Vreman2015; Ghosh & Goswami Reference Ghosh and Goswami2022a,Reference Ghosh and Goswamib).

In the present study, the point-particle approach (Bagchi & Balachandar Reference Bagchi and Balachandar2003; Mehrabadi et al. Reference Mehrabadi, Horwitz, Subramaniam and Mani2018) is considered with drag and lift corrections as discussed in detail by Muramulla et al. (Reference Muramulla, Tyagi, Goswami and Kumaran2020). The point particles are tracked in the Lagrangian frame, and Newton's second law describes their motion. The particle–wall and particle–particle collisions have also been considered. The particle motion is described by

(2.4) \begin{equation} m_p \frac{{\rm d} \boldsymbol{v}_{i,I}}{{\rm d}t} = \boldsymbol{F}_{i,I}^{D} + \boldsymbol{F}_{i,I}^{L} + \sum_{I\neq J} \boldsymbol{F}_{i,IJ} + \boldsymbol{F}_{i,Iw} + m_p \boldsymbol{g},\end{equation}

where $m_p$ is the mass, $\boldsymbol {v}_{i,I}$ is the velocity of the $I$th particle, $\boldsymbol {F}_{i,I}^{D}$ is the drag force and $\boldsymbol {F}_{i,I}^{L}$ is the lift force exerted on the particle. In (2.4), $\boldsymbol {g}$ is the gravitational acceleration, $\boldsymbol {F}_{i,IJ}$ is the force due to interaction between the $I$th and $J$th particles, and $\boldsymbol {F}_{i,Iw}$ is the force due to interaction between the $I$th particle and wall. In the present study, the effect of gravity and lift are included as it is reported that the implementation of the lift force improves the particle statistics in the near-wall region (Marchioli, Picciotto & Soldati Reference Marchioli, Picciotto and Soldati2007; Costa et al. Reference Costa, Brandt and Picano2020). Marchioli et al. (Reference Marchioli, Picciotto and Soldati2007) explored the effect of gravity and lift on the particle distribution in wall-bounded flows. The hard-sphere approach is taken to account the the particle–particle and particle–wall elastic collisions. The inertia corrected drag law (Naumann & Schiller Reference Naumann and Schiller1935) is used to calculate the drag force as

(2.5)\begin{equation} \boldsymbol{F}_{i,I}^D = 3 {\rm \pi}\mu d_p (\tilde{\boldsymbol{u}}_{i,I}(\boldsymbol{x},t) - \boldsymbol{v}_{i,I}) (1+ 0.15Re_p^{0.687}).\end{equation}

Although a point-particle approximation has been used, the grid size very near the wall in the wall-normal direction may be smaller than the particle size in the near-wall region. Thus, the fraction of particle surface present in the cell is used to calculate the force on the grid. The fluid velocity is interpolated at the particle location to calculate the drag and lift. The detailed simulation procedure for the calculation of the feedback force, corrections to obtain the undisturbed fluid velocity field at the particle location, and the near-wall corrections in lift and drag force have been discussed in our earlier work (Muramulla et al. Reference Muramulla, Tyagi, Goswami and Kumaran2020). The implementation of Saffmann lift and correction for undisturbed velocity are important to accurately predict the spatial particle distribution (Mehrabadi et al. Reference Mehrabadi, Horwitz, Subramaniam and Mani2018; Wang et al. Reference Wang, Fong, Coletti, Capecelatro and Richter2019; Costa et al. Reference Costa, Brandt and Picano2020; Brandt & Coletti Reference Brandt and Coletti2022), which have been included in this study. The particle-to-fluid density ratios considered in the present study are $\approx$1000 or higher. Therefore, the buoyancy and Basset history effects are neglected in the particle's equation of motion.

The simulations have been performed in a vertical channel with length $8{\rm \pi} \delta$ in the streamwise ($x$) direction, $2\delta$ in the wall normal ($y$) direction and $(4/3){\rm \pi} \delta$ in spanwise ($z$) direction, where $\delta$ is the half-channel width. Wall normal direction is wall bounded, streamwise and spanwise directions are considered as periodic. No-slip boundary conditions are applied on the walls. The bulk Reynolds numbers ($Re_b = \rho _f \times \bar {u} \times 2\delta /\mu _f$) are fixed at 3300 and 5600 based on the channel width ($2\delta$) and average fluid velocity ($\bar {u}$). For Reynolds numbers of 3300 and 5600, 128 and 64, and 192 and 160 Fourier modes are used in the streamwise and spanwise directions, respectively (Muramulla et al. Reference Muramulla, Tyagi, Goswami and Kumaran2020). For wall-normal direction, 65 and 129 Chebyshev modes are used for $Re_b = 3300$ and 5600, respectively. The corresponding Reynolds numbers ($Re_c$) based on the centreline turbulent velocity and half-channel width are 2000 and 3360, and the Reynolds numbers ($Re_\tau$) based on unladen frictional velocity and half-channel width are 115 and 180. It is worth mentioning that the present investigation is performed under the fully developed turbulent flow condition as the transitional regime for the channel flow occurs in the range $1300\leq Re_b \leq 1800$ and $62.5\leq Re_\tau \leq 73.5$ (Patel & Head Reference Patel and Head1969; Carlson, Widnall & Peeters Reference Carlson, Widnall and Peeters1982; Sano & Tamai Reference Sano and Tamai2016; Zhang Reference Zhang2017). Thus, the simulations in our work are carried out at $Re_b$ of 1.8 and 3 times the upper limit of the transitional regime. The domain lengths ($L_x^+ \times L_y^+ \times L_z^+$) are $2921\times 232 \times 486$ and $4562 \times 363 \times 760$ for Reynolds numbers of 3300 and 5600, respectively. The ($^+$) symbol is used to indicate that the quantities are normalized with viscous scales. A time step of $0.0033h/\bar {u}$ is used in the simulations, where $h(=2\delta )$ is the channel width and $\bar {u}$ is the bulk velocity of the fluid phase. The pressure gradient is adjusted to maintain a constant bulk flow rate. The range of particle Stokes numbers is shown in table 1. In the present work, the particle diameter ($d_p$) is taken as $39\ \mathrm {\mu }$m and the particle density ($\rho _p$) is changed to have a range of Stokes number. The fluid phase is considered as air, where the fluid density ($\rho _f$) and dynamic viscosity ($\mu _f$) are 1.179 kg m$^{-3}$ and $1.75\times 10^{-5}\ {\rm kg}\ {\rm m}\ {\rm s}^{-1}$, respectively. The ratio of channel width ($2\delta$) to particle diameter ($d_p$) is 54. Simulations are performed for a range of Stokes numbers ($St = 52.73\unicode{x2013}210.93$). Various authors have reported simulations and experiments for Stokes numbers in this range (Li et al. Reference Li, McLaughlin, Kontomaris and Portela2001; Hwang & Eaton Reference Hwang and Eaton2006; Vreman et al. Reference Vreman, Geurts, Deen, Kuipers and Kuerten2009; Goswami & Kumaran Reference Goswami and Kumaran2011; Capecelatro, Desjardins & Fox Reference Capecelatro, Desjardins and Fox2015; Vreman Reference Vreman2015; Muramulla et al. Reference Muramulla, Tyagi, Goswami and Kumaran2020).

Table 1. The Stokes number ($St$) is defined as $St = \tau _p/\tau _f$, where $\tau _p = \rho _p d_p^2/18\mu _f$, $\tau _f = 2\delta /\bar {u}$, $\rho _p$ is the particle density, $d_p$ is the particle diameter, $\mu _f$ is the fluid dynamic viscosity, $\delta$ is the half channel width and $\bar {u}$ is the average fluid velocity. Here, $Re_b$ is the fluid bulk Reynolds number.

4. New modelling approach based on modified energetics

From the above analysis of the compensated spectra and second-order velocity structure function, it is observed that there is a non-monotonic decrease in the apparent Kolmogorov constant across the channel in the presence of dispersed particles. However, a significant variation of the apparent Kolmogorov constant, especially in the channel centre, happens constantly with an increase in particle volume loading. This observation can be used to model the effect of particles on the fluid phase without actually adding the discrete solid phase. The use of the Kolmogorov constant appears in the classical Smagorinsky model (Smagorinsky Reference Smagorinsky1963). In the Smagorinsky model (LES approach), the eddy viscosity term models the subgrid-scale dissipation. This eddy viscosity is represented in terms of Smagorinsky coefficient which is calculated from the Kolmogorov constant (Sagaut Reference Sagaut2006). The filtered continuity and momentum equations for LES are written as

(4.1)$$\begin{gather} \frac{\partial \tilde{\boldsymbol{u}}_i}{\partial \boldsymbol{x}_i}=0, \end{gather}$$

(4.2) $$\begin{gather} \frac{\partial \tilde{\boldsymbol{u}}_i}{\partial t} + \frac{\partial \tilde{\boldsymbol{u}}_i\tilde{\boldsymbol{u}}_j}{\partial \boldsymbol{x}_j} ={-} \frac{1}{\rho_f} \frac{\partial \tilde{p}}{\partial \boldsymbol{x}_i} + {\nu} \frac{\partial^2 \tilde{\boldsymbol{u}}_i}{\partial \boldsymbol{x}_j \partial \boldsymbol{x}_j } + \frac{\partial (\tilde{\boldsymbol{u}}_i \tilde{\boldsymbol{u}}_j - \widetilde{\boldsymbol{u}_i \boldsymbol{u}_j})}{\partial \boldsymbol{x}_j}. \end{gather}$$

Here, $\tilde {p}$ is the filtered pressure, $\tilde {\boldsymbol {u}}_i$ is the filtered velocity, $\nu$ is the kinematic viscosity and $\rho _f$ is the fluid density. The subgrid scale (SGS) stress term, ($\tilde{\boldsymbol {u}}_i \tilde {\boldsymbol {u}}_j - \widetilde {\boldsymbol {u}_i \boldsymbol {u}_j}$), in the Smagorinksy model is expressed as

(4.3)$$\begin{gather} -\boldsymbol{\tau}_{ij} = \tilde{\boldsymbol{u}}_i \tilde{\boldsymbol{u}}_j - \widetilde{\boldsymbol{u}_i \boldsymbol{u}_j} = 2 \nu_t \tilde{\boldsymbol{S}}_{ij}, \end{gather}$$

(4.4)$$\begin{gather}\tilde{\boldsymbol{S}}_{ij} = \frac{1}{2} \left[ \frac{\partial \tilde{\boldsymbol{u}}_i }{ \partial \boldsymbol{x}_j} + \frac{\partial \tilde{\boldsymbol{u}}_j }{ \partial \boldsymbol{x}_i} \right], \end{gather}$$

where $\nu _t$ is the eddy viscosity and $\tilde {\boldsymbol {S}}_{ij}$ is the filtered strain rate tensor. The eddy viscosity is written as

(4.5)$$\begin{gather} \nu_t = (C_s \tilde{\triangle})^2 |\tilde{S}| , \end{gather}$$

(4.6)$$\begin{gather}|\tilde{S}| = \sqrt{2 \tilde{\boldsymbol{S}}_{ij} \tilde{\boldsymbol{S}}_{ij}}. \end{gather}$$

In (4.5), $C_s$ is the Smagorinsky coefficient, $|\tilde {S}|$ is the magnitude of the strain rate and $\tilde {\triangle } = (\tilde {\triangle }_1 \tilde {\triangle }_2 \tilde {\triangle }_3 )^{1/3}$ is the cube-root volume of the grid size. Here, $\tilde {\triangle }_1$, $\tilde {\triangle }_2$ and $\tilde {\triangle }_3$ are the grid spacing in the $x$, $y$ and $z$ directions, respectively. The relation between the Kolmogorov constant and the Smagorinsky coefficient is given as (Sagaut Reference Sagaut2006)

(4.7)\begin{equation} C_s = \frac{0.23}{C^{3/4}}.\end{equation}

Therefore, the variation of the Kolmogorov constant ($C$) causes a modification in the Smagorinsky coefficient $C_s$. The apparent Kolmogorov constant decreases with an increase in the particle volume loading, which results in an increased $C_s$. Thus, the simulations are performed for single-phase flow with higher $C_s$ values for the Smagorinsky model without adding the particles. In this context, it is worth noting the earlier work of Yeo et al. (Reference Yeo, Dong, Climent and Maxey2010) who performed the fully resolved simulation for bubbles, neutrally buoyant and inertial particles in homogenous turbulence, and commented on the possibility of representing the particle feedback effect with an additional effective viscosity.

The present exercise aims to capture the effect of particles using a single-phase simulation by varying the Smagorinsky coefficient without adding the particles to the system. Although we have observed that the variation of the apparent Kolmogorov constant is a function of wall-normal location, as a first approximation, we consider a single Kolmogorov constant approximation across the channel as a function of particle concentration. Use of a uniform concentration across the channel would be a good assumption for the range of Stokes number used here, which is also evident from the DNS results, shown later. It should be noted that the particles increase the anisotropy across the fluid fluctuations and an anisotropic inhomogeneous modelling approach should be taken to model the fluid phase accurately, which is left as a future scope. However, the present method of modelling the dynamics of the fluid phase in a two-phase flow using a single-phase simulation is computationally less expensive. It provides insights into the mechanism of turbulence modulation. Simulations are performed with the Smagorinksy model which has been used in our earlier work (Rohilla et al. Reference Rohilla, Muramulla and Goswami2022). The number of grids used is $128\times 65\times 64$ and $64\times 65\times 32$ in the streamwise ($x$), wall-normal ($y$) and spanwise ($z$) directions for bulk Reynolds numbers of 5600 and 3300. Here, $C_s = 0.125$ is used in the unladen flow simulations and the Van Driest damping is implemented to avoid high dissipation in the near-wall region (Rohilla et al. Reference Rohilla, Muramulla and Goswami2022).

First, we have increased the value of $C_s$ in LES to determine the effect of $C_s$ on the fluid phase fluctuations. Then, we have presented the variation of second moments of fluctuation simultaneously as a function of $C_s$ and particle volume loading ($\phi _{av}$). In figure 9, the temporal evolution of the normalized sum of the fluid fluctuations is plotted for $Re_b = 5600$ with different $C_s$ values. The simulation with a lower $C_s$ value reaches a stationary state quickly. For high $C_s$ such as 0.30, the sum decreases initially and reaches the stationary state after a long time. For $C_s = 0.35$ and higher values, the sum of fluid fluctuations continuously decreases, representing a decay of the intensity of fluid turbulence. In figure 10, the profiles of fluid fluctuations are plotted along the wall-normal direction for a range of $C_s$. It is observed that the streamwise fluid fluctuations initially increase with an increase in $C_s$, and a sudden decrease happens at $C_s = 0.35$. However, a continuous decrease is observed for Reynolds stress, wall-normal and spanwise fluid velocity fluctuations, and a complete collapse of turbulence is observed at $C_s = 0.35$. As the transverse fluctuations are decreased, there is a decrease in momentum flux. The non-monotonic behaviour in the streamwise fluid fluctuations has also been reported by Zhou et al. (Reference Zhou, Zhao, Huang and Xu2020) where the authors performed the DNS of particle-laden channel flow. An increase in the streamwise fluctuations was observed at a solid volume fraction of $1.8\times 10^{-4}$, and a decrease was observed with a further increase in volume fraction. The authors mentioned two competitive phenomena. First, the near-wall vortices become weaker, leading to the larger spacing between the streaks. This effect causes a decrease in streamwise fluctuations. However, in the second case, the streaks become more organized and aligned in the streamwise direction. This phenomenon increases the streamwise fluctuations. A continuous decrease in the wall-normal and spanwise fluctuations was observed by Zhou et al. (Reference Zhou, Zhao, Huang and Xu2020). In the present study, an increase in $C_s$ leads to an initial increase in the streamwise fluid fluctuations. Then, a sudden turbulence collapse is observed at a higher $C_s = 0.35$. It is to be noted that we have used a constant value of $C_s$ across the channel. A further analysis using inhomogeneous $C_s$ is expected to provide more quantitative modulation of fluid turbulence. The present simulation shows a monotonic decrease in the Reynolds stress, wall-normal and spanwise fluid fluctuations (figure 10b–d), as reported by DNS studies (Li et al. Reference Li, McLaughlin, Kontomaris and Portela2001; Vreman Reference Vreman2015; Muramulla et al. Reference Muramulla, Tyagi, Goswami and Kumaran2020; Zhou et al. Reference Zhou, Zhao, Huang and Xu2020).

Figure 9. The average of fluid fluctuations (normalized with the $t^* = 200$ value) across the channel evolved over the time for a range of $C_s$. The $t^*$ is normalized with channel width and fluid bulk velocity.

Figure 10. The profiles of fluid fluctuations in the wall-normal direction. The simulations are performed with Smagorinksy model with increasing Smagorinsky coefficient ($C_s$) for $Re_b = 5600$. (a) Streamwise fluctuations, (b) wall-normal fluctuations, (c) Reynolds stress and (d) spanwise fluctuations. The symbols in panels (b–d) are the same as those in panel (a).

Figure 10 shows that the Reynolds stress and other components of fluctuations are a strong function of wall-normal distance. Therefore, we define the channel-averaged second moments of velocity fluctuations (4.8) and present the effect of variation of $C_s$:

(4.8)\begin{equation} \langle \star \rangle_s = \frac{1}{\delta} \int_{0}^{\delta} {{\rm d} y} \langle \star \rangle.\end{equation}

Here, $\langle \star \rangle _s$ is the averaged quantity over half-channel width ($\delta$). The average fluid fluctuations across the channel as a function of $C_s$ are plotted in figure 11(b,d) for both the Reynolds numbers. For DNS, the average fluid fluctuations as a function of $\phi _{av}$ for both the Reynolds numbers are shown in figure 11(a,c). The Stokes numbers considered are 105.47 and 210.93 for $Re_b = 3300$ and 5600, respectively, for DNS. The average fluid fluctuations decreases with an increase in particle volume fraction, and a complete turbulence collapse is observed at $\phi _{av} = 10^{-3}$ and $2.8\times 10^{-3}$ for $Re_b = 3300$ and 5600, respectively, shown in figure 11(a,c). It is observed that the fluid fluctuations, except the streamwise fluctuations, decrease with an increase in $C_s$ value, and a drastic collapse of more than one order of magnitude is observed at $C_s = 0.20$ for $Re_b = 3300$ and at $C_s = 0.35$ for $Re_b = 5600$, shown in figure 11(b,d). The decrease in Reynolds stress, wall-normal and spanwise fluctuations (figure 11b,d) with an increase in $C_s$ is similar to the decrease observed for DNS (figure 11a,c). The complete turbulence collapse observed in fluid fluctuations is similar to the collapse observed for particle-laden cases by Mito & Hanratty (Reference Mito and Hanratty2006), Shringarpure et al. (Reference Shringarpure, Cantero and Balachandar2012), Capecelatro et al. (Reference Capecelatro, Desjardins and Fox2018), Kumaran et al. (Reference Kumaran, Muramulla, Tyagi and Goswami2020), Duque-Daza et al. (Reference Duque-Daza, Ramirez-Pastran and Lain2021), Muramulla et al. (Reference Muramulla, Tyagi, Goswami and Kumaran2020) and Rohilla et al. (Reference Rohilla, Muramulla and Goswami2022). However, once turbulence is collapsed, the fluid fluctuations do not remain constant or increase as happens in the case of particle-laden cases due to particle-induced fluctuations (Capecelatro et al. Reference Capecelatro, Desjardins and Fox2018; Kumaran et al. Reference Kumaran, Muramulla, Tyagi and Goswami2020; Yu et al. Reference Yu, Xia, Guo and Lin2021). The volume fractions ($\phi _{eq}$) which correspond to different $C_s$ values are plotted in figure 12, where fluid fluctuations (shown in figure 11) predicted by Smagorinksy model are compared with DNS. The Kolmogorov constant, computed using the second-order velocity structure function and compensated spectra, is found to decrease linearly in the channel centre location, shown in figures 6(b) and 7(d). The equivalence between the $C_s$ (used as the Smagorinsky coefficient) and the particle volume loading has been derived as follows. From the value of $C_s$ (0.125) used for the unladen flow, $C$ is calculated using (4.7). It is assumed that the variation of $C$ as a function of $\phi _{av}$ follows the same functional form as that of $C_2$ versus $\phi _{av}$, shown in figure 6(b). Therefore, we can estimate $C$ and $C_2$ for a range of particle volume loading ($\phi _{av}$). The equivalent volume fraction ($\phi _{eq}$) which corresponds to $C_s$ is plotted in figure 12 along with the average volume fraction ($\phi _{av}$). For $Re_b = 3300$ and 5600, the turbulence collapse was observed at $\phi _{av} = 10^{-3}$ and $2.8\times 10^{-3}$ , respectively, for particle-laden DNS (Muramulla et al. Reference Muramulla, Tyagi, Goswami and Kumaran2020). However, the turbulence collapse for the present simulations is observed at $\phi _{eq} = 1.4\times 10^{-3}$ and $3.2\times 10^{-3}$ for $Re_b = 3300$ and 5600, which is very close the critical loading predicted by DNS.

Figure 11. The average fluid fluctuations normalized with fluid bulk velocity ($\bar {u}$) in (a,b) $Re_b = 3300$ and (c,d) $Re_b = 5600$. (a,c) The average fluid fluctuations are shown for particle-laden DNS (Muramulla et al. Reference Muramulla, Tyagi, Goswami and Kumaran2020). (b,d) The average fluid fluctuations are shown for the Smagorinksy model over a range of Smagorinsky coefficient ($C_s$). The dashed lines with closed symbols are for the Smagorinksy model and the dashed lines with open symbols are for the DNS. The legend is as follows: black line, streamwise; red line, Reynolds stress; blue line: wall-normal; and magenta line, spanwise fluctuations.

Figure 12. The average fluid fluctuations (shown in figure 11) for the Smagorinsky model and DNS are compared in (a,b) $Re_b = 3300$ and (c,d) $Re_b = 5600$. Here, $C_s$ is represented in terms of equivalent volume fraction ($\phi _{eq}$). (a,c) Streamwise fluid fluctuations and (b,d): cross-stream, wall-normal and spanwise fluid fluctuations. The dashed lines with closed symbols are for the Smagorinksy model and the dashed lines with open symbols are for the DNS. The legend is as follows: black line, streamwise; red line, Reynolds stress; blue line, wall-normal; and magenta line, spanwise fluctuations.

In the present study, an attempt is made to capture the effect of change in particle volume loading in a two-phase flow by a fluid phase only simulation with a modified Smagorinsky coefficient. Therefore, it is expected that the total dissipation caused by the viscous term and modified eddy viscosity term in the present simulation should be similar to the total dissipation by the mean viscous term and dissipation due to feedback force exerted by the particles. Therefore, we have compared the total dissipation predicted in these two cases. The kinetic energy balance equation of the mean fluid flow in the stationary state for the particle-laden flow can be described as

(4.9)\begin{align} &- {U}_x \frac{1}{\rho } \frac{\partial {P}}{\partial x} - \frac{\partial (\overline{{u_x^{'}} {u_y^{'}}} {U}_x)}{\partial y} + \overline{{u_x^{'}} {u_y^{'}}} \frac{\partial {U}_x}{\partial y} \nonumber\\ &\quad +\nu \frac{\partial}{\partial y} \left({U}_x \frac{\partial {U_x}}{\partial y}\right) - \nu \frac{\partial {U}_x}{\partial y} \frac{\partial {U}_x}{\partial y} - {U}_x \overline{ \frac{\rho_p f \phi_c}{ \rho_f \tau_p} ({u}_x - {v}_x) } = 0. \end{align}

Here, ${P}$ is the mean pressure, ${u_i}$ is the instantaneous velocity, ${U}_x$ is the mean velocity and ${u'_i}$ is the fluctuating velocity. Additionally, $\rho _f \overline {{u^{'}_i} {u^{'}_j}}$ is the Reynold stress, $\tau _p$ is the particle relaxation time, $\phi _c$ is the volume fraction in the respective grid, $f$ is the drag coefficient, $\rho _p$ is the material density of the particle and $\rho _f$ is the fluid density. In (4.9), the first term is the energy due to pressure work, the second term is the divergence of energy fluxes due to Reynolds stress, the third term is the energy used for the turbulence production, the fourth term is the divergence of energy flux due to fluid viscous stress, the fifth term is the viscous dissipation ($\epsilon _m$) due to mean shear and the sixth term is the loss of energy due to the particle drag ($F_p$). The filtered mean kinetic energy equation for unladen fluid flow is written as

(4.10)\begin{align} &- \tilde{U}_x \frac{1}{\rho } \frac{\partial \tilde{P}}{\partial x} - \frac{\partial (\overline{\widetilde{u_x^{'}} \widetilde{u_y^{'}}} \tilde{U}_x)}{\partial y} + \overline{\widetilde{u_x^{'}} \widetilde{u_y^{'}}} \frac{\partial ( \tilde{U}_x)}{\partial y}\nonumber\\ &\quad +(\nu +\nu_t) \frac{\partial}{\partial y}\left(\tilde{U}_x \frac{\partial \widetilde{U_x}}{\partial y}\right) - (\nu +\nu_t) \frac{\partial \tilde{U}_x}{\partial y} \frac{\partial \tilde{U}_x}{\partial y} = 0.\end{align}

Here, $\tilde {P}$ is the mean filtered pressure, $\tilde {u}_i$ is the instantaneous filtered velocity, $\tilde {U}_x$ is the mean filtered velocity, $\widetilde {u'_i}$ is the filtered fluctuating velocity and $\rho _f \overline {\widetilde {u^{'}_i} \widetilde {u^{'}_j}}$ is the filtered Reynold stress. The terms in (4.10) are similar to those in (4.9) with filtered quantities. The dissipation in the mean kinetic energy equation (4.9) is caused by the last two terms, mean viscous dissipation and particle-induced dissipation. However, the dissipation in (4.10) is due to the last term which is due to molecular ($\epsilon _m$) and eddy viscosity ($\epsilon _{eddy}$). Thus, these two dissipation terms should be comparable to predict fluid phase dynamics accurately when a single-phase simulation supplants the two-phase simulation.

Figures 13(a) and 13(c) show the terms due to mean viscous dissipation, dissipation due to particle drag and input energy due to pressure work for particle-laden DNS cases (Rohilla et al. Reference Rohilla, Muramulla and Goswami2022). In figure 13(b,d), the Smagorinsky coefficient ($C_s$) is varied, and input energy due to pressure work, the dissipation due to mean viscous and eddy viscosity are plotted using the Smagorinsky model for the unladen flow. It is observed that with an increase in $C_s$, the mean viscous dissipation decreases. However, the dissipation due to eddy viscosity increases. The dissipation due to particle feedback in figure 13(a,c) is comparable and shows a similar trend with the dissipation by the modified eddy viscous term in figure 13(b,d). Also, the total dissipation in the present unladen simulations is comparable to the total dissipation caused by the particle-laden DNS. It is to be noted that the total dissipation (mean viscous and particle induced dissipation) is almost constant in all the cases, which has also been reported previously (Rohilla et al. Reference Rohilla, Muramulla and Goswami2022). In the case of homogenous isotropic turbulence, it was observed by Squires & Eaton (Reference Squires and Eaton1990) that the total dissipation was constant for all the cases. For the mass loading of one, the decrease in the viscous dissipation was nearly $50\,\%$, and another $50\,\%$ dissipation was caused by the particles. A decrease in the input energy is observed for the particle-laden DNS, figure 13(a,c), and in the case of single-phase Smagorinsky simulations, figure 13(b). In figure 13(d), a decrease in the input energy is observed at a low value of $C_s$, and a non-monotonic variation is observed before the turbulence collapse. The difference between input energy and the total dissipation, which is used in the turbulence production, is plotted in figure 14, and is compared with DNS (Rohilla et al. Reference Rohilla, Muramulla and Goswami2022). Here, $C_s$ is represented in terms of equivalent volume fraction ($\phi _{eq}$). A decrease in the difference ($diff$) is observed with an increase in $\phi _{eq}$ for the low Reynolds number ($Re_b = 3300$). However, it is not monotonic at moderate Reynolds number ($Re_b = 5600$). A sudden collapse occurs at $\phi _{eq} = 0.0014$ and 0.0032 for $Re_b = 3300$ and 5600, respectively. The turbulence collapse happens at $\phi = 0.001$ and 0.0028 for DNS at $Re_b = 3300$ and 5600, respectively. The critical volume fraction predicted by the Smagorinsky model closely matches with DNS at both Reynolds numbers.

Figure 13. The terms from the mean kinetic energy equation are presented in (a) $Re_b = 3300$, $St = 105.47$, (b) $Re_b = 3300$, (c) $Re_b = 5600$, $St$ = 210.93 and (d) $Re_b = 5600$. Panels (a,c) are particle-laden DNS cases where $\epsilon _m$ is the mean viscous dissipation, $F_{p,DNS}$ is the dissipation due to particles and $sum$ is the addition of feedback and mean viscous dissipation. Panels (b,d) are single-phase simulations where the Smagorinsky coefficient is varied, $\epsilon _m$ is the viscous dissipation at mean flow, $\epsilon _{eddy}$ is the dissipation due to eddy viscosity and $sum$ is the addition of eddy and viscous dissipation. Here, $IE$ is the input energy from pressure work.

Figure 14. Difference of input energy ($IE$) and total dissipation ($sum$), which are plotted in figure 13, $diff = IE - sum$. (a) $Re_b = 3300$ and (b) $Re_b = 5600$. The subscripts are as S, Smagorinsky model and DNS, direct numerical simulation.

The feedback term from particle-laden DNS and eddy viscosity-based dissipation due to modified $C_s$ from the Smagorinksy model is compared in figure 15. Here, the dissipation due to eddy viscosity for the unladen simulation ($\nu _t ({\rm d}\tilde {U}/{{\rm d} y})^2$ for the $C_s = 0.125$ case) is subtracted from the $F_{p}$ to depict the equivalent particle dissipation only. The feedback term in DNS is taken for $St = 105.47$ and 210.93 for $Re_b = 3300$ and 5600, respectively. The feedback term as a function of volume fraction ($\phi _{av}$) from particle-laden DNS and as a function of equivalent volume fraction ($\phi _{eq}$) is plotted from the Smagorinsky model. As the turbulence collapse happens at different $\phi _{av}$ and $\phi _{eq}$ values for particle-laden DNS and Smagorinsky coefficient (figure 12), the $\phi _{av}$ and $\phi _{eq}$ are divided by the $\phi _{av,cr}$ and $\phi _{eq,cr}$, which are the critical loadings for DNS and the Smagorinsky model, respectively, where turbulence collapse is observed. The Smagorinsky model predicts the feedback term and the trend with reasonable accuracy for both Reynolds numbers. For $Re_b = 3300$ in figure 15(a), the Smagorinsky model accurately predict the dissipation at low $\phi _{eq}$, but underpredicts particle dissipation near the $\phi _{eq,cr}$. In the case of moderate Reynolds number ($Re_b = 5600$), the prediction by the Smagorinksy model matches with DNS, figure 15(b). In the case of the Smagorinsky model, it is to be noted that particles act as a source/sink in particle-laden cases depending on the local relative velocity. However, in the case of modified average $C_s$, it will be a dissipative effect only. Therefore, further analysis is required so that $C_s$ can be expressed as a function of Reynolds number, Stokes number and wall-normal distance to capture the effect of particles more accurately. This will be an interesting future scope.

Figure 15. The dissipation term in DNS due to particles and an equivalent dissipation in the Smagorinksy model due to modified $C_s$ are compared. Here, $C_s$ is represented in terms of equivalent volume fraction ($\phi _{eq}$). (a) $Re_b = 3300$ and (b) $Re_b = 5600$. The volume fraction ($\phi _{av}$ or $\phi _{eq}$) is normalized with a critical volume fraction at which turbulence collapse is observed in individual cases.

The mean velocity, the gradient of mean velocity and the pressure gradient across the channel are plotted in figure 16. The plots are shown for the simulations performed with the Smagorinsky model with varying $C_s$ for $Re_b = 5600$, and the particle-laden DNS for $Re_b = 5600$ and $St = 210.93$. The mean velocity profiles in figure 16(a,b) show that the mean velocity decreases in the buffer region and increases in the channel centre with an increase in $C_s$. However, the extent of decreases predicted by modified $C_s$ is more than that predicted by DNS, shown in figure 16(b). In the case of the mean velocity gradient, a significant decrease is observed in the near-wall region for the Smagorinsky model (figure 16c) compared with the DNS case (figure 16d). The Smagorinksy model captures the qualitative behaviour of the pressure gradient, which is observed in particle-laden DNS, figure 16(e). Here, $C_s$ is denoted in terms of the equivalent volume fraction ($\phi _{eq}$). Thus, the variation of $C_s$ in the Smagorinksy model captures the effect of particles qualitatively, and further analysis will lead to new LES models in the future.

Figure 16. The streamwise mean velocity (a) for $Re_b = 5600$, (b) $Re_b = 5600$, $St = 210.93$, mean velocity gradient (c) for $Re_b = 5600$, (d) $Re_b = 5600$, $St = 210.93$, and pressure gradient (e) for $Re_b = 5600$. The streamwise velocity in panels (a,b) are normalized with unladen frictional velocity ($u_\tau$). The mean velocity and pressure gradient in panels (c–e) are normalized with fluid bulk velocity and channel width. In panels (a,c), the simulations are performed with different $C_s$ for the Smagorinsky model. Panels (b,d) show the results from particle-laden DNS for $St = 210.93$. The pressure gradient for a range of $\phi _{av}$ and $\phi _{eq}$ is compared in panel (e) where the legend is S, Smagorinsky model and DNS, direct numerical simulation.

The particle statistics are also presented for both the Reynolds numbers, as shown in figures 17 and 18. It is worth noting that the simulations to predict the particle properties are performed as one-way coupling. Modulation of fluid properties due to the presence of particles is embedded in our new modelling technique through the modified eddy viscosity term. The particle mean velocities predicted by LES match well with the DNS results for a range of volume fractions and the velocity profiles are flat across the channel width, as shown in figure 17(a) for $Re_b = 3300$. The unladen mean fluid velocity for DNS is also plotted. It is observed that the particles lead the fluid in the near-wall region while lag in the channel centre. The particle fluctuations scaled with averaged fluid velocity fluctuation are plotted in figure 17(b,c). The streamwise particle fluctuations are higher near the wall and decrease away from the wall. The wall-normal particle fluctuations are flat across the channel width. The streamwise velocity fluctuation decreases, while wall-normal fluctuations increases as the solid loading is increased. This happens due to increased particle collision frequency and the transfer of momentum from the streamwise direction to the other directions. The relative particle fluctuations predicted by the modified Smagorinsky model (present model) deviate by a maximum of 20 % compared with the DNS results for lower volume fractions. Near the critical loading, the deviation increases as the DNS predicted fluid fluctuation is much lower than that predicted by the LES. The particle concentration scaled with the average value is plotted along the wall-normal direction. The normalized particle concentration profiles are almost flat across the channel width and the profiles predicted by the LES model agree well with the DNS results for a range of volume fractions, as shown in figure 17(d). For $Re_b = 5600$, the particle mean velocity profiles are also flat for a range of volume fractions, and the modified-Smagorinsky model predicts the mean velocities with a deviation of less than 3 % compared with the DNS results. For $Re_b = 5600$, a similar trend is observed for the particle velocity fluctuations as for $Re_b = 3300$ (figure 18). A maximum deviation of $20\,\%$ is observed. The scaled particle concentration profiles are almost flat and are in good agreement with DNS. The results for the prediction of second moments of particle velocity fluctuations may be improved if an inhomogeneity and Stokes number dependence are included in the $C_s$ prediction.

Figure 17. The particle properties are plotted in the wall-normal direction for $Re_b = 3300$ and $St = 105.47$. The particle mean velocities ($\bar {v}_x$) are scaled with fluid bulk velocity ($\bar {u}$), the particle fluctuations are scaled by average fluid fluctuations and the wall-normal distance is scaled with channel half-width ($\delta$). The $\phi _{eq}$ for the Smagorinsky model is represented by $\phi$ itself in the legends. (a) Mean velocity, (b) root mean square of streamwise, (c) wall-normal fluctuations and (d) particle concentration. The particle concentration is normalized with an average concentration across the channel width. Unladen DNS in panel (a) shows the unladen fluid mean velocity for DNS. The legend in panel (c) is the same as that in panel (b).

Figure 18. Particle properties plotted in the wall-normal direction for $Re_b = 5600$ and St = 210.93. The particle mean velocities ($\bar {v}_x$) are scaled with fluid bulk velocity ($\bar {u}$), the particle fluctuations are scaled by average fluid fluctuations and wall-normal distance is scaled with channel half-width ($\delta$). The $\phi _{eq}$ for the Smagorinsky model is represented by $\phi$ itself in the legends. (a) Mean velocity, (b) root mean square of streamwise,(c) wall-normal fluctuations and (d) particle concentration. The particle concentration is normalized with an average concentration across the channel width. Unladen DNS in panel (a) shows the unladen fluid mean velocity for DNS.