The response of small-scale shear layers to perturbations in turbulence

Abstract

The perturbation response of small-scale shear layers in turbulence is investigated with direct numerical simulations (DNS). The analysis of shear layers in isotropic turbulence suggests that the typical layer thickness is about four times the Kolmogorov scale $\eta$. Response for sinusoidal perturbations is investigated for an isolated shear layer, which models a mean flow around the shear layers in turbulence. The vortex formation in the shear layer is optimally promoted by the perturbation whose wavelength divided by the layer thickness is about 7. These results indicate that the small-scale shear instability in turbulence is efficiently promoted by velocity fluctuations with a wavelength of about $30\eta$. Furthermore, DNS are carried out for decaying turbulence initialised by the artificially modified velocity field of isotropic turbulence. The vortex formation from shear layers is accelerated under the influence of external perturbations with the efficient wavelength to promote the instability. When velocity fluctuations with this wavelength are eliminated by a band-cut filter, the shear layers tend to persist for a long time without producing vortices. These behaviours affect the number of vortices in turbulence, which increases and decreases when velocity perturbations with the unstable wavelength of the instability are artificially amplified and damped, respectively. The increase in the number of vortices results in the enhancement of kinetic energy dissipation, enstrophy production and strain self-amplification. These results indicate that the perturbation response of shear layers is important in the small-scale dynamics of turbulence as well as the modulation of small-scale turbulent motions by external disturbance.

1. Introduction

Turbulent flows are characterised by fluid motions with a wide range of length scales. Large-scale turbulent motions possess a large part of kinetic energy, which is transferred to small scales by inertial processes and eventually dissipated by viscous effects (Kolmogorov 1941; Pope 2000; Davidson 2004). The scale dependence of turbulence has been studied extensively in previous studies. The characteristics of large-scale turbulent motions depend strongly on the process by which turbulence is generated. On the other hand, small-scale motions are considered to have universal statistical properties which hardly depend on flows (Monin & Yaglom 1975; Sreenivasan & Antonia 1997). Turbulence is also studied in terms of flow structures, which are often identified as coherent patterns in flow visualisation. As also found for the statistical properties, large-scale structures are flow-dependent whereas small-scale structures are more universal (Brown & Roshko 1974; Yule 1978; Jiménez et al. 1993; da Silva, Dos Reis & Pereira 2011; Philip & Marusic 2012).

Small-scale properties of turbulence are often studied with a velocity gradient tensor $\boldsymbol {\nabla } \boldsymbol {u}$, whose component is denoted by $\partial u_i/\partial x_j$. Hereafter, a component of vectors and matrices is denoted by subscripts, e.g. $(\boldsymbol {\nabla } {\boldsymbol {u}})_{ij}=\partial u_i/\partial x_j$. The decomposition of $\boldsymbol {\nabla } \boldsymbol {u}$ into symmetric and antisymmetric parts as $(\boldsymbol {\nabla } {\boldsymbol {u}})_{ij}={\mathsf{S}}_{ij}+\textsf{$\mathit{\Omega}$} _{ij}$ is widely used to study incompressible turbulence, where ${\mathsf{S}}_{ij}=(\partial u_i/\partial x_j+\partial u_j/\partial x_i)/2$ is the rate-of-strain tensor and $\textsf{$\mathit{\Omega}$} _{ij}=(\partial u_i/\partial x_j-\partial u_j/\partial x_i)/2$ is the rate-of-rotation tensor. One of the most famous small-scale structures in turbulence is a vortex tube with rotating motion. The radius of vortices is about four to six times the Kolmogorov scale in various turbulent flows (Jiménez et al. 1993; Jiménez & Wray 1998; Kida & Miura 1998; Tanahashi, Iwase & Miyauchi 2001; Kang, Tanahashi & Miyauchi 2007; Mouri, Hori & Kawashima 2007; Ganapathisubramani, Lakshminarasimhan & Clemens 2008; da Silva et al. 2011; Jahanbakhshi, Vaghefi & Madnia 2015; Ghira, Elsinga & da Silva 2022). Another small-scale structure observed in turbulence is a vortex sheet. The vortex sheet is a thin layer with intense shear. Therefore, the vortex sheets are also called shear layers in recent studies because the structures are identified with shear intensity or shear vorticity (Eisma et al. 2015; Nagata et al. 2020). The small-scale shear layers were originally identified in turbulence as regions with moderately large enstrophy $\omega ^2/2=\textsf{$\mathit{\Omega}$} _{ij}\textsf{$\mathit{\Omega}$} _{ij}$ when the enstrophy field was visualised to detect vortical structures (Jiménez et al. 1993). Because of a strong correlation between vorticity and strain in shear layers, they are also identified with an eigenvalue of the second-order velocity gradient tensor ${\mathsf{A}}_{ij}={\mathsf{S}}_{ik}\textsf{$\mathit{\Omega}$} _{kj}+{\mathsf{S}}_{jk}\textsf{$\mathit{\Omega}$} _{ki}$ (Horiuti & Takagi 2005). The eigenvalue of ${\mathsf{A}}_{ij}$ is useful in identifying the location of shear layers and is used to investigate the flow statistics in the regions occupied by the shear layers, which are shown to be dynamically important structures related to the budget of turbulent kinetic energy and enstrophy (Buxton & Ganapathisubramani 2010; Pirozzoli, Bernardini & Grasso 2010; Nakamura, Watanabe & Nagata 2023).

Both vortex tubes and shear layers are perceivable in an enstrophy profile because both motions of rigid-body rotation and shear contribute to a vorticity vector ${\boldsymbol {\omega }}=\boldsymbol {\nabla } \times {\boldsymbol {u}}$. Recent studies of vortex identification schemes have developed new decompositions of $\boldsymbol {\nabla } \boldsymbol {u}$ by which the contribution of shearing motion is extracted from $\boldsymbol {\nabla } {\boldsymbol {u}}$ or ${\boldsymbol {\omega }}$. Such examples are the triple decomposition and the Rortex-based decomposition (Kolář 2007; Liu et al. 2018). To distinguish vortex tubes and shear layers, the triple decomposition extracts a component of shearing motion, $\boldsymbol {\nabla } \boldsymbol {u}_S$, from $\boldsymbol {\nabla } \boldsymbol {u}$. The vorticity vector defined with the residual component $\boldsymbol {\nabla } \boldsymbol {u}-\boldsymbol {\nabla } \boldsymbol {u}_S$ contains only rigid-body rotation. Thus, the magnitudes of vorticity vectors defined for $\boldsymbol {\nabla } \boldsymbol {u}_S$ and $\boldsymbol {\nabla } \boldsymbol {u}-\boldsymbol {\nabla } \boldsymbol {u}_S$ can identify shear layers and vortex tubes, respectively. The details of the triple decomposition can be found in Kolář (2007), which is briefly explained in § 2.2. The shear-layer identification in turbulent flows has been reported with criteria based on the norm or vorticity of $\boldsymbol {\nabla } \boldsymbol {u}_S$ (Eisma et al. 2015; Watanabe, Tanaka & Nagata 2020; Fiscaletti, Buxton & Attili 2021). Gul, Elsinga & Westerweel (2020) have compared shear layers detected with the triple decomposition and the eigenvalue of ${\mathsf{A}}_{ij}$ in a turbulent pipe flow and have found that the same flow regions with intense shear can be identified with both quantities. However, unlike the eigenvalue of ${\mathsf{A}}_{ij}$, the components of $\boldsymbol {\nabla } \boldsymbol {u}_S$ provide the information of shear orientation in addition to the location of shear layers. The mean flow field around shear layers in isotropic turbulence has been assessed in the reference frame defined with the shear orientation, which is identified with the triple decomposition (Watanabe & Nagata 2022). The shear layers in turbulence are formed in a region of biaxial strain. The shear layers are small-scale structures whose mean length scales are characterised by the Kolmogorov scale in isotropic turbulence and turbulent free shear flows (Watanabe et al. 2020; Fiscaletti et al. 2021; Hayashi, Watanabe & Nagata 2021a). The enstrophy production and strain self-amplification actively occur in the shear layers because of the interaction between the shear and biaxial strain (Watanabe et al. 2020; Hayashi et al. 2021a). Similar layer structures with intense shear have been observed between uniform momentum zones in wall-bounded shear flows (Eisma et al. 2015; Fan et al. 2019; Gul et al. 2020; Chen, Chung & Wan 2021). An experimental study on a turbulent boundary layer has shown that the thickness of these layers scales with the Taylor microscale (Eisma et al. 2015) while the thickness of small-scale shear layers found in isotropic turbulence, jets and mixing layers scales with the Kolmogorov scale. Therefore, it is not clear at this moment if the shear layers that separate the uniform momentum zones are the same as the small-scale shear layers found in turbulence far from the wall although both layer structures are characterised by intense shear.

