3.1. Definitions

We wish to model turbulent velocity fields $\boldsymbol {U}(\boldsymbol {x}) = \langle \boldsymbol {U}(\boldsymbol {x})\rangle + \boldsymbol {u}(\boldsymbol {x}) \in \mathbb {R}^3$. Here, $\langle \boldsymbol {U}\rangle = (\langle U_1\rangle ,\langle U_2\rangle ,\langle U_3\rangle )$ is the mean velocity field and $\boldsymbol {u} = (u_1,u_2,u_3)$ (sometimes also written $(u,v,w)$) are the zero-mean turbulent fluctuations. All of the models we choose to consider for $\boldsymbol {u}$ are Gaussian. That is, they are determined entirely from the two-point correlation tensor

(3.1)\begin{equation} R_{ij}(\boldsymbol{r},\boldsymbol{x},t) = \langle u_i(\boldsymbol{x},t)u_j(\boldsymbol{x}+\boldsymbol{r},t)\rangle . \end{equation}

When $R(\boldsymbol {r},\boldsymbol {x},t) = R(\boldsymbol {r},t)$ depends only on the separation vector $\boldsymbol {r}$, the model is said to be spatially homogeneous. Alternatively, when $R(\boldsymbol {r},\boldsymbol {x},t) = R(\boldsymbol {r},\boldsymbol {x})$ is independent of the time variable $t$, the model is said to be temporally stationary.

Frequently, it is convenient to consider the Fourier transform of the velocity field $\boldsymbol {u}$. In such cases, we express the field in terms of a generalized Fourier–Stieltjes integral,

(3.2)\begin{equation} \boldsymbol{u}(\boldsymbol{x}) = \int_{\mathbb{R}^3} \operatorname{e}^{\operatorname{i}\boldsymbol{k}\boldsymbol{\cdot}\boldsymbol{x}} \,\mathrm{d} \boldsymbol{Z}(\boldsymbol{k}), \end{equation}

where $\boldsymbol {Z}(\boldsymbol {k})$ is a three-component measure on $\mathbb {R}^3$. The validity of this expression follows from the Wiener–Khinchin theorem (Lord, Powell & Shardlow Reference Lord, Powell and Shardlow2014). Likewise, in the homogeneous setting, we may consider the Fourier transform of the covariance tensor, otherwise known as the velocity-spectrum tensor,

(3.3)\begin{equation} \varPhi_{ij}(\boldsymbol{k},t) = \frac{1}{(2{\rm \pi})^3}\int_{\mathbb{R}^3} \operatorname{e}^{-\operatorname{i}\boldsymbol{k}\boldsymbol{\cdot}\boldsymbol{r}}R_{ij}(\boldsymbol{r},t)\,\mathrm{d} \boldsymbol{r} . \end{equation}

Consider three-dimensional additive white Gaussian noise (Hida et al. Reference Hida, Kuo, Potthoff and Streit2013; Kuo Reference Kuo2018) in the physical and frequency domains, denoted $\boldsymbol {\xi }(\boldsymbol {x})$ and $\hat {\boldsymbol {\xi }}(\boldsymbol {k})$, respectively, such that

(3.4)\begin{equation} \boldsymbol{\xi}(\boldsymbol{x}) = \int_{\mathbb{R}^3} \operatorname{e}^{\operatorname{i}\boldsymbol{k}\boldsymbol{\cdot}\boldsymbol{x}} \hat{\boldsymbol{\xi}}(\boldsymbol{k}) \,\mathrm{d} \boldsymbol{k} = \int_{\mathbb{R}^3} \operatorname{e}^{\operatorname{i}\boldsymbol{k}\boldsymbol{\cdot}\boldsymbol{x}} \,\mathrm{d} \boldsymbol{W}(\boldsymbol{k}), \end{equation}

where $\boldsymbol {W}(\boldsymbol {k})$ is three-dimensional Brownian motion. We assume $\mathrm {d} \boldsymbol {Z}(\boldsymbol {k}) = \boldsymbol{\mathsf{G}}(\boldsymbol {k}) \,\mathrm {d} \boldsymbol {W}(\boldsymbol {k}) = \boldsymbol{\mathsf{G}}(\boldsymbol {k}) \hat {\boldsymbol {\xi }}(\boldsymbol {k}) \,\mathrm {d} \boldsymbol {k}$, where $\boldsymbol{\mathsf{G}}(\boldsymbol {k})^\ast \boldsymbol{\mathsf{G}}(\boldsymbol {k}) = \varPhi (\boldsymbol {k})$.

This section and the next are devoted to deriving FPDE models for homogeneous turbulence. The approach we follow involves a commonly used definition of fractional differential operators facilitated by the spectral theorem (Reed Reference Reed2012). Note that, for an abstract closed normal operator $A\colon \mathcal {D}(A)\subseteq H\to H$ on a complex Hilbert space $H$, $AA^\ast = A^\ast A$, there exists a finite measure space $(Y, \mu )$, together with a complex-valued measurable function $\lambda (y)$, defined on $Y$, and a unitary map $U : H\to L_2(Y, \mu )$, such that

(3.5)\begin{equation} UA\phi= \lambda U\phi \quad \text{for all}\ \phi\in H . \end{equation}

In this case, one may define the $\alpha$-fractional power of $A$ as follows:

(3.6)\begin{equation} A^\alpha= U^\ast \lambda^\alpha U . \end{equation}

For an operator $A\colon \mathcal {D}(A)\subseteq L^2(\varOmega )\to L^2(\varOmega )$ with a discrete spectrum, we may simply write

(3.7)\begin{equation} A^{\alpha}\phi = \sum_{j=1}^\infty \lambda_j^{\alpha} (\phi,e_j)_\varOmega\, e_j . \end{equation}

Here, $e_j$ and $\lambda _j$ denote the corresponding eigenmodes and eigenvalues of $A$ and $(\phi ,\chi )_\varOmega = \int _\varOmega \phi \boldsymbol {\cdot }\chi \,\mathrm {d} x$ denotes the $L^2$-inner product on the domain $\varOmega \subseteq \mathbb {R}^3$.

For example, consider the vector Laplacian operator $A = -\varDelta$ on $\varOmega = \mathbb {R}^d$. Letting $k = |\boldsymbol {k}|$ denote the magnitude of the wavenumber vector $\boldsymbol {k} = (k_1,k_2,k_3)$ in Fourier space and $\mathcal {F}$ and $\mathcal {F}^{-1}$ denote the Fourier and inverse Fourier transforms, respectively, we have

(3.8)\begin{equation} (-\varDelta)^{\alpha}\boldsymbol{\phi}(\boldsymbol{x}) = \frac{1}{(2{\rm \pi})^d} \int_{\mathbb{R}^d} k^{2\alpha} (\boldsymbol{\phi},\operatorname{e}^{-\operatorname{i}\boldsymbol{k}\boldsymbol{\cdot} \boldsymbol{x}})_{\mathbb{R}^d}\,\operatorname{e}^{\operatorname{i}\boldsymbol{k}\boldsymbol{\cdot} \boldsymbol{x}} \,\mathrm{d} \boldsymbol{k} = \mathcal{F}^{{-}1}\{k^{2\alpha}\mathcal{F}\{\boldsymbol{\phi}\}(\boldsymbol{k})\}(\boldsymbol{x}). \end{equation}

Evidently, in this setting, $\mathcal {F}$ is the analogue of the unitary operator $U$ present in the abstract expression (3.6). On the other hand, when $\varOmega = (0,1)^d$ is a periodic domain, it is well known that $A = -\varDelta$ has a discrete spectrum. Here, recall that

(3.9)\begin{equation} (-\varDelta)^{\alpha}\boldsymbol{\phi}(\boldsymbol{x}) = \frac{1}{(2{\rm \pi})^d} \sum_{\boldsymbol{j}\in\mathbb{Z}^d} k_{\boldsymbol{j}}^{2\alpha} {(\boldsymbol{\phi},\operatorname{e}^{-\operatorname{i}\boldsymbol{k}_{\boldsymbol{j}}\boldsymbol{\cdot} \boldsymbol{x}})}_{\mathbb{R}^d}\,\operatorname{e}^{\operatorname{i}\boldsymbol{k}_{\boldsymbol{j}}\boldsymbol{\cdot} \boldsymbol{x}}. \end{equation}

For further details on the spectral representation of closed operators, we refer the interested reader to de Dormale & Gautrin (Reference de Dormale and Gautrin1975), Weidmann (Reference Weidmann2012) and Kowalski (Reference Kowalski2009).

3.2. The von Kármán model

Let us begin with a standard form of the spectral tensor used in isotropic stationary and homogeneous turbulence models, namely,

(3.10)\begin{equation} \varPhi_{ij}(\boldsymbol{k}) = (4{\rm \pi})^{{-}1}k^{{-}2} E(k) P_{ij}(\boldsymbol{k}). \end{equation}

Here, $E(k)$ is called the energy spectrum function and $P_{ij}(\boldsymbol {k}) = \delta _{ij} - {k_ik_j}/{k^2}$ is commonly referred to as the projection tensor. One common empirical model for $E(k)$, suggested by Von Kármán (Reference Von Kármán1948), is given by the expression

(3.11)\begin{equation} E(k) = c_0^2 \varepsilon^{2/3}k^{{-}5/3} \left(\frac{k L}{(1 + (k L)^2)^{1/2}}\right)^{17/3} . \end{equation}

Here, $\varepsilon$ is the viscous dissipation of the turbulent kinetic energy, $L$ is a length scale parameter and $c_0^2\approx 1.7$ is an empirical constant.

Recall that the Fourier transform of the scalar Laplacian is simply $-k^2$. Likewise, consider the Fourier transform $\boldsymbol{\mathsf{Q}}(\boldsymbol {k})$ of the $\mathrm {curl}$ operator, $\int _{\mathbb {R}^3} \boldsymbol {\nabla } \times \boldsymbol {v}(\boldsymbol {r}) \operatorname {e}^{-\operatorname {i}\boldsymbol {k}\boldsymbol {\cdot }\boldsymbol {r}}\,\mathrm {d}\boldsymbol {r} = \boldsymbol{\mathsf{Q}}(\boldsymbol {k})\hat {\boldsymbol {v}}(\boldsymbol {k})$, where $\hat {\boldsymbol {v}}(\boldsymbol {k}) = \int _{\mathbb {R}^3} \boldsymbol {v}(\boldsymbol {r}) \operatorname {e}^{-\operatorname {i}\boldsymbol {k}\boldsymbol {\cdot }\boldsymbol {r}}\,\mathrm {d}\boldsymbol {r}$. Observe that

