Modified Method of Characteristics Variational Multiscale Finite Element Method for Time Dependent Navier-Stokes Problems

AbstractIn this paper, a modified method of characteristics variational multiscale (MMOCVMS) finite element method is presented for the time dependent Navier-Stokes problems, which is leaded by combining the characteristics time discretization with the variational multiscale (VMS) finite element method in space. The theoretical analysis shows that this method has a good convergence property. In order to show the efficiency of the MMOCVMS finite element method, some numerical results of analytical solution problems are presented. First, we give some numerical results of lid-driven cavity flow with Re=5000 and 7500 as the time is sufficient long. From the numerical results, we can see that the steady state numerical solutions of the time-dependent Navier-Stokes equations are obtained. Then, we choose Re=10000, and we find that the steady state numerical solution is not stable from t=200 to 300. Moreover, we also investigate numerically the flow around a cylinder problems. The numerical results show that our...


Introduction
In this paper, we consider the time-dependent Navier-Stokes problems x ∈ Ω, u(x, t) = 0, x ∈ ∂Ω × [0, T ], (1.1) where Ω is a bounded domain in R 2 assumed to have a Lipschitz continuous boundary ∂Ω. u = (u 1 (x, t), u 2 (x, t)) T represents the velocity vector, p(x, t) represents the pressure, f (x, t) is the body force, ν = 1/Re the viscosity number, Re is the Reynolds number. Developing efficient finite element methods for the Stokes and Navier-Stokes equations is a key component in the incompressible flow simulation. There are some iterative methods for solving the stationary Navier-Stokes equations under some strong uniqueness conditions were presented by some authors [15,18,19]. The VMS method was first introduced by Hughes and his coworkers in [22,23]. The basic idea is splitting the solution into resolved and unresolved scale, representing the unresolved scales in terms of the resolved scales, and using this representation in the variational equation for the resolved scales. And there are many works devoted to this method, e.g. VMS method for the Navier-Stokes equations [26]; a two-level VMS method for convectiondominated convection-diffusion problems [29]; VMS methods for turbulent flow [9,30,31]; large-eddy simulation (LES) [24,25,32]; subgrid-scale models for the incompressible flow [22,44]. In [27,28], John et al. presented the error analysis of the two different kinds of VMS for the Navier-Stokes equations. The main difference is the definition of the large scales projection (L 2 -projection in [28] and elliptic projection in [27]). There is another class of VMS method which rely on a three-scale decomposition of the flow field into large, resolved small and unresolved scales [8].
In the characteristics method, which is a highly effective method for advection dominated problems, the hyperbolic part (the temporal and advection term) is treated by a characteristic tracking scheme. In 1982, Douglas and Russell [11] presented the modified method of characteristics finite element method (MMOCFEM) firstly. Russell [36] extended it to nonlinear coupled systems in two and three spatial dimensions. The second order in time method for linear convection diffusion problems had been given by Ewing and Russell [13]. A characteristics mixed finite element method for advection-dominated transport problems was presented by Arbogast [2] and for Navier-stokes equations by Buscagkia and Dari [7], respectively. In [34], a detail theory analysis for Navier-Stokes equations have been done by Pironneau, who obtained suboptimal convergence rates of the form O(h m + t + h m+1 / t) and improved by Dawson et al [10]. In [4], Boukir et al. presented a second-order time scheme based on the characteristics method and spatial discretization of finite element type for the incompressible Navier-Stokes equations. An optimal error estimate for the Lagrange-Galerkin mixed finite element approximation of the Navier-Stokes equations was given by Süli in [43]. In [38,39], one order and second order MMOC mixed defect-correction finite element methods for timedependent Navier-Stokes problems were given. In [41], second order in time MMOC variational multiscale finite element method for the time-dependent Navier-Stokes equation was shown. In [40], we presented modified characteristics gauge-Uzawa finite element method for the conduction-convection equation. In [42], the modified characteristics finite element method for the Navier-Stokes/Darcy problem was shown.
In this paper, we present a modified method of characteristics VMS finite element method based on the L 2 projection for the time dependent Navier-Stokes equations. In our method, the hyperbolic part (the temporal and advection term) is treated by a characteristic tracking scheme. The VMS based on projection finite element method is used for spatial discretization. The error analysis shows that this method has a good convergence property. In order to show the efficiency of the MMOCVMS finite element method, we firstly present some numerical results of analytical solution problems. The numerical results show that the convergence rates are O(h 3 ) in the L 2 -norm for u, O(h 2 ) in the semi H 1 -norm for u and O(h 2 ) in the L 2 -norm for p by using the Taylor-Hood element, which agrees very well with our theoretical results. Then, some numerical results of the lid-driven cavity flow with Re = 5000 and 7500 are given. We present the numerical results as the time are sufficient long enough, then the solution of the Navier-Stokes should approximate the solution of the steady state cases. From the numerical results, we can see that a steady state numerical solutions of the time-dependent Navier-Stokes equations are obtained. And the numerical solutions are in good agreement with that of the steady Navier-Stokes equations shown by Ghia et al. [14] and Erturk et al. [12]. Finally, we present some numerical results for Re = 10000. It shows that the solution is mainly periodic with small variations in the amplitude of the time evolution at the monitoring points as the time is sufficient long. And the phase portraits of the monitoring points show that the variations in amplitude yield a solution which is quasi-periodic. Furthermore, we present numerical simulations for the two-dimensional fluid flow around a cylinder. It is observed from these numerical results that the scheme can result in good accuracy, which shows that this method is highly efficient.
The rest of this paper is organized as follows. In Section 2, the functional settings of the Navier-Stokes equations are presented. In Section 3 the MMOCVMS finite element methods is proposed. The following Section 4 we derive optimal error estimate for the new algorithm. Numerical experiments will be carried out in Section 5 and some concluding remarks will be given in the final section.

