A Subgrid Model for the Time-Dependent Navier-Stokes Equations

We propose a stabilized subgrid finite-element method for the two-dimensional 2D nonstationary incompressible Naver-Stokes equation NSE . This method yields a subgrid eddy viscosity which does not act on the large flow structures. The proposed eddy viscous term is constructed by a fluctuation operator based on an L2-projection. The fluctuation operator can be implemented by the L2-projection from high-order interpolation finite-element spaces to the low-order interpolation finite-element spaces. In this paper, P2/P1 mixed finite-element spaces are adopted to implement the calculation and the analysis. The error analysis is given based on some regular assumptions. Finally, in the part of numerical tests, the numerical computations show that the numerical results agree with theoretical analysis very well. Meanwhile, the numerical investigations demonstrate that the proposed subgrid is very effective for high Reynolds number fluid flows and superior to other referred subgrid models.


Introduction
In this paper, we focus on formulating a subgrid eddy viscosity method for the nonstationary incompressible Navier-Stokes equations.For the subgrid method, we must admit that there exists a scale separation between large and small scales and this model can be viewed as a viscous correction for large-scale fluid flows.Meanwhile, for laminal fluid flows, the added subgrid viscosity term should not affect the large-scale structures of fluid flow field and should tend to vanish.This kind of subgrid method is just a flexible and effective method for high Reynolds number fluid flows.
It is well known that for most problems of fluid flows, the numerical algorithms capturing all scales of fluid flows are impossible and there exist several scales that span from the large-scales to the small Kolmogorov scales which can hardly be resolved by state-of-theart computers for most engineering problems very efficiently.For the convection dominating fluid flows, we need to consider the dispersive effects of unresolved scales on resolved We introduce the following functional settings: Ω q dx 0 .

2.2
We denote by For convenience, we introduce the following bilinear form a •, • on X × X: a u, u u, u , ∀u ∈ X, 2.4 and d •, • on X × Q defined by d v, q − v, ∇q q, div v , ∀v ∈ X, q ∈ Q.

2.5
The trilinear term is defined by which is the skew-symmetric form of the convective term.It is easy to gain b u; v, w −b u; w, v .

2.7
Also, we have the following estimates 9 : where C is a positive constant depending only on the domain Ω, which stands for different values at different occurrences.The weak formulation of 2.1 reads as follows.Find u, p ∈ X, Q such that and u 0, x u 0 x ∈ X.

Advances in Numerical Analysis
For the finite-element analysis, we need some regularity assumptions of the Navier-Stokes equations and the solution.
Theorem 2.1 7, 10 .One assumes that f ∈ L 2 0, T; L 2 2 , u 0 ∈ X, 2.10 and that 2.9 possesses a solution u, p with Note.These assumptions imply that the solution of 2.9 is unique.For simplicity, let f f h .In addition, we assume that Ω has a polygonal boundary such that no boundary approximation in the application of the finite-element method becomes necessary.

Discretization of Navier-Stokes Equations and Subgrid Model
We give a family τ h , which is a partition of Ω into triangles or quadrilaterals, assumed to be regular in the usual sense 11 .The diameter of the cell K is denoted by h K .The mesh parameter h describes the maximum diameter of the cells K ∈ τ h .
We introduce the finite-dimensional subspace X h and Q h :

3.1
We define the space of discretely divergence free functions denoted by V h : and let the following approximation properties hold 12 :

3.4
Meanwhile, the velocity-pressure pair in X h , Q h satisfies the following discrete inf-sup condition 13 : inf Remark 3.1 14, 15 .Let Π : Q → R 0 be the standard L 2 -projection with the following properties: where For a tensor T, ΠT denotes that Π acts on each component of this tensor.
We know that for high Reynolds number fluid flows, when the fluid convection dominates fluid flow fields, under the finite-resolution of meshes, the flow becomes very instable.When the mesh scales cannot resolve the smallest scale in fluid flows, we must add some term into Navier-Stokes equations to smear out the effect from the unresolve scales.Here, we chose the following subgrid stabilization term to control the effect from the unresolved scales: where α is the user-selected stabilization parameter and typically, α h s s is a real number .
The analogous stabilization is used to circumvent the pressure stabilization LBB condition for Stokes problems 16 .Since Π is an L 2 -projection, it follows for v h ∈ X h and ∇v h 0 > 0 that