(3.12)\begin{equation} \boldsymbol{\mathsf{Q}}(\boldsymbol{k}) = \operatorname{i} \begin{bmatrix} 0 & -k_3 & k_2\\ k_3 & 0 & -k_1\\ -k_2 & k_1 & 0 \end{bmatrix} \end{equation}

and, moreover, $P(\boldsymbol {k}) = k^{-2}\boldsymbol{\mathsf{Q}}(\boldsymbol {k})^\ast \boldsymbol{\mathsf{Q}}(\boldsymbol {k})$. Motivated by the decomposition $\varPhi (\boldsymbol {k}) = \boldsymbol{\mathsf{G}}(\boldsymbol {k})^\ast \boldsymbol{\mathsf{G}}(\boldsymbol {k})$, we choose to simply write $\boldsymbol{\mathsf{G}}(\boldsymbol {k}) = ({1}/{\sqrt {4{\rm \pi} }}) k^{-2} E^{1/2}(k) \boldsymbol{\mathsf{Q}}(\boldsymbol {k})$. Next, recalling $\mathrm {d} \boldsymbol {Z}(\boldsymbol {k}) = \boldsymbol{\mathsf{G}}(\boldsymbol {k})\, \mathrm {d} \boldsymbol {W}(\boldsymbol {k})$, it immediately follows that

(3.13)\begin{equation} \mathrm{d} \boldsymbol{Z}(\boldsymbol{k}) = \boldsymbol{\mathsf{Q}}(\boldsymbol{k}) \left(\frac{1}{\sqrt{4{\rm \pi}}k^2} E^{1/2}(k) \,\mathrm{d}\boldsymbol{W}(\boldsymbol{k})\right) . \end{equation}

Integrating both sides with respect to $\boldsymbol {k}$, we arrive at the expression $\boldsymbol {u} = \boldsymbol {\nabla } \times \boldsymbol {\psi }$, with a vector potential defined

(3.14)\begin{equation} \boldsymbol{\psi}(\boldsymbol{x}) = \frac{1}{\sqrt{4{\rm \pi}}} \int_{\mathbb{R}^3} k^{{-}2} E^{1/2}(k) \operatorname{e}^{\operatorname{i}\boldsymbol{k}\boldsymbol{\cdot}\boldsymbol{x}} \,\mathrm{d} \boldsymbol{W}(\boldsymbol{k}). \end{equation}

We now proceed to relate the vector potential $\boldsymbol {\psi }(\boldsymbol {x})$ to the solution of a FPDE. Writing $\boldsymbol {\psi }(\boldsymbol {x}) = \int _{\mathbb {R}^3} \operatorname {e}^{\operatorname {i}\boldsymbol {k}\boldsymbol {\cdot }\boldsymbol {x}} \,\mathrm {d} \boldsymbol {Y}(\boldsymbol {k})$, similar to (3.2), and rearranging the factors in (3.14), leads to

(3.15)\begin{equation} (1 + (k L)^2)^{17/12} \,\mathrm{d} \boldsymbol{Y}(\boldsymbol{k}) = c_0\varepsilon^{1/3}L^{17/6}\,\mathrm{d} \boldsymbol{W}(\boldsymbol{k}) . \end{equation}

Then, upon integrating both sides with respect to $\boldsymbol {k}$, we arrive at the FPDE

(3.16)\begin{equation} (I-L^2\varDelta)^{17/12} \boldsymbol{\psi} = c_0 \varepsilon^{1/3} L^{17/6} \boldsymbol{\xi} . \end{equation}

This and all future differential equations are only properly understood in the sense of distributions, yet we continue to use the ‘strong form’ for readability.

Let $I$ denote the identity operator, $A = I-L^2\varDelta$, $\mu = c_0 \varepsilon ^{1/3}$, and $\alpha = 17/12$. With these symbols in hand, the derivation above can be summarized as follows:

(3.17)\begin{equation} \boldsymbol{u} = \operatorname{\boldsymbol{\nabla}\times}\boldsymbol{\psi}, \quad \text{where}\ A^{\alpha} \boldsymbol{\psi} = \mu L^{2\alpha} \boldsymbol{\xi} . \end{equation}

In the next section, we extend the simple FPDE model above in order to describe inhomogeneous turbulence on bounded domains. This is achieved by both generalizing the definition of the length scale $L$ and the fractional operator $A^{\alpha }$ as well as introducing a physical notion of boundary conditions.

Remark 3.1 Note that the vector potential $\boldsymbol {\psi }(\boldsymbol {x})$, defined in (3.14), is not divergence-free. In an alternative model, one may seek to enforce this condition. In this case, one would naturally arrive at the Stokes-type system

(3.18)\begin{equation} A^{\alpha} \boldsymbol{\psi} + \boldsymbol{\nabla} \phi = \mu L^{2\alpha} \boldsymbol{\xi} , \quad \operatorname{\boldsymbol{\nabla}\hspace{1.25pt}\boldsymbol{\cdot}} \boldsymbol{\psi} = 0 . \end{equation}

Here, $\phi$ plays the role of an additional pressure-like Lagrange multiplier. Note that by taking the curl of the first equation above, the turbulence $\boldsymbol {u}(\boldsymbol {x})$ can be characterized by just one equation; namely,

(3.19)\begin{equation} A^{\alpha} \boldsymbol{u} = \mu L^{2\alpha} \operatorname{\boldsymbol{\nabla}\times} \boldsymbol{\xi}. \end{equation}

For the sake of completeness, note that we may also define a generalized vorticity field $\boldsymbol {w} = -\varDelta \boldsymbol {\psi }$. One may show that $\boldsymbol {w}(\boldsymbol {x}) = ({1}/{\sqrt {4{\rm \pi} }}) \int E^{1/2}(k) \operatorname {e}^{\operatorname {i}\boldsymbol {k}\boldsymbol {\cdot }\boldsymbol {x}} \,\mathrm {d} \boldsymbol {W}(\boldsymbol {k})$. This expression, in combination with the partial differential equation (PDE)

(3.20)\begin{equation} -\varDelta \boldsymbol{u} = \operatorname{\boldsymbol{\nabla}\times} \boldsymbol{w}, \end{equation}

can also be used to characterize $\boldsymbol {u}(\boldsymbol {x})$.

Both (3.19) and (3.20) are perfectly valid and equivalent characterizations of the homogeneous turbulent velocity field considered above, $\boldsymbol {u}(\boldsymbol {x})$, on the free space domain $\mathbb {R}^3$. More importantly, they will likely lead to alternative turbulence models on more complicated domains, once appropriate boundary conditions are selected. We have chosen not to use (3.19) because it is not valid in the presence of non-homogeneous length scales $L = L(\boldsymbol {x})$; a modelling consideration we wish to incorporate. The non-homogeneous setting still requires the saddle-point problem (3.18) in order to enforce volume conservation in $\boldsymbol {\psi }(\boldsymbol {x})$. Because $\boldsymbol {u} = \operatorname {\boldsymbol {\nabla }\times } \boldsymbol {\psi }$ does not depend on the irrotational part of $\boldsymbol {\psi }(\boldsymbol {x})$, equation (3.18) appears to be a valid alternative model which we leave open for future investigation. Finally, we have chosen to avoid (3.20) because of the low regularity of the solution variable $\boldsymbol {w}(\boldsymbol {x})$; cf. figure 1.

Figure 1. Normalized magnitudes of (a) $\boldsymbol {\psi }$, (b) $\boldsymbol {u} = \operatorname {\boldsymbol {\nabla }\times }\boldsymbol {\psi }$ and (c) $\boldsymbol {w} = -\varDelta \boldsymbol {\psi }$. Observe the decrease of regularity, from left to right, with higher-order derivatives of the vector potential. The fields are computed using a discrete Fourier transform.

4.1. A general class of inhomogeneous models

Equation (3.16) was derived from a very specific form of the energy spectrum function $E(k)$. Under the same decomposition of the spectral tensor $\varPhi (\boldsymbol {x})$ given in (3.10), a much more general family of homogeneous and stationary random field models derive from the following ansatz on the energy spectrum function:

(4.1)\begin{equation} k^{{-}4}E(\boldsymbol{k}) = \mu^2 \det(\bar{\boldsymbol{\varTheta}})^{2/3\gamma}(1 +\boldsymbol{k}^\top\bar{\boldsymbol{\varTheta}}\boldsymbol{k})^{{-}2\alpha_1}(\boldsymbol{k}^\top\bar{\boldsymbol{\varTheta}}\boldsymbol{k})^{{-}2\alpha_2} . \end{equation}

Here, $\bar {\boldsymbol {\varTheta }} \in \mathbb {R}^{3\times 3}$ is a fixed symmetric positive-definite matrix and $\alpha _2$, $\alpha _1$, $\gamma$ and $\mu$ are additional scalar parameters.

Equation (4.1) is a broad generalization of (3.11) which replaces $E$ as function of $k = |\boldsymbol {k}|$ by $E$ as function of $\boldsymbol {k}$. This flexibility allows us, for example, to consider anisotropic effects. Indeed, just as $L$ played the role of a length scale in (3.11), here, $\bar {\boldsymbol {\varTheta }}$ plays the role of a metric in Fourier space. In addition, observe that if $\bar {\boldsymbol {\varTheta }}=L^2\boldsymbol{\mathsf{I}}$, where $\boldsymbol{\mathsf{I}}$ denotes the identity matrix, $4\alpha _2=4-p_0$, $4\alpha _1=5/3+p_0$, $\gamma = \alpha _1+\alpha _2$ and $\mu ^2=C \varepsilon ^{2/3}$, then (4.1) reproduces the following common one-parameter homogeneous energy spectrum model (see, e.g., Pope Reference Pope2001, p. 232):

(4.2)\begin{equation} E(k) = C \varepsilon^{2/3}k^{{-}5/3} \left(\frac{k L}{((k L)^2 + 1)^{1/2}} \right)^{5/3 + p_0}. \end{equation}

In this scenario, $p_0=4$ corresponds exactly to the von Kármán spectrum (3.11) considered previously, i.e. $\alpha _1 = \gamma = {17}/{12}$ and $\alpha _2 = 0$.

As in (3.14), the vector potential $\boldsymbol {\psi }(\boldsymbol {x}) = \int \operatorname {e}^{\operatorname {i}\boldsymbol {k}\boldsymbol {\cdot }\boldsymbol {x}} \,\mathrm {d} \boldsymbol {Y}(\boldsymbol {k})$ can also be written in terms of a Fourier–Stieltjes integral, weighted by $k^{-2} E^{1/2}(\boldsymbol {k})$. After rearranging factors, equation (4.1) characterizes the vector potential $\boldsymbol {\psi }$ as the solution of the following fractional stochastic PDE on $\mathbb {R}^3$:

