A VIRTUAL ELEMENT DISCRETIZATION FOR THE TIME DEPENDENT NAVIER–STOKES EQUATIONS IN STREAM-FUNCTION FORMULATION

. In this work, a new Virtual Element Method (VEM) of arbitrary order 𝑘 ≥ 2 for the time dependent Navier–Stokes equations in stream-function form is proposed and analyzed. Using suitable projection operators, the bilinear and trilinear terms are discretized by only using the proposed degrees of freedom associated with the virtual space. Under certain assumptions on the computational domain, error estimations are derived and shown that the method is optimally convergent in both space and time variables. Finally, to justify the theoretical analysis, four benchmark examples are examined numerically


Introduction
In this work, we study a Virtual Element Method (VEM) for a fourth order nonlinear problem arising in the numerical discretization of the Navier-Stokes problem. The VEM, introduced in [6], is a generalization of the Finite Element Method (FEM) which is characterized by the capability of dealing with very general polygonal/polyhedral meshes, and it also permits to construct in a straightforward way highly regular discrete spaces. Indeed, by avoiding the explicit construction of the local basis functions, the VEM can easily handle general polygons/polyhedrons without complex integrations on the element (see [7] for details on the coding aspects of the method). The VEM has been applied successfully for problems in fluid mechanics; see for instance [1,8,17,18,24,26,30,33,34,39,41], where Stokes, Brinkman, Stokes-Darcy and Navier-Stokes equations have been recently developed.
The Navier-Stokes system is a paradigm of fluid flow problems. Usually, the variables and denote the velocity and the pressure field, respectively. It is proved that if the body force f and the initial data 0 are smooth enough and the boundary of domain Ω is locally Lipschitz continuous, then the two dimensional non stationary Navier-Stokes problem has weak solution. In [38], Temam showed that ∈ ∞ (︀ 0, ; [ 2 (Ω)] 2 )︀ with the assumption that f and initial data 0 are suitably smooth. Since the model problem consists of nonlinear term, it is not straightforward to find analytical solution. Therefore, numerical approximation is Keywords and phrases. Virtual element method, Navier-Stokes equations, time dependent problem, stream-function.

Preliminaries and the model problem
Let Ω ⊂ R 2 be a simply connected polygonal domain with boundary Γ := Ω. We denote by 2 (Ω), the space of square integrable scalar functions with the standard inner product ( , ) 0,Ω := ∫︀ Ω . For each positive integer ∈ N, we define (Ω) [2], the Sobolev space with standard norm where is multi index and denotes th partial derivative of . Let denote the time variable taking values in the interval := (0, ], where is a given final time. Moreover, the function space 2 (0, ; (Ω)) consists of scalar functions such that for almost all ∈ [0, ], (·, ) ∈ (Ω) [15] with the norm, In addition, given any Hilbert space , we will denote by [ ] 2 the space of vectors functions with entries in (see [2]). Further, we define := d d , , ≤2 denotes the Hessian matrix of , and := ∇ · , where is outward normal vector. For second order tensor fields , : Ω → R 2×2 , we define scalar product : : R 2×2 × R 2×2 → R by , where and are the entries at ( , )-th position of and , respectively.
Then, equation (2.2) can be rewritten in the following form: find ( ) ∈ Z such that Now, we reformulate the above problem as follows: since Ω is a simply connected domain, a well known result states that a vector function ∈ Z if and only if there exists a scalar function ∈ 2 (Ω) (called stream-function) such that = curl ∈ H.

Virtual element method
Let { ℎ } ℎ be a sequence of decompositions of Ω into general polygonal elements . Let ℎ denote the diameter of the element and ℎ the maximum of the diameters of all the elements of the mesh , i.e., ℎ := max ∈ ℎ ℎ . In what follows, we denote by the number of vertices of , by a generic edge of ℎ and for all ∈ , we define a unit normal vector that points outside of and a unit tangent vector to . Further, we denote by ℎ the length of the edge . For each vertex , we associate a characteristic length ℎ which is the average of the diameter of the elements having as a vertex. For the theoretical analysis, we will make the following assumptions: there exists a real number > 0 such that, for every ℎ and every ∈ ℎ , Assumption 3.1.
(a) The ratio between the shortest edge and the diameter ℎ of is larger than ; (b) ∈ ℎ is star-shaped with respect to every point of a ball of radius ℎ .