Flow visualisation has revealed that a roll-up of small-scale shear layers results in the formation of a vortex tube by a similar process to Kelvin–Helmholtz instability (Vincent & Meneguzzi 1994; Passot et al. 1995; Watanabe et al. 2020). The dynamics of the roll-up process in a shear flow are often studied with an ideal model of a shear layer, which is initially uniform in the flow direction. These studies mostly concern large-scale shear flows, which are observed in plane mixing layers, jets and separated flows although they might be relevant to the instability of small-scale shear layers in turbulence. The instability of a parallel uniform shear flow has been studied with a linear stability theory, which relies on a linearised equation for small pertubations (Miles 1961; Betchov & Szewczyk 1963; Maslowe & Thompson 1971). Although the linear model of the shear layer is usually applied to an early stage of the layer development, the predicted initial growth rate of the instability is consistent with numerical simulations of Navier–Stokes equations (Patnaik, Sherman & Corcos 1976). Other studies also consider a parallel shear flow in a biaxial strain field, which is also known as Burgers’ vortex (Lin & Corcos 1984; Beronov & Kida 1996). Both linear and nonlinear analyses of these uniform shear flows have confirmed that the growth of the instability strongly depends on the wavelength of perturbations. When the wavelength matches the most unstable mode of the instability, the perturbed shear layer rapidly develops and results in a fast roll-up. This wavelength dependence is often taken into account in numerical simulations of turbulent shear flows, which rapidly develop from a laminar state with the perturbations of the most unstable mode (Moser & Rogers 1991; Caulfield & Peltier 1994; Vreman, Geurts & Kuerten 1995; Abraham & Magi 1997; Le Ribault, Sarkar & Stanley 1999; Okong'o & Bellan 2002; Stanley, Sarkar & Mellado 2002). The dependence of the instability on the perturbation wavelength has also proved useful in flow control, by which a flow is forced by actuators or passive objects. The flow is often efficiently modulated when an external disturbance matches the frequency band of instabilities (Crow & Champagne 1971; Oster & Wygnanski 1982). The application of this flow control strategy has been reported for various turbulent shear flows, which include canonical free shear flows and more practical flows such as a flow over an airfoil (Ghoniem & Ng 1987; Gutmark, Schadow & Yu 1995; Soteriou & Ghoniem 1995; Brehm, Gross & Fasel 2013; Samimy, Webb & Crawley 2018)

One of the interesting questions is whether small-scale shear layers in turbulence also have a wavenumber dependence in the response to perturbations, as also observed for uniform shear layers considered as numerical and analytical models and large-scale shear layers in laboratory experiments. Considering the application of the results of stability analysis in active flow control of large-scale shear flows, one may expect that the response of small-scale shear layers to perturbations will be important in the studies of modulation of turbulence. For example, small-scale perturbations can be introduced in particle-laden flows by the interaction between small particles and turbulence (Maxey 2017). As also considered for active flow control of turbulent shear flows, the perturbation response of small-scale shear layers can be important in establishing the strategy of active flow control of small-scale turbulent motions (Buzzicotti, Biferale & Toschi 2020; Cazaubiel et al. 2021). In addition, it can also be important in understanding the interaction of turbulent motions with different scales. The roll-up of the shear layers may be understood as the interaction between small-scale shearing motion and velocity fluctuations at larger scales. The interaction between different scales of motions is often investigated in wavenumber space (Yeung & Brasseur 1991; Ohkitani & Kida 1992; Aluie & Eyink 2009). However, some turbulent structures are difficult to identify in wavenumber space because the structures are not defined solely by their scales. For this reason, a large number of publications have been devoted to the investigation of turbulent structures in physical space. It is sometimes considered that the structural approach to turbulence is less rigorous than the statistical approach because the former mostly relies on flow visualisation and thresholds applied to flow variables. However, both approaches are important in understanding turbulence from different aspects as discussed in Tsinober (2009). For example, small-scale intermittency can be discussed with the scaling exponents for structure functions (Kholmyansky, Tsinober & Yorish 2001) while the spatiotemporal distributions of variables defined with a velocity gradient tensor often explain these statistical behaviour from the viewpoint of turbulent structures (Siggia 1981; Jiménez et al. 1993). The perturbation response of small-scale shear layers observed in physical space may shed light on a physical process that underlies the interaction of different scales in wavenumber space.

Most fundamental studies on shear instability assume a uniform shear layer, whose thickness is much smaller than the length scale in the other directions. However, visualisation of shear layers in turbulence has suggested that this assumption is not valid as they have a pancake or ribbon shape with a finite aspect ratio (Buxton & Ganapathisubramani 2010; Bhatt & Tsuji 2021; Fiscaletti et al. 2021; Hayashi et al. 2021a). For this reason, it is not clear how the stability analysis of shear flows with an infinite aspect ratio is related to the evolution of small-scale shear layers in turbulence. This might be because early studies of small-scale shear layers have relied on flow visualisation, which does not provide important parameters in the stability analysis, such as the probability distribution of thickness and velocity jump of the shear layers. However, recent studies have overcome these difficulties with the new statistical analysis of small-scale shear layers based on the conditional averaging procedure in a local reference frame defined for each shear layer (Eisma et al. 2015; Fiscaletti et al. 2021; Watanabe et al. 2020).

This study aims to investigate the response of small-scale shear layers in turbulence to perturbations, which may be relevant to vortex formation due to the instability in turbulence. Three kinds of numerical simulations are carried out for this purpose. First, the characteristics of shear layers are investigated with direct numerical simulations (DNS) of statistical steady homogeneous isotropic turbulence to reveal the variations of length and velocity scales and shear Reynolds number among shear layers. The analysis of shear layers in homogeneous isotropic turbulence is presented in § 2. Based on the results for isotropic turbulence, numerical simulations for an isolated shear layer with a finite aspect ratio subject to a biaxial strain are performed in § 3 to examine the development of the perturbed shear layer. The model follows Burgers’ vortex layer except for the aspect ratio and approximates a mean flow observed around the shear layers in turbulence. Furthermore, § 4 presents DNS of decaying turbulence which is initialised with an artificially modified velocity field of isotropic turbulence. The DNS is carried out to investigate the response of shear layers to external perturbations or velocity fluctuations that internally exist in turbulence. For the former case, the initial field of DNS is obtained by superimposing perturbations with a certain length scale on isotropic turbulence. The latter considers decaying turbulence initialised with a band-cut filtered velocity field of isotropic turbulence. The length scale of the perturbations and filter is determined based on the other two simulations carried out to examine the shear layer characteristics and the response of the isolated shear layer to perturbations. These simulations are compared with decaying turbulence of original isotropic turbulence to confirm whether vortex formation in turbulence can be promoted by velocity fluctuations with a wavelength within the shear instability band. Finally, the paper is summarised in § 5.

2. Characteristics of small-scale shear layers in turbulence

2.1. DNS of statistically steady homogeneous isotropic turbulence

DNS databases of statistically steady homogeneous isotropic turbulence are analysed to relate the small-scale shear layers in turbulence to a model of an isolated shear layer in § 3, where some parameters are determined based on the results presented in this section. The DNS databases are the same as those used in our previous studies (Watanabe et al. 2020; Watanabe, Tanaka & Nagata 2021; Watanabe & Nagata 2022). Turbulence is sustained at a statistically steady state with a linear forcing scheme (Carroll & Blanquart 2013) The governing equations are the Navier–Stokes equations for an incompressible fluid, which are written as

(2.1)$$\begin{gather} \frac{\partial u_j}{\partial x_j}=0, \end{gather}$$

(2.2)$$\begin{gather}\frac{\partial u_i}{\partial t }+u_j\frac{\partial u_i}{\partial x_j}=-\frac{1}{\rho}\frac{\partial p}{\partial x_i} +\nu\frac{\partial^2 u_i}{\partial x_{j}^2}+f_i, \end{gather}$$

where $u_{i}$ is the velocity vector, $p$ is the pressure, $\rho$ is the constant density, $\nu$ is the kinematic viscosity and $f_i$ is the forcing term (Carroll & Blanquart 2013). The DNS code is based on the fractional step method and employs a fourth-order fully conservative finite-difference scheme for spatial discretisation (Morinishi et al. 1998) and a third-order and low-storage Runge–Kutta method for temporal advancement (Spalart, Moser & Rogers 1991). The Poisson equation for pressure is solved with the biconjugate gradient stabilised method.

Table 1 summarises the parameters of DNS databases. Here, the statistics are defined with volume averages in the computational domain and ensemble averages of different time instances. An average of a variable $f$ is denoted by $\langle\, f\rangle$. The spatial resolution $\varDelta$ is about 0.8 times the Kolmogorov scale $\eta =(\nu ^3/\varepsilon _0)^{1/4}$, where $\varepsilon _0=\langle 2\nu {\mathsf{S}}_{ij}{\mathsf{S}}_{ij}\rangle$ is the averaged kinetic energy dissipation rate. Table 1 includes the Reynolds numbers based on the integral scale and Taylor microscale, which are defined as $Re_L=u_0L_{0}/\nu$ and $Re_\lambda =u_0\lambda /\nu$, respectively. When the linear forcing scheme is used, the integral scale is often evaluated as $L_0=u_0^3/\varepsilon _0$ with the root-mean-squared (r.m.s.) velocity fluctuations $u_0=\sqrt { ( \langle u^2\rangle +\langle v^2\rangle +\langle w^2\rangle )/3}$ because $L_0$ normalised by the domain size is determined as the property of the forcing scheme (Rosales & Meneveau 2005; Carroll & Blanquart 2013; Watanabe et al. 2021). On the other hand, the Taylor microscale is evaluated as $\lambda =\sqrt {15\nu u_0^{2}/\varepsilon _0}$. For the present DNS databases, $Re_\lambda$ ranges from 43 to 202.

Table 1. Parameters of DNS of statistically steady homogeneous isotropic turbulence: the number of grid points $N^3$, the spatial resolution $\varDelta /\eta$ and Reynolds numbers based on the integral scale and Taylor microscale ($Re_L$ and $Re_\lambda$).

2.2. Analysis of small-scale shear layers in turbulence

The small-scale shear layers in isotropic turbulence are analysed by detecting the shear layers with the triple decomposition of a velocity gradient tensor, as also reported in previous studies of shear layers in turbulent flows (Eisma et al. 2015; Watanabe et al. 2020; Fiscaletti et al. 2021). The DNS databases are analysed with the procedures briefly described below. Previous studies can be referred to for further details of the algorithms of the triple decomposition (Kolář 2007; Nagata et al. 2020; Fiscaletti et al. 2021) and the shear layer analysis (Watanabe & Nagata 2022) .

The shear layers are detected with the triple decomposition (Kolář 2007), which considers three local fluid motions: shear, rigid-body rotation and elongation (irrotational strain). The velocity gradient tensor $\boldsymbol {\nabla } {\boldsymbol {u}}$ is decomposed into the components associated with these motions as $\boldsymbol {\nabla } {\boldsymbol {u}}= \boldsymbol {\nabla } {\boldsymbol {u}}_S +\boldsymbol {\nabla } {\boldsymbol {u}}_R +\boldsymbol {\nabla } {\boldsymbol {u}}_E$, where the subscripts, $S$, $R$ and $E$, represent shear, rotation and elongation, respectively. The triple decomposition has to be applied in a so-called basic reference frame, which is different depending on position and time. The basic reference frame is the reference frame where the decomposition formula can most effectively extract the shear contribution from $\boldsymbol {\nabla } {\boldsymbol {u}}$ among all the reference frames. The basic reference frame can be specified from reference frames defined with three sequential rotational transformations ${\boldsymbol{\mathsf{Q}}}(\theta _{1},\theta _{2},\theta _{3})$ with $0^{\circ }\leq \theta _{1}\leq 180^{\circ }$, $0^{\circ }\leq \theta _{2}\leq 180^{\circ }$ and $0^{\circ }\leq \theta _{3}\leq 90^{\circ }$. The rotational transformation tensor ${\boldsymbol{\mathsf{Q}}}(\theta _1, \theta _2, \theta _3)$ is written as

(2.3) \begin{align} {\boldsymbol{\mathsf{Q}}} = \left( \begin{array}{@{}ccc@{}} \cos\theta_{1}\cos\theta_{2}\cos\theta_{3} - \sin\theta_{1}\sin\theta_{3} & \sin\theta_{1}\cos\theta_{2}\cos\theta_{3} + \cos\theta_{1}\sin\theta_{3} & -\sin\theta_{2}\cos\theta_{3}\\ -\cos\theta_{1}\cos\theta_{2}\sin\theta_{3} - \sin\theta_{1}\cos\theta_{3} & -\sin\theta_{1}\cos\theta_{2}\sin\theta_{3} + \cos\theta_{1}\cos\theta_{3} & \sin\theta_{2}\sin\theta_{3} \\ \cos\theta_{1}\sin\theta_{2} & \sin\theta_{1}\sin\theta_{2} & \cos\theta_{2} \end{array} \right). \end{align}

The velocity gradient tensor in a rotated reference frame is calculated as $(\boldsymbol {\nabla } \boldsymbol {u})^{*}={\boldsymbol{\mathsf{Q}}}(\boldsymbol {\nabla }\boldsymbol {u} ){\boldsymbol{\mathsf{Q}}}^{T}$, where $*$ represents a quantity evaluated in the rotated reference frame. The basic reference frame can be identified with a quantity called an interaction scalar, $I^{*}= |\textsf{$\mathit{\Omega}$} ^{*}_{12}{\mathsf{S}}^{*}_{12}| +|\textsf{$\mathit{\Omega}$} ^{*}_{23}{\mathsf{S}}^{*}_{23}| +|\textsf{$\mathit{\Omega}$} ^{*}_{31}{\mathsf{S}}^{*}_{31}|$. The basic reference frame assumes that $I^{*}$ has the largest value among all reference frames defined with ${\boldsymbol{\mathsf{Q}}}$ for $0^{\circ }\leq \theta _{1}\leq 180^{\circ }$, $0^{\circ }\leq \theta _{2}\leq 180^{\circ }$ and $0^{\circ }\leq \theta _{3}\leq 90^{\circ }$. By discretely changing the angles, $I^{*}$ is calculated for many reference frames. Then, $(\theta _{1},\theta _{2},\theta _{3})$ which gives the largest $I_S^*$ is used to define the basic reference frame. The angles are changed by $5^{\circ }$, for which the basic reference frame is accurately determined (Nagata et al. 2020; Fiscaletti et al. 2021). In the basic reference frame, the following decomposition is applied to extract $(\boldsymbol {\nabla } {\boldsymbol {u}}_{S})^*$:

(2.4)$$\begin{gather} (\boldsymbol{\nabla} {\boldsymbol{u}}_{RES})^{*}_{ij}= \mathrm{sgn}[(\boldsymbol{\nabla} {\boldsymbol{u}})^{*}_{ij}] \min[|(\boldsymbol{\nabla} {\boldsymbol{u}})^{*}_{ij}|,|(\boldsymbol{\nabla} {\boldsymbol{u}})^{*}_{ji}|], \end{gather}$$

(2.5)$$\begin{gather}(\boldsymbol{\nabla} {\boldsymbol{u}}_{S})^{*}_{ij}= (\boldsymbol{\nabla} {\boldsymbol{u}})^{*}_{ij}-(\boldsymbol{\nabla} {\boldsymbol{u}}_{RES})^{*}_{ij}, \end{gather}$$

for $i,j=1,2$ and 3 and sgn is a sign function. In addition, $(\boldsymbol {\nabla } {\boldsymbol {u}}_{RES})^{*}$ is further decomposed into the components of rotation and elongation as $(\boldsymbol {\nabla } {\boldsymbol {u}}_{R})^{*}_{ij}=[(\boldsymbol {\nabla } {\boldsymbol {u}}_{RES})^{*}_{ij}-(\boldsymbol {\nabla } {\boldsymbol {u}}_{RES})^{*}_{ji}]/2$ and $(\boldsymbol {\nabla } {\boldsymbol {u}}_{E})^{*}_{ij}=[(\boldsymbol {\nabla } {\boldsymbol {u}}_{RES})^{*}_{ij}+(\boldsymbol {\nabla } {\boldsymbol {u}}_{RES})^{*}_{ji}]/2$. Finally, $(\boldsymbol {\nabla } {\boldsymbol {u}}_{S})$, $(\boldsymbol {\nabla } {\boldsymbol {u}}_{R})$ and $(\boldsymbol {\nabla } {\boldsymbol {u}}_{E})$ in the original coordinate are obtained by applying the inverse transformation of ${\boldsymbol{\mathsf{Q}}}$ to $(\boldsymbol {\nabla } {\boldsymbol {u}}_{S})^{*}$, $(\boldsymbol {\nabla } {\boldsymbol {u}}_{R})^{*}$ and $(\boldsymbol {\nabla } {\boldsymbol {u}}_{E})^{*}$. The intensities of motions of shear and rotation are defined as $I_S=\sqrt {2(\boldsymbol {\nabla } {\boldsymbol {u}}_{S})_{ij}(\boldsymbol {\nabla } {\boldsymbol {u}}_{S})_{ij}}$ and $I_R=\sqrt {2(\boldsymbol {\nabla } {\boldsymbol {u}}_{R})_{ij}(\boldsymbol {\nabla } {\boldsymbol {u}}_{R})_{ij}}$. Figure 1 visualises the isosurfaces of $I_R$ and $I_S$ in NS2. These quantities are often used to detect vortices and shear layers in a three-dimensional flow. The shear layers identified with $I_S$ (white) have a flat shape whereas vortex tubes are detected with $I_R$ (orange).

Figure 1. Visualisation of shear layers and vortices, which are identified with the isosurfaces of $I_S=2\langle I_S\rangle$ (white) and $I_R=4\langle I_R\rangle$ (orange), respectively (NS2). Only a small part of the computational domain is shown here.

The small-scale shear layers in turbulence can be identified as flow regions with large shear intensity $I_S$. Therefore, following Watanabe et al. (2020), the present study examines a flow field around the local maxima of $I_S$. Here, the Hessian matrix of $I_S$, $\partial ^2 I_S/\partial x_i\partial x_j$, uniquely determines the locations of the local maxima without thresholds (Hayashi, Watanabe & Nagata 2021b). For each detected local maximum, a local shear coordinate $(\zeta _1, \zeta _2, \zeta _3)$ is introduced to observe the flow in a reference frame associated with shearing motion. The velocity vector in the shear coordinate is denoted by $(u_1, u_2, u_3)$. The shear coordinate assumes that the shear is expressed with $\partial u_3/\partial \zeta _2$. In this case, the shear contributes to the vorticity in the $\zeta _1$ direction. Therefore, the $\zeta _1$ direction is given by the unit vorticity vector of $\boldsymbol {\nabla } {\boldsymbol {u}}_{S}$, $\boldsymbol {\omega }_{S}/|\boldsymbol {\omega }_{S}|$, where $(\boldsymbol {\omega }_{S})_i=\epsilon _{ijk}(\boldsymbol {\nabla } {\boldsymbol {u}}_{S})_{jk}$ is a component of shear vorticity $\boldsymbol {\omega }_{S}$ and $\epsilon _{ijk}$ is the Levi-Civita symbol. Then, $\zeta _2$ and $\zeta _3$ are determined such that shearing motion is expressed with $\partial u_3/\partial \zeta _2$. The identification of $\zeta _2$ and $\zeta _3$ is based on a similar algorithm used to identify the basic reference frame. A large number of different reference frames are examined to identify the reference frame in which $(\boldsymbol {\nabla } {\boldsymbol {u}}_{S})_{32}$ becomes the largest among all the reference frames under the condition that the $\zeta _1$ direction is given by $\boldsymbol {\omega }_{S}/|\boldsymbol {\omega }_{S}|$. The numerical algorithm to find the shear coordinate is described in detail in Watanabe et al. (2020) and Hayashi et al. (2021a).

Once the shear coordinate is identified for a location of a local maximum of $I_S$, flow variables on the DNS grid, namely $f(x,y,z)$, are interpolated onto the shear coordinate with a third-order polynomial interpolation, providing $f(\zeta _1,\zeta _2,\zeta _3)$. Here, the shear coordinate is also discretised with a spacing smaller than the Kolmogorov scale. The interpolation on the shear coordinate is applied for all local maxima of $I_S$. Figure 2(a,b) visualise one of the detected shear layers in NS2 in the reference frame of DNS $(x,y,z)$ and in the shear coordinate $(\zeta _1,\zeta _2,\zeta _3)$. The shear intensity $I_S$ and the velocity vectors with respect to the velocity at the centre are shown on the two-dimensional planes. As the shear layers have no preference in the alignment with $(x,y,z)$, a shear flow is not always observed in the reference frame of DNS in figure 2(a). The shearing motion in the shear coordinate is represented by the flows in the $\pm \zeta _3$ directions. In figure 2(b), these parallel flows in opposite directions are observed on the sides of the thin shear layer. Therefore, the velocity in the $\zeta _3$ direction, $u_3$, significantly varies across the shear layer. The characteristics of each shear layer are assessed with the profiles of $u_3$ and $I_S$ along the $\zeta _2$ axis at $(\zeta _1,\zeta _3)=(0,0)$. Figure 2($c$) presents the plots of $u_3$ and $I_S$ across the shear layer visualised in figure 2($b$). The Kolmogorov length and velocity scales, $\eta$ and $u_\eta =(\nu \varepsilon _0)^{1/4}$, are used to normalise $\zeta _2$ and $u_3$, respectively, whereas $I_S$ is normalised by the shear intensity at $\zeta _2=0$, $I_{S,max}$. As the shear layer is thin in the $\zeta _2$ direction, $I_{S}(\zeta _2)$ decreases rapidly with $\zeta _2$ from the peak at $\zeta _2=0$. The layer thickness $\delta _S$ is evaluated as the half-width of $I_S$ on the $\zeta _2$ axis, namely the distance between two locations with $I_{S}(\zeta _2)/I_{S,max}=0.5$ for $\zeta _2>0$ and $\zeta _2<0$. In addition, the velocity jump across the shear layer, $u_S$, is calculated with $\partial u_3/\partial \zeta _2$ at $\zeta _2=0$ as $u_S=(\partial u_3/\partial \zeta _2)\delta _S$. These length and velocity scales are also shown in figure 2(c). Then, the shear Reynolds number is defined as $Re_S=u_S\delta _S/\nu$. These definitions of $u_S$, $\delta _S$ and $Re_S$ are applied for an instantaneous flow field around each shear layer whereas the length and velocity scales estimated from a mean flow field around the shear layers have been discussed in previous studies (Watanabe et al. 2020; Fiscaletti et al. 2021; Hayashi et al. 2021a).

Figure 2. A local flow field around a shear layer in ($a$) the coordinate used in DNS and ($b$) the shear coordinate. The local maximum of $I_S$ is located at the centre of each figure. A colour contour of shear intensity $I_S$ normalised by the Kolmogorov time scale $\tau _{\eta }=(\nu /\varepsilon _0)^{1/2}$ and two-dimensional velocity vectors relative to the fluid motion at $(y,z)=(0,0)$ or $(\zeta _2, \zeta _3)=(0,0)$ are shown on the two-dimensional planes. The length of the arrows represents the vector magnitude. ($c$) Variations of $I_S$ and velocity in the $\zeta _3$ direction, $u_3$, across the shear layer visualised in ($b$) along the $\zeta _2$ axis at $(\zeta _1, \zeta _3)=(0,0)$.

The mean velocity profile around shear layers has shown that the shear layer appears in a biaxial strain field. The enstrophy production $P_\omega = \omega_i {\mathsf{S}}_{ij} \omega_j$ actively occurs in the shear layers because of the interaction between the shear and biaxial strain (Buxton & Ganapathisubramani 2010; Elsinga & Marusic 2010; Eisma et al. 2015; Elsinga et al. 2017; Watanabe et al. 2020). However, some shear layers may have negative $P_\omega$. The number of the identified local maxima of $I_S$ is evaluated separately for $P_\omega >0$ and $P_\omega <0$, for which the numbers are denoted by $N_{P_\omega >0}$ and $N_{P_\omega <0}$, respectively. The fractions of shear layers with $P_\omega >0$ and $P_\omega <0$ are given by $F_{P_\omega >0}=N_{P_\omega >0}/N_{SH}$ and $F_{P_\omega <0}=N_{P_\omega <0}/N_{SH}$ with $N_{SH}=N_{P_\omega >0}+N_{P_\omega <0}$. Figure 3 plots $F_{P_\omega >0}$ and $F_{P_\omega <0}$ as functions of $Re_\lambda$. For all cases, $F_{P_\omega >0}\approx 0.94$ and $F_{P_\omega <0}\approx 0.06$ are obtained regardless of the Reynolds number. Almost all shear layers have positive enstrophy production. The primary interest of this paper is in the characteristics of strained shear layers with $P_\omega >0$ for comparisons with Burgers’ vortex layer and numerical simulations of an isolated shear layer in a biaxial strain field, both of which also have $P_\omega >0$. Therefore, the statistics of $u_S$, $\delta _S$ and $Re_S$ are calculated with ensemble averages of the shear layers with $P_\omega >0$.

Figure 3. Fractions of shear layers with positive and negative enstrophy production $P_\omega =\omega _i {\mathsf{S}}_{ij}\omega _j$, $F_{P_\omega >0}$ and $F_{P_\omega <0}$, as functions of the turbulent Reynolds number $Re_\lambda$.

2.3. Thickness, velocity jump and Reynolds number of small-scale shear layers

Variations of the shear layer characteristics are examined with the probability density functions (p.d.f.s) of $\delta _S$, $u_S$ and $Re_S$. Figure 4 shows the p.d.f.s of $\delta _S$ and $u_S$, which are normalised by the Kolmogorov length and velocity scales ($\eta$ and $u_\eta$), respectively. The p.d.f.s hardly depend on $Re_\lambda$. Because of $u_\eta \eta /\nu =1$, $Re_S$ can be written as $Re_S=(u_S/u_\eta )(\delta _S/\eta )$. Therefore, the p.d.f.s of $\delta _S/\eta$ and $u_S/u_\eta$ suggest that $Re_S$ hardly depends on $Re_\lambda$ and the behaviour of typical small-scale shear layers in turbulence can be insensitive to the flow Reynolds number. However, the Reynolds number dependence can be important for intense shear layers because an extremely large velocity gradient tends to be observed at a very high Reynolds number (Elsinga et al. 2017; Buaria et al. 2019; Das & Girimaji 2020; Ghira et al. 2022).

Figure 4. Probability density functions of ($a$) thickness $\delta _S$ and ($b$) velocity jump $u_S$ of shear layers. The thickness and velocity jump are normalised by the Kolmogorov length and velocity scales, $\eta$ and $u_\eta$, respectively.

The $Re_S$ dependence of $\delta _S$ and $u_S$ is examined with the joint p.d.f.s of $(\delta _S/\eta, Re_S)$ and $(u_S/u_\eta, Re_S)$ in figure 5(a,b), respectively. Here, the results for NS2 are presented because the $Re_\lambda$ dependence is weak when the Kolmogorov scales are used for normalisation. The Reynolds number $Re_S$ is positively correlated with $\delta _S$ and $u_S$, as expected from the definition of $Re_S=u_S \delta _S/\nu$. The highest probability is observed for $(\delta _S/\eta,Re_S)\approx (4,20)$ and $(u_S/u_\eta,Re_S)\approx (5,20)$, suggesting that the shear layers in turbulence typically have $(\delta _S/\eta,u_S/\eta, Re_S)\approx (4,5,20)$ although variations from these values are confirmed from the distributions of the joint p.d.f.s. Specifically, the shear Reynolds number $Re_S$ ranges approximately between 10 and 70.

Figure 5. Joint probability density functions of ($a$) $(\delta _S/\eta, Re_S)$, ($b$) $(u_S/u_{\eta }, Re_S)$ and ($c$) $(\delta _S/\delta _{BV}, Re_S)$ of shear layers in NS2. Here, the shear layer thickness $\delta _S$ is normalised with the Kolmogorov scale $\eta$ or the thickness of Burgers’ vortex layer $\delta _{BV}=1.67\sqrt {2\nu /\alpha }$ calculated with the strain rate $\alpha =\omega _i {\mathsf{S}}_{ij}\omega _j/\omega _k\omega _k$ while the velocity jump $u_S$ is normalised by the Kolmogorov velocity scale $u_\eta$.

The relation between the shear layer and the nearby strain field in turbulence agrees with Burgers’ vortex layer (Watanabe et al. 2020), which is one of the exact steady solutions of Navier–Stokes equations for a uniform shear layer subject to a biaxial strain (Burgers 1948; Davidson 2004). In both shear layers in turbulence and Burgers’ vortex layer, a compressive strain acts in the normal direction of the shear layer while an extensive strain acts in the vorticity direction of the shear. For Burgers’ vortex layer, the half-width of $I_S$ is related to the strain rate $\alpha$ as $\delta _{BV}=1.67\sqrt {2\nu /\alpha }$. The strain rate acting on vortices is often estimated as $\alpha =\omega _i {\mathsf{S}}_{ij}\omega _j/\omega _k\omega _k$, for which the diameter of Burgers’ vortex has been compared with the actual diameter of vortex tubes (Jiménez & Wray 1998; da Silva et al. 2011; Jahanbakhshi et al. 2015; Watanabe et al. 2017b; Ghira et al. 2022). Similarly, $\delta _{BV}$ is evaluated with $\alpha$ at $(\zeta _1, \zeta _2, \zeta _3)=(0,0,0)$ for each shear layer and is compared with the half-width of $I_S$, $\delta _S$. Here, the $Re_S$ dependence of the relation between $\delta _S$ and $\delta _{BV}$ is examined with the present DNS database, whereas a general comparison between $\delta _S$ and $\delta _{BV}$ has been provided for isotropic turbulence and planar jets in Watanabe et al. (2020) and Hayashi et al. (2021a). Figure 5(c) shows the joint p.d.f. between $\delta _S/\delta _{BV}$ and $Re_S$. For the typical shear Reynolds number $Re_S\approx 20$, the joint p.d.f. reaches the maximum for $\delta _S/\delta _{BV}\approx 1$, implying that the relation between the layer thickness and the strain rate is consistent with Burgers’ vortex layer. The p.d.f. also shows that $\delta _S/\delta _{BV}$ is positively correlated with $Re_S$. Thus, the layer thickness tends to be larger than that of Burgers’ vortex layer at large $Re_S$, and vice versa. An unsteady analysis of Burgers’ vortex layer predicts that the layer thickness asymptotically approaches $\delta _{BV}$ with time when the initial thickness differs from $\delta _{BV}$ (Davidson 2004). Therefore, the shear layers with large $Re_S$ are possibly being flattened with time due to the compressible strain while those with small $Re_S$ are thickened via the viscous effect. In addition, the p.d.f. is distributed for $0.5\lesssim \delta _S/\delta _{BV}\lesssim 1.5$ and most shear layers in isotropic turbulence are close to an equilibrium state in terms of the relation between the thickness and the strain rate.

These probability distributions of length and velocity scales are consistent with previous studies of a mean flow in a reference frame defined with the eigenvectors of the rate-of-strain tensor (Elsinga & Marusic 2010; Elsinga et al. 2017). This reference frame is often referred to as a strain eigenframe. These studies introduce the strain eigenframe at each point of the flow and take averages of various quantities as functions of the position in the strain eigenframe. This approach is similar to the present method, which relies on the local reference frame defined with $\boldsymbol {\nabla } {\boldsymbol {u}}_{S}$. Even though the averages in the strain eigenframe are not taken for shear layers, a mean flow profile in the strain eigenframe also exhibits a shear-layer pattern (Elsinga & Marusic 2010), which is similar to the mean flow around the shear layers (Watanabe & Nagata 2022). The p.d.f. of shear intensity $I_S$ in isotropic turbulence and turbulent jets indicates that shearing motions appear anywhere in turbulence (Nagata et al. 2020; Hayashi et al. 2021a). For this reason, the statistics of turbulence can be strongly influenced by shearing motions, and therefore, the averages in the strain eigenframe taken for an entire turbulent region is consistent with those taken solely for shear layers. The present study examines the length and velocity scales of each shear layer with their probability distribution and confirms that the most typical values of the half-width and velocity jump of the shear layers are $\delta _S/\eta \approx 4$ and $u_S/u_\eta \approx 5$, respectively. Similarly, the shear layer of the mean flow in the strain eigenframe has a full width of about $10\eta$ and a velocity jump of about $6u_\eta$. The p.d.f. in figure 4 confirms that some shear layers have different velocity and length scales from the typical values. However, the mean flow in the strain eigenframe is strongly dominated by the most typical shear layers with $\delta _S/\eta \approx 4$ and $u_S/u_\eta \approx 5$.

3. Numerical simulations of an isolated shear layer in a biaxial strain field

3.1. The model of an isolated shear layer

Numerical simulations are performed for an isolated shear layer in a biaxial strain field illustrated in figure 6. This simple model is an approximation of the mean-flow pattern observed around the small-scale shear layers in turbulence (Watanabe et al. 2020; Hayashi et al. 2021a; Watanabe & Nagata 2022), which also resembles the shear-layer pattern observed in a mean flow with a strain eigenframe taken regardless of the presence of shear layers (Elsinga & Marusic 2010; Elsinga et al. 2017). The model assumes that the aspect ratio of the length scales in the $x$ and $y$ directions is finite. The simulations are carried out to investigate the response of the shear layer to perturbations. Some physical parameters are determined based on the shear layer characteristics in isotropic turbulence presented in § 2 to relate the behaviour of the isolated shear layer to that of small-scale shear layers in turbulence. The shear coordinates, $\zeta _1$, $\zeta _2$ and $\zeta _3$, used in § 2 are respectively related to $z$, $y$ and $x$ in figure 6. The model follows the numerical study of a row of alternating vortex sheets in a biaxial strain field (Lin & Corcos 1984) although the present study considers a single isolated shear layer. In this section, a dimensional variable is denoted with a symbol with a tilde, e.g. $\tilde {f}$, whose non-dimensional counterpart is $f$. The model assumes an incompressible fluid. The transverse and vorticity directions of the shear layer are denoted by $\tilde {y}$ and $\tilde {z}$, respectively. The shear layer is initially parallel to the $\tilde {x}$ direction. The shear layer is subject to a biaxial strain with compression in the $\tilde {y}$ direction and stretching in the $\tilde {z}$ direction. The velocity components in the $\tilde {x}$, $\tilde {y}$ and $\tilde {z}$ directions are denoted by $\tilde {u}$, $\tilde {v}$ and $\tilde {w}$, respectively. The origin of the coordinate system is located at the centre of the shear layer.

Figure 6. An isolated shear layer with a finite aspect ratio in a biaxial strain field.

The initial velocity field $\tilde {{\boldsymbol {u}}}=(\tilde {u},\tilde {v},\tilde {w})$ can be decomposed into the components of a base flow and perturbations as $\tilde {\boldsymbol {u}}=\tilde {\boldsymbol {u}}_B+\tilde {\boldsymbol {u}}'$. In addition, the base flow, $\tilde {\boldsymbol {u}}_B$, is further decomposed into the components associated with the strain and vorticity as $\tilde {\boldsymbol {u}}_B=\tilde {\boldsymbol {u}}_S+\tilde {\boldsymbol {u}}_V$, the latter of which can be defined with the vorticity profile. The present model assumes that an initial profile of vorticity is given by

(3.1a–c)\begin{align} \tilde{\omega}_x(\tilde{x},\tilde{y},\tilde{z})=0; \quad \tilde{\omega}_y(\tilde{x},\tilde{y},\tilde{z})=0; \quad \tilde{\omega}_z(\tilde{x},\tilde{y},\tilde{z}) = \omega_0 \exp\left(-\frac{\tilde{y}^2}{\delta^2}\right) \exp\left[-\frac{\tilde{x}^2}{(A_R \delta)^2}\right]. \end{align}

Here, $\omega _0$ is the peak vorticity at the centre of the shear layer, $A_R$ is the aspect ratio and $\delta$ is the shear layer thickness. The velocity components for the time-independent biaxial strain, $\tilde {\boldsymbol {u}}_S=(\tilde {u}_S, \tilde {v}_S, \tilde {w}_S)$, are given by

(3.2a–c)\begin{equation} \tilde{u}_S(\tilde{x},\tilde{y},\tilde{z})=0, \quad \tilde{v}_S(\tilde{x},\tilde{y},\tilde{z})=-\gamma \tilde{y}, \quad \tilde{w}_S(\tilde{x},\tilde{y},\tilde{z})=\gamma \tilde{z} \end{equation}

with a constant strain rate $\gamma$, which is assumed to satisfy the equilibrium relation $\gamma =2\nu /\delta ^2$ of Burgers’ vortex layer. The DNS results in § 2 have suggested that the majority of shear layers in turbulence approximately satisfy this relation. For the vorticity profile of (3.1a–c), $w_V=0$ is satisfied anywhere. The remaining components $(\tilde {u}_V,\tilde {v}_V)$ are determined as $(\tilde {u}_V,\tilde {v}_V)=(-\partial \tilde {\phi }/\partial \tilde {y}, \partial \tilde {\phi }/\partial \tilde {x})$ with a vector potential $\tilde {\phi }$ which satisfies $\tilde {\omega }_z=\nabla ^2 \tilde {\phi }$. For $A_R\rightarrow \infty$, (3.1a–c) and (3.2a–c) express Burgers’ vortex layer, which is infinitely long in the $x$ direction. A similar vorticity profile was also used in Lin & Corcos (1984), where a cosine function was used instead of $\exp [-\tilde {x}^2/(A_R \delta )^2]$ to model the row of shear layers with positive and negative spanwise vorticities. Their model allows the interaction between adjacent shear layers, which is observed in the development of a plane mixing layer, whereas this interaction does not occur in the present model of the isolated shear layer because the simulations are designed to investigate the development of a single shear layer under perturbations.

The perturbations $\tilde {\boldsymbol {u}}'$ are defined based on the spanwise-vorticity perturbations $\tilde {\omega }_z'$, whose profile is given by

(3.3)$$\begin{gather} \tilde{\omega}_z'(\tilde{x},\tilde{y},\tilde{z}) = \omega'_0 \mathrm{sin}\left[2{\rm \pi} \left(\frac{\tilde{x}}{\varLambda}+\varphi_x\right)\right] \mathrm{sin}\left[2{\rm \pi} \left(\frac{\tilde{y}}{\varLambda}+\varphi_y\right)\right] G(\tilde{x}, \tilde{y}), \end{gather}$$

(3.4)$$\begin{gather}G(\tilde{x}, \tilde{y})= \left[ \frac{1}{2} + \frac{1}{2}\mathrm{tanh}\left( \frac{L_{xf}-2|\tilde{x}|}{12\delta} \right) \right] \left[ \frac{1}{2} + \frac{1}{2}\mathrm{tanh}\left( \frac{L_{yf}-2|\tilde{y}|}{12\delta} \right) \right], \end{gather}$$

where $\varLambda$ and $\omega '_0$ are the wavelength and amplitude, respectively, whereas $\varphi _x$ and $\varphi _y$ are the phases of the perturbations. The top-hat function $G(\tilde {x}, \tilde {y})$ is equal to 1 near the shear layer and becomes 0 in the far field. Therefore, the sinusoidal perturbations exist only for $|\tilde {x}|\lesssim L_{xf}/2$ and $|\tilde {y}|\lesssim L_{yf}/2$. Once the profile of the vorticity perturbations is determined, $(\tilde {u}',\tilde {v}',\tilde {w}')=(-\partial \tilde {\phi }'/\partial \tilde {y}, \partial \tilde {\phi }'/\partial \tilde {x}, 0)$ is calculated with a vector potential $\tilde {\phi }'$ which satisfies $\tilde {\omega }_z'=\nabla ^2 \tilde {\phi }'$. The profile of $(\tilde {u}',\tilde {v}')$ exhibits fluctuations around $(\tilde {u}',\tilde {v}')=(0,0)$, where the maximum amplitude is denoted by $u'_0$. Because $\omega '_0$ is directly related to $u'_0$, the numerical parameter is determined by specifying $u'_0$. This can be done by generating $(\tilde {u}',\tilde {v}')$ from (3.3) and (3.4) with given values of $\varLambda$ and $\omega '_0$. Then, $(\tilde {u}',\tilde {v}')$ are normalised such that the maximum value of $\tilde {u}'$ is equal to $u'_0$.