(4.3)\begin{equation} (I-\boldsymbol{\nabla}\boldsymbol{\cdot}(\bar{\boldsymbol{\varTheta}}\boldsymbol{\nabla}))^{\alpha_1} (-\boldsymbol{\nabla}\boldsymbol{\cdot}(\bar{\boldsymbol{\varTheta}}\boldsymbol{\nabla}))^{\alpha_2} \boldsymbol{\psi} = \mu\det(\bar{\boldsymbol{\varTheta}})^{\gamma/3}\boldsymbol{\xi}. \end{equation}

Two immediate modifications of (4.3) are now in order. First, we may replace the constant matrix $\bar{\boldsymbol{\varTheta}}$ by a spatially varying metric tensor $\boldsymbol {\varTheta }(\boldsymbol {x})$. This change immediately induces an inhomogeneous turbulence model. Second, we may consider substituting the white noise random variable $\boldsymbol {\xi }$ for a well-chosen coloured noise variable denoted ${\boldsymbol {\eta }}$. Together, these two generalizations lead to a family of random field models written

(4.4)\begin{equation} \left(I-\boldsymbol{\nabla}\boldsymbol{\cdot}(\boldsymbol{\varTheta}(\boldsymbol{x})\boldsymbol{\nabla})\right)^{\alpha_1}\left(-\boldsymbol{\nabla}\boldsymbol{\cdot}(\boldsymbol{\varTheta}(\boldsymbol{x})\boldsymbol{\nabla})\right)^{\alpha_2}\boldsymbol{\psi} = \mu\det(\boldsymbol{\varTheta}(\boldsymbol{x}))^{\gamma/3}{\boldsymbol{\eta}} . \end{equation}

Physically, the metric tensor $\boldsymbol {\varTheta }(\boldsymbol {x})$ introduces inhomogeneous and anisotropic diffusion; this corresponds to local changes of the turbulence length scales which may result from complicated dynamics of interacting eddies. Statistically, it incorporates the possibility for spatially varying correlation lengths and also may contain distortion.

In order to motivate one possible choice of the stochastic forcing term ${\boldsymbol {\eta }}$, note that (4.2) can adequately characterize both the energy-containing and inertial subranges, however, it fails in the dissipative range; namely, where $k$ is large. In order to fit the dissipative range, one approach is to scale the energy spectrum with a decaying exponential function (see, e.g., Pope Reference Pope2001, p. 233):

(4.5)\begin{equation} E_{\beta}(k) = E(k)\operatorname{e}^{-\beta k}, \end{equation}

where $\beta >0$ is a positive constant, usually close to the Kolmogorov length scale. In such scenarios, we suggest using the following definition for ${\boldsymbol {\eta }}$ in (4.4) (see Stein & Weiss Reference Stein and Weiss1971, Theorem 1.14):

(4.6)\begin{equation} \boldsymbol{\xi}_\beta(\boldsymbol{x}) = \int_{\mathbb{R}^3} \operatorname{e}^{\operatorname{i}\boldsymbol{k}\boldsymbol{\cdot}\boldsymbol{x}} \operatorname{e}^{-\beta k} \,\mathrm{d} \boldsymbol{W}(\boldsymbol{k}) \propto \boldsymbol{\xi}(\boldsymbol{x})\ast\frac{\beta}{(\beta^2 + |\boldsymbol{x}|^2)^{2}}, \end{equation}

which converges to (3.4) as $\beta \to 0$. In the presence of shear, a different time-dependent modification is also natural to consider from the point of view of rapid distortion theory. That is the subject of the following subsection.

Remark 4.1 When $\alpha _2$ and $\alpha _1$ are chosen to match the energy spectrum model (4.2), it is clear that $\alpha _2+\alpha _1 = 17/12$ is independent of $p_0$. Under this constraint, $\alpha _2$ mainly affect the behaviour of the power spectrum at the origin and, likewise, the large-scale structure of $\boldsymbol {u}$. In other words, the shape of the spectrum in the inertial subrange is unaffected by the precise choice of $\alpha _2$ and $\alpha _1 = 17/12-\alpha _2$; only the shape of the spectrum in the energy-containing range is affected (see figure 2). An illustration of some energy spectra possibilities is included in figure 2.

Figure 2. The energy spectrum $E_{\beta }(k)/\mu ^2$ corresponding to (4.1) with $\bar {\boldsymbol {\varTheta }}=L^2\boldsymbol{\mathsf{I}}$. The sum $\alpha _2+\alpha _1$ is fixed to $17/12$, which guarantees the slope $k^{-5/3}$ in the inertial subrange. Different values of $\alpha _2$ control the energy-containing range. And the exponent $\beta /L=10^{-3}$ defines the exponential decay in the dissipative subrange.

4.2. A model for shear flows

Consider the velocity field $\boldsymbol {U} = \langle \boldsymbol {U}\rangle + \boldsymbol {u}$ and define the average total derivative of the turbulent fluctuations $\boldsymbol {u} = (u_1,u_2,u_3)$ as follows:

(4.7)\begin{equation} \frac{\bar{D} u_i}{\bar{D} t} = \frac{\partial u_i}{\partial t} + \langle U_j\rangle\frac{\partial u_i}{\partial x_j}. \end{equation}

The rapid distortion equations (see, e.g. Townsend Reference Townsend1980; Maxey Reference Maxey1982; Hunt & Carruthers Reference Hunt and Carruthers1990) are a linearization of the Navier–Stokes equations in free space when the turbulence-to-mean-shear time scale ratio is arbitrarily large. They can be written

(4.8a,b)\begin{equation} \frac{\bar{D} u_i}{\bar{D} t} ={-}u_i \frac{\partial \langle U_j\rangle}{\partial x_i} - \frac{1}{\rho}\frac{\partial p}{\partial x_i}, \quad \frac{1}{\rho} \varDelta p ={-}2\frac{\partial \langle U_i\rangle}{\partial x_j}\frac{\partial u_j}{\partial x_i}. \end{equation}

Under a uniform shear mean velocity gradient, $\langle U_i(\boldsymbol {x})\rangle = x_j\partial \langle U_i\rangle / \partial x_j$, where $\partial \langle U_i\rangle / \partial x_j$ is a constant tensor, a well-known form of these equations can be written out in Fourier space. In this case, the rate of change of each frequency $\boldsymbol {k}(t) = (k_1(t),k_2(t),k_3(t))$ is defined ${\mathrm {d} k_i}/{\mathrm {d} t} = -k_j{\partial \langle U_j\rangle }/{\partial x_i}$. We then have the following Fourier representation of the average total derivative of $\boldsymbol {u}$:

(4.9)\begin{equation} \frac{\bar{D} u_i}{\bar{D} t} = \int_{\mathbb{R}^3} \operatorname{e}^{\operatorname{i}\boldsymbol{k}\boldsymbol{\cdot}\boldsymbol{x}}\left( \left(\frac{\partial }{\partial t} + \frac{\mathrm{d} k_j}{\mathrm{d} t}\frac{\partial }{\partial k_j}\right) \,\mathrm{d} Z_i(\boldsymbol{k},t) \right) = \int_{\mathbb{R}^3} \operatorname{e}^{\operatorname{i}\boldsymbol{k}\boldsymbol{\cdot}\boldsymbol{x}}\left(\frac{\bar{D} \,\mathrm{d} Z_i(\boldsymbol{k},t)}{\bar{D} t}\right) . \end{equation}

With this expression, the Fourier representation of (4.8a,b) can be written

(4.10)\begin{equation} \frac{\bar{D} \,\mathrm{d} Z_j(\boldsymbol{k},t)}{\bar{D} t} = \frac{\partial U_\ell}{\partial x_k} \left( 2 \frac{k_j k_\ell}{k^2} - \delta_{j\ell} \right) \,\mathrm{d} Z_k(\boldsymbol{k},t). \end{equation}

Exact solutions to (4.10) are well-known (see, e.g., Townsend Reference Townsend1980; Mann Reference Mann1994), given the initial conditions $\boldsymbol {k}_0 = (k_{10},k_{20},k_{30})$ and $\mathrm {d} \boldsymbol {Z}(\boldsymbol {k}_0,0)$. In the scenario

(4.11)\begin{equation} \langle \boldsymbol{U}(\boldsymbol{x})\rangle = (U_0 + S x_3) \boldsymbol{e}_1 , \end{equation}

the solution can be written in terms of the evolving Fourier modes $\boldsymbol {k}(t)$ and non-dimensional time $\tau = S t$, as follows:

(4.12)\begin{equation} \mathrm{d} \boldsymbol{Z}(\boldsymbol{k},t) = \boldsymbol{\mathsf{D}}_\tau(\boldsymbol{k}) \mathrm{d} \boldsymbol{Z}(\boldsymbol{k}_0,0) , \end{equation}

where

(4.13a–c)\begin{equation} \boldsymbol{\mathsf{D}}_\tau(\boldsymbol{k}) = \begin{bmatrix} 1 & 0 & \zeta_1\\ 0 & 1 & \zeta_2\\ 0 & 0 & \zeta_3 \end{bmatrix} , \quad \boldsymbol{k}_0 = \boldsymbol{\mathsf{T}}_\tau \boldsymbol{k} , \quad \boldsymbol{\mathsf{T}}_\tau = \begin{bmatrix} 1 & 0 & 0 \\ 0 & 1 & 0 \\ \tau & 0 & 1 \end{bmatrix} . \end{equation}

In the expression for $\boldsymbol{\mathsf{D}}_\tau (\boldsymbol {k})$, the non-dimensional coefficients $\zeta _i = \zeta _i(\boldsymbol {k},\tau )$, $i=1,2,3$, are defined

(4.14a–c)\begin{equation} \zeta_1 = C_1 - C_2 k_2/k_1, \quad \zeta_2 = C_1 k_2/k_1 + C_2, \quad \zeta_3 = k_0^2/k^2, \end{equation}

where $k_0 = |\boldsymbol {k}_0|$ and

(4.15a,b)\begin{align} C_1 = \frac{\tau k_1^2(k_0^2 -2k_{30}^2 + \tau k_1 k_{30})}{k^2(k_1^2+k_2^2)}, \quad C_2 = \frac{k_2 k_0^2}{(k_1^2+k_2^2)^{3/2}} \arctan\left( \frac{\tau k_1(k_1^2 + k_2^2)^{1/2}}{k_0^2 - \tau k_{30}k_1} \right) . \end{align}

One may observe that