Functional settings of the Navier-Stokes equations
In this section, we aim to describe some notations and results which will be frequently used in this paper. Firstly we introduce the standard Sobolev spaces Here Φ denotes a Hilbert space with the norm · Φ . Especially, when n = 0 we note And, we define For the mathematical setting of the Navier-Stokes problems (1.1), we introduce the following Hilbert spaces The following assumptions and results are recalled (see [18]). (A 1 ) There exists a positive constant C 0 which only depends on Ω, such that u 0 ≤ C 0 ∇u 0 ∀u ∈ X.
(A 2 ) Assume that Ω is smooth, hence the unique solution (v, q) ∈ (X, M ) of the steady Stokes problem for any prescribed g ∈ L 2 (Ω) 2 exits and satisfies where c > 0 is a generic constant depending on Ω, which may stand for different values at its different occurrences. · 2 means the H 2 −norm, · 1 means the H 1 −norm, and · 0 means the L 2 −norm, respectively.

MMOC-mixed finite element methods for the time-dependent Navier-Stokes equations
For each positive integer N , let {J n : 1 ≤ n ≤ N } be a partition of [0, T ] into subintervals J n = (t n−1 , t n ], with t n = n t, t = T /N . Set u n = u(·, t n ) be the u(x, t) at t = t n . The characteristic trace-back along the field u n−1 of a point x ∈ Ω at time t n to t n−1 is approximatelyx(x, t n−1 ) = x − tu n−1 .
Consequently, the hyperbolic part in the first equation of (1.1) at time t n is approximated by for any function w.
Let h be a quasi-uniform partition ofΩ into non-overlapping triangles, indexed by a parameter h = max K∈ h {h K ; h K = diam(K)}. We introduce the finite element subspace X h ⊂ X, M h ⊂ M as follows where P (K) is the space of polynomials of degree on K. ≥ 1, k ≥ 1 are two integers. Assume that (X h , M h ) satisfies the discrete LBB condition where d(ϕ, v) = (ϕ, ∇ · v), β > 0 is a constant independent of h. In this paper, V h denotes the kernel of the discrete divergence operator With the previous notations, we present the MMOC finite element method for the Navier-Stokes (1.1): find (u n , p n ) : otherwise.

