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 -projection. The fluctuation operator can be implemented by the -projection from high-order interpolation finite-element spaces to the low-order interpolation finite-element spaces. In this paper, 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.

1. 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-the-art 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 scales. The eddy viscosity models are often utilized to solve the convection-dominating or high Reynolds number NSE by engineers, which have been achieved many successes in engineering applications [1]. These kinds of models are firstly proposed by Frisch and Orszag [2], developed by Iliescu and Layton [3], and introduced a dissipative mechanism by Smagorinsky [4]. At present, these models have been further improved by various numerical methods [57]. In some recent models, the scale-separation subgrid terms are constructed by two different finite-element spaces based on a two-hierarchy mesh structure. Recently, Hughes et al. have proposed a variational multscale method (VMM) in which the diffusion acts only on the finest resolved scales. Generally, there exist different ways to define coarse and small scales under VMM frameworks [8]. According to the idea of VMM, the referred subgrid methods are variational multiscale methods.

In this paper, we will implement a novel way to circumvent the dispersive effects from small resolve scales. The subgrid term is established in a complement space of a low-order interpolation finite-element space in a high-order interpolation finite-element space. The complement space is established by an -projection decomposition of flow strain tensors. This method is easy to be implemented on a one-level mesh. This subgrid term does not act on the large-scale flow structures. For laminar low-Reynolds number flows, the action of this subgrid term tends to vanish. Numerical investigations demonstrate that the proposed model is effective and flexible for fluid flows of high Reynolds numbers ( and ).

The outline of the paper is organized as follows. In the following section, we introduce the nonstationary incompressible Navier-Stokes equations and the corresponding function settings. In Section 3, the subgrid viscous term is introduced into NSEs and the standard Galerkin discretization of the Navier-Stokes problem is given. In Section 4, we show the results of the error estimates. Some numerical results are presented in Section 5, which show the correctness and efficiency of the methods. Finally, we give some conclusions.

Let be a polygonal domain with Lipschitz continuous boundary . We consider the time-dependent Navier-Stokes equations for incompressible flow as follows: where a finite-time interval, represents the velocity vector, is the pressure, is the body force, and is the viscosity.

We introduce the following functional settings: We denote by and the inner product and norm in or . The space or denotes the standard Sobolev spaces within norm and seminorm . The space or is equipped with the following scalar product and norm: For convenience, we introduce the following bilinear form on : and on defined by The trilinear term is defined by which is the skew-symmetric form of the convective term. It is easy to gain Also, we have the following estimates [9]: where 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 such that and .

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 and that (2.9) possesses a solution with

Note. These assumptions imply that the solution of (2.9) is unique. For simplicity, let . In addition, we assume that has a polygonal boundary such that no boundary approximation in the application of the finite-element method becomes necessary.

3. Discretization of Navier-Stokes Equations and Subgrid Model

We give a family , which is a partition of into triangles or quadrilaterals, assumed to be regular in the usual sense [11]. The diameter of the cell is denoted by . The mesh parameter describes the maximum diameter of the cells .

We introduce the finite-dimensional subspace and : We define the space of discretely divergence free functions denoted by : Let , be arbitrary, let be a projection of , and let the following approximation properties hold [12]: here, , the constants depend only on . The regularity assumptions (2.11) imply that Meanwhile, the velocity-pressure pair in satisfies the following discrete inf-sup condition [13]:

Remark 3.1 ([14, 15]). Let be the standard -projection with the following properties: where . For a tensor , 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, ( is a real number). The analogous stabilization is used to circumvent the pressure stabilization LBB condition for Stokes problems [16]. Since is an -projection, it follows for and that In addition, from , it follows that Note. If depends on , then too. From (3.9) it follows that if is bounded almost everywhere in the time interval. If , then since . In this case, we set . The analogous formula was proposed to define a reduced Reynolds number for a variational multiscale method of the Navier-Stokes equations [7].

Using the stabilization term , we give the following stabilization finite-element discretization form of the variational form (2.9).

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

Find such that

4. 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 and , that is, certain norms of and are bounded by the data of the problem: .(2)Derive an error equation by subtracting (3.11) from (2.9) for test functions from . Split the error into an approximation term and a remainder : where is a projection of which fulfills the estimates (3.3). Then take 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 .(4)Apply Gronwall's lemma to (4.2), that is, derive all the functions in (4.2) belong to and get an estimate of .(5)Derive the error estimate of by applying the triangle inequality to (4.1).

Along these lines, the estimate will be proved [7]. This error estimate uses the parameter 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 and .

Lemma 4.1. The solution of the finite-element problem (3.10) fulfills and . The velocity solution of the continuous problem (2.9) fulfills and .

Proof. The proof for the and is very similar. We will show the result for . Setting in (3.11), use the skew-symmetric form of , (3.8), and integrate over with , we get Here is the value of at .