(4.16)\begin{equation} \begin{bmatrix} 1 & 0 & \zeta_1\\ 0 & 1 & \zeta_2\\ 0 & 0 & \zeta_3 \end{bmatrix} \begin{bmatrix} 0 & -k_{30} & k_2\\ k_{30} & 0 & -k_1\\ -k_2 & k_1 & 0 \end{bmatrix} = \begin{bmatrix} 0 & -k_{3} & k_2\\ k_3 & 0 & -k_1\\ -k_2 & k_1 & 0 \end{bmatrix} \begin{bmatrix} \zeta_3 & 0 & 0\\ 0 & \zeta_3 & 0\\ -\zeta_1 & -\zeta_2 & 1 \end{bmatrix} , \end{equation}

or, equivalently, $\boldsymbol{\mathsf{D}}_\tau (\boldsymbol {k}) k_0^{-2}\boldsymbol{\mathsf{Q}}(\boldsymbol {k}_0) = k^{-2}\boldsymbol{\mathsf{Q}}(\boldsymbol {k}) \boldsymbol{\mathsf{D}}_\tau ^{-\top }(\boldsymbol {k})$. Moreover, $\mathrm {d}\boldsymbol {W}(\boldsymbol {k}_0) = \mathrm {d}\boldsymbol {W}(\boldsymbol {k})$, owing to translational invariance. Therefore, taking $\mathrm {d} \boldsymbol {Z}(\boldsymbol {k}_0,0) = \boldsymbol{\mathsf{Q}}(\boldsymbol {k}_0) (({1}/{\sqrt {4{\rm \pi} }k_0^2}) E^{1/2} (\boldsymbol {k}_0) \,\mathrm {d}\boldsymbol {W}(\boldsymbol {k}_0)) ,$ it holds that

(4.17)\begin{equation} \mathrm{d} \boldsymbol{Z}(\boldsymbol{k},t) = \boldsymbol{\mathsf{Q}}(\boldsymbol{k}) \left(\frac{1}{\sqrt{4{\rm \pi}}k^2} E^{1/2}(\boldsymbol{\mathsf{T}}_\tau\boldsymbol{k})\, \boldsymbol{\mathsf{D}}_\tau^{-\top}(\boldsymbol{k})\,\mathrm{d}\boldsymbol{W}(\boldsymbol{k})\right) . \end{equation}

Finally, invoking the general expression for $E(\boldsymbol {k})$ written in (4.1), one arrives at the rapid distortion equation FPDE

(4.18)\begin{equation} \left(I-\boldsymbol{\nabla}\boldsymbol{\cdot}(\bar{\boldsymbol{\varTheta}}_\tau\boldsymbol{\nabla})\right)^{\alpha_1}\left(-\boldsymbol{\nabla}\boldsymbol{\cdot}(\bar{\boldsymbol{\varTheta}}_\tau\boldsymbol{\nabla})\right)^{\alpha_2}\boldsymbol{\psi} = \mu\det(\bar{\boldsymbol{\varTheta}}_\tau)^{\gamma/3}{\boldsymbol{\eta}}_\tau \end{equation}

where $\bar {\boldsymbol {\varTheta }}_\tau = \boldsymbol{\mathsf{T}}_\tau ^\top \bar {\boldsymbol {\varTheta }} \boldsymbol{\mathsf{T}}_\tau$ and ${\boldsymbol {\eta }}_\tau (\boldsymbol {x}) = \int _{\mathbb {R}^3} \operatorname {e}^{\operatorname {i}\boldsymbol {k}\boldsymbol {\cdot }\boldsymbol {x}} (\boldsymbol{\mathsf{D}}_\tau ^{-\top }(\boldsymbol {k})\,\mathrm {d}\boldsymbol {W}(\boldsymbol {k}) )$. Note that $\det (\bar {\boldsymbol {\varTheta }}_\tau ) = \det (\bar {\boldsymbol {\varTheta }})$.

Remark 4.3 An important extension of the rapid distortion model above involves replacing the constant $\tau$ by a wavenumber-dependent ‘eddy lifetime’ $\tau (k)$ (see, e.g., Mann Reference Mann1994). Such models are considered more realistic because, at some point, the shear from the mean velocity gradient will cause the eddies to stretch and eventually they will breakup within a size-dependent timescale. In this case, the generalization of ${\boldsymbol {\eta }}_\tau$ above is straightforward. Meanwhile, at least when $\bar {\boldsymbol {\varTheta }} = L^2 \boldsymbol{\mathsf{I}}$, one may consider replacing the operator $\bar {\boldsymbol {\varTheta }}_\tau$ in (4.18) by

(4.19)\begin{equation} L^2\, \mathcal{F}^{-1} \begin{bmatrix} 1+\tau(k)^2 & 0 & \tau(k) \\ 0 & 1 & 0 \\ \tau(k) & 0 & 1 \end{bmatrix} \mathcal{F}. \end{equation}

To solve such an equation numerically, one does not need to construct the closed form of the linear operator, but may instead choose to use a matrix-free Krylov method (Saad Reference Saad2003).

4.3. Boundary conditions

There are a number of different, equivalent, definitions of fractional operators on $\mathbb {R}^3$. However, moving from the free-space equation (4.4) to a boundary value problem relies on heuristics and can be done in a wide variety of ways; each of which may also differ by the specific definition of the fractional operator being used (Lischke et al. Reference Lischke2020). As stated previously, in this work, we choose to only deal with the spectral definition. In this setting, boundary conditions are applied to the corresponding integer-order operator and then incorporated implicitly by modifying the spectrum; cf. (3.6) and (3.7).

Assume that (4.4) is posed on a three-dimensional simply connected domain $\varOmega \subsetneq \mathbb {R}^3$ with boundary $\varGamma = \partial \varOmega$. We begin with the following heuristically chosen impermeability condition for the velocity field:

(4.20)\begin{equation} \boldsymbol{u} = \operatorname{\boldsymbol{\nabla}\times} \boldsymbol{\psi} \quad \text{in } \varOmega , \quad \boldsymbol{u}\boldsymbol{\cdot}\boldsymbol{n} = 0 \quad \text{on } \varGamma . \end{equation}

Although more relaxed boundary conditions are also possible, we choose to enforce (4.20) via a no-slip condition on the vector potential $\boldsymbol {\psi }$; specifically,

(4.21)\begin{equation} \boldsymbol{\psi} - (\boldsymbol{\psi}\boldsymbol{\cdot}\boldsymbol{n})\boldsymbol{n} = \boldsymbol{0} \quad \text{on } \varGamma . \end{equation}

It turns out that (4.21) is not enough to uniquely define $\boldsymbol {\psi }$ on $\varOmega$. In fact, the remaining boundary condition must restrict $\boldsymbol {\psi }$ normal to $\varGamma$.

We are somewhat free to select what the remaining boundary condition will be. Both the Dirichlet-type boundary condition $\boldsymbol {\psi }\boldsymbol {\cdot }\boldsymbol {n} = 0$ and the Neumann-type boundary condition $(\boldsymbol {\varTheta }(\boldsymbol {x})\boldsymbol {\nabla }\boldsymbol {\psi })\boldsymbol {n}\boldsymbol {\cdot } \boldsymbol {n} = 0$ are possible candidates which would close the equations. Another option is to enforce a weighted average of those two boundary conditions. To be more specific, we may also consider the generalized (homogeneous) Robin boundary condition

(4.22)\begin{equation} \kappa\boldsymbol{\psi}\boldsymbol{\cdot}\boldsymbol{n} + (\boldsymbol{\varTheta}(\boldsymbol{x})\boldsymbol{\nabla}\boldsymbol{\psi})\boldsymbol{n}\boldsymbol{\cdot} \boldsymbol{n} = 0 \quad \text{on } \varGamma , \end{equation}

where the new model parameter $\kappa \geq 0$ could be inferred from available data.

In this work, we choose to close the equations with (4.22) because it is flexible enough to fit a wide variety of data and simple to implement alongside (4.21). We note that $\kappa$ affects the horizontal velocity near the surface because of its control over the normal component of $\boldsymbol {\psi }$. Thus, in the proposed model, $\kappa$ may be a parameter which distinguishes between different types of surfaces. We also note that, in the limit $\kappa \to \infty$, we recover the boundary condition $\boldsymbol {\psi }\boldsymbol {\cdot }\boldsymbol {n} = 0$. Together with (4.21), it implies the complete Dirichlet boundary condition, $\boldsymbol {\psi } = \boldsymbol {0}$ on $\varGamma$. Hereon, we use the notation $\kappa =\infty$ to indicate this special limiting scenario.

Note that (4.4) can be written $\mathcal {L}\boldsymbol {\psi } = \boldsymbol {b}$, where

(4.23a,b)\begin{align} \mathcal{L} := \left(I-\boldsymbol{\nabla}\boldsymbol{\cdot}(\boldsymbol{\varTheta}(\boldsymbol{x})\boldsymbol{\nabla})\right)^{\alpha_1}\left(-\boldsymbol{\nabla}\boldsymbol{\cdot}(\boldsymbol{\varTheta}(\boldsymbol{x})\boldsymbol{\nabla})\right)^{\alpha_2} \quad\text{and}\quad \boldsymbol{b} := \mu\det(\boldsymbol{\varTheta}(\boldsymbol{x}))^{\gamma/3}{\boldsymbol{\eta}} . \end{align}

In order to define the domain $\mathcal {D}(\mathcal {L})$ of the multi-fractional operator $\mathcal {L}\colon \mathcal {D}(\mathcal {L}) \subseteq [L^2(\varOmega )]^3 \to [L^2(\varOmega )]^3$, we start by letting $A := (I\,{-}\,\boldsymbol {\nabla }\boldsymbol {\cdot }(\boldsymbol {\varTheta }(\boldsymbol {x})\boldsymbol {\nabla })) \colon \mathcal {D}(A) \,{\subseteq}\, [L^2(\varOmega )]^3 \to [L^2(\varOmega )]^3$. For notational convenience, we assume that $A$ has a discrete spectrum.

In the spectral definition of $A^{\alpha _1}$, the domain $\mathcal {D}(A)$ characterizes the boundary conditions on $\varGamma$. In this work, assuming that $\det (\boldsymbol {\varTheta }(\boldsymbol {x}))$ is uniformly bounded from above and below by positive constants, we define

(4.24)\begin{equation} \mathcal{D}(A) = \{ \boldsymbol{\psi} \in [H^2(\varOmega)]^3 \colon \eqref{eqn48}\hbox{ and }\eqref{eqn48}\hbox{ hold in the sense of traces} \} . \end{equation}

