Diffusive dynamics and stochastic models of turbulent axisymmetric wakes

Abstract

A modelling methodology to reproduce the experimental measurements of a turbulent flow in the presence of symmetry is presented. The flow is a three-dimensional wake generated by an axisymmetric body. We show that the dynamics of the turbulent wake flow can be assimilated by a nonlinear two-dimensional Langevin equation, the deterministic part of which accounts for the broken symmetries that occur in the laminar and transitional regimes at low Reynolds numbers and the stochastic part of which accounts for the turbulent fluctuations. Comparison between theoretical and experimental results allows the extraction of the model parameters.

1. Introduction

Turbulent flows are ubiquitous in natural phenomena and engineering applications (Pope 2000), and therefore a mathematically tractable description of them is desirable for their prediction and control. At low Reynolds numbers, corresponding to laminar regimes, bifurcation theory has aided understanding of the dynamic behaviour of fluid flows in the presence of symmetry: under certain assumptions, the infinite-dimensional nonlinear fluid system governed by the Navier–Stokes equations can be reduced to a simple and finite-dimensional system close to the threshold of bifurcation (Crawford & Knobloch 1991; Holmes et al. 2012). However, departure from the laminar regimes and the critical bifurcating points renders the flow chaotic and finally turbulent, increasing the order of the system and the complexity for a mathematical description of it.

Here we analyse the turbulent wake generated by an axisymmetric bluff body at high Reynolds numbers. This type of flow undergoes a steady bifurcation followed by an unsteady one at low Reynolds numbers (Bohorquez et al. 2011; Bury & Jardin 2012), the dynamics of which can be accurately captured using simple, weakly nonlinear models (Fabre, Auguste & Magnaudet 2008; Meliga, Chomaz & Sipp 2009) prior to the emergence of chaos and turbulent behaviour. During the transitional regime, continuous spatial and temporal symmetries are spontaneously broken. Although bifurcations break the symmetries en route to turbulence, fully developed turbulence is known to restore the possible symmetries in a statistical sense at very high Reynolds numbers (Frisch 1996).

In this paper we show that the dynamics of the turbulent wake can be approximated by a simple nonlinear Langevin equation, the deterministic part of which accounts for the broken symmetries that occur in the laminar and transitional regimes and the stochastic part of which accounts for the turbulent fluctuations.

The paper first presents an overview of axisymmetric wake flows and models in the laminar (§ 1.1) and turbulent (§ 1.2) regimes, before presenting the stochastic model in § 2. The model is validated using experimental data (§ 3), and the results are presented in § 4.

1.1. Laminar regime and modelling

The laminar and linearly stable axisymmetric wake loses spatial rotational symmetry, $\mathscr{R}_{{\rm\pi}}$ , in the azimuthal direction (rotation of angle ${\rm\pi}$ around any radial axis passing through the centre of the body) due to a supercritical pitchfork bifurcation (Fabre et al. 2008; Meliga et al. 2009; Bohorquez et al. 2011; Bobinski, Goujon-Durand & Wesfreid 2014), the normal form of which reads

(1.1) $$\begin{eqnarray}\dot{\boldsymbol{x}}={\it\alpha}\boldsymbol{x}+{\it\lambda}\boldsymbol{x}|\boldsymbol{x}|^{2}\end{eqnarray}$$

( ${\it\alpha},{\it\lambda}\in \mathbb{R},{\it\lambda}<0$ ). Above the threshold of instability ( ${\it\alpha}>0$ ), (1.1) is associated with symmetry breaking since the steady-state solutions are not invariant under the $\boldsymbol{x}\rightarrow -\boldsymbol{x}$ symmetry. This model can be directly obtained from the Navier–Stokes equations through a weakly nonlinear expansion around the critical bifurcating point (Meliga et al. 2009) and has been used extensively for the description of laminar flows undergoing supercritical bifurcations (Drazin & Reid 2004). In the context of the weakly nonlinear analysis, $\boldsymbol{x}=\boldsymbol{x}(t)$ describes the complex amplitude evolution of the bifurcated state (global mode). Then, the model coefficient ${\it\alpha}$ denotes the linear growth rate of the unstable bifurcated mode and ${\it\lambda}$ the coefficient of the nonlinear saturating term.

On increasing the Reynolds number, a subsequent Hopf bifurcation is responsible for the loss of continuous temporal symmetry and the emergence of periodic shedding of vortices before the turbulent axisymmetric wake reaches a chaotic regime (Mittal, Wilson & Najjar 2002; Bury & Jardin 2012).

