Research Article  Open Access
A Subgrid Model for the TimeDependent NavierStokes Equations
Abstract
We propose a stabilized subgrid finiteelement method for the twodimensional (2D) nonstationary incompressible NaverStokes 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 highorder interpolation finiteelement spaces to the loworder interpolation finiteelement spaces. In this paper, mixed finiteelement 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 NavierStokes 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 largescale fluid flows. Meanwhile, for laminal fluid flows, the added subgrid viscosity term should not affect the largescale 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 largescales to the small Kolmogorov scales which can hardly be resolved by stateoftheart 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 convectiondominating 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 [5–7]. In some recent models, the scaleseparation subgrid terms are constructed by two different finiteelement spaces based on a twohierarchy 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 loworder interpolation finiteelement space in a highorder interpolation finiteelement space. The complement space is established by an projection decomposition of flow strain tensors. This method is easy to be implemented on a onelevel mesh. This subgrid term does not act on the largescale flow structures. For laminar lowReynolds 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 NavierStokes equations and the corresponding function settings. In Section 3, the subgrid viscous term is introduced into NSEs and the standard Galerkin discretization of the NavierStokes 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.
2. NavierStokes Equations and Functional Settings
Let be a polygonal domain with Lipschitz continuous boundary . We consider the timedependent NavierStokes equations for incompressible flow as follows: where a finitetime 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 skewsymmetric 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 finiteelement analysis, we need some regularity assumptions of the NavierStokes 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 finiteelement method becomes necessary.
3. Discretization of NavierStokes 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 finitedimensional 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 velocitypressure pair in satisfies the following discrete infsup 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 finiteresolution 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 NavierStokes 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 userselected 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 NavierStokes equations [7].
Using the stabilization term , we give the following stabilization finiteelement discretization form of the variational form (2.9).
Find satisfying Under the infsup condition (3.5), formulation (3.10) is equivalent to the following problem.
Find such that
4. Error Analysis
The proof of the finiteelement 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 finiteelement 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 skewsymmetric 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 righthand side of (4.4) have to be estimated. Using the CauchySchwarz 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 CauchySchwarz 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 finiteelement error estimate for the error in the pressure can also be derived following Heywood and Rannacher [12]. Using the infsup 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 finiteelement 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 loglog 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 secondorder and the firstorder, 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. LidDriven Cavity Flows
In this part, the comparisons among three kinds of subgrid models (Current model, GuermondMarraQuartapelle (GMQ) model [5], and KayaLaytonRiviere (KLR) model [6, 7]) are investigated for liddriven cavity flows. These investigations will address the actions of subgrid models on largescale flow structures. Now, we introduce the GMQ model and KLR model. The GMQ model is based on the twolevel Lagrange finiteelement 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 twolevel finiteelement space, which is defined by [6, 7] where is the projection of However, is a userselected 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 liddriven 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 timedependent NavierStokes 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 bottomleft 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 errors. It is known that the liddriven flows of is laminar and the “ideal" subgrid models should not affect flow structures. Under the sufficient mesh resolution, the twogrid subgrid models introduce the extraaction (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 bottomright vortex. So, the other two subgrid models tend to be “unideal."
(a)
(b)
(a) Galerkin method
(b) Current model
(c) GMQ model
(d) KLR model
5.3. Fluid Flows around a Cylinder by High Reynolds Numbers
In this part, we investigate the twodimensional underresolved 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 timedependent 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 (finiteelement 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.
(a) Galerkin method (GM)
(b) Current model
(c) GMQ model
(d) KLR model.
(a) Galerkin method (GM)
(b) Current model
(c) GMQ model
(d) KLR model
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 twodimensional turbulence. The vortex filaments are clearly visible. These results demonstrate that the proposed subgrid model is effective and flexible.
(a)
(b)
(a)
(b)
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.
(a) Velocity
(b) Velocity
(c) Pressure
(d) Vorticity
(a) Velocity
(b) Velocity
(c) Pressure
(d) Vorticity
6. Conclusion
In this paper, a subgrid model is proposed for the timedependent NavierStokes 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 smallscale fluid flows and does not affect the largescale 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 finiteelement 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.
Acknowledgments
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).
References
 V. John, Large Eddy Simulation of Turbulent Incompressible Flows, vol. 34 of Lecture Notes in Computational Science and Engineering, Springer, Berlin, Germany, 2004. View at: Zentralblatt MATH  MathSciNet
 U. Frisch and S. A. Orszag, “Turbulence: challenges for theory and experiment,” Physics Today, vol. 43, no. 1, pp. 24–32, 1990. View at: Publisher Site  Google Scholar
 T. Iliescu and W. J. Layton, “Approximating the larger eddies in fluid motion. III. The Boussinesq model for turbulent fluctuations,” Analele Ştiinţifice ale Universitătii “Al. I. Cuza” din Iaşi. Serie Nouă. Matematică, vol. 44, no. 2, pp. 245–261, 1998. View at: Google Scholar  Zentralblatt MATH  MathSciNet
 J. Smagorinsky, “General circulation experiments with the primitive equation—I: the basic experiment,” Monthly Weather Review, vol. 91, pp. 99–164, 1963. View at: Publisher Site  Google Scholar
 J.L. Guermond, A. Marra, and L. Quartapelle, “Subgrid stabilized projection method for 2D unsteady flows at high Reynolds numbers,” Computer Methods in Applied Mechanics and Engineering, vol. 195, no. 44–47, pp. 5857–5876, 2006. View at: Publisher Site  Google Scholar  Zentralblatt MATH  MathSciNet
 S. Kaya, W. Layton, and B. Rivière, “Subgrid stabilized defect correction methods for the NavierStokes equations,” SIAM Journal on Numerical Analysis, vol. 44, no. 4, pp. 1639–1654, 2006. View at: Publisher Site  Google Scholar  Zentralblatt MATH  MathSciNet
 V. John and S. Kaya, “Finite element error analysis of a variational multiscale method for the NavierStokes equations,” Advances in Computational Mathematics, vol. 28, no. 1, pp. 43–61, 2008. View at: Publisher Site  Google Scholar  Zentralblatt MATH  MathSciNet
 T. J. R. Hughes, L. Mazzei, and K. E. Jansen, “Large eddy simulation and the variational multiscale method,” Computing and Visualization in Science, vol. 3, no. 12, pp. 47–59, 2000. View at: Publisher Site  Google Scholar  Zentralblatt MATH
 W. Layton and L. Tobiska, “A twolevel method with backtracking for the NavierStokes equations,” SIAM Journal on Numerical Analysis, vol. 35, no. 5, pp. 2035–2054, 1998. View at: Publisher Site  Google Scholar  Zentralblatt MATH  MathSciNet
 R. Temam, NavierStokes Equations, Theory and Numerical Analysis, NorthHolland, Amsterdam, The Netherlands, 1983.
 V. Girault and P.A. Raviart, Finite Element Methods for NavierStokes Equations, vol. 5 of Springer Series in Computational Mathematics, Springer, Berlin, Germany, 1986. View at: MathSciNet
 J. G. Heywood and R. Rannacher, “Finite element approximation of the nonstationary NavierStokes problem. III. Smoothing property and higher order error estimates for spatial discretization,” SIAM Journal on Numerical Analysis, vol. 25, no. 3, pp. 489–512, 1988. View at: Publisher Site  Google Scholar  Zentralblatt MATH  MathSciNet
 A. Quarteroni and A. Valli, Numerical Approximation of Partial Differential Equations, Springer, Berlin, Germany, 1997.
 J. G. Heywood and R. Rannacher, “Finite element approximation of the nonstationary NavierStokes problem. I. Regularity of solutions and secondorder error estimates for spatial discretization,” SIAM Journal on Numerical Analysis, vol. 19, no. 2, pp. 275–311, 1982. View at: Publisher Site  Google Scholar  Zentralblatt MATH  MathSciNet
 P. B. Bochev, C. R. Dohrmann, and M. D. Gunzburger, “Stabilization of loworder mixed finite elements for the Stokes equations,” SIAM Journal on Numerical Analysis, vol. 44, no. 1, pp. 82–101, 2006. View at: Publisher Site  Google Scholar  Zentralblatt MATH  MathSciNet
 S. Ganesan, G. Matthies, and L. Tobiska, “Local projection stabilization of equal order interpolation applied to the Stokes problem,” Mathematics of Computation, vol. 77, no. 264, pp. 2039–2060, 2008. View at: Publisher Site  Google Scholar  MathSciNet
 Y. He and A. Wang, “A simplified twolevel method for the steady NavierStokes equations,” Computer Methods in Applied Mechanics and Engineering, vol. 197, no. 1718, pp. 1568–1576, 2008. View at: Publisher Site  Google Scholar  MathSciNet
Copyright
Copyright © 2009 Yan Zhang and Yinnian He. This is an open access article distributed under the Creative Commons Attribution License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.