3.1. Spatial distribution of particle and fluid velocity

Figure 3.

Time evolution of fluid velocity field and particle motion. Results are shown for ( $a$ ) $t = 0.61$ , ( $b$ ) $t = 1.21$ , ( $c$ ) $t = 1.52$ , ( $d$ ) $t = 1.82$ , ( $e$ ) $t = 2.42$ , ( $f$ ) $t = 3.03$ , ( $g$ ) $t = 4.24$ , ( $h$ ) $t = 5.45$ , ( $i$ ) $t = 6.67$ and ( $j$ ) $t = 7.88$ . The slice indicates the spatial distribution of the spanwise-averaged fluid streamwise velocity $u_{\!f}$ with the grey lines representing the streamlines of the spanwise-averaged fluid velocity $\boldsymbol{u}_f$ . Particles are coloured with particle streamwise velocity $u_{{\kern-1pt}p}$ .

This section discusses the result of case 3 ( ${Re}=3030$ ) as an example. Figure 3 illustrates the time evolution of the fluid velocity field and particle behaviour. It can be seen from the figure that under the dominance of particle gravity, the particle–fluid mixture gradually collapses and evolves into a TC downstream. This development process aligns with findings from previous studies (Cantero et al. Reference Cantero, Balachandar and García2008; Espath et al. Reference Espath, Pinto, Laizet and Silvestrini2015; Xie et al. Reference Xie, Hu, Zhu, Yu and Pähtz2023a ; Lu et al. Reference Lu, Liu, Xie, Sun, Yang and Guo2024). During the entire TC evolution, the interaction between particles and fluids causes the two phases to exhibit similar motion characteristics. At the very early stage of the collapse (figure 3 $a$ ), both the fluid and the particles in the turbidity front exhibit positive flow velocities. Due to the conservation of mass, the upper ambient fluid is pumped in the opposite direction into the lock (Cantero et al. Reference Cantero, Lee, Balachandar and Garcia2007), resulting in a counterclockwise flow in the ambient fluid region. This phenomenon explains why the fluid and particle streamwise velocities in the upper TC region are initially negative.

As the collapse progresses, the local positive particle and fluid streamwise velocities ( $u_{{\kern-1pt}p}$ and $u_{\!f}$ ) at some locations near the interface within the TC region increase, as shown in figure 3( $b$ – $d$ ). This is accompanied by an increase in the negative streamwise velocity of the ambient fluid on the upper side of the TC interface. It suggests that there is a local increase in shear between the TC and the ambient fluid. Consequently, a local flow trend emerges, moving away from the collapsing body, leading to undulations in the interface morphology and a loss of its initial smoothness. At $t = 2.42$ , a distinct KH billow composed of a series of anticlockwise KH vortex pairs appears as shear instability intensifies, as seen in figure 3( $e$ ). This observation contrasts with the finding that KH billows do not form in TCs at low Reynolds numbers (Xie et al. Reference Xie, Hu, Pähtz, He and Cheng2022, Reference Xie, Hu, Zhu, Yu and Pähtz2023a ). It suggests that the larger flow Reynolds number ${\textit{Re}}$ , representing the ratio of inertial forces to viscous diffusion, promotes the development of interfacial instability of the lock-exchange TC (Ottolenghi et al. Reference Ottolenghi, Adduce, Inghilesi, Roman and Armenio2016). The KH vortex connects with the TC and entertains particles from it through the downstream side of the vortices. Some KH vortices merge and are eventually destroyed due to vortex–vortex interactions at $t = 3.03$ (figure 3 $f$ ). The primary reason is the intense shearing action brought about by the opposing vertical velocities of adjacent KH vortices (Cerretelli & Williamson Reference Cerretelli and Williamson2003). Simultaneously, as the current evolves, the reflected bore generated by the backflowing ambient fluid impacting the wall propagates forward (Shin, Dalziel & Linden Reference Shin, Dalziel and Linden2004) and eventually interacts with each KH vortex. It eventually rises upward under the obstruction of the downstream TC, as shown in figure 3( $g$ ), generating a counterclockwise vortex structure that exists for a long time – this is consistent with Lu et al. (Reference Lu, Liu, Xie, Sun, Yang and Guo2024). In the vicinity of $x=2{\;{\sim}\;} 2.5$ , it induces a clockwise vortex by forcing the particles moving upward and forward, as depicted in figure 3( $h$ , $i$ ), which was also reported in Cantero et al. (Reference Cantero, Lee, Balachandar and Garcia2007). Particles keep leaking out of the vortex in the centre of the simulation region in figure 3( $i$ , $j$ ) because of the attenuation of vortex intensity caused by energy dissipation. Furthermore, an anticlockwise vortex is still generated close to the TC head (figure 3 $g$ ), where the reflected bore is not able to catch up, leading to particle ejection in figure 3( $h$ ).

Figure 4.

Time evolution of streamwise fluid–particle slip velocity $u_{\!f} - u_{{\kern-1pt}p}$ . Results are shown for ( $a$ ) $t = 0.61$ , ( $b$ ) $t = 1.21$ , ( $c$ ) $t = 1.52$ , ( $d$ ) $t = 1.82$ , ( $e$ ) $t = 2.42$ and ( $f$ ) $t = 3.03$ .

Figure 4 displays the spatial distribution of streamwise interphase slip velocity $u_{\!f}-u_{{\kern-1pt}p}$ , which directly reflects the direction of the drag force exerted by the fluid on the particles. Using the definition of the streamwise particle Reynolds number ${Re}_{p,x} = \rho _f {\textrm{d}}_p |u_{\!f} - u_{{\kern-1pt}p}|/{\mu} _f$ of the order of $O(10^{-4}$ ), figure 4 also provides insight into the strength of the fluid–particle interaction in the streamwise direction. Streamwise slip velocity $u_{\!f} - u_{{\kern-1pt}p}$ , of less than $O(10^{-4})$ , reflects a weak deviation in streamwise velocity between the two phases, as the particle Stokes number $St$ is very small ( $St \ll 1$ ). Most areas in the TC body exhibit a low positive slip velocity, consistent with Eulerian–Eulerian results in Cantero et al. (Reference Cantero, Balachandar and García2008). In the early stages, a relatively significant positive slip velocity is observed near the left wall, while in the middle and later stages, it is primarily located near the upper interface of the middle and tail regions of the flow, which is linked to the downstream propagation of the reflected bore. The negative slip velocity is predominantly found near the upper interface, especially within the interfacial instability vortices, due to the shearing effect of the ambient fluid. This indicates that particles within these vortices experience substantial two-phase interactions.

Figure 5.

Time evolution of fluid relative pressure $p'$ and vortical coherent structures. Results are shown for ( $a$ ) $t = 0.61$ , ( $b$ ) $t = 1.21$ , ( $c$ ) $t = 1.52$ , ( $d$ ) $t = 1.82$ , ( $e$ ) $t = 2.42$ and ( $f$ ) $t = 3.03$ . The slice shows the spatial distribution of the spanwise-averaged fluid relative pressure $p'$ . Translucent grey areas indicate TCs. Vortical coherent structures are captured by $Q$ -criterion ( $Q = 5$ ) that are marked in cyan.

Figure 5 depicts the spatial distribution of the fluid relative pressure $p'$ , which is defined as the difference between the fluid total pressure $p$ and the local theoretical hydrostatic pressure. During the early stages of the TC, particle sedimentation causes the lower fluids to exhibit a high relative pressure $p'$ , while the upper layer displays a low value (Si, Shi & Yu Reference Si, Shi and Yu2018). It also affects the relative pressure of the nearby ambient fluid, which exhibits comparable distribution features. As time goes by, regions of high relative pressure shift to the TC front and adjacent areas, as well as to the bottom of the initial fluid–particle collapse region. Affixed to the upper interface of the TC, the low relative pressure centre follows the TC in the longitudinal direction (Lee Reference Lee2019). In the early stages (figure 3 $a$ – $e$ ), the low-pressure centre’s relative pressure steadily falls, and in the latter stages (figure 3 $e$ , $f$ ), it gradually increases.

Figure 5 also depicts the coherent vortical structures captured by the $Q$ -criterion. The $Q$ value is given by the symmetric tensor and antisymmetric tensor of the fluid velocity gradient $\boldsymbol\nabla \boldsymbol{u}_f$ (Hunt, Wray & Moin Reference Hunt, Wray and Moin1988; da Silva & Pereira Reference De Silva, Fernando, Eaton and Hebert2008), i.e.

(3.1) \begin{equation} Q = \frac {1}{2}\left ( {\Omega _{ij}}{\Omega _{ij}} - {{S _{ij}}{S _{ij}}} \right ), \end{equation}

in which the symmetric tensor $\mathcal{S}_{ij}$ (strain rate tensor) and anti-symmetric tensor $\Omega _{ij}$ (rotation rate tensor) are determined by