1.2. Turbulent regime

At high Reynolds numbers, a large number of studies have shown that the bifurcating states observed in the laminar wakes persist and manifest as large-scale structures. For the axisymmetric wake, it was recently shown that the steady bifurcated state (responsible for the loss of rotational symmetry) and the unsteady bifurcated state (responsible for the vortex shedding) persist at high Reynolds numbers and fully turbulent regimes (Grandemange, Gohlke & Cadot 2014; Rigas et al. 2014). Similar behaviour has been observed in wakes with other types of symmetries (see, for example, the rectilinear three-dimensional Ahmed body (Grandemange, Cadot & Gohlke 2012; Grandemange, Gohlke & Cadot 2013)).

Here we extend the laminar, weakly nonlinear modelling approach to the fully turbulent regime, and we show that the derived model captures the dynamic evolution of the turbulent three-dimensional wake. This is achieved by modelling the effect of the turbulent fluctuations acting on the deterministic dynamics of the system, reminiscent of the transitional instabilities observed in the laminar regimes, as stochastic forcing. The stochastic modelling approach of the turbulent dynamics has been successfully applied to various turbulent flows exhibiting global and large-scale dynamics including the turbulent Rayleigh–Bénard convection (Sreenivasan, Bershadskii & Niemela 2002; Brown & Ahlers 2007) and the turbulent von Kármán swirling flow (de la Torre & Burguete 2007). Hereafter we demonstrate our modelling approach considering the dynamics associated with the first steady bifurcation given by (1.1).

2. The stochastic model

We consider the deterministic system given by (1.1) and apply a stochastic modelling approach for the turbulent background fluctuations. In Cartesian coordinates, $\boldsymbol{x}=(x,y)$ , the deterministic system with independent additive white noise becomes

(2.1) $$\begin{eqnarray}\displaystyle \left.\begin{array}{@{}lcl@{}}{\dot{x}}\hspace{3.0pt} & =\hspace{3.0pt} & {\it\alpha}x+{\it\lambda}x(x^{2}+y^{2})+{\it\sigma}{\it\xi}_{x}(t),\\ {\dot{y}}\hspace{3.0pt} & =\hspace{3.0pt} & {\it\alpha}y+{\it\lambda}y(x^{2}+y^{2})+{\it\sigma}{\it\xi}_{y}(t),\end{array}\right\} & & \displaystyle\end{eqnarray}$$

where $\boldsymbol{{\it\xi}}$ accounts for the random forcing of turbulence and ${\it\sigma}^{2}$ the variance of it. Here, ${\it\sigma}\equiv {\it\sigma}_{x}={\it\sigma}_{y}$ and the stochastic process in (2.1) is rotationally symmetric. The coefficients ${\it\alpha}$ and ${\it\lambda}$ have the same meaning as those in the laminar regime. We transform the system from $(x,y)\rightarrow (r,{\it\phi})$ , where $r$ is the radial distance from the centre (amplitude) and ${\it\phi}$ the angle (phase), using the Ito interpretation (Gardiner 1985). In polar variables, the Langevin system given by (2.1) becomes

(2.2a,b ) $$\begin{eqnarray}\displaystyle {\dot{r}}={\it\alpha}r+{\it\lambda}r^{3}+\frac{{\it\sigma}^{2}}{2r}+{\it\sigma}{\it\xi}_{r},\quad \dot{{\it\phi}}=\frac{{\it\sigma}}{r}{\it\xi}_{{\it\phi}}, & & \displaystyle\end{eqnarray}$$

where $r=\sqrt{x^{2}+y^{2}}$ and ${\it\phi}=\tan ^{-1}(y/x)$ . Notice that in polar variables the radial component is independent of the angular position.

The stationary probability density function (PDF) of the above system can be found from the steady-state Fokker–Planck equation and is given by

(2.3) $$\begin{eqnarray}P_{s}(r,{\it\phi})=C\exp \left[-\frac{U(r,{\it\phi})}{K}\right],\end{eqnarray}$$

where $C$ is a normalisation constant, $K={\it\sigma}^{2}/2$ is the noise intensity (diffusivity) and $U$ the potential. The potential $U$ is

(2.4) $$\begin{eqnarray}U(r,{\it\phi})=-\!\left[\frac{{\it\alpha}r^{2}}{2}+\frac{{\it\lambda}r^{4}}{4}+K\ln r\right].\end{eqnarray}$$