For this operator domain, there exists an orthonormal basis of eigenvectors $\{\boldsymbol {a}_j\}_{j=1}^\infty \subseteq \mathcal {D}(A)$, with corresponding eigenvalues $\{a_j\}_{j=1}^\infty$ in non-increasing order (cf. Bolin, Kirchner & Kovács Reference Bolin, Kirchner and Kovács2020). Then, following (3.7), the fractional differential operator $A^{\alpha _1} \colon \mathcal {D}(A^{\alpha _1}) \subseteq [L^2(\varOmega )]^3 \to [L^2(\varOmega )]^3$ is defined

(4.25)\begin{equation} A^{\alpha_1}\boldsymbol{\psi} = \sum_{j=1}^\infty a_j^{\alpha_1} (\boldsymbol{\psi},\boldsymbol{a}_j)_\varOmega\, \boldsymbol{a}_j \end{equation}

and $\mathcal {D}(A^{\alpha _1}) = \{ \boldsymbol {\psi } \in [L^2(\varOmega )]^3\, \colon \sum _{j=1}^\infty a_j^{2\alpha _1} (\boldsymbol {\psi },\boldsymbol {a}_j)_\varOmega ^2 < \infty \}$.

Now consider $A-I \colon \mathcal {D}(A) \to [L^2(\varOmega )]^3$ and note that $\mathcal {L} = A^{\alpha _1} (A-I)^{\alpha _2}$. In this case, $A^{\alpha _1}$ and $(A-I)^{\alpha _2}$ commute because they share the same eigenmodes:

(4.26)\begin{equation} A^{\alpha_1} (A-I)^{\alpha_2} \boldsymbol{\psi} = \sum_{j=1}^\infty a_j^{\alpha_1} (a_j-1)^{\alpha_2} (\boldsymbol{\psi},\boldsymbol{a}_j)_\varOmega\, \boldsymbol{a}_j = (A-I)^{\alpha_2} A^{\alpha_1}\boldsymbol{\psi} . \end{equation}

Accordingly, we define the domain of the operator $\mathcal {L}$ as follows:

(4.27)\begin{equation} \mathcal{D}(\mathcal{L}) = \left\{ \boldsymbol{\psi} \in [L^2(\varOmega)]^3\, \colon \sum_{j=1}^\infty a_j^{2\alpha_1}(a_j-1)^{2\alpha_2} (\boldsymbol{\psi},\boldsymbol{a}_j)_\varOmega^2 < \infty\right\} . \end{equation}

We may now write the boundary value problem given by (4.4), (4.21) and (4.22) as the abstract operator equation $\mathcal {L}\boldsymbol {\psi } = \boldsymbol {b}$, with $\mathcal {D}(\mathcal {L})$ defined in (4.27). Nevertheless, we will still usually refer to this problem in the ‘strong form’

(4.28)\begin{equation} \left. \begin{array}{cl} \left(I-\boldsymbol{\nabla}\boldsymbol{\cdot}(\boldsymbol{\varTheta}(\boldsymbol{x})\boldsymbol{\nabla})\right)^{\alpha_1}\left(-\boldsymbol{\nabla}\boldsymbol{\cdot}(\boldsymbol{\varTheta}(\boldsymbol{x})\boldsymbol{\nabla})\right)^{\alpha_2}\boldsymbol{\psi} = \mu\det(\boldsymbol{\varTheta}(\boldsymbol{x}))^{\gamma/3}{\boldsymbol{\eta}}, & \text{in}\ \varOmega,\\ \boldsymbol{\psi} - (\boldsymbol{\psi}\boldsymbol{\cdot}\boldsymbol{n})\boldsymbol{n} = 0, & \text{on}\ \varGamma,\\ \kappa\boldsymbol{\psi}\boldsymbol{\cdot}\boldsymbol{n} + (\boldsymbol{\varTheta}(\boldsymbol{x})\boldsymbol{\nabla}\boldsymbol{\psi})\boldsymbol{n}\boldsymbol{\cdot}\boldsymbol{n} = 0, & \text{on}\ \varGamma, \end{array} \right\} \end{equation}

because it is much more physically illustrative.

5.1. Shear-free boundary layers

There are many different examples of turbulence confined by a solid boundary, without any significant mean shear (Hunt Reference Hunt1984). In such flows, the rate of turbulent kinetic energy dissipation $\epsilon$ can be assumed to be approximately constant with height. This setting has been studied in detail by various authors (see, e.g., Hunt Reference Hunt1984; Hunt et al. Reference Hunt, Moin, Lee, Moser, Spalart, Mansour, Kaimal and Gaynor1989; Perot & Moin Reference Perot and Moin1995a,Reference Perot and Moinb; Aronson, Johansson & Löfdahl Reference Aronson, Johansson and Löfdahl1997 and references therein) and so forms a solid proving ground to validate (4.28).

5.1.1. A von Kármán-type model

We begin with the inhomogeneous turbulence model (4.28), with fractional coefficients corresponding to the von Kármán energy spectrum (3.11), on the open half space domain $\mathbb {R}_+^3 = \{ (x,y,z) \in \mathbb {R}^3 \colon z >0 \}$. Based on the supposed absence of shear, we also consider the following simple diagonal form for the diffusion tensor, in Cartesian coordinates:

(5.1)\begin{equation} \boldsymbol{\varTheta}(z) = \begin{bmatrix} L_{1}(z)^2 & 0 & 0 \\ 0 & L_{2}(z)^2 & 0 \\ 0 & 0 & L_{3}(z)^2 \end{bmatrix} . \end{equation}

Defining $L(z) = \sqrt [3]{L_1(z)L_2(z)L_3(z)}$, the appropriate form of (4.28) can be written as follows:

(5.2)\begin{equation} \left. \begin{array}{cl} \left(I-\operatorname{\boldsymbol{\nabla}\hspace{1.25pt}\boldsymbol{\cdot}}(\boldsymbol{\varTheta}(z)\boldsymbol{\nabla})\right)^{17/12} \boldsymbol{\psi}= \mu L(z)^{17/6} \boldsymbol{\xi}, & \text{in}\ \mathbb{R}_+^3, \\ \kappa\psi_3 + L_{3}(z)^2 \dfrac{\partial\psi_3}{\partial z} = \psi_1 = \psi_2 = 0, & \text{at}\ z=0 . \end{array} \right\} \end{equation}

Both the Robin coefficient $\kappa$ and an explicit parametric expression for each $L_i(z)$ give rise to a model design parameter vector, say ${\boldsymbol {\theta }}$. This vector ${\boldsymbol {\theta }}$ may then be subject to calibration with respect to experimental data, e.g. using the technique described in § 6.2. This process of model calibration is important because wall roughness, Reynolds number and the nature of the turbulence may affect the near-wall statistics (Pope Reference Pope2001) and may be incorporated through proper parameter selection. For instance, let us consider the following exponential expansion

(5.3)\begin{equation} L_i(z) = L_\infty \cdot\left(1 + \sum_{k=1}^K c_{i,k}\operatorname{e}^{{-}d_{i,k} ({z}/{L_\infty})} \right) , \end{equation}

with each $d_{i,k}\ge 0$, $c_{1,k} = c_{2,k}$ and $d_{1,k} = d_{2,k}$. Taking only two terms in each expansion above ($K=2$), we arrive through calibration at a statistical model which closely matches the experimental data found in Thomas & Hancock (Reference Thomas and Hancock1977). Note that with such a model, $L_1(z) = L_2(z)$ and each $L_i(z)$ exponentially converges to the homogeneous length scale $L_\infty$, as $z\to \infty$, as illustrated in figure 3.

Figure 3. Optimal diffusion coefficients $L_i(z)$ and Robin constant $\kappa$ determined by fitting the Reynolds stress data in figure 4. Note that $L_1(z) = L_2(z)$.

The prescribed boundary conditions will affect the physical length scales of the random velocity field $\boldsymbol {u} = \boldsymbol {\nabla } \times \boldsymbol {\psi }$. Therefore, the diffusion coefficients $L_i(z)$ do not necessarily correspond to the physical length scales. For this reason, we follow Lee & Hunt (Reference Lee and Hunt1991) and define the (physical) so-called integral length scales

(5.4)\begin{equation} \ell_{ij}^{(x_m)}(z) = \frac{\displaystyle\int_{\mathbb{R}}\langle u_i(\boldsymbol{x} + r\boldsymbol{e}_m) u_j(\boldsymbol{x})\rangle\,\mathrm{d} r}{\langle u_i(\boldsymbol{x}) u_j(\boldsymbol{x})\rangle} = \frac{\displaystyle\int_{\mathbb{R}} R_{ij}(r\boldsymbol{e}_m,z) \,\mathrm{d} r}{R_{ij}(\boldsymbol{0},z)} . \end{equation}

In these expressions, we have accounted for the fact that all solutions of (5.2) are temporary stationary and statistically homogeneous in the $x$- and $y$-directions, i.e. $R(\boldsymbol {r},\boldsymbol {x},t) = R(\boldsymbol {r},z)$.

In § 6, we explain how to solve this problem numerically and to calibrate its solutions to Reynolds stress data. The difference between the Reynolds stress profiles in the calibrated model and the corresponding experimental data is depicted in figure 4, alongside the resulting integral length scales $\ell _{ij}^{(x_m)}(z)$. Because this model has many free parameters which can be calibrated to experimental data, it is much more flexible than the classical theory proposed by Hunt (Reference Hunt1973, Reference Hunt1984) and Hunt & Graham (Reference Hunt and Graham1978). Indeed, a comparison between the two theories, which highlights this flexibility, is given in the next subsection. Note that the exact definitions of the optimized model parameters used in the results above are stated explicitly in the table in figure 3.

Figure 4. Reynolds stress data from Thomas & Hancock (Reference Thomas and Hancock1977) compared with Reynolds stresses from (a) the calibrated SFBL turbulence model (5.2) and (b) corresponding integral length scales. Observe that the model is able to closely fit the experimental data.

5.1.2. Comparison with the classical theory

It is important to consider the special case of (5.2) where each $L_i(z)$ is constant in $z$. In Hunt's idealized shear-free boundary layer (SFBL) theory (Hunt & Graham Reference Hunt and Graham1978; Hunt Reference Hunt1984), derived from the energy spectrum ansatz (3.11) and briefly summarized in § 2, one can show that

(5.5)\begin{equation} \frac{\langle u^2 \rangle}{\langle u^2_\infty \rangle} = \frac{\langle v^2 \rangle}{\langle v^2_\infty \rangle} \to 1.5 \quad\text{and} \quad \frac{\langle w^2 \rangle}{\langle w^2_\infty \rangle} = O\left(\left(\frac{z}{L_\infty}\right)^{2/3}\right) \quad\text{as} \ \frac{z}{L_\infty} \to 0, \end{equation}