Local and global virtual spaces
Now, for any subset ⊆ R 2 and non negative integer , we will denote by ( ) the space of polynomials of degree up to defined on . Then, for ( 1 , 2 ) ∈ N × N, we define the set of scale monomials as where ( , ) denotes centroid of . Then, we define ℳ ( ) := ∪ ≤ ℳ * ( ) as a basis of ( ). Analogously, we consider the set of the scaled monomials defined on each edge : where is the midpoint of . Then, for any ≥ 2 and for every polygon ∈ ℎ , we introduce the following preliminary local virtual space:︀ where := max{3, } and := − 1.
Next, for a given ℎ ∈̃︀ ℎ ( ), we introduce five sets D 1 − D 5 of linear operators from the local virtual spacẽ︀ ℎ ( ) into R.
-D 1 : contains linear operators evaluating ℎ at the vertices of ; -D 2 : contains linear operators evaluating ℎ ∇ ℎ ( ) for all vertices of , where 1 ≤ ≤ ; In order to construct the discrete scheme, we first observe that Then, we decompose the bilinear form (2.7) in the following element by element contribution: In what follows, we are going to build the discrete version of the local bilinear forms listed above. With this aim, we define the following projector operator Π ,Δ :̃︀ ℎ ( ) −→ ( ) ⊆̃︀ ℎ ( ) for each ℎ ∈̃︀ ℎ ( ), as the solution of the local problems (on each polygon ): wherê︁ (·) is defined as follows:̂︁ and , 1 ≤ ≤ , are the vertices of . The following result establishes that the projector Π ,Δ is computable using the output values of the sets Proof. For detail proof, we refer to [21].
On the other hand, we observe that ( ) ⊆ ℎ ( ) which will guarantee the good approximation properties for the space. Moreover, following similar arguments presented in [3,16] (see also [21]) we obtain that the sets of linear operators D 1 − D 5 constitutes a set of degrees of freedom for ℎ ( ). Now, we introduce the global virtual space by combining the local spaces ℎ ( ) and incorporating the homogeneous Dirichlet boundary conditions. For every decomposition ℎ of Ω into polygons , we define

Construction of bilinear forms and the force term
In order to build the discrete local and global forms, we observe that the particular condition appearing in the definition of the local virtual space ℎ ( ) will be useful to construct an 2 -projection which will be employed to build the discrete bilinear forms. In particular, we consider the 2 ( )-projection onto −2 ( ). For each The following lemma establishes that Π −2 is computable on ℎ ( ). The proof follows from the definition of the local virtual space and the set of degrees of freedom. Proof. For a detail proof, we refer to [21]. Now, for ≥ 2, we will introduce some additional projectors which will be used to write the virtual scheme. First, we define Π ,∇ ⊥ : ℎ ( ) −→ ( ) ⊆ ℎ ( ) for each ℎ ∈ ℎ ( ) as the solution of the following local problem.
wherê︁ (·) has been defined in (3.2). The following result states that this operator is fully computable.
using only the information of the set of degrees freedom D 1 − D 5 . Proof. First we note that (3.5b) is computable using the information of the set D 1 . On the other hand, we integrate by parts on the right hand side of (3.5a) to obtain: where we have used the fact that ∆ ∈ −2 ( ) and the definition of the projection Π −2 (cf. (3.4)). The previous equality allows us to conclude that the polynomial Π ,∇ ⊥ ℎ can be explicitly computed from the degrees of freedom D 1 − D 5 . Now, we will introduce an additional projection operator onto the polynomial space [ −1 ( )] 2 , which will be used to construct a local approximation of (·, ·) and (·; ·, ·). For ∈ ℎ , and ∈ [ 2 ( )] 2 , we define We observe that for any ℎ ∈ ℎ ( ), the vector functions Π −1 curl ℎ ∈ [ −1 ( )] 2 and Π −1 ∇ ℎ ∈ [ −1 ( )] 2 can be explicitly computed from the degrees of freedom D 1 − D 5 . In fact, for all ∈ ℎ and for all ℎ ∈ ℎ ( ), using integration by parts on the right-hand side of (3.6) (with curl ℎ instead of ), we have (see [28]) where we have used (3.4). The first term on the right-hand side above depends only on Π −2 ℎ and this depends only on the values of the degrees of freedom (see Lem. 3.4). The second term is an integral on the boundary of the element , which is fully computable.
Next, we use the above projection operators to construct computable approximations of the continuous bilinear and trilinear forms, and for the right-hand side. First, let Δ (·, ·) and curl (·, ·) be any symmetric positive definite bilinear forms to be chosen as to satisfy: with 0 , 1 , 2 and 3 are positive constants independent of ℎ and . From (3.7), we deduce that Δ (·, ·), and curl (·, ·) scale same as (·, ·) and (·, ·), respectively. On each element , we define the local discrete bilinear forms It can be observed that the forms curl (·, ·) and Δ (·, ·) reduce to zero when one of the two entries ℎ or ℎ is a polynomial function. Different computable form of the stabilizers are available in the literature [5,16,36]. However, we choose the following representation where curl and Δ are the suitable constants, and dof denotes the number of degrees freedom of ℎ ( ). Adding the local contribution, the global forms are defined as The following result establishes the usual consistency and stability properties for the discrete local bilinear forms.
Lemma 3.7. There exists a constant > 0, independent of ℎ, such that Now, we proceed to discretize the force function as follows Globally, the force function ℎ is defined as follows