The stationary PDF for the angle is uniform, $\int _{0}^{\infty }P_{s}(r,{\it\phi})\,\text{d}r=1/(2{\rm\pi})$ . In the case of a rotationally symmetric experimental set-up and inflow conditions, the drifts and diffusivities are independent of ${\it\phi}$ . The stationary PDF, figure 1, peaks around the minimum of the Mexican-hat-shaped potential.

3. Experimental set-up

The turbulent wake is generated by an axisymmetric bluff body with a blunt trailing edge, a schematic of which is shown in figure 2. The body diameter $D$ is 196.5 mm and the length-to-diameter ratio $L/D$ is 6.48. More details of the experimental set-up and procedure can be found in Rigas et al. (2014) and Oxlade et al. (2015). Experiments were performed at a constant free stream velocity, $U_{\infty }=15~\text{m}~\text{s}^{-1}$ . The Reynolds numbers based on diameter and boundary layer momentum thickness at separation are $\mathit{Re}_{D}=1.88\times 10^{5}$ and $\mathit{Re}_{{\it\theta}}=2050$ respectively. The boundary layer separation and the ensuing wake are turbulent.

Figure 1. Stationary probability density function (a) and potential (b) of the nonlinear Langevin model.

Figure 2. Schematic representation of the model. (a) Side view; (b) rear view. Pressure measurements are taken on the base of the body from 64 static pressure tappings, based on which the centre-of-pressure location is calculated.

Pressure measurements are taken on the base of the body from 64 static taps equally spaced on a uniform polar grid with ${\it\delta}r=0.056D$ and ${\it\delta}{\it\phi}=45^{\circ }$ in the radial and azimuthal directions, respectively. The sampling frequency was 225 Hz. A total of 19 200 s of data was acquired over 16 independent experiments, providing approximately 2000 independent measurements with 95 % uncertainties of approximately 0.45 % and 1 % in time average and r.m.s. pressure respectively.

4. Results

Let $\boldsymbol{x}$ of (1.1) represent the centre-of-pressure (CoP) coordinates. The CoP is used as a macroscopic variable quantifying the symmetry of the turbulent wake. Here it is calculated from the space-averaged pressure, defined on the Cartesian coordinate system of the base, $\boldsymbol{x}^{\prime }=(x^{\prime },y^{\prime })$ , and non-dimensionalised with the body diameter $D$ as

(4.1) $$\begin{eqnarray}\boldsymbol{x}(t)=\frac{1}{D\displaystyle \int p(t)\,\text{d}A}\int _{A}p(t)\boldsymbol{x}^{\prime }\,\text{d}A,\end{eqnarray}$$

where $A$ is the area of the base of the body. A zero CoP value corresponds to an $\mathscr{R}_{{\rm\pi}}$ -symmetric flow, whereas departure from this value indicates an increased asymmetry of the flow (Rigas et al. 2014). In polar coordinates, the CoP components are ( $r,{\it\phi}$ ). In the subsequent analysis, $r$ is the radial distance of the CoP from the centre of the body normalised with the body diameter and ${\it\phi}$ represents the orientation of the wake expressed in radian units.

Figure 3. A time series and probability density of the centre of pressure: angular (a) and radial (b) position. Symbols in the probability density (right) correspond to experimental data and solid lines to theoretical model predictions. Dashed line: mean radial value indicating the broken $\mathscr{R}_{{\rm\pi}}$ symmetry.

The temporal evolution and PDF of the radial and azimuthal location of the CoP for the highly turbulent regime are shown in figure 3. Highly erratic motion with time is observed in both components. Statistically, the wake spends most of the time in a non-zero radial location, indicating broken $\mathscr{R}_{{\rm\pi}}$ symmetry. Due to the uniform PDF in the azimuthal location, rotational symmetry is recovered in the long time average. For the laminar regime and after the first bifurcation, which is described by (1.1), the stable equilibrium corresponds to a non-zero CoP, the angle of which is determined from the initial conditions and is unique. During the turbulent regime, the laminar stable equilibrium point explores a continuum of states around its mean radial value that is uniformly distributed in the angular direction: hence, the wake explores an infinite number of metastable states, restoring the lost $\mathscr{R}_{{\rm\pi}}$ symmetry.