Numerical simulations are performed with governing equations non-dimensionalised with reference velocity and length scales, $u_r$ and $l_r$, e.g. $u=\tilde {u}/u_r$ and ${\omega }_z=\tilde {\omega }_z/(u_r/l_r)$. Time $\tilde {t}$ can also be non-dimensionalised as $t=\tilde {t}/t_r$ with the reference time scale $t_r=l_r/u_r$. Here, $u_r$ and $l_r$ are given by $u_r=\sqrt {{\rm \pi} }\omega _0\delta$ and $l_r=\delta$. The choice of $u_r$ is based on Burgers’ vortex layer, where the velocity jump across the layer is given by $\sqrt {{\rm \pi} }\omega _0\delta$.

3.2. Governing equations, numerical methods and parameters

The present model assumes that the initial profiles of $u$ and $v$ are independent of $z$. In this case, $u$ and $v$ remain independent of $z$ even after the shear layer development (Lin & Corcos 1984). Therefore, the continuity and momentum equations non-dimensionalised by $u_r$ and $l_r$ are written as

(3.5)$$\begin{gather} \frac{\partial u}{\partial x} + \frac{\partial v}{\partial y}=-\frac{2}{ Re}, \end{gather}$$

(3.6)$$\begin{gather}\frac{\partial u}{\partial t}+u\frac{\partial u}{\partial x} +v\frac{\partial u}{\partial y} =-\frac{\partial p}{\partial x} +\frac{1}{Re}\left( \frac{\partial^2 u}{\partial x^2} +\frac{\partial^2 u}{\partial y^2} \right), \end{gather}$$