Discretization of the nonlinear term
In this section, we would like to discretize the trilinear term associated with the problem (2.3). Using the projection operators Π ,Δ and Π −1 , we discretize the trilinear term.
The term is fully computable from the degrees of freedom. Globally, the trilinear term is defined as Moreover, it can be shown that the discrete trilinear form ℎ (·; ·, ·) is uniformly bounded on ℎ .
Proof. An application of boundedness of the projection operators Π ,Δ and Π −1 , Hölder inequality, and Sobolev's embedding theorem yields the proof.
Remark 3.9. The discrete trilinear form ℎ (·; ·, ·) does not contain non-polynomial part or stabilizer. It is defined using the projection operators Π ,Δ and Π −1 which are computable from the information provided by the degrees of freedom. With this definition, we will show that the semi-discrete and fully-discrete schemes are well-posed and we will obtain the corresponding error estimates.

Discrete schemes and their well posedness
In this section, we will introduce the semi-discrete and fully-discrete virtual element schemes for problem (2.5), by using the discrete forms introduced in Sections 3.2 and 3.3. We will also prove that under some assumptions on , the fully-discrete scheme is well posed.
First, we observe that the matrix representation corresponding to the discrete bilinear form ℎ (·, ·) is positive definite and hence inverse exists. Further, let us assume that A, B be the matrix representation corresponding to the discrete forms ℎ (·, ·), ℎ (·, ·), respectively. Therefore, problem (4.1) reduces to a system of nonlinear differential equations as follows where ℎ denotes the vector whose entries are the components in the basis of ℎ . Moreover, C( ℎ ) is the matrix corresponding to the nonlinear term and F be the right hand side load vector corresponding to the basis ℎ . Before going into further details, we would like to prove that the nonlinear term , i.e., ℎ ( ℎ ; ℎ , ℎ ) satisfies Lipschitz's continuity condition. In this direction, let 1 ℎ , 2 ℎ be two elements in ℎ . Then, we can write as An application of Hölder inequality and using the continuity of Π −1 with respect to 4 -norm and stability of Π ,Δ , we obtain ∑︁ using Sobolev's embedding theorem, we obtain Using analogous arguments, we derive that Inserting (4.5) and (4.6) into (4.4), we can claim that the nonlinear term ℎ ( ℎ ; ℎ , ℎ ) is Lipschitz continuous. Therefore, from Picard's Theorem of existence and uniqueness of system of differential equation, we can deduce that (4.2) and (4.3) has a unique solution.

