(This content is originally written by ‪Kyongmin Yeo’s manual)


To maintain stationary turbulence, external force term is added to Navier-Stokes equation \[ \dfrac{d \hat{u}_i}{dt} = - i \kappa_{i} \hat{P} + \hat{H}_i - \nu \kappa^2 \hat{u}_i + \hat{f}_i \]

where \(\hat{f}_i\) is a external forcing term.

\(\hat{f}_i\) is defined as the projection of a vector \(\hat{\mathbf{b}}\) onto the plane normal to the wavenumber vector \(\mathbf{\kappa}\) to ensure divergence-free condition.

\[ \hat{f}_{i} = \hat{b}_{i} - \dfrac{\kappa_i}{\kappa^2} \kappa_j \hat{b}_{j} \]

So, how we define vector \(\hat{\mathbf{b}}\)? Eswaran & Pope suggested stochastic forcing. They define 3D complex vector \(\hat{\mathbf{b}}\) which is non-zero in the range \(0 < \kappa < \kappa_f \), in which \(\kappa_f\) is the maximum forcing wavenumber. This can be interpreted as forcing to sphere in wavenumber space.

\(\hat{\mathbf{b}}\) is composed of six independent Uhlenbeck-Ornstein process.

\[ \hat{\mathbf{b}} = \begin{bmatrix} UO1 \newline UO3 \newline UO5 \newline \end{bmatrix} + i \begin{bmatrix} UO2 \newline UO4 \newline UO6 \newline \end{bmatrix} \]

Solving Uhlenbeck-Ornstein process

Each stochastic process, \(UO1 ~ UO6\), is chosen so as to satisfy hte Langevin equation with a time scale \(T^f_L\) and stadnard deviation \(\sigma_f\),

\[ dUO = - \dfrac{UO}{T^f_L} \Delta t + \left( \dfrac{2\sigma^2_f}{T^f_L} \right)^{1/2} dW_t \]

in which \(W_t\) denotes a Wiener process satisfying \[ dW_t \sim \mathcal{N} (0, \Delta t) \]

Analytical solution to Langevin equation is given by

\[ UO(t) = UO(0) e^{-t/T^f_L} + e^{-t/T^f_L} \int^t_{0} e^{s/T^f_L} (2\sigma^2_f/T^f_L)^{1/2}dW_s \]

Above solution can be solved discretely by applying Itô integral. With RK3 method, UO process discretized solution is

\[ UO^{n+1} = e^{-(a_n+b_n)\Delta t / T^f_L}\left[ UO^{n} + e^{s/T^f_L} (2\sigma^2_f/T^f_L)^{1/2}dW_s dW^n \right] \]

in which discretized Wiener proces is

\[ dW^n \sim \mathcal{N} (0, (a_n + b_n) \Delta t) \]

Estimating Reynolds Number

Given parameters

  • \(\nu\) : Fluid viscosity
  • \(\beta\) : constant (\(=0.8\))
  • \(\kappa_0\) : smallest wavenumber
  • \(\kappa_f\) : maximum forcing wavenumber
  • \(T_L\) : Forcing time scale
  • \(\epsilon^* \equiv \sigma^2 T_L \) where \(\sigma\) is a forcing amplitude, usually just given by constant


  • \(\epsilon \propto N_f \epsilon^*\) \begin{align} \epsilon &= \dfrac{4\epsilon^* T_e N_f}{T_L + T_e} \newline T_e &= \dfrac{\beta}{(N_f \epsilon^* \kappa^2_{0})^{1/3}} \end{align}
  • \(\kappa^{-1}_{0}\) : Integral length scales

Computed parameters

  • \(N_f\) : The number of forced modes, \(\kappa < \kappa_f>\), counted manually
  • Predicted value of the energy dissipation, replace \(T_e\) by \(\beta / (N_f \epsilon^* \kappa^2_{0})^{1/3}\) \[ \epsilon^*_{T} \equiv \dfrac{4 \epsilon^* N_f}{1 + T_{L} N^{1/3}_{F}/\beta} \]
  • Predicted Kolmogorov microscale \[ \eta_T \equiv (\nu^3 / \epsilon_T) \]

Predicted \(Re\)

Using above parameters Taylor Reynolds number is estimated by

\[ Re \simeq \dfrac{8.5}{(\eta_T \kappa_0)^{5/6} N^{2/9}_{F}} \]


  • More detail explanation of Uhlenbeck-Ornstein process
  • Add code for computing \(Re\)