(3.7)$$\begin{gather}\frac{\partial v}{\partial t}+u\frac{\partial v}{\partial x} +v\frac{\partial v}{\partial y} =-\frac{\partial p}{\partial y} +\frac{1}{Re}\left( \frac{\partial^2 v}{\partial x^2} +\frac{\partial^2 v}{\partial y^2} \right), \end{gather}$$

(3.8)$$\begin{gather}\frac{\partial w}{\partial t}=0, \end{gather}$$

with the Reynolds number $Re=u_rl_r/\nu =\sqrt {{\rm \pi} }\omega _0\delta ^2/\nu$. Here, $2/Re$ on the right-hand side of (3.5) is obtained from $\partial w/\partial z= \gamma /(u_r/l_r)$ with $\gamma =2\nu /\delta ^2$.

Solving these equations yields the temporal evolution of the shear layer initialised with the velocity profile explained above. As $u$ and $v$ are independent of $z$ and $w$ remains unchanged with time, the simulations can be conducted in a two-dimensional computational domain. The governing equations are solved with an in-house finite-difference code, which is developed by reducing the dimension of our DNS code from three to two. The original DNS code was used for various turbulent flows in our previous studies, where the code was validated by comparisons with experiments and other numerical simulations (e.g. Watanabe & Nagata 2018; Watanabe, Zhang & Nagata 2019; Katagiri, Watanabe & Nagata 2021; Watanabe, Zheng & Nagata 2022). The code is the same as that used for isotropic turbulence except for spatial discretisation, for which a second-order fully conservative finite-difference scheme is employed (Morinishi et al. 1998).

The isolated shear layer is simulated in a rectangular domain with a size of $(L_x,L_y)=(250\delta, 250\delta )$, which is much larger than the shear layer. Therefore, the boundary conditions assume the far-field condition of the base flow. Here, $u=0$ is applied at all boundaries. The zero-gradient condition is applied to $v$ at the boundaries in the $x$ direction while $\tilde {v}=-\gamma \tilde {y}$ is used for the $y$ direction. The computational domain is discretised with $(N_x\times N_y)=(300\times 300)$ grid points. A non-uniform grid is employed to achieve a high spatial resolution near the shear layer. Here, $x$ and $y$ positions of the grid are determined with integers $i$ and $j$ by the following mapping functions:

(3.9)$$\begin{gather} \tilde{x}(i)=-\frac{L_{x}}{2\alpha_x} {\mathrm{atanh}} \left[\left( {\mathrm{tanh}}\,\alpha_{x}\right) \left(1- \frac{2i}{N_x} \right) \right]\quad \mathrm{with}\ i=0,\ldots,N_{x}, \end{gather}$$

(3.10)$$\begin{gather}\tilde{y}(\,j)=-\frac{L_{y}}{2\alpha_y}{\mathrm{atanh}}\left[\left({\mathrm{tanh}}\,\alpha_{y}\right) \left(1-\frac{2j}{N_y} \right) \right]\quad \mathrm{with}\ j=0,\ldots,N_{y}, \end{gather}$$

with $\alpha _x=\alpha _y=3.5$. The grid size becomes large as being away from the centre of the domain. The distribution of the grid size for these functions can be found in Watanabe et al. (2018). A non-dimensional time increment is $\Delta t=0.025$, and time is advanced until $t=50$ in each simulation.

The numerical simulations are performed for various sets of the Reynolds number $Re$, the aspect ratio $A_R$ and the perturbation wavelength $\varLambda$. Table 2 summarises $Re$ and $A_R$ considered in this study. Here, these simulations consider $\varphi _x=\varphi _y=0$ while the influence of the phases is investigated with additional simulations, as explained below. The Reynolds number is determined based on the DNS data while $A_R=10$ and 30 are chosen to discuss two different types of the evolution of the shear layer. In turbulent flows, $A_R$ differs for each shear layer as confirmed from visualisations of shear layers in figure 1 and previous studies (Nagata et al. 2020; Fiscaletti et al. 2021; Hayashi et al. 2021b). We have noticed that it is difficult to quantify the aspect ratio which can be directly compared with $A_R$ of the ideal numerical model because of the curved shape of shear layers in turbulence. In figure 2($b$), the shear layer extends approximately from $\zeta _3/\eta =-15$ to 30 whereas the thickness is about $4\eta$. Thus, the aspect ratio is close to 10. Visualisation of shear intensity $I_S$ or shear vorticity $|\boldsymbol {\omega }_{S}|$ in isotropic turbulence and turbulent free shear flows has suggested that the aspect ratio of some shear layers is close to or smaller than about $10$ (Fiscaletti et al. 2021; Hayashi et al. 2021a; Watanabe & Nagata 2022). However, these studies have also observed shear layers with a larger aspect ratio than 10. A long shear layer is broken during the roll-up and generates more than one shear layer with small $A_R$ around a vortex tube (Horiuti & Takagi 2005; Watanabe et al. 2020). This might be the reason why the shear layers with a wide range of $A_R$ are observed in turbulence. Therefore, instead of matching the aspect ratio, two cases of $A_R=10$ and 30 are considered because the shear layer behaves differently for large and small $A_R$. In the numerical simulations of the isolated shear layer, $Re$ is defined with the layer thickness $\delta$ in (3.1a–c) and the velocity difference $u_r=\sqrt {{\rm \pi} }\omega _0\delta$ across the layer. On the other hand, the DNS results of isotropic turbulence have been examined with the shear Reynolds number $Re_S$ defined with the half-width of the shear intensity, $\delta _S$, and the velocity jump $u_S$ defined with $\delta _S$ and the local velocity gradient in the shear layer. Table 2 also lists $\delta _S/\delta$, $u_S/u_r$ and the shear Reynolds number $Re_S$ evaluated with the same definitions of $\delta _S$ and $u_S$ as in the DNS of isotropic turbulence. The present simulations consider $20\lesssim Re_S\lesssim 70$, which is determined from the p.d.f. of $Re_S$ in the isotropic turbulence. For each set of $(Re, A_R)$, the simulations are performed for the perturbation wavelength $\varLambda /\delta =8, 12,\ldots$ or 40 to examine the wavelength dependence. The amplitude of velocity fluctuations is fixed in all simulations as $u'_0/u_r=0.1$, for which the results are compared with unperturbed cases $(u'_0=0)$. For $G(x,y)$ in (3.4), $L_{xf}=125\delta$ and $L_{yf}=100\delta$ are used for the perturbations to cover a large area around the shear layer.

Table 2. Parameters of numerical simulations of an isolated shear layer in a biaxial strain field: the Reynolds number $Re=u_rl_r/\nu$ defined with reference scales of the simulations, the aspect ratio of the shear layer $A_R$, the ratio between the half-width of $I_S$ and the layer thickness $\delta _S/\delta$, the velocity scale ratio $u_S/u_r$ and the shear Reynolds number $Re_S=u_S\delta _S/\nu$.

A uniform shear layer with $A_R=\infty$ is also investigated for a comparison with the finite $A_R$ cases. The initial vorticity profile is (3.1a–c) with $A_R\rightarrow \infty$. The periodic boundary condition is applied in the $x$ direction. The region of the perturbations is given by (3.4) with $L_{xf}\rightarrow \infty$. The uniform shear layer is simulated for $Re=30$ with perturbations with $u'_0/u_r=0.1$ and $\varLambda /\delta =8, 12,\ldots$ or 40. As the perturbations should also be periodic in the $x$ direction, the domain size is determined with $\varLambda$ as $(L_x,L_y)=(4\varLambda, 250\delta )$. The grid spacing in the $x$ direction is uniform, and $N_x$ is determined such that the spacing is $0.2\delta$. Other parameters and conditions are the same for both infinite- and finite-$A_R$ cases.