Fully-discrete formulation
A classical backward Euler integration method is employed for the time discretization of (4.1) with time step ∆ = / , where is a positive integer. In addition, we introduce ℎ := ℎ ( ) for = 0, 1, 2, . . . , . This results in the following fully discrete method: find ℎ ∈ ℎ such that where 0 ℎ ∈ ℎ is an initial approximation of at = 0. Next, we prove the well posedness of the fullydiscrete scheme (4.7). In this direction, we first recollect Brouwer's fixed point theorem. Then, under certain assumption, we will show that the fully-discrete scheme has unique solution ℎ and the solution is bounded, i.e., ‖ ℎ ‖ 2,Ω ≤ ℛ, where ℛ is a positive constant which will be specified later.
Let be a Banach space and let ℬ ⊂ be a compact and convex subset. If ℒ : ℬ → ℬ is continuous, then ℒ has a fixed point.
Then, for 1 ≤ ≤ and for sufficiently small values of ∆ , there exists a unique solution of the fully-discrete problem (4.7) and the solution ℎ satisfies the condition ‖ ℎ ‖ 2,Ω ≤ ℛ, with The proof of the result will be divided in three steps. We first define a mapping ℱ from ℎ to ℎ and prove that the mapping is well-defined and maps a ball ℬ ℛ to a ball ℬ ℛ . In second stage, we prove that the mapping is continuous. Then, from Brouwer's Theorem, we deduce that ℱ has a fixed point inside the ball ℬ ℛ which is the solution of the fully-discrete scheme (4.7). Finally, using assumption (4.8), we prove that the solution is unique.
Uniqueness of solution: Let 1 ℎ , 2 ℎ ∈ ℬ ℛ be two solutions of (4.7). Then from (4.9), we have To bound the nonlinear term in (4.19), we add and substract ℎ (︀ 2 ℎ ; 1 ℎ , ℎ )︀ in the above equality, to get (4.20) Using Hölder's inequality and Sobolev embedding theorem, the first two terms in (4.20) can be bounded as follows The last two terms on the right hand side of (4.20) can be written as follows, (4.19), we have that the right hand side above vanish (cf. (3.18)), and using the stability property of the discrete bilinear forms ℎ (·, ·) and ℎ (·, ·), and Young's inequality, we obtain An application of kick back arguments, we obtain  )︃ > 0.
Then for sufficiently small values of ∆ , we have

Convergence analysis
In this section, we will derive a priori error estimation for the virtual element semi-discrete and fully-discrete schemes. With this aim, first we introduce a discrete energy projection operator ℎ : 2 0 (Ω) → ℎ , which is defined as follows: Upon exploiting the energy projection operator ℎ , we will split the error of the stream-function as Next, we define In what follows, we will prove approximations properties for ℎ , thus the first term ℎ will be easily bounded. Then, using the continuous problem (2.3) and the polynomial approximation properties, we bound the term ℎ . In this regard, we introduce the polynomial approximation property and the interpolation operator on the virtual element space ℎ . Further, to derive the a priori error estimations of the semi-discrete and fully-discrete schemes, some additional results are needed which will be presented in the next subsection.

Preliminary results
We start with the following approximation result, on star-shaped polygons, which is derived by interpolation between Sobolev spaces (see for instance [28], Thm. I.1.4 from the analogous result for integer values of ). We mention that this result has been stated in Proposition 4.2 of [6] for integer values and follows from the classical Scott-Dupont theory (see [15] and [5], Prop. 3.1): Proposition 5.1. There exists a constant > 0, independent of mesh size ℎ but depends on mesh regularity parameter (Assumption 3.1) such that for every ∈ ( ) there exists ∈ ( ), ≥ 0 such that with [ ] denoting largest integer equal or smaller than ∈ R. Now, we present an interpolation result in the virtual space ℎ (see [5,10]).
Proposition 5.2. Assume A1 and A2 are satisfied, then for all ∈ ( ) there exist ∈ ℎ and > 0 independent of ℎ such that where is independent of mesh size ℎ but depends on mesh regularity parameter (Assumption 3.1).