where $\langle u^2_\infty \rangle = \langle v^2_\infty \rangle = \langle w^2_\infty \rangle$ denotes the far-field limit $z\to \infty$ of the non-zero Reynolds stresses. The limit $\langle u^2 \rangle /\langle u^2_\infty \rangle \to 1.5$ is not always achieved in experiments (cf. figure 4), however, the limiting behaviour $\langle w^2 \rangle /\langle w^2_\infty \rangle = O((z/L_\infty )^{2/3})$ is well-established in the literature (Priestley Reference Priestley1959; Kaimal et al. Reference Kaimal, Wyngaard, Haugen, Coté, Izumi, Caughey and Readings1976).

The corresponding scenario in our class of models is exactly (5.2) with each $L_i = L_\infty$. In this setting, the non-zero Reynolds stresses, $\langle u^2 \rangle = \langle v^2 \rangle$ and $\langle w^2 \rangle$, can be derived analytically, at least for certain values of $\kappa \geq 0$. These exact analytical solutions are summarized in Lemmas 5.1–5.3, the proofs of which can be found in Appendix A. Exact analytical solutions for the integral length scales $\ell _{ij}^{(x_m)}(z)$ can also be derived by a similar technique, but we do not include their derivation in this work for the sake of brevity. Plots of the analytical Reynolds stresses and integral length scales are depicted in figure 5.

Lemma 5.1 Given $\boldsymbol {u} = (u,v,w) = \boldsymbol {\nabla }\times \boldsymbol {\psi }$, where $\boldsymbol {\psi }$ is any solution of (5.2) with constant $L_1 = L_2 = L_\infty$, it holds that

(5.6)\begin{equation} \frac{\langle w^2 \rangle}{\langle w^2_\infty \rangle} = 1 - \mathcal{M}_{1/3}\left(\frac{2z}{L_\infty}\right) , \end{equation}

where $\mathcal {M}_{\nu }(x)$ is the Matérn kernel (Matérn Reference Matérn1986; Stein Reference Stein1999; Khristenko et al. Reference Khristenko, Scarabosio, Swierczynski, Ullmann and Wohlmuth2019) given by

(5.7)\begin{equation} \mathcal{M}_{\nu}(x) = \frac{x^\nu\mathrm{K}_{\nu}(x)}{2^{\nu-1}\varGamma(\nu)}, \quad \nu \geq 0, \end{equation}

and $\mathrm {K}_{\nu }(x)$ denotes the modified Bessel function of the second kind (Abramowitz & Stegun Reference Abramowitz and Stegun1948; Bateman Reference Bateman1953; Watson Reference Watson1995). Moreover, near the boundary the following asymptotic holds:

(5.8)\begin{equation} \frac{\langle w^2 \rangle}{\langle w^2_\infty \rangle} \sim \frac{\varGamma(2/3)}{\varGamma(4/3)} \left(\frac{z}{L_\infty}\right)^{2/3} \quad\text{as} \ \frac{z}{L_\infty} \to 0 . \end{equation}

Figure 5. (a) The analytically derived non-zero Reynolds stresses stated in Lemmas 5.1–5.3 and (b) corresponding integral length scales.

Lemma 5.2 Given $\boldsymbol {u} = (u,v,w) = \boldsymbol {\nabla }\times \boldsymbol {\psi }$, where $\boldsymbol {\psi }$ is the solution of (5.2) with constant $\boldsymbol {\varTheta } = L_\infty ^2\boldsymbol{\mathsf{I}}$ and $\kappa =0$, it holds that

(5.9)\begin{equation} \frac{\langle u^2 \rangle}{\langle u^2_\infty \rangle} = \frac{\langle v^2 \rangle}{\langle v^2_\infty \rangle} = 1 + (\nu+1)\mathcal{M}_{\nu}\left(\frac{2z}{L_\infty}\right) - \nu\mathcal{M}_{\nu+1}\left(\frac{2z}{L_\infty}\right) . \end{equation}

Hence, near the boundary,

(5.10)\begin{equation} \frac{\langle u^2 \rangle}{\langle u^2_\infty \rangle} = \frac{\langle v^2 \rangle}{\langle v^2_\infty \rangle} \to 2 \quad \textrm{as}\ \frac{z}{L_\infty} \to 0. \end{equation}

Lemma 5.3 Given $\boldsymbol {u} = (u,v,w) = \boldsymbol {\nabla }\times \boldsymbol {\psi }$, where $\boldsymbol {\psi }$ is any solution of (5.2) with constant $\boldsymbol {\varTheta } = L_\infty ^2\boldsymbol{\mathsf{I}}$ and $\kappa =\infty$, it holds that

(5.11)\begin{equation} \frac{\langle u^2 \rangle}{\langle u^2_\infty \rangle} = \frac{\langle v^2 \rangle}{\langle v^2_\infty \rangle} = 1 + \nu\mathcal{M}_{\nu}\left(\frac{2z}{L_\infty}\right) - \nu\mathcal{M}_{\nu+1}\left(\frac{2z}{L_\infty}\right). \end{equation}

Hence, near the boundary,

(5.12)\begin{equation} \frac{\langle u^2 \rangle}{\langle u^2_\infty \rangle} = \frac{\langle v^2 \rangle}{\langle v^2_\infty \rangle} \to 1 \quad \textrm{as}\ \frac{z}{L_\infty} \to 0. \end{equation}

Remark 5.4 The Robin boundary condition $\kappa \psi _3 + L_\infty ^2 ({\partial \psi _3}/{\partial z}) = 0$ has no effect on $\langle w^2 \rangle$. Therefore, the asymptotic expansion of the well-known (Priestley Reference Priestley1959; Kaimal et al. Reference Kaimal, Wyngaard, Haugen, Coté, Izumi, Caughey and Readings1976; Hunt Reference Hunt1984; Hunt et al. Reference Hunt, Moin, Lee, Moser, Spalart, Mansour, Kaimal and Gaynor1989) asymptotic behaviour $\langle w^2 \rangle /\langle w^2_\infty \rangle = O((z/L_\infty )^{2/3})$ as $z/L_\infty \to 0$ always holds when $L_1 = L_2 = L_\infty$.

Remark 5.5 The limit $\langle u^2 \rangle /\langle u^2_\infty \rangle \to 1.5$ from Hunt's theory lies exactly in between the range of analogous limits, $\langle u^2 \rangle /\langle u^2_\infty \rangle \to 1$ and $\langle u^2 \rangle /\langle u^2_\infty \rangle \to 2$, coming from the exact solutions of (5.2) when $\kappa = \infty$ and $\kappa = 0$, respectively. Numerical experiments show that $\langle u^2 \rangle /\langle u^2_\infty \rangle = \langle v^2 \rangle /\langle v^2_\infty \rangle$ always limits to a value in the interval $(1,2)$ when $\kappa \in (0,\infty )$ and $\boldsymbol {\varTheta } = L_\infty ^2\boldsymbol{\mathsf{I}}$.

5.1.3. A more general energy spectrum

In order to illustrate the dependence of (4.28) on the parameter $\alpha _2 = 17/12 -\alpha _1$, we may consider an alternative form of (5.2) which corresponds to the energy spectrum (4.2) with $p_0 = 2$. Here, for additional complexity, we also consider the load ${\boldsymbol {\eta }} = \boldsymbol {\xi }_\beta$ with $\beta /L_\infty =10^{-2}$:

(5.13)\begin{equation} \left. \begin{array}{cl} \left(I-\operatorname{\boldsymbol{\nabla}\hspace{1.25pt}\boldsymbol{\cdot}}(\boldsymbol{\varTheta}(z)\boldsymbol{\nabla})\right)^{11/12} \left(-\operatorname{\boldsymbol{\nabla}\hspace{1.25pt}\boldsymbol{\cdot}}(\boldsymbol{\varTheta}(z)\boldsymbol{\nabla})\right)^{1/2} \boldsymbol{\psi} = \mu L(z)^{17/6} \boldsymbol{\xi}_\beta, & \text{in}\ \mathbb{R}_+^3,\\ \kappa\psi_3 + L_3(z)^2 \dfrac{\partial\psi_3}{\partial z} = \psi_1 = \psi_2 = 0 , & \text{at}\ z=0 . \end{array} \right\} \end{equation}

We do not analyze these equations in detail here, however, we present a single realization of their solution figure 6 for visual comparison. Observe that the velocity field coming from (5.13) is visibly smoother than its counterpart coming from (5.2). This is due to the high regularity load $\boldsymbol {\xi }_\beta$.

Figure 6. Normalized magnitudes of $\boldsymbol {u} = \operatorname {\boldsymbol {\nabla }\times }\boldsymbol {\psi }$ from (a) (5.2), (b) (5.13) with $\beta /L_\infty = 10^{-2}$ and (c) (5.15) with $\tau = 1.0$. The additional model parameters are specified in figure 3. Observe that the central field is visibly smoother than its counterpart on the left owing to the high regularity load $\boldsymbol {\xi }_\beta$. The field on the right, issued from the same noise, presents distortion.

5.2. Uniform shear boundary layers

Classically, rapid distortion theory is used to describe the short time evolution of isotropic turbulence. As pointed out in, e.g., Lee & Hunt (Reference Lee and Hunt1991), it is also possible to extend its use to some examples of inhomogeneous turbulence. In this example, we follow Lee & Hunt (Reference Lee and Hunt1991) in considering a uniform shear boundary layer (USBL) model where the only effect of the wall is to block velocity fluctuations in the normal direction. Our derivation begins from the assumption $\langle \boldsymbol {U}(\boldsymbol {x})\rangle = (U_0 + S x_3) \boldsymbol {e}_1$ taken in § 4.2, but we also allow for a $z$-dependent inhomogeneous diffusion tensor,

(5.14)\begin{equation} \boldsymbol{\varTheta}_\tau(z) = \begin{bmatrix} L_1(z)^2+\tau^2 L_3(z)^2 & 0 & \tau L_3(z)^2\\ 0 & L_2(z)^2 & 0 \\ \tau L_3(z)^2 & 0 & L_3(z)^2 \end{bmatrix} . \end{equation}

With this expression in hand, we may consider the following inhomogeneous version of (4.4) with $\tau = 1.0$:

(5.15) \begin{equation} \left. \begin{array}{cl} \left(I-\boldsymbol{\nabla}\boldsymbol{\cdot}(\boldsymbol{\varTheta}_\tau(z)\boldsymbol{\nabla})\right)^{17/12}\boldsymbol{\psi} = \mu L(z)^{17/6}{\boldsymbol{\eta}}_\tau, & \text{in}\ \mathbb{R}_+^3, \\ \kappa\psi_3 + L_3(z)^2\left(\dfrac{\partial\psi_3}{\partial z} + \tau\dfrac{\partial\psi_3}{\partial x}\right)= \psi_1 = \psi_2 = 0, & \text{at}\ z=0 . \end{array} \right\} \end{equation}

It is possible that the inhomogeneous length scales in this tensor, $L_i(z)$, may be tuned to compensate for the presence of small non-zero Reynolds stress gradients, however, we do not seek to verify that hypothesis here. Instead, we settle for a visual comparison between the solutions of the various models.

Figure 6 depicts a reference velocity field coming from a single realization of (5.2), (5.13) and (5.15). In order to demonstrate the flexibility of the models, we have taken the same calibrated model parameters used in § 5.1.1. For a fair reference, we have also used the same additive white Gaussian noise vector to generate the load for each realization.

5.3. Turbulent inlet generation for numerical wind tunnel simulations

The mean profile $\langle \boldsymbol {U}(z) \rangle$ in many wall-bounded shear flows is often assumed to follow a logarithmic curve, sometimes with a Reynolds number modification (see, e.g., Barenblatt & Chorin Reference Barenblatt and Chorin2004). In the atmospheric boundary layer, one such model for the mean velocity, $\langle \boldsymbol {U}(\boldsymbol {x})\rangle = U(z)\boldsymbol {e}_1$, found in the wind engineering community is written in terms of the height above ground, $z$, as follows (Mendis et al. Reference Mendis, Ngo, Haritos, Hira, Samali and Cheung2007; Kareem & Tamura Reference Kareem and Tamura2013):

(5.16)\begin{equation} U(z) = \frac{u_\ast}{\kappa} \ln\left(\frac{z-d}{z_0}\right). \end{equation}

Here, $u_\ast$ is the friction velocity, $z_0$ is the roughness length and $d$ is the zero-plane displacement. Although all such models violate the uniform shear assumption made in deriving (4.18), it has been argued that the assumption is still valid for describing eddies of ‘linear dimension smaller than the length over which the shear changes appreciably’ (Mann Reference Mann1994, p. 145). For this reason, turbulence models similar to those presented in the previous subsections (see, e.g. Mann Reference Mann1994, Reference Mann1998; Chougule et al. Reference Chougule, Mann, Kelly and Larsen2018), have established themselves in wind engineering (IEC 61400-1:2005). An account of some physical violations of such models is given in detail in Hunt (Reference Hunt1984) and Hunt et al. (Reference Hunt, Moin, Lee, Moser, Spalart, Mansour, Kaimal and Gaynor1989). It remains to be demonstrated whether the non-homogenous diffusion coefficient in, e.g., (5.15) may ameliorate some of these issues.

Our final application involves using (5.15) to generate synthetic turbulent inlet conditions, which is an important application in computational fluid dynamics as a whole (Tabor & Baba-Ahmadi Reference Tabor and Baba-Ahmadi2010). We choose to follow an established approach used in the wind engineering industry (see Michalski et al. Reference Michalski, Kermel, Haug, Löhner, Wüchner and Bletzinger2011; Andre, Mier-Torrecilla & Wüchner Reference Andre, Mier-Torrecilla and Wüchner2015 and references therein). Here, a contiguous section of spatially correlated turbulence is transformed into a stationary Gaussian process by identifying the $x$-component of the turbulent velocity field with a time axis via the transformation $x = U_{m} t$. Then, at each time step $t= t_k$, the turbulent fluctuations $\boldsymbol {U}(\boldsymbol {x})|_{x = t_k/U_{m}}$ are projected onto the inflow boundary of a numerical wind tunnel; see the depiction in figure 7. Here, $U_{m} > 0$ is a mean velocity parameter which directly affects the spatial-to-temporal correlation of the synthetic turbulent inlet boundary conditions. With this application, we highlight the potential of calibrated FPDE models to improve the accuracy of numerical wind tunnel simulations.

Remark 5.6 The physical justification for the transformation $x = U_{m} t$ derives from a manipulated Taylor's hypothesis, as described in Mann (Reference Mann1994, § 2.3).

Figure 7. Snapshots of synthetic wind, $\boldsymbol {U}(\boldsymbol {x}) = \langle \boldsymbol {U}(\boldsymbol {x})\rangle + \boldsymbol {u}(\boldsymbol {x})$, mapped onto the inlet boundary in a numerical wind tunnel test of a modern high rise building. Turbulent fluctuations $\boldsymbol {u}(\boldsymbol {x})$ generated using model (5.15) with $\tau = 1.0$, $\kappa = 0$, and $L_1(z) = L_2(z) = L_3(z) = L_\infty$. The LES was performed with the finite element software Kratos Multiphysics (Dadvand, Rossi & Oñate Reference Dadvand, Rossi and Oñate2010).

6.1. Numerical solution of FPDEs

Numerical solution of boundary value problems involving fractional powers of elliptic operators is challenging and computationally expensive, owing in part to the non-locality of the resulting operator. Methods based on direct diagonalization of the elliptic operator (Ilic et al. Reference Ilic, Liu, Turner and Anh2005; Yang et al. Reference Yang, Turner, Liu and Ilić2011) are generally too expensive for practical applications. Common methods for solving fractional differential equations involve introducing an additional dimension, provided the extension operator is local, which defines the increased computational complexity. A typical example is Caffarelli–Silvestre extension (Caffarelli & Silvestre Reference Caffarelli and Silvestre2007; Nochetto, Otárola & Salgado Reference Nochetto, Otárola and Salgado2015), where the FPDE solution is viewed as a restriction onto a subspace of the solution to an associated integer-order boundary value problem on a higher-dimensional semi-infinite cylinder. The resulting higher-dimensional integer-order PDEs can thus be solved using Galerkin approach (see, e.g., Meidner et al. Reference Meidner, Pfefferer, Schuürholz and Vexler2018; Banjai et al. Reference Banjai, Melenk, Nochetto, Otarola, Salgado and Schwab2019; Banjai, Melenk & Schwab Reference Banjai, Melenk and Schwab2020). In Balakrishnan et al. (Reference Balakrishnan1960), Bonito & Pasciak (Reference Bonito and Pasciak2015) and Bonito, Lei & Pasciak (Reference Bonito, Lei and Pasciak2019), the integral representation of the inverse fractional operator is considered, introducing a special quadrature rule, which yields a series of integer-order PDEs. Alternatively, Vabishchevich (Reference Vabishchevich2015) and Lazarov & Vabishchevich (Reference Lazarov and Vabishchevich2017) propose to represent the FPDE as a transient pseudo-parabolic problem. Then, numerical integration in ‘time’ (transition parameter) provides a sequence of integer-order PDEs. In this work, we follow the rational approximation approach (see, e.g., Harizanov & Margenov Reference Harizanov and Margenov2018; Bolin & Kirchner Reference Bolin and Kirchner2019). It consists of approximating of the operator's spectrum with a rational function. Multipole expansion of this function defines a family of independent integer-order PDEs, which can be solved in parallel. The method is briefly summarized in the following. The interested reader is referred to Bonito et al. (Reference Bonito, Borthagaray, Nochetto, Otárola and Salgado2018) and Lischke et al. (Reference Lischke2020) and the references therein for further information on fractional diffusion problems.

Let $A$ be an abstract bounded elliptic symmetric positive definite operator with spectrum $\sigma (A)\subseteq [\lambda _{min}, \lambda _{max}]$, $0<\lambda _{min}< \lambda _{max}$. For illustration, consider the associated fractional problem

(6.1)\begin{equation} A^\alpha\boldsymbol{\psi} = \boldsymbol{b}, \end{equation}

for some $\alpha >0$. If the rational function $r_N(\lambda ) = \sum _{n=1}^{N}({c_n}/{\lambda + d_n})$ approximates the function $f(\lambda )=\lambda ^{-\alpha }$ on the interval $[\lambda _{min}, \lambda _{max}]$, then the solution $\boldsymbol {\psi }$ can be approximated as the weighted average of solutions of $N$ other elliptic problems; namely,

(6.2a,b)\begin{equation} \boldsymbol{\psi} \approx \sum_{n=1}^{N} c_n\boldsymbol{\psi}_n, \qquad \left(d_nI + A\right)\boldsymbol{\psi}_n = \boldsymbol{b} . \end{equation}

If $A$ is an integer-order differential operator, e.g. $A = I-\boldsymbol {\nabla }\boldsymbol {\cdot }(\boldsymbol {\varTheta }(\boldsymbol {x})\boldsymbol {\nabla })$, then each of these $N$ problems can be solved using standard discretization methods for integer-order operators, e.g. finite elements. Remark 6.2 contains a number of general comments about such discretizations. For the reader's interest, an example of the numerical method we used for the problems in § 5 is described in brief in Appendix B.

The rational approximation technique above can be extended to the solution of (4.28), which, notably, has two fractional powers, $\alpha _1$ and $\alpha _2$. Indeed, in this case, we need to construct a rational approximation $r_N(\lambda )$ for the function $f(\lambda ) = \lambda ^{-\alpha _1}(\lambda -1)^{-\alpha _2}$. With this alternative rational approximation in hand, the approximate vector potential $\tilde {\boldsymbol {\psi }}$ is again given by (6.2a,b).

Remark 6.2 Note that the load $\boldsymbol {b} = \mu \det (\boldsymbol {\varTheta }(\boldsymbol {x}))^{\gamma /3}{\boldsymbol {\eta }}$ in (4.28) is a random variable. The reader is referred to Lindgren, Rue & Lindström (Reference Lindgren, Rue and Lindström2011), Du & Zhang (Reference Du and Zhang2002) and Croci et al. (Reference Croci, Giles, Rognes and Farrell2018) for details of numerical solution to stochastic PDEs and approximation of additive white Gaussian noise. Typically, a discretization of the integer-order operator equation $(d_nI + A)\boldsymbol {\psi }_n = \boldsymbol {b}$ results in a linear system

(6.3)\begin{equation} (d_n\boldsymbol{\mathsf{M}} + \boldsymbol{\mathsf{A}})\boldsymbol{\mathsf{p}}_n = \boldsymbol{\mathsf{b}}, \quad \text{with}\ \boldsymbol{\mathsf{b}}\sim \mathcal{N}(0,\boldsymbol{\mathsf{B}}), \end{equation}

