The Upwind Finite Volume Element Method for Two-Dimensional Burgers Equation

and Applied Analysis 3 where P l (V) (l = 0, 1) denotes the set of polynomials on V with a degree of not more than l. Multiplying (6a) by test function z ∈ V h and integrating by parts yield ( ∂θ ∂t , z) + B 1 (θ; θ, z) + B 2 (θ; θ, z) + A (θ, z) = 0, ∀z ∈ V h , (9)


Introduction
We consider the following two-dimensional Burgers equation [1][2][3]: for the unknown functions  and  in a bounded spatial domain Ω ⊂ R 2 , over a time interval [0, ].The coefficient  is a positive number.Burgers equation is the simplest nonlinear convectiondiffusion model [1].It is often used in modeling such physical phenomena as turbulence, shocks, and so forth.The study of Burgers equation has been a very active area because of its importance.
It is well known that strictly parabolic discretization schemes applied to Burgers equation do not work well when it is advection dominated.Effective discretization schemes recognize to some extent the hyperbolic nature of the equation.
The finite volume element method (FVEM) [4][5][6][7][8][9][10][11][12] is an important discretization technique for partial differential equations, especially those that arise from physical conservation laws.FVEM has ability to be faithful to the physics in general and conservation in particular, to produce simple stencils, and to treat effectively Neumann boundary conditions and nonuniform grids, and so forth.
Liang [11,12] combined the upwind technique and the FVEM to handle the linear convection-dominated problems.In this paper, we will consider upwind finite volume element method for the approximation of (1).Upwind approximation is applied to handle the nonlinear convection term.The semidiscrete and fully discrete schemes are defined, respectively.We prove that they are both convergent to order one in space.Numerical experiments are presented finally to validate the theoretical analysis.
An outline of the paper follows.In the next section we define the upwind finite volume element schemes for (1).Some lemmas are presented in Section 3. We derive the  2norm error estimates for the semi-discrete scheme and the fully discrete scheme in Sections 4 and 5, respectively.Finally in Section 6, we give some numerical experiments.
Throughout the paper we will denote by  and   ( = 1, 2, . ..) generic constants independent of the mesh parameters, which may take different values in different occurrences.

The Approximation Schemes
In order to rewrite (1) as the vector form we define some vector notations.The gradient of a vector function  = ( 1 ,  2 ) : R 2 → R 2 is a matrix, and the divergence of a matrix function  = (  ) 1≤,≤2 : R 2 → R 2 × 2 is a vector Consequently, we have for a vector function  = ( 1 ,  2 ) Let  = (, ),  0 () = ( 0 (),  0 ()), and let  = ( 1 ,  2 ); then the system (1) can be written as the following vector form: where Let T ℎ = {} be a triangulation of the domain Ω, and as usual, we assume the triangles  to be shape regular.Denote by Ω ℎ = {  } the set of the vertices of all the triangles , and let Ω ℎ = Ω ℎ \ Ω.For a given triangulation T ℎ , we construct a dual mesh T * ℎ whose elements are called control volumes.Each triangle  ∈ T ℎ can be divided into three subdomains by connecting an inner point of the triangle to the midpoints of the three edges.Around each   ∈ Ω ℎ , we associate a control volume  *  =  *   , which consists of the union of subregions having   as a vertex.For a vertex   ∈ Ω, we can define its control volume in a similar way.Then we define the dual partition T * ℎ = { *   ,   ∈ Ω ℎ } to be the union of all the control volumes.Usually we can choose the inner points as the barycenters or the circum centers, and in the later case we assume that all the inner angles of each triangle are not larger than /2.We will use the barycenters duel mesh in this paper, while, with some trivial changes, our analysis can be also applied to the case when the circum centers are used.
We now characterize the finite-dimensional spaces which will be employed in approximating (6).For the sake of simplicity, we will assume that  1 =  2 = 0. We define the following finite dimensional spaces: where P  () ( = 0, 1) denotes the set of polynomials on  with a degree of not more than . where here  is the unit outward normal vector of  *  .Now we approximate  2 (; , ) by using the upwind technique.
Let Λ  = { :   is adjoint with   }.Assuming that  ∈ Λ  , let Γ  =  *  ∩  *  and   is the length of Γ  .Denote by   the unit outward normal vector of Γ  when Γ  is regarded as the boundary of  *  .Define The upwind discretization of the nonlinear term  2 (; , ) is defined by the form Using the heaviside function The semi-discrete upwind finite volume scheme of ( 6) is as follows: find  ℎ : [0, ] →  ℎ such that where  0ℎ () is the interpolation projection of  0 , that is, Our analysis is valid for variable time steps, but we drop the superscript from  for convenience.For functions  on Ω × , we write   () for (,   ).By approximating  ℎ / at the time  =   with the backward difference     ℎ = (  ℎ − −1 ℎ )/, we define the fully discrete upwind finite volume scheme for (6) as follows: find