Error estimation for semi-discrete scheme
In this section, we will derive the error estimation for the semi-discrete scheme (cf. (4.1)). With this end, we state the following lemma.
where we have used the Sobolev embedding (Ω) ˓→ 4 (Ω) for ∈ (1/2, 1]. On the other hand, for the case 1 ≤ ≤ − 1, we proceed as follows. Using the orthogonality property of the projection operator Π −1 on the polynomial function of degree − 1 and Cauchy-Schwarz inequality, we bound the first term on the right hand side of (5.13).
An application of Hölder's inequality and Sobolev embedding theorem Lemma 4.2 of [9] yields, Collecting the above inequalities, we obtain for > 1/2 that Using Hölder's inequality, the term 2 in (5.13) can be bounded as follows Using the continuity of Π −1 on the space 4 ( ) and optimal approximation property of the polynomial projection operator, we have Also, the term can be bounded as where we have used that 1 (Ω) ˓→ 4 (Ω). Hence, by the Sobolev embedding 1+ (Ω) ˓→ ,4 (Ω), we can write Using the Hölder's inequality and suming over each element , the third term of (5.13) can be bounded as follows Therefore, from the fact that Π ,Δ is the projector defined by (3.1a), the continuity of Π −1 on the space 4 ( ) and Sobolev's embedding theorem, we obtain Finally, the proof follows by collecting all the estimations (5.14)-(5.16) and inserting into (5.13).
The following theorem provides the error estimates for the semi-discrete virtual element scheme.
Theorem 5.5. Let ( ) ∈ be the solution of problem (2.5) and let ℎ ( ) ∈ ℎ be the solution of problem (4.1). Assume that ( ) ∈ 2+ (Ω), ( ) ∈ 1+ (Ω), for 1 2 < ≤ − 1 and for almost all ∈ [0, ]. In addition, assume that , ℎ ∈ := {︁ ∈ : Then, there exists a positive generic constant , which could be depended on mesh regularity Assumption 3.1, Sobolev regularity of the solution , , but independents of mesh size ℎ such that the following estimation holds Proof. Upon applying the projection operator ℎ , we split the error as (see (5.2)): Since the estimation for ℎ is known from the Lemma 5.3, we attempt to bound ℎ . Using the semi-discrete scheme (4.1), definition of the projection operator ℎ (5.1) and the continuous weak formulation (2.3), we obtain Now, we will bound all the terms on the right hand side above. We start with the nonlinear terms ( ; , ℎ ) − ℎ ( ℎ ; ℎ , ℎ ). We rewrite the term as follows: The term 1 has been bounded in Lemma 5.4. Thus, we will bound the term 2 . We rewrite the term 2 as follows Using Hölder's inequality, the first term on the right hand side of (5.18) can be bounded as follows, where we have used the boundedness of Π ,Δ , Sobolev's embedding theorem and the approximation property of the operator ℎ (cf. (5.3)). Using Hölder's inequality, we can bound the second term on the right hand side of (5.18) as follows, where we have exploited Sobolev inequality and approximation property of ℎ (cf. (5.3)).
The proof is complete.

Error estimation for fully-discrete scheme
In this section, we would like to derive the error estimation for the fully-discrete scheme (cf. (4.7)). To derive the estimates, we first split the error as follows: where ℎ := ℎ − ℎ ( ) and ℎ := ℎ ( ) − ( ). Since the estimation for ℎ is known from the Lemma 5.3, we will focus in bounding the term ℎ . , and force function f but independent of mesh size ℎ and time steps ∆ such that the following estimation holds Proof. Using the fully discrete scheme (4.7), weak formulation (2.3), and the biharmonic projection operator In order to derive the desired estimation, we will put ℎ = ℎ into (5.23).
The nonlinear term can be rewritten as follows: Proceeding same as semi-discrete case, another term of (5.25) can be estimated as follows Using the projection operator Π ,Δ , and polynomial consistency property of the discrete bilinear form ℎ (·, ·), we obtain (5.28) Following the analogous arguments as Theorem 3.3 of [40], we can write the term Σ 1 as follows, Further, in order to bound Σ 2 , we derive the following estimate Adding and subtracting curl ( ( ) − ( −1 )), the term Σ 3 (cf. (5.28)) can be written as Using the approximation properties of the operator ℎ , we obtain Analogously, we have Inserting the above estimates in (5.28), and then using it together with (5.26) and (5.27) into (5.24), and multiplying by ∆ , we obtain that An application of the Young's inequality and kick-back argument yields Since ( ), ℎ ∈˜, from the assumption of the theorem, and using Young's inequality and iterating = 1, . . . , , we have )︁ .
Using the approximation properties for ℎ , and Proposition 5.2, we get )︁ .
Finally, using the fact that ℎ − ( ) = ℎ + ℎ , from the above estimation and approximation properties for ℎ (see Lem. 5.3), we have the desired thesis. Remark 5.7. In Theorems 5.5 and 5.6, we have chosen that the analytical solution at time = and fully discrete solution ℎ belong to a bounded subset of ℬ ℛ . To satisfy this condition, an additional condition has to be imposed on viscosity , which leads For sufficiently small values of ∆ and satisfying the above mentioned condition, we advocate a numerical approximation of (2.3) that converges optimally in both space and time variable. The primary advantage of this scheme is that the condition imposed on is independent of 1/∆ which ensures the robustness of the scheme for very small values of ∆ .

