1. Introduction
Turbulent entrainment is the process by which a localised turbulent region draws in an outer non-turbulent (irrotational) fluid, leading to the expansion of the turbulent region. Entrainment plays a crucial role in various flows encountered in engineering and physics, and is also evident in many environmental flows, such as atmospheric boundary layers (Schols Reference Schols1984; Mahrt Reference Mahrt1999) and clouds (Mellado Reference Mellado2017). At a cloud edge, dry air from an external region is often mixed into the cloud (Sardina et al. Reference Sardina, Picano, Brandt and Caballero2015). Entrainment significantly influences turbulent transport processes, especially in scalar mixing (e.g. heat and substances) (Westerweel et al. Reference Westerweel, Fukushima, Pedersen and Hunt2009). Consequently, the control and enhancement of turbulent entrainment remain important topics in both science and engineering (da Silva et al. Reference da Silva, Hunt, Eames and Westerweel2014).
Over the past decades, extensive research has been conducted to understand the entrainment process, emphasising the turbulent/non-turbulent interface (TNTI) layer that separates turbulent from non-turbulent regions (Westerweel et al. Reference Westerweel, Fukushima, Pedersen and Hunt2005; Holzner & Lüthi Reference Holzner and Lüthi2011; Taveira et al. Reference Taveira, Diogo, Lopes and da Silva2013; van Reeuwijk & Holzner Reference van Reeuwijk and Holzner2014; Krug et al. Reference Krug, Holzner, Lüthi, Wolf, Kinzelbach and Tsinober2015; Jahanbakhshi & Madnia Reference Jahanbakhshi and Madnia2016; Watanabe et al. Reference Watanabe, da Silva, Sakai, Nagata and Hayase2016b), as summarised in da Silva et al. (Reference da Silva, Hunt, Eames and Westerweel2014). The TNTI layer appears in canonical turbulent flows, such as jets, wakes and mixing layers. The local entrainment process in many such flows is dominated by small-scale turbulent motions located at or near the TNTI layer. Larger-scale motions predominantly dictate the total entrainment rate by influencing the area of the interface (Mistry et al. Reference Mistry, Philip, Dawson and Marusic2016; Krug et al. Reference Krug, Holzner, Marusic and van Reeuwijk2017). Notably, the TNTI layer forms near small-scale vortices (vortex tubes) and shear layers (da Silva, Dos Reis & Pereira Reference da Silva, Dos Reis and Pereira2011; Elsinga & da Silva Reference Elsinga and da Silva2019; Hayashi, Watanabe & Nagata Reference Hayashi, Watanabe and Nagata2021b; Neamtu-Halic et al. Reference Neamtu-Halic, Mollicone, Van Reeuwijk and Holzner2021). Shear layers, originally termed vortex sheets due to their inherent vorticity, are sheet-like structures exhibiting shearing motion (Vincent & Meneguzzi Reference Vincent and Meneguzzi1994). Recent studies use the term shear layers to emphasise that their essential characteristic is shear rather than vorticity (Eisma et al. Reference Eisma, Westerweel, Ooms and Elsinga2015; Gul, Elsinga & Westerweel Reference Gul, Elsinga and Westerweel2020; Watanabe, Tanaka & Nagata Reference Watanabe, Tanaka and Nagata2020b; Fiscaletti, Buxton & Attili Reference Fiscaletti, Buxton and Attili2021). This study employs ‘shear layers’ to denote small-scale shear layers arising from velocity fluctuations, distinguishing them from those resulting from the mean velocity gradient, the latter being large-scale structures. The properties of the small-scale vortices and shear layers provide reasonable explanations for the entrainment process described by fluid particle movements and vorticity transport near the TNTI layer (Taveira & da Silva Reference Taveira and da Silva2014; Jahanbakhshi, Vaghefi & Madnia Reference Jahanbakhshi, Vaghefi and Madnia2015; Watanabe et al. Reference Watanabe, Jaulino, Taveira, da Silva, Nagata and Sakai2017; Elsinga & da Silva Reference Elsinga and da Silva2019; Hayashi et al. Reference Hayashi, Watanabe and Nagata2021b; Neamtu-Halic et al. Reference Neamtu-Halic, Mollicone, Van Reeuwijk and Holzner2021). These studies imply that vortices and shear layers within the TNTI layer play dominant roles in the entrainment process.
Compared to vortices, small-scale shear layers have received limited investigation due to challenges in their identification. Recent advancements in identifying shear layers based on novel mathematical treatments of velocity gradient tensors have enabled researchers to scrutinise their characteristics, as summarised herein (Kolář & Šístek Reference Kolář and Šístek2014). The shear layer is one of the smallest turbulent structures, characterised using the Kolmogorov scales, and exhibits biaxial straining motion (Elsinga & Marusic Reference Elsinga and Marusic2010; Watanabe et al. Reference Watanabe, Tanaka and Nagata2020b; Fiscaletti et al. Reference Fiscaletti, Buxton and Attili2021). The interplay between shear and biaxial strain leads to substantial enstrophy production and self-amplification of strain. Moreover, the velocity field induced by shear layers makes a more pronounced contribution to the energy cascade compared to that induced by vortex tubes (Enoki, Watanabe & Nagata Reference Enoki, Watanabe and Nagata2023). The stability of small-scale shear layers has also been investigated, given their role in vortex formation in turbulent regions. Prior to the investigation of small-scale shear layers, parallel shear flows had been studied theoretically as approximations of large-scale shear layers induced by a mean flow. The flow is inherently unstable against weak perturbations, leading to layer roll-up due to the Kelvin–Helmholtz instability. Linear stability theory suggests that parallel shear flows are most unstable in the presence of perturbations with a specific wavelength (Betchov & Szewczyk Reference Betchov and Szewczyk1963). Similarly, small-scale shear layers in turbulence generate vortices in a process reminiscent of the Kelvin–Helmholtz instability (Vincent & Meneguzzi Reference Vincent and Meneguzzi1994; Watanabe et al. Reference Watanabe, Tanaka and Nagata2020b). Recent studies of small-scale shear layers confirmed the optimal wavelength for the small-scale shear instability, which is approximately 30 times the Kolmogorov scale 
$\eta$ (Watanabe & Nagata Reference Watanabe and Nagata2023). When turbulence is subject to perturbations of this wavelength, the shear layers promptly collapse, resulting in a higher population of vortices arising from the shear instability. The evolution of shear layers is largely unaffected by perturbations with wavelengths much larger or smaller than 
$30\eta$. The optimal perturbation wavelength 
$30\eta$ was determined using two simulations (Watanabe & Nagata Reference Watanabe and Nagata2023). One is a direct numerical simulation (DNS) of turbulent flows, which examines the length and velocity scales of the shear layers. The other is a numerical simulation of a modelled shear layer based on a conditionally averaged flow field observed around the shear layers in DNS. The first has confirmed that shear layers have a typical thickness 
$\delta$ of 
$4\eta$ regardless of Reynolds numbers. Simulating modelled shear layers suggests that the optimal wavelength for shear instability is 
$7\delta$, which is consistent with the stability analysis of uniform shear layers (Lin & Corcos Reference Lin and Corcos1984). Because of the probability distribution of 
$\delta$, perturbations with wavelength approximately 
$30\eta$ effectively promote the instability of many shear layers.
Previous studies concerning the entrainment process and small-scale shear instability suggest the potential for enhancing entrainment by modulating small-scale shear layers. Given the efficacy of flow control strategies that exploit various flow instabilities to amplify the effects of weak disturbances (Cattafesta & Sheplak Reference Cattafesta III and Sheplak2011), it is plausible that manipulating small-scale shear layers can similarly lead to efficient flow control. In this study, numerical experiments are carried out to explore entrainment enhancement by stimulating small-scale shear layers, employing the DNS of a turbulent front perturbed by artificial velocity fluctuations with a wavelength comparable to the unstable mode of shear layers.
2. The DNS of a turbulent front evolving into a non-turbulent fluid
The DNS is conducted for a turbulent front evolving into a non-turbulent fluid without mean shear. The flow configuration aligns with previous works (Cimarelli et al. Reference Cimarelli, Cocconi, Frohnapfel and De Angelis2015; Silva, Zecchetto & da Silva Reference Silva, Zecchetto and da Silva2018; Watanabe, da Silva & Nagata Reference Watanabe, da Silva and Nagata2020a), which are outlined briefly herein. The governing equations are the incompressible Navier–Stokes equations, which are solved using an in-house finite-difference code based on the fractional step method (Watanabe et al. Reference Watanabe, Tanaka and Nagata2020b). The code employs a fourth-order fully conservative central difference and third-order Runge–Kutta schemes for spatial and temporal discretisation (Morinishi et al. Reference Morinishi, Lund, Vasilyev and Moin1998). The initial field is generated by embedding homogeneous isotropic turbulence (HIT) within a quiescent fluid. Snapshots from DNS databases of statistically steady HIT (Watanabe et al. Reference Watanabe, Tanaka and Nagata2020b) are used. The computational domain is set as a cube with side length 
$L$. Periodic boundary conditions are applied in three directions. Within the coordinate system 
$(x,y,z)$, 
$y=0$ is positioned at the domain centre, and 
$(x,z)=(0,0)$ is anchored at the corner of the 
$x$–
$z$ plane. During the initialisation phase, the velocity field 
$u_i$ of HIT is modified by the top-hat function 
$C(y)$:
\begin{equation} C(y)=0.5 + 0.5\tanh\left[\frac{4}{\varDelta_I}\left(1-\frac{2\,|y-L/2|}{L_T}\right)\right], \end{equation}
resulting in 
$C(y)\,u_i(x,y,z)$, with 
$L_T=L/3$ and 
$\varDelta _I=10\eta$, where 
$\eta$ is the Kolmogorov scale of HIT. Turbulence, centred at 
$y=0$ with initial width 
$L_T$, evolves into the surrounding non-turbulent fluid. By definition, 
$C=1$ in the turbulent region and 
$C=0$ in the non-turbulent region. Figure 1 shows a two-dimensional profile of vorticity magnitude 
$\omega =|\boldsymbol {\nabla } \times \boldsymbol {u}|$ after the evolution of turbulence from its initial state. The turbulent region exhibits large values of 
$\omega$, and grows by entraining non-turbulent fluid.
Figure 1. Distribution of vorticity magnitude 
$\omega$ in a turbulent front at 
$t=14\tau _{\eta 0}$ for R4. The white lines represent the irrotational boundary. The normalised vorticity magnitude 
$\hat {\omega }$ is defined as 
$\omega$ divided by the mean value at 
$y=0$.
The flow is statistically homogeneous in the 
$x$ and 
$y$ directions. Averages of flow variables, represented by 
$\bar {f}$, are determined as spatial averages over the 
$x$–
$z$ planes as functions of 
$y$ and time 
$t$. A fluctuating component of a variable 
$f$ is expressed as 
$f'=f-\bar {f}$. The turbulent Reynolds number is given by 
$Re_\lambda =u_{rms}\lambda /\nu$. Here, 
$\nu$ is the kinematic viscosity, 
$u_{rms}=\sqrt {\overline {u_{i}'u_{i}'/3}}$ is the root mean square (rms) of velocity fluctuations, and 
$\lambda =\sqrt {15\nu u_{rms}^2/\varepsilon }$ is the Taylor microscale. Throughout this paper, successive indices imply summation. The average kinetic energy dissipation rate is 
$\varepsilon =\overline {2\nu \boldsymbol{\mathsf{S}}_{ij}\boldsymbol{\mathsf{S}}_{ij}}$, and the rate-of-strain tensor is 
$\boldsymbol{\mathsf{S}}_{ij}=(\partial u_i/\partial x_j+\partial u_j/\partial x_i)/2$. The Kolmogorov length, time and velocity scales are defined as 
$\eta =(\nu ^3/\varepsilon )^{1/4}$, 
$\tau _{\eta }=(\nu /\varepsilon )^{1/2}$ and 
$u_\eta =(\nu \varepsilon )^{1/4}$, respectively. Hereafter, a subscript 0 indicates statistics at 
$y=0$ of the initial condition without perturbations, such as 
$\eta _0$ and 
$u_{\eta 0}$.
The evolution of the turbulent front is analysed under two distinct initial conditions. In the first case, the velocity field of HIT, denoted as 
$\boldsymbol {u}_{HIT}$, serves as the initial field, i.e. 
$\boldsymbol {u}=C\boldsymbol {u}_{HIT}$. In the second case, solenoidal velocity perturbations 
$\boldsymbol {u}_P$ are introduced to the velocity field in the turbulent region. The perturbed velocity field is formulated as 
$\boldsymbol {u}=C(\boldsymbol {u}_{HIT}+\boldsymbol {u}_P)$. A prior study addressed the impacts of 
$\boldsymbol {u}_P$ on small-scale shear instability in HIT, where the changes due to the promoted instability were discussed in terms of the number of vortices, energy dissipation, enstrophy production and strain self-amplification (Watanabe & Nagata Reference Watanabe and Nagata2023). Adopting this same methodology, 
$\boldsymbol {u}_P$ is described using sinusoidal functions: 
$\boldsymbol {u}_P=[ u_f\sin (2{\rm \pi} y/\lambda _f), u_f\sin (2{\rm \pi} z/\lambda _f), u_f\sin (2{\rm \pi} x/\lambda _f) ]$, where 
$u_f$ is the amplitude, and 
$\lambda _f$ is the wavelength. Different perturbation forms, such as random disturbances with a single length scale, were tested for HIT (Watanabe & Nagata Reference Watanabe and Nagata2023), showing that the evolution of perturbed HIT remains independent of perturbation type.
The parameters for DNS are detailed as follows. Five DNS databases of HIT with 
$Re_\lambda =43$, 72, 128, 202 and 296 are labelled as R1, R2, R3, R4 and R5, respectively. The numbers of grid points 
$N^3$ for these cases are 
$256^3$, 
$512^3$, 
$1024^3$, 
$2048^3$ and 
$4096^3$, ensuring a grid spacing smaller than 
$0.8\eta$. Cases without perturbations, i.e. 
$\boldsymbol {u}_P=(0,0,0)$, are simply designated as R1 to R5. Cases incorporating perturbations are designated as R
$n\varLambda a$ with 
$n=1,\ldots,5$, indicating the Reynolds number, where 
$\varLambda ={\rm S}$ or L distinguishes the wavelength 
$\lambda _f$, and 
$a=u_f/u_{\eta 0}$ is the normalised amplitude. The designation 
$\varLambda ={\rm S}$ corresponds to a perturbation of 
$\lambda _f=30\eta _0$, namely the optimal wavelength of small-scale shear instability (Watanabe & Nagata Reference Watanabe and Nagata2023), whereas 
$\varLambda ={\rm L}$ assumes a larger wavelength, 
$\lambda _f=140\eta _0$. The normalised amplitude 
$a$ ranges between 1.0 and 2.6. For instance, R2S14 considers turbulence at 
$Re_\lambda =72$ with a perturbation defined by 
$\lambda _f=30\eta _0$ and 
$u_f/u_{\eta 0}=1.4$. Most simulations adapt perturbations with 
$\lambda _f=30\eta _0$ (
$\varLambda ={\rm S}$). One additional case for R2 with 
$\varLambda ={\rm L}$ (
$\lambda _f=140\eta _0$), in which the perturbation is unlikely to influence the evolution of small-scale shear layers (Watanabe & Nagata Reference Watanabe and Nagata2023), is designated as R2L14. Each simulation is allowed to progress until 
$t=15\tau _{\eta 0}$. Three-dimensional data files for post-processes are saved at intervals close to 
$\tau _{\eta 0}$. Ensemble averages are obtained from 
$N_S$ simulations, each initialised with different snapshots of HIT. The values of 
$N_S$ are 10, 5, 3, 1 and 1, applied sequentially from case R1 to case R5. The perturbations are introduced to the initial velocity field of shear-free turbulent fronts. The Appendix examines the dependence of the flow evolution on the time at which perturbations are introduced.
The ratios between the integral scale 
$L_I$ and Kolmogorov scale 
$\eta$ in the HIT used for the initial condition are 37, 79, 191, 373 and 664, ascending from the lowest to the highest 
$Re_\lambda$ cases. Here, 
$L_I=(2k_T/3)^{3/2}/\varepsilon$ is evaluated with the turbulent kinetic energy (TKE) 
$k_T=\overline {u_iu_i}/2$ and its dissipation rate 
$\varepsilon$ (Rosales & Meneveau Reference Rosales and Meneveau2005). The integral scale consistently exceeds the thickness of the shear layer, which is approximately 
$4\eta$. The same DNS databases were used to explore the statistical properties of small-scale shear layers (Watanabe et al. Reference Watanabe, Tanaka and Nagata2020b). Notably, even at the smallest 
$Re_\lambda$ value with 
$L_I/\eta =37$, the characteristics of small-scale shear layers, when normalised by Kolmogorov scales, remain congruent with those observed at higher 
$Re_\lambda$ (Watanabe et al. Reference Watanabe, Tanaka and Nagata2020b). A key parameter considered herein is the ratio of the perturbation wavelength 
$\lambda _f$ to the Kolmogorov scale 
$\eta$. This study focuses on entrainment in shear-free turbulent fronts, chosen for their initial statistical homogeneity, which assures constant 
$\lambda _f/\eta$ in the turbulent region at a fixed value of 
$\lambda _f$. In contrast, turbulent shear flows, such as jets and mixing layers, display statistical inhomogeneity, leading to spatial variations in 
$\lambda _f/\eta$. Previous studies indicate that mean shear effects exert minimal influence on the scaling of the TNTI and small-scale shear layers (Silva et al. Reference Silva, Zecchetto and da Silva2018; Hayashi et al. Reference Hayashi, Watanabe and Nagata2021b), indicating that the findings from the shear-free turbulent front can be extended to more complex flows with mean shear.
The present DNS considers internal perturbations in the turbulent region. Both internal and external perturbations are applicable in laboratory experiments. Internal perturbations can be generated in turbulence by objects, such as cylinders and spheres in a mean flow, and settling particles, whose wakes produce velocity fluctuations at a certain scale (Britter, Hunt & Mumford Reference Britter, Hunt and Mumford1979; Nagata et al. Reference Nagata, Noguchi, Kusama, Nonomura, Komuro, Ando and Asai2020b; Kato, Takamure & Uchiyama Reference Kato, Takamure and Uchiyama2022). External perturbations can be introduced by fluidic actuators outside the turbulent region (Smith & Glezer Reference Smith and Glezer1998).
3. Entrainment analysis
The entrainment process is explored by examining the isosurface of vorticity magnitude 
$\omega$ at the outer edge of the turbulent front. As the vorticity decays over time with the evolution of the shear-free turbulent front, the normalised vorticity magnitude 
$\hat {\omega }$, defined as 
$\omega$ divided by the average 
$\bar {\omega }$ at 
$y=0$, is utilised to identify the turbulent region. Similar normalisations are commonly employed to examine the TNTI in turbulent shear flows, where the mean vorticity magnitude decreases as the flow evolves (Bisset, Hunt & Rogers Reference Bisset, Hunt and Rogers2002; Attili, Cristancho & Bisetti Reference Attili, Cristancho and Bisetti2014; Watanabe, Zhang & Nagata Reference Watanabe, Zhang and Nagata2018). Turbulent and non-turbulent regions are distinguishable by 
$\hat {\omega }>\omega _{th}$ and 
$\hat {\omega }\leq \omega _{th}$, respectively, using a threshold 
$\omega _{th}$. Given that 
$\omega$ rapidly diminishes to zero from the turbulent region across a thin TNTI layer, the identified turbulent region remains largely unaffected by the specific choice of 
$\omega _{th}$, provided that 
$\omega _{th}$ is chosen from a range for which the detected turbulent volume does not depend on 
$\omega _{th}$, as outlined in Taveira et al. (Reference Taveira, Diogo, Lopes and da Silva2013) and da Silva et al. (Reference da Silva, Hunt, Eames and Westerweel2014). This threshold insensitivity in the TNTI analysis has been documented for various flows (da Silva et al. Reference da Silva, Hunt, Eames and Westerweel2014; Jahanbakhshi et al. Reference Jahanbakhshi, Vaghefi and Madnia2015; Watanabe et al. Reference Watanabe, Zhang and Nagata2018). The present investigation adopts 
$\omega _{th}=0.01$, determined by assessing the 
$\omega _{th}$ dependence of the turbulent volume. Consequently, the surface bounding on the irrotational region, termed the irrotational boundary (Watanabe et al. Reference Watanabe, Sakai, Nagata, Ito and Hayase2015b; Zecchetto & da Silva Reference Zecchetto and da Silva2021), is demarcated as an isosurface of 
$\hat {\omega }=\omega _{th}$. The irrotational boundary defines the outer edge of the TNTI layer, which has a finite thickness. Figure 1 also shows the irrotational boundary as an iso-vorticity line with 
$\hat {\omega }=0.01$. The iso-vorticity line effectively separates the turbulent front characterised by large vorticity magnitude from the non-turbulent region with negligible vorticity.
The entrainment of the turbulent front is examined through analysis of the irrotational boundary using isosurface area density 
$\varSigma '$, as described in Blakeley, Olson & Riley (Reference Blakeley, Olson and Riley2022). The post-processing procedure is outlined briefly below, with a more in-depth exploration of 
$\varSigma '$ available in Blakeley et al. (Reference Blakeley, Olson and Riley2022), Blakeley, Olson & Riley (Reference Blakeley, Olson and Riley2023) and works referenced therein. In each snapshot, the detection function 
$X(x,y,z)$ of the turbulent region is defined such that 
$X=1$ when 
$\hat {\omega }(x,y,z)\leq \omega _{th}$, and 
$X=0$ when 
$\hat {\omega }>\omega _{th}$. The isosurface area density 
$\varSigma '$ is then calculated as 
$\varSigma '(x,y,z)=-|\boldsymbol {\nabla } \hat {\omega }|^{-1}\,\boldsymbol {\nabla } X\boldsymbol {\cdot } \boldsymbol {\nabla } \hat {\omega }$. The isosurface area 
$A$ is written as 
$A=\iiint _{\mathcal {V}} \varSigma ' \,{\rm d}\kern0.06em x\,{\rm d} y\,{\rm d} z$, where the integration spans the entire computational domain 
$\mathcal {V}$, while the volume of the turbulent region 
$V_T$ is expressed using 
$X$ as 
$V_T=\iiint _{\mathcal {V}} (1-X) \,{\rm d}\kern0.06em x\,{\rm d} y\,{\rm d} z$. The entrainment rate is defined as the temporal change in turbulent volume, represented by 
$\dot {V}_T={\rm d} V_T/{\rm d} t$. The mean propagation velocity of the irrotational boundary is then evaluated as 
$V_P=\dot {V}_T/A$. These metrics related to entrainment may be compared between turbulent fronts originating from perturbed and unperturbed initial fields. For a given quantity 
$Q$, the relative disparity between the two cases – for example, R2S14 and R2 – is evaluated as 
$\Delta Q= (Q_P-Q_U)/Q_U$, where 
$Q_P$ and 
$Q_U$ denote 
$Q$ in the perturbed and unperturbed cases, respectively.
4. Results and discussion
The impact of perturbations on the turbulent front is first explored through visualisation of the shear layers. The triple decomposition of the velocity gradient tensor (Kolář Reference Kolář2007; Nagata et al. Reference Nagata, Watanabe, Nagata and da Silva2020a) is applied to obtain the local intensities of rigid-body rotation 
$I_R$ and shearing motion 
$I_S$ (Watanabe & Nagata Reference Watanabe and Nagata2023). Vortex tubes and shear layers can be identified using 
$I_R$ and 
$I_S$, respectively. Figure 2(a) shows the temporal evolution of shear layers (white) and vortices (orange) near the irrotational boundary for R2, for which the initial field is derived from the original HIT without perturbations. Figure 2(b) shows the equivalent flow region and time range for R2S14, for which the perturbation wavelength corresponds to the unstable mode of small-scale shear instability. Both R2 and R2S14 initially observe the same shear layer, identified as ‘S’, at 
$t/\tau _{\eta 0}=2$. However, in the perturbed case, a vortex, designated as ‘V’, forms within the shear layer by 
$t/\tau _{\eta 0}=4$. Eventually, vortex formation occurs together with the breakdown of the shear layer at 
$t/\tau _{\eta 0}=6$. In contrast, vortex formation occurs more gradually in the unperturbed case. Weak perturbations with 
$\lambda _f=30\eta _0$ efficiently trigger shear instability at small scales near the TNTI layer, as previously observed for HIT (Watanabe & Nagata Reference Watanabe and Nagata2023).
Figure 2. Temporal evolution of shear layers near the irrotational boundary for (a) R2 (unperturbed case) and (b) R2S14 (perturbed case). The isosurfaces of 
$I_R/\overline {I_R}=6$ and 
$I_S/\overline {I_S}=1.5$ represent vortices (orange) and shear layers (white), respectively. The mean intensities are evaluated at 
$y=0$. From left to right, the time instances are 
$t/\tau _{\eta 0}=2$, 
$4$ and 6. The colour contour represents 
$I_S$.
In the case of HIT, the shear instability promoted by perturbations results in an increased number of vortices (Watanabe & Nagata Reference Watanabe and Nagata2023). This phenomenon is further confirmed for the shear-free turbulent front. Vortices are identified following the methodology presented in Watanabe & Nagata (Reference Watanabe and Nagata2023). Within this framework, the current study assumes that a grid point belongs to a vortex when the intensity of rigid-body rotation 
$I_R$ normalised by its time-dependent average at 
$y=0$, denoted by 
$\overline {I_R}$, surpasses a threshold 
$I_{th}$. Each separate region of vortex points constitutes a single vortex structure. Accordingly, the number of vortices, represented as 
$N_V$, is obtained as a function of time. Vortices are identified as regions with 
$I_R/\overline {I_R} > I_{th}$, with 
$I_{th}=3$. Figure 3 presents the temporal variations of the relative difference between perturbed and unperturbed cases, 
$\Delta N_V$, for the R2 series. Positive 
$\Delta N_V$ values signify an increase in vortices due to perturbation effects, negative values indicate a decrease, and 
$\Delta N_V=0$ suggests that perturbations do not influence the number of vortices. The plots of 
$\Delta N_V$ are therefore intended to examine the qualitative impacts of perturbations on vortices rather than to provide an exact number of vortices. Watanabe & Nagata (Reference Watanabe and Nagata2023) evaluated the threshold (
$I_{th}$) dependency of 
$N_V$, focusing on the perturbation response in HIT with the same methodology as in the present DNS. By plotting the number of vortices against 
$I_{th}$, they observed identical perturbation effects, namely an increase in 
$N_V$ due to promoted shear instability, regardless of 
$I_{th}$. They also evaluated the potential of different variables to detect vortices, e.g. a second invariant of 
$\partial u_i/\partial x_j$, confirming that an increase in vortices is observed irrespective of the vortex identification method used. For the turbulent front in figure 3, perturbations initially have no influence on vortices, as shown by 
$\Delta N_V=0$ at 
$t=0$. When perturbations exhibit the optimal wavelength of shear instability (R
$n$S
$a$), an increase in 
$\Delta N_V$ is observed, suggesting that more vortices are generated when shear instability is promoted. No such increase is apparent in R2L14, where the perturbation wavelength considerably exceeds 
$30\eta$. The R2S
$a$ cases with 
$a$ representing a normalised amplitude demonstrate that perturbations with a larger amplitude have more pronounced effects on the number of vortices. The decaying HIT corresponding to R2S10 in Watanabe & Nagata (Reference Watanabe and Nagata2023) may be compared with the shear-free turbulent front, showing that variations in 
$\Delta N_V$ are similar in both flows. The responses of small-scale shear layers to perturbations are consistent in both HIT and shear-free turbulent fronts. Watanabe & Nagata (Reference Watanabe and Nagata2023) also demonstrated that the increment in vortices remains similar for different Reynolds numbers in HIT.
Figure 3. Changes in the number of vortices 
$N_V$ by perturbations, represented by 
$\Delta N_V$. HIT refers to the simulation of decaying HIT with a Reynolds number and perturbations corresponding to R2S10 (Watanabe & Nagata Reference Watanabe and Nagata2023).
The ramifications of enhanced shear instability are explored for the entrainment rate, represented as 
$\dot {V}_T$. Figure 4(a) presents the temporal variations of the entrainment rate 
$\dot {V}_T$ for three cases. A positive value of 
$\dot {V}_T$ indicates the growth of the turbulent region. The entrainment rate is calculated using a second-order central difference on discrete data sets of 
$V_T(t)$. The entrainment rate in R2S14 becomes higher than in the other cases after 
$t/\tau _{\eta 0}\approx 2$. This amplified entrainment from 
$t/\tau _{\eta 0}\approx 2$ is virtually absent in R2L14 despite having the same perturbation amplitude as R2S14.
Figure 4. (a) Temporal variations of entrainment rate 
$\dot {V}_T$. (b) Changes due to perturbations in the entrainment rate 
$\dot {V}_T$, represented by 
$\Delta \dot {V}_T$.
Figure 4(b) assesses the relative difference in the entrainment rate between perturbed and unperturbed cases by depicting 
$\Delta \dot {V}_T(t)$. Some simulation cases are omitted from the figure for clarity, although they are presented later to assess the Reynolds number dependence. Introducing perturbations that stimulate shear instability amplifies the entrainment rate, as shown by 
$\Delta \dot {V}_T>0$ in the R
$n$S
$a$ cases. Notably, 
$\Delta \dot {V}_T$ reaches its maximum at approximately 
$t/\tau _{\eta 0}\approx 7$ irrespective of the Reynolds number. The increase in vortices due to the shear instability promoted by perturbations is also most prominent at 
$t/\tau _{\eta 0}\approx 5$, as shown in figure 3. This consistency suggests that the activity of small-scale vortices plays a significant role in entrainment. After 
$t/\tau _{\eta 0}\approx 7$, 
$\Delta \dot {V}_T$ begins to decline and the influence of perturbations diminishes with time. When the perturbation wavelength is much larger than the unstable mode wavelength, the entrainment rate exhibits minimal change, as evidenced by the small 
$\Delta \dot {V}_T$ for R2L14. Consequently, the entrainment is enhanced when the perturbations are characterised by a wavelength similar in scale to the unstable mode of small-scale shear layers.
After imposing impulsive perturbations at a given time, the flow evolves freely. Although these perturbations initially cause an increase in the entrainment rate, indicated by 
$\Delta \dot {V}_T>0$ in figure 4(b), a subsequent decrease of 
$\Delta \dot {V}_T$ occurs after reaching a peak, and the difference between the perturbed and unperturbed cases becomes small. This indicates that the perturbation effect is transient. Similarly, the impact on the number of vortices is also temporary, as evidenced by the analogous trend of 
$\Delta N_V$ in figure 3. Continuous or repetitive application of perturbations would more effectively sustain this enhanced entrainment. The transient enhancement of entrainment achieved by impulsive perturbations has important applications for chemically reacting flows. In situations where two reactants are introduced into separate flows, their mixing and subsequent reaction occur rapidly near the TNTI layer, especially when the reaction time scale is shorter than that of turbulence (Watanabe et al. Reference Watanabe, Sakai, Nagata, Terashima, Suzuki, Hayase and Ito2013, Reference Watanabe, Sakai, Nagata, Ito and Hayase2014, Reference Watanabe, Naito, Sakai, Nagata and Ito2015a; Jahanbakhshi & Madnia Reference Jahanbakhshi and Madnia2018; Ren et al. Reference Ren, Wang, Luo and Fan2023). Thus enhancing entrainment even for a short time could be leveraged to increase the reaction rate near the TNTI layer.
Entrainment is the process by which the non-turbulent fluid gains vorticity near the TNTI layer (Westerweel et al. Reference Westerweel, Fukushima, Pedersen and Hunt2005; Holzner & Lüthi Reference Holzner and Lüthi2011). A comparison of the enstrophy budget within the TNTI layer and around vortex tubes and shear layers suggests that these small-scale structures are relevant to the entrainment process described by vorticity transport (Watanabe et al. Reference Watanabe, Jaulino, Taveira, da Silva, Nagata and Sakai2017; Hayashi et al. Reference Hayashi, Watanabe and Nagata2021b; Neamtu-Halic et al. Reference Neamtu-Halic, Mollicone, Van Reeuwijk and Holzner2021). However, the mechanisms by which these structures cause entrainment locally has not yet been determined. A more straightforward way to understand entrainment is provided by Lagrangian viewpoints (Mathew & Basu Reference Mathew and Basu2002; Holzner & Lüthi Reference Holzner and Lüthi2011; Taveira et al. Reference Taveira, Diogo, Lopes and da Silva2013; Watanabe, da Silva & Nagata Reference Watanabe, da Silva and Nagata2016a; Watanabe et al. Reference Watanabe, da Silva, Sakai, Nagata and Hayase2016b). Entrained fluid particles pass across the TNTI layers to reach the turbulent core region. The shear layers within the TNTI layer are unstable, and their roll-up generates vortex tubes as shown in figure 2. Layer roll-up causes a significant mixing of fluids on both sides of the shear layer (Winant & Browand Reference Winant and Browand1974; Rogers & Moser Reference Rogers and Moser1992). One of the possible mechanisms by which small-scale structures cause entrainment is small-scale mixing associated with the roll-up of shear layers. The shear layers within the TNTI layer separate the turbulent and non-turbulent fluids (Elsinga & da Silva Reference Elsinga and da Silva2019; Hayashi et al. Reference Hayashi, Watanabe and Nagata2021b), which are expected to be mixed by shear instability. Another entrainment mechanism is directly related to the velocity distribution around vortex tubes and shear layers. Vortex tubes in turbulence form within an axial straining flow (Jiménez & Wray Reference Jiménez and Wray1998; da Silva et al. Reference da Silva, Dos Reis and Pereira2011), which induces an inward flow to the vortex core (Davidson Reference Davidson2004). It has been noted that this inward velocity of vortex tubes can initiate entrainment (Watanabe et al. Reference Watanabe, Jaulino, Taveira, da Silva, Nagata and Sakai2017) when the non-turbulent fluid is drawn into vortex tubes within the TNTI layer by the inward flow and is then further transferred into the turbulent core region by the rotating motion of the vortex. This explanation, based on the inward velocity of the straining flow, is also valid for shear layers within the TNTI layer because the shear layers undergo a biaxial strain with an inward velocity in the layer normal direction. The entrained fluid gains vorticity by viscous diffusion when it is transferred by the inward flow, and this entrainment process is consistent with Eulerian investigations of entrainment with vorticity transport.
The perturbation effect on the entrainment rate, described by 
$\Delta \dot {V}_T$, is similar to that of the number of vortices, 
$\Delta N_V$. This correlation between 
$\Delta \dot {V}_T$ and 
$\Delta N_V$ has an important implication for the dominant entrainment mechanism. The perturbations promote shear instability, resulting in a greater number of vortices in the flow. The increase ‘rate’ of 
$\Delta N_V$, 
$\partial \Delta N_V/\partial t$, provides a measure of the promotion of shear instability due to perturbations. As shown in figure 3, an increase in 
$\Delta N_V$ is observed for 
$t/\tau _\eta \leq 5$. The roll-up of shear layers in the perturbed cases is expected to occur more frequently than in the unperturbed cases for 
$t/\tau _\eta \leq 5$. If mixing due to the roll-up of shear layers dominates the entrainment process, then 
$\Delta \dot {V}_T$ should be maximised until 
$t/\tau _\eta \leq 5$. However, this contradicts the results for 
$\Delta \dot {V}_T$ shown in figure 4(b). Instead, 
$\Delta \dot {V}_T$ and 
$\Delta N_V$ reach their peaks at almost the same time. This suggests that the entrainment rate is related to the number of vortices 
$N_V$ rather than the increase rate of 
$N_V$ due to perturbations. It should be noted that 
$N_V$ is also correlated with the number of shear layers because the latter forms in the proximity of the former (Watanabe & Nagata Reference Watanabe and Nagata2023). Here, 
$N_V$ concerns the vortices throughout the entire turbulent region; however, 
$N_V$ can also be related to the number of vortices near the TNTI layer because vortices appear anywhere within the turbulent region. Indeed, the probability for vortex tubes and shear layers to appear at a given location does not change between the vicinity of the TNTI layer and the turbulent core region except within the viscous superlayer (da Silva et al. Reference da Silva, Dos Reis and Pereira2011; Hayashi et al. Reference Hayashi, Watanabe and Nagata2021b). The correlation between 
$\Delta \dot {V}_T$ and 
$\Delta N_V$ suggests that the presence of small-scale vortical structures near the TNTI layer plays an important role in entrainment, supporting a mechanism of entrainment linked to the inward flows of vortex tubes and shear layers as discussed above.
The entrainment rate is linked to two main factors: the surface area of the irrotational boundary 
$A$, and the mean propagation velocity 
$V_P$. Figure 5 examines the perturbation-induced variations in 
$A$ and 
$V_P$ with 
$\Delta A$ and 
$\Delta V_P$, respectively, both of which increase from 0 with time. The combined effect of these increments culminates in an overall enhancement of the entrainment rate 
$\dot {V}_T=V_PA$. An increase in the area is expected from surface stretching due to vortices generated by shear instability (Neamtu-Halic et al. Reference Neamtu-Halic, Krug, Mollicone, van Reeuwijk, Haller and Holzner2020).
Figure 5. Changes due to perturbations in the isosurface area 
$A$ and mean propagation velocity 
$V_P$, represented by 
$\Delta A$ and 
$\Delta V_P$, respectively.
Figure 6(a) shows the peak values of 
$\Delta \dot {V}_T$ in figure 4(b) plotted against the perturbation amplitude normalised by 
$u_{\eta 0}$ for all Reynolds number cases with 
$\varLambda =\mathrm {S}$ (
$\lambda _f=30\eta _0$). A larger amplitude results in greater enhancement of the entrainment. Notably, when normalised by the Kolmogorov velocity, the perturbation effects remain largely independent of the Reynolds numbers. Values of 
$\Delta \dot {V}_T$ peak at approximately 
$7\tau _{\eta 0}$, which is also identical for all Reynolds numbers considered here. Thus the enhancement of entrainment by small-scale shear instability occurs at the Kolmogorov scales.
Figure 6. (a) Maximum values of the increase ratio of entrainment rate due to perturbations, represented by 
$\Delta \dot {V}_T$, plotted against the normalised perturbation amplitude 
$a=u_f/u_{\eta 0}$. (b) Maximum values of 
$\Delta \dot {V}_T$ plotted against the relative increase in TKE due to the perturbations, represented by 
$\Delta k_T$.
The perturbations applied in the DNS presented herein result in an increase in the initial TKE, denoted as 
$k_T=\overline {u_iu_i}/2$. The relative increase is evaluated as 
$\Delta k_T=(k_{TP}-k_{T0})/k_{T0}$, wherein 
$k_{TP}$ and 
$k_{T0}$ correspond to the initial TKE at 
$y=0$ for perturbed and unperturbed cases. Figure 6(b) displays the maximal 
$\Delta \dot {V}_T$ values for each case against 
$\Delta k_T$. The same entrainment enhancement can be achieved with weaker perturbations, i.e. smaller 
$\Delta k_T$, at a higher Reynolds number. An approximate 10 % increase in the entrainment rate (
$a=2.2$) occurs with increases in the TKE of 2 % and 12 % for R5 and R1, respectively. The line in figure 6(a) illustrates an empirical fit using a quadratic function, represented by 
$\Delta \dot {V}_T=C_1a^2$, with 
$C_1=1.75$, determined using a least squares method. This fitting function is introduced only for discussing the Reynolds number dependence for the range of 
$a$ considered in this study, and not for establishing the precise relationship between 
$\Delta \dot {V}_T$ and 
$a$. In fully developed turbulent regions, a scaling 
$u_{\eta 0}^2/k_{T0}\sim Re_\lambda ^{-1}$ is anticipated (Pope Reference Pope2000); hence this empirical relation for 
$\Delta \dot {V}_T$ suggests a scaling 
$\Delta \dot {V}_T\sim (u_f^{2}/k_{T0})\,Re_\lambda$. For the present test case, the energy introduced by the perturbation required to achieve equivalent entrainment enhancement (
$\Delta \dot {V}_T$) decreases as 
$Re_\lambda ^{-1}$ with increasing Reynolds numbers. Case R2L14 investigates the influence of perturbations with a wavelength significantly larger than the unstable mode wavelength of shear layers. Under such conditions, perturbations exert minimal influence on small-scale vortices and entrainment rate, as shown in figures 3 and 4(b). A weak influence of large-scale perturbations was also reported for vortices in HIT (Watanabe & Nagata Reference Watanabe and Nagata2023). The perturbations in R2L14 lead to an increase in the initial kinetic energy by 
$\Delta k_T=1.3\,\%$. This increment is similar to that in R5S18, where 
$\Delta k_T=1.5\,\%$. For R5S18, the entrainment rate increases by approximately 6 %. In the numerical set-up used for this study, the original HIT exhibits TKE approximately 
$k_{T0}\approx 1.5$ (arbitrary units), irrespective of the Reynolds number. Therefore, the similarity between the 
$\Delta k_T$ values in these two cases indicates that their initial kinetic energies are also comparable. However, the influence of perturbations on the vortices and entrainment is evident when the wavelength aligns with the unstable mode of shear layers. This comparison suggests that the kinetic energy added by perturbations alone is not a significant factor in enhancing entrainment; rather, the wavelength of the perturbations plays a crucial role. Additionally, no clear correlation is observed between the entrainment rate and 
$\Delta k_T$ across different Reynolds numbers in figure 6(b). The kinetic energy increment alone does not dictate the enhancement in entrainment rate. Instead, the perturbation amplitude (analogous to the kinetic energy of perturbations) normalised by the Kolmogorov velocity determines the increase in entrainment rate. This may be related to the scaling of shear layers. The velocity jump across the shear layers is also determined by the Kolmogorov velocity (Watanabe et al. Reference Watanabe, Tanaka and Nagata2020b; Fiscaletti et al. Reference Fiscaletti, Buxton and Attili2021; Hayashi, Watanabe & Nagata Reference Hayashi, Watanabe and Nagata2021a).
5. Conclusions
The numerical experiments reported in this study demonstrate that entrainment rate can be amplified significantly by introducing perturbations with wavelengths corresponding to the unstable mode of small-scale shear layers, i.e. at approximately 
$30$ times the Kolmogorov scale. These perturbations efficiently stimulate the shear layers near the TNTI layer, and accelerate vortex formations due to instability. The effects of impulsive perturbations on entrainment and vortices are transient, diminishing over time. The Kolmogorov scales primarily govern the process underlying the entrainment enhancement. Even weak perturbations can yield substantial enhancements in entrainment via small-scale shear instability, especially at higher Reynolds numbers, due to the smaller value of the corresponding Kolmogorov velocity. The statistical behaviour of small-scale shear layers is minimally influenced by flow types (Elsinga & Marusic Reference Elsinga and Marusic2010; Watanabe et al. Reference Watanabe, Tanaka and Nagata2020b; Fiscaletti et al. Reference Fiscaletti, Buxton and Attili2021; Hayashi et al. Reference Hayashi, Watanabe and Nagata2021a), suggesting that this method of enhancement is applicable to various turbulent flows. Furthermore, turbulent flows are filled with many small-scale shear layers, which exist anywhere in turbulent regions (Horiuti & Takagi Reference Horiuti and Takagi2005; Buxton & Ganapathisubramani Reference Buxton and Ganapathisubramani2010; Pirozzoli, Bernardini & Grasso Reference Pirozzoli, Bernardini and Grasso2010; Nagata et al. Reference Nagata, Watanabe, Nagata and da Silva2020a; Fiscaletti et al. Reference Fiscaletti, Buxton and Attili2021; Hayashi et al. Reference Hayashi, Watanabe and Nagata2021a). Intricate sensing mechanisms are not required to locate and manipulate these shear layers, unlike other active flow control techniques that rely on identifying specific turbulent structures. Although the results reported herein are derived from idealised numerical experiments conducted under controlled conditions, the resulting features of small-scale shear instability have broad practical applications, from engineering to environmental flows. Although additional research, including experimental validation, is required for future applications, the present study introduces a new and widely applicable framework for flow control that leverages small-scale shear instabilities triggered by weak and small-scale disturbances.
Acknowledgements
Numerical simulations presented in this paper were performed using the high-performance computing systems at the Japan Agency for Marine-Earth Science and Technology and Nagoya University. This work was also supported by a Collaborative Research Project on Computer Science with High-Performance Computing in Nagoya University.
Funding
This work was supported by JSPS KAKENHI grant no. JP22K03903.
Declaration of interests
The author reports no conflict of interest.
Data availability statement
The data supporting the findings of this study are available from the corresponding author upon reasonable request.
Appendix. Effects of the initial transition of a shear-free turbulent front on entrainment enhancement
The DNS of the shear-free turbulent front conducted in this study was initialised using a velocity field featuring an artificial interface generated by a truncated velocity field defined by (2.1). The flow evolution at an early time is influenced by the initial conditions, which may affect the observed entrainment enhancement. This appendix addresses the dependence of entrainment enhancement on the time at which perturbations are imposed. To this end, the DNS is performed for two additional cases involving the transfer of a passive scalar 
$\phi$, whose evolution is governed by an advection–diffusion equation (Watanabe et al. Reference Watanabe, Sakai, Nagata, Ito and Hayase2015b). The Schmidt number 
$Sc = \nu /D$ is set to 1, where 
$D$ represents the diffusivity coefficient. The initial distribution of the passive scalar 
$\phi$ is specified as 
$\phi (y) = C(y)$, which is 1 in the turbulent region and 0 in the non-turbulent region. The passive scalar with 
$Sc=1$ serves as the marker of the turbulent region. The interfaces identified using vorticity and passive scalar are consistent, exhibiting minimal differences in both the location and the local curvature of the isosurfaces (Gampert et al. Reference Gampert, Boschung, Hennig, Gauding and Peters2014; Watanabe et al. Reference Watanabe, Zhang and Nagata2018). Perturbations are introduced into the turbulent core region with 
$\phi >0.5$ after the flow has evolved from the initial condition for a time 
$T$. The Kolmogorov scales corresponding to the time at which the perturbations are added are denoted by 
$\eta _P$, 
$\tau _{\eta P}$ and 
$u_{\eta P}$. These scales are evaluated by considering a volume average of the energy dissipation rate in the turbulent region identified with the vorticity magnitude. The simulations in this appendix are conducted for R2, 
$\lambda _f=30\eta _P$ and 
$u_f=1.4u_{\eta P}$. Based on the integral time scale 
$\tau _I = L_I / \sqrt {2k_T/3}$ defined using the TKE 
$k_T$ and the characteristic length scale of large-scale motions 
$L_I$, 
$T$ is determined to be either 0 or 
$0.51\tau _I$ (
$9.4\tau _\eta$). Here, 
$\tau _I$ and 
$\tau _\eta$ are evaluated for the original HIT used for the initial conditions of the shear-free turbulent front. Figure 7 visualises the shear-free turbulent front at 
$t=9.4\tau _\eta$. The flow evolves until this time instance, and then the sinusoidal perturbations described in § 2 are introduced into the region with 
$\phi >0.5$. In this appendix, the results of a single simulation for each case are presented without taking ensemble averages of different simulations.
Figure 7. The shear-free turbulent front visualised by passive scalar 
$\phi$ at 
$t= 9.4\tau _\eta$ in R2. The coordinates are normalised by the Kolmogorov scale 
$\eta$ at 
$y=0$.
Figure 8(a) plots the entrainment rate 
$\dot {V}_T$ as a function of 
$t-T$, which represents the time after adding the perturbation. The results are compared between the perturbed and unperturbed cases. Regardless of the time 
$T$ at which the perturbations are added, the perturbed cases exhibit larger 
$\dot {V}_T$ values than the corresponding unperturbed cases. The choice of 
$T$ influences the temporal evolution of 
$\dot {V}_T$ because of the transient behaviour of the shear-free turbulent front at early times. Figure 8(b) illustrates changes in entrainment rate due to perturbations, denoted as 
$\Delta \dot {V}_T$. The temporal variation of 
$\Delta \dot {V}_T$ is quantitatively similar for all cases. Entrainment enhancement is observed for both 
$T$ values, and is not influenced by the initial transient behaviour of the shear-free turbulent front.
Figure 8. (a) Entrainment rate 
$\dot {V}_T$ plotted as a function of time after adding perturbations, 
$t-T$. The results are compared between the perturbed (
$u_f=1.4u_{\eta P}$) and unperturbed (
$u_f=0$) cases. (b) Changes due to perturbations in the entrainment rate, represented by 
$\Delta \dot {V}_T$.
