1 Introduction
The flow of a viscous, incompressible Newtonian fluid is described by the (incompressible) NavierStokes equations:
Here is the velocity, is the pressure, which plays the role of a Lagrange multiplier to enforce the divergencefree constraint on the velocity, is the kinematic viscosity and represents the effects of external forces, for example, gravity and buoyancy. This system also needs to be supplied with initial and boundary conditions.
NavierStokes equations play a fundamental role in many applications and have been studied extensively for more than a century. The existence of globalintime weak solutions, both in two and three spatial dimensions, can be traced back to the pioneering works of Leray Le1934 and Hopf Ho1951 . In two spatial dimensions, the uniqueness of weak solutions of NSE has been established. However, in three spatial dimensions, the question of uniqueness remains unanswered and has been designated as a Millennium Prize Problem by the Clay Mathematics Institute.
Keeping the uniqueness question aside, it is known that, for small viscosity (or, equivalently, large Reynolds number
), i.e. when inertial forces are much stronger than viscous diffusion, the fluid flow is very sensitive to initial conditions and characterized by chaotic motions. Therefore, in a deterministic framework, measurement errors in the problem data could have a drastic affect on the solutions of NSE, and so, it is not physically meaningful to describe these turbulent flows as individual solutions. On the other hand, there is ample experimental evidence that the statistical observables, for example, mean and variance, can be inferred reliably for turbulent flows. Hence, for practical applications, it is more useful to consider the evolution of NSE under uncertainties on the problem data in a suitable probabilistic framework. In their seminal work
Fo1972 ; FPr1976, Foiaş and Prodi proposed the framework of statistical solutions of incompressible NavierStokes equations, in which, given a probability distribution on the initial velocity, velocity ensembles are evolved according to the NSE and described by a timeparameterized family of probability measures on the function space of initial velocity. The existence of such statistical solutions, both in two and three spatial dimensions, has been well established
Fo1972 ; FMRTe2001 . The global uniqueness in three spatial dimensions is an open problem. In two dimensions, statistical solutions are unique and they are defined as the pushforward of the probability measure on the initial velocity data. The authors in FMRTe2001 show that several results of the conventional theory of turbulence for NSE attributed to the groundbreaking work of Kolmogorov Ko1941a ; Ko1941b can be recovered with statistical solutions, thus, providing evidence of the importance of this solution framework.The computation of statistical solutions of NSE is a special case of Uncertainty Quantification (UQ) in Computational Fluid Dynamics (CFD), for which different methods are available in the literature (see BLMSc2013 and the references therein). Loosely speaking, the different methods used to solve UQ problems can be categorized into two classes: stochastic Galerkin methods and stochastic collocation methods. Stochastic Galerkin methods consider the RitzGalerkin formulation of the underlying PDEs in the stochastic parameter space, thus, these methods are highly intrusive and not amenable for implementation from a practical viewpoint (Su2015, , Chap. 12). On the other hand, stochastic collocation methods are wellsuited for large scale applications as they are nonintrusive (Su2015, , Chap. 13). For applications, Monte Carlo (MC) and QuasiMonte Carlo sampling can be seen as a subclass of stochastic collocation methods. Multilevel Monte Carlo methods have also received a lot of attention in research in recent times (see LMSc2016 and the references therein).
In the discussion above, the nonintrusive methods require a deterministic solver that numerically approximates the solution of the underlying system of PDEs (in our case the NSE). There are many numerical methods available in the literature that can efficiently approximate the NSE. For problems with periodic boundary, spectral methods are appealing because of their efficiency and high resolution Le2018 . Finite difference methods with Leray projection were first introduced in KLa1966 ; Ch1967 ; Ch1969 . Later, the authors in BCGl1989 proposed a finitevolume scheme that was an efficient variant of Ch1969 . In LMSc2016 , the authors develop a finitedifference scheme for the vorticity formulation of the NSE with uniformly stable bounds with respect to viscosity; this scheme is a variant of the one proposed in LTa1997 for incompressible Euler equations. Finite element methods are another class of methods to solve NSE GRa1986 . While standard mixed FEM are not pressurerobust, i.e. pressure approximation influences the velocity approximation, (weakly) divergencefree mixed methods are (see JLMNRe2017 and the references therein). In these weakly divergencefree mixed methods, if the velocity field belongs to the Hilbert space
, then the a priori error bounds are nonuniform with respect to viscosity. However, if H(div)conforming spaces are used instead, then the a priori velocity error estimates are robust with respect to viscosity
SLe2018 . Moreover, the divergencefree constraint on the velocity field is satisfied pointwise in the spatial domain. The authors in LSc2016 proposed a computationally efficient hybridized variant of H(div) FEM, called the H(div)HDG method, where they introduce hybridization for the tangential velocity components. This method was further improved in LLSc2018 and used to study turbulent flows in FKLLSc2019 ; SJLLLSc2019 .1.1 Contributions
We restrict ourselves to twodimensional open, bounded spatial domains and choose the framework of statistical solutions for NSE. Given the discussion above, to numerically approximate these solutions, we use Monte Carlo sampling for the stochastic space and H(div)FEM as the deterministic solver. We can easily use QMC instead of MC, but it is not necessary as our main focus is on the presence or absence of convergence of the approximate solutions and their statistics, not the rates. We choose the H(div) scheme and not the more efficient H(div)HDG scheme because of the ease of implementation. In the deterministic solver, we use the implicit Euler method for time integration as there is no CFL restriction on the timestep size (see MPRo1990 ; LSc2016 for higherorder time integration methods). To the best of my knowledge, this is the first computational effort in this direction.
Structure functions are an important observable in turbulence literature Fr1995 . In the context of incompressible Euler equations, a uniform decay of the structure functions ensures the convergence of approximate statistical solutions in the very recent work LMPa2021a . The behaviour of structure functions is also crucial in the study of the vanishing viscosity limit of incompressible NavierStokes in LMPa2021b . In this light, a novel contribution of this work is the development of an algorithm for approximating structure functions on unstructured meshes. We compute the structure functions only at the final time, which is a good indicator of their behaviour over the whole time interval as seen in LMPa2021a .
1.2 Outline
After covering some preliminaries, we state the initial boundary value problem (IBVP) and introduce the concept of weak solutions for NSE. In §4, we describe statistical solutions of NSE and define structure functions and Wasserstein distances. We present the deterministic H(div)conforming scheme in §5. In §6, we describe the algorithm for approximating structure functions on unstructured meshes. We present the results of our numerical experiments in §7 and conclusions in §8.
2 Preliminaries
Consider an open time interval , where denotes a finite time horizon, and let be an open and bounded polygon with Lipschitz boundary and with outward unit normal . Define the spacetime cylinder as the Cartesian product of the spatial domain and the time domain , i.e. . Let the scalar
and the vector
denote the time coordinate and the spatial Cartesian coordinates, respectively, then the spacetime Cartesian coordinates are .We highlight vector fields by bold face fonts and tensor fields by bold face fonts with an underline.
Given and a vector , the norm , we also denote the Euclidean norm by .
Given and two matrices with entries , and , we denote by the Frobenius inner product of the two matrices.
2.1 Function spaces
For any open, bounded domain , we use the Lebesgue spaces for scalarvalued functions with the associated norm , for . We also use the standard Hilbert spaces with associated norms and seminorms , for . For vectorvalued functions of size , we indicate these spaces as and . Spaces and norms for tensorvalued functions are indicated with bold font. We use the standard notation for the space of continuous functions on and , , for the space of functions on that are times differentiable with continuous th derivative. The space
(2.1) 
For a Hilbert space and , we use the standard notation for Bochner spaces and define .
2.2 Probability
Let be a probability space Sc2005
. For any random variable
, its mean and variance are given by(2.2a)  
(2.2b) 
We assume that is complete.
Given a topological space , we denote by the Borel algebra of and by the space of all probability measures on .
2.3 Meshes, mesh faces, averages and jumps
We choose a mesh of domain such that it is shaperegular and does not allow hanging nodes. We allow both simplicial meshes and meshes with quadrilateral elements, which we also refer to as triangular and quadrilateral meshes, respectively. The meshsize of is denoted by .
We define the following sets and unions for mesh faces:

