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. Reference Watanabe, Tanaka and Nagata2020; Watanabe, Tanaka & Nagata Reference Watanabe, Tanaka and Nagata2021; Watanabe & Nagata Reference Watanabe and Nagata2022). Turbulence is sustained at a statistically steady state with a linear forcing scheme (Carroll & Blanquart Reference Carroll and Blanquart2013) 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 Reference Carroll and Blanquart2013). 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. Reference Morinishi, Lund, Vasilyev and Moin1998) and a third-order and low-storage Runge–Kutta method for temporal advancement (Spalart, Moser & Rogers Reference Spalart, Moser and Rogers1991). 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 Reference Rosales and Meneveau2005; Carroll & Blanquart Reference Carroll and Blanquart2013; Watanabe et al. Reference Watanabe, Tanaka and Nagata2021). 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. Reference Eisma, Westerweel, Ooms and Elsinga2015; Watanabe et al. Reference Watanabe, Tanaka and Nagata2020; Fiscaletti et al. Reference Fiscaletti, Buxton and Attili2021). 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ář Reference Kolář2007; Nagata et al. Reference Nagata, Watanabe, Nagata and da Silva2020; Fiscaletti et al. Reference Fiscaletti, Buxton and Attili2021) and the shear layer analysis (Watanabe & Nagata Reference Watanabe and Nagata2022) .

The shear layers are detected with the triple decomposition (Kolář Reference 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. Reference Nagata, Watanabe, Nagata and da Silva2020; Fiscaletti et al. Reference Fiscaletti, Buxton and Attili2021). 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. (Reference Watanabe, Tanaka and Nagata2020), 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 Reference Hayashi, Watanabe and Nagata2021b). 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. (Reference Watanabe, Tanaka and Nagata2020) and Hayashi et al. (Reference Hayashi, Watanabe and Nagata2021a).

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. Reference Watanabe, Tanaka and Nagata2020; Fiscaletti et al. Reference Fiscaletti, Buxton and Attili2021; Hayashi et al. Reference Hayashi, Watanabe and Nagata2021a).

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 Reference Buxton and Ganapathisubramani2010; Elsinga & Marusic Reference Elsinga and Marusic2010; Eisma et al. Reference Eisma, Westerweel, Ooms and Elsinga2015; Elsinga et al. Reference Elsinga, Ishihara, Goudar, da Silva and Hunt2017; Watanabe et al. Reference Watanabe, Tanaka and Nagata2020). 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. Reference Elsinga, Ishihara, Goudar, da Silva and Hunt2017; Buaria et al. Reference Buaria, Pumir, Bodenschatz and Yeung2019; Das & Girimaji Reference Das and Girimaji2020; Ghira et al. Reference Ghira, Elsinga and da Silva2022).

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. Reference Watanabe, Tanaka and Nagata2020), which is one of the exact steady solutions of Navier–Stokes equations for a uniform shear layer subject to a biaxial strain (Burgers Reference Burgers1948; Davidson Reference Davidson2004). 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 Reference Jiménez and Wray1998; da Silva et al. Reference da Silva, Dos Reis and Pereira2011; Jahanbakhshi et al. Reference Jahanbakhshi, Vaghefi and Madnia2015; Watanabe et al. Reference Watanabe, da Silva, Nagata and Sakai2017b; Ghira et al. Reference Ghira, Elsinga and da Silva2022). 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. (Reference Watanabe, Tanaka and Nagata2020) and Hayashi et al. (Reference Hayashi, Watanabe and Nagata2021a). 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 Reference Davidson2004). 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 Reference Elsinga and Marusic2010; Elsinga et al. Reference Elsinga, Ishihara, Goudar, da Silva and Hunt2017). 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 Reference Elsinga and Marusic2010), which is similar to the mean flow around the shear layers (Watanabe & Nagata Reference Watanabe and Nagata2022). 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. Reference Nagata, Watanabe, Nagata and da Silva2020; Hayashi et al. Reference Hayashi, Watanabe and Nagata2021a). 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.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. Reference Watanabe, Tanaka and Nagata2020; Hayashi et al. Reference Hayashi, Watanabe and Nagata2021a; Watanabe & Nagata Reference Watanabe and Nagata2022), 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 Reference Elsinga and Marusic2010; Elsinga et al. Reference Elsinga, Ishihara, Goudar, da Silva and Hunt2017). 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 Reference Lin and Corcos1984) 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 (Reference Lin and Corcos1984), 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 Reference Lin and Corcos1984). 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 Reference Watanabe and Nagata2018; Watanabe, Zhang & Nagata Reference Watanabe, Zhang and Nagata2019; Katagiri, Watanabe & Nagata Reference Katagiri, Watanabe and Nagata2021; Watanabe, Zheng & Nagata Reference Watanabe, Zheng and Nagata2022). 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. Reference Morinishi, Lund, Vasilyev and Moin1998).

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. (Reference Watanabe, Riley, Nagata, Onishi and Matsuda2018). 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. Reference Nagata, Watanabe, Nagata and da Silva2020; Fiscaletti et al. Reference Fiscaletti, Buxton and Attili2021; Hayashi et al. Reference Hayashi, Watanabe and Nagata2021b). 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. Reference Fiscaletti, Buxton and Attili2021; Hayashi et al. Reference Hayashi, Watanabe and Nagata2021a; Watanabe & Nagata Reference Watanabe and Nagata2022). 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 Reference Horiuti and Takagi2005; Watanabe et al. Reference Watanabe, Tanaka and Nagata2020). 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 Reference Lin and Corcos1984). 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. Reference Patnaik, Sherman and Corcos1976) 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 Reference Lin and Corcos1984). 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 Reference Lin and Corcos1984). 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.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. (Reference Johnsen2010). 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. Reference Hayashi, Watanabe and Nagata2021a). 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 (Reference Yeung and Brasseur1991) 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 (Reference Goto and Vassilicos2016). Another example is a study of a turbulent mixing layer in Takamure et al. (Reference Takamure, Sakai, Ito, Iwano and Hayase2019), 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 (Reference Cimbala, Nagib and Roshko1988), Jiménez & Moin (Reference Jiménez and Moin1991) and Jiménez (Reference Jiménez2018), 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 (Reference Leung, Swaminathan and Davidson2012) and Doan, Swaminathan & Chakraborty (Reference Doan, Swaminathan and Chakraborty2017), 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 Reference Cao, Chen and Doolen1999; Baerenzung et al. Reference Baerenzung, Mininni, Pouquet, Politano and Ponty2010; Goto, Saito & Kawahara Reference Goto, Saito and Kawahara2017; Hirota et al. Reference Hirota, Nishio, Izawa and Fukunishi2017), 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.