1. Introduction
Synchronisation of turbulent flows has been a topic of great interest in recent years. In the simplest form, the technique seeks to reproduce the evolution of one flow in another, by feeding incomplete data obtained from the former into the latter through some coupling mechanism. The technique is at the core of applications such as data assimilation (Kalnay Reference Kalnay2003), where there is a need to enhance the control, construction or prediction of instantaneous turbulence field. Though the interest in the topic has seen resurgence recently, the earliest research can be traced back to two decades ago. Henshaw et al. (Reference Henshaw, Kreiss and Ystróm2003) derive theoretical estimate for the number of Fourier modes required to synchronise two turbulent solutions of Burgers’ equation and those of unforced Navier–Stokes equations (NSEs). Numerical simulations are conducted which verify the theoretical estimate but also show the latter tends to over estimate the bound. Soon after, Yoshida et al. (Reference Yoshida, Yamaguchi and Kaneda2005) investigate the complete synchronisation of two isotropic turbulent flows in three-dimensional (3-D) periodic boxes and find that it is achieved if all Fourier modes with wavenumber less than at least
$k_c$
are coupled, where
$k_c \eta \approx 0.2$
and
$\eta$
is the Kolmogorov length scale (see also Lalescu et al. (Reference Lalescu, Meneveau and Eyink2013)). More recent work on the synchronisation of two direct numerical simulations (DNSs) of isotropic turbulence includes those using the nudging coupling (Di Leoni et al. Reference Di Leoni, Mazzino and Biferale2018, Reference Di Leoni, Mazzino and Biferale2020) and the investigation on partial synchronisation of intense vortices when the number of coupled modes is smaller than the threshold value (Vela-Martin Reference Vela-Martin2021). Nikolaidis & Ioannou (Reference Nikolaidis and Ioannou2022) investigate the synchronisation of Couette flows by coupling selected streamwise Fourier modes. They compute the conditional leading Lyapunov exponents (LLEs) (Boccaletti et al. Reference Boccaletti, Kurths, Osipov, Valladares and Zhou2002) and show that synchronisation happens when the conditional LLE is negative. Channel flows are investigated by Wang & Zaki (Reference Wang and Zaki2022), where the synchronisability of the flows by coupling velocity from different spatial domains is investigated. This work appears to be the first to investigate coupling in physical space. The new scaling relation of the coupling domain is obtained. The behaviour of synchronisation is also investigated when the coupling is applied at the interface between two domains. Inubushi et al. (Reference Inubushi, Saiki, Kobayashi and Goto2023) revisit the synchronisation of turbulence in a periodic box at Reynolds numbers significantly higher than attempted before. Though they also find
$k_c \eta$
to be approximately
$0.2$
based on the results for the conditional LLEs (the transverse Lyapunov exponents in their terminology), it is noted that the conditional LLEs, and hence the exact threshold value for
$k_c \eta$
, depend on the Reynolds number. Wang et al. (Reference Wang, Yuan and Wang2023) investigate the cases when the large-scale data are not sufficient to dominate the small scales. Several techniques are employed to improve the correlation between the small scales when only partial synchronisation is achieved. Using nudging coupling, Buzzicotti & Di Leoni (Reference Buzzicotti and Di Leoni2020) identify the optimal parameters in subgrid-scale (SGS) stress models that minimise the synchronisation error between synchronised large eddy simulation (LES) and DNS. Li et al. (Reference Li, Tian and Li2022) study the synchronisation between LES and DNS using the master–slave coupling, focusing on the threshold wavenumber and the behaviours of different SGS stress models. Most recently, Li et al. (Reference Li, Tian, Li, Si and Mohammed2024) report an investigation into the synchronisation of rotating turbulence in a 3-D periodic box. They show that the forcing scheme has significant impacts on the threshold wavenumber
$k_c$
. More interestingly, they find that the energy spectra of the leading Lyapunov vectors (LLVs) have peaks which always locate approximately at the threshold wavenumber
$k_c$
, regardless the rotation rates and the forcing mechanisms.
The present research is mainly concerned with the last observation mentioned above, namely there appears to be an intrinsic relationship between the peak in the energy spectrum of the LLV and the threshold wavenumber for synchronisation. Though Li et al. (Reference Li, Tian, Li, Si and Mohammed2024) are the first to highlight the link between the two wavenumbers, they are not the first to investigate the energy spectrum of the LLV. Budanur & Kantz (Reference Budanur and Kantz2022) calculate the energy spectrum of the LLV for turbulence in a periodic box for simulations with up to
$432^3$
grid points. The spectrum is found to peak at
$k \eta \approx 0.2$
at all Reynolds numbers. The energy spectrum of small velocity perturbation in box turbulence is investigated by Ge et al. (Reference Ge, Rolland and Vassilicos2023) with up to
$512^3$
grid points. When the perturbation experiences self-similar growth (at which stage it is analogous to the LLV), the peak of its energy spectrum is found to locate at
$kL_\Delta \approx 2$
, where
$L_\Delta$
is the integral length scale of the perturbation. It can be deduced from the data in the paper that this also corresponds to
$k\eta \approx 0.2$
. In other words, the peak wavenumbers found in these two works are both the same as the threshold wavenumber for synchronisation for box turbulence. In rotating turbulence, Li et al. (Reference Li, Tian, Li, Si and Mohammed2024) show that the two wavenumbers display similar dependence on the rotation rate. In particular, for flows with constant power forcing, both decrease as rotation rate is increased. Furthermore, there also seems to be a relationship between the synchronisation threshold and the peaks of the LLV in Couette flow. Nikolaidis & Ioannou (Reference Nikolaidis and Ioannou2022) note in passing that the threshold wavenumber they find coincides with the peak of the spectra of both the vertical and spanwise velocity components at a distance of
$14$
wall units from the wall for the LLV reported by Nikitin (Reference Nikitin2018). The authors stop short of suggesting that there is an intrinsic link between the two parameters.
The evidence summarised above is highly suggestive. However, an analytical theory explaining this relationship is yet to be developed. Therefore, one may justifiably question whether it is only coincidental and further numerical evidence is needed to substantiate the relationship. One obvious approach is to demonstrate this relationship in DNS with much higher Reynolds numbers. However, arguably, the evidence is more convincing if it shows that this relationship also exists in other chaotic systems, thus indicating that the observation in DNS turbulence is an example of a more general principle. Guided by this view, we address in this paper the question of whether the relationship is genuine with numerical experiments for the synchronisation between LES (rather than DNS), and supplement these experiments with synchronisation experiments based on shell models. Specifically, we seek to show that the peak wavenumber of the energy spectrum of the LLV and the threshold wavenumber for synchronisation are the same for both LES turbulence and shell models. We contend that the positive outcome from these experiments will serve to substantiate the authenticity of the relationship in DNS and real turbulence. It is worth emphasising that we do not expect, and do not intend to show, that the threshold or the peak wavenumbers in LES or shell models are the same as their counterparts in the DNS of box turbulence. In essence, LES is not chosen as a model for DNS, but rather a nonlinear dynamical system that is different from DNS. Nevertheless, given the widespread use of LES in turbulence modelling, the findings of this investigation may have yet to be seen implications for turbulent simulations too.
Shell models are a class of reduced-order models which capture important features of turbulence. A detailed discussion about shell models can be found in the monograph by Ditlevsen (Reference Ditlevsen2011) (see also Biferale (Reference Biferale2003) and Bohr et al. (Reference Bohr, Jensen, Paladin and Vulpiani1998)). The shell model we use in the present research is the SABRA model. The state of the model is described by a vector of complex numbers, which mimics the Fourier modes of turbulent velocity. The system evolves according to a set of coupled nonlinear ordinary differential equations. Two copies of the system can be synchronised in a way similar to the synchronisation of two DNSs or two LESs, by coupling some components of the state vectors together during the evolution. Therefore, we may investigate and characterise the synchronisation of the SABRA model in a similar way. However, the simplicity of the SABRA model allows one to simulate a much wider range of length scales that are not reachable even with LES and that is not interfered by the SGS model. These benefits prompt us to supplement LES with simulations of the SABRA model. As we will show later, the SABRA model does display features that are not present in LES.
To examine the relationship between the LLV and the threshold wavenumber more fully, the synchronisation experiments for LES are conducted for three different SGS models: the standard Smagorinsky model, the dynamic Smagorinsky model and the dynamic mixed model. Our hope is to show that the relationship is universally observed for different SGS stress models. However, as will be explained below, the picture is more complicated than this, with the dynamic Smagorinsky model (DSM) exhibiting somewhat different behaviours.
In the course of addressing the main objective laid out above, it becomes clear that two auxiliary questions need to be answered so that we can interpret the main results of the investigation in a meaningful way. The first one is on the interpretation of the LLE of an LES. The LLE of an LES has been calculated previously (Nastac et al. Reference Nastac, Labahn, Magri and Ihme2017; Budanur & Kantz Reference Budanur and Kantz2022). The general observation is that the LLE of an LES is lower than that of the DNS of the same flow and that it decreases as the filter length scale increases. However, the physics behind these results is unclear. We will address this question as the second objective of the current research. We report a scaling law for the LLE of an LES that has not been discussed before and show that the LLE of an LES should be understood as an approximation to the LLE of a filtered DNS. Though this result is auxiliary to our main objective, we consider it a main result of this paper given its foundational nature.
The second auxiliary question is the relationship between the synchronisation threshold and the master flow. This question arises because, in practice, the master system one uses to drive synchronisation is usually different from the slave system. For example, in the studies by both Buzzicotti & Di Leoni (Reference Buzzicotti and Di Leoni2020) and Li et al. (Reference Li, Tian and Li2022), the slave is an LES whereas the master is a DNS. How the threshold and the conditional LLEs depend on the master system is a question one needs to answer to gauge if the results we obtain will be relevant to different practices. We address this question as the third objective of this paper. We do so by comparing synchronisation experiments with different master–slave coupling, which will be described below. We present the results related to this objective in the appendices.
This paper documents our attempt at answering the above questions. We start with a summary of the governing equations in § 2, where we also explain the various coupling used in the synchronisation experiments and the definitions of the conditional LLEs in the context of LES and the SABRA model. The details of the numerical experiments are explained in § 3. The results are presented and discussed in § 4. The results related to the first auxiliary question above are considered fundamental for the interpretation of results related to the LLEs of LES, and therefore are presented first in this section. We then focus on the relationship between the peak wavenumbers for the LLVs and the threshold wavenumber in the synchronisation for LES. The results for the SABRA model are discussed next. The conclusions are summarised in § 5. The results related to the second auxiliary question are summarised in the appendices, as a complement to the main findings of the investigation.
2. Governing equations
2.1. Flow field simulations
We consider incompressible turbulent flows in a 3-D periodic box
$[0, 2\pi ]^3$
. DNS is conducted for flows with relatively low Reynolds numbers, where we integrate the NSEs
\begin{equation} \partial _t {\boldsymbol {u}} + ({\boldsymbol {u}} \cdot \nabla ) {\boldsymbol {u}} = - \nabla p + \nu \nabla ^2 {\boldsymbol {u}} + {\boldsymbol {f}}, \end{equation}
in which
${\boldsymbol {u}} \equiv (u_1, u_2, u_3)$
is the velocity field,
$p$
is the pressure divided by the constant density,
$\nu$
is the viscosity and
$\boldsymbol {f}$
is the forcing term. The velocity is assumed to be incompressible so that
\begin{equation} \nabla \cdot {\boldsymbol {u}} = 0. \end{equation}
LESs are conducted for flows with high Reynolds numbers, which are based on the filtered NSEs (fNSEs)
\begin{equation} \partial _t \overline {\boldsymbol {u}} + (\overline {\boldsymbol {u}} \cdot \nabla ) \overline {\boldsymbol {u}} = - \nabla {\overline {p}}+ \nabla \cdot (-\boldsymbol {\tau }) + \nu \nabla ^2 \overline {\boldsymbol {u}} + \overline {\boldsymbol {f}} + {\boldsymbol {f}}_s, \end{equation}
and the filtered continuity equation
$\nabla \cdot \overline {\boldsymbol {u}} = 0$
. In these equations, the line
$\overline {\phantom {a}}$
represents filtering with filter length
$\Delta$
. Therefore,
$\overline {\boldsymbol {u}}$
and
$\overline {p}$
are the filtered velocity and pressure, respectively. The SGS stress tensor
$\boldsymbol {\tau }$
in (2.3) is defined as
$\tau _{ij} = \overline {u_iu}_j - \overline {u}_i \overline {u}_j,$
whereas
$\overline {\boldsymbol {f}}$
is the filtered forcing term and
${\boldsymbol {f}}_s$
symbolically represents the effects of coupling in synchronisation experiments, which will be further explained below.
We let
${\boldsymbol {f}} = (a_f\cos k_f x_2,0,0)$
with
$a_f = 0.15$
and
$k_f = 1$
in all simulations. The flow resulting from this type of deterministic sinusoidal forcing is often called the Kolmogorov flow (Borue & Orszag Reference Borue and Orszag1996). The flow is inhomogeneous with a sinusoidal mean velocity profile. Previous research (Yoshida et al. Reference Yoshida, Yamaguchi and Kaneda2005; Inubushi et al. Reference Inubushi, Saiki, Kobayashi and Goto2023; Li et al. Reference Li, Tian, Li, Si and Mohammed2024) shows that the forcing only has at most secondary effects on the synchronisability of flows in periodic boxes for non-rotating turbulence. Therefore, the conclusions we draw can be compared with previous research where different forcing is used.
Simulations with different SGS stress models are conducted to assess the reproducibility of key observations. Three canonical models are considered, namely the standard Smagorinsky model (SSM), the dynamic Smagorinsky model (DSM) and the dynamic mixed model (DMM). Details of these models can be found from, e.g., Meneveau & Katz (Reference Meneveau and Katz2000). We summarise the formulae here for references. The SSM is defined with
\begin{equation} \tau _{ij} = - 2(c_s \Delta )^2 \vert \overline {\boldsymbol {s}}\vert {\overline {s}}_{ij}, \end{equation}
where
$c_s = 0.16$
is the Smagorinsky coefficient,
${\overline {s}}_{ij}$
is the filtered strain rate tensor given by
\begin{equation} {\overline {s}}_{ij} = \frac {1}{2}(\partial _j {\overline {u}}_i + \partial _i {\overline {u}}_j) \end{equation}
and
$\vert \overline {\boldsymbol {s}}\vert \equiv (2 {\overline {s}}_{ij}{\overline {s}}_{ij})^{1/2}$
is the magnitude of the filtered strain rate. The DSM is also defined by (2.4), but the coefficient
$c_s^2$
is calculated dynamically using the Germano identity (Germano Reference Germano1992). Following the procedure as described by, e.g., Germano (Reference Germano1992) and Meneveau & Katz (Reference Meneveau and Katz2000), one can derive
\begin{equation} c_s^2 = \frac {\langle L_{ij} M_{ij}\rangle _v}{\langle M_{ij} M_{ij}\rangle _v}, \end{equation}
with
\begin{equation} M_{ij} = -2 \Delta ^2 \left [4 \vert \widetilde {\overline {\boldsymbol {s}}} \vert \widetilde {{\overline {s}}}_{ij} - \widetilde {\vert \overline {\boldsymbol {s}}\vert {\overline {s}}_{ij}}\right ], \quad L_{ij} = \widetilde {{\overline {u}}_i {\overline {u}}_j} - \widetilde {{\overline {u}}}_i \widetilde {{\overline {u}}}_j, \end{equation}
where
$\langle \phantom {\,}\rangle _v$
represents volume average and
$\widetilde {\phantom {u}}$
denotes test-filtering with filter scale
$2\Delta$
using the cut-off filter (Pope Reference Pope2000). The model expression for DMM is given by
\begin{equation} \tau _{ij} = - 2(c_s \Delta )^2 \vert \overline {\boldsymbol {s}}\vert {\overline {s}}_{ij} + c_{nl}\Delta ^2 \overline {A}_{ik} \overline {A}_{jk} \,, \end{equation}
where
$\overline {A}_{ij} \equiv \partial _j \overline {u}_i$
is the filtered velocity gradient and
$c_{nl}$
is the nonlinear coefficient, which is to be determined dynamically. Following the dynamic procedure, the expressions for
$c_s^2$
and
$c_{nl}$
are found to be
\begin{align} c_s^2 = \frac {\langle L_{ij} M_{ij} \rangle _v \langle N_{ij} N_{ij} \rangle _v - \langle L_{ij} N_{ij} \rangle _v \langle M_{ij} N_{ij}\rangle _v}{\langle M_{ij} M_{ij} \rangle _v \langle N_{ij} N_{ij} \rangle _v - \langle M_{ij} N_{ij} \rangle _v^2 }, \end{align}
\begin{align} c_{nl} = \frac {\langle L_{ij} N_{ij} \rangle _v \langle M_{ij} M_{ij} \rangle _v - \langle L_{ij} M_{ij} \rangle _v \langle M_{ij} N_{ij}\rangle _v}{\langle M_{ij} M_{ij} \rangle _v \langle N_{ij} N_{ij} \rangle _v - \langle M_{ij} N_{ij} \rangle _v^2 }, \end{align}
where
\begin{equation} N_{ij} = \Delta ^2 \left [4 \widetilde {\overline {A}}_{ik}\widetilde {\overline {A}}_{jk} - \widetilde {\overline {A}_{ik} \overline {A}}_{jk}\right ]. \end{equation}
The three models are compared in depth by Kang et al. (Reference Kang, Chester and Meneveau2003) with experimental data. DMM produces better agreement in most statistics and correlates better with the real SGS stress. SSM is known to be too dissipative. DSM has the same expression as SSM, but
$c_s^2$
, obtained dynamically, fluctuates slightly over time and, in general, is somewhat smaller. As a consequence, DSM is less dissipative compared with SSM, and, in general, agrees with experimental data better than SSM.
Figure 1. Set-up of the three-body experiments and the twin experiments showing differing coupling mechanisms. The master in the former is a DNS in Group I and an LES in Group II. It may have a different number of slave modes from the slave systems, and its slave modes do not converge to those in the slave systems
$A$
and
$B$
.
2.2. Coupling between the masters and slaves
We run two groups of synchronisation experiments, and the experiments in each group are categorised into two sets with differing coupling methods. The set-up of the experiments is illustrated in figure 1. Group I consists of numerical experiments with lower Reynolds numbers so that the results can be compared with DNS. Group II involves simulations with higher Reynolds numbers that are not reachable by DNS. Therefore, DNS is replaced by LES in Group II. In the first set of experiments, a master system
$M$
and two slave systems
$A$
and
$B$
are simulated concurrently. We thus call them three-body experiments. In three-body experiments, systems
$A$
and
$B$
are identical LES with the same SGS stress model, but they have different initial conditions. The two systems are synchronised indirectly through their coupling with the same master
$M$
. System
$M$
is a DNS in Group I, but it is an LES in Group II that has a SGS stress model different from that in systems
$A$
and
$B$
. In the second set of experiments, systems
$A$
and
$B$
are coupled directly together, with
$A$
now serving as the master. These experiments are called twin experiments.
The three-body experiments are similar to those reported by Li et al. (Reference Li, Tian and Li2022). However, for such three-body experiments, the properties of the conditional LLEs and the LLVs, and their relationships with the synchronisation threshold, have not been investigated, which are addressed here. The twin experiments are similar to those reported by, e.g., Yoshida et al. (Reference Yoshida, Yamaguchi and Kaneda2005), Inubushi et al. (Reference Inubushi, Saiki, Kobayashi and Goto2023) and Li et al. (Reference Li, Tian, Li, Si and Mohammed2024), in the sense that they all deal with synchronisation between identical systems. However, here, the systems are LES as opposed to DNS. How the conditional LLEs and the LLVs behave differently is the question that motivates this set of experiments. Furthermore, the contrast between the two sets of experiments allows us to assess the impacts of the coupling mechanism: indirect coupling via a third (different) master system versus direct coupling. Note that the results addressing the questions outlined above are mostly summarised in the appendices. In the main text, the discussion centres around the relationship between the threshold wavenumber and the LLV, using the data gathered from both the three-body and the twin experiments.
We now provide the mathematical representation of the coupling and introduce several definitions. We use
${\boldsymbol {u}}_{M}({\boldsymbol {x}},t)$
to denote the velocity of system
$M$
and
$\hat {{\boldsymbol {u}}}_{M}({\boldsymbol {k}},t)$
its Fourier mode, where
$\boldsymbol {k}$
is the wavenumber. Similarly, the duos
$\overline {\boldsymbol {u}}_{A}$
and
$\hat {\overline {\boldsymbol {u}}}_{A}$
, and
$\overline {\boldsymbol {u}}_{B}$
and
$\hat {\overline {\boldsymbol {u}}}_{B}$
are used to denote the velocity and its Fourier modes for systems
$A$
and
$B$
, respectively. We use
$\vert {\boldsymbol {k}}\vert$
to denote the Euclidean norm of
$\boldsymbol {k}$
, i.e.
$\vert {\boldsymbol {k}} \vert \equiv (k_1^2 + k_2^2 + k_3^2)^{1/2}$
.
Let
$k_m$
be a constant and suppose we intend to couple the Fourier modes with wavenumbers
$\vert {\boldsymbol {k}} \vert \leqslant k_m$
. In the three-body experiments, the coupling is implemented by letting
\begin{equation} \hat {\overline {\boldsymbol {u}}}_{A}({\boldsymbol {k}},t) = \hat {{\boldsymbol {u}}}_{M}({\boldsymbol {k}},t), \qquad \hat {\overline {\boldsymbol {u}}}_{B}({\boldsymbol {k}},t) = \hat {{\boldsymbol {u}}}_{M}({\boldsymbol {k}},t), \end{equation}
for
$\vert {\boldsymbol {k}} \vert \leqslant k_m$
at each time step. The effects of (2.12) are represented by
${\boldsymbol {f}}_s$
in (2.3) for both systems
$A$
and
$B$
. In the twin experiments, the coupling is implemented by letting
\begin{equation} \hat {\overline {\boldsymbol {u}}}_{B}({\boldsymbol {k}},t) = \hat {\overline {\boldsymbol {u}}}_{A}({\boldsymbol {k}},t) \end{equation}
for
$\vert {\boldsymbol {k}} \vert \leqslant k_m$
. This coupling introduces the forcing
${\boldsymbol {f}}_s$
into (2.3) for system
$B$
, but does not impact system
$A$
.
The constant
$k_m$
is called the coupling wavenumber. The Fourier modes with
$\vert {\boldsymbol {k}} \vert \leqslant k_m$
are called the master modes, whereas those with
$\vert {\boldsymbol {k}} \vert \gt k_m$
are called the slave modes. Systems
$A$
and
$B$
synchronise asymptotically when
$k_m$
is sufficiently large, i.e. the difference between
${\boldsymbol {u}}^A$
and
${\boldsymbol {u}}^B$
decays towards zero over time. The smallest
$k_m$
value for which such synchronisation happens is the threshold wavenumber, denoted as
$k_c$
. The wavenumber
$k_c$
might be different for different experiments. The dependence of
$k_c$
on the coupling, flow parameters and its relationship with the LLV is the main focus of this paper.
2.3. Lyapunov exponents and Lyapunov vectors
The synchronisation experiments are analysed via the LLEs (i.e. the leading Lyapunov exponents), the conditional LLEs and the LLVs (i.e. the leading Lyapunov vectors) of the LES of the flows. We introduce the relevant definitions in this subsection.
As we have mentioned in § 1, a physical interpretation of the LLE of an LES is lacking, and we will show that it should be interpreted as an approximation to the LLE of filtered DNS of the same flow. Therefore, we start with the definition of the LLE of filtered DNS. Let
$\boldsymbol {u}$
be a DNS velocity and
${\boldsymbol {u}}_\Delta ^\delta$
be the filtered infinitesimal perturbation, which is zero for wavenumber larger than the cut-off wavenumber
$k_\Delta \equiv \pi /\Delta$
, where
$\Delta$
is the filter length. The Fourier mode of
${\boldsymbol {u}}_\Delta ^\delta$
is denoted by
$\hat {{\boldsymbol {u}}}_\Delta ^\delta$
. By definition,
\begin{equation} \hat {{\boldsymbol {u}}}_\Delta ^\delta ({\boldsymbol {k}}, t) = 0 \end{equation}
for
$\vert {\boldsymbol {k}} \vert \gt k_\Delta$
. The perturbation
${\boldsymbol {u}}_\Delta ^\delta$
is divergence free and obeys the linearised NSE:
\begin{equation} \partial _t {\boldsymbol {u}}_\Delta ^\delta + ({\boldsymbol {u}} \cdot \nabla ) {\boldsymbol {u}}_\Delta ^\delta + ({\boldsymbol {u}}_\Delta ^\delta \cdot \nabla ) {\boldsymbol {u}} = - \nabla p^\delta + \nu \nabla ^2 {\boldsymbol {u}}_\Delta ^\delta + {\boldsymbol {f}}^\delta , \end{equation}
where
$p^\delta$
is the pressure perturbation and
${\boldsymbol {f}}^\delta$
is the perturbation to the forcing term. The LLE for the filtered DNS, denoted by
$\lambda _\Delta$
, is defined as
\begin{equation} \lambda _\Delta = \lim _{t\to \infty } \frac {1}{t}\log \frac {\Vert {\boldsymbol {u}}_\Delta ^\delta ({\boldsymbol {x}}, t+t_0)\Vert }{\Vert {\boldsymbol {u}}_\Delta ^\delta ({\boldsymbol {x}}, t_0)\Vert }, \end{equation}
where
$t_0$
is an arbitrary initial time and
$\Vert \cdot \Vert$
represents the 2-norm defined by (for a generic vector field
$\boldsymbol {w}$
)
\begin{equation} \Vert \boldsymbol {w} \Vert \equiv \langle \boldsymbol {w} \cdot \boldsymbol {w} \rangle _v^{1/2}. \end{equation}
Obviously,
$\lambda _\Delta$
depends on the filter scale
$\Delta$
. It measures the average growth rate of the perturbations only applied to the Fourier modes of
$\boldsymbol {u}$
with length scales larger than
$\Delta$
.
The conditional LLE of a slave system in a synchronisation experiment is defined in a similar way. The conditional LLE measures the averaged growth rate of an infinitesimal perturbation to the slave modes of a slave system. Let
$\overline {\boldsymbol {u}}$
be the velocity of a slave, i.e. systems
$A$
or
$B$
in a three-body experiment or system
$B$
in a twin experiment. We let
$\overline {\boldsymbol {u}}^\delta$
be an infinitesimal perturbation to the slave modes of
$\overline {\boldsymbol {u}}$
. As the master modes are not perturbed, we have
\begin{equation} \hat {\overline {\boldsymbol {u}}}^\delta ({\boldsymbol {k}},t) = 0 \quad \text {for}\ \vert {\boldsymbol {k}} \vert \leqslant k_m. \end{equation}
By definition, the perturbation
$\overline {\boldsymbol {u}}^\delta$
is solenoidal and obeys the linearised fNSE
\begin{equation} \partial _t \overline {\boldsymbol {u}}^\delta + (\overline {\boldsymbol {u}} \cdot \nabla ) \overline {\boldsymbol {u}}^\delta + (\overline {\boldsymbol {u}}^\delta \cdot \nabla ) \overline {\boldsymbol {u}} = - \nabla {\overline {p}}^\delta + \nabla \cdot (-\boldsymbol {\tau }^\delta ) + \nu \nabla ^2 \overline {\boldsymbol {u}}^\delta + \overline {\boldsymbol {f}}^\delta , \end{equation}
where
${\overline {p}}^\delta$
,
$\overline {\boldsymbol {f}}^\delta$
and
$\boldsymbol {\tau }^\delta$
are the pressure perturbation, the perturbation to the filtered forcing and the perturbation to the SGS stress
$\boldsymbol {\tau }$
, respectively. Here,
$\overline {\boldsymbol {f}}^\delta$
is included for completeness, but it is actually always zero for the forcing we are using.
The conditional LLE for the slave is defined in terms of
$\overline {\boldsymbol {u}}^\delta$
in a way similar to (2.16). Let
$\lambda (k_m)$
denote the conditional LLE with coupling wavenumber
$k_m$
. Then,
$\lambda (k_m)$
is defined by (Boccaletti et al. Reference Boccaletti, Kurths, Osipov, Valladares and Zhou2002; Nikolaidis & Ioannou Reference Nikolaidis and Ioannou2022; Inubushi et al. Reference Inubushi, Saiki, Kobayashi and Goto2023)
\begin{equation} \lambda (k_m) = \lim _{t\to \infty } \frac {1}{t}\log \frac {\Vert \overline {\boldsymbol {u}}^\delta ({\boldsymbol {x}}, t+t_0)\Vert }{\Vert \overline {\boldsymbol {u}}^\delta ({\boldsymbol {x}}, t_0)\Vert }. \end{equation}
The conditional LLE
$\lambda (k_m)$
is a function of
$k_m$
. It becomes the usual (unconditional) LLE when
$k_m=0$
. Synchronisation is successful only when the conditional LLE is negative, as is shown previously in twin experiments for turbulent channel flows (Nikolaidis & Ioannou Reference Nikolaidis and Ioannou2022), turbulence in periodic boxes (Inubushi et al. Reference Inubushi, Saiki, Kobayashi and Goto2023) and rotating turbulence (Li et al. Reference Li, Tian, Li, Si and Mohammed2024). As we show in the appendices, the observation also holds for the three-body experiments. Therefore, the threshold wavenumber
$k_c$
for
$k_m$
can be found from
$\lambda (k_m)$
as the root of the equation
$\lambda (k_m) = 0$
.
For sufficiently large
$t$
, the perturbation velocity
$\overline {\boldsymbol {u}}^\delta$
approaches the conditional LLV or the (unconditional) LLV if
$k_m=0$
(Bohr et al. Reference Bohr, Jensen, Paladin and Vulpiani1998; Kuptsov & Parlitz Reference Kuptsov and Parlitz2012; Nikitin Reference Nikitin2018). The conditional LLV represents the most unstable perturbation to the slave modes on long-time average. The properties of the (unconditional) LLV for a turbulent channel flow are analysed by Nikitin (Reference Nikitin2018). In our previous work (Li et al. Reference Li, Tian, Li, Si and Mohammed2024), we find that the threshold coupling wavenumber
$k_c$
correlates closely with the peak of the mean energy spectrum of the LLVs. One of the objectives of this work is to look into the properties of the LLVs in relation to the synchronisability of LES.
2.4. SABRA model
The investigation based on LES is supplemented with experiments conducted with the SABRA model. The state of the SABRA model is a vector
$\boldsymbol {z}\equiv (u_1, u_2, \ldots , u_{N_{SB}})$
, where
$u_n$
is a complex number referred to as the
$n$
th shell velocity and
$N_{SB}$
is the number of shells in the model. The equation for
$u_n$
(
$n=1,2,\ldots , N_{SB}$
) reads
\begin{equation} \frac {{\rm d}u_n}{{\rm d}t} + \nu k_n^2 u_n = i(k_n u_{n+2} u^*_{n+1} - b k_{n-1} u_{n+1} u^*_{n-1} - c k_{n-2} u_{n-1} u_{n-2}) + f \delta _{n1}, \end{equation}
where
$\nu$
is the viscosity,
$i$
is the imaginary unit,
$k_n$
is referred to as the wavenumber for the
$n$
th shell,
$^*$
denotes complex conjugate,
$b$
and
$c\equiv b-1$
are constant parameters,
$f$
is the forcing amplitude, and
$\delta _{n1}$
is the Kronecker delta
$\delta _{ij}$
with
$(i,j)=(n, 1)$
, which indicates that the forcing is applied only to the first shell.
We use the following standard values for the parameters. The wavenumber
$k_n$
is given by
$ k_n= k_0 q^n$
with
$k_0 = 0.05$
and
$q = 2$
. The forcing amplitude
$f = (0.005, 0.005)$
. We let
$b = 1/2$
so that the system has a helicity-like inviscid invariant. This leads to
$c = -1/2$
due to the relation
$c = b-1$
, which guarantees that the total energy of the system is also an inviscid invariant of the model. The viscosity
$\nu$
determines the Reynolds number, which will be adjusted in the synchronisation experiments.
For suitable parameters, the evolution of the SABRA model becomes chaotic. In this regime, the energy spectrum of the system displays an inertial range where a
$k_n^{-5/3}$
distribution is observed. That is, the slope of the spectrum is the same as that for isotropic turbulence. The energy flux also becomes a constant across the inertial range (Ditlevsen Reference Ditlevsen2011). Some of these features will be shown below. These and many other properties make the SABRA model a useful simplified model for real turbulence, which is also the motivation for us to choose it for our synchronisation experiments.
We only conduct twin experiments using the SABRA model, and only use it to investigate the relationship between the peak wavenumber of its LLV and the threshold wavenumber for synchronisation. The set-up of these twin experiments is as illustrated in figure 1, except that systems
$A$
and
$B$
are now both the SABRA model, with different initial states. The states of the two systems are denoted by
$\boldsymbol {z}_A$
and
$\boldsymbol {z}_B$
. The shell velocities
$u_n$
with wavenumber
$k_n\leqslant k_m$
are the master modes. The master modes in
$\boldsymbol {z}_A$
are copied into
$\boldsymbol {z}_B$
and replace those in
$\boldsymbol {z}_B$
at every time step in the twin experiments (cf. (2.13)). For notational simplicity, the conditional LLE for the SABRA model with coupling wavenumber
$k_m$
is also denoted by
$\lambda (k_m)$
, and it is defined by (2.20) with
$\overline {\boldsymbol {u}}^\delta ({\boldsymbol {x}}, t)$
replaced by
$\boldsymbol {z}^\delta (t)$
, where
$\boldsymbol {z}^\delta (t)$
is the infinitesimal perturbation to
$\boldsymbol {z}_A(t)$
. The limit of
$\boldsymbol {z}^\delta (t)$
as
$t\to \infty$
is the LLV for the SABRA model.
3. Numerical experiments
3.1. Summary of the numerical experiments
Tables 1 and 2 summarise the key parameters for the LES in Group I and II, respectively, with the definitions of the parameters given below in § 3.2. We label a case with a code that consists of three segments. The first segment from the left is in the form of ‘R
$a$
’ with
$a$
being an integer between
$1$
and
$5$
. Different
$a$
corresponds to different kinetic viscosity
$\nu$
. The middle segment is the acronym of the SGS stress model. The third segment is the value of
$N_{LES}$
, with
$N_{LES}^3$
being the number of grid points of the simulation. Each case (i.e. each code) represents a series of simulations with different coupling wavenumber
$k_m$
. Where necessary, we append ‘K
$b$
’ to the end of the code to distinguish such simulations, with
$b$
being the value of
$k_m$
. Also, the cases with codes starting with the same ‘R
$a$
’ are sometimes collectively referred to as ‘subgroup R
$a$
’. The cases in such a subgroup may have different SGS stress models,
$N_{LES}$
or
$k_m$
.
Table 1. Parameters for low-Reynolds-number LESs in Group I.
Table 2. Parameters for the high-Reynolds-number LESs in Group II.
The parameters for the DNS in the three-body experiments are summarised in table 3. They correspond to the LES in table 1. We use ‘R
$a$
DNS’ with
$N_{DNS}$
appended to refer to the DNS in subgroup R
$a$
,
$N_{DNS}^3$
being the number of grid points for the DNS. Note that each subgroup R
$a$
is associated with at most a single DNS (there is no DNS for subgroups R4 and R5). In total, we conduct approximately 300 simulations where two or three velocity fields are solved concurrently (cf. § 3.3), with up to
$256^3$
grid points for each field.
Table 3. Parameters for the DNSs corresponding to the LESs in Group I.
The parameters for the twin experiments for the SABRA model are summarised in table 4, with the definitions of the parameters given in § 3.2. Ten cases with different Reynolds numbers are computed and they are named cases
$Re_0$
,
$Re_1$
, …,
$Re_{9}$
, respectively.
Table 4. Parameters for the SABRA model.
3.2. Definitions of main parameters
We now summarise briefly the main parameters of the simulations. We introduce the energy spectrum first. Using
$\boldsymbol {v}({\boldsymbol {x}},t)$
and
$\hat {\boldsymbol {v}}({\boldsymbol {k}},t)$
to represent a generic DNS or LES velocity and its Fourier transform, respectively, the energy spectrum of
$\boldsymbol {v}$
is denoted by
$E(k)$
with
\begin{equation} E(k) = \frac {1}{2}\sum _{k \leqslant \vert {\boldsymbol {k}} \vert \leqslant k+1} \langle \hat {\boldsymbol {v}}({\boldsymbol {k}},t) \cdot \hat {\boldsymbol {v}}^*({\boldsymbol {k}},t)\rangle, \end{equation}
where
$\langle \, \rangle$
represents ensemble average, i.e. average over space and time as well as different realisations where applicable. The turbulent kinetic energy
$K$
and the viscous energy dissipation rate
$\epsilon$
then can be calculated from
$E(k)$
with
\begin{equation} K = \int _0^\infty E(k) {\rm d}k, \quad {\epsilon } = 2 \nu \int _0^\infty k^2 E(k){\rm d}k. \end{equation}
The above equations apply to both LES and DNS, with the understanding that
$E(k) = 0$
when
$k\ge k_\Delta$
for LES.
For DNS,
$\epsilon$
is a key parameter as it represents the mean energy flux across the inertial range in stationary turbulence (Pope Reference Pope2000). For LES, however,
$\epsilon$
is usually not pertinent, because the energy flux is the sum of
$\epsilon$
and the SGS energy dissipation
$\Pi _\Delta \equiv -\langle \tau _{ij} {\overline {s}}_{ij}\rangle$
and
$\epsilon$
is usually much smaller than
$\Pi _\Delta$
(Meneveau & Katz Reference Meneveau and Katz2000; Pope Reference Pope2000). We refer to the sum as the total dissipation rate and denote it by
${\epsilon }_t$
, i.e.
\begin{equation} {\epsilon }_t = {\epsilon } + \Pi _\Delta = {\epsilon } - \langle \tau _{ij} {\overline {s}}_{ij}\rangle . \end{equation}
Here,
${\epsilon }_t$
is the suitable parameter for the characterisation of LES velocity and it can be calculated from the LES velocity with
$\tau _{ij}$
given by a SGS stress model. As
${\epsilon }_t = {\epsilon }$
in DNS, we can use
${\epsilon }_t$
uniformly to introduce relevant length and time scales, and other parameters that we will use to normalise the conditional LLEs and the threshold wavenumbers.
We thus define the root-mean-square (r.m.s.) velocity
$u_{\textit{rms}}$
, the Taylor length scale
$\lambda _a$
and the Taylor scale Reynolds number via
\begin{equation} u_{\textit{rms}} = (2K/3)^{1/2}, \quad \lambda _a = (15 \nu u^2_{\textit{rms}}/{\epsilon }_t)^{1/2}, \quad Re_\lambda = \frac {u_{\textit{rms}} \lambda _a }{\nu }, \end{equation}
and the Kolmogorov length scale
$\eta$
and time scale
$\tau _k$
by
\begin{equation} \eta = (\nu ^3/{\epsilon }_t)^{1/4}, \quad \tau _k = (\nu /{\epsilon }_t)^{1/2}. \end{equation}
The eddy turnover time at the filter scale
$\Delta$
, denoted by
$\tau _\Delta$
, is defined as
\begin{equation} \tau _\Delta = \Delta ^{2/3}{\epsilon }_t^{-1/3}. \end{equation}
For DNS, the definitions introduced in (3.4) and (3.5) are the same as the usual definitions. For LES,
$\lambda _a$
,
$\eta$
and
$\tau _k$
do not usually represent real scales in the velocity field. Rather they should be understood as approximations to those of the true velocity which the LES is simulating. Though LES results are very often normalised by parameters at the filter scale (e.g.
$\tau _\Delta$
), we also use
$\eta$
and
$\tau _k$
defined above to compare the results with those in the literature and to present the results in a unified way at both low and high Reynolds numbers.
As will be shown later, it is instructive to correlate the LES results (in particular, the threshold wavenumber
$k_c$
) using an effective dissipation length scale
$\Delta _t$
. We first define a (constant) effective viscosity
$\nu _t$
by
\begin{equation} \nu _t \equiv \frac { {\epsilon }_t}{2 \langle {\overline {s}}_{ij} {\overline {s}}_{ij}\rangle }. \end{equation}
Since the above equation implies
\begin{equation} {\epsilon }_t = 2 \nu _t \langle {\overline {s}}_{ij}{\overline {s}}_{ij}\rangle , \end{equation}
it is clear that
$\nu _t$
is defined in such a way that the total energy dissipation can be calculated from the eddy viscosity model
$\tau _{ij}= -2 \nu _t {\overline {s}}_{ij}$
. The effective dissipation length scale
$\Delta _t$
is then defined as
\begin{equation} \Delta _t = (\nu _t^3/{\epsilon }_t)^{1/4}. \end{equation}
According to the Kolmogorov phenomenology,
$\Delta _t$
should scale with
$\Delta$
at high Reynolds numbers for
$\Delta$
in the inertial range. Nevertheless, our numerical experiments later show that
$\Delta _t$
correlates better with data than
$\Delta$
.
For the SABRA model, similar quantities can be defined. The turbulent kinetic energy
$K$
is defined by (Biferale Reference Biferale2003; Ditlevsen Reference Ditlevsen2011)
\begin{equation} K = \frac {1}{2} \sum _{n=1}^{N_{SB}} \langle \vert u_n \vert ^2 \rangle , \end{equation}
whereas the viscous energy dissipation rate
$\epsilon$
is defined by
\begin{equation} {\epsilon } = \nu \sum _{n=1}^{N_{SB}} k_n^2\langle \vert u_n \vert ^2\rangle . \end{equation}
We then define
$u_{\textit{rms}}$
using (3.4), and
$\eta$
and
$\tau _k$
using (3.5), with the understanding that
${\epsilon }_t = {\epsilon }$
for the SABRA model. The Taylor scale Reynolds number is not often used in the SABRA model. Instead, we introduce a crude ‘integral’ length scale
$L \equiv \pi /k_1 \approx 31.42$
, where
$k_1$
is the wavenumber for the first shell, and define the Reynolds number as
\begin{equation} Re = \frac {u_{\textit{rms}} L}{\nu }. \end{equation}
The values of these parameters for the LES are summarised in tables 1 and 2, those for DNS in table 3, and the SABRA model in table 4.
3.3. Numerical methods
Information for the numerical schemes is given in what follows. The NSE and the fNSE (see (2.1) and (2.3)), together with the continuity equations, are integrated with the pseudo-spectral method. The domain is discretised uniformly with
$N^3_{DNS}$
and
$N^3_{LES}$
grid points for DNS and LES, respectively. The two-thirds rule (Pope Reference Pope2000) is applied to de-alias the advection term so that the maximum effective wavenumber is
$N_{DNS}/3$
for DNS and
$N_{LES}/3$
for LES. Time stepping uses an explicit second-order Euler scheme (Li et al. Reference Li, Zhang, Dong and Abdullah2020). The viscous diffusion term is treated with an integration factor prior to time discretisation. The step-size
$\delta t$
is chosen in such a way that the Courant–Friedrichs–Lewy number
$ u_{\textit{rms}} \delta t/\delta x$
is less than
$0.1$
, where
$\delta x$
is the grid size.
The conditional LLEs are calculated using the algorithm reported by, e.g., Boffetta & Musacchio (Reference Boffetta and Musacchio2017) and Budanur & Kantz (Reference Budanur and Kantz2022), with some minor modification. We run concurrent simulations for
$\overline {\boldsymbol {u}}_{A}$
and
$\overline {\boldsymbol {u}}_{B}$
in twin experiments, and for
$\boldsymbol {u}_M$
,
$\overline {\boldsymbol {u}}_A$
and
$\overline {\boldsymbol {u}}_B$
in three-body experiments. Here,
$\overline {\boldsymbol {u}}_A$
and
$\overline {\boldsymbol {u}}_B$
initially differ by a small perturbation such that
\begin{equation} e \equiv \Vert \overline {\boldsymbol {u}}_A({\boldsymbol {x}},0) - \overline {\boldsymbol {u}}_B({\boldsymbol {x}},0)\Vert \end{equation}
is small, and the perturbation is only introduced to the slave modes in system
$B$
. During the simulation,
$\overline {\boldsymbol {u}}_{B}$
is periodically re-initialised at
$t_n = n \triangle t$
(
$n=1,2,\ldots$
) for some short time interval
$\triangle t ( \ne \delta t)$
so that
\begin{equation} \overline {\boldsymbol {u}}_{B}({\boldsymbol {x}}, t_n^+) = \overline {\boldsymbol {u}}_{A}({\boldsymbol {x}}, t_n) + g_n^{-1} \left [\overline {\boldsymbol {u}}_{B}({\boldsymbol {x}}, t_n^-) - \overline {\boldsymbol {u}}_{A}({\boldsymbol {x}}, t_n)\right ], \end{equation}
where
$t_n^\pm$
are the times immediately after/before
$t_n$
and
\begin{equation} g_n = e^{-1}\Vert \overline {\boldsymbol {u}}_{B}({\boldsymbol {x}}, t_n^-)- \overline {\boldsymbol {u}}_{A}({\boldsymbol {x}}, t_n)\Vert \end{equation}
is the amplification factor for the perturbation. The conditional LLE is given by
\begin{equation} \lambda (k_m) \approx \lim _{M\to \infty } \frac {1}{M \triangle t} \sum _{n=1}^M \log g_n. \end{equation}
In our simulations,
$\triangle t \approx 0.1 \tau _k$
for cases in Group I. It varies somewhat for cases in Group II, but is never larger than
$0.05 \tau _\Delta$
.
Figure 2. Energy spectra: (a) for cases in R3; (b) for cases in R5. Dash-dotted lines without symbols, the
$k^{-5/3}$
power law; the dashed line without symbol, DNS for R3.
Note that because of the coupling between
$\overline {\boldsymbol {u}}_A$
and
$\overline {\boldsymbol {u}}_B$
(cf. (2.12) and (2.13)), the perturbation
$\overline {\boldsymbol {u}}_{B}-\overline {\boldsymbol {u}}_{A}$
is always zero for Fourier modes with wavenumber
$\vert {\boldsymbol {k}} \vert \leqslant k_m$
, as required by (2.18). The LLE for filtered DNS, namely
$\lambda _\Delta$
, is calculated in much the same way. Here, we run two uncoupled DNS for
${\boldsymbol {u}}_A$
and
${\boldsymbol {u}}_B$
. When
${\boldsymbol {u}}_B$
is re-initialised, the perturbation in Fourier modes with
$\vert {\boldsymbol {k}} \vert \gt k_\Delta$
is set to zero by letting
\begin{equation} \hat {{\boldsymbol {u}}}_B({\boldsymbol {k}}, t) = \hat {{\boldsymbol {u}}}_A({\boldsymbol {k}},t) \quad (\vert {\boldsymbol {k}} \vert \gt k_\Delta ), \end{equation}
hence enforcing (2.14).
The SABRA model is integrated in time with an explicit fourth-order Runge–Kutta method with step-size
$\delta t_{SB} = 10^{-5}$
. The conditional LLEs for the SABRA model are calculated using the same algorithm described above, with
$\overline {\boldsymbol {u}}_A$
and
$\overline {\boldsymbol {u}}_B$
replaced by respective copies of the SABRA model.
4. Results and analyses
We begin this section with a preamble to present several supplementary results. First, we present the energy spectra for selected cases as a validation of the numerics. Figure 2 shows the energy spectra, with low-Reynolds-number results in panel (a) and high-Reynolds-number results in panel (b). The DNS spectrum included in figure 2(a) displays features that are consistent with the
$k^{-5/3}$
Kolmogorov spectrum. The results from LES reproduce the main features of the DNS spectrum. DSM and DMM, in particular, appear to capture the tail of the spectrum better. No DNS data are available for comparison in the high-Reynolds-number cases, but an extended range with a slope of
$-5/3$
is clearly visible in the spectra, reproducing the expected behaviours from these canonical SGS stress models.
As explained previously, several questions auxiliary to our main objective are addressed in this investigation. One of them is how different coupling methods impact the synchronisation process. Specifically, one may ask whether the synchronisation threshold is independent of the coupling mechanism, or more pointedly, whether synchronisation in three-body experiments can be described by the conditional LLEs of the slave system at all, noting that the actual measure of synchronisation is the synchronisation error, not the conditional LLE. To keep the focus of the paper clear, we present the discussion of these questions in the appendices, and only relate the conclusions as follows.
-
(i) First, the conditional LLEs of the slave system still characterise the synchronisation process in three-body experiments. Specifically, the decay rate of the synchronisation error is still approximately equal to the conditional LLE.
-
(ii) Second, the threshold wavenumbers for the three-body and the twin experiments are essentially the same. Therefore, the synchronisation process is insensitive to the nature of the master flow.
The above conclusions show that we may bypass the results on synchronisation error and focus on the conditional LLEs. Also, we do not need to distinguish the three-body from the twin experiments as far as the synchronisation threshold and the conditional LLEs are concerned. We will follow these principles in what follows.
4.1. Relationship between the unconditional LLEs for LES and filtered DNS
This subsection is devoted to the unconditional LLE for LES, i.e.
$\lambda (k_m=0)$
, and its relationship with the LLE for filtered DNS,
$\lambda _\Delta$
. To the best of our knowledge, the LLE for filtered DNS has not been reported before. Only unconditional LLE is considered in this subsection so we will simply use
$\lambda$
to represent
$\lambda (k_m = 0)$
.
Figure 3. Normalised unconditional LLEs
$\lambda \tau _k$
as functions of
$\eta /\Delta$
for the cases in (a) Group I and (b) Group II. The horizontal dash-dotted line indicates
$\lambda \tau _k = 0.12$
corresponding to the value for DNS with moderate Reynolds number (Boffetta & Musacchio Reference Boffetta and Musacchio2017; Berera & Ho Reference Berera and Ho2018; Inubushi et al. Reference Inubushi, Saiki, Kobayashi and Goto2023).
The main result is shown in figure 3, which compares the LLEs of LES with those of filtered DNS, with the cases in Group I plotted in panel (a) and those in Group II in panel (b). The normalisation parameters are given in tables 1–3. To compare with the results shown in figure 1 of Budanur & Kantz (Reference Budanur and Kantz2022) and figure 7 of Nastac et al. (Reference Nastac, Labahn, Magri and Ihme2017), we use
$\tau _k$
to normalise
$\lambda$
and plot
$\lambda \tau _k$
against
$\eta /\Delta$
. There are no DNS data with which to compare the high-Reynolds-number LES results. We plot the same set of low-Reynolds-number DNS results in panel (b), merely to highlight the scaling relation.
The most salient feature in the figure is the power-law behaviour emerging for smaller
$\eta /\Delta$
. For the low-Reynolds-number data, the scaling law can be observed for
$\eta /\Delta$
being less than a value between
$0.1$
and
$0.2$
for the three SGS stress models. For filtered DNS, it is observed for
$\eta /\Delta \lesssim 0.08$
. For the high-Reynolds-number data, the scaling behaviour appears to be present for all
$\eta /\Delta$
values covered in our simulations.
No scaling law can be discerned for larger
$\eta /\Delta$
values for either LES or filtered DNS results. It appears that this range of values for
$\eta /\Delta$
coincide with those of Nastac et al. (Reference Nastac, Labahn, Magri and Ihme2017) and Budanur & Kantz (Reference Budanur and Kantz2022), which might explain why no scaling regime was shown in their results. In this range, the LLE for the LES generally increases with
$\eta /\Delta$
. This trend is consistent with Nastac et al. (Reference Nastac, Labahn, Magri and Ihme2017) and Budanur & Kantz (Reference Budanur and Kantz2022). In addition, our results also show that the LLE is different for different models, with SSM and DMM producing the smallest and largest values, respectively. In contrast to the LLEs for LES, the LLE of the filtered DNS actually remains roughly a constant, a trend that is captured better by the DMM.
Returning back to the power-law region, the data suggest that the three models display the same scaling exponent and the scaling exponent is close to that of the filtered DNS. This observation can be quantified by fitting a power law in the form of
\begin{equation} \lambda \tau _k = a (\eta /\Delta )^b \end{equation}
to the data. The parameters
$a$
and
$b$
are found by least-squares fitting applied to the transformed data
$(\ln \eta /\Delta , \ln \lambda \tau _k)$
. The values for
$a$
and
$b$
and 4.1 are also shown in figure 3. Equation (4.1) appears as a straight line with slope
$b$
in the figure. The data in panel (a) show that the slopes for the LES results are slightly steeper than that for filtered DNS. However, the slopes for the LES results at higher Reynolds numbers, shown in panel (b), are slightly smaller, though we should note this is compared to data for filtered DNS obtained at lower Reynolds numbers. In any case, it is justified to say that the differences in the slopes between different cases are small.
Figure 4. Normalised unconditional LLE
$\lambda \tau _\Delta$
as functions of
$\eta /\Delta$
for cases in subgroups R4 and R5. Symbols, simulations results; lines, least-squares fitting based on (4.4).
There are conspicuous differences for coefficient
$a$
between SGS models and the filtered DNS. However, given the close agreement for the scaling exponent
$b$
, we argue that the LLE for an LES should be interpreted as an approximation to the LLE for the filtered DNS that the LES is simulating. This relationship, as intuitive as it may seem, has not been reported before. In this view, the LLE for an LES is a measure of the growth rate of the infinitesimal perturbation confined in the scales larger than
$\Delta$
in the DNS.
The above claim can be supported by the following phenomenological analysis. In fact, Kolmogorov phenomenology suggests that
\begin{equation} \lambda \sim \tau _\Delta ^{-1}, \qquad \lambda _\Delta \sim \tau _\Delta ^{-1}, \end{equation}
so that
$\lambda \tau _\Delta \sim \text {constant}$
for
$\Delta$
in the inertial range. In light of (3.6), this would imply
$\lambda \tau _k \sim \Delta ^{-2/3} \sim (\eta /\Delta )^{2/3}$
. Therefore, according to the Kolmogorov phenomenology,
$\lambda \tau _k$
should decrease as
$\Delta$
increases and the slope for
$\lambda \tau _k$
should be
$2/3$
. The empirical values of the exponent
$b$
again are close to this semi-analytical value. This phenomenological argument applies equally to
$\lambda$
and
$\lambda _\Delta$
, which shows that there is a physical basis for the two to be the same.
There is non-trivial discrepancy between the empirical values for
$b$
and
$2/3$
(observed for both filtered DNS and LES results). We explore this issue further with the high-Reynolds-number simulations in Group II, and demonstrate below that the deviation can be plausibly explained with current theories. We start by plotting the results for
$\lambda \tau _\Delta$
, which are shown with the symbols in figure 4. To quantify the deviation, we note (3.6) and (4.1) imply
\begin{equation} \lambda \tau _\Delta = a \left (\frac {\eta }{\Delta }\right )^{\alpha }, \quad \alpha = b - \frac {2}{3}. \end{equation}
Thus,
$\alpha$
can be found from the values of
$b$
we obtained previously (cf. figure 3
b). Here,
$\alpha$
obviously is smaller than
$b$
, indicating weaker dependence on
$\eta /\Delta$
. This is consistent with figure 4. Figure 4 also shows a systematic difference between results for subgroups R4 and R5 which have different Reynolds numbers (which thus have been plotted separately). This observation suggests that
$a$
in (4.3) depends on the Reynolds number, and we may write empirically
\begin{equation} \lambda \tau _\Delta = c_1 \left (\frac {\eta }{\Delta }\right )^{\alpha } Re_\lambda ^\beta \end{equation}
for some constants
$c_1$
and
$\beta$
. The values for
$c_1$
and
$\beta$
are found with linear regression, and are given in table 5. Equation (4.4) with coefficients given in table 5 is shown with lines in figure 4.
To compare the above scaling relation with results in the literature, we rewrite (4.4) in terms of the Reynolds number defined with an integral length scale
$L$
. Here,
$L$
can be defined via the energy spectrum as shown by, e.g., Batchelor (Reference Batchelor1953), though its precise formula is not important for the scaling argument below. Introducing the dissipation constant
$C_{\epsilon } \equiv L {\epsilon }_t/ u_{\textit{rms}}^3$
, it is easy to show that (4.4) implies
\begin{equation} \lambda \tau _\Delta \sim C_{\epsilon }^{-\beta /2-\alpha /4} Re_L^{\beta /2 - 3\alpha /4} \left (\frac {\Delta }{L}\right )^{-\alpha }, \end{equation}
where
$Re_L \equiv u_{\textit{rms}} L/\nu$
is the integral scale Reynolds number. As is reviewed by Vassilicos (Reference Vassilicos2015), the dissipation constant
$C_{\epsilon }$
decreases with increasing
$Re_L$
though perhaps only weakly at high Reynolds numbers. It also displays interesting scaling behaviours reported in non-equilibrium turbulence (Valente & Vassilicos Reference Valente and Vassilicos2012). With the values of
$\alpha$
and
$\beta$
obtained above, we then find
\begin{equation} \lambda \tau _\Delta \sim \begin{cases} \,\, C_{\epsilon }^{-0.12} Re_L^{0.06}\displaystyle \left (\frac {\Delta }{L}\right )^{-0.05} & \text {for SSM,}\\ \,\, C_{\epsilon }^{-0.06} Re_L^{0.01}\displaystyle \left (\frac {\Delta }{L}\right )^{-0.05} & \text {for DSM,} \\ \,\, C_{\epsilon }^{-0.10} Re_L^{0.02}\displaystyle \left (\frac {\Delta }{L}\right )^{-0.08} & \text {for DMM.} \end{cases} \end{equation}
Note that
$C_{\epsilon }$
tends to decrease as
$Re_L$
increases. Therefore,
$\lambda \tau _\Delta$
increases with the Reynolds number and increases when
$\Delta /L$
is decreased. Both are corrections to the Kolmogorov scaling. As one may see, they are both relatively small, not unexpected.
The corrections are reminiscent of what is observed for the Lyapunov exponents of unfiltered DNS (Boffetta & Musacchio Reference Boffetta and Musacchio2017; Berera & Ho Reference Berera and Ho2018), where small deviation from Kolmogorov scaling is also found. The quantity corresponding to
$\lambda \tau _\Delta$
in those investigations is
$\lambda _{DNS}\tau _k$
, where
$\lambda _{DNS}$
denotes the LLE for the unfiltered DNS. It is found that
$\lambda _{DNS} \tau _k \sim Re_L^{0.14}$
(Boffetta & Musacchio Reference Boffetta and Musacchio2017; Berera & Ho Reference Berera and Ho2018; Inubushi et al. Reference Inubushi, Saiki, Kobayashi and Goto2023), whereas Kolmogorov phenomenology suggests
$\lambda _{DNS} \tau _k \sim$
constant.
The origin of the deviation is debated. One might attribute it to small-scale intermittency (for reviews of small-scale intermittency in turbulence, see, e.g., Frisch (Reference Frisch1995) and Sreenivasan & Antonia (Reference Sreenivasan and Antonia1997)). One might reasonably expect the intermittency correction to increase with
$Re_L$
, which thus might account for the
$Re_L$
factor in (4.6). Meanwhile, intermittency is stronger for smaller filter scale
$\Delta$
. This is also consistent with (4.6). However, it should be noted that the multi-fractal formalism predicts opposite dependence on
$Re_L$
(Boffetta & Musacchio Reference Boffetta and Musacchio2017).
Another possible origin of the deviation is the sweeping effect, as argued recently by Ge et al. (Reference Ge, Rolland and Vassilicos2023). The sweeping effect introduces the sweeping time scale
$\eta /u_{\textit{rms}}$
(also known as the Eulerian time scale) into the parametrisation of
$\lambda _{DNS}\tau _k$
. The same argument, applied to
$\lambda \tau _\Delta$
, suggests that we should introduce a sweeping time scale
$\tau ^E_\Delta = \Delta /u_{\textit{rms}}$
into the parametrisation of
$\lambda \tau _\Delta$
, which thus implies
$\lambda \tau _\Delta$
is a function of
$\tau ^E_\Delta /\tau _\Delta$
. As a demonstration, let us assume
$\lambda \tau _\Delta \sim (\tau ^E_\Delta / \tau _\Delta )^c$
for some constant
$c$
. We then obtain
\begin{equation} \lambda \tau _\Delta \sim C_{\epsilon }^{c/3} \left (\frac {\Delta }{L}\right )^{c/3}. \end{equation}
As
$C_{\epsilon }$
decreases with
$Re_L$
(Goto & Vassilicos Reference Goto and Vassilicos2009; Vassilicos Reference Vassilicos2015), a negative
$c$
can potentially provide an explanation for the behaviour of
$\lambda \tau _\Delta$
expressed in (4.6).
Finally, the dependence of
$\lambda \tau _\Delta$
on
$Re_L$
can also be a finite-Reynolds-number effect. If the ratio
$\tau _k/\tau _\Delta$
is not negligible at our Reynolds numbers,
$\lambda \tau _\Delta$
should be a function of
$\tau _k/\tau _\Delta$
, and hence depends on
$Re_L$
and
$\Delta /L$
. To sum up, there are competing theories for the deviation from Kolmogorov phenomenology observed in
$\lambda \tau _\Delta$
. However, the fact that the deviation can be explained plausibly by these theories lends significant support to the contention that the scaling behaviour observed in
$\lambda$
is not incidental, but rather originates from the fact that
$\lambda$
is an approximation to
$\lambda _\Delta$
.
The knowledge about the relationship between
$\lambda$
and
$\lambda _\Delta$
can be beneficial. In this regard, one is tempted to compare
$\lambda$
(or
$\lambda _\Delta$
) with the finite size Lyapunov exponents (FSLEs) (Boffetta et al. Reference Boffetta, Cencini, Falcioni and Vulpiani2002; Boffetta & Musacchio Reference Boffetta and Musacchio2017). The FSLE measures the growth rate of a finite size velocity perturbation
$\delta v$
by calculating the doubling time of the perturbation. It is found that the FSLE is
$\sim (\delta v)^{-2}$
in isotropic turbulence for
$\delta v$
corresponding to fluctuations over length scales in the inertial range. There is significant difference between the LLEs of an LES velocity field and the FSLE, as the LLEs measures the growth of infinitesimal perturbation. Nevertheless, if we assume that the most unstable velocity perturbation scales with
$\Delta$
according to Kolmogorov scaling, we have
$\Vert \overline {\boldsymbol {u}}^\delta \Vert \sim \Delta ^{1/3}$
. Ignoring the dependence on
$Re_L$
and
$\Delta /L$
shown in (4.6), we have
$\lambda \sim \tau _\Delta ^{-1}$
, and thus,
\begin{equation} \lambda \sim \tau _\Delta ^{-1} \sim \frac {\Vert \overline {\boldsymbol {u}}^\delta \Vert }{\Delta } \sim \Vert \overline {\boldsymbol {u}}^\delta \Vert ^{-2}, \end{equation}
which mimics the scaling law for the FSLEs. This argument hints at the possibility of using LES to gain insight into the FSLEs.
4.2. Conditional LLEs and the threshold coupling wavenumber
We now move on to the results for the conditional LLEs
$\lambda (k_m)$
. Two flows coupled with coupling wavenumber
$k_m$
can be synchronised only when
$\lambda (k_m) \lt 0$
. One main parameter of interest is the threshold value
$k_c$
, i.e. the
$k_m$
value for which
$\lambda (k_m) = 0$
. We will present the results for
$\lambda (k_m)$
first, and then discuss those for
$k_c$
.
Figure 5. Normalised conditional LLEs
$\lambda \tau _\Delta$
as functions of
$k_m \Delta _t$
: (a) R1; (b) R3; (c) R4; (d) R5.
The data for
$\lambda (k_m) \tau _\Delta$
are plotted against
$k_m \Delta _t$
in figure 5, where
$\Delta _t$
is the effective dissipation length defined in (3.9). The reason why
$\Delta _t$
is used to normalise
$k_m$
will become clear when we discuss the results for
$k_c$
. Only the data from subgroups R1, R3, R4 and R5 are shown in figure 5. The results for R2 have been omitted as they fall somewhere between those of R1 and R3.
In all cases,
$\lambda (k_m) \tau _\Delta$
is a decreasing function of
$k_m \Delta _t$
. For the low-Reynolds-number cases in R1 (panel a) and R3 (panel b), the results for SSM and DSM only change slightly with
$Re_\lambda$
or the filter scale
$\Delta$
. In contrast, those for DMM depend on the two parameters quite significantly. For DMM,
$\lambda \tau _\Delta$
is larger for larger
$\Delta$
and this effect is the strongest for cases in R3 which have larger Reynolds numbers. Data for high-Reynolds-number cases are shown in figure 5(c,d). The results depend on
$Re_\lambda$
and
$\Delta$
to some extent, but the variation is moderate and it is reduced as
$k_m\Delta _t$
increases. The curves collapse well as
$\lambda (k_m)\tau _\Delta$
approaches zero.
Figure 6. Threshold wavenumber
$k_c \Delta _t$
for all subgroups R1–R5.
Figure 7. (a) The threshold wavenumber
$k_c \eta$
for all subgroups R1–R5. (b)
$\Delta _t/ \Delta$
. The slope of the dashed line is
$1$
.
The dimensionless threshold wavenumber
$k_c \Delta _t$
is the value of
$k_m\Delta _t$
where the curve for
$\lambda (k_m)$
crosses the horizontal axis. The values of
$k_c\Delta _t$
are read from figure 5 and plotted as a function of
$\eta /\Delta$
in figure 6. The interesting observation is that
$k_c \Delta _t$
for SSM and DSM are almost constant across the whole range of values for
$\eta /\Delta$
. For SSM,
$k_c\Delta _t \approx 0.18$
, and
$k_c \Delta _t\approx 0.20$
for DSM. These values, incidentally, are close to the DNS value of
$k_c \eta \approx 0.2$
. For DMM,
$k_c\Delta _t$
remains approximately a constant
$0.36$
for
$\eta /\Delta$
up to approximately
$0.08$
, and decreases towards
$0.2$
as
$\eta /\Delta$
increases further. Therefore, for large
$\eta /\Delta$
,
$k_c \Delta _t$
approaches the DNS value
$0.2$
for all models, which is to be expected since LES approaches DNS when
$\Delta \to \eta$
. At a given
$\eta /\Delta$
,
$k_c \Delta _t$
is the smallest for SSM, somewhat larger for DSM and the largest for DMM.
We may contrast the results for
$k_c \Delta _t$
with those for
$k_c \eta$
. The values of
$k_c\eta$
are plotted against
$\eta /\Delta$
in figure 7(a). The main observation is
$k_c\eta$
increases monotonically with
$\eta /\Delta$
for SSM and DSM, while for DMM, it increases monotonically up to a value close to
$0.2$
at
$\eta /\Delta \approx 0.13$
and then undulates for larger
$\eta /\Delta$
.
The normalisation
$k_c \eta$
has the advantage of being more intuitive, as the behaviour of
$\eta$
is well known. However,
$k_c\eta$
approaches zero for all models as
$\eta /\Delta \to 0$
, which obscures the differences between the models at small
$\eta /\Delta$
. Meanwhile,
$k_c \Delta _t$
is able to capture the differences and it displays better asymptotic behaviours as
$\eta /\Delta \to 0$
. Based on figure 6, we can extrapolate the results to much smaller
$\eta /\Delta$
, but it is difficult to do so based on figure 7.
To gain more insight into the above results, it is of interest to look into the behaviours of
$\Delta _t$
. Here,
$\Delta _t/\Delta$
is shown in figure 7(b) as a function of
$\eta /\Delta$
. We can see
$\Delta _t/\Delta$
approaches constants as
$\eta /\Delta \to 0$
. When
$\eta /\Delta \to 1$
, the curves tend towards
$\Delta _t/\Delta = \eta /\Delta$
, as can be seen from the dashed line. Therefore, roughly speaking,
$\triangle _t$
is an interpolation between
$\eta$
and
$\Delta$
. The asymptotic values for
$\Delta _t/\Delta$
as
$\eta /\Delta \to 0$
can be read from the figure. Using these values, we can find that, as
$\eta /\Delta \to 0$
,
$k_c\Delta$
approaches approximately
$0.16$
,
$0.13$
and
$0.22$
for SSM, DSM and DMM, respectively.
The data for
$k_c \Delta _t$
obtained in this subsection will be used to examine their relationship with the peak wavenumbers of the LLVs in the next subsection, which is our main objective. These data by themselves mainly serve to reveal the difference in the synchronisability of DNS and LES. The data could also supplement previous discussions related to the predictability of LES and the dynamic contents in LES (Nastac et al. Reference Nastac, Labahn, Magri and Ihme2017; Budanur & Kantz Reference Budanur and Kantz2022), which are important questions from the perspective of turbulent modelling.
4.3. Threshold wavenumbers and the leading Lyapunov vectors
As is reviewed in § 1, evidence suggests that the synchronisation threshold wavenumber obtained in the previous subsection should be the same as the peak wavenumber for the energy spectrum of the unconditional LLV for the slave system. Though the amount of evidence is significant, in this subsection, we further examine this relationship using LES synchronisation experiments (and later with the SABRA model). As we explain in § 1, here, LES is treated as a chaotic dynamical system that is different from DNS. Our thesis is that the relationship can be substantiated by showing it can also be observed in systems different from DNS. For brevity, we use LLVs to refer to the unconditional LLVs.
Figure 8. Average energy spectra of the unconditional LLV from the DNS data in subgroups R1, R2 and R3. The spectra are normalised so that the total energy is unity. The vertical line marks the value
$k\eta = 0.2$
.
We denote the average energy spectrum of the LLV by
$E_\Delta (k)$
. To give a first impression of the relationship, we plot in figure 8
$E_\Delta (k)$
from our DNS data. This figure is essentially the same as figure 2(d) of Budanur & Kantz (Reference Budanur and Kantz2022), although here we use a linear scale for
$E_\Delta (k)$
to highlight the profiles around the peaks of the spectra. The figure shows that the spectra do peak approximately at
$k\eta \approx 0.2$
, which coincides with the threshold wavenumber found in synchronisation experiments for isotropic turbulence.
Figure 9. Energy spectrum
$E_\Delta (k)$
for LES with (a,b) SSM, (c,d) DSM and (e,f) DMM. Results for lower Reynolds numbers are shown in panels (a,c,e) whereas results for higher Reynolds numbers are shown in panels (b,d,f).
The results for
$E_\Delta (k)$
calculated from LES data are given in figure 9. Following the normalisation we adopted previously,
$E_\Delta (k)$
is plotted against
$k \Delta _t$
, where
$\Delta _t$
can be different for different cases. We may comment that, in most cases,
$E_\Delta (k)$
peaks at an intermediate wavenumber. That is, the perturbations with energy localised on intermediate wavenumbers are the most unstable. However, what we are most interested in is the wavenumber
$k_p$
where
$E_\Delta (k)$
finds its maximum. To compare
$k_p$
with
$k_c$
quantitatively, we need to extract the precise values of
$k_p$
from the spectra. There are challenges for this since the gap between two data points can be fairly large around the peaks. To reduce the uncertainty, we use the method of Li et al. (Reference Li, Tian, Li, Si and Mohammed2024). Namely, we fit a smooth curve to the spectrum using cubic splines. The peak wavenumber of the fitted curve is taken to be
$k_p$
. More details of the method can be found from Li et al. (Reference Li, Tian, Li, Si and Mohammed2024).
We consider the results for lower Reynolds numbers first. The normalised peak wavenumbers
$k_p\Delta _t$
for the LES in subgroups R1–R3 are plotted in figure 10(a) together with the data for
$k_c\Delta _t$
in these cases extracted from figure 6. These data can be cross-checked with the distributions given in figure 9. There are some discrepancies between the two quantities, but the overall agreement is evident. The discrepancy is somewhat larger for DMM at low
$\eta /\Delta$
, where
$k_p$
is smaller than
$k_c$
. However, it is significant to observe that both
$k_p$
and
$k_c$
decrease when
$\eta /\Delta$
increases, despite the fact that the two quantities are obtained in entirely different ways.
Figure 10. Comparison between the threshold coupling wavenumber
$k_c$
and the peak wavenumber
$k_p$
of the LLV: (a) data from cases in R1, R2 and R3 (lower Reynolds numbers); (b) data from cases in R4 and R5 (high Reynolds numbers), and for SSM and DMM only.
For the high-Reynolds-number cases, the physical picture is more complicated. The peak wavenumber
$k_p$
in this case is shown in figure 10(b), which corresponds to figure 9(b,d,f). As we can see from the latter, for SSM and DMM,
$E_\Delta (k)$
displays a clear peak at a wavenumber that is roughly constant for different Reynolds numbers and different
$N_{LES}$
. The peak wavenumber
$k_p$
is plotted in figure 10(b). We can see that the agreement between
$k_p$
and
$k_c$
is excellent for these models over a wide range of values of
$\eta /\Delta$
. However, the results for DSM show unexpected trends. On coarse grids, such as in the cases with
$N_{LES} = 32$
and
$64$
,
$E_\Delta (k)$
for DSM increases with
$k \Delta _t$
monotonically. Therefore, it does not peak at an intermediate wavenumber, which means that the relationship between
$k_p$
and
$k_c$
breaks down in these cases. However, it is interesting to observe that the middle part of the spectrum appears to be amplified relative to the rest of the spectrum when
$N_{LES}$
is increased. This trend can be seen from the diamonds and circles in figure 9(d). As a matter of fact, the spectrum for the case R4DSM256 develops a local maximum at
$k \Delta _t \approx 0.2$
. This value again is close to the value for
$k_c \Delta _t$
found for the same case. This local maximum, however, is not the global maximum. The latter is found at the tail of the spectrum. Nevertheless, this observation suggests that the same mechanism creating the peak in
$E_\Delta (k)$
for other models is also at play here, and a peak may emerge at higher
$N_{LES}$
for the spectrum for DSM too. It would be interesting to verify this conjecture in the future with larger simulations.
To summarise, we find strong, though not perfect, evidence to show that the relationship between the threshold wavenumber for synchronisation and the peak wavenumber of the LLV also holds in LES. This finding adds to the evidence already found in DNS of rotating and non-rotating turbulence, thus suggesting that the relationship is universal among a class of turbulent flows. We further strengthen this claim with experiments conducted with the SABRA model below.
4.4. Synchronisation threshold for the SABRA model
This section focuses on the results of the twin experiments for the SABRA model. The computational aspect of the experiments has been described before. Some additional measures that are taken to ensure the accuracy of the results are described here.
When integrating a twin system of the SABRA model, we calculate simultaneously the cumulative average of the energy spectrum of the LLV
$E_\Delta (k_n)$
. The convergence of
$E_\Delta (k_n)$
and the conditional LLE
$\lambda (k_m)$
is monitored and checked at fixed time intervals. Let
$t_e$
represent the relative change in the core part of
$E_\Delta (k_n)$
, which, more precisely, is defined as
\begin{equation} t_e \equiv \max _n \left \{\left .\frac {\vert \delta E_\Delta (k_n)\vert }{ E_\Delta (k_n)} \,\right \vert \,E_\Delta (k_n)\ge 10^{-3} \right \}, \end{equation}
where
$\delta E_\Delta (k_n)$
represents the change in
$E_\Delta (k_n)$
between current and previous checks. Similarly,
$t_\ell \equiv \vert \delta \lambda (k_m)\vert /\vert \lambda (k_m)\vert$
represents the relative error in
$\lambda (k_m)$
between consecutive checks. The time integration of the twin systems is conducted for at least
$9\times 10^7$
time steps, and is terminated afterwards if both
$t_e$
and
$t_\ell$
are found to be less than
$5 \times 10^{-4}$
for five consecutive checks.
The final result for
$E_\Delta (k_n)$
is averaged over forty realisations simulated in the manner described above. Only one realisation is used to calculate
$\lambda (k_m)$
. As we will see later,
$E_\Delta (k_n)$
displays features unseen in LES (cf. figure 12). Our additional measures make sure that the differences are not numerical artefacts.
We now begin with some basic results about the energy spectrum
$E(k_n)$
and the energy flux
$\Pi (k_n )$
. Here,
$E(k_n)$
is defined by
\begin{equation} E(k_n) = \frac {1}{2} k_n^{-1} \langle \vert u_n\vert ^2 \rangle . \end{equation}
For the energy flux, we use a definition suggested by Biferale (Reference Biferale2003):
\begin{equation} \Pi (k_n ) = \Im \left \{ \langle u_n^* u_{n+1}^* u_{n+2} \rangle + \frac {1-b}{q} \langle u_{n-1}^* u_n^* u_{n+1}\rangle \right \}, \end{equation}
where
$\Im$
indicates the imaginary part,
$b$
is defined in (2.21) and
$q$
is the ratio between two consecutive wavenumbers. For more discussions about these two quantities, including alternative definitions for
$\Pi (k_n)$
, see e.g. Ditlevsen (Reference Ditlevsen2011) and Biferale (Reference Biferale2003).
The spectrum, normalised by the parameters given in table 4, is shown in figure 11(a), for all ten Reynolds numbers. The data collapse on one single curve that exhibits a
$-5/3$
slope over a very wide range of scales. The normalised energy flux
$\Pi (k_n )$
is shown in figure 11(b). It plateaus at a value of
$1$
for
$k_n \eta \lesssim 10^{-2}$
. This behaviour mimics that of the energy flux in energy cascade for high-Reynolds-number turbulence, which is expected to be a constant over the scales in the inertial ranges.
Figure 11. (a) The (normalised) energy spectrum
$E(k_n)$
. (b) The (normalised) energy flux
$\Pi (k_n)$
. (c) Synchronisation error
${\mathcal E}_{AB }(t)$
for case
$Re_0$
, with coupling wavenumber
$k_m = k_7$
,
$k_8$
,
$k_9$
and
$k_{10}$
.
Figure 12. (a) Normalised conditional LLE
$\lambda (k_m)$
for the SABRA model. (b) Spectra
$E_\Delta (k_n)$
for the LLV.
Figure 13. Threshold coupling wavenumber
$k_c$
obtained from the conditional LLEs
$\lambda (k_m )$
versus the peak wavenumber
$k_p$
of the energy spectrum of the LLV
$E_\Delta (k_n)$
.
Figure 11(c) plots the synchronisation error
${\mathcal E}_{AB}(t ) \equiv \Vert \boldsymbol {z }_A - \boldsymbol {z}_B\Vert$
for the case
$Re_0$
for several coupling wavenumber
$k_m$
, which illustrates the general features of the synchronisation process. The figure demonstrates that the SABRA model can be synchronised in the usual twin experiments. Here,
${\mathcal E}_{AB}(t)$
also decays exponentially for sufficiently large
$k_m$
, but it displays stronger fluctuations than what is observed in LES (cf. figure 14 in Appendix A). Figure 11 shows that synchronisation fails for small
$k_m$
, e.g. for
$k_m = k_7$
, and that the threshold wavenumber
$k_c$
is somewhere between
$k_8$
and
$k_9$
in this case.
Figures 12 and 13 are the main results we intend to obtain from the synchronisation experiments for the SABRA model. Figures 12(a) and 12(b) show the conditional LLEs
$\lambda (k_m)$
and
$E_\Delta (k_n)$
, respectively. The threshold wavenumbers
$k_c$
obtained from these two plots are compared in figure 13. Before focusing on
$k_c$
, we comment briefly on the features of
$\lambda (k_m)$
and
$E_\triangle (k_n )$
. The distribution of
$\lambda (k_m)$
, as can be seen in figure 12(a), follows the general trend of the conditional LLEs for LES, with some small differences. Specifically,
$\lambda (k_m)$
is not necessarily a decreasing function of
$k_m$
; for small
$k_m \eta$
, it may increase with
$k_m \eta$
. However,
$\lambda (k_m)$
eventually becomes negative for sufficiently large
$k_m\eta$
, and its distributions seem to collapse onto one curve as
$k_m\eta$
becomes bigger.
The spectrum
$E_\Delta (k_n)$
, however, displays more significant differences. Most importantly, the distributions might have multiple peaks. For the Reynolds numbers we consider, there may be one to three peaks. The existence of multiple peaks might be related to the fact that each shell
$u_n$
in the SABRA model is only coupled to the neighbouring two shells on each side through the nonlinear term, but we will not explore this further.
The fact that
$E_\triangle (k_n)$
may have multiple peaks indicates that there might be multiple mechanisms behind the amplification of the most unstable perturbation. This observation raises the question of whether the maximum of
$E_\Delta (k_n)$
plays the same role as that observed in LES. Nevertheless, we can still follow our practice with LES results and identify the wavenumber
$k_p$
where the maximum of
$E_\Delta (k_n)$
is located. Note that
$k_p$
has to be one of the discrete wavenumbers
$k_n$
, because we do not apply smoothing here when we identify the peak of
$E_\Delta (k_n)$
. On the other hand, the threshold wavenumber
$k_c$
can be identified from
$\lambda (k_m)$
, as it is the root of the equation
$\lambda (k_m) = 0$
. We find
$k_c$
by interpolating between the two wavenumbers bracketing the root of
$\lambda (k_m)$
(cf. figure 12
a), so it does not have to be one of the discrete wavenumber
$k_n$
. The two nearest neighbouring wavenumbers are also recorded as indicators for the uncertainty in determining the threshold wavenumber.
To test our hypothesis that
$k_p$
is the same as the threshold wavenumber
$k_c$
, the two wavenumbers are compared in figure 13. The two wavenumbers neighbouring
$k_c$
are represented by the error bars in the figure. We observe that
$k_c \eta$
is almost a constant over a wide range of Reynolds numbers and its average value is
$0.035$
. Meanwhile,
$k_p$
oscillates somewhat around this value, but it never exceeds the error bars. Actually, it always lands on one of the two error bars. That is, the difference between
$k_p$
and
$k_c$
falls within the uncertainty introduced by the discretisation of the wavenumber
$k_n$
. This observation again supports our hypothesis that the peak of the energy spectrum of the LLV is the same as the synchronisation threshold wavenumber, which corroborates the observations we made previously based on LES with SSM and DMM as the SGS stress models.
5. Conclusions
The relationship between the threshold coupling wavenumber for two turbulent flows and the leading Lyapunov vector (LLV) for the flows is investigated via a series of large eddy simulations (LESs) covering different subgrid-scale (SGS) stress models, filter scales, Reynolds numbers, coupling methods, and different combinations between the master and the slave flows. The large eddy simulations are complemented by a series of synchronisation experiments based on the SABRA model. Our main conclusions can be summarised as follows.
First, we report a new scaling relation for the leading Lyapunov exponents (LLEs) of LES with canonical SGS stress models. We introduce the concept of the LLE of filtered DNS, and show that it displays same scaling behaviour observed in the LLEs of LES. Phenomenological arguments are presented to give plausible explanations to the empirical data, which thus show that the LLE of LES can be seen as an approximation to the LLE of the filtered DNS. It appears that this intuitive interpretation of the LLEs of LES has not been expounded before. The interpretation strengthens the physical foundation of the results related to the LLEs of LES reported in the literature.
Second, we show that synchronisation experiments conducted with LES and the SABRA model strongly support the hypothesis that, for a class of chaotic systems, the peak wavenumber of the energy spectrum of the LLV is the same as the threshold coupling wavenumber for synchronisation. Evidence shows that DNS of box turbulence and the SABRA model are among this class. LES for box turbulence with the standard Smagorinsky model or the dynamic mixed model also belongs to the class. However, LES with the dynamic Smagorinsky model generally does not display such a relationship.
Finally, we show that the conditional LLEs characterise the synchronisation process in indirectly coupled systems, as they do in directly coupled systems. The synchronisation process is also insensitive to the nature of the master flow. These additional results help to establish the domain of applicability of our results.
The relationship between the threshold coupling wavenumber and the peak wavenumber of the LLV extends the finding in rotating turbulence reported by Li et al. (Reference Li, Tian, Li, Si and Mohammed2024), where a qualitative argument was also put forward as an explanation for the relationship. The argument hypothesises that the peak of the energy spectrum of the LLV corresponds to the most unstable Fourier modes in the velocity perturbation, and synchronisation can be achieved when the growth of these modes is suppressed by coupling. Therefore, the coupling wavenumber should be approximately the same as the peak wavenumber of the energy spectrum of the LLV. We suspect that this argument is similarly applicable here. However, an analytical theory for the relationship is still lacking and will be the subject of our future research.
Funding.
J.L. acknowledges the support of the National Natural Science Foundation of China (No. 12102391).
Declaration of interests.
The authors report no conflict of interest.
Data availability statement.
The data that support the findings of this study are available from the corresponding author upon reasonable request.
Appendix A. Decay rates of synchronisation error and conditional Lyapunov exponents
The synchronisation between
$\overline {\boldsymbol {u}}_{A}$
and
$\overline {\boldsymbol {u}}_{B}$
is characterised by the synchronisation error
\begin{equation} {\mathcal E}_{AB}(t) = \Vert \overline {\boldsymbol {u}}_{A}-\overline {\boldsymbol {u}}_{B}\Vert . \end{equation}
Synchronisation happens when
${\mathcal E}_{AB}(t) \to 0$
as
$t\to \infty$
. Empirically, it has been found
${\mathcal E}_{AB}(t) \sim \exp (\alpha t)$
, i.e.
${\mathcal E}_{AB}$
decays exponentially when synchronisation happens, where
$\alpha$
is the decay rate which depends on
$k_m$
. In twin experiments, it has been shown (Nikolaidis & Ioannou Reference Nikolaidis and Ioannou2022; Li et al. Reference Li, Tian, Li, Si and Mohammed2024) that
$\alpha$
is equal to the conditional LLE
$\lambda (k_m)$
. The question we address here is whether they are also the same for three-body experiments.
Figure 14. Comparison between the decay rates of the synchronisation error
${\mathcal E}_{AB}(t)$
and the conditional LLEs of the slave flow.
Figure 14 shows the results for
${\mathcal E}_{AB}(t)$
with symbols for selected three-body experiments. It is clear from the figure that the decay of
${\mathcal E}_{AB}(t)$
on average is exponential for sufficiently large
$k_m$
. The synchronisation can happen with very different rates for different models and different
$k_m$
. The lines in the figure represent exponential decay given by
$\exp [\lambda (k_m) t ]$
, where
$\lambda (k_m)$
is the conditional LLE for the slave
$\overline {\boldsymbol {u}}_A$
. The agreement between the symbols and the lines shows that the decay rate of the synchronisation error is also the same as
$\lambda (k_m)$
in these three-body experiments. Therefore, the synchronisation in three-body experiments can also be fully characterised by the conditional LLEs.
Appendix B. The impacts of coupling on synchronisation threshold
We now present data on the threshold wavenumber
$k_c \Delta _t$
to show that it is insensitive to the coupling mechanism or the master flows. Figure 15 compares the thresholds obtained from the three-body experiments with those obtained from the twin experiments. The latter is the same as those shown in figure 6. For the data shown in panel (a) with lower Reynolds numbers, the master in the three-body experiments is a DNS. For the high-Reynolds-number data shown in panel (b), the master is an LES with an SGS stress model different from that in the slaves. The three cases shown in the figure have a
$256^3$
DMM as the master and two identical SSM as the slaves, a
$256^3$
DMM as the master and two identical DSM as the slaves, and a
$256^3$
DSM as the master and two identical DMM as the slaves, respectively. For the sake of brevity, we do not conduct further tests (e.g. LES with DSM as the master driving two identical LES with SSM).
Figure 15. Threshold wavenumber
$k_c \Delta _t$
. Dashed lines with empty symbols, three-body experiments; solid lines with solid symbols, twin experiments (same as those in figure 6). (a) Low-Reynolds-number cases (R1, R2 and R3). (b) High-Reynolds-number cases (R4 and R5).
Figure 15 shows that the threshold
$k_c \Delta _t$
for the twin experiments are very close to those obtained in the three-body experiments. This figure is the evidence for the claim made at the beginning of this subsection. Given the close agreement between the threshold wavenumbers, it is reasonable to infer that the conditional LLEs obtained from the three-body experiments would be close to those from the twin experiments too. We omit the numerical tests for the latter.
Open access