where the vector $\boldsymbol{\mathsf{p}}_n$ denotes the coefficients of the discrete solution $\boldsymbol {\psi }^h_n$ in a preselected basis, say $\varPhi$. Here, $\boldsymbol{\mathsf{M}}$ is a discretization of the identity operator $I$, $\boldsymbol{\mathsf{A}}$ is a discretization of the integer order differential operator $A$, and $\boldsymbol{\mathsf{B}} = \langle \boldsymbol{\mathsf{b}}\boldsymbol{\mathsf{b}}^\top \rangle$ is a given covariance matrix. Via a change of variables, the random load may also be written $\boldsymbol{\mathsf{b}} = \boldsymbol{\mathsf{H}}\boldsymbol {\xi }$, where $\boldsymbol{\mathsf{H}}\boldsymbol{\mathsf{H}}^\top = \boldsymbol{\mathsf{B}}$ and $\boldsymbol {\xi } \sim \mathcal {N}(0, {\mathsf{I}})$ is a standard Gaussian vector $\langle \boldsymbol {\xi }\boldsymbol {\xi }^\top \rangle = {\mathsf{I}}$, with ${\mathsf{I}}$ denoting the identity matrix. One particular form of $\boldsymbol{\mathsf{H}}$ comes from the Cholesky decomposition, although many other are factorizations are also possible (Croci et al. Reference Croci, Giles, Rognes and Farrell2018; Kessy, Lewin & Strimmer Reference Kessy, Lewin and Strimmer2018). Finally, note that if the same basis $\varPhi$ is used the solve for each $\boldsymbol {\psi }_n^h$, then the discrete solution $\boldsymbol {\psi }^h = \sum _{n=1}^{N} c_n\boldsymbol {\psi }_n^h \approx \boldsymbol {\psi }$ can also be expressed using $\varPhi$, with the coefficient vector $\boldsymbol{\mathsf{p}} = \sum _{n=1}^N c_n \boldsymbol{\mathsf{p}}_n$.

6.2. Fitting Reynolds stress data

Various statistical quantities of a turbulent flow field can be measured experimentally. Near a solid boundary, some of the most important of these quantities are the Reynolds stresses $\tau _{ij} = \langle u_{i}u_{j} \rangle$. In order to calibrate the parameters in (5.2) to Reynolds stress data $\tau _{ij}^\text {data}(\boldsymbol {x}_l)$, collected at a number of locations in the flow domain $\boldsymbol {x}_l\in S$, we propose the following optimization problem:

(6.4)\begin{equation} \min_{\boldsymbol{\theta}} \mathcal{J}({\boldsymbol{\theta}}), \quad \text{where} \ \mathcal{J}({\boldsymbol{\theta}}) = \frac{1}{|S|} \sum_{\boldsymbol{x}_l\in S} \sum_{i,j=1}^3 \left( \tau_{ij}(\boldsymbol{x}_l;{\boldsymbol{\theta}}) - \tau_{ij}^\text{data} (\boldsymbol{x}_l) \right)^2 . \end{equation}

Here, the design variable ${\boldsymbol {\theta }}$ denotes a coefficient vector taking accounting for all of the undetermined model parameters present in (4.28). For instance, in (5.1.1) we used

(6.5)\begin{equation} {\boldsymbol{\theta}} = (c_{1,1}, d_{1,1}, c_{3,1}, d_{3,1},\ldots,c_{1,K}, d_{1,K}, c_{3,K}, d_{3,K},\kappa) \in \mathbb{R}^{4 K+1} , \end{equation}

where $c_{i,k}$ and $d_{i,k}$, $i=1,3$, $k=1,\ldots ,K$, appear in the representation of each $L_i(z)$ with ${K=2}$ terms; cf. (5.3).

Remark 6.4 It turns out that (6.4) can be rewritten as a deterministic optimization problem. To see this, recall Remark 6.2 and consider the common basis $\varPhi = \{\phi _m\boldsymbol {e}_i \,\colon m = 1,\ldots , M,~ i = 1,2,3\} \subseteq [H^1(\varOmega )]^3$ for the discretization (6.3) of each sub-problem (6.2a,b). We may then write $\boldsymbol{\mathsf{p}} = ({\mathsf{p}}_1, \ldots , {\mathsf{p}}_{3M})\in \mathbb {R}^{3M}$ and $\boldsymbol {\psi }^h = \sum _{i=1}^3\sum _{m=1}^M {\mathsf{p}}_{m + (i-1)\cdot M}\phi _{m}\boldsymbol {e}_i$. Likewise, we may also write $\boldsymbol {u}^h = \sum _{i=1}^3\sum _{m=1}^M {\mathsf{p}}_{m + (i-1)\cdot M}\operatorname {\boldsymbol {\nabla }\times }(\phi _{m}\boldsymbol {e}_i)$. As remarked previously, $\boldsymbol{\mathsf{p}} = \sum _{n=1}^{N}(d_n\boldsymbol{\mathsf{M}} + \boldsymbol{\mathsf{A}})^{-1}c_n\boldsymbol{\mathsf{b}}$, where $\boldsymbol{\mathsf{b}} \sim \mathcal {N}(0,\boldsymbol{\mathsf{B}})$. Note that both matrices $\boldsymbol{\mathsf{A}}$ and $\boldsymbol{\mathsf{B}}$ generally depend on ${\boldsymbol {\theta }}$. Throughout the rest of this section, we use the shorthand $\boldsymbol{\mathsf{L}}^{-1}$ to denote the linear operator $\sum _{n=1}^{N}(d_n\boldsymbol{\mathsf{M}} + \boldsymbol{\mathsf{A}})^{-1}c_n$. With this notation at our disposal, we may simply write $\boldsymbol{\mathsf{p}} = \boldsymbol{\mathsf{L}}^{-1}\boldsymbol{\mathsf{b}}$ or, equivalently, $\boldsymbol{\mathsf{L}}\boldsymbol{\mathsf{p}} = \boldsymbol{\mathsf{b}}$. An associated adjoint problem can be used to approximate $\tau _{ij}$ at any location $\boldsymbol {x}_l$.

Suppose that we wish to evaluate the covariance tensor $\langle u_i(\boldsymbol {x}) u_j(\boldsymbol {y})\rangle$ at a point, say $\boldsymbol {x}_l$. This may be approximated by applying the delta function (or some approximation thereof) in both $\boldsymbol {x}$- and $\boldsymbol {y}$-coordinates to $\langle u_i^h(\boldsymbol {x}) u_{j}^h(\boldsymbol {y})\rangle$:

(6.6)\begin{equation} \langle u_i^h(\boldsymbol{x}_l) u_{j}^h(\boldsymbol{x}_l)\rangle = \int_\varOmega \int_\varOmega \delta(\boldsymbol{x}-\boldsymbol{x}_l) \langle u_i^h(\boldsymbol{x}) u_{j}^h(\boldsymbol{y})\rangle \delta(\boldsymbol{y}-\boldsymbol{x}_l)\, \mathrm{d}\boldsymbol{x}\,\mathrm{d}\boldsymbol{y} . \end{equation}

Upon substitution of the expression $u_j^h = \sum _{i=1}^3\sum _{m=1}^M {\mathsf{p}}_{m + (i-1)\cdot M}\operatorname {\boldsymbol {\nabla }\times }(\phi _{m}\boldsymbol {e}_i)\boldsymbol {\cdot } \boldsymbol {e}_j$, we find that

(6.7)\begin{equation} \langle u_i^h(\boldsymbol{x}_l) u_{j}^h(\boldsymbol{x}_l)\rangle = \boldsymbol{\mathsf{d}}_{i,l}^\top \langle\boldsymbol{\mathsf{p}}\boldsymbol{\mathsf{p}}^\top\rangle \boldsymbol{\mathsf{d}}_{j,l} = \boldsymbol{\mathsf{d}}_{i,l}^\top \boldsymbol{\mathsf{L}}^{{-}1} \langle\boldsymbol{\mathsf{b}}\boldsymbol{\mathsf{b}}^\top\rangle \boldsymbol{\mathsf{L}}^{{-}1} \boldsymbol{\mathsf{d}}_{j,l} = \boldsymbol{\mathsf{d}}_{i,l}^\top \boldsymbol{\mathsf{L}}^{{-}1} \boldsymbol{\mathsf{B}} \boldsymbol{\mathsf{L}}^{{-}1} \boldsymbol{\mathsf{d}}_{j,l} , \end{equation}

where each vector $\boldsymbol{\mathsf{d}}_{j,l} = ({\mathsf{d}}_{j,l,1},\ldots , {\mathsf{d}}_{j,l,3M})\in \mathbb {R}^{3M}$ is defined component-wise as ${\mathsf{d}}_{j,l,m+(i-1)\cdot M} = \int _\varOmega \delta (\boldsymbol {x}-\boldsymbol {x}_l)\boldsymbol {e}_j\boldsymbol {\cdot } \operatorname {\boldsymbol {\nabla }\times } (\phi _m(\boldsymbol {x})\boldsymbol {e}_i)$ for $m=1,\ldots , M$ and $i=1,2,3$. Hence, upon discretization, we may rewrite

(6.8)\begin{equation} \mathcal{J}({\boldsymbol{\theta}}) = \frac{1}{|S|} \sum_{\boldsymbol{x}_l\in S} \sum_{i,j=1}^3 ( \boldsymbol{\mathsf{f}}_{i,l}^\top\boldsymbol{\mathsf{B}}\boldsymbol{\mathsf{f}}_{j,l} - \tau_{ij}^\text{data} (\boldsymbol{x}_l) )^2 , \quad \text{where each}\ \boldsymbol{\mathsf{L}}\boldsymbol{\mathsf{f}}_{i,l} = \boldsymbol{\mathsf{d}}_{i,l} . \end{equation}

Because expression (6.8) is deterministic, equation (6.4) can be solved accurately and efficiently using a very wide variety of standard optimization software.

Remark 6.5 Owing to the fact that the loss function $\mathcal {J}({\boldsymbol {\theta }})$ may simply be written

(6.9)\begin{equation} \mathcal{J}({\boldsymbol{\theta}}) = \frac{1}{|S|} \sum_{\boldsymbol{x}_l\in S} \sum_{i,j=1}^3 (\mathbb{E}[ u_{i}({\boldsymbol{\theta}})u_{j}({\boldsymbol{\theta}})|_{\boldsymbol{x}_l} - \tau_{ij}^\text{data} (\boldsymbol{x}_l) ])^2 , \end{equation}

the optimization problem (6.4) can be solved with many stochastic optimization techniques commonly used in, e.g., the machine learning community. However, it is much more efficient to proceed by rewriting (6.4) as the deterministic optimization problem (6.8).

Alternatively, the optimization problem can be posed in the abstract setting of Bayesian inference. In this framework, the parameters are defined as random distributions (Stuart Reference Stuart2010).