Numerical experiments
In this section, we report the results of four numerical tests carried out with the fully-discrete virtual element scheme proposed in Section 4, in order to validate the theoretical results presented in Section 5. We have developed a MATLAB code that implements the fully-discrete scheme for = 2 and = 3. For the time discretization, we employ a backward Euler scheme and for each time step, we use the Newton-Raphson method to solve the resulting nonlinear system, with maximum 10 iterations, a user specified tolerance tol:= 10 −8 , and taking in ℎ = 0 as initial guess. For our numerical tests, we have used different families of polygonal meshes (see Fig. 1 ℎ : Distorted concave rhombic quadrilaterals. In order to test the convergence properties of the VEM method, we measure the errors as the difference between the exact solution and the suitable projections of the numerical solution ℎ . More precisely, for the norm 2 (︀ 0, ; 2 (Ω) )︀ , we consider the following quantity:

Test 1. Homogeneous Dirichlet boundary conditions and initial data
In this numerical test, we solve the Navier-Stokes problem (2.1) on the square domain Ω := (0, 1) 2 . We take the load term f , boundary and initial conditions in such a way that the analytical solution is given by: )︂ , Table 1. Test 1. Errors for the stream-function ℎ in the discrete 2 (︀ 0, ; 2 (Ω) )︀ norm obtained with = 2, = 1 and 2 ℎ .  We report in Table 3 the errors E 2 ( ) for the family of meshes 4 ℎ and different refinement levels and time steps. In this example, the maximum number of iterations that are required for the Newton method is 5. Once again, it can be seen along the diagonal of Table 3 that the error E 2 ( ) reduces linearly with respect to ℎ, which is the expected order of convergence for = 2. In addition, we have highlighted the errors which are dominated by space in Table 3 for small values of time-step ∆ .
We consider the time interval [0, 1], the viscosity = 1 and we start the process with ℎ 0 = 1/8 and ∆ 0 = 1/8.  We report in Table 4 the errors E 2 ( ) obtained for = 2 and using the family of meshes 1 ℎ with different refinement levels and time steps. In this experiment, differently to the first and second tests, we can observe that, for small values of ℎ, the error E 2 ( ) reduces linearly with respect to ∆ (see the last row of Tab. 4), which is the expected order of convergence in time according to Theorem 5.6. We also observe that for big values of ∆ the error E 2 ( ) is almost constant with respect to ℎ. In this numerical test, the maximum number of iterations that are required for the Newton method is 3.
Further, in Figure 3, we have posted a convergence graph where errors E 2 ( ) are dominated by time and virtual element space of order = 3 is chosen. We have used all family of meshes. We deduce that the rate of convergence is closer to 1 (expected order of convergence for time variable) for very small values of ℎ.

Test 4. Chorin problem
In this example, we consider the well-known Chorin problem for incompressible Navier-Stokes equations [22]. For this test the load term is f = 0, and the initial and boundary conditions correspond to the analytical solution: ( , , ) =  It has been observed that some finite element methods for a velocity-pressure formulation, the 2 error of the velocity converge suboptimally or even lock for small values of the viscosity (see [32], Sect. 4.2). It can be seen that the 2 error of the velocity is related with the 1 error for the stream-function. Thus, in order to assess the performance of the virtual scheme for this numerical example, we introduce the following discrete quantity: We observe that an additional order of convergence is expected in this discrete error. To show this fact, we will compute experimental rates of convergence for each individual error as follows: r ( ) := log(E ( )/E ′ ( )) log(ℎ/ℎ ′ ) , = 1, 2, where ℎ, ℎ ′ denote two consecutive mesh sizes with their respective errors E and E ′ . We report in Table 5 the discrete errors E 1 ( ) and E 2 ( ), for the family of meshes 3 ℎ . The results were obtained by considering the final time = 0.01 and time stepping ∆ = 0.001. For the viscosity , we take two values: = 10 −3 and = 10 −6 . The maximum number of iterations that are required for the Newton method in this example is 4 for = 10 −3 and 5 for = 10 −6 .
It can be clearly observed from Table 5 a linear order of convergence in the norm E 2 ( ) and a quadratic order in the norm E 1 ( ) (which has not been proved). Thus, we conclude that our virtual scheme does not suffer from a suboptimal order of convergence or locking phenomenon.
Exact and approximate solutions (including a postprocessed velocity field) are illustrated in Figure 4.

Conclusion
In this article, we have proposed a 1 VEM to discretize the time-dependent Navier-Stokes problem in the stream-function form. Exploiting enhanced virtual element spaces, we have approximated the spatial variables and we have discretized the time variable using the backward Euler scheme. We have derived a priori error estimations for semi-discrete and fully-discrete schemes and the theoretical results are verified by four numerical experiments. Moreover, the fourth numerical experiment allows us to conclude that our virtual scheme does not suffer from a suboptimal order of convergence when the diffusion coefficient is small, in contrast to some finite element methods in velocity-pressure formulation, where a suboptimal convergence is observed in 2 -norm of the velocity for the Chorin problem with small values of (see [32]).