# Analysis and approximation of a fractional Laplacian-based closure model for turbulent flows and its connection to Richardson pair dispersion^{†}^{†}thanks: Supported by the US National Science Foundation grant DMS-1315259 and the US Air Force Office of Scientific Research grant FA9550-15-1-0001.

###### Abstract

We study a turbulence closure model in which the fractional Laplacian of the velocity field represents the turbulence diffusivity. We investigate the energy spectrum of the model by applying Pao’s energy transfer theory. For the case , the corresponding power law of the energy spectrum in the inertial range has a correction exponent on the regular Kolmogorov -5/3 scaling exponent. For this case, this model represents Richardson’s particle pair-distance superdiffusion of a fully developed homogeneous turbulent flow as well as Lévy jumps that lead to the superdiffusion. For other values of , the power law of the energy spectrum is consistent with the regular Kolmogorov -5/3 scaling exponent. We also propose and study a modular time-stepping algorithm in semi-discretized form. The algorithm is minimally intrusive to a given legacy code for solving Navier-Stokes equations by decoupling the local part and nonlocal part of the equations for the unknowns. We prove the algorithm is unconditionally stable and unconditionally, first-order convergent. We also derive error estimates for full discretizations of the model which, in addition to the time stepping algorithm, involves a finite element spatial discretization and a domain truncation approximation to the range of the fractional Laplacian.

Key words. turbulence modeling, fractional Laplacians, nonlocal closure, Navier-Stokes equations, Richardson pair dispersion, finite elements

## 1 Introduction

Nonlocal models have attracted intensive research interests in recent years due to their ability to model phenomena that cannot be correctly described by classical partial differential equation models. Many advances have been made in various scientific and engineering areas including continuum mechanics [34], graph theory [25], image denoising [6], machine learning [32], and phase transitions [5]. In particular, fractional derivative models have been found to be effective in modeling anomalous diffusion processes [27, 26]. In this work, we study a new closure model based on the fractional Laplacian operator that accounts for the anomalous diffusion (superdiffusion in this case) that arises in fully-developed turbulent fluid flows.

Turbulence modeling remains one of the most challenging scientific problems. Despite the fact that the governing equations for turbulence have been known since 1845, a full understanding of turbulence is still far from complete due to its extremely complex behavior and chaotic nature. The wide range of scales present in turbulence results in a very high computational complexity and renders direct numerical simulations infeasible even with modern supercomputers. Thus, turbulence models are introduced to predict the mean flow and coherent structures with the effects of the turbulence on the mean flow being modeled. The mean flow and the smaller scales of turbulence interact through a quantity referred to as the Reynolds stress that appears in the evolution equations of the mean flow and which, to close the system, must be replaced by terms that are solely dependent on the mean flow; this is the closure problems of turbulence. It is worth noting that all turbulence models inevitably invoke additional heuristic hypotheses and thus tend to work only for a narrow class of problems. In this paper, we consider a class of nonlocal operators, namely fractional Laplacian operators, as a turbulent closure model [7]. These operators have a deep connection with Lévy jump processes in probability theory and corresponding superdiffusion behavior in turbulent flows.

The closure model we consider is given by

\hb@xt@.01(1.1) |

for , where denotes a bounded, open domain and a temporal interval of interest. The fractional Laplacian operator is most often defined in terms of Fourier transforms as

where denotes the Fourier transform and its inverse. For , an equivalent definition [3] is

\hb@xt@.01(1.2) |

where is a normalizing constant.

We first explore, in §LABEL:energytransfer, the energy spectrum of the model (LABEL:model) by applying Pao’s transfer theory [29]. Following the methodology developed in [28] for the analysis of a family of approximate deconvolution models of turbulence, we derive an expression for the long-time averaged energy distribution. The results show that if , the corresponding power law of the energy spectrum in the inertial range has a correction to the well-known Kolmogorov scaling exponent. For this case, this model (LABEL:model) corresponds to the Richardson particle pair-distance superdiffusion of a fully developed homogeneous turbulent flow (see [31] and also [22, 21, 15]) and also to a Lévy jump process that lead to superdiffusion. For other values of in , the power law of the energy spectrum is consistent with the standard Kolmogorov scaling exponent.