3.8
In addition, from 0 ≤ Π∇v h 0 ≤ ∇v h 0 , it follows that 0 ≤ α add v h ≤ α. 3.9 Note.If v h depends on t, then α add v h too.From 3.9 it follows that α add v h t, • ∈ L ∞ 0, T if α is bounded almost everywhere in the time interval.If ∇v h 0 0, then v h 0 since v h ∈ X h .In this case, we set α add v h 0. The analogous formula was proposed to define a reduced Reynolds number for a variational multiscale method of the Navier-Stokes equations 7 .

Advances in Numerical Analysis
Using the stabilization term M u h , v h , we give the following stabilization finiteelement discretization form of the variational form 2.9 .
Find u h , q h ∈ X h , Q h satisfying

3.10
Under the inf-sup condition 3.5 , formulation 3.10 is equivalent to the following problem. Find

Error Analysis
The proof of the finite-element error estimate uses an approach by John and Kaya 7 and Heywood and Rannacher 12, 14 .The theoretical results of error analysis are classical 7 .
We first show an outline 7 .
1 Prove stability of u and u h , that is, certain norms of u and u h are bounded by the data of the problem: f, u 0 , ν.
2 Derive an error equation by subtracting 3.11 from 2.9 for test functions from V h .Split the error e into an approximation term η and a remainder φ h : where u h ∈ V h is a projection of u which fulfills the estimates 3.3 .Then take φ h ∈ V h as test function in the error equation.
3 Estimate the right hand side of the error equation which has the following form where all functions are nonnegative for all t ∈ 0, T .
4 Apply Gronwall's lemma to 4.2 , that is, derive all the functions in 4.2 belong to L 1 0, T and get an estimate of φ h .
5 Derive the error estimate of e by applying the triangle inequality to 4.1 .
Along these lines, the estimate will be proved 7 .This error estimate uses the parameter α add defined in 3.8 .The proving method is based on the classical scheme 7 .However, we still give the details of theoretical analysis for completeness.Firstly, we present the stability of u and u h .Proof.The proof for the u h and u is very similar.We will show the result for u h .Setting v h u h in 3.11 , use the skew-symmetric form of b •, •, • , 3.8 , and integrate over 0, t with t ≤ T , we get

4.3
Here u 0,h is the value of u h at t 0.
Subtraction of the last term and the regularity 2.10 gives Now, Step 2 of the proof is carried out by a strightforward calculation.We obtain 1 2

Advances in Numerical Analysis
To get the form of 4.2 , the terms on the right-hand side of 4.4 have to be estimated.Using the Cauchy-Schwarz inequality, Young's inequality, and 3.8 , we get

4.5
The trilinear term is decomposed into three terms as follows: Using 2.8 and Young's inequality, we have

4.7
Advances in Numerical Analysis 9 Collecting all the above terms, we obtain 1 2

4.8
We define Re red :

4.9
It follows that Re red is smaller or equal than Re.Using 3.9 finishes Step 3 of the proof:

4.10
To perform Step 4 of the proof, we have to study the L 1 0, T -regularity of the terms in 4.10 .Let t ∈ 0, T be arbitrary.Using Poincaré's inequality, the Cauchy-Schwarz inequality, then 2.11 and 3.4 , we get

Advances in Numerical Analysis
Similarly by H ölders inequality, Lemma 4.1, and 3.4 , we have

4.12
The L 1 0, T -regularity of other terms is a direct consequence of 2.11 , 3.3 , and 3.4 .The application of Gronwall's inequality and the last step of the proof give the following theorem.
Theorem 4.2.Let u, p ∈ X × Q be the solution of 2.9 and let u h ∈ V h be the solution of 3.11 .Let the regularity assumptions 2.11 be fulfilled, let u h be a projection of u into V h such that η fulfills 3.3 .Let the reduced Reynolds number Re red be defined in 4.9 .Then, the error e u − u h satisfies, for T ≥ 0, In particular,

4.15
Remark 4.4.A finite-element error estimate for the L 2 Ω -error in the pressure can also be derived following Heywood and Rannacher 12 .Using the inf-sup condition 3.5 and the estimates of 3.3 , the pressure error can be estimated by approximation errors and the velocity error u − u h t 0 .Theorem 4.2 finishes the error estimate for the pressure.Since the analysis is lengthy and follows closely 12 , we will not present it here.