3.3. Temporal development of an isolated shear layer

Figure 7 shows the profiles of vorticity $\omega _z$ to visualise the temporal evolution of a shear layer with a low aspect ratio ($A_R=10$). The unperturbed case with $u_0'=0$ is shown in figure 7($a$) while the perturbed cases with the wavelength of $\varLambda /\delta =16$ and 32 are shown in figure 7(b,c). In all cases, the initially flat shear layer is distorted with time and forms a circular region with large vorticity. The formation of the vortex occurs even in the unperturbed case. This is caused by the self-induced velocity of the shear layer with the low aspect ratio. The sheet-like distribution of vorticity accompanies the upward and downward velocities on the right and left of the shear layer, respectively (Lin & Corcos 1984). Consequently, the shear layer is distorted by itself and generates the vortex centred at $(x,y)=(0,0)$. A similar formation of the vortex is also observed in the perturbed cases in figure 7(b,c). However, the time required for the vortex formation is different among the three cases. At $t=40$, the value of $\omega _z$ in the vortex for the unperturbed case is smaller and larger than those for the perturbed cases with $\varLambda /\delta =16$ and $32$, respectively. In addition, as the vortex forms at $t=30$ and 40, the shape of the large-$\omega _z$ region changes from an ellipse to a circular shape. Compared with the unperturbed case, the vortex formation is accelerated for the perturbations with $\varLambda /\delta =16$ whereas it is delayed for $\varLambda /\delta =32$.