Figure 4. Mean-square displacement as a function of the time interval ${\it\tau}$ . (a) Angular and (b) radial components of the CoP location. Symbols: experimental data. Solid line: model. Thick dashed lines: power-law fits. Insets show power spectral density of the experimental CoP components. An exponent of $-2$ , consistent with free diffusive motion, is obtained for the azimuthal and radial components. At low frequencies, the PSD of the radial components is constant due to the spatial confinement imposed by the potential well, in accordance with the MSD results.

Insight into the random dynamics of the turbulent wake is provided from the calculation of the time-averaged mean-square displacement (MSD), defined as $\langle [{\rm\Delta}\boldsymbol{x}({\it\tau})]^{2}\rangle =\langle (\boldsymbol{x}(t+{\it\tau})-\boldsymbol{x}(t))^{2}\rangle$ . The MSD of the angular and radial components is plotted in figure 4 for different sampling times. In the azimuthal direction, the MSD increases linearly with time, $\langle [{\rm\Delta}{\it\phi}({\it\tau})]^{2}\rangle \propto {\it\tau}$ , consistent with free diffusive motion (Einstein 1905). In the radial direction, the linear relation holds only for short time scales below a threshold $t_{s}$ , $\langle [{\rm\Delta}r({\it\tau})]^{2}\rangle \propto {\it\tau},{\it\tau}<t_{s}$ , and reaches a saturation plateau at larger time scales, $\lim _{{\it\tau}\rightarrow \infty }\langle [{\rm\Delta}r({\it\tau})]^{2}\rangle =H^{2}$ .

The above results for the CoP are consistent with the predictions of the model given by (2.2). The coefficients of the Langevin equation, ${\it\alpha}$ , ${\it\lambda}$ and $K$ , are obtained from the experimental MSD and PDF and are given in table 1. The slope of the radial MSD relation is directly correlated with the diffusion coefficient $K$ , which is obtained through linear fitting, $\langle [{\rm\Delta}r({\it\tau})]^{2}\rangle =2K{\it\tau},~{\it\tau}<t_{s}$ . Knowing the diffusion coefficient $K$ , the model coefficients ${\it\alpha}$ and ${\it\lambda}$ are uniquely defined from (2.3) and here are obtained through least-squares fitting.

In figure 3, we compare the measured probability density function with the one predicted by the model from (2.3). In figure 4 we compare the MSD of the experimental results with the numerically calculated MSD from the model. Direct numerical integration of (2.2) was performed using an Euler–Maruyama scheme. We notice that the dynamics of the CoP can be described over all the time scales from the Langevin model.

Table 1. Coefficient values and associated errors for the Langevin model obtained experimentally at $\mathit{Re}_{D}=1.88\times 10^{5}$ . The units of the coefficients are $\text{s}^{-1}$ since lengths are non-dimensionalised with the body diameter.

A qualitative explanation of the above results is provided by the potential well shown in figure 1. The turbulent wake, the state of which is quantified by the CoP, meanders in the Mexican-hat-shaped potential and explores an infinite number of states through a random walk (diffusive motion). Specifically, in the azimuthal direction, which determines the orientation of the wake, it can explore freely any azimuthal location, resulting in unbounded reorientations. In the radial direction, the motion is restricted due to spatial constraints (confinement imposed by the potential well), resulting in a constant MSD at large time scales.

The diffusive dynamics of the turbulent wake are also depicted in the power spectral density of the CoP, plotted as inset in figure 4. The spectral density is closely related to the MSD: in general, for a power-law behaviour of the MSD, $\langle [{\rm\Delta}\boldsymbol{x}({\it\tau})]^{2}\rangle \propto {\it\tau}^{{\it\kappa}}$ , the asymptotic form of the power spectral density is ${\it\Phi}(f)\propto f^{-(1+{\it\kappa})}$ . In accordance with the linear increase of the MSD with time ( ${\it\kappa}=1$ ), an exponent of $-2$ is observed in the spectrum of the angular component consistent with Brownian motion. A similar decay is observed for the radial component when $f>1/t_{s}$ . However, at low frequencies corresponding to $f\rightarrow 0$ , or equivalently large time scales, it reaches a plateau and levels off, in accordance with the MSD measurements for ${\it\tau}\rightarrow \infty$ .

The dynamics associated with the radial and angular motion of the CoP have been analysed independently so far. Now, we analyse the coupled dynamics predicted by the model. Specifically the coupling arises as an inverse relationship between $\dot{{\it\phi}}={\rm\Delta}{\it\phi}/{\rm\Delta}t$ and $r$ for the angular component, as described by (2.2). The model suggests that the conditional PDF of $\dot{{\it\phi}}$ for a given $r$ follows a Gaussian distribution with variance inversely proportional to $r^{2}$ , that is

