Numerical simulations of a stochastic dynamics leading to cascades and loss of regularity: applications to fluid turbulence and generation of fractional Gaussian fields

TL;DR

This work introduces a linear stochastic PDE in Fourier space, driven by a delta-correlated forcing, to capture the energy cascade and the emergence of small-scale structure in turbulence. A novel finite-volume discretization in the Fourier domain, coupled with a splitting time integrator, yields accurate, isotropic solutions and reproduces key second-order statistics, including a viscosity-independent velocity variance in the inviscid limit and a Kolmogorov-like PSD $E(t,k) \propto |k|^{-(2H+d)}$. The method is validated in $d=1,2,3$, demonstrating correct PSD scaling and structure-function behavior $\mathbb{E}[|\delta_\ell u|^2] \propto |\ell|^{2H}$ with $H=1/3$, matching turbulence phenomenology and providing a robust computational framework for analyzing stochastic cascade dynamics. The work also outlines avenues for rigorous convergence analysis, extension to vector fields, inverse-transform definitions, intermittency modeling, and kinetic-energy budgeting, highlighting the model’s potential impact on turbulence theory and numerical simulation.

Abstract

Motivated by the modeling of the spatial structure of the velocity field of three-dimensional turbulent flows, and the phenomenology of cascade phenomena, a linear dynamics has been recently proposed able to generate high velocity gradients from a smooth-in-space forcing term. It is based on a linear Partial Differential Equation (PDE) stirred by an additive random forcing term which is delta-correlated in time. The underlying proposed deterministic mechanism corresponds to a transport in Fourier space which aims at transferring energy injected at large scales towards small scales. The key role of the random forcing is to realize these transfers in a statistically homogeneous way. Whereas at finite times and positive viscosity the solutions are smooth, a loss of regularity is observed for the statistically stationary state in the inviscid limit. We here present novel simulations, based on finite volume methods in the Fourier domain and a splitting method in time, which are more accurate than the pseudo-spectral simulations. We show that the novel algorithm is able to reproduce accurately the expected local and statistical structure of the predicted solutions. We conduct numerical simulations in one, two and three spatial dimensions, and we display the solutions both in physical and Fourier spaces. We additionally display key statistical quantities such as second-order structure functions and power spectral densities at various viscosities.

Figures (7)