The use of the fractional Laplacian operator results in a dense matrix that requires different types of linear solvers than the ones employed for the sparse matrices that arise when solving discretized local turbulent models. Because many industrial flow codes are highly optimized and extensively calibrated, a direct implementation of the proposed model encounters many technical issues and requires substantial code modifications. To ease the implementation process for the model we consider and reduce the required effort, we propose, in §LABEL:imexmethod, a novel modular algorithm that splits the local and nonlocal parts of the equations so that only minimal changes need be done to legacy codes. The algorithm consists of two steps. The first is to solve the Navier-Stokes equations with a modified right-hand side so that a legacy code can be easily modified without changing linear solvers used or the manner in which matrices are stored. The second step is a post-processing step that involves solving a discretized fractional Laplacian problem. This step can be added to the legacy code as a separate routine. We study our modular algorithm based on a first-order time-stepping method also given in §LABEL:imexmethod, proving that the algorithm is unconditionally stable and unconditionally first-order convergent. This modular algorithm can be extended to higher-order time-stepping methods.

The model we study employs a standard definition, i.e., (LABEL:ffll), of the fractional Laplacian operator, instead of some variants that are defined on the bounded domain for which, e.g., the integral appearing in (LABEL:ffll) is replaced by an integral over . Thus, in our model, although the domain of the fractional Laplacian operator appearing in (LABEL:model) is the bounded domain , its range is the infinite domain . In practice, however, the integral in (LABEL:ffll) is approximated by an integral over a finite domain strictly containing the given domain . In §LABEL:truncatedinteractions, for a fixed bounded domain , we derive an estimate for the error incurred by truncation as the extent of the truncated containment domain increases, in particular showing that solutions of the truncated domain problem converge to those of (LABEL:model).

In §LABEL:sec:femfem, we complete our study by defining and analyzing full discretizations of the problem (LABEL:model) for which finite element spatial discretizations are added to the time-stepping methods of §LABEL:imexmethod and the domain truncation of §LABEL:truncatedinteractions.

## 2 Energy transfer

We investigate the energy spectrum of the new model based on the energy transfer theory of Pao [29, 30]. We thus consider the Navier-Stokes equations in a periodic box in :

\hb@xt@.01(2.1) |

where, for , , , or ,

with , , denoting the Cartesian unit vectors. Then, the fluid velocity can be expanded in Fourier series as

and its associated kinetic energy is given by

where with , being non-negative integers and . We partition the kinetic energy into wave number shells given by

where , so that the total energy is then given by

An evolution equation for can be derived by taking inner product of (LABEL:periodicNSE) with a single Fourier mode and then summing over all modes; see Davidson [8] or Pope [9]. Using the Kronecker delta, satisfies

where denotes the complex conjugate of . Note that for each , is a complex vector that satisfies conjugate symmetry, i.e. . Define the energy transfer function by

where

Because fully developed, homogeneous, isotropic turbulence is characterized by a wide range of persistent scales, in transfer theory is considered as a continuous variable. We then redefine the energy transfer function as

We further assume that for all the energy is input into the modes by smooth, persistent body forces, i.e., , where is fixed, which is a representative large scale velocity. Based on all the assumptions made above, we have the evolution equation for the kinetic energy for a given given by

This system is not closed. So far the most successful closure model is Pao’s model [30] given by

where and is the Kolmogorov constant. Based on Pao’s closure model, we now consider the problem

We define the long-time averaged energy distribution as

which then satisfies

along with . This equation can be easily solved and the solution is given by

\hb@xt@.01(2.2) |

where

### 2.1 The inertial range energy spectrum

In the inertial range, the viscous dissipation effect is negligible because is small. Then, over this range, the expression (LABEL:Spectrum) for reduces to

This shows that if the exponent of the fractional Laplacian is equal to , the corresponding power law of the energy spectrum in the inertial range has a deviation from the exponent of the regular Kolmogorov scaling exponent, whereas for other values of the power law of the energy spectrum is consistent with the Kolmogorov theory. The Kolmogorov scaling theory (often referred to as the “K41 theory”) is the most celebrated turbulence theory and is supported by much experimental evidence from atmospheric and oceanographic turbulence at sufficiently high Reynolds number [14, 8]. However, small deviations from the scaling exponent have also been observed in various turbulence experiments [2, 24, 17, 14]. These deviations, although small in the spectrum, considerably affect higher-order statistics. There have been many theoretical attempts to modify the exponent in the power law. Actually, Kolmogorov himself first proposed a modification of the exponent [20]. Most of the theories developed concern the intermittency in the inertial range and various intermittency models have been built to try to fit experimental data such as the -model; see [14] for a review. However, as far as we know, there exists no partial differential equation turbulence model that is able to correct the exponent of the Kolmogorov scaling exponent as does the model we consider.