(3.2) \begin{align}& \mathcal{S}_{ij} = \frac {1}{2} \left ( \frac {\partial u_{f,i}}{ \partial x_j} + \frac {\partial u_{f,j}}{\partial x_i} \right )\!, \end{align}

(3.3) \begin{align}& {\Omega _{ij}} = \frac {1}{2} \left ( \frac {\partial u_{f,i}}{\partial x_j} - \frac {\partial u_{f,j}}{\partial x_i} \right )\!. \end{align}

The figure illustrates that the coherent vortical structure mostly originates on the upper TC interface and that it significantly coincides with the low-pressure centre (Hunt et al. Reference Hunt, Wray and Moin1988; Jeong & Hussain Reference Jeong and Hussain1995), both of which are likely to be related to the KH vortices close to the interface. Thus, the KH vortex’s dynamic features can be well described by the coherent vortical structure’s attributes, such as particle count, vortex size and energy. This serves as the foundation for the analysis of the KH vortex dynamics in the text that follows.

Figure 6.

Spatial distribution of particle vertical velocity divided by particle Stokes settling velocity $w_p / u_{\tau }$ at six selected times. Results are shown for ( $a$ ) $t = 0.6$ , ( $b$ ) $t = 1.8$ , ( $c$ ) $t = 2.4$ , ( $d$ ) $t = 3.0$ , ( $e$ ) $t = 3.6$ , ( $f$ ) $t = 4.2$ , ( $g$ ) $t = 5.5$ and ( $h$ ) $t = 6.7$ . The isosurface represents the vortical coherent structures captured at $Q=5$ with the colour indicating the spanwise fluid velocity $v_f$ .

The impact of the KH vortex on particle transport and settling behaviour is preliminarily addressed here. Figure 3 shows the impact of the KH vortex on horizontal particle transport, and the spatial distribution of the vertical particle velocity during the TC evolution is displayed in figure 6. As figure 3 illustrates, the interfacial KH vortices promote the streamwise transport of particles near the TC interface. In TCs with low Reynolds numbers (Xie et al. Reference Xie, Hu, Zhu, Yu and Pähtz2023a ), this vortex-related local particle transport enhancement is not visible. This can be explained in two ways: first, by the existence of the KH vortex, which prevents the ambient fluid and current from directly interacting, thereby lessening the ambient fluid’s shear on the TC interface; second, by the fluid force produced by the vortex structure acting as an accelerator.

The particle vertical velocity, shown in figure 6( $a$ ), is initially generally higher than the particle Stokes velocity, i.e. $|w_p|\gt u_\tau$ , due to the influence of upper backflow. Subsequently, the particles on the upstream side of the vortices continue to settle at a speed greater than the Stokes velocity, while the particles near the downstream side of the vortex – especially those close to the head – tend to move upward due to the anticlockwise interfacial vortices (figure 6 $c,\!d$ ). The reflected wave is the primary cause of the extremely rapid particle settling within $x = 0 - 0.5$ in figure 6( $b$ – $d$ ). In regions unaffected by the reflected wave and vortices, particles in the TC mostly settle at a speed of approximately $u_\tau$ .

The figure also presents a three-dimensional visualisation of the vortex coherent structures extracted by the $Q$ -criterion. In the early and middle stages, the vortex structures at the TC interface exhibit prominent two-dimensional characteristics. In the later stages (figure 6 $f$ ), the development of three-dimensional instability leads to the destruction of spanwise vortices and the breakdown of KH billows, consistent with the previous findings (Winant & Browand Reference Winant and Browand1974; Caulfield & Peltier Reference Caulfield and Peltier2000). Notably, in figure 6( $g$ , $h$ ) a series of long coherent structures appear, connecting the TC with the resident vortex in the intermediate region, with opposite spanwise fluid velocities $v_f$ .

We investigate the relationship between the relative settling velocity of particles and the $Q$ value at the particle position in the flow, as depicted in figure 7. The relative particle settling velocity, $w_p-w_f$ , is defined as the difference between the particle vertical velocity and the fluid vertical velocity, which here is used to eliminate the influence of the TC background flow. As seen in the figure, there is a general positive correlation between the relative settling velocity $w_p-w_f$ and the $Q$ value, denoted as $w_p-w_f \;{\sim}\; Q$ . Specifically, negative $Q$ values are associated with an increase in particle settling velocity, while positive $Q$ values correlate with a decrease in settling velocity, which aligns with previous cognition (Maxey Reference Maxey1987; Tom & Bragg Reference Tom and Bragg2019). Maxey (Reference Maxey1987) explained this relationship by the preferential sweeping of settling inertial particles through negative $Q$ vortex regions under centrifugal force, which causes downward-moving vortices to intensify their downward motion. Furthermore, the relative settling velocity is mostly close to $u_\tau$ , as illustrated by $ (w_p-w_f ) / u_\tau \approx -1$ in figure 7.

Figure 7.

Relationship between the vertical velocity of the particles after subtracting the fluid velocity ( $w_p-w_f$ ) and the $Q$ values at the particle positions for ${Re} = 3030$ . Results are shown for ( $a$ ) $t = 1.2$ , ( $b$ ) $t = 1.8$ , ( $c$ ) $t = 2.4$ , ( $d$ ) $t = 3.0$ , ( $e$ ) $t = 3.6$ and ( $f$ ) $t = 4.2$ . The black line shows the linear fitting result, and the fitting slope is shown in the figure.

Figure 8.

Time evolution of the current front $x_f$ .

Figure 9.

Spanwise-averaged TC profiles for different ${\textit{Re}}$ values at various times, which are distinguished by the particle volumetric fraction $\alpha _p = 10^{-5}$ . Results are shown for ( $a$ ) $t = 0.6$ , ( $b$ ) $t = 1.8$ , ( $c$ ) $t = 2.4$ , ( $d$ ) $t = 3.0$ , ( $e$ ) $t = 4.2$ and ( $f$ ) $t = 6.7$ .

3.2. Effect of Reynolds number on TC evolution

Figure 8 plots the time evolution of the front position, $x_f$ , of the TC under various ${\textit{Re}}$ conditions. The front position $x_f$ is defined as the distance between the farthest moving particle and the gate, i.e. $x_f = \max (x_{p,i}) - L_g$ , where $x_{p,i}$ is the particle position along the $x$ axis. As can be seen from the figure, the simulation results almost align with the previous results (Cantero et al. Reference Cantero, Lee, Balachandar and Garcia2007; Xie et al. Reference Xie, Zhu, Hu, Yu and Pan2023b ). However, the simulated TC front decelerates earlier, primarily due to the greater particle settling effect (Bonnecaze, Huppert & Lister Reference Bonnecaze, Huppert and Lister1993; Gladstone & Woods Reference Gladstone and Woods2000; Meiburg et al. Reference Meiburg, Radhakrishnan and Nasr-Azadani2015) and the aspect ratio of the initial fluid–particle column (Cantero et al. Reference Cantero, Lee, Balachandar and Garcia2007; Sun et al. Reference Sun, Zhang, Wang and Liu2020; Lu et al. Reference Lu, Liu, Xie, Sun, Yang and Guo2024). The time evolution of the front includes the initial collapse acceleration stage, the slumping stage characterised by a constant front velocity, and the later deceleration stage. The curves indicate that the constant velocity during the slumping stage is roughly 0.55, consistent with earlier simulation and experimental findings (Rottman & Simpson Reference Rottman and Simpson1983; Gladstone, Phillips & Sparks Reference Gladstone, Phillips and Sparks1998; Shin et al. Reference Shin, Dalziel and Linden2004; Steenhauer, Tokyay & Constantinescu Reference Steenhauer, Tokyay and Constantinescu2017). With increasing flow Reynolds number, the TC movement distance increases; additionally, the constant velocity of the slumping stage slightly increases, and its duration is extended. These findings are in line with the established understanding of how the Reynolds number influences TCs (Marino, Thomas & Linden Reference Marino, Thomas and Linden2005; Cantero et al. Reference Cantero, Lee, Balachandar and Garcia2007; Hu et al. Reference Hu, Tao, Li and He2020).

The TC profiles, captured by the particle volumetric fraction $\alpha _p = 10^{-5}$ , at six selected instants for various Reynolds numbers are displayed in figure 9. As can be observed in figure 9( $a$ ), the initial collapse event of the fluid–particle mixture is independent of the Reynolds number during the early stages of the TC. Later on (figure 9 $b$ – $f$ ), however, the front position becomes positively correlated with the Reynolds number, which is consistent with the observations regarding the front position in figure 8. At this time, TCs with higher Reynolds numbers typically exhibit greater current thickness, demonstrating a higher ability of particle suspension. Furthermore, more KH vortices can be excited at the interface at high Reynolds numbers, aligning with the findings from simulations of KH billow development in a horizontally periodic domain (Smyth Reference Smyth2003; Fontane & Joly Reference Fontane and Joly2008). This corresponds to a decrease in wavelength of KH instability. Given their ability to promote particle transport and obstruct settling (figures 3 and 6), these KH vortices near the TC interface may be the cause of the higher transport levels of high Reynolds number TCs.