• Figure 1: (a) Representation of the discretization of the plane spanned by the wave vector $k$ in dimension $d=2$ using finite volumes method, each cell corresponds to a finite volume $\mathcal{K}_{i,a}$ (Eq. \ref{['Meshelements']}) with radial step $\Delta \rho_i$. (b) Representation of a unit cell $\mathcal{K}_{i,a}$ with our notation. We superimpose the non-vanishing radial fluxes (Eq. \ref{['eq:RadialFluxes']}).
• Figure 2: Solution to the dynamics in Fourier space and physical space in the statistically stationary regime. All simulations are conducted with $H=1/3$, $c=1$, $\Delta \rho = h=\kappa=2^{-3}$ and $\widehat{C}_f(k)= \mathds{1}_{\kappa\;\leqslant\; |k| \leqslant k_f}$ (see \ref{['corr-fia']}), with $k_f=4\kappa$. (a) Volume-averaged Fourier mode amplitude $|\hat{u}_{\mathcal{K}_{n,a}}(t)|$ in the statistically stationary regime (i.e. for $t>T^\star$), choosing $\nu=10^{-9}$ and using $N= 2^{12}$ collocation points in the radial direction, as a function of the radial coordinate $\rho_i$\ref{['eq:midpointapproximation']}. (b), (c) Physical space representation of the solution in the statistically stationary regime for $\nu=10^{-5}$ and $10^{-9}$ using correspondingly $N= 2^7$ and $2^{12}$, at a given time in the statistically stationary regime, as a function of the non dimensional variable $xk_f$. The physical space representations are obtained using the inversion formula \ref{['eq:DeftildeUDeltaPhysSpace1D']} over $|x| \leqslant L_{\text{tot}}/2$, with the spatial resolution $\Delta x= 1/k_{\text{max}}$ with $k_{\text{max}}=\kappa+N\Delta \rho$, and the total length $L_{\text{tot}}$ of the physical domain chosen to be $L_{\text{tot}}=1/\Delta \rho$.
• Figure 3: Numerical estimation of second-order statistical quantities for the one dimensional $d=1$ fields: (a) Periodograms $E^{\Theta}_\mathcal{K}$\ref{['eq:DefEKPSDTheta']} as a function of the radial coordinate $\rho_i$ and (b) Second-order structure functions $\left\langle \left(\delta_\ell \widetilde{u}_{\Delta}\right)^2\right\rangle$\ref{['eq:DefS2tildeUDelta']}. In both figures, the representation is made in a logarithmic fashion, and the darker the curve, the lower the viscosity. These quantities have been estimated while averaging over $10^3$ instances in the statistically stationary regime, every $10$ time units. All the simulations are conducted with the same parameters $h$, $c$, $H$, $h=\kappa$, $k_f$ and $\widehat{C}_f(k)$ as they are given in the caption of Fig. \ref{['fig:Phys1D']}. Values of viscosity correspond to, from lighter to darker, $\nu=10^{-5},10^{-6},10^{-7},10^{-8},10^{-9}$, with corresponding number of collocation points along the radial direction $N=2^{7},2^{8},2^{9},2^{10},2^{11}$. With dashed lines, we superimpose the theoretical predictions of the power-law behaviors in (a) based on \ref{['eq:ExpEKPSDGivenCf']} with the particular value $d=1$, and in (b) based on \ref{['eq:ComputSecondOrderSFModelDoubleLimitSmallScales']}, with geometrical factor $c_1$\ref{['eq:PredContCdS2D1']}.
• Figure 4: Snapshot of the solution for the 2D dynamics, in (a) Fourier and (b) physical spaces, at a time pertaining to the statistically stationary regime. In both cases, we have used the following parameters: $\nu=10^{-5}$, $N=2^{7}$, $H=1/3$, $c=1$, $h=2^{-7}$, $N_\vartheta=2^9$, $\kappa=1$, and $\widehat{C}_f(k)= \mathds{1}_{\kappa\;\leqslant\; |k| \leqslant k_f}$ (see \ref{['corr-fia']}), with $k_f=\kappa+3h$. Notice that to get the inverse Fourier transform $\widetilde{u}_{\Delta}(t,x,y)$, represented in (b), based on the modes $\hat{u}_{\mathcal{K}}(t)$ displayed in (a), we have used the inversion formula based on \ref{['eq:DeftildeUDeltaPhysSpace2D']}.
• Figure 5: Numerical estimations of the second-order statistical quantities for the two dimensional case ($d=2$). In both figures, darker the curves, lower viscosities. These quantities are computed in the statistically stationary regime with $500$ instances of the corresponding fields, every $10$ units of time. All the simulations are conducted with $H=1/3$, $c=1$, $h=2^{-7}$, $N_\vartheta=2^9$, $\kappa=1$, $k_f=\kappa+3h$ and the same $\widehat{C}_f(k)$ used in Fig. \ref{['fig:Phys2D']}. Chosen values of viscosities are $\nu=10^{-5},10^{-6},10^{-7},10^{-8},10^{-9}$ with, respectively $N=2^{10},2^{11},2^{12},2^{13},2^{14}$ collocation points in the radial direction. (a) Angle averaged periodograms $E_\mathcal{K}^\theta$\ref{['eq:DefEKPSDTheta']}, as a function of the radial coordinate $\rho_i$, weighted by the corresponding volume of unit cells $|\mathcal{K}_{i,a}|$ (which is independent of the angle coordinate $a$). (b) Second-order structure functions $\left\langle \left(\delta_\ell \widetilde{u}_{\Delta}(t,x) \right)^2\right\rangle$\ref{['eq:DefS2tildeUDelta']} for different values of the viscosity (solid line). We superimpose with dotted lines in (a) the precise asymptotic power-law behavior given in \ref{['eq:ExpEKPSDGivenCf']}, and in (b) the predicted asymptotic power-law based on \ref{['eq:ComputSecondOrderSFModelDoubleLimitSmallScales']} with corresponding geometrical factor $c_2$\ref{['eq:PredContCdS2D2']}. In (b), we also indicate as a guide to the eyes the dissipative behavior $\ell^2$.