Numerical Tests
Firstly, we give the algorithm used to deal with the nonlinear term and the subgrid eddy viscosity term.Given where and Δt is the time interval.Also, u n h is the finite-element approximate solutions of u h t n nΔt 1 < n < N .The fixed point iteration is adopted to solve 5.1 based on the Newtonian method 17 and the iteration tolerance is set to equal 10 −8 .For calculating the subgrid term α I − Π ∇u n h , I − Π ∇v h , the subgrid term is rewritten as follows: In order to reduce the computational work, the following linear time extrapolation method is used to calculate Π∇u n h :

Validate Convergence Rate by Taylor Vortex
It is essential to investigate the subgrid model 3.7 for low viscosity fluid flow and validate the flexibility and convergence orders of this model.So, we need to choose a true solution.
We consider 2.1 on the domain Ω 0, 1 × 0, 1 , with a body force obtained such that the following true solution is given by u u 1 , u 2 :  mean h denotes the averaging mesh scale.The time step 0.001 and the final time t 2.
In Table 1, the relative errors and absolute errors are given by different mesh scale h.In Figures 1 and 2, L 2 H 1 convergence orders are given by two log-log plots.From these two figures, L 2 and H 1 convergence orders are equal to 2.31194 and 1.34385, respectively.The calculated convergence orders coincide with the theoretical analysis in Corollary The expected convergence orders for the velocity in L 2 and H 1 are the second-order and the first-order, respectively.The computational convergence orders are a little higher than the expected orders.These results are mainly attributed to a specific test case which has a good regularity of the solutions.

Lid-Driven Cavity Flows
In this part, the comparisons among three kinds of subgrid models Current model, Guermond-Marra-Quartapelle GMQ 5 , and Kaya-Layton-Riviere KLR model 6, 7 are investigated for lid-driven cavity flows.These investigations will address the actions of subgrid models on large-scale flow structures.Now, we introduce the GMQ model  and KLR model.The GMQ model is based on the two-level Lagrange finite-element setting, which is defined by 5 where |K| is the measure of K, u H h 1 − P H u h , and v H h 1 − P H v 5 for the details of the functional settings and the definition of the operator P H The KLR model is also based on two-level finite-element space, which is defined by 6, 7 where g H is the L 2 projection of ∇u h g H − ∇u h , l H 0.

5.7
However, α is a user-selected stabilization parameter and typically α O h see 6, 7 for the details of the functional settings .In the numerical computations, the iterative scheme is adopted to implement these two subgrid models 6 .
There exists a fundamental requirement for subgrid models, that is, effects from subgrid should tend to vanish for laminar flows the suitable mesh resolutions.In order to implement the investigations, Reynolds number Re 400 is adopted to implement lid-driven cavity flows.The computational domain Ω 0, 1 × 0, 1 and the top boundary velocity u, v 1, 0 , as well as the other three boundaries are nonslip boundary conditions.The mesh scales are 1/49.Under this Reynolds number, the time-dependent Navier-Stokes equations will converge to stationary laminar solutions.The convergence tolerance is set to  equal 10 −8 .In Figure 3, the comparisons among four methods are given based on the profiles of horizontal and vertical velocity in mid planes of x and y axes.It is obvious that the three kinds of subgrid models give the almost same results compared with the Galerkin method.In Figure 4, we show the streamline patterns by the four different computational methods.From the streamline patterns, it is observed that the center position of the bottom-left vortex by the current subgrid model coincides with that by Galerkin method very well.However, GMQ and KLR models do not show the good coincidence.Numerically, the proposed subgrid model yields 0.3% deviation from the Galerkin solution for the center position of the bottomleft vortex, but GMQ and KLR models introduce about 2% ∼ 10% errors.It is known that the lid-driven flows of Re 400 is laminar and the "ideal" subgrid models should not affect flow structures.Under the sufficient mesh resolution, the two-grid subgrid models introduce the extra-action dissipation to flow fields according to this simple investigation, though three subgrid models can give the almost same results about the center positions of the center vortex and the bottom-right vortex.So, the other two subgrid models tend to be "unideal."