5.1. Formation of KH billows

The instability of the TC upper interface is affected by velocity shear and density stratification. The stratification strength of the flow is characterised by the mixture density gradient in the numerator of (1.2), which can be carefully deduced to yield the following relationship:

(5.1) \begin{equation} \frac {\partial \rho _m}{\partial z} = \frac {\partial (\alpha _f \rho _f + \alpha _p \rho _p)}{\partial z} = \frac {\partial [(1-\alpha _p) \rho _f + \alpha _p \rho _p]}{\partial z} = (\rho _p - \rho _f)\frac {\partial \alpha _p}{\partial z}. \end{equation}

This implies that $\partial \rho _m / \partial z \varpropto \partial \alpha _p / \partial z$ , suggesting that the stratification can also be quantified by the particle volume fraction gradient $\partial \alpha _p / \partial z$ . Figure 12 displays the fluid streamwise velocity profile ( $u_{\!f}$ ) and particle concentration profile ( $\alpha _p/\alpha _{p0}$ ) of four cases at $t = 1.7$ , which is close to the critical moment of KH excitation. In terms of fluid velocity, the profiles for all cases are generally similar. The primary difference occurs in the region of strong shear near the upper interface of the TC flow – an increase in Reynolds number complicates the interface velocity profile at the upstream position (figure 12 $a$ , $c$ ). This may be a manifestation of greater susceptibility to instability. As for the particle concentration profiles, particles are suspended to a higher position at larger Reynolds numbers. Notably, the velocity profile does not exactly coincide with the particle concentration profile, especially at low Reynolds numbers ( ${Re}=1716$ ), as shown by the red line in figure 12( $g$ , $h$ ). The offset between the velocity and particle concentration interfaces is indicative of asymmetric KH instability, which is related to the hybrid mode of symmetric KH instability and Holmboe instability (Olsthoorn et al. Reference Olsthoorn, Kaminski and Robb2023; Yang et al. Reference Yang, Timmermans and Lawrence2024a ). Additionally, this offset may be related to the different dominant wavenumber and amplitude caused by different boundary conditions above and below the TC interface at different Reynolds numbers (Carpenter et al. Reference Carpenter, Tedford, Rahmani and Lawrence2010), which in turn lead to different interface instability modes. This instability ejects particles away from the TC upper interface, resulting in a sharp interface (Yang et al. Reference Yang, Timmermans and Lawrence2024a ), as demonstrated by the profile at ${Re}=1716$ in figure 9( $f$ ). Figure 13 shows the interfacial shear Reynolds number ${Re}_I$ of the velocity profile at different streamwise positions at $t=1.7$ . In a TC event, ${Re}_I$ is generally larger at the TC front, as it is the first to experience ambient fluid shear. As the TC flow Reynolds number ${\textit{Re}}$ increases, the interfacial shear Reynolds number ${Re}_I$ also increases, with a scaling relationship of ${Re}_I\;{\sim}\; {Re}^{1.25}$ typically holding between the two.

Figure 12.

Profiles of spanwise-averaged fluid streamwise velocity and particle volume fraction for varying ${\textit{Re}}$ cases at $t = 1.7$ . ( $a$ , $c$ , $e$ , $g$ ) Fluid streamwise velocity $u_{\!f}$ and ( $b$ , $d$ , $f$ , $h$ ) particle volume fraction that is divided by the initial particle volume fraction $\alpha _p / \alpha _{p0}$ . Results are shown for ( $a$ , $b$ ) $x = 0.44$ , ( $c$ , $d$ ) $x = 0.59$ , ( $e$ , $f$ ) $x = 0.73$ and ( $g$ , $h$ ) $x = 1.20$ .

Figure 13.

Interfacial shear Reynolds number ${Re}_I$ versus flow Reynolds number ${\textit{Re}}$ at various positions at $t = 1.7$ .

Figure 14.

Spatial distribution of gradient Richardson number $Ri_g$ of case ${Re} = 3836$ at ( $a$ ) $t = 1.7$ and ( $b$ ) $t = 3.7$ . Black lines depict the envelope of the TC given by the particle volume fraction ( $\alpha _p = 10^{-5}$ ). The high- $Ri_g$ region beneath the low- $Ri_g$ region, where the KH instability is excited, is outlined by a red dashed ellipse. The gradient Richardson number profiles at the selected vortex excitation positions are displayed in the inset in panel ( $a$ ) and the two right insets in panel ( $b$ ). The left inset in panel ( $b$ ) shows the statistical results of the gradient Richardson numbers $Ri_g$ near the TC interface where the interfacial vortices are excited at different ${\textit{Re}}$ values.

Figure 15.

Time evolution of ( $a$ ) interfacially averaged fluid streamwise velocity gradient $\langle -\textrm{d} u_{\!f} / \textrm{d} z\rangle _{\textit{IF}}$ , ( $b$ ) interfacially averaged particle volume fraction gradient $\langle -\textrm{d} \alpha _p / \textrm{d} z\rangle _{\textit{IF}}$ , and ( $c$ ) interfacially averaged gradient Richardson number $\langle Ri_g\rangle _{\textit{IF}}$ .