denotes the set of all faces of ,

denotes the set of all interior faces of ,

denotes the set of all boundary faces of ,

the mesh skeleton ,

the interior mesh skeleton ,

the boundary mesh skeleton .
Further, we divide into two disjoint sets and that correspond to the boundaries with Dirichlet and outflow boundary conditions, respectively, in NSE (3.1). Using these sets, we also define the unions and .
For , we denote the outwardpointing unit normal vector on its boundary by . Let be distinct. Then, for a vector field and a tensor field , which are elementwise continuous on , we define the averages and jumps on the interior mesh face as follows:
(2.3)  
Here is the unit normal at face and we set by convention. We extend the definitions of the jumps and averages for boundary faces. Let such that , then we define
(2.4)  
We frequently drop the subscript in for convenience.
2.4 RaviartThomas spaces
We revisit the definition of RaviartThomas spaces and some important results associated with them, which are presented in (Ga2014, , Chapter 3) and (BBFo2013, , Chapter 2). We use these spaces to construct elementwise polynomial subspaces of in our numerical scheme.
Given , we denote by the space of polynomials of total degree at most , and we define the discrete polynomial space
(2.5) 
We denote RaviartThomas spaces of degree by , they are defined as follows:
Definition 2.1 (RaviartThomas spaces).
For any element , the local RaviartThomas space of degree is defined as
(2.6) 
Let denote the global RaviartThomas space, it is defined as
(2.7) 
Lemma 2.2 ((BBFo2013, , §2.5.2)).
Any function is continuous in the normal direction at all interior faces of the mesh, i.e.
(2.8) 
Moreover, the following holds true:
(2.9) 
3 NavierStokes equations
We assume that the boundary can be divided into Dirichlet boundary and outflow boundary , where either or may be empty. We also assume that and consist of entire segments of .
The IBVP for the timedependent incompressible NavierStokes equations in the divergence form reads as follows:
(3.1a)  
(3.1b)  
(3.1c)  
(3.1d)  
(3.1e) 
where is the velocity, is the pressure, is the kinematic viscosity, is the initial velocity, is the forcing, is the boundary data and is the identity tensor.
Remark 3.1.
For very large Reynolds numbers, a stabilized outflow boundary condition is needed to ensure numerical stability in the case of backflow at the outlet, we refer to BCBBr2018 for a review of the different stabilization methods available in the literature and in particular to DKCh2014 ; Do2015 for some popular choices. In the present work, see §7, our numerical experiments for problems with outflow preclude any backwardflow at the outlet. Therefore, the outflow condition (3.1e) is sufficient for our purpose.
3.1 Weak solutions
We assume homogeneous Dirichlet BCs, i.e. and . We choose the functional setting used in (FMRTe2001, , §II.5) and define the following Hilbert spaces:
(3.2)  
(3.3) 
On these spaces, we define the inner products
(3.4)  
(3.5) 
and the associated norms
(3.6) 
The weak formulation reads as: find such that, ,
(3.7) 
with the viscosity , the initial velocity and the forcing . This weak formulation can be traced back to the seminal work of Leray Le1933 ; Le1934 .
It is wellknown that in twodimensional spatial domains with sufficiently smooth boundaries, for any , any initial velocity and any forcing , weak solutions of NSE exist and they are unique. We state this result in the following (cf. (FMRTe2001, , Chapter II, Theorem 7.3 and Remark 7.2)):
Theorem 3.2 (Existence and uniqueness of weak solutions of NSE).
Assume that the spatial domain has boundary , and that
(3.8) 
Then, for every , there exists a unique solution of (3.7) such that
(3.9) 
and, for all , it satisfies the energy equation
Moreover, there exists a continuous solution operator such that .
4 Statistical solutions
The discussion in this section is based on the monograph (FMRTe2001, , Chapter V) and it partially uses the notation prescribed in (LMSc2016, , §2.2). We consider statistical solutions of NSE (3.1) with fixed, noslip boundary, i.e. and , where the boundary is of class .
Given a probability distribution on the initial velocity data, i.e.
the main idea of statistical solutions, as introduced by FoiaşProdi Fo1972 ; FPr1976 , is to describe the evolution of velocity ensembles by a timeparameterized family of probability measures , such that (up to null sets), and for every , . This solution framework includes the individual weak solutions of NSE as a special case with the probability measure , where is the Dirac measure (also known as the unit mass) at , cf. (Sc2005, , §4.7).
From Theorem 3.2, we know that the individual weak solutions of NSE are unique and the solution operator exists. As a result, for any timeindependent forcing , unique solutions are given by transporting the initial probability measure under the solution operator , i.e.
(4.1) 
Concretely, we state this existence and uniqueness result in the following (cf. (FMRTe2001, , Chapter V, Theorems 1.1 and 1.2)):
Theorem 4.1.
Let be a probability measure on with finite kinetic energy, i.e.
Assume that the forcing term . Then, for every and homogeneous Dirichlet boundary conditions, there exists a statistical solution of the incompressible NavierStokes equations (3.1) on .
Moreover, if has bounded support in and is independent of time, then the statistical solution is unique and is given by , where is the solution operator of the incompressible NavierStokes equations, cf. Theorem 3.2.
4.1 Structure functions
For a solution , other than its common statistics like mean and variance, its structure functions are of great interest to us. The structure functions are defined as follows (we refer to Ly2020 ):
Definition 4.2 (Structure functions).
Let , and . Let be a random field such that , for any , and be the probability distribution of at time . Then, the structure function associated with is given by:
(4.2) 
where denotes the ball of radius centered at a point .
We use the structure functions as a measure of the regularity of the velocity field , and we investigate the scaling behaviour of with respect to the offset in our numerical experiments.
4.2 Wasserstein distances
In the next section, we describe a numerical method to approximate statistical solutions with ensembles of approximate individual velocity solutions. To study the convergence of these velocity ensembles, it is important to define a notion of distance between two ensembles. For this purpose, we use Wasserstein distance (cf. (FLMi2017, , Definition 2.2)):
Definition 4.4.
Let be a separable Banach space, and let be probability measures on with finite
th moments, i.e.
and . Then, the Wasserstein distance between and is defined aswhere the infimum is taken over the set of all transport plans from to , i.e. those that satisfy
Here denotes the space of bounded, continuous, realvalued functionals on .
Given , we can define Wasserstein distance for the point distribution at and the point correlation at , respectively, as:
(4.3)  
(4.4) 
For the whole domain , we can define
(4.5)  
(4.6) 
In this work, we always use Wasserstein distances and , so we simply refer to them as Wasserstein distances and drop the subscripts.
5 Numerical approximation of statistical solutions
In this section, our goal is to approximate the statistical solutions described in the previous section. To this end, we assume that initial probability measure is given as the law of a random field defined on the underlying probability space . This allows us to draw random samples from the initial distribution, which can then be evolved using a numerical solution operator. We estimate the statistical solutions using (5.1). Concretely, this Monte Carlo type method is described in Algorithm 1.
(5.1) 
Here is a numerical approximation of the operator such that .
To approximate an individual solution, we choose an H(div)conforming method based on the work of Sequeira et al. GSSe2017 and Cockburn et al. CKSc2007 ; this scheme has also been considered in SLe2018 ; SLu2018 .
For spatial discretisation §5.1, we use the H(div)conforming scheme designed in GSSe2017 for the incompressible Euler equations along with the symmetric interior penalty discretisation for the viscous diffusion given in CKSc2007 and (DEr2012, , §4.2.2). This semidiscrete scheme is pressurerobust and semirobust, we refer to SLe2018 for details, which means that the velocity error bounds have no explicit dependence on the pressure approximation and the Reynolds number. To obtain a fully discrete scheme §5.2, we combine the H(div) spatial discretisation with implicit Euler timestepping as formulated in GSSe2017 .
5.1 H(div)conforming spatial discretisation
We use the function space , which is defined as
(5.2) 
for the velocity field. Here the space is given by (2.7). In the case of purely Dirichlet boundary, we denote by .We use the function space (2.5) for the pressure field.
Remark 5.1.
For purely Dirichlet boundary, i.e. , if is a solution of (3.1) such that , then , for any constant , is also a solution. Therefore, an additional constraint is needed to ensure a unique pressure solution. A common choice is to impose vanishing mean for the pressure (cf., for example, GSSe2017 ; JLMNRe2017 ) using the function space
(5.3a)  
(5.3b) 
Given initial velocity data , let denote the orthogonal projection of into . Then, the H(div)conforming spatial discretisation of the IBVP (3.1) reads as follows:
Given the initial data , the forcing and the boundary data , find and such that, for all and for all ,
(5.4a)  
(5.4b) 
where the mass bilinear form
(5.5) 
the convection trilinear form with upwind flux
(5.6)  
where (cf. (2.8)), the symmetric interior penalty diffusion bilinear form
(5.7)  
where is the jump penalization parameter and is the size of face , the divergence bilinear form
(5.8) 
and the linear form
(5.9)  
Remark 5.2.
In the formulation (5.4), the outflow boundary conditions are imposed weakly through the boundary integrals on that appear in the bilinear forms and .
Remark 5.3.
For homogeneous Dirichlet boundary conditions, the scheme (5.4) is the same as the one analyzed in SLe2018 , except for the slight difference in the convection trilinear term that does not affect the error estimates. Therefore, for , from Theorem 5.3, Theorem 5.6 and Corollary 5.9 in SLe2018 , we deduce that the error for the discrete solution of (5.4) is bounded by , , for with sufficient regularity.
5.2 Fullydiscrete scheme
As mentioned, we combine the semidiscrete formulation (5.4) with implicit Euler timestepping to obtain a fully discrete scheme §5.2, by following the procedure described in GSSe2017 .
For , let be a set of points in the time interval such that with . Given the solution of the IBVP (3.1), we define the notation and , and we denote by and the numerical approximations of and , respectively.
The Taylor series expansion of at is given by
(5.10) 
which yields the time derivative
(5.11) 
with the truncation error
(5.12) 
Substituting (5.11) in (3.1), at any discrete time , we obtain
which can be written equivalently as
(5.13)  
Using (5.2) and (5.12), we deduce that the term . Dropping the truncation term in (5.13), the remaining equation is linear with respect to the variable . We introduce H(div)conforming spatial discretisation in the truncated system (5.13) to arrive at a fully discrete method, which reads as follows:
Given the initial data , the forcing and the boundary data , find and such that, for all and for all ,
(5.14a)  
(5.14b) 
for , with .
Remark 5.4.
The discrete formulation (5.14) satisfies the following:
for . Then, for the forcing , we have
Therefore, the scheme is stable.
Remark 5.5.
For viscosity , to the best of my knowledge, error estimates of our numerical scheme are not available, and it is beyond the scope of this work to derive them. However, based on the reasoning in the proofs of (SLu2018, , Lemma 5.5, Theorem 5.6) and (SLe2018, , Lemma 3.4, Theorem 3.6), we expect the velocity error to be bounded by .
6 Implementation details
In this section, we discuss important details about our implementation of the Monte Carlo algorithm 1, which we refer to as MCFEM solver; here all the samples are evolved in parallel using the selfscheduling algorithm (GLSk2014, , §3.6). In particular, we describe the deterministic solver that is used for a single Monte Carlo sample and present the algorithms used to approximate structure functions and Wasserstein distances.
We denote vectors by a lowercase font with an underline and matrices by uppercase font or bold font with an underline, only for the discussion on the deterministic solver. Given bases of and , we denote the corresponding vectors of basis coefficients by and . Then, the finite element assembly for the discrete formulation (5.14) yields the following blockstructured linear system of equations:
(6.1) 
where the matrices , , and are assembled from the mass form , the convection form , the diffusion form and the divergence form , respectively. Here the righthand side vector is assembled from the linear form . Note that we must use the Piola transformation in our finite element assembly.
We need to assemble the matrices , and only once, but the matrix has to be reassembled for each discrete time because of its dependence on the velocity solution of the previous time step . Using direct solvers for the linear system (6.1) would involve factorizing the reassembled system matrix at each time step, which is computationally expensive. Therefore, we use suitable iterative solvers with preconditioning. Rescaling the pressure solution as
Comments
There are no comments yet.