The fractional Laplacian is the generator of -stable Lévy processes in probability theory. The special case of that leads to a correction exponent in the power law of energy spectrum corresponds to the -stable Lévy process. This has an interesting connection with the Richardson’s particle pair-distance superdiffusion in a fully developed homogeneous turbulence for which so that the displacement increment also obeys the -stable Lévy distribution. Thus, the fractional Laplacian with which we use in our model actually introduces the corresponding Lévy flight mechanism into the system and represents Richardson’s turbulence superdiffusion.

### 2.2 The dissipation range energy spectrum

In the dissipation range, viscous dissipation is dominant and removes energy from the system. We rewrite (LABEL:Spectrum) as

decays exponentially in the dissipation range, which is consistent with the Kolmogorov scaling theory. In the dissipation range, . Then, if , we have and and thus . Similarly, if , we have and and . So for all , our closure model results in enhanced exponential decay in the dissipation range.

## 3 First-order IMEX time-stepping methods

Most turbulent flows of engineering and scientific interest occur in bounded flow regions. Thus, we are more interested in computing turbulent flows on bounded domains. Accordingly, , , denotes an open, bounded domain and consider the problem

\hb@xt@.01(3.1) |

Here we do not change the definition of the fractional Laplacian as defined in (LABEL:ffll) but merely restrict its range to the bounded domain . Note that because the domain of integration in (LABEL:ffll) is , we impose, in (LABEL:eq:model), the constraint on the complement domain instead of on the boundary of as is done in the PDE setting.

We first consider the simple first-order implicit-explicit (IMEX) Euler time-stepping scheme given as follows. Given , find and satifying

\hb@xt@.01(3.2) |

We prove, in §LABEL:imex1, that this time-stepping scheme is unconditionally stable and first-order accurate.

Because of its implicit treatment of the fractional Laplacian term, the scheme (LABEL:IMEX-Euler) requires the solution of a dense linear system at each time step. Having to also handle the Navier-Stokes terms makes for an even greater computational challenge. Thus, it is tempting to lag the fractional Laplacian term to the previous time step; however, this leads to serious stability issues so that that term has to be treated implicitly. However, there does exist a way to split the equations so that one can still solve a (modified) local momentum equation involving the usual sparse matrices and subsequently correct the solution by solving a linear nonlocal fractional Laplacian equation. Specifically, we propose to modify the IMEX Euler scheme (LABEL:IMEX-Euler) into the following modular algorithm.

###### Algorithm 3.1 (Modular IMEX Euler)

Stage 1: Given in , find in satisfying

\hb@xt@.01(3.3) |

Stage 2: Given and in , find in satisfying

\hb@xt@.01(3.4) |

In §LABEL:imex2, we prove that this time-stepping scheme is also unconditionally stable and first-order accurate. However, we can now solve the Stage 1 problem using a legacy Navier-Stokes code with the only modification necessary being in the construction of the right-hand side. Then, in the second stage, one solves a “Poisson” problem for the fractional Laplacian operator which involves a symmetric, positive definite, albeit dense linear system. This two-stage algorithm, although involving two linear system solves per time step, requires, compared to the algorithm given in (LABEL:IMEX-Euler), much less coding effort and introduces efficiencies not possible for the scheme (LABEL:IMEX-Euler).

### 3.1 Preliminaries

We first recall that for , the fractional Sobolev space is defined as

which is an intermediary Banach space between and , equipped with the natural norm

For , we have and

We denote the Gagliardo (semi)-norm of by

Define

Let denote the velocity space and be the pressure space, defined by

###### Remark 3.2

Generally, and thus . This seems to cause some difficulties in finding a suitable finite element space that is a subset of to approximate solutions. However, in practical numerical simulations, one cannot integrate over all of . So the domain of integration must be restricted. Here, for a , we assume is the interaction domain defined by

\hb@xt@.01(3.5) |

and one only integrates over . Let

Then, to use finite element methods, one only needs to find a finite element space that is a subset of . Thus, the usual finite element spaces, such as continuous piecewise-quadratic elements, can be employed. It is shown in §LABEL:truncatedinteractions that the error incurred by domain truncation is of .

### 3.2 Analysis for the algorithm (LABEL:IMEX-Euler)