(4.2) $$\begin{eqnarray}P(\dot{{\it\phi}}\mid r)=\mathscr{N}\left(0,\frac{{\it\sigma}^{2}}{r^{2}}\right)=\frac{r}{{\it\sigma}\sqrt{2{\rm\pi}}}\exp \left(-\frac{1}{2}\left(\frac{r\dot{{\it\phi}}}{{\it\sigma}}\right)^{2}\right).\end{eqnarray}$$

Indications of the inverse relationship can be found in the time series of the CoP in figure 3. One observes that for small values of the radial component, abrupt changes of the azimuthal orientation occur. The inverse relation is further validated from the joint PDF of the angular variation and radial location, $p(r,{\rm\Delta}{\it\phi})$ , shown in figure 5, where large reorientations are more probable at small radii. The conditional PDFs of the reorientations ${\rm\Delta}{\it\phi}$ , $p({\rm\Delta}{\it\phi}\mid r)=p(r,{\rm\Delta}{\it\phi})/P(r)$ , collapse to a zero-mean Gaussian distribution with variance ${\it\sigma}^{2}$ when scaled by $1/r$ and plotted against $r\dot{{\it\phi}}$ . The latter becomes evident in the inset of figure 5, where we compare the scaled Gaussian PDF with identified variance ${\it\sigma}^{2}$ , obtained from the model given in table 1, and the normalised conditional distributions of the experimental data. In the limit of small radii values, $r\rightarrow 0$ , angular rotations with an almost uniform PDF are observed in the joint PDF. This situation of small $r$ corresponds to large variance and therefore a very broad Gaussian that, over the finite measurable window of $-{\rm\pi}$ to ${\rm\pi}$ , will appear uniform. This is analogous to the cessation events observed by Brown, Nikolaenko & Ahlers (2005) during the reorientation of the large scales of Rayleigh–Bénard convection and can be understood without specifying arbitrary thresholds for a small  $r$ .

Figure 5. Reorientation characteristics: joint PDF of angular variation and radial position of the CoP as measured from the experiment. The inset shows normalised conditional probabilities: experimental data (symbols) collapse on a Gaussian distribution (line) with identified variance ${\it\sigma}^{2}/r^{2}$ , as described by the model.

5. Conclusions

The physical picture that can be drawn for the turbulent axisymmetric wake based on the above is as follows. The laminar large-scale structures, associated with spatially broken symmetries, persist at high Reynolds numbers. In the turbulent regime, these structures undergo diffusive motion (random walk) in a two-dimensional Mexican-hat-shaped potential well (see figure 1), restoring statistically the broken symmetries. In this paper we have shown that this behaviour can be captured by a simple nonlinear model consisting of two coupled stochastic differential equations. The deterministic part of the model accounts for the broken symmetries observed in the laminar regime, and the stochastic part models in a phenomenological sense the turbulent fluctuations acting on the large-scale structures.

The unsteady bifurcation, responsible for the vortex shedding (limit cycle), has not been considered in the present analysis. This can be done by adding an extra equation accounting for the temporal broken symmetries due to a Hopf bifurcation (Fabre et al. 2008; Meliga et al. 2009). Although the stochastic model was demonstrated for a case with high Reynolds number ( $\mathit{Re}_{D}=188\,000$ ), we expect it to be valid over a wide range. An investigation of the Reynolds dependence of the model coefficients will contribute towards a deeper understanding of the dynamic evolution of the large-scale structures observed in chaotic/turbulent wakes and their modelling.

The diffusive dynamics of the large-scale structures presented here and in Rigas et al. (2014) show close similarities with the diffusive dynamics of the large-scale circulation observed in the turbulent Rayleigh–Bénard convection (Brown et al. 2005; Brown & Ahlers 2006) and the turbulent von Kármán swirling flow (de la Torre & Burguete 2007; Berhanu et al. 2009; Pétrélis et al. 2009). However, this similarity is not surprising if the symmetries of the problem are accounted for: spatial rotational symmetry that breaks en route to turbulence.

In summary, we have shown that turbulent dynamics can be described by deterministic equations derived from symmetry arguments with stochastic forcing terms that give rise to turbulent dynamics. The modelling process we have presented here should be applicable to other turbulent flows, provided that their specific symmetries are taken into account.