Figure 14 displays the spatial distribution of the gradient Richardson number $Ri_g$ , defined by (1.2), at two moments ( $t=1.7$ and $3.7$ ) for the ${Re}=3836$ case. The stratification effect is dominant in the green areas where $Ri_g \gt 1$ , whereas the shear effect is dominant in the pink areas where $0\lt Ri_g\lt 1$ . The TC front, the region where the TC body borders the bottom wall, and the area where the TC interacts with the ambient fluid are the three main locations where shear dominance occurs. As shown in the insets of figure 14, the gradient Richardson number ( $Ri_g$ ) near the interface, particularly at the location of the KH vortex initiation, is less than 1, with most values falling in the region of $Ri_g \lt 1/4$ . This observation aligns with previous studies (Turner Reference Turner1986; Kneller et al. Reference Kneller, Nasr-Azadani, Radhakrishnan and Meiburg2016). Notably, this threshold is less stringent than the critical value of $Ri_g = 1/4$ , a `rule of thumb’, which has been discussed by many studies (Harang et al. Reference Harang, Thual, Brancher and Bonometti2014; Howland, Taylor & Caulfield Reference Howland, Taylor and Caulfield2018; Parker, Caulfield & Kerswell Reference Parker, Caulfield and Kerswell2019). The offset between the particle concentration interface and the velocity interface induces asymmetry, which drives the upward extension of KH instability (Yang et al. Reference Yang, Timmermans and Lawrence2024a ) and contributes to the reduction of $Ri_g$ . Additionally, the unsteadiness of TC flow (Caulfield Reference Caulfield2021) and the reduction in velocity gradients due to particle settling dynamics (Kneller et al. Reference Kneller, Nasr-Azadani, Radhakrishnan and Meiburg2016) may enable KH instability to occur at slightly higher $Ri_g$ values. The critical threshold of $Ri_g$ in TCs thus remains a subject of ongoing investigation. Here $Ri_g$ values below the critical threshold do not always result in interface instability. In other words, a $Ri_g$ value below the threshold is only a prerequisite that must be met in order for the KH vortex to start.

It is worth mentioning that in the early stages of the lock-exchange TC, such as at $t=1.21$ in figure 3( $b$ ), slight signs of instability can be observed, which is primarily attributed to the vertical shear term ( $\partial w_f / \partial x$ ) at the interface of the collapsing fluid–particle column. During this phase, the Miles–Howard criterion, typically used to predict the stability of stratified shear flows, may no longer be valid. Subsequently, the particle column collapse progresses rapidly into the horizontal development stage (see $t = 1.82$ in figure 3 $d$ ). At this stage, most interface vortex structures are governed by the horizontal shear term ( $\partial u_{\!f}/\partial z$ ), with the exception of vortices near the sidewalls that remain strongly influenced by reflected bore. The empirical $Ri_g$ criterion can then be used to analyse these vortex structures.

The gradient Richardson number profile does not exhibit a discernible pattern of increase or decrease. Interestingly, though, a stratification-dominated region dyed with dark green ( $Ri_g \gt 10^2$ ) always lies directly below the $Ri_g\lt 1$ region where the interfacial vortices are excited. Because of the strong stratification effect, the cost of vertical momentum transfer and mixing between the upper and lower fluids increases to a certain extent (Noh & Fernando Reference Noh and Fernando1995; Caulfield Reference Caulfield2021). This creates a scenario where the lower fluids act as a barrier, preventing the upper fluids from penetrating downward, thereby, in a sense fostering an environment conducive to the upward excitation of instabilities.

To further investigate the relationship between interface instability and Reynolds number, we select a region extending 0.025 in both the upward and downward directions (i.e. a vertical thickness of 0.05 in dimensionless units) relative to the three-dimensional upper interface of the TC, and the physical quantities in this region are then averaged to obtain the interfacially averaged quantities $\langle \boldsymbol\cdot \rangle _{\textit{IF}}$ . The time evolution of the interfacially averaged fluid streamwise velocity gradient $\langle -\textrm{d} u_{\!f} / \textrm{d} z\rangle _{\textit{IF}}$ , particle volume fraction gradient $\langle -\textrm{d} \alpha _p / \textrm{d} z\rangle _{\textit{IF}}$ and gradient Richardson number $\langle Ri_g\rangle _{\textit{IF}}$ is plotted in figure 15. As can be observed in figure 15( $b$ ) and in figure 12, an increase in Reynolds number slightly increases the interface particle volume fraction gradient $\langle -\textrm{d} \alpha _p / \textrm{d} z\rangle _{\textit{IF}}$ , with the effect remaining within one order of magnitude. Due to the greater TC height at large ${\textit{Re}}$ , a larger fluid velocity gradient is observed in the interface (figure 12). The positive contribution of the increasing Reynolds number to $\langle -\textrm{d} u_{\!f} / \textrm{d} z\rangle _{\textit{IF}}$ is more than an order of magnitude as illustrated in figure 15( $a$ ), exhibiting greater than its effect on $\langle -\textrm{d} \alpha _p / \textrm{d} z\rangle _{\textit{IF}}$ . Accordingly, this produces a smaller gradient Richardson number at a larger Reynolds number, as reported in previous studies (Salehipour et al. Reference Salehipour, Peltier, Whalen and MacKinnon2016b ; Caulfield Reference Caulfield2021), which means that the TC interface is more susceptible to shear instability. This explains why higher Reynolds number TCs develop larger and more numerous KH billows, as shown in figure 9.

Figure 16.

Spatial distribution of forces acting on the particles of case ${Re} = 3030$ at $t = 1.7$ . ( $a$ , $b$ ) Total forces $\boldsymbol{F}^T$ , ( $c$ , $d$ ) effective drag force $\boldsymbol{F}^{Ed}$ , ( $e$ , $f$ ) lift force $\boldsymbol{F}^l$ , ( $g$ , $h$ ) added mass force $\boldsymbol{F}^{\textit{add}}$ , ( $i$ , $j$ ) pressure gradient force $\boldsymbol{F}^{q}$ and ( $k$ , $l$ ) viscous force $\boldsymbol{F}^{\boldsymbol\nabla \boldsymbol\cdot \boldsymbol{\tau }}$ . ( $a$ , $c$ , $e$ , $g$ , $i$ , $k$ ) Streamwise component $\boldsymbol{F}_x$ and ( $b$ , $d$ , $f$ , $h$ , $j$ , $l$ ) vertical component $\boldsymbol{F}_z$ .

Figure 17.

Spatial distribution of forces acting on the interfacial particles of various cases at $t = 1.7$ . ( $a$ , $b$ ) Total forces $\boldsymbol{F}^T$ , ( $c$ , $d$ ) effective drag force $\boldsymbol{F}^{Ed}$ , ( $e$ , $f$ ) lift force $\boldsymbol{F}^l$ , ( $g$ , $h$ ) added mass force $\boldsymbol{F}^{\textit{add}}$ and ( $i$ , $j$ ) pressure gradient force $\boldsymbol{F}^{q}$ . ( $a$ , $c$ , $e$ , $g$ , $i$ ) Streamwise component $\boldsymbol{F}_x$ and ( $b$ , $d$ , $f$ , $h$ , $j$ ) vertical component $\boldsymbol{F}_z$ . The scatter represents the forces on the interfacial particles (within 0.05 dimensionless thickness from the TC interface), and the curve represents the spatial variation of average force on the interfacial particles $\langle \boldsymbol{F} \rangle _{\textit{IF}}$ .

Figure 16 depicts the spatial distribution of the forces acting on the particles in the vicinity of the KH excitation moment ( $t=1.7$ ) at ${Re} = 3030$ . These forces include the total force on the particles $\boldsymbol{F}^T$ , the effective drag force $\boldsymbol{F}^{Ed}$ (the sum of the drag force and the effective gravity), the lift force $\boldsymbol{F}^l$ , the added mass force $\boldsymbol{F}^{\textit{add}}$ , the pressure gradient force $\boldsymbol{F}^{q}$ and the viscous force $\boldsymbol{F}^{\boldsymbol\nabla \boldsymbol\cdot \boldsymbol{\tau }}$ . Given the low particle volume fraction of the TCs in this study ( $O(10^{-4} \;{\sim}\; 10^{-3})$ ) and accordingly the low particle collision frequency, the contact force is negligible on average. Thus, the contact force is not shown in the figure.

The components of each force in figure 16 are all smaller than the predominant vertical effective gravity of the suspended particles driving the TC process, which is associated with the dilution regime associated at low volume fraction (Balachandar & Eaton Reference Balachandar and Eaton2010). Nevertheless, despite its small magnitude, this interaction force plays a subtle yet significant role in the evolution of TCs through complex interphase interactions (Xie et al. Reference Xie, Hu, Pähtz, He and Cheng2022). With regard to the total force acting on the particles of KH vortices, it is observed that in the streamwise direction (figure 16 $a$ ) the particles in the front part of the vortex experience a negative total force, whereas the particles in the rear undergo a positive total force; in the vertical direction (figure 16 $b$ ) the particles in the upper and lower regions experience, respectively, a negative total force and a positive total force. This suggests that the fluid force applied by the KH vortex on the particles is convergent, which is in line with Lee (Reference Lee2019). This effect facilitates the particle accumulation within the vortex. Conversely, the particles’ response to the fluid is divergent, meaning that they will, in a way, impede the formation of the KH structure. The forces on particles in the TC are also affected by the KH vortices. The majority of particles, with the exception of those in the head region, receive positive $\boldsymbol{F}^T_x$ and $\boldsymbol{F}^T_z$ , which facilitate streamwise transport and vertical suspension of particles, respectively. This explains the spatial distribution of particle velocity shown in figures 3 and 6.

The spatial distribution of the total force on particles near the KH billows exhibits significant similarity to the spatial distribution of the added mass force (figure 16 $g$ , $h$ ) and the pressure gradient force (figure 16 $i$ , $j$ ). Given that the fluid acceleration and particle acceleration are relatively synchronised when $St \ll 1$ , the strong correlation observed between the total force and the added mass force makes sense. The pressure gradient force is primarily closely related to the spatial distribution of pressure in figure 5. It is evident that the low-pressure centre formed by the interfacial vortex plays a crucial role in entraining particles.

Furthermore, as figure 16( $c$ ) illustrates, shear from the ambient fluids causes particles close to the upper TC interface, including those in KH billows, to experience a negative streamwise effective drag force. There is a negative vertical effective drag on average (figure 16 $d$ ) because the vertical drag is insufficient to resist the particle gravity as they settle. Positive streamwise lift force and positive vertical lift force are produced on the interfacial particles by the negative $u_{\!f}-u_{{\kern-1pt}p}$ and positive $w_f-w_p$ under the control of the anticlockwise flow near the interface, which is represented as positive vorticity $\Omega _{xz}$ . It validates the lift’s beneficial contribution to particle suspension and transport in TCs, as concluded by Xie et al. (Reference Xie, Hu, Pähtz, He and Cheng2022), especially to rolling up the interfacial particles. On average, the fluid’s viscous force on the particles inhibits their streamwise movement at the TC interface but helps to their suspension.

Akin to the distribution of the added mass force and the pressure gradient force in figure 16( $g$ – $j$ ), the fluid force on the particles away from the TC interface is also evidently influenced by the KH billows. Due to the large-scale start-up vortex created by rolling up the edges of the vorticity sheet, the particles in the head are subject to a negative total streamwise force. It is explained by the dominance of the added mass force and the pressure gradient force. The back particles experience a positive streamwise force from the fluids.

Figure 17 shows the various forces acting on the TC interface particles at different flow Reynolds numbers, and plots the distribution of the average forces applied to the particles along the $x$ direction. The viscous force $\boldsymbol{F}^{\boldsymbol\nabla \boldsymbol\cdot \boldsymbol{\tau }}$ is not depicted in the figure since it is of a comparatively negligible order as compared with other predominant forces (figure 16). As can be observed, with the Reynolds number of the TC increasing, the total force and the component forces also increase, though to different degrees. In particular, the lift force, added mass force and pressure gradient force – crucial for particles rolling up – increase with Reynolds number, facilitating the rolling up of interface particles in instability billows. This explains why, as the Reynolds number rises, the behaviour of particles rolling up becomes more pronounced as shown in figure 9.

5.2. Evolution mechanism of KH billows

Figure 18.

( $a$ ) Relative particle volumetric fraction $\langle \alpha _{p} \rangle _V / \alpha _{p0}$ and ( $b$ ) average particle Reynolds number $\langle {Re}_p \rangle _V$ in the typical KH vortices of case ${Re} = 3030$ . Vortices A–D are shown in figure 5.

Figure 19.

Spatial distribution of finite-time Lyapunov exponent (FTLE) when ${Re} = 3030$ at six instants. Results are shown for ( $a$ ) $t = 0.61$ , ( $b$ ) $t = 1.21$ , ( $c$ ) $t = 1.52$ , ( $d$ ) $t = 1.82$ , ( $e$ ) $t = 2.42$ and ( $f$ ) $t = 3.03$ . Vortical coherent structures are captured by $Q$ -criterion ( $Q = 5$ ) that are marked in cyan.

For the KH billow evolution, we characterise the particle-carrying capacity and the strength of the fluid–particle interaction by using two physical quantities. They are the average particle concentration $\langle \alpha _{p0}\rangle _V/\alpha _{p0}$ and average particle Reynolds number $\langle {Re}_p \rangle _V$ of the vortex structures identified by $Q = 5$ . The particle Reynolds number ${Re}_p$ is determined by (2.15). The average particle volume fraction and average particle Reynolds number of four typical vortex structures in the case of ${Re}=3030$ (vortices A–D, marked in figure 5) are displayed over time in figure 18. Figure 18( $a$ ) shows $\langle \alpha _{p0} \rangle _V/\alpha _{p0} \lt 1$ , suggesting that the average particle volume fraction of the vortices is always less than the initial particle concentration. The average particle volume fraction continues to drop as the interfacial vortex develops due to the expansion of the vortex structure (figure 5) and the rapid decline in the number of particles within it. The reflected bore is the cause of the increase in vortex A’s particle volume fraction at $t=2$ . With regard to the particle Reynolds number, figure 18( $b$ ) reveals a progressive increase until the vortex merges or dissipates. It implies the enhancement of fluid–particle interactions as the vortices evolve.

The finite-time Lyapunov exponent (FTLE), which quantifies the local fluid stretching level, is often used to extract Lagrangian coherent structures (Lipinski & Mohseni Reference Lipinski and Mohseni2010; Balasuriya, Kalampattel & Ouellette Reference Balasuriya, Kalampattel and Ouellette2016; Rist, Weinschenk & Wenzel Reference Rist, Weinschenk and Wenzel2023). It helps in the assessment of the flow instability. Here, we compute the FTLE using the forward method in order to understand the evolution of interfacial KH billows in TCs. The spatial distribution of the FTLE at six chosen times during the TC process with ${Re}=3030$ is displayed in figure 19. Darker region in the image correspond to larger FTLE values where the fluid stretching is stronger.

A prominent feature of the large-scale FTLE structure with high stretching is its close correlation with the TC interface and the vortical coherent structures in the flow. The KH billows, which are part of the interface morphology, also exhibit a correlation with the FTLE structure. There is only a thin, sizable FTLE region close to the interface in the very early stage shown in figure 19( $a$ ). The region progressively widens and the FTLE value dramatically rises as the current develops (figure 19 $b$ ), indicating an improvement in fluid stretching level. At $t = 1.52$ (figure 19 $c$ ), the large FTLE streak begins to become distorted, coinciding with the onset of interfacial particle flow instability (see figure 3 $c$ ). As the KH vortex is excited, the large FTLE streak evolves into a spatially continuous braided structure at $t=1.82$ in figure 19( $d$ ). Three large FTLE streaks are visible in the typical vortex structure at $(x,z) = (0.6,0.6)$ in figure 19( $d$ ). Two of these streaks originate from the downstream and upstream vortex structures, respectively, and they wrap around the middle largeFTLE streak up and down, which is located at the vortex centre. The particles are rolled up and firmly attached to the large FTLE streak from the downstream. The TC is connected to the ambient fluid by the low FTLE fluid in this braided structure, which destroys the initial spatially continuous large FTLE streak and facilitates the two fluids to mix. As the mixing progresses, the fluid near the TC interface exhibits a decreasing stretching level in figure 19( $e$ , $f$ ), and some vortex structures are engulfed by other vortices.

The pairing and merging of KH vortices play a critical role in the development of interfacial instability, which significantly influences the mixing of the two fluids (TC and ambient fluid). This phenomenon represents the dominant two-dimensional instability in the evolution of KH billows (Yang et al. Reference Yang, Timmermans and Lawrence2024a ). In this context, the phase difference between the primary KH mode and the subharmonic mode, which governs vortex pairing (Dong et al. Reference Dong, Tedford, Rahmani and Lawrence2019; Liu et al. Reference Liu, Kaminski and Smyth2022), is influenced by the lock-exchange TC flow, particularly by the reflected bore that propagates downstream. We focus on the pairing process shown in figure 19( $e$ , $f$ ), in which the three interfacial vortices pair and merge into two. This may correspond to the case of the phase of the subharmonic mode relative to the KH mode $\theta _{sub} = -\pi /2$ proposed in figure 1( $b$ ) of Dong et al. (Reference Dong, Tedford, Rahmani and Lawrence2019), where the middle vortex is weakened and the vortices on both sides are strengthened.

Figure 20.

Merging process of KH vortices during the TC evolution. Spatial distribution of ( $a$ – $e$ ) fluid rotation rate tensor $\Omega _{xz}$ , ( $f$ – $j$ ) fluid strain rate tensor $\mathcal{S}_{xz}$ , ( $k$ – $o$ ) FTLE, ( $p$ – $t$ ) particle streamwise velocity $u_{{\kern-1pt}p}$ and ( $u$ – $y$ ) total streamwise force acting on particles $\boldsymbol{F}^T_x$ . The first to fifth columns represent $t=2.42, 2.55, 2.67, 2.79$ and $3.03$ , respectively.

Here, we investigate this evolution mechanism from both the fluid and the particle phases. The spatial distribution of fluid rotation rate tensor $\Omega _{xz}$ , fluid strain rate tensor $\mathcal{S}_{xz}$ , FTLE, particle streamwise velocity $u_{{\kern-1pt}p}$ and total streamwise force acting on particles $\boldsymbol{F}^T_x$ during the vortex merging process of case ${Re}=3030$ are displayed in figure 20. Of note is the wall-reflected bore produced by the lock-exchange set-up. As shown in the figure, part of the flow structure formed by the development of the reflected bore gradually merges into vortex A, while the remainder travels beneath the KH billow. Although this process slightly affects the engulfment of vortex B, it does not qualitatively disrupt the merging mechanism. Therefore, the understanding of vortex merging presented here remains largely applicable to other similar scenarios.

Strong positive vorticity $\Omega _{xz}$ dominates the core of the KH vortices from the perspective of vorticity distribution (Parker et al. Reference Parker, Caulfield and Kerswell2019). Continuous vortex filaments develop on the front and rear sides of vortex B, connecting it with adjacent vortices (vortices C, D and A) (Parker et al. Reference Parker, Caulfield and Kerswell2019) before vortex B is engulfed and disappears. The fluid strain rate $\mathcal{S}_{xz}$ in the vortex filament area is very high as shown in figure 20( $f$ ), particularly at the location where the vortex filaments are incredibly thin (Corcos & Sherman Reference Corcos and Sherman1976; Dong et al. Reference Dong, Tedford, Rahmani and Lawrence2019). Vortex B, as well as the front and rear vortices, are traversed by the braided structure of the large FTLE (figure 20 $k$ , $l$ ). The combined regime of strong vorticity and strong strain rate has a significant impact on the distribution of the FTLE. Vortex B flattens out as it shrinks due to the stretching force from the front and rear vortices. At the same time, the high fluid strain rate areas corresponding to the front and rear vortex filaments gradually merge in figure 20( $f$ – $i$ ), and the braided structure of the large FTLE gradually becomes blurred. Finally, vortex B is eliminated. The regions of high FTLE, high strain rate $\mathcal{S}_{xz}$ and high positive vorticity $\Omega _{xz}$ gradually overlap during this merging process.

From the perspective of the particle movement, the root of vortex A moves forward more quickly due to the greater positive fluid streamwise force at the roots (figure 20 $u$ , $v$ ). As a result, vortex A in figure 20( $r$ , $w$ ) catches up with the particles entrained in vortex B. When $t=2.55{\;{\sim}\;} 2.79$ , vortex A begins to progressively compress vortex B’s spatial position, leading to the eventual merging of vortex B into vortex A’s root. Regarding the different component forces illustrated in figure 31, the drag force, the added mass force and the pressure gradient force in the streamwise direction all work together to promote the merging process of this particle flow. This is largely affected by the strong pressure gradient between vortices, as shown in figure 5. Additionally, from figure 20( $f$ – $j$ , $p$ – $t$ ), it can be observed that particles mainly accumulate or have a tendency to accumulate locally in regions with high strain rate $\mathcal{S}_{xz}$ rather than high vorticity $\Omega _{xz}$ , which accords with the general understandings of particle flows (Wang & Maxey Reference Wang and Maxey1993; Cantero et al. Reference Cantero, Balachandar and García2008).

7.1. Energy budget of TC evolution

We examine the influence of interfacial vortex structures on the evolution of TCs by analysing changes in energy components. First, we assess the total energy budget. The energy of the lock-exchange TC originates from the particle gravitational potential energy $E_p^p$ that is generated by particle settling (Parker, Fukushima & Pantin Reference Parker, Fukushima and Pantin1986; Necker et al. Reference Necker, Härtel, Kleiser and Meiburg2002; Nasr-Azadani & Meiburg Reference Nasr-Azadani and Meiburg2014; Zgheib, Bonometti & Balachandar Reference Zgheib, Bonometti and Balachandar2015). As particles settle, the fluid potential energy $E_p^f$ increases. The available potential energy, $E_{p,\textit{avail}}$ , is defined as the difference between the decrease in particle gravitational potential energy and the increase in fluid potential energy, i.e. $E_{p,\textit{avail}} = -\Delta E^p_p - \Delta E^f_p$ , in which the symbol $\Delta$ represents the variation in energy with respect to the initial time. Within the TC system, this available potential energy is converted into particle kinetic energy $E_k^p$ , fluid kinetic energy $E_k^f$ and dissipated energy $E_{\textit{Diss}}$ . This results in the following relationship (Xie et al. Reference Xie, Zhu, Hu, Yu and Pan2023b ):

(7.1) \begin{equation} \Delta {E_{\textit{Diss}}} = - \big[ {\Delta E_p^p + \Delta E_p^f + \Delta E_k^p + \Delta E_k^f} \big]. \end{equation}

Particle gravitational potential energy $E_p^p$ , particle kinetic energy $E_k^p$ , fluid potential energy $E_p^f$ and fluid kinetic energy $E_k^f$ are defined as follows in accordance with earlier research (Xie et al. Reference Xie, Hu, Pähtz, He and Cheng2022, Reference Xie, Zhu, Hu, Yu and Pan2023b ; Lu et al. Reference Lu, Liu, Xie, Sun, Yang and Guo2024; Yang et al. Reference Yang, Chen and Wang2024b ):

(7.2) \begin{align} E_p^p\left (t\right ) &= \sum \limits _{i=1}^{N_p} m_i \left |\boldsymbol{g}\right | z_{p,i}, \end{align}

(7.3) \begin{align} E_k^p\left (t\right ) &= \sum \limits _{i=1}^{N_p} \left ( \frac {1}{2} m_i |\boldsymbol{u}_{p,i}|^2 + \frac {1}{2} I |\boldsymbol{\omega }_{p,i}|^2 \right )\!, \end{align}

(7.4) \begin{align} E_p^f\left (t\right ) &= \int _{\mathcal{D}} \alpha _f \rho _f \left |\boldsymbol{g}\right | z \textrm{d} V, \end{align}

(7.5) \begin{align} E_k^f\left (t\right ) &= \int _{\mathcal{D}} \frac {1}{2} \alpha _f \rho _f |\boldsymbol{u}_f|^2 \textrm{d} V. \end{align}

Here $N_p$ denotes the particle number, $z_{p,i}$ is the elevation of the particle $i$ and $\mathcal{D}$ represents the whole simulation domain. Figure 26 illustrates how the various energy components alter as the TC evolves. With the initial particle gravitational potential energy $E^p_p(0)$ , the energy components in the figure are dimensionless.

Figure 26.

Time evolution of the variation in energy components for varying cases. ( $a$ ) Particle gravitational potential energy $\Delta E^p_p$ , ( $b$ ) fluid potential energy $\Delta E^f_p$ , ( $c$ ) available potential energy $\Delta E_{p,\textit{avail}}$ , ( $d$ ) particle kinetic energy $\Delta E^p_k$ , ( $e$ ) fluid kinetic energy $\Delta E^f_k$ and ( $f$ ) dissipated energy $\Delta E_{\textit{Diss}}$ .

The figure illustrates that, as the currents evolve, the fluid potential energy $E_p^f$ , available potential energy $E_{p,\textit{avail}}$ and dissipated energy $E_{\textit{Diss}}$ all continue to rise, while the particle gravity potential energy $E_p^p$ consistently decreases. Particle gravity potential energy, fluid potential energy and available potential energy all progress (increase or decrease) at much higher rates at $t=0{\;{\sim}\;} 2$ than they do in the middle and later stages, as previously reported by other researchers (Necker et al. Reference Necker, Härtel, Kleiser and Meiburg2002, Reference Necker, Härtel, Kleiser and Meiburg2005; Ooi, Constantinescu & Weber Reference Ooi, Constantinescu and Weber2007; Nasr-Azadani & Meiburg Reference Nasr-Azadani and Meiburg2014). This rapid change is attributed to the initial collapse of the fluid–particle mixture. Both the fluid and particle kinetic energies ( $E_k^f$ and $E_k^p$ ) increase prior to $t=2$ and decrease thereafter. This pattern is due to the substantial supply of available potential energy during the early stage, which drives the kinetic energies, and the subsequent dissipation that reduces them in the later stage (Espath et al. Reference Espath, Pinto, Laizet and Silvestrini2015). The crucial time node, $t=2$ , marks when approximately $70\,\% -80\,\%$ of the initial particle potential energy $E^p_p(0)$ is transformed, yielding $E_{p,\textit{avail}}(t=2) = (\rho _p - \rho _f )/\rho _p |\Delta E^p_p (t=2) | \approx 0.3$ . Around $60\,\% -70\,\%$ of this available energy is crucial for the kinetic energies of particles and fluids, while the remaining $30\,\% -40\,\%$ of the energy is dissipated due to viscosity. This dissipation is slightly greater than that reported in previous studies (Espath et al. Reference Espath, Pinto, Laizet and Silvestrini2014; Nasr-Azadani & Meiburg Reference Nasr-Azadani and Meiburg2014), likely attributed to the increased terminal velocity of the particles, $u_\tau$ , in our TC simulations (Espath et al. Reference Espath, Pinto, Laizet and Silvestrini2015).

Both the potential energy change ( $\Delta E_p^p$ in figure 26 $a$ , $\Delta E_p^f$ in figure 26 $b$ and $\Delta E_{p,\textit{avail}}$ in figure 26 $c$ ) and dissipated energy (figure 26 $f$ ) decrease with increasing Reynolds numbers ${\textit{Re}}$ . However, the kinetic energies of the particles and the fluid both rise (Espath et al. Reference Espath, Pinto, Laizet and Silvestrini2014), which is consistent with the observation that the TC front moves more quickly at higher Reynolds number in figure 8 (Cantero et al. Reference Cantero, Lee, Balachandar and Garcia2007). This may be related to the enhanced large-sized vortex structures and KH billows at higher Reynolds numbers (Necker et al. Reference Necker, Härtel, Kleiser and Meiburg2002; Espath et al. Reference Espath, Pinto, Laizet and Silvestrini2014; Francisco et al. Reference Francisco, Espath and Silvestrini2017). An intriguing implication is that a rise in Reynolds number will result in a decrease in the conversion ratio from the system’s available potential energy to dissipated energy, while increasing the conversion ratio to the kinetic energies. At the end of the simulation ( $t=10$ ), the potential energies ( $E_p^p$ , $E_p^f$ and $E_{p,\textit{avail}}$ ) are quickly converted and stabilise, and the particle kinetic energy $E_k^p$ approaches zero. However, the fluid kinetic energy $E_k^f$ is still in the process of being converted into dissipated energy and has not yet fully dissipated. This indicates the need for a longer simulation time to reach the final equilibrium state.

7.2. Energy dissipation route analysis

It is necessary to investigate how the dissipation process is weakened by an increase in Reynolds number and whether this effect is attained by varying the interfacial vortex structures. The dissipated energy is separated into three categories: dissipation from fluid viscosity $E_{\textit{Diss}}^f$ (called fluid dissipation), dissipation from particle collisions $E_{\textit{Diss}}^{p \leftrightarrow p}$ (called particle dissipation) and dissipation from fluid–particle interactions $E^{f \leftrightarrow p}_{\textit{Diss}}$ (called interaction dissipation). Thus, the following relationship is given:

(7.6) \begin{equation} E_{\textit{Diss}} = E_{\textit{Diss}}^f + E_{\textit{Diss}}^{p \leftrightarrow p} + E^{f \leftrightarrow p}_{\textit{Diss}}. \end{equation}

The dissipation resulting from particle collisions can be disregarded in this work, i.e. $E_{\textit{Diss}}^{p \leftrightarrow p} \approx 0$ , because the particle volume fraction $\alpha _{p0}$ of the TC is very low ( $O(10^{-4} \;{\sim}\; 10^{-3})$ ). Thus, (7.6) can be rewritten as

(7.7) \begin{equation} E_{\textit{Diss}} = E_{\textit{Diss}}^f + E^{f \leftrightarrow p}_{\textit{Diss}}. \end{equation}

The energy dissipated at time $t$ due to fluid viscosity, $E_{\textit{Diss}}^f$ , can be computed as

(7.8) \begin{equation} E_{\textit{Diss}}^f (t) = \int _0 ^t \int _{\mathcal{D}} \alpha _f \rho _f \epsilon \textrm{d} V {\textrm{d}}t, \end{equation}

where $\epsilon$ denotes the fluid dissipation rate, which is determined by (e.g. Rahmani et al. Reference Rahmani, Lawrence and Seymour2014a ; Vassilicos Reference Vassilicos2015; He et al. Reference He, Zhao, Zhu and Hu2019)

(7.9) \begin{equation} \epsilon = \nu _f \langle \mathcal{S}_{ij}\mathcal{S}_{ij} \rangle _{\mathcal{F\,V}}, \end{equation}

where $\langle \boldsymbol\cdot \rangle _{\mathcal{F\,V}}$ denotes an average over the fluid volume. Figure 27 shows the time evolution of the dissipated energy ( $\Delta E_{\textit{Diss}}^f$ , $\Delta E_{\textit{Diss}}^V$ , $\Delta E_{\textit{Diss}}^{f \leftrightarrow p}$ ), energy dissipation rate ( $\textrm{d} E_{\textit{Diss}}^f/\textrm{d}t$ , $\textrm{d} E_{\textit{Diss}}^V/\textrm{d}t$ , $\textrm{d} E_{\textit{Diss}}^{f \leftrightarrow p}/\textrm{d}t$ ) and the proportion of the dissipation component in the total dissipated energy ( $\Delta E_{\textit{Diss}}^f / \Delta E_{\textit{Diss}}$ , $\Delta E_{\textit{Diss}}^V / \Delta E_{\textit{Diss}}$ , $\Delta E_{\textit{Diss}}^{f \leftrightarrow p} / \Delta E_{\textit{Diss}}$ ) of different components, where $\Delta E_{\textit{Diss}}^V$ is the fluid dissipation of the vortical coherent region. The vortical coherent structures are captured by the $Q$ -criterion ( $Q=5$ ) as well, which essentially characterises the KH billows as seen in figure 5. The curves indicate that the dissipation process can be divided into two stages based on the dissipation rate, $t\lt 2$ and $t\gt 2$ , which shows the impact of collapse during the initial stage. The change rate of each dissipation component is rapidly increasing in the earlier phase, as illustrated in figure 27( $d$ – $f$ ), and is generally decreasing in the subsequent stage. As a result, the early stage of dissipated energy in figure 27( $a$ – $c$ ) exhibits a faster growth rate than in the later stage.

Figure 27.

Process analysis of the energy dissipation. Time evolution of variation in ( $a$ ) fluid dissipation for the whole fluid domain $\Delta E^f_{\textit{Diss}}$ , ( $b$ ) fluid dissipation for the vortex domain $\Delta E^V_{\textit{Diss}}$ and ( $c$ ) interaction dissipation $\Delta E^{f \leftrightarrow p}_{\textit{Diss}}$ . Time evolution of the change rate of the above three dissipation components: ( $d$ ) $\textrm{d} E^f_{\textit{Diss}}/\textrm{d}t$ , ( $e$ ) $\textrm{d} E^V_{\textit{Diss}}/\textrm{d}t$ and ( $f$ ) $\textrm{d} E^{f \leftrightarrow p}_{\textit{Diss}}/\textrm{d}t$ . Time evolution of variation in the above dissipation components normalised by the total dissipation variation: ( $g$ ) $\Delta E^f_{\textit{Diss}} /\Delta E_{\textit{Diss}}$ , ( $h$ ) $\Delta E^V_{\textit{Diss}} /\Delta E_{\textit{Diss}}$ and ( $i$ ) $\Delta E^{f \leftrightarrow p}_{\textit{Diss}} /\Delta E_{\textit{Diss}}$ .

Generally speaking, an increase in Reynolds number at $t=0{\;{\sim}\;} 1$ reduces interaction dissipation (figure 27 $f$ ) but has a limited impact on fluid dissipation (figure 27 $d$ , $e$ ). The increase in ${\textit{Re}}$ at $t=1-4$ results in a decrease in the energy dissipation of each component. This has to do with the KH billows that form at the interface after $t=1$ . The increase in fluid dissipation at lower Reynolds numbers also may be attributed to the significant asymmetry between the fluid velocity and density profiles (Yang et al. Reference Yang, Timmermans and Lawrence2024a ), as shown in figure 12. The larger the Reynolds number, the stronger the interfacial vortex, which, acting as a transition zone, plays a role in reducing the shear between the ambient fluid and the TC. The curves for various ${\textit{Re}}$ cases in figure 27( $a$ – $c$ ) have similar increments because, after $t = 4$ , the dissipation energy change rates become essentially independent of the Reynolds number (with the exception of ${Re} = 1716$ in figure 27 $f$ ). Because of the ongoing advancement of the reflected bore, the critical moment $t=4$ , where the Reynolds number shows a distinct effect, occurs slightly later than the moment $t=2$ , when the dissipation change rate reaches its maximum. After $t = 4$ , the interaction dissipation at low Reynolds number ( ${Re} = 1716$ ) drastically drops and almost reaches zero. The reason for this is that the current has essentially stopped moving at this time.

The fluid-particle interactions are responsible for a substantial portion of the energy dissipation in the TC process, as shown in figure 27( $g$ – $i$ ). While the percentage of dissipation resulting from fluid–particle interaction in figure 27( $i$ ) steadily declines, the percentage of dissipation due to the fluid viscosity in figure 27( $g$ ) keeps rising. When $t\lt 4$ , the proportion of the dissipation brought on by fluid viscosity in the vortex structure (figure 27 $h$ ) initially rises and subsequently falls. Interestingly, after $t\gt 4$ , this proportion $\Delta E_{\textit{Diss}}^V / \Delta E_{\textit{Diss}}$ essentially stays constant and is independent of the Reynolds number, which is about 0.01.

Figure 28.

( $a$ ) Time evolution of average dissipation rates in the vortex region and in the fluid simulation domain $\langle \epsilon \rangle _V$ and $\langle \epsilon \rangle _f$ , and ( $b$ ) time evolution of the ratio between them $\langle \epsilon \rangle _V/ \langle \epsilon \rangle _f$ .

As illustrated in figure 28, we then go into great detail about the relationship between the average dissipation rate of the fluid region, $\langle \epsilon \rangle _f$ , and that of the vortex region, $\langle \epsilon \rangle _V$ . Both $\langle \epsilon \rangle _f$ and $\langle \epsilon \rangle _V$ rapidly increase in the early stages before progressively declining, essentially in time with the change in fluid kinetic energy in figure 26( $e$ ). The peak of $\langle \epsilon \rangle _V$ occurs earlier than that of $\langle \epsilon \rangle _f$ . The value of $\langle \epsilon \rangle _V/\langle \epsilon \rangle _f$ in figure 28( $b$ ), which always stays greater than 1, corresponds to the average dissipation rate in the vortex being always higher than the globally averaged fluid dissipation rate. Notably, its ratio exceeds 10 during the initial collapse stage of $t\lt 2$ , highlighting a significantly higher dissipation level for the fluids within the vortex structure. The ratio keeps going down as the TC evolves, mostly since $\langle \epsilon \rangle _V$ decays more quickly than $\langle \epsilon \rangle _f$ . The Reynolds number has a negligible effect on the ratio at $t\lt 2$ ; however, at $t\gt 2$ , the ratio increases with increasing Reynolds number, which is indicative of a higher dissipation rate of vortex fluids. This observation further supports the notion that vortex structures are stronger at higher Reynolds numbers, consistent with the interpretation of the preceding figure 5.

Figure 29.

Relationship between the average dissipation rate components of the interfacial vortices as the TC evolves: ( $a$ ) $\langle \epsilon _{\textit{NS}}\rangle _V$ vs $\langle \epsilon _{\textit{SD}}\rangle _V$ and ( $b$ ) $\langle \epsilon _{xx}\rangle _V$ vs $\langle \epsilon _{zz}\rangle _V$ . The colours of the symbols in panels ( $a$ , $b$ ) represent time $t$ . Snapshots of different stages of case ${Re} = 3030$ : ( $c$ ) $t = 1.21$ (stage I), ( $d$ ) $t = 1.33$ (the end of stage I), ( $e$ ) $t = 1.82$ (stage II), ( $f$ ) $t = 2.18$ (the end of stage II) and ( $g$ ) $t = 2.78$ (stage III). In panels ( $c$ – $g$ ) the slices show the spatial distribution of fluid vorticity $\Omega _{xz}$ , the black contour represents the TC profiles determined by the particle volume fraction $\alpha _p=10^{-5}$ , and the cyan contours are vortex structures captured by $Q$ -criterion ( $Q=5$ ).

Here we discuss the different components of average dissipation rate in the vortex structure. We decompose the dissipation rate $\epsilon$ into a part related to the normal strain $\epsilon _{\textit{NS}}$ and a remaining part related to the shear deformation $\epsilon _{\textit{SD}}$ :

(7.10) \begin{equation} \epsilon = \epsilon _{\textit{NS}} + \epsilon _{\textit{SD}}. \end{equation}

Given that the TCs in our work mainly exhibit noticeable two-dimensional characteristics before the onset of three-dimensional instability, the spanwise dissipation component at this stage is relatively negligible. Thus, $\epsilon _{\textit{NS}}$ and $\epsilon _{\textit{SD}}$ can be expressed by

(7.11) \begin{align} \epsilon _{\textit{NS}} &= \underbrace {\nu _f \langle \mathcal{S}_{xx}\mathcal{S}_{xx}\rangle _{\mathcal{F\,V}}}_{\epsilon _{xx}} + \underbrace { \nu _f \langle \mathcal{S}_{zz}\mathcal{S}_{zz}\rangle _{\mathcal{F\,V}}}_{\epsilon _{zz}}, \end{align}

(7.12) \begin{align} \epsilon _{\textit{SD}} &= \nu _f \left ( \langle \mathcal{S}_{xz}\mathcal{S}_{xz}\rangle _{\mathcal{F\,V}} + \langle \mathcal{S}_{zx}\mathcal{S}_{zx}\rangle _{\mathcal{F\,V}} \right ), \end{align}

where $\mathcal{S}_{xx}$ , $\mathcal{S}_{zz}$ , $\mathcal{S}_{xz}$ and $\mathcal{S}_{zx}$ can be calculated by (3.2). Figure 29 shows the relationship between the different average dissipation rate components of interfacial vortices and displays the fluid field snapshots at different stages. The dissipation process of the interfacial vortices can be divided into three stages in figure 29( $a$ ) by using the approximate equalisation of $\langle \epsilon _{\textit{NS}}\rangle _V$ and $\langle \epsilon _{\textit{SD}}\rangle _V$ as a basis for stage division.

1. (i) Stage I: the dissipation rate from shear deformation $\langle \epsilon _{\textit{SD}}\rangle _V$ is higher than that from normal strain $\langle \epsilon _{\textit{NS}}\rangle _V$ , with $\langle \epsilon _{\textit{SD}}\rangle _V$ gradually decreasing while $\langle \epsilon _{\textit{NS}}\rangle _V$ increases. The relationship between the two dissipation rates can be described by

   (7.13) \begin{equation} \frac {\textrm{d} \langle \epsilon _{\textit{NS}}\rangle _V}{\textrm{d} \langle \epsilon _{\textit{SD}}\rangle _V} \approx -1. \end{equation}
   This indicates a mutual restriction between the dissipation intensity brought on by shear deformation and that brought on by normal strain. At this time, the vortex structures at the interface appear as slender strips (figure 29 $c$ ). This dissipation relationship persists until the two dissipation rates approach equilibrium, i.e. $\langle \epsilon _{\textit{NS}}\rangle _V \approx \langle \epsilon _{\textit{SD}}\rangle _V$ , marking the end of the first stage.

2. (ii) Stage II: in this stage, $\langle \epsilon _{\textit{SD}}\rangle _V$ is less than $\langle \epsilon _{\textit{NS}}\rangle _V$ , and the interface vortex size keeps increasing. The difference between the major and minor axes of the elliptical vortex narrows and then expands (figure 29 $d$ – $f$ ), indicating that the vortex oblateness decreases and then increases. Here $\langle \epsilon _{\textit{SD}}\rangle _V$ first decreases and then increases, exhibiting a positive correlation with vortex oblateness. In contrast, $\langle \epsilon _{\textit{NS}}\rangle _V$ increases first and then decreases, which is negatively correlated with vortex oblateness. The time when stage II starts or ends ( $\langle \epsilon _{\textit{SD}}\rangle _V = \langle \epsilon _{\textit{NS}}\rangle _V$ ) may be related to the oblateness of the vortex. At the end of this stage, both $\langle \epsilon _{\textit{SD}}\rangle _V$ and $\langle \epsilon _{\textit{NS}}\rangle _V$ are lower than at the start, consistent with the variation of $\langle \epsilon \rangle _V$ shown in figure 28( $a$ ).

3. (iii) Stage III: the vortex elongates progressively (figure 29 $g$ ), with the shear deformation dissipation rate gradually increasing and the normal strain dissipation rate decreasing – opposite to the variation in stage I. Notably, the relationship given by (7.13) in stage I remains valid here.

Note that the relationship of (7.13) in stages I and III yields

(7.14) \begin{equation} \langle \epsilon _{\textit{NS}}\rangle _V + \langle \epsilon _{\textit{SD}}\rangle _V = const = \langle \epsilon \rangle _V. \end{equation}

It suggests the conservation of dissipation rate of the interfacial vortices. In addition, it can be seen from figure 29( $b$ ) that, during the development of the interfacial vortices, the streamwise and vertical components of normal strain dissipation rate are nearly equal, i.e.

(7.15) \begin{equation} \langle \epsilon _{xx}\rangle _V \approx \langle \epsilon _{zz}\rangle _V. \end{equation}

By combining (7.11), (7.13) and (7.15), it yields

(7.16) \begin{equation} \frac {\textrm{d} \langle \epsilon _{xx}\rangle _V}{\textrm{d} \langle \epsilon _{\textit{SD}}\rangle _V} \approx \frac {\textrm{d} \langle \epsilon _{zz}\rangle _V}{\textrm{d} \langle \epsilon _{\textit{SD}}\rangle _V} \approx -\frac {1}{2}. \end{equation}

In general, fluid shear deformation is the primary cause of fluid dissipation in the early and later stages of interfacial vortices, whereas normal strain dominates dissipation in the middle stage when the vortex size is relatively larger.

Temporal evolution of the average spanwise dissipation rate $\langle \epsilon _{\textit{SD}}\rangle _V$ of various selected interfacial vortices is presented in figure 30. To align all interfacial vortices on the same time axis, the horizontal axis origin is set at the moment when the average spanwise dissipation rate reaches its minimum value. It is noteworthy that the maximum vortex size coincides approximately with the point at which the average spanwise dissipation rate is at its lowest. Figure 30 shows that the average shear deformation dissipation rate $\langle \epsilon _{\textit{SD}}\rangle _V$ of vortices A, B and C first decreases and then increases. The $\langle \epsilon _{\textit{SD}}\rangle _V$ curves of these vortices exhibit a high degree of overlap, particularly in the vicinity of $t = t|_{\langle \epsilon _{\textit{SD}}\rangle _{V,min }}$ . This indicates that the majority of the fluids in the interfacial KH vortices have comparable fluid energy dissipation processes, with the exception of vortex A at the front. The start-up vortex A, the ending of the vortex sheet, undergoes a fluctuating change in the average shear deformation dissipation rate. Moreover, as can be seen in the inset of figure 30( $b$ ), the change rate of dissipated energy of each vortex $\textrm{d} E_{\epsilon ,SD} / {\textrm{d}}t$ is roughly synchronised with $\langle \epsilon _{\textit{SD}}\rangle _V$ , initially decreasing and subsequently increasing. This suggests a negative correlation between the vortex size and the dissipation of the interfacial vortices due to the shear deformation.

Figure 30.

Time evolution of average dissipation rate of turbulent kinetic energy $\langle \epsilon _{\textit{SD}}\rangle _{V}$ for typical KH vortices. Results are shown for ( $a$ ) ${Re} = 2426$ , ( $b$ ) ${Re} = 3030$ , ( $c$ ) ${Re} = 3431$ and ( $d$ ) ${Re} = 3836$ . Vortices A, B, C and D are four typical vortices along the streamwise direction, respectively, and vortex D develops at the TC front. The inset in figure 30( $b$ ) shows the time evolution of shear deformation dissipation change rate $\textrm{d}E_{\epsilon ,SD}/\textrm{d}t$ of the vortices. Vortices A–D in case ${Re} = 3030$ , shown in figure 5, correspond to vortices A–D at similar positions in other cases.