Subtraction of the last term and the regularity (2.10) gives . Taking the supremum of gives the statement .

Now, Step (2) of the proof is carried out by a strightforward calculation. We obtain with arbitrary .

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 The trilinear term is decomposed into three terms as follows: Using (2.8) and Young's inequality, we have Collecting all the above terms, we obtain We define It follows that is smaller or equal than Re. Using (3.9) finishes Step (3) of the proof: To perform Step (4) of the proof, we have to study the -regularity of the terms in (4.10). Let be arbitrary. Using Poincaré's inequality, the Cauchy-Schwarz inequality, then (2.11) and (3.4), we get Similarly by Hölders inequality, Lemma 4.1, and (3.4), we have The -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 be the solution of (2.9) and let be the solution of (3.11). Let the regularity assumptions (2.11) be fulfilled, let be a projection of into such that fulfills (3.3). Let the reduced Reynolds number be defined in (4.9). Then, the error satisfies, for ,

Corollary 4.3. From Theorem 4.2, the approximation result (3.3) and Remark 3.1, one has In particular,

Remark 4.4. A finite-element error estimate for the -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 . 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.

5. Numerical Tests

Firstly, we give the algorithm used to deal with the nonlinear term and the subgrid eddy viscosity term. Given , we find satisfying where and is the time interval. Also, is the finite-element approximate solutions of . The fixed point iteration is adopted to solve (5.1) based on the Newtonian method [17] and the iteration tolerance is set to equal .

For calculating the subgrid term , the subgrid term is rewritten as follows: In order to reduce the computational work, the following linear time extrapolation method is used to calculate :

5.1. 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 , with a body force obtained such that the following true solution is given by : The viscosity is and the corresponding Reynolds number is . The mesh scales we choose are (). Let and , where denotes the averaging mesh scale. The time step and the final time . In Table 1, the relative errors and absolute errors are given by different mesh scale . In Figures 1 and 2, and convergence orders are given by two log-log plots. From these two figures, and convergence orders are equal to and , respectively. The calculated convergence orders coincide with the theoretical analysis in Corollary 4.3. The expected convergence orders for the velocity in and 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.

5.2. Lid-Driven Cavity Flows

In this part, the comparisons among three kinds of subgrid models (Current model, Guermond-Marra-Quartapelle (GMQ) model [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 is the measure of , , and (see [5] for the details of the functional settings and the definition of the operator ). The KLR model is also based on two-level finite-element space, which is defined by [6, 7] where is the projection of However, is a user-selected stabilization parameter and typically (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 under the suitable mesh resolutions. In order to implement the investigations, Reynolds number is adopted to implement lid-driven cavity flows. The computational domain and the top boundary velocity , as well as the other three boundaries are nonslip boundary conditions. The mesh scales are . Under this Reynolds number, the time-dependent Navier-Stokes equations will converge to stationary laminar solutions. The convergence tolerance is set to equal . In Figure 3, the comparisons among four methods are given based on the profiles of horizontal and vertical velocity in mid planes of and 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 bottom-left vortex, but GMQ and KLR models introduce about errors. It is known that the lid-driven flows of 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."

5.3. 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 and width of the channel are equal to and , respectively. The origin of the cylinder is at and the radius is equal to . The time-dependent inlet flow velocity profiles are given by and the boundary condition of the outlet is set as The boundary conditions of the two parallel planes and the cylinder surface are set as the nonslip boundary conditions. If the diameter of the cylinder is chosen as the characteristic length and the mean velocity of inlet as the characteristic velocity, the Reynolds number is defined by . 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 . By this Reynolds number, we give the comparison results among referred four different methods. The time interval . The mesh scale for Galerkin method and for the three kinds of subgrid models. The mesh scale of Galerkin method is enough to resolve the small scale of current flow structures () [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 ( and ) 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.

5.4. Flow around an NACA0012 Airfoil at Zero Incidence

In this part, the proposed subgrid model will be used to simulate the flow around the NACA0012 airfoil at zero incidence [5]. Two different Reynolds numbers are used to simulate fluid flows: and . 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 and the upstream velocity is enforced on . Natural boundary conditions are set at the downstream boundary. The mesh scale on the surface of airfoil is . The mesh includes 45894 triangle cells. The mesh partition is given in Figure 10. The time interval is employed for implementing computations of and .

The effects of the proposed subgrid term are shown in Figures 11 and 12. The proposed subgrid model suppresses the spurious oscillations plaguing the Galerkin solution [5]. The proposed subgrid model possesses the same function as the other two referred subgrid 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.

6. 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 finite-element spaces and the corresponding scale separation is achieved by -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.


The authors thank the reviewers for their suggestions and critical comments which helped us to improve this paper. This work is supported by Chinese NSF (Grant no. 10671154) and the National Basic Research Program (no. 2005CB321703).