Figure 7. Temporal evolutions of an isolated shear layer in a biaxial strain field with the aspect ratio $A_R=10$ and the Reynolds number $Re=30$: ($a$) unperturbed case ($u_0'=0$); perturbed cases with ($b$) $u_0'/u_r=0.1$ and $\varLambda /\delta =16$ and ($c$) $u_0'/u_r=0.1$ and $\varLambda /\delta =32$. From top to bottom, time increases as $t=0$, 5, 10, 20, 30 and 40. Only a small part of the computational domain is shown here.

Figure 8 shows the temporal evolution of a shear layer with a high aspect ratio ($A_R=30$) without perturbations whereas the perturbed cases with $\varLambda /\delta =16$ and 32 are visualised in figures 9 and 10, respectively. Unlike the shear layer with small $A_R$, the vortex formation does not occur for the unperturbed case with large $A_R$ in figure 8. The vortices form in the shear layer with perturbations with $\varLambda /\delta =16$ in figure 9, where three vortices can be found at $t=40$. However, the perturbations with a twice-as-large wavelength, $\varLambda /\delta =32$, are not effective to trigger the formation of vortices, as shown in figure 10. The vortex formation for $\varLambda /\delta =16$ is initiated with the distortion of the shear layer by the perturbations, and its process is more similar to Kelvin–Helmholtz instability (Patnaik et al. 1976) than that for the low $A_R$ case in figure 7.

Figure 8. Temporal evolution of an unperturbed isolated shear layer in a biaxial strain field with the aspect ratio $A_R=30$, the Reynolds number $Re=30$ and $u_0'/u_r=0$. From top to bottom, time increases as $t=0$, 5, 10, 20, 30 and 40. Only a small part of the computational domain is shown here.

Figure 9. The same as figure 8 but for a perturbed case with $u_0'/u_r=0.1$ and $\varLambda /\delta =16$.

Figure 10. The same as figure 8 but for a perturbed case with $u_0'/u_r=0.1$ and $\varLambda /\delta =32$.

When the vortices are successfully formed in the shear layer, the maximum vorticity $\omega _{max}$ in the shear layer increases with time. Figure 11 shows the temporal variation of $\omega _{max}$ of the shear layers with $Re=30$ for different $\varLambda$ and $A_R$. The results are not shown for all values of $\varLambda$ for visibility of the figure while $\varLambda$ dependence is examined in detail at several time instances below. For the low aspect ratio with $A_R=10$ in figure 11(a), $\omega _{max}$ increases for both unperturbed and perturbed cases with all wavelengths because of the vortex formation. The growth of $\omega _{max}$ occurs rapidly by the perturbations for $A_R/\delta =12$ and 16 whereas it is delayed for larger wavelengths. For the high aspect ratio in figure 11(b), $\omega _{max}$ without perturbations does not vary until $t\approx 30$ and then begins to decrease. The increase of $\omega _{max}$ can be seen for the perturbed case with $\varLambda /\delta =12$, 16 and 24. The perturbations with larger $\varLambda$ do not cause a continuous increase of $\omega _{max}$ and result in the decay at a late time for $\varLambda /\delta =32$ and 40. For both aspect ratios, the perturbations with $\varLambda /\delta =12$ and 16 are effective to promote the vortex formation, which does not occur or is delayed for larger $\varLambda$.

Figure 11. Temporal variations of the maximum vorticity $\omega _{max}$ of the isolated shear layer with $Re=30$ and ($a$) $A_R=10$ and ($b$) $A_R=30$.

The above results are obtained for $Re=30$ whereas the Reynolds number dependence is examined in figures 12 and 13, where the vorticity ratio $\omega _r$ between $\omega _{max}(t)$ and its initial value is plotted as a function of $\varLambda /\delta$ for three Reynolds numbers at $t=20$, 30 and 50. When the vortex formation occurs in the shear layer and the spanwise vorticity is amplified, $\omega _r$ increases with time. On the other hand, the increase of $\omega _r$ does not occur when the shear layer is simply distorted without vortex formation. The results for $A_R=10$ at $t=20, 30$ and 50 are shown in figures 12($a$), 12($b$) and 12($c$), respectively. For comparison, $\omega _r$ for the unperturbed case is plotted at $\varLambda =0$. As the vortex formation occurs without perturbations for $A_R=10$, $\omega _r$ gradually increases with time even without perturbations. However, the rapid growth of $\omega _r$ is observed for the perturbations with $\varLambda /\delta \approx 12$, for which $\omega _r$ tends to be large compared with other values of $\varLambda /\delta$. At each time instance, $\omega _r$ decreases from the peak of $\omega _r$ at $\varLambda /\delta =8$–$16$ as $\varLambda /\delta$ increases. In addition, $\omega _r$ for $\varLambda /\delta =8$ is smaller than those for $\varLambda /\delta =12$ at $t=30$ and 50 in figure 12(b,c) except for $Re=45$ at $t=30$. Therefore, the perturbations with $\varLambda /\delta \approx 12$ more effectively cause the amplification of vorticity than those with larger or smaller $\varLambda /\delta$. This wavelength dependence is similar for all $Re$. Figure 13 shows the results for $A_R=30$, for which perturbations are necessary for the vortex formation. The results for the uniform shear layer ($A_R\rightarrow \infty$) with $Re=30$ are also shown for comparison. For $Re=30$, the $\varLambda$ dependence of $\omega _r$ is similar for the isolated shear layer with $A_R=30$ and the uniform shear layer. Therefore, the effects of the finite aspect ratio are not significant for $A_R=30$, except for a late time. The wavelength dependence of the growth of $\omega _r$ for $A_R=30$ is stronger than that for $A_R=10$, especially at $t=50$. As also found for $A_R=10$, $\omega _r$ is large for $\varLambda /\delta \approx 12$, and the perturbations with larger or smaller $\varLambda /\delta$ do not cause the rapid increase of $\omega _r$. A different trend is observed for $Re=15$, for which $\omega _r$ decreases with time for the unperturbed case shown at $\varLambda /\delta =0$. In figure 11(b), $\omega _{max}$ for $u'_0=0$ does not vary with time until $t=30$ and decreases for a later time. This decrease is also observed for $Re=15$, resulting in $\omega _r<1$ in figure 13. However, we note that the vortex formation with the amplification of $\omega _z$ occurs for the perturbed cases with $12\leq \varLambda /\delta \leq 20$ and $Re=15$. Therefore, $\omega _r$ for $12\leq \varLambda /\delta \leq 20$ is larger than those for the unperturbed case and for larger or smaller $\varLambda /\delta$. In summary, the perturbations with $\varLambda /\delta \approx 12$ tend to amplify the spanwise vorticity for all $Re$ and $A_R$ considered here. This wavelength of $\varLambda /\delta =12$ corresponds to $\varLambda /\delta _S\approx 7.2$, where $\delta _S$ is the half-width of shear intensity used in the shear layer analysis for isotropic turbulence and is related to $\delta$ as presented in table 2.

Figure 12. Wavenumber dependence of maximum vorticity normalised by its initial value, $\omega _r=\omega _{max}(t)/\omega _{max}(t=0)$ for the aspect ratio $A_R=10$. The results taken at ($a$) $t=20$, ($b$) $t=30$ and ($c$) $t=50$ are plotted as functions of the perturbation wavelength $\varLambda /\delta$.

Figure 13. The same as figure 12 but for $A_R=30$: ($a$) $t=20$, ($b$) $t=30$ and ($c$) $t=50$. The results for the uniform shear layer ($A_R\rightarrow \infty$) with $Re=30$ are also shown for comparison.

Simulations with $(\varphi _x,\varphi _y)\neq (0,0)$ are also performed for the shear layer with $Re=30$ and $A_R=10$ or 30 and the perturbations with $\varLambda /\delta =12$. All computational parameters except the phases are the same as in the other simulations. Figure 14 shows the temporal variations of the maximum vorticity $\omega _{max}$ for various values of $\varphi _x$ and $\varphi _y$. The result for the unperturbed case ($u'_0=0$) is also shown for comparison. For $\varphi _y=0.25$, the spanwise vorticity of the perturbations is not zero at the centre of the shear layer and affects the initial value of $\omega _{max}$. Except for $(\varphi _x, \varphi _y)=(0,0.5)$ and $A_R=10$, $\omega _{max}$ in the perturbed cases increases with time. For $A_R=10$, the perturbations with $(\varphi _x, \varphi _y)=(0,0.5)$ result in a slow growth of $\omega _{max}$. The formation of the vortex is initially caused by the self-induced velocity of the shear layer with the low aspect ratio. The perturbations given with $(\varphi _x, \varphi _y)=(0,0.5)$ negate the self-induced velocity associated with a positive spanwise vorticity because of the matched phase of velocity distribution. Therefore, the self-induced velocity does not effectively distort the shear layer. The growth of $\omega _{max}$ for $A_R=30$ is observed regardless of the phases because the vortex formation begins with the distortion of the shear layer by the perturbations rather than the self-induced velocity.

Figure 14. The effects of the phases of perturbations on the shear layer development for ($a$) $A_R=10$ and ($b$) $A_R=30$.

The wavenumber dependence of the isolated shear layer with a finite aspect ratio $A_R$ agrees with previous studies of a parallel shear flow with $A_R\rightarrow \infty$. A comparison between parallel shear flows with and without a biaxial strain suggests that the most unstable mode of the shear instability is hardly influenced by the strain in Burgers’ vortex layer (Lin & Corcos 1984). The non-dimensional wavenumber in the linear stability analysis is often written as $\alpha =2{\rm \pi} \delta /\varLambda$. The maximum initial growth rate of the instability of Burgers’ vortex layer was found for $\alpha \approx 0.40\unicode{x2013}0.43$, which slightly increases with the Reynolds number (Lin & Corcos 1984). These values of $\alpha$ correspond to $\varLambda /\delta \approx 15$, which is close to the wavelength for which the vorticity in the isolated shear layer rapidly grows ($\varLambda /\delta \approx 12$) in figures 12 and 13. Therefore, the wavelength dependence observed for the shear layer with a finite $A_R$ can be related to the instability of the uniform shear layer with $A_R\rightarrow \infty$. One of the important differences is the vortex formation due to the self-induced velocity observed for the shear layer with low $A_R$. However, the present results have confirmed that the vortex formation by this mechanism is also promoted by perturbations with $\varLambda /\delta \approx 12$.

4. The response of small-scale shear layers to velocity perturbations in turbulence

Turbulence contains many small-scale shear layers with a finite aspect ratio that are subject to a biaxial strain. Therefore, the wavelength dependence of the perturbed isolated shear layer may be relevant to the evolution of shear layers in turbulence. Here, DNS is carried out for decaying turbulence to investigate the response of small-scale shear layers to velocity fluctuations with the efficient wavelength to promote shear instability. The results in § 3 suggest that the vortex formation from small-scale shear layers is promoted when they are subject to perturbations with $\varLambda /\delta \approx 12$ ($\varLambda /\delta _S\approx 7.2$). The visualisations of shear layers in turbulence indicate that the aspect ratio $A_R$ is different for each shear layer (Nagata et al. 2020; Fiscaletti et al. 2021; Hayashi et al. 2021a). However, the rapid growth of spanwise vorticity by perturbations with $\varLambda /\delta \approx 12$ is observed for both large and small $A_R$ in § 3. Furthermore, the p.d.f. of the layer thickness has a peak at $\delta _S\approx 4\eta$. Therefore, perturbations with a wavelength determined with $\delta _S\approx 4\eta$ and $\varLambda /\delta \approx 12$ ($\varLambda /\delta _S\approx 7.2$) are expected to affect the development of many shear layers in turbulence. This hypothesis is examined by simulating decaying turbulence initialised with an artificially modified velocity field of isotropic turbulence. Here, velocity fluctuations with the unstable wavelength for shear layers with a thickness of about $4\eta$ are artificially amplified or damped in isotropic turbulence. The response of shear layers to velocity fluctuations with a certain length scale is examined by comparing decaying turbulence initialised by the modified and original velocity profiles of isotropic turbulence.

4.1. Decaying turbulence with artificial velocity perturbations

DNS of decaying turbulence with artificial velocity perturbations is carried out to investigate the effects of external perturbations on the shear layers. The boundary conditions are the same as in the DNS of isotropic turbulence in § 2. The governing equations are (2.1) and (2.2) without external forcing, namely $f_i=0$. Non-dimensionalisation considered in § 3 is not adapted for DNS of isotropic turbulence, and appropriate normalisation is applied when the results are presented. The DNS code is also the same as in § 2. The initial velocity field is generated by superimposing solenoidal velocity perturbations ${\boldsymbol {u}}_{P}(x,y,z)$ on an instantaneous flow field ${\boldsymbol {u}}_{HIT}(x,y,z)$ of forced homogeneous isotropic turbulence. Two types of perturbations are tested here. One is defined in a deterministic way with a sinusoidal function whereas another is given by a random velocity field with a prescribed energy spectrum. Both velocity fields are divergence-free and are characterised by a single length scale. For the sinusoidal velocity perturbations, ${\boldsymbol {u}}_{P}(x,y,z)$ is given by

(4.1)\begin{equation} {\boldsymbol{u}}_{P}(x,y,z) =\left(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)\right), \end{equation}

where the amplitude $u_f$ and wavelength $\lambda _f$ are the computational parameters. As the velocity and length scales of small-scale shear layers in turbulence are characterised by the Kolmogorov scales, $u_f$ and $\lambda _f$ are determined in terms of $u_{\eta }$ and $\eta$.

Random velocity perturbations are generated by applying the inverse Fourier transform to velocity vectors in wavenumber space, $\hat {\boldsymbol {u}}_{P}(k_x,k_y,k_z)$, whose three-dimensional energy spectrum $E(k)$ is given by a Gaussian function as

(4.2)\begin{equation} E(k)= A_{E}\exp\left[ -\frac{(\log k-\log k_p)^2}{2\delta_E^2} \right], \end{equation}

with $k=\sqrt {k_x^2+k_y^2+k_z^2}$. The spectrum has a peak at $k=k_p=2{\rm \pi} /\lambda _f$ whereas the width of the distribution is $\delta _E=C_{\eta } \log k_\eta$ with $k_\eta =2{\rm \pi} /\eta$ and $C_\eta =0.005$. With two uniform random numbers that determine the phases, $\hat {\boldsymbol {u}}_{P}(k_x,k_y,k_z)$ can be calculated with (4.2) and the divergence-free condition in the wavenumber space following Johnsen et al. (2010). Then, the inverse Fourier transform of $\hat {\boldsymbol {u}}_{P}(k_x,k_y,k_z)$ yields ${\boldsymbol {u}}_{P}(x,y,z)$. The solenoidal velocity vector generated in this method has random phases in the wavenumber space and the energy spectrum with a single peak for $k_p=2{\rm \pi} /\lambda _f$. The constant $A_{E}$ is determined such that the increment of r.m.s. velocity fluctuations due to the perturbations is equal to $\Delta u_{rms}$, which is a parameter of the simulation. Figure 15 shows $E(k)$ and $u_{P}$ used in one of the simulations with $\lambda _f/\eta =30$. A large peak appears at a wavenumber of $k_p$. The length corresponding to $k_p$ characterises the profile of $u_{P}$, which does not have coherent patterns because of the random phases.

Figure 15. Random velocity perturbations used for the initial condition of D3L30R: ($a$) the energy spectrum $E(k)$ and ($b$) velocity in the $x$ direction, $u_P$, with $\lambda _f/\eta \approx 30$. In ($a$), $E(k)$ is normalised for the peak to be 1.

Table 3 summarises the parameters of DNS. For statistical analyses, DNS for each case is repeated $N_S$ times with different initial velocity profiles of ${\boldsymbol {u}}_{HIT}$, which are obtained from different time instances of statistically steady homogeneous isotropic turbulence. Here, $N_S$ is larger for a lower Reynolds number because the number of small-scale structures contained in the computational domain decreases as the Reynolds number becomes small and more simulations have to be carried out for statistical convergence (Hayashi et al. 2021a). The perturbations are generated with sinusoidal functions (4.1) for ‘Sin’ and random velocity for ‘Random’. The results of NS1, NS2 and NS3 in § 2 are used as the initial conditions. Therefore, the computational domain size and the number of grid points $N^3$ in the DNS of decaying turbulence are the same as those for the statistically steady isotropic turbulence. The table also includes the initial turbulent Reynolds number $Re_{\lambda }$. Cases D1, D2 and D3 consider decaying turbulence without perturbations, namely $u_{f}=0$. Cases D$n$L$m$ apply sinusoidal perturbations with $\lambda _f/\eta \approx m$ to the velocity profile taken from NS$n$ whereas cases D$n$L$m$R adapt random velocity perturbations. As the periodic boundary conditions are applied in the DNS, $\lambda _f$ is chosen such that ${\boldsymbol {u}}_{P}$ of (4.1) also satisfies the periodic boundary conditions. DNS is performed for $\lambda _f/\eta \approx 30$ and 70, which are determined with the shear layer thickness examined in § 2 and the wavenumber dependence of the isolated shear layer in § 3 as explained below. In each of the $N_S$ simulations, the initial value of $\eta$ is slightly different because different instances of isotropic turbulence are used in the initialisation. Therefore, the average of $\lambda _f/\eta$ taken for $N_S$ simulations is listed in table 3 whereas the r.m.s. fluctuations of $\lambda _f/\eta$ among $N_S$ simulations are less than 1 % of the averages. For all Reynolds numbers, $u_f/u_{\eta }=0$ or $1.0$ is assumed for the initial perturbations. In addition, the simulations are also performed with $u_f/u_{\eta }=1.2$ and 1.4 (cases D1L30u1 and D1L30u2) for the series of D1, which examine the influence of the perturbation amplitude. For random velocity in D$n$L$m$R, $A_{E}$ in (4.2) is determined for $\Delta u_{rms}/u_{rms0}$ to match the corresponding cases with sinusoidal perturbations, D$n$L$m$. The unperturbed and perturbed cases, D$n$ and D$n$L30, are compared for each Reynolds number. Time is advanced until $t=15\tau _{\eta 0}$, where $\tau _{\eta 0}$ is the Kolmogorov time scale at $t=0$.

Table 3. Parameters of DNS of decaying turbulence with and without initial perturbations: the turbulent Reynolds number $Re_\lambda$ of the original isotropic turbulence used for the initial condition, the number of grid points $N^3$, the number of simulations repeated for ensemble averages, $N_S$, the amplitude $u_f$ and wavelength $\lambda _f$ of the perturbations and the increase in r.m.s. velocity fluctuations due to the perturbations normalised by r.m.s. velocity fluctuations of the original isotropic turbulence, $\Delta u_{rms}/u_{rms0}$.

The parameters of the perturbations are compared with the characteristics of shear layers in § 2 and the isolated shear layer in § 3. For isotropic turbulence, the typical thickness of shear layers is $\delta _S/\eta \approx 4$, which is estimated from the peak of the joint p.d.f. in figure 5(a). The numerical simulations in § 3 suggest that the vortex formation is effectively promoted by perturbations with a wavelength of $\varLambda /\delta \approx 12$, which is equivalent to $\varLambda /\delta _S\approx 7.2$ because of $\delta _S/\delta =1.67$. Thus, $\varLambda /\delta \approx 12$ corresponds to $\varLambda /\eta =(\varLambda /\delta _S)(\delta _S/\eta )\approx 29$ for small-scale shear layers in isotropic turbulence. The wavelength of the perturbations is close to $\varLambda /\eta \approx 29$ in cases D$n$L30, for which the perturbations are expected to have significant influences on the evolution of the shear layers. In addition, case D1L70 considers the perturbations with $\varLambda /\eta = 71.7$, which corresponds to $\varLambda /\delta \approx 30$. The results in § 3 suggest that the perturbations with this wavelength are not as conducive to the vortex formation as those with $\varLambda /\eta \approx 30$ ($\varLambda /\delta \approx 12$).

The typical velocity jump across the shear layers in turbulence is about $u_S/u_{\eta }=5$, which is estimated from the joint p.d.f. in figure 5(b). The velocity jump $u_S$ is related to the reference velocity scale $u_r$ used in the simulations of the isolated shear layer by $u_S/u_r\approx 0.9$. Thus, $u_f/u_{\eta }=1$, 1.2 and 1.4 used for the perturbations of decaying turbulence are equivalent to $u_f/u_r=0.18$, 0.22 and 0.25, respectively. These amplitudes are larger than those considered for the isolated shear layer in § 3, and the perturbations are expected to be strong enough to affect the shear layer development in turbulence. The artificial velocity fluctuations affect the initial r.m.s. velocity fluctuations. The increase of the r.m.s. velocity fluctuations, $\Delta u_{rms}$, normalised by the original r.m.s. velocity fluctuations of isotropic turbulence, $u_{rms0}$, is also listed in table 3. For all cases, $\Delta u_{rms}/u_{rms0}$ is smaller than 5 %. Because the amplitude of perturbations is determined with the Kolmogorov velocity scale, $\Delta u_{rms}/u_{rms0}$ decreases with the Reynolds number. For the highest $Re$ case, D3L30, $\Delta u_{rms}/u_{rms0}$ is less than 1 %. These small values confirm that the perturbations are much weaker than the velocity fluctuations of isotropic turbulence.

The perturbations with a length of $\lambda _f$ cause a spike in an energy spectrum. Turbulence with a spectral spike has also been considered in other studies. These studies mostly concern a fully developed state of turbulence with a spectral spike, whose effects are already felt throughout turbulent motions of all length scales. Yeung & Brasseur (1991) investigated turbulence generated by anisotropic vortical forcing, which also causes a spectral spike for a wavelength of forced vortices. Their analysis concerns fully developed turbulence generated by vortical forcing, by which large-scale anisotropy affects small-scale turbulent motions. Similar anisotropic turbulence with a spectral spike was also studied in Goto & Vassilicos (2016). Another example is a study of a turbulent mixing layer in Takamure et al. (2019), where large-scale coherent vortices cause spikes in the energy spectrum. They have shown the relevance of the spectral spikes to the non-equilibrium scaling of the energy dissipation rate. These studies have discussed the influences of spectral spikes on the statistics of turbulence. The perturbation response of small-scale shear layers in the present DNS can be interpreted as the transient regime of the interaction between the spectral spike and small-scale velocity fluctuations of isotropic turbulence, where the influence of the perturbations gradually prevails. In addition, the present study mainly discusses small-scale vortices and shear layers, which are usually difficult to identify in wavenumber space. However, the influence of the perturbations observed for shear layers in physical space should affect the flow evolution in wavenumber space, and the insight from the present analysis can be useful for future studies of the scale-by-scale analysis in wavenumber space.

4.2. Decaying turbulence initialised with a band-cut filtered velocity field

Fluid motions with the efficient wavelength to promote small-scale shear instability also exist in turbulence because this length, namely about $30\eta$, is greater than the Kolmogorov scale. The role of velocity fluctuations of turbulence in the shear layer development is also investigated with DNS of decaying turbulence. Although the interaction between turbulent motions with a length scale of about $30\eta$ and small-scale shear layers may be understood as the interaction of different scales studied with Fourier analysis, it is difficult to explicitly examine the influence of velocity fluctuations at a particular length scale on the shear layers because the shear layers cannot be identified in wavenumber space. Therefore, the present study takes an alternative approach to examine the response of shear layers to velocity fluctuations of turbulence. The DNS of decaying turbulence is initialised with a velocity field which is obtained by applying a band-cut filter to a three-dimensional velocity field of isotropic turbulence. The filter damps velocity fluctuations of a certain length scale. Then, flow evolution is compared for the initial conditions given by the filtered and original velocity profiles to determine whether or not the velocity fluctuations damped by the filter are important in the shear layer development. A similar approach was used in Cimbala, Nagib & Roshko (1988), Jiménez & Moin (1991) and Jiménez (2018), where large-scale turbulent structures or vorticity fluctuations within a region of interest were artificially damped or eliminated to examine their influences on flow evolution.

Following Leung, Swaminathan & Davidson (2012) and Doan, Swaminathan & Chakraborty (2017), the filter is applied in wavenumber space. First, the three-dimensional Fourier transform is applied to the velocity field of isotropic turbulence to obtain the velocity vector in the wavenumber space, $\hat {\boldsymbol {u}}(k_x,k_y,k_z)$. Then, $\hat {\boldsymbol {u}}$ is multiplied by the following band-cut filter $T(k)$ defined with $k=\sqrt {k_x^2+k_y^2+k_z^2}$:

(4.3)\begin{equation} T(k)= 1-\exp\left[ -\frac{(\log k-\log k_p)^2}{2\delta_E^2} \right]. \end{equation}

This function is similar to the energy spectrum in (4.2), and $T(k)\approx 0$ for $k\approx k_p$ and $T(k)=1$ otherwise. Then, the band-cut filtered velocity field is obtained by the inverse Fourier transform applied to $T\hat {\boldsymbol {u}}$ and is used as an initial condition of DNS of decaying turbulence. The wavelength $k_p$ is determined as $k_p=2{\rm \pi} /\lambda _f$ while $\delta _E$ is $\delta _E=C_{\eta } \log k_\eta$ with $C_\eta =0.005$, as also assumed for the perturbations given by (4.2). The filter damps velocity fluctuations with a length scale of $\lambda _f$. We have also tested different filters, such as a band-cut filter defined with a top-hat band-pass filter in Fourier space (Cao, Chen & Doolen 1999; Baerenzung et al. 2010; Goto, Saito & Kawahara 2017; Hirota et al. 2017), and have confirmed that the results are not sensitive to the choice of filters.

Filtered velocity fields for the initial condition are prepared with the results of NS3. The same snapshots of NS3 have been used in D3, D3L30 and D3L30R presented above. The cutoff wavelength of the filter, $k_p=2{\rm \pi} /\lambda _f$ in (4.3), is given by $\lambda _f=30\eta$ or $140\eta$. The filter with $\lambda _f=140\eta$ damps large-scale velocity fluctuations because $140\eta$ corresponds to 0.6 times the integral scale. On the other hand, velocity fluctuations of the unstable wavelength for shear layers are damped by the filter with $\lambda _f=30\eta$. Cases with the filters of $\lambda _f=30\eta$ and $140\eta$ are referred to as D3L30F and D3L140F, respectively. Computational parameters other than the initial conditions are the same as in D3. The results are compared for D3, D3L30F and D3L140F to reveal the role of the filtered velocity component in the shear layer development.

4.3. Results and discussion

The results of DNS with artificial perturbations in § 4.1 are presented first whereas the role of internal velocity fluctuations in turbulence is discussed later. Figure 16(a) shows the temporal evolution of a one-dimensional energy spectrum of $v$, $E_v$, for a wavenumber in the $x$ direction, $k_x$, in D3L30. The results are shown from $t/\tau _{\eta 0}=0$ to 6. The initial perturbations for $v$ are given in a form of $\sin (x)$, which contributes to the peak in $E_v(k_x)$. This peak decays with time and is no longer identified after $t/\tau _{\eta 0}=4.8$. As the perturbations added to the initial field is weak, it hardly influences the overall spectral shape. Figure 16(b) presents a comparison between D3 and D3L30 for the decay of r.m.s. velocity fluctuations normalised by the value of isotropic turbulence. At $t=0$, the perturbations in D3L30 cause $u_{rms}$ to increase by about 0.6 % compared with the unperturbed case (D3). This increase is not large enough to significantly alter the decay rate of $u_{rms}$. However, $u_{rms}$ in D3L30 seems to gradually approach that in D3, implying a larger energy dissipation rate in D3L30.

Figure 16. ($a$) Energy spectrum of $v$, $E_v$, calculated for the wavenumber in the $x$ direction, $k_x$, in D3L30. ($b$) The decay of r.m.s. velocity fluctuations $u_{rms}$ in D3 and D3L30. The results are non-dimensionalised with the r.m.s. velocity fluctuations $u_{rms0}$, integral scale $L_{0}$ and the Kolmogorov time scale $\tau _{\eta 0}$ of statistically stationary homogeneous isotropic turbulence used for the initial condition.

The number of vortices (vortex tubes) is compared for decaying turbulence with and without initial perturbations. As also done in previous studies (Jiménez & Wray 1998), the vortices are detected by thresholding a quantity related to rotating motion. The intensity of rotating motion $I_R$ defined with the triple decomposition of a velocity gradient tensor is used in the vortex identification in this study. The vortices identified with $I_R$ are visualised in figure 1 for isotropic turbulence while those in other flows, such as wakes, jets and boundary layers, can be found in Kolář & Šístek (2014) and Hayashi et al. (2021a). Vortices in turbulence can also be detected with so-called $Q$ and $\lambda _2$ criteria (Chong, Perry & Cantwell 1990; Jeong & Hussain 1995). Vortices identified with $I_R$ were compared with those with $Q$ and $\lambda _2$ criteria in Šístek et al. (2012) and Kolář & Šístek (2014), where the difference among these identification methods was shown to be small for various flows, such as a transitional boundary layer and a turbulent wake of a flat plate. This is because positive $Q$ is strongly correlated with $I_R$ in turbulence (Nagata et al. 2020). The results presented in this paper do not change even if the vortices are identified with the $Q$ criterion, as shown below.

The algorithm to identify vortices with $I_R$ is the same as that used with enstrophy for various turbulent flows (Jiménez et al. 1993; Jiménez & Wray 1998; da Silva et al. 2011; Jahanbakhshi et al. 2015; Ghira et al. 2022) although the enstrophy is not used in this study because it cannot distinguish between vortex tubes and shear layers. The computational grid points that belong to vortices are identified with a given threshold $I_{Rth}$ as points with $I_{R}>I_{Rth}$. Then, the connectivity of the vortex points is examined to define the vortices as connected points with $I_{R}>I_{Rth}$, which belong to the same vortex. Then, the number of points belonging to each vortex is counted. Following Jiménez & Wray (1998) and da Silva et al. (2011), the vortices for which less than 30 points are assigned are discarded from the subsequent analysis because these vortices have $I_R$ close to the threshold and appear as noise-like patterns rather than coherent structures. Finally, the number of vortices detected with $I_{R}$, $N_{I_R}$, is counted. As each simulation is repeated $N_S$ times, $N_{I_R}$ is averaged for $N_S$ simulations.

The number of detected vortices $N_{I_R}$ depends on the choice of the threshold $I_{Rth}$. Instead of using a specific threshold, $N_{I_R}$ is evaluated for a wide range of $I_{Rth}$. Figure 17 plots $N_{I_R}$ with the normalised threshold $I_{Rth}/\langle I_R\rangle$ at $t/\tau _{\eta 0}=0$, 2 and 6 for D1, D1L30 and D1L70. Here, $\langle I_R\rangle$ is used for normalisation because $\langle I_R\rangle$ decreases with time and becomes large with the Reynolds number. Except for very small $I_{Rth}/\langle I_R\rangle$, $N_{I_R}$ decreases with an increase of $I_{Rth}/\langle I_R\rangle$ because the vortices with moderately large $I_R$ are not identified with a large threshold. However, $N_{I_R}$ also decreases as $I_{Rth}/\langle I_R\rangle$ becomes small for $I_{Rth}/\langle I_R\rangle \lesssim 1.5$. Vortex regions identified with a very small threshold are connected even if these regions belong to different vortices. Therefore, the threshold should be chosen from the range where $N_{I_R}$ decreases with increasing $I_{Rth}$, namely $I_{Rth}/\langle I_R\rangle \gtrsim 1.5$, for which the effects of the initial perturbations are discussed by comparing D1 with D1L30 and D1L70. At $t=0$ in figure 17(a), $N_{I_R}$ hardly differs for the three cases, and the perturbations do not alter the initial distribution of vortices. However, a difference in $N_{I_R}$ between D1L30 and the others becomes clearer with time. For all thresholds, $N_{I_R}$ for D1L30 is larger than those for D1 and D1L70 at $t/\tau _{\eta 0}=2$ and 6, and more vortices are found in the flow when the initial velocity profile is perturbed by velocity fluctuations with $\lambda _f/\eta \approx 30$. The perturbations are also added in D1L70, and the increase of initial $u_{rms}$ caused by the perturbations is the same for D1L30 and D1L70. However, the increase in the number of vortices due to the perturbations does not occur for D1L70. These results are qualitatively the same for all thresholds.

Figure 17. Threshold dependence of the number of vortices, $N_{I_R}$, in D1, D1L30 and D1L70 at ($a$) $t=0$, ($b$) $t/\tau _{\eta 0}=2$ and ($c$) $t/\tau_{\eta 0}=6$. The vortices are detected by the intensity of rotating motion $I_R$ with a threshold $I_{Rth}$. $N_{I_R}$ is plotted as a function of $I_{Rth}$ normalised by the average of $I_{R}$, $\langle I_R\rangle$.

The temporal variations of $N_{I_R}$ are examined with several thresholds for D1, D1L30 and D1L70. As the threshold increases, $N_{I_R}$ decreases because the vortices with moderately strong rotation are disregarded for a large threshold. However, as discussed above, a very small threshold results in a decrease of $N_{I_R}$ because multiple vortices are identified as one large structure. These tendencies suggest that the threshold corresponding to the maximum $N_{I_R}$ is appropriate to detect vortices with both large and small $I_R$ without the problem of vortex connection caused by a small threshold. For this case, the threshold value normalised by $\langle I_R\rangle$ also varies with time. Figure 18(a) plots the temporal variations of $N_{I_R}$ calculated with the threshold corresponding to the maximum $N_{I_R}$. When the perturbations are not added to the initial velocity field (D1), $N_{I_R}$ monotonically decreases with time. However, $N_{I_R}$ increases until $t/\tau _{\eta 0}\approx 4$ in D1L30, for which more vortices exist in the flow compared with the unperturbed case. This influence of the perturbations is observed only for D1L30, where the perturbation wavelength matches the anticipated unstable wavelength of typical small-scale shear layers. When the perturbations have a wavelength of $\lambda _f/\eta \approx 70$ (D1L70), the temporal variation of $N_{I_R}$ hardly differs from the unperturbed case (D1) even though the perturbations imposed in D1L30 and D1L70 have the same turbulent kinetic energy as confirmed from the same value of $\Delta u_{rms}/u_{rms0}$ in table 3.

Figure 18. Temporal variations of the number of detected vortices $N_{I_R}$ in D1, D1L30 and D1L70. The vortices are detected with the intensity of rotating motion $I_R$ with ($a$) the threshold that yields the maximum number of vortices at each time step and ($b$) those determined with the average of $I_R$ as $C_{th}\langle I_R\rangle$ with $C_{th}=2$, 3 and 4 whereas the results obtained with the $Q$ criterion are presented in ($c$), where the threshold is determined such that the number of detected vortices is the maximum at each time step. Panel ($a$) also presents the number of shear layers, $N_{I_S}$, detected with $I_S$, for which the threshold is determined based on the threshold dependence of $N_{I_S}$.

Shear layers can be detected with the same method used for vortices by using $I_S$ instead of $I_R$. The threshold $I_{Sth}$ is also determined by the same method as $I_{Rth}$ based on the $I_{Sth}$ dependence of the number of detected shear layers. Figure 18($a$) shows the number of shear layers, $N_{I_S}$, for the three simulations. Both $N_{I_S}$ and $N_{I_R}$ in D1L30 are larger than those in D1 and D1L70, and the perturbations in D1L30 result in the increase of both vortices and shear layers. The kinetic energy budget around shear layers suggests that the energetic flow that contributes to the shear is generated and sustained by vortices (Watanabe & Nagata 2022). Furthermore, when the roll-up of a shear layer occurs, it splits into more than one shear layer around the generated vortex (Corcos & Lin 1984). Therefore, more shear layers appear when the number of vortices is increased by the influence of perturbations.

Figure 18($b$) presents $N_{I_R}$ obtained with the threshold determined with the average of $I_R$ as $I_{Rth}=C_{th}\langle I_R\rangle$. The temporal variations of $N_{I_R}$ in D1, D1L30 and D1L70 are compared for three thresholds of $C_{th}=2$, 3 and 4. Although $N_{I_R}$ varies with $C_{th}$, the temporal variations of $N_{I_R}$ are qualitatively similar for all $C_{th}$. Therefore, the increase in the number of vortices due to the perturbations in D1L30 can be observed regardless of the threshold. Finally, the effects of the perturbations are also examined for the vortices detected with the second invariant of $\boldsymbol {\nabla } \boldsymbol {u}$, $Q=(\textsf{$\mathit{\Omega}$} _{ij}\textsf{$\mathit{\Omega}$} _{ij}-{\mathsf{S}}_{ij}{\mathsf{S}}_{ij})/2$, for which the region of $Q>Q_{th}$ is identified as vortices. Here, $Q_{th}$ is determined such that the number of detected vortices, $N_Q$, is the maximum. Figure 18(c) plots $N_Q$ with time. As also found for the vortices identified with $I_R$, more vortices appear in D1L30 whereas the number of vortices is not influenced by the perturbations in D1L70. These results confirm that the increase in the number of vortices in D1L30 is a physically meaningful result and is not an artefact due to the vortex identification methods.

Figure 19 visualises the evolution of one of the shear layers. The shear layers (white) and vortices (orange) are visualised with the isosurfaces of $I_S$ and $I_R$, respectively. Figures 19(a–d) and 19(e–h) show the results of D1 and D1L30, respectively. A large shear layer is observed at $t=0$ in both cases (figure 19a,e). As this shear layer evolves with time, vortex A appears in D1L30 in figure 19(g). The shear layer which generates vortex A is disappearing in figure 19(h). However, this shear layer keeps its flat shape for a long time in D1, where no perturbations are added to the initial velocity. In figure 19(d,h), vortex B also fully develops from the shear layer in D1L30 whereas this vortex is still developing in the shear layer in D1. A comparison between figures 19(c,d) and 19(g,h) suggests that the lifetime of the shear layers is shortened by the perturbations. The vortex formation from the shear layers naturally occurs in turbulence (Vincent & Meneguzzi 1994; Passot et al. 1995; Watanabe et al. 2020). This process visualised in Watanabe et al. (2020) looks similar to the evolution of the shear layers and vortices in figure 19(e–h). Therefore, the perturbations in D1L30 do not directly destroy the shear layers at $t=0$ but promote the vortex formation from the shear layers.

Figure 19. Temporal evolutions of shear layers and vortices, which are visualised with the isosurfaces of $I_S/\langle I_S\rangle =2$ (white) and $I_R/\langle I_R\rangle =6$ (orange), respectively. The colour contour of $I_S$ is shown on the surface of the visualised domain. (a–d) Unperturbed case D1 and (e–h) perturbed case D1L30. Time increases from left to right as (a,e) $t/\tau _{\eta 0}=0$, (b,f) $t/\tau _{\eta 0}=2$, (c,g) $t/\tau _{\eta 0}=4$ and (d,h) $t/\tau _{\eta 0}=6$.

Figure 20(a) shows the variations of the number of vortices identified with $I_R$ for D2, D3, D2L30 and D3L30. Here, the threshold is determined as $I_{Rth}=C_{th}\langle I_R\rangle$ with $C_{th}=3$. Because the vortices are small-scale structures characterised by the Kolmogorov scale, which becomes small with increasing $Re_\lambda$, $N_{I_R}$ increases from D1 to D3. However, the effect of the initial perturbations is similar for all Reynolds numbers, and $N_{I_R}$ in the perturbed cases with $\lambda _f/\eta \approx 30$ becomes large compared with the unperturbed cases.

Figure 20. Temporal variations of ($a$) the number of vortices, $N_{I_R}$, in D2, D2L30, D3 and D3L30 and ($b$) a ratio of $N_{I_R}$ between the perturbed and unperturbed cases.

Figure 20(b) shows the variations of the ratio of $N_{I_R}$ between the perturbed and unperturbed cases. The results for sinusoidal and random perturbations are also compared in the figure. For both perturbations, the ratio exceeds 1 for $\lambda _f/\eta \approx 30$ (D$n$L30 and D$n$L30R), and the increase in the number of vortices is observed for this wavelength regardless of the Reynolds number and the types of perturbations. Because the shear layers are located almost randomly in the flow, the relative positions of the shear layers to the peak of sinusoidal perturbations are not constant. Therefore, the effects of sinusoidal perturbations are similar to those of random velocity perturbations. As the number of vortices increases, the ratio reaches a peak at $t/\tau _{\eta 0}\approx 7$ and then decays with time. The time corresponding to the peak hardly depends on the Reynolds number when the Kolmogorov time scale is used for normalisation. Therefore, how the perturbations affect the shear layers is a small-scale process whose time scale is characterised by the Kolmogorov time scale. Visualisation of the temporal evolution of small-scale shear layers in isotropic turbulence has shown that the natural roll-up of the shear layer is also completed within about five times the Kolmogorov time scale (Watanabe et al. 2020). This time scale of the roll-up is close to $t/\tau _{\eta 0}\approx 7$, at which the ratio of $N_{I_R}$ reaches the peak in figure 20(b). The peak value at $t/\tau _{\eta 0}\approx 7$ is about 1.07 for D1L30, D2L30 and D3L30 despite the difference in the Reynolds number. For these cases, the amplitude of perturbations, $u_f$, is determined based on the Kolmogorov velocity scale $u_{\eta }$ as $1.0u_{\eta }$ whereas $u_{\eta }/u_{rms}$ decreases with the Reynolds number. Therefore, the increment in the number of vortices for a given wavelength depends on $u_f/u_{\eta }$. Because $u_{\eta }$ decreases with the Reynolds number, less turbulent kinetic energy is required for the perturbations to cause an increase in the number of vortices at a higher Reynolds number. A comparison among the series of D1 suggests that $N_{I_R}$ increases as $u_f/u_{\eta }$ becomes large. For $u_f/u_{\eta }=1.4$ in D1L30u2, $N_{I_R}$ increases by more than 10 % from the unperturbed case. The increase in the number of vortices does not occur in D1L70, for which the ratio stays around 1. These results are consistent with the response of the isolated shear layer to perturbations in § 3: the vortex formation rapidly occurs with the perturbations with $\lambda _f/\eta \approx 30$ estimated with $\lambda _f/\delta _S\approx 7.2$ and $\delta _S/\eta \approx 4$ while those with a larger wavelength are not efficient in promoting the vortex formation.

The above results have shown that the external perturbations with a wavelength of $30\eta$ promote the vortex formation from shear layers. The results for decaying turbulence initialised with a filtered velocity field explained in § 4.2 are presented to examine whether the shear instability in turbulence is promoted by internal velocity fluctuations with this wavelength. Figure 21 compares the variations of the number of vortices $N_{I_R}$ for D3, D3L30F and D3L140F, together with the ratio of $N_{I_R}$ between the filtered and original cases. The initial velocity fields of D3L30F and D3L140F do not contain velocity fluctuations with wavelengths of $30\eta$ and $140\eta$, respectively. A comparison of three cases shows that $N_{I_R}$ in D3L30F is smaller than those in the other cases, suggesting that fewer vortices are generated by the shear instability when the velocity fluctuations with a length of $30\eta$ are damped by the filter. The reduction of $N_{I_R}$ is not significant in D3L140F. Figures 21(b) and 21(c) visualise one of the shear layers and nearby vortices in D3 and D3L30F, respectively. Both figures visualise the same region of the computational domain at $t/\tau _{\eta 0}=6$, at which the reduction of the number of vortices is observed in D3L30F. The shear layer in figure 21(b) is being broken and a vortex has been generated at the centre of the shear layer. This shear layer has a flat shape in D3L30F without generating the vortex. Thus, the shear layer can survive a long time when the velocity fluctuations with $30\eta$ are damped in the initial field. These results confirm that turbulent motions with a length of $30\eta$ play a dynamically important role in vortex formation by small-scale shear instability.

Figure 21. ($a$) The number of vortices $N_{I_R}$ in decaying turbulence initialised with a band-cut filtered velocity field and an original velocity field of isotropic turbulence. The ratio of $N_{I_R}$ between the filtered and original cases is also shown for comparison. The vortices are detected as regions of $I_R\geqslant C_{th}\langle I_R\rangle$ with $C_{th}=3$, where $\langle I_R\rangle$ is evaluated with D3. Visualisation of shear layers (white) and vortices (orange) at $t/\tau _{\eta 0}=6$ in ($b$) D3 and ($c$) D3L30F. The isosurface values are the same as in figure 19.

The flow field around shear layers in the initial filtered velocity field of D3L30F is investigated by taking conditional averages around the shear layers. For each shear layer detected by the method described in § 2, variables on the DNS grid are interpolated on the shear coordinate $(\zeta _1, \zeta _2, \zeta _3)$. Then, the ensemble averages of many shear layers are calculated as functions of $(\zeta _1, \zeta _2, \zeta _3)$. This average is denoted by an overline, e.g. $\overline {I_S}$. Figure 22 shows the mean intensities of shear and rotation, $\overline {I_S}$ and $\overline {I_R}$, around the shear layers in the initial field of D3L30F. The averages are shown on the $\zeta _2$–$\zeta _3$ plane at $\zeta _1=0$. Figure 22(a) also shows the mean velocity vectors. The results for original isotropic turbulence were presented in Watanabe & Nagata (2022). The shear is associated with the mean flows in $\pm \zeta _3$ directions on the sides of the shear layer with large $\overline {I_S}$ centred at $(\zeta _2, \zeta _3)=(0,0)$. In figure 22(b), four approximately circular regions with large $\overline {I_R}$ are found around $(\zeta _2, \zeta _3)=(0,0)$. Two of them are located within the shear layer on $\zeta _2=0$ whereas the others appear on the sides of the shear layer. Similar rotating motions were also identified in the mean flow pattern in the strain eigenframe (Elsinga & Marusic 2010). The profiles of $\overline {I_S}$, $\overline {I_R}$ and mean velocity hardly differ from isotropic turbulence (Watanabe & Nagata 2022) even though the present results are obtained for the filtered velocity field, for which the instability of shear layers is inhibited. It is implied that rotating motions around the shear layers are not related to the velocity fluctuations which promote small-scale shear instability.

Figure 22. A mean flow field around shear layers in the initial velocity field of D3L30F, which is generated by applying the band-cut filter with the cutoff length of $\lambda _f=30\eta$ to isotropic turbulence. ($a$) Mean shear intensity $\overline {I_S}$ and mean velocity vectors and ($b$) mean intensity of rotating motion $\overline {I_R}$ on the $\zeta _2$–$\zeta _3$ plane at the centre of shear layers ($\zeta _1=0$).

The effects of the increase or decrease in the number of vortices are examined for the kinetic energy dissipation rate $\varepsilon =2\nu {\mathsf{S}}_{ij}{\mathsf{S}}_{ij}$. Figure 23(a) compares the temporal variation of the averaged dissipation rate $\langle \varepsilon \rangle$ for D3, D3L30 and D3L30F. The shear instability is promoted and inhibited in D3L30 and D3L30F, respectively. For $t/\tau _{\eta 0}\leq 10$, the dissipation rates in D3L30 and D3L30F are larger and smaller than that of the original case, D3. Thus, when the initial velocity field is modified such that the shear instability is promoted, the dissipation rate also increases whereas the opposite effect is found when the instability is inhibited. The increase in the number of vortices results in more small-scale shear layers. Although for the local fluid motions considered in the triple decomposition, both shear and elongation contribute to the rate-of-strain tensor, the shear has a greater contribution to ${\mathsf{S}}_{ij}{\mathsf{S}}_{ij}$ than the elongation (Nagata et al. 2020). Therefore, a large dissipation rate in D3L30 can be explained by the increased number of shear layers. For $t/\tau _{\eta 0}\geqslant 10$, an opposite trend for $\langle \varepsilon \rangle$ is observed: $\langle \varepsilon \rangle$ for the perturbed and filtered cases is smaller and larger than that for the original case. The enhanced dissipation in the perturbed case results in a rapid decay of turbulent kinetic energy at small scales. Then, there is less kinetic energy to be dissipated at a late time when the perturbations are added, and therefore $\langle \varepsilon \rangle$ can be small. Figure 23(b,c) compare the p.d.f. of $\varepsilon$ at $t/\tau _{\eta 0}=5$ for the three cases. In figure 23($b$), a logarithmic scale for $\varepsilon$ is used to compare the probability distribution for small $\varepsilon$. The p.d.f. suggests that small $\varepsilon$ is less frequently observed when the perturbations are added in D3L30. The peak of the p.d.f. for the filtered velocity in D3L30F is larger than that in D3, and the probability to observe small $\varepsilon$ increases in D3L30F. The p.d.f. for relatively large $\varepsilon$ is shown in figure 23($c$). The probability to observe large $\varepsilon$ increases and decreases when the velocity fluctuations with a length of $30\eta$ are added and damped, respectively.

Figure 23. ($a$) Temporal variations of averaged kinetic energy dissipation rate $\langle \varepsilon \rangle$ normalised by the initial value $\varepsilon _0$ in D3, D3L30 and D3L30F. ($b,c$) Probability density function of $\varepsilon$ at $t/\tau _{\eta 0}=5$.

Figure 24 shows the p.d.f. of the production rates of enstrophy and strain at $t/\tau _{\eta 0}=5$ in D3, D3L30 and D3L30F. The enstrophy production and strain self-amplification are written as $P_\omega =\omega _{i}{\mathsf{S}}_{ij}\omega _{j}$ and $P_S=-{\mathsf{S}}_{ij}{\mathsf{S}}_{jk}{\mathsf{S}}_{ki}$, respectively. The overall distribution of the p.d.f. is similar for all cases. However, the profile becomes slightly wider and narrower for D3L30 and D3L30F, respectively, than for D3. The insets confirm that the p.d.f. for large $P_\omega$ and $P_S$ increases when the shear instability is promoted by the perturbations whereas the opposite effect is found when the instability is inhibited. Consequently, the averages of $P_\omega$ and $P_S$ increase by 20 % by the perturbations in D3L30 and decrease by 4 % by the filter in D3L30F. Shearing motions have dominant contributions to $P_\omega$ and $P_S$ (Watanabe et al. 2020). The increases of $\langle P_\omega \rangle$ and $\langle P_S\rangle$ by the perturbations are also explained by the increased number of shear layers, which can be generated by vortices.

Figure 24. Probability density functions of ($a$) enstrophy production $P_\omega =\omega _i{\mathsf{S}}_{ij}\omega _j$ and ($b$) strain self-amplification $P_S=-{\mathsf{S}}_{ij}{\mathsf{S}}_{jk}{\mathsf{S}}_{ki}$ at $t/\tau _{\eta 0}=5$ in D3, D3L30 and D3F30F. The production rates, $P_\omega$ and $P_{S}$, are normalised by their initial mean values denoted by $P_{\omega 0}$ and $P_{S0}$.

5. Conclusions

The response of small-scale shear layers to velocity perturbations has been studied with the analysis of shear layers in isotropic turbulence, numerical simulations of a perturbed isolated shear layer in a biaxial strain field and DNS of decaying turbulence. As summarised below, vortex formation from the shear layers in turbulence is efficiently promoted by velocity fluctuations with a wavelength of approximately $30\eta$, which is the most unstable wavelength for the instability of small-scale shear layers.

The thickness $\delta _S$ and velocity jump $u_S$ of each shear layer have been examined with DNS of isotropic turbulence. The typical shear layers have $\delta _S/\eta \approx 4$ and $u_S/u_\eta \approx 5$ regardless of $Re_\lambda$. Many shear layers have the shear Reynolds number of $Re_S=\delta _Su_S/\nu \approx 20$ although $Re_S$ ranges approximately from 10 to 70. Furthermore, the equilibrium relation between a strain rate and the layer thickness of Burgers’ vortex layer is valid for most shear layers with $Re_S\approx 20$.

Numerical simulations have been carried out to investigate the development of an isolated shear layer with a finite aspect ratio $A_R$. The model considers a shear layer subject to a biaxial strain, which is observed in the mean flow field around shear layers in turbulence (Watanabe et al. 2020; Hayashi et al. 2021a). When the aspect ratio is small, the self-induced velocity of the shear layer causes the roll-up of the layer and generates a vortex. On the other hand, the shear layer with large $A_R$ is simply distorted without vortex formation if no perturbations are added. Perturbations can trigger the vortex formation in the large-$A_R$ shear layer, as also observed for a parallel shear flow with $A_R\rightarrow \infty$ in previous studies. For both cases of large and small $A_R$, the vortex formation is promoted by the perturbations whose wavelength $\varLambda$ divided by the layer thickness $\delta _S$ is $\varLambda /\delta _S\approx 7$, where $\delta _S$ is the half-width of shear intensity. The perturbations with smaller or larger wavelengths are not efficient in promoting the instability and sometimes delay or inhibit the vortex formation.

Many small-scale shear layers in turbulence are expected to be significantly influenced by perturbations with a wavelength of about $30\eta$, which is estimated from the typical thickness of small-scale shear layers in turbulence and the perturbation wavelength which is efficient in promoting vortex formation of the isolated shear layer. The response of turbulence to artificial velocity perturbations has been investigated with DNS of decaying turbulence initialised by superimposing solenoidal velocity perturbations with a length scale of $\lambda _f$ on isotropic turbulence. The number of vortices increases under the influence of perturbations with $\lambda _f/\eta \approx 30$, which promote the vortex formation from shear layers while the shear layers can survive for a longer time without the perturbations. The number of vortices remains unchanged for perturbations with $\lambda _f/\eta \approx 70$. How the perturbations affect shear layer development is a small-scale process whose time scale is given by the Kolmogorov time scale. These results are consistent with numerical simulations of the isolated shear layer, where perturbations only with a certain wavelength can promote vortex formation. The amplitude of velocity perturbations $u_f$ affects the increment in the number of vortices, which is more relevant to the ratio of $u_f$ to $u_{\eta }$ rather than that to $u_{rms}$. As $u_{\eta }/u_{rms}$ decreases with a Reynolds number, less turbulent kinetic energy is required for the perturbations to cause an increase in the number of vortices at a higher Reynolds number. For the present DNS, the number of vortices increases by more than 10 % when $u_f$ is $1.4u_{\eta }$. The shear layers are generated or sustained by energetic motions induced by the vortices (Watanabe & Nagata 2022). In addition, shear layer is split into more than one shear layer during the roll-up (Passot et al. 1995; Watanabe et al. 2020). Therefore, more shear layers appear in turbulence when the perturbations increase the number of vortices. The external perturbations which promote the shear instability also enhance the kinetic energy dissipation, enstrophy production and strain self-amplification because these are significantly enhanced by shearing motion.

DNS of decaying turbulence has also been carried out with an initial condition obtained by applying a band-cut filter to isotropic turbulence. When the velocity fluctuations with a wavelength of $30\eta$ are eliminated, the shear layers tend to persist for a long time, resulting in fewer vortices. The reduction in the number of vortices does not occur when the filter is applied to an energy-containing length scale. Thus, turbulent motions with a length of $30\eta$ play a dynamically important role in vortex formation by small-scale shear instability. Vortices frequently appear near and within shear layers (Horiuti & Takagi 2005; Elsinga & Marusic 2010; Watanabe & Nagata 2022; Nakamura et al. 2023). However, the filter that leads to the suppression of the shear instability does not affect these vortices. Therefore, rotating motions observed around the shear layers are not directly related to the velocity fluctuations which efficiently promote the shear instability.

The present study has confirmed that perturbations with the unstable wavelength of small-scale shear instability can efficiently modulate small-scale turbulent structures. The effects observed for the number of vortices can be important in many phenomena related to vortices in turbulence. One of the examples is turbulent entrainment in intermittent flows such as jets, mixing layers and boundary layers, where a turbulent/non-turbulent interface appears between turbulent and outer non-turbulent regions (da Silva et al. 2014). Previous studies of the interface have revealed that small-scale vortices near the interface are relevant to the local entrainment of non-turbulent fluids (Mathew & Basu 2002; Holzner & Lüthi 2011; Taveira & da Silva 2014; Watanabe et al. 2016, 2017a; Neamtu-Halic et al. 2021). The increase in the number of vortices by perturbations will be possibly important in the enhancement of entrainment by external disturbance. Further investigations are required for future applications of small-scale shear instability to turbulence modulation.