Some Lemmas
Now we present several Lemmas.Let From [4] we know that for  1 ,  2 , w1 , w2 ∈  ℎ , where  1 and  2 are some positive constants that are independent of ℎ.Thus we obtain ( 19) and (20) immediately.
Lemma 2. For the bilinear form (⋅, Π * ℎ ⋅), one has the following conclusions: (ii) There exists a positive constant  such that (iii) There exists a positive constant  such that Proof.For ,  ∈  ℎ , define the bilinear form Then, we get By combining the above results and the corresponding conclusions for (⋅, Π * ℎ ⋅) in [4], we can obtain ( 23)-(25).
Proof.First we have Noting that (  ) = Π ℎ (  ),   ∈ Ω ℎ , we can easily deduce by the definition of  2ℎ (⋅; ⋅, Π * ℎ ⋅).Now we only need to bound We denote the last two terms on the right-hand side of (31) by  1 and  2 , respectively.We now turn to analyze the two terms.Noting that   = −  , we rewrite  1 as here Λ  is the set of vertex of .From the Taylor's Formula and the linear property of  ℎ = ( where  depends on principally , and Proof.We derive the following error equation from ( 6) and ( 17): Let  =  − Π ℎ ,  = Π ℎ  −  ℎ .We rewrite the previously mentioned equation as We choose  ℎ =  in (41) to get Using Lemmas 1, 2 and Young's inequality, we have Now we bound the last two terms on the right-hand side of (43).We need the following induction hypothesis: We know from [13] that where , and ‖‖  1 (( 2 ) 2 ) .We now complete the proof of the theorem.

Theorem 5. Assume that 𝜃 satisfies the necessary regularities and the discretization parameters obey the relation 𝜏 = 𝑂(ℎ).
Then the error of the approximation (18) of (6) satisfies where  depends on , and Proof.Subtract (18) from ( 9) to obtain that Choose  ℎ =   to obtain that For the left-hand side of (62), from Lemmas 1 and 2, we have Now we estimates the terms  1 , . . .,  5 one by one.From the Taylor's formula, we have  (67) We make the following induction hypothesis: For  4 , using the similar argument as (48) and noting (68), we deduce that         Now we prove the induction hypothesis (68).Noting that  0 ℎ = Π ℎ  0 , we know that  0 = 0. From (75) and the assumption  = (ℎ), we get that (log 1 ℎ )

Numerical Example
In this section, we will show the affectivity of our method by numerical experiments.The exact solutions to problem  where  = 1/.We present numerical results for the  2 -norm estimates of  −  ℎ and  −  ℎ .In Tables 1 and 2, we present the numerical results for  = 1 and  = 0.01, respectively.In all runs, we use the uniform mesh step ℎ = Δ and choose the time  = 1.As seen in these tables, in all cases the errors decrease by a factor of about two as ℎ decreases by the factor of two.This indicates that all  2 -norm error estimates are of first-order convergence, which is consistent with our theoretical analysis.When  = 1.0 and  = 0.01, the figures of the exact solutions  and the numerical solutions  ℎ at  = 1 for ℎ = 1/32 are given in Figures 1, 2, 3, and 4. In order to show that our method keeps stable when  is smaller, we also give the comparison figures of exact solution  and numerical solution  ℎ for  = 0.001 in Figures 5 and 6.The comparison figure of numerical solution by using finite volume element method (FVEM) without upwinding is given in Figure 7, which show that the approximation produces unacceptable nonphysical oscillations.