The VMS finite element method and MMOCVMS finite element method for the time-dependent Navier-Stokes equations
For the VMS finite element method, a large scale space L H ⊂ L = {L ∈ (L 2 (Ω)) d×d , L = L T } and a so-called turbulent viscosity ν T ≥ 0 are introduced. The semi-discrete problems reads as follows (please see [31]): By 0 ≤ P L H ∇v 0 ≤ ∇v 0 , we get 0 ≤ ν add ≤ ν T . By straightforward calculation, we deduce Hence, system (3.2) can be reformulated as: The MMOC time discretization, combined with the VMS finite element method in space, leads to the following MMOCVMS finite element method. otherwise.

Remark 1.
Since L H has been distinguish between resolved small scales and large scales, with L H representing the large scales, it must be in some sense a coarse finite element space. One way of achieving this is choosing it to be a lower order finite element space than V h on the same grid (e.g. [27] for details). The other way is defining L H on a coarser grid. This two methods are called one-level and two-level projection-based FEVMS methods respectively (see [29,31]). The so-called turbulent viscosity ν T can be chosen in many ways. In this paper we choose as follows (Smagorinsky-type [26] then we have the following properties

Error estimate
In order to obtain the error analysis, we give some lemmas firstly.
. Then, we define the following Galerkin projection

Error estimate for the velocity
2Ln on the time step it is an easy matter to verify that this mapping has a positive Jacobian, since u l h vanishes on ∂Ω; this mapping is one-to-one and this it is a change of variables from Ω onto Ω. This yields for any positive function φ on Ω the estimate (please see [4] for details)

Remark 4.
Here, we assume that the initial condition u 0 has a strong regularity.
Or there are some restriction on the time step t [17].
Proof. We prove this lemma by induction. By the definition ofx, we can see that (4.5) and (4.4) hold for n = 0. We assume that (4.5) and (4.4) hold for 1 ≤ n ≤ l − 1, then L n < +∞. Now we will prove them hold for n = l.
Letting ϕ h = 0 in (3.6) and using the definition of V h , we have Define η l = u l − R l h , we arrive at where Π i denotes the i-th term of the right-hand side of (4.8).
Next we give the estimate for each term. By the definition ofx andx, we arrive atx Using Taylor formula, we obtain Therefore, we get Now, we estimate the bound of η l −η l−1 t 0 . By Cauchy-Schwarz inequality, we get By the definition ofX l x (t l−1 ) , we have J(X l x (t l−1 )) = Hence, det J(X l x (t l−1 )) = 1 + O( t). Then, we deduce Then, we deduce By (4.11) and (4.13), we get . (4.14) Similarly, we obtain Setting v h = ξ l h in (4.9), we obtain By (4.1), (4.10), (4.14) and (4.15), we get Summing inequality (4.16) from i = 1 to l, we obtain Using discrete Gronwall's Lemma, we have By triangle inequality, we deduce Using (4.17) and Lemma 4.4 we finish the proof. Theorem 1. Let u n h be defined by (3.1) and u be the solution of (1.1). Under the assumptions of Lemma 4, for all 1 ≤ n ≤ N we have Proof. Using triangle inequality, (4.3) and (4.4), we obtain the desired result.

Error estimate for pressure
The following result on the pressure is a consequence of the previous theorem on the velocity.
where C is a positive constant independent of t and h.
Proof. Multiplying (1.1) by v h ∈ X h and subtracting the result form (3.6) with ϕ h = 0, we obtain By the inf-sup condition and Cauchy-Schwarz inequality, we get Using (4.1), (4.10) and (4.18), we arrive at By triangle inequality, we have Therefore, we finish the proof.

Numerical Results
In this section, we present some numerical results to verify the theoretical results obtained in the previous section and show the effect of our methods.
Here, we use the software package FreeFEM++ [21] for our program.

Analytical solution problems
In this subsection, we present some numerical results of the Navier-Stokes problems with the analytical solution The boundary and initial conditions in (1.1) are set equal to the analytical solution and f is given by evaluating the momentum equation of the problem (1.1) for the analytical solution. We choose t = 0.001, T = 1, δ = 0.5 and Re = 2000. We present the numerical results with different h respectively by using Taylor-Hood element, and L H = {L ∈ (L 2 (Ω)) 2×2 , L| K ∈ P 2×2 1dc ∀K ∈ h }, P 1dc means the piecewise linear discontinuous finite element space [16].     The numerical results in Table 1 and 2 show that the convergence rates are O(h 3 ) of the L 2 -norm for u, O(h 2 ) of the semi H 1 -norm for u and O(h 2 ) of the L 2 -norm for p, which agree very well with our theoretical results by using P 2 − P 1 finite element spaces. In Table 3 and 4, we present the convergence test results for time. We can see that the time convergence order is O( t), which agrees well with our theory analysis. Comparing our numerical results with the MMOCFEM's, we can see that the error and the convergence rates are very similar. But our method can save much time.

The lid driven cavity problem
In this subsection we show the numerical results of lid driven cavity problem. The two-dimensional lid driven was formulated as in Ω = (0, 1) 2 , the boundary conditions are u 1 = 1, u 2 = 0 on the top lid and u 1 = 0, u 2 = 0 on the other lids. In a former work [5], it was suggested that the first Hopf bifurcation occurs around Reynolds number Re = 7500. Since then various results were given in the literature [5,12,14]. In this numerical experiments, h = 1 60 , t = 0.01, ν T = δh 2 ∇u h 0 , δ = 0.5 are chosen.We choose the Taylor-Hood element and L H = {L ∈ (L 2 (Ω)) 2×2 , L| K ∈ P 2×2 1dc ∀K ∈ h }. Firstly, we choose Re = 5000, which is a good choice as there are some comparisons available in the literature and as the steady solution is still stable but not too far from the first Hopf bifurcation. Figure 1           cavity agree very well.
Finally, we show some numerical results for Re = 10000, which is probably the most famous value and for quite a long time the question was to known if the steady solution was stable or not for this Reynolds number. In [6], Bruneau and Saad shown that the steady is not stable and the first Hopf bifurcation occurs around Re = 8000. Figure 3 presents the streamlines at different time. We choose two monitoring points (1/8, 13/16) and (7/8, 13/16). Figures 4  and 5 show u-velocity history (a), v-velocity history (b) and the phase (c) portrait at monitoring points (1/8, 13/16) and (7/8, 13/16). It shows that the stable solution in mainly periodic with small variations in the amplitude of the time evolution at the monitoring points. And the phase portraits show that the variations in amplitude yield a solution which is quasi-periodic. Figure 6 shows the evolution of the kinetic energy in time (a) and evolution of u n h − u n−1 h 0 / u n h 0 in time (b). We can see that the kinetic energy dose not change as the time change long enough and the errors don't change smaller. The results are very close the those shown in [6].

The flow around a cylinder problems
In this part, we investigate the two-dimensional under-resolved fluid flow around a cylinder [37], the physical model given by Figure 7. From Figure 7 we can see that, the height of the channel H = 0.41 and the width W = 2.2. The origin of the cylinder is at (0.2, 0.2), and the radius R is equal to 0.05. The time-dependent inlet flow velocity profiles are given by The boundary conditions of the two parallel planes and the cylinder surface are set as the non-slip boundary conditions. The Reynolds number is defined by Re = 2Ū R/ν andŪ = 2u(0, H/2)/3. Here, we choose Re = 100, t = 0.01. The finite element grid is presented by Figure 7. We choose the Taylor-Hood element spaces. Figures 8 and 9 give the numerical results by the MMOCVMSFEM. Figures 10 and 11

Conclusions
In this work, we proposed a modified algorithm of characteristics variational multiscale finite element method for the time dependent Navier-Stokes equations. The theoretical analysis is discussed and illustrated. Error estimate shows that this method has a good convergence property. Numerical experiments are presented to demonstrate the high efficiency of the new algorithm. Also, further developments can extend these techniques and ideas to the general nonlinear problems and high-dimensional problems.