We prove that the time-stepping scheme (LABEL:IMEX-Euler) is unconditionally stable and first-order convergent. We do not provide the proofs because they are similar to those for the theorems considered in §LABEL:imex2.

###### Theorem 3.3 (Unconditional stability of time-stepping scheme (LABEL:IMEX-Euler))

The IMEX Euler scheme (LABEL:IMEX-Euler) is unconditionally, long-time stable. Specifically, for any , we have that

To analyze the rate of convergence of the approximation we assume the following regularity on the exact solution and the body force of (LABEL:eq:model):

Let and . Denote . We introduce the discrete norms

Let denote the error at the time between the exact solution of (LABEL:eq:model) and the approximation obtained using the time-stepping scheme (LABEL:IMEX-Euler). Then, we have the following error estimate.

###### Theorem 3.4 (Unconditional convergence of algorithm (LABEL:IMEX-Euler))

For any , there exits a positive constant independent of the time step such that

\hb@xt@.01(3.6) | ||||

### 3.3 Analysis for the modular algorithm (LABEL:step1)-(LABEL:step2)

We first prove that the modular time-stepping scheme (LABEL:step1)-(LABEL:step2) is unconditionally stable. We emphasize that the coefficient multiplying the fractional Laplacian term in (LABEL:step2) is essential for the unconditional stability of the modular algorithm.

###### Theorem 3.5 (Unconditional stability of the modular algorithm (LABEL:step1)-(LABEL:step2))

The modular algorithm (LABEL:step1)-(LABEL:step2) is unconditionally, long-time stable. Specifically, for any , we have that

\hb@xt@.01(3.7) | ||||

Proof. The proof is provided in Section LABEL:thm3.5proof.

We next prove that the modular algorithm (LABEL:step1)-(LABEL:step2) is also a first-order in time scheme. To analyze the rate of convergence of the approximation we assume that the true solution and and the body force of (LABEL:eq:model) have regularity given by

Let and denote the solution to Stages 1 and 2 in algorithm (LABEL:step1)-(LABEL:step2), respectively. Denote and . Then, we have the following error estimate.

###### Theorem 3.6 (Unconditional first-order convergence of the modular algorithm (LABEL:step1)-(LABEL:step2))

Consider the modular algorithm (LABEL:step1)-(LABEL:step2). Assuming , then for any , there is a positive constant independent of time step and mesh size such that

\hb@xt@.01(3.8) | ||||

Proof. The proof is provided in Section LABEL:thm3.6proof.

## 4 Truncation of interactions

As already stated in Remark LABEL:rem2, it is impractical to integrate over all of . Furthermore, it is reasonable to assume that the velocity and pressure fields at points that are far away from a point in the bounded domain of interest have negligible effect on the velocity and pressure fields at points close to . This is clearly true for the integrand kernels we consider in our model. Thus, we consider limiting the extent of nonlocal interactions for a point to the ball of radius centered at denoted by

and define the interaction domain as in (LABEL:lambda). We then consider the truncated variational problem

\hb@xt@.01(4.1) |

The nonlocal term in (LABEL:prob) now involves a finite integration domain, in contrast to that in (LABEL:eq:model) that features an infinite integration domain.

We recall that for , the fractional Sobolev space is defined as

equipped with the natural norm

For , we have

We define the following norm of by

We define the velocity and pressure spaces

respectively, where

and

Let denote divergence free velocity space

In the next theorem, we analyze the error due to domain truncation and the rate of convergence with respect to . We assume that the true solution of (LABEL:eq:model) has regularity given by

Also, we denote the error due to truncation by , where and denote the solutions of (LABEL:eq:model) and (LABEL:prob), respectively.^{1}^{1}1For the simpler setting of the fractional Laplacian Poisson problem, the error due to truncation was considered in [10].

###### Theorem 4.1 (Error due to domain truncation)

Let . Then, for any , there exists a positive constant independent of the time step and the truncation radius such that

\hb@xt@.01(4.2) | ||||

In particular, if , we have

\hb@xt@.01(4.3) | ||||

Proof. The proof is provided in Section LABEL:thm4.1proof.

## 5 Finite element approximations

We denote by and conforming velocity and pressure finite element spaces, respectively, based on an edge-to-edge triangulation of with maximum triangle diameter and with consisting of triangle vertices and/or edges. We assume that and satisfy the usual discrete inf-sup condition and the approximation properties, [23]