Fluid Flows around a Cylinder by High Reynolds Numbers
In this part, we investigate the two-dimensional under-resolved fluid flow around a cylinder.In this kind of fluid flows, the flow patterns are affected by the interaction of the fluid flows with two parallel planes and the surface of the cylinder.This problem is very useful to validate the subgrid model by vortex street patterns.The success and failure of the subgrid model simulations are useful for real fluid flows and engineering applications.
The geometry of fluid flows is shown in Figure 5.The height H and width W of the channel are equal to 1 and 6, respectively.The origin of the cylinder is at x, y 1, 0.   the radius R is equal to 0.15.The time-dependent inlet flow velocity profiles are given by and the boundary condition of the outlet is set as ∂u ∂n 0. 5.9 The boundary conditions of the two parallel planes and the cylinder surface are set as the nonslip boundary conditions.If the diameter D of the cylinder is chosen as the characteristic length and the mean velocity U mean of inlet as the characteristic velocity, the Reynolds number is defined by Re U mean D/ν.In order to implement the Galerkin method finite-element direct numerical simulation FEDNS under the current computer capability and limited hardware resources, the Reynolds number is set to equal 1000.By this Reynolds number, we give the comparison results among referred four different methods.The time interval Δt 0.001.The mesh scale h 1/50 for Galerkin method and h 1/16 for the three kinds of subgrid models.The mesh scale of Galerkin method is enough to resolve the small scale of current flow structures h ∝ Re −1/2 2 .From Figures 6 and 7, it is very clear that the proposed subgrid model in this paper predicts the flow structures best, compared with the other two subgrid models based on the Galerkin benchmark solutions.From these results, it is proved that the GMQ and KLR subgrid models introduce a too strong local dissipation into the flow structures, but the current subgrid model presents a suitable local dissipation behavior.
In order to assess performance of the current subgrid model, the very high Reynolds numbers Re 10 5 and Re 10 6 are chosen to implement the computations of complex flow phenomena.In Figures 8 and 9, the snapshots of vortex shedding are given.It is obvious that under these very high Reynolds numbers, the flow structures are complex and develop into the two-dimensional turbulence.The vortex filaments are clearly visible.These results demonstrate that the proposed subgrid model is effective and flexible.

Flow around an NACA0012 Airfoil at Zero Incidence
In this part, the proposed subgrid model will be used to simulate the flow around the  models.For the sake of brevity, we only show the comparisons between the proposed subgrid model and literature results 5 .By this, it is demonstrated that the proposed model can remove the spurious oscillations reasonably.

Conclusion
In this paper, a subgrid model is proposed for the time-dependent Navier-Stokes equations.
The corresponding error estimates are given.The numerical investigations validate that the proposed model is effective and flexible.At the same time, the proposed model only acts on the small-scale fluid flows and does not affect the large-scale flow structures.Numerical investigations also demonstrate that the proposed subgrid model is superior to recent proposed subgrid models and is effective for the simulations of very high Reynolds number fluid flows.This model is established based P 2 /P 1 finite-element spaces and the corresponding scale separation is achieved by L 2 -projection.Computationally, the implementation of the proposed model is very easy for some legacy codes of fluid flows and does not need the background coarse mesh.In future, this subgrid method will be attempted to implement simulations for 3D high Reynolds turbulence flows.

Lemma 4 . 1 .Figure 1 :
Figure 1: L 2 Convergence rate of u h by a log-log plot.

Figure 5 :
Figure 5: The schematic geometry of fluid flow past a cylinder.

Figure 6 :
Figure 6: The contour patterns of pressure fields t 1 .

Figure 7 : 1 . a t 6 b t 7 Figure 8 :t 6 b t 7 Figure 9 :
Figure 7: The contour patterns of horizontal velocity fields t 1 .

Figure 10 :
Figure 10: Mesh near the edge of airfoil.

NACA0012 airfoil at zero incidence 5 .
Two different Reynolds numbers are used to simulate fluid flows: Re 10 5 and Re 10 6 .The velocity of the incoming flow at infinity is the reference velocity scale and the chord of the airfoil is the reference length scale 5 .The flow domain is −5, 6 × −5, 5 and the upstream velocity is enforced on {x −5; −5 ≤ y ≤ 5} ∪ {−5 ≤ x ≤ 6; y ±5}.Natural boundary conditions are set at the downstream boundary.The mesh scale on the surface of airfoil is h 1/120.The mesh includes 45894 triangle cells.The mesh partition is given in Figure 10.The time interval Δt 10 −3 is employed for implementing computations of Re 10 5 and Re 10 6 .The effects of the proposed subgrid term are shown in Figures 11 and 12

Table 1 :
L 2 and H 1 relative errors and absolute errors of numerical velocity u h .