The spectral implementation of the nonstationary Stokes problem with nonstandard boundary conditions

Abstract

This paper deals with the implementation of the spectral discretization of the vorticity–velocity–pressure formulation of the nonstationary Stokes problem. We use the implicit Euler scheme for the discretization in time. We propose a global method for the resolution of the linear system resulting from the discrete problem. This method is implemented, and some numerical results are presented which confirm the good convergence for the three unknowns.

1 Introduction

The low velocity flow of an incompressible viscous fluid can be modeled by the nonstationary Stokes equations. Such problem is typically provided with a Dirichlet boundary condition on the velocity [1–3]. In this paper, we deal with this problem in two- or three-dimensional bounded domain using nonstandard boundary conditions. In dimension 2, the problem is supplied with boundary conditions on the normal component of the velocity and the Dirichlet condition on the vorticity (since it is a scalar). However, in dimension 3, it is supplied with boundary conditions on the normal component of the velocity and on the tangential components of the vorticity vector field. This permits us to prove that the above problem is equivalent to a time-dependent variational formulation with unknown vorticity, velocity, and pressure [4–7]. For the discretization of this problem, we combined the implicit Euler’s scheme in time and the spectral method in space. One can see [8] for the analysis of this discretization in the case of the stationary Stokes problem.

In [9], the Generalized Minimal Residual method (GMRES) [10] is used to solve the stationary Stokes problem. However, this method is not easy to implement since the matrix is not symmetric and needs a huge number of iterations to converge. In [11], a direct method was proposed permitting to simplify and to improve the resolution.

For the nonstationary case, which is the subject of this paper, we will use the global resolution based on the good results established in [11]. The discretization in time by the implicit Euler method permits us to stabilize the discrete problem. The resulting global matrix is now symmetric and positive definite, making possible the use of the gradient conjugate method. As a result, we obtained high accuracy with an optimized number of iterations. Some numerical results are presented to show the time and space convergence.

The paper is organized as follows:

• In Sect. 2 the continuous problem and some regularity results are presented.

• Sect. 3 deals with the discrete problem and error estimates.

• The nonhomogeneous case is handled in Sect. 4.

• The implementation details of the discrete problem are developed in Sect. 5.

• Sect. 6 is dedicated to describe thoroughly the numerical results.

2 The continuous problem

Let Ω be a bounded simply connected domain of $\mathbb{R}^{d}$ ($d=2$ or $d=3$), and Γ be its boundary that we suppose connected and Lipschitz; $\mathbf{x}=(x,y)$ in $d=2$ or $\mathbf{x}=(x,y,z)$ in $d=3$ is the generic point in Ω. We introduce $[0,T]$ an interval of $\mathbb{R}$ where T is a positive constant and n is the unit outward vector to Ω on Γ. The nonstationary Stokes problem is presented as follows:

$$ \textstyle\begin{cases} \frac{{\partial \mathbf{u}}}{{\partial t}}(\mathbf{x},t)- \nu \Delta \mathbf{u}(\mathbf{x},t)+\nabla p(\mathbf{x},t)=\mathbf{f}(\mathbf{x},t) &\text{in } \varOmega \times [0,T], \\ \operatorname{div}\, \mathbf{u}(\mathbf{x},t)=0 &\text{in } \varOmega \times [0,T], \\ \mathbf{u}(\mathbf{x},t) \cdot \mathbf{n}(\mathbf{x})=0 &\text{on } \varGamma \times [0,T], \\ \zeta (\mathbf{curl}\, \mathbf{u}) (\mathbf{x},t)=\mathbf{0} & \text{on } \varGamma \times [0,T], \\ \mathbf{u}(\mathbf{x},0)=\mathbf{u}_{0} & \text{in } \varOmega , \end{cases} $$

(1)

where f is the data function and ν is the positive constant viscosity. The unknowns are the velocity u and pressure p. The boundary operator ζ is defined according to the dimension as follows:

• In dimension 2 ($d=2$), ζ represents the trace operator on the boundary Γ.

• In dimension 3 ($d=3$), if the vector field $\mathbf{w}=(w_{x},w_{y},w_{z})$, then $\zeta (w)=\mathbf{curl}\, w \times \mathbf{n}$ on the boundary Γ.

Using the formula

$$ -\Delta \mathbf{w}= \mathbf{curl}(\mathbf{curl}\, \mathbf{w}) - \nabla (\operatorname{div}\, \mathbf{w}), $$

and the vorticity ${{\tau }}=\mathbf{curl}\, \mathbf{u}$, we prove that (1) is equivalent to the following system:

$$ \textstyle\begin{cases} \frac{{\partial \mathbf{u}}}{{\partial t}}(\mathbf{x},t)+\nu\, \mathbf{curl}\, {{\tau }}( \mathbf{x},t)+\nabla p(\mathbf{x},t)=\mathbf{f}(\mathbf{x},t) &\text{in } \varOmega \times [0,T], \\ \operatorname{div}\, \mathbf{u}(\mathbf{x},t)=0 &\text{in } \varOmega \times [0,T], \\ {{\tau }}(\mathbf{x},t)=\mathbf{curl}\, \mathbf{u}(\mathbf{x},t) &\text{in } \varOmega \times [0,T], \\ \mathbf{u}(\mathbf{x},t) \cdot \mathbf{n}(\mathbf{x})=0 &\text{on } \varGamma \times [0,T], \\ \zeta ({{\tau }})(\mathbf{x},t)=\mathbf{0} & \text{on } \varGamma \times [0,T], \\ \mathbf{u}(\mathbf{x},0)=\mathbf{u}_{0} &\text{in } \varOmega . \end{cases} $$

(2)

We suppose that the initial velocity and vorticity satisfy

$$ \operatorname{div}\, \mathbf{u}_{0}=0 \quad \text{and}\quad {{\tau }}(\mathbf{x},0)={{\tau }}_{0}=\mathbf{curl}\, \mathbf{u}_{0}\quad \text{in } \varOmega . $$

(3)

We need to introduce the following Sobolev spaces:

$$ W^{m,p}(\varOmega )=\bigl\{ \varphi \in L^{p}(\varOmega ): \partial ^{\alpha } \varphi \in L^{p}(\varOmega ), \forall \vert \alpha \vert \leq m\bigr\} $$

associated with the following norm and seminorm:

$$ \Vert \varphi \Vert _{m,p,\varOmega }= \biggl(\sum _{ \vert \alpha \vert \leq m} \int _{\varOmega } \bigl\vert \partial ^{\alpha }\varphi (\mathbf{x}) \bigr\vert ^{p} \biggr)^{\frac{1}{p}}\quad \text{and}\quad \vert \varphi \vert _{m,p,\varOmega }= \biggl(\sum_{ \vert \alpha \vert = m} \int _{\varOmega } \bigl\vert \partial ^{ \alpha }\varphi (\mathbf{x}) \bigr\vert ^{p} \biggr)^{\frac{1}{p}}. $$

Also $H^{m}(\varOmega )=W^{m,2}(\varOmega )$ is a Hilbert space provided with the scalar product

$$ (\varphi ,\psi )_{m,\varOmega }= \biggl(\sum_{| \alpha | \leq m} \bigl( \partial ^{\alpha }\varphi ,\partial ^{\alpha }\psi \bigr)^{2} \biggr)^{\frac{1 }{2}}, $$

where $(\cdot ,\cdot )$ is the $L^{2}(\varOmega )$ scalar product. We denote by $L_{0}^{2}(\varOmega )$ the space of functions in $L^{2}(\varOmega )$ which have a null integral on Ω, and $\mathcal{D}(\varOmega )$ the space of infinitely many times differentiable functions with a compact support in Ω. Then we introduce the space $H({\mathrm{div}},\varOmega )$ as

$$ H({\mathrm{div}},\varOmega )=\bigl\{ \boldsymbol {\varphi }\in \bigl(L^{2}(\varOmega ) \bigr)^{d}: \operatorname{div}\, \boldsymbol {\varphi }\in L^{2}(\varOmega )\bigr\} $$

with the norm

$$ \Vert \boldsymbol {\varphi }\Vert _{H({\mathrm{div}},\varOmega )}= \bigl( \Vert \boldsymbol {\varphi }\Vert ^{2}_{L^{2}(\varOmega )^{d}} + \Vert \operatorname{div}\, \boldsymbol {\varphi }\Vert ^{2}_{L^{2}( \varOmega )} \bigr)^{\frac{1}{2}} $$

and

$$ H_{0}({\mathrm{div}},\varOmega )=\bigl\{ \boldsymbol {\varphi }\in H({\mathrm{div}}, \varOmega ): \boldsymbol {\varphi }\cdot \mathbf{n}=0 \text{ on } \partial \varOmega \bigr\} . $$

In the same way, we introduce the space $H({\mathbf{curl}},\varOmega )$ by

$$ H({\mathbf{curl}},\varOmega )=\bigl\{ \boldsymbol {\varphi }\in L^{2}(\varOmega )^{\frac{d(d-1)}{2}}: \mathbf{curl}\, \boldsymbol {\varphi }\in L^{2}(\varOmega )^{d} \bigr\} $$

with the norm

$$ \Vert \boldsymbol {\varphi }\Vert _{H({\mathbf{curl}},\varOmega )}= \bigl( \Vert \boldsymbol {\varphi }\Vert ^{2}_{L^{2}(\varOmega )^{\frac{d(d-1)}{2}}} + \Vert \mathbf{curl}\, \boldsymbol {\varphi }\Vert ^{2}_{L^{2}(\varOmega )^{d}} \bigr)^{\frac{1}{2}} $$

and

$$ H_{0}(\mathbf{curl},\varOmega )=\bigl\{ \boldsymbol {\varphi }\in H({\mathbf{curl}},\varOmega ): \boldsymbol {\varphi }\times \mathbf{n}=0 \text{ on } \partial \varOmega \bigr\} . $$

Remark 1

Notice that the spaces $H({\mathbf{curl}},\varOmega )$ and $H_{0}({\mathbf{curl}},\varOmega )$ coincide with the spaces $H^{1}(\varOmega )$ and $H^{1}_{0}(\varOmega )$ in dimension $d = 2$.

Now we introduce the following spaces to deal with time issues. Let B be a separable Banach space. We define $\mathcal{C}^{m}(0,T;B)$ as the set of $\mathcal{C}^{m}$-class time functions with values in B. Then $\mathcal{C}^{m}(0,T;B)$ is a Banach space with the norm

$$ \Vert \boldsymbol {\varphi }\Vert _{\mathcal{C}^{m}(0,T;B)}=\sup_{0\leq t\leq T}\sum _{k=0}^{m} \bigl\Vert \partial _{t}^{k} \boldsymbol {\varphi }\bigr\Vert _{B}, $$

where $\partial _{t}^{k} \boldsymbol {\varphi }$ is the partial derivative of order k in time of function φ. Consider also the spaces

$$ L^{p}(0,T;B)=\biggl\{ \boldsymbol {\varphi }:\boldsymbol {\varphi }\mbox{ measurable on } ]0,T[ \mbox{ such that } \int _{0}^{T} \bigl\Vert \boldsymbol {\varphi }(t,\cdot ) \bigr\Vert ^{p}_{B} \,dt< \infty \biggr\} $$

and

$$ H^{s}(0,T;B)=\bigl\{ \boldsymbol {\varphi }\in L^{2}(0,T;B):\partial ^{m} \boldsymbol {\varphi }\in L^{2}(0,T;B), m\leq s\bigr\} . $$

Then $L^{p}(0,T;B)$ is a Banach space equipped with the norm

$$ \Vert \boldsymbol {\varphi }\Vert _{L^{p}(0,T;B)}= \textstyle\begin{cases} ( \int _{0}^{T} \Vert \boldsymbol {\varphi }(t,\cdot ) \Vert ^{p}_{B} \,dt)^{\frac{1 }{p}} &\mbox{for } 1 \leq p < + \infty , \\ \sup_{0\leq t\leq T} \Vert \boldsymbol {\varphi }(t,\cdot ) \Vert _{B}& \mbox{for } p = + \infty , \end{cases} $$

and $H^{s}(0,T;B)$ is a Hilbert space with the scalar product

$$ (\boldsymbol {\varphi },\psi )_{H^{s}(0,T;B)}= \Biggl((\boldsymbol {\varphi },\psi )^{2}_{L^{2}(0,T;B)} + \sum_{m=0}^{s}\bigl(\partial ^{m} \boldsymbol {\varphi },\partial ^{m} \psi \bigr)^{2}_{L^{2}(0,T;B)} \Biggr)^{\frac{1}{2}}. $$

If the data $\mathbf{f}\in L^{2}(0,T;(H_{0}({\mathrm{div}},\varOmega ))^{\prime })$, where $(H_{0}({\mathrm{div}},\varOmega ))^{\prime }$ is the dual space of $H_{0}({\mathrm{div}},\varOmega )$ (see [12] for some properties of $(H_{0}({\mathrm{div}},\varOmega ))^{\prime }$), then we show, using the density of $(\mathcal{D}(\varOmega ))^{d}$ in $H_{0}(\operatorname{div},\varOmega )$ and the density of $\mathcal{D}(\varOmega )^{\frac{{d(d-1)}}{2}}$ in $H_{0}(\mathbf{curl},\varOmega )$ (see [2, Chap. I, Sect. 2]), that problem (2) is equivalent to the following variational formulation:

Find $({{\tau }},\mathbf{u},p) \in L^{2}(0,T;H_{0}(\mathbf{curl},\varOmega ))\times L^{2}(0,T;H_{0}( \operatorname{div},\varOmega ))\times L^{2}(0,T;L^{2}_{0}(\varOmega ))$

$$ \textstyle\begin{cases} \forall \mathbf{v}\in H_{0}(\operatorname{div},\varOmega ), \quad (\frac{{\partial \mathbf{u}}}{{\partial t}},\mathbf{v}) + a({{\tau }},\mathbf{u};\mathbf{v}) + b(\mathbf{v},p) =\prec \mathbf{f}, \mathbf{v}\succ , \\ \forall q\in L_{0}^{2}(\varOmega ), \quad b(\mathbf{u},q) =0, \\ \forall \boldsymbol {\vartheta }\in H_{0}(\mathbf{curl},\varOmega ), \quad c({{\tau }},\mathbf{u};\boldsymbol {\vartheta }) =0, \end{cases} $$

(4)

where ≺⋅, ⋅≻ is the duality product between $(H_{0}({\mathrm{div}},\varOmega ))^{\prime }$ and $H_{0}({\mathrm{div}},\varOmega )$. The bilinear forms $a(\cdot ,\cdot ;\cdot )$, $b(\cdot ,\cdot )$ and $c(\cdot ,\cdot ;\cdot )$ are defined as follows:

$$\begin{aligned}& a({{\tau }},\mathbf{u};\mathbf{v})=\nu \int _{ \varOmega }\,\mathbf{curl}({{\tau }}) (\mathbf{x},t)\cdot \mathbf{v}(\mathbf{x})\,d\mathbf{x},\qquad b(\mathbf{u},q)=- \int _{\varOmega }\operatorname{div}\, \mathbf{u}(\mathbf{x},t) q(\mathbf{x})\,d\mathbf{x}\quad \mbox{and} \\& c({{\tau }},\mathbf{u};\boldsymbol {\vartheta })= \int _{\varOmega }{{\tau }}( \mathbf{x},t)\cdot \boldsymbol {\vartheta }(\mathbf{x})\,d\mathbf{x}- \int _{\varOmega }\mathbf{u}(\mathbf{x},t) \cdot \mathbf{curl}\boldsymbol {\vartheta }(\mathbf{x})\,d\mathbf{x}. \end{aligned}$$

To prove the existence and uniqueness of the solution of problem (4), we need to define the kernel of the bilinear form $b(\cdot ,\cdot )$:

$$ V = \bigl\{ \boldsymbol {\varphi }\in H_{0}({\mathrm{div}},\varOmega ): \forall q \in L^{2}_{0}( \varOmega ), b(\boldsymbol {\varphi },q) = 0 \bigr\} , $$

which coincides with the space of divergence-free functions in $H_{0}(\operatorname{div},\varOmega )$.

We also introduce the kernel of the bilinear form $c(\cdot ,\cdot ;\cdot )$:

$$ \begin{aligned} U &= \bigl\{ (\boldsymbol {\vartheta },\boldsymbol {\varphi }) \in H_{0}({ \mathbf{curl}},\varOmega ) \times V: \forall \boldsymbol {\psi }\in H_{0}({ \mathbf{curl}},\varOmega ), c(\boldsymbol {\vartheta }, \boldsymbol {\varphi };\boldsymbol {\psi }) = 0 \bigr\} \\ &= \bigl\{ (\boldsymbol {\vartheta },\boldsymbol {\varphi }) \in H_{0}({\mathbf{curl}},\varOmega ) \times V: \boldsymbol {\varphi }= \mathbf{curl}\, \boldsymbol {\vartheta }\bigr\} . \end{aligned} $$

Then $({{\tau }},\mathbf{u})$ is a solution of the following reduced problem:

Find $({{\tau }},\mathbf{u}) \in L^{2}(0,T;U)$ such that

$$ \forall \mathbf{v}\in V,\quad \biggl(\frac{{\partial \mathbf{u}}}{{\partial t}},\mathbf{v}\biggr) + a( {{\tau }},\mathbf{u};\mathbf{v}) = \prec \mathbf{f},\mathbf{v}\succ . $$

(5)

We refer to [9, Lemma 2.3], [3, Chap. III, Theorem 1.1] and ([8, Proposition 1] regarding the arguments of the demonstration of the following proposition.

Proposition 1

For any data function$\mathbf{f}\in L^{2}(0,T;(H_{0}(\operatorname{div},\varOmega ))^{\prime })$, and$\mathbf{u}_{0}\in L^{2}(\varOmega )^{d}$which satisfies condition (3), problem (5) has a unique solution$({{\tau }},\mathbf{u})\in L^{2}(0,T;U)$satisfying

$$ \begin{aligned} &\Vert {{\tau }}\Vert ^{2}_{L^{\infty }(0,t;L^{2}( \varOmega )^{\frac{{d(d-1)}}{{2}}})} + \Vert \mathbf{u}\Vert ^{2}_{L^{\infty }(0,t;L^{2}(\varOmega )^{d})} \\ &\quad \leq c \bigl( \Vert {{\tau }}_{0} \Vert ^{2}_{L^{2}(\varOmega )^{\frac{{d(d-1)} }{{2}}}} + \Vert \mathbf{u}_{0} \Vert ^{2}_{L^{2}(\varOmega )^{d}} + \Vert \mathbf{f}\Vert ^{2}_{L^{2}(0,t;(H_{0}(\operatorname{div},\varOmega ))^{\prime })} \bigr), \end{aligned} $$

wherecis a positive constant that depends only onΩandT.

We recall the following inf–sup condition on bilinear form $b(\cdot ,\cdot )$ (see [2, Chap. I, Lemma 4.1]).

There exists a constant $\beta >0$ such that

$$ \forall q \in L^{2}_{0}(\varOmega ),\quad \sup_{\varphi \in H_{0}({ \mathrm{div}},\varOmega )} \frac{b(\varphi ,q)}{ \Vert \varphi \Vert _{H({ \mathrm{div}},\varOmega )}} \ge \beta \Vert q \Vert _{L^{2}(\varOmega )}. $$

(6)

Using Proposition 1 and (6), we obtain the proof of the following theorem (see [8, Sect. 1]):

Theorem 1

For a function$\mathbf{f}\in L^{2}(0,T;(H_{0}(\operatorname{div},\varOmega ))^{\prime })$, and$\mathbf{u}_{0}\in L^{2}(\varOmega )^{d}$which satisfies condition (3), problem (4) has a unique solution$({{\tau }},\mathbf{u},p)\in L^{2}(0,T;H_{0}(\mathbf{curl},\varOmega ))\times L^{2}(0,T;H_{0}( \operatorname{div},\varOmega ))\times L^{2}(0,T;L^{2}_{0}(\varOmega ))$. This solution satisfies

$$ \begin{aligned}[b] &\Vert {{\tau }}\Vert ^{2}_{L^{\infty }(0,t;L^{2}( \varOmega )^{\frac{{d(d-1)}}{{2}}})} + \Vert \mathbf{u}\Vert ^{2}_{L^{\infty }(0,t;L^{2}(\varOmega )^{d})}+ \Vert p \Vert ^{2}_{L^{2}(0,t;L_{0}^{2}( \varOmega ))} \\ &\quad \leq c \bigl( \Vert {{\tau }}_{0} \Vert ^{2}_{L^{2}(\varOmega )^{\frac{{d(d-1)} }{{2}}}} + \Vert \mathbf{u}_{0} \Vert ^{2}_{L^{2}(\varOmega )^{d}} + \Vert \mathbf{f}\Vert ^{2}_{L^{2}(0,t;(H_{0}(\operatorname{div},\varOmega ))^{\prime })} \bigr), \end{aligned} $$

(7)

wherecis a positive constant that depends onΩandT.

From [13, Sect. 2], [14], and [15], we derive the following regularity result.

Theorem 2

For a data functionfbelonging to$L^{2}(0,T;H^{r-1}(\varOmega )^{d})$, the solution$({{\tau }},\mathbf{u},p)$of problem (4) belongs to$L^{2}(0,T;H^{r}(\varOmega )^{\frac{d(d-1)}{2}}) \times L^{2}(0,T;H^{r}( \varOmega )^{d}) \times L^{2}(0,T;H^{r}(\varOmega ))$for any$r>0$such that:

• (i) $r\le \frac{1}{2}$in the general case,

• (ii) $r\le 1$whenΩis convex,

• (iii) $r < \frac{\pi }{\alpha }$in dimension$d = 2$whenΩis a polygon with largest angle equal toα.

3 Discrete problems and error estimates

3.1 The time discrete problem

In this section, we start by the discretization in time of problem (4) using implicit Euler method. We propose a partition of the interval $[0,T]$ into subintervals $[t_{k-1},t_{k}]$, for $1\leq k\leq K$, where K is a positive integer, and where

$$ 0=t_{0}< t_{1}< \cdots < t_{K}=T. $$

Let $h_{k}=t_{k}-t_{k-1}$, $h=(h_{1},h_{2},\ldots ,h_{K})$ and $|h|=\max_{1\leq k\leq K} h_{k}$.

For any function $\mathbf{f}\in {L^{2}(0,T;(H_{0}(\operatorname{div},\varOmega ))^{\prime })}$ and $u_{0}\in H_{0}(\operatorname{div},\varOmega )$ satisfying condition (3), we are considering the following scheme:

For all $k, 1 \leq k \leq K$, find $({{\tau }}^{k})_{0\leq k\leq K}\in (H_{0}(\mathbf{curl},\varOmega ))^{K+1}$, $(\mathbf{u}^{k})_{0\leq k\leq K}\in (H_{0}(\operatorname{div},\varOmega ))^{K+1}$ and $(p^{k})_{1\leq k\leq K}\in (L^{2}_{0}(\varOmega ))^{K}$ such that

$$\begin{aligned}& \mathbf{u}^{0}=\mathbf{u}_{0} \quad \text{and}\quad {{\tau }}^{0}=\mathbf{curl}\, \mathbf{u}^{0} \quad \text{in } {\varOmega }, \end{aligned}$$

(8)

$$\begin{aligned}& \textstyle\begin{cases} \forall \mathbf{v}\in H_{0}(\operatorname{div},\varOmega ), \quad (\mathbf{u}^{k} ,\mathbf{v}) + h_{k} a({{\tau }}^{k},\mathbf{u}^{k};\mathbf{v}) + h_{k} b(\mathbf{v},p^{k}) =(\mathbf{u}^{k-1},\mathbf{v})+h_{k} \prec \mathbf{f}^{k},\mathbf{v}\succ , \\ \forall q\in L_{0}^{2}(\varOmega ), \quad b(\mathbf{u}^{k},q) =0, \\ \forall \boldsymbol {\vartheta }\in H_{0}(\mathbf{curl},\varOmega ), \quad c({{\tau }}^{k},\mathbf{u}^{k};\boldsymbol {\vartheta }) =0, \end{cases}\displaystyle \end{aligned}$$

(9)

where $\mathbf{f}^{k}=\mathbf{f}(\cdot ,t_{k})$.

For any k, $1 \leq k \leq K$, let

$$ \hat{a}\bigl({{\tau }}^{k},\mathbf{u}^{k};\mathbf{v}\bigr)=\bigl(\mathbf{u}^{k} ,\mathbf{v}\bigr) + h_{k}a\bigl({{\tau }}^{k},\mathbf{u}^{k}; \mathbf{v}\bigr)\quad \text{and}\quad L(\mathbf{v})=\bigl(\mathbf{u}^{k-1}, \mathbf{v}\bigr)+h_{k}\prec \mathbf{f}^{k},\mathbf{v}\succ . $$

The functional L is linear and continuous on V. Thus if $({{\tau }}^{k},\mathbf{u}^{k},p^{k})$ is the solution of problem (9) then $({{\tau }}^{k},\mathbf{u}^{k})$ belongs to U and is the solution of the following reduced problem:

$$ \forall \mathbf{v}\in V,\quad \hat{a}\bigl({{\tau }}^{k},\mathbf{u}^{k};\mathbf{v}\bigr)=L(\mathbf{v}). $$

(10)

Based on the property of the bilinear form $\hat{a}(\cdot ,\cdot ,\cdot )$, proved in [8, Sect. 2], we state the following proposition.

Proposition 2

Suppose that the data function$\mathbf{f}\in L^{2}(0,t;H_{0}(\operatorname{div},\varOmega )^{\prime })$and that the initial velocity$\mathbf{u}_{0} \in H_{0}(\operatorname{div},\varOmega )$satisfies condition (3). Then, for any$k, 1 \leq k \leq K$, problem (8), (10) has a unique solution$({{\tau }}^{k},\mathbf{u}^{k})$inUsuch that

$$ \bigl\Vert \mathbf{u}^{k} \bigr\Vert ^{2}_{L^{2}(\varOmega )^{d}}\leq c \Biggl( \Vert \mathbf{u}_{0} \Vert ^{2}_{L^{2}(\varOmega )^{d}} + \sum_{j=1}^{k} h_{j} \bigl\Vert \mathbf{f}^{j} \bigr\Vert ^{2}_{(H_{0}(\operatorname{div},\varOmega ))^{\prime }} \Biggr), $$

(11)

wherecis a positive constant independent ofk.

Then, according to the inf–sup condition (6), we can state the following theorem (see [8, Sect. 2] for its proof).

Theorem 3

Suppose that the data function$\mathbf{f}\in L^{2}(0,t;H_{0}(\operatorname{div},\varOmega )^{\prime })$and that the initial velocity$\mathbf{u}_{0} \in H_{0}(\operatorname{div},\varOmega )$satisfies condition (3). Then, for any$k, 1 \leq k \leq K$, problem (8), (9) has a unique solution$({{\tau }}^{k},\mathbf{u}^{k},p^{k})$in$H_{0}(\mathbf{curl},\varOmega )\times H_{0}(\operatorname{div},\varOmega )\times L^{2}(\varOmega )$.

3.2 The spectral discrete problem

Hereinafter, we assume that the domain Ω is a square in dimension $d=2$ or a cube in dimension $d=3$. We adopt the Nédélec’s finite element method analog for the spectral discretization [16, Sect. 2].

To define the discrete spaces, we introduce $\mathbb{P}_{n,m}(\varOmega )$ as the space of polynomials with degree ≤n with respect to x and ≤m with respect to y, and $\mathbb{P}_{n,m,p}(\varOmega )$ as the space of polynomials with degree ≤n with respect to x, ≤m with respect to y, and ≤p with respect to z. When $n=m=p$, these spaces are equal to $\mathbb{P}_{n}(\varOmega )$.

Let N be an integer ≥2. The space ${\mathbb{D}}_{N}$ which approximates $H_{0}({\mathrm{div}},\varOmega )$ is defined by

$$ {\mathbb{D}}_{N} = H_{0}({ \mathrm{div}},\varOmega ) \cap \textstyle\begin{cases} \mathbb{P}_{N,N-1}(\varOmega ) \times \mathbb{P}_{N-1,N}(\varOmega ) & \text{if $d = 2$}, \\ \mathbb{P}_{N,N-1,N-1}(\varOmega ) \times \mathbb{P}_{N-1,N,N-1}(\varOmega ) \times \mathbb{P}_{N-1,N-1,N}( \varOmega ) &\text{if $d = 3$}. \end{cases} $$

(12)

The space ${\mathbb{C}}_{N}$ approximating $H_{0}({\mathbf{curl}},\varOmega )$ depends on the dimension (see Remark 1) and is defined by

$$ {\mathbb{C}}_{N} = \textstyle\begin{cases} H^{1}_{0}(\varOmega ) \cap \mathbb{P}_{N}(\varOmega ) &\text{if $d = 2$}, \\ H_{0}({\mathbf{curl}},\varOmega ) \cap (\mathbb{P}_{N-1,N,N}(\varOmega ) \times \mathbb{P}_{N,N-1,N}(\varOmega ) \times \mathbb{P}_{N,N,N-1}(\varOmega ) ) &\text{if $d = 3$}. \end{cases} $$

(13)

Finally, we consider the space ${\mathbb{M}}_{N}$ for the approximation of $L^{2}_{0}(\varOmega )$ given by

$$ {\mathbb{M}}_{N} = L^{2}_{0}( \varOmega ) \cap \mathbb{P}_{N-1}(\varOmega ). $$

(14)

If $\xi _{0} = -1$ and $\xi _{N} = 1$, we consider the nodes $\xi _{i}$, $1 \le i \le N-1$, and the set of $N+1$ weights $\rho _{i}$, $0 \le i \le N$, of the Gauss–Lobatto quadrature formula. We denote $\mathbb{P}_{n}(-1,1)$ the space of restrictions to $[-1,1]$ of polynomials with degree ≤n, and then

$$ \forall \varPhi \in \mathbb{P}_{2N-1}(-1,1),\quad \int _{-1}^{1} \varPhi (\zeta ) \,d\zeta = \sum _{j = 0}^{N} \varPhi (\xi _{j}) \rho _{j}. $$

(15)

We also recall a property from [17, (13.20)], which is useful in what follows:

$$ \forall \varphi _{N} \in \mathbb{P}_{N}(-1,1),\quad \Vert \varphi _{N} \Vert _{L^{2}(-1,1)}^{2} \le \sum_{j = 0}^{N} \varphi _{N}^{2}(\xi _{j}) \rho _{j} \le 3 \Vert \varphi _{N} \Vert _{L^{2}(-1,1)}^{2}. $$

(16)

According to (16), we define the discrete product on $\mathbb{P}_{N}(\varOmega )$, for two continuous functions u and v by

$$ (u,v)_{N} = \textstyle\begin{cases} \sum_{i = 0}^{N} \sum_{j = 0}^{N} u(\xi _{i},\xi _{j})v(\xi _{i}, \xi _{j}) \rho _{i}\rho _{j} &\text{if $d = 2$}, \\ \sum_{i = 0}^{N} \sum_{j = 0}^{N} \sum_{k = 0}^{N} u(\xi _{i},\xi _{j}, \xi _{k})v(\xi _{i},\xi _{j},\xi _{k}) \rho _{i}\rho _{j}\rho _{k} &\text{if $d = 3$}. \end{cases} $$

(17)

Finally, let $\mathfrak{I}_{N}$ denote the Lagrange interpolation operator on the nodes $(\xi _{i},\xi _{j})$, $0\le i$, $j \le N$, in dimension $d = 2$, and on the nodes $(\xi _{i},\xi _{j},\xi _{k})$, $0 \le i$, $j,k \le N$, in dimension $d = 3$, with values in $\mathbb{P}_{N}(\varOmega )$.

Henceforth, we suppose that the data f is continuous on $\overline{\varOmega }\times [0,T]$. The spectral discrete problem is constructed from the time semidiscrete problem (8), (9) using the Galerkin method combined with numerical integration. It reads as follows:

If $\mathbf{u}_{N}^{0}=\mathfrak{I}_{N}(\mathbf{u}_{0})$, knowing $\mathbf{u}_{N}^{k-1}$, we seek $({{\tau }}^{k}_{N},\mathbf{u}^{k}_{N},p^{k}_{N})$ in ${\mathbb{C}_{N}}\times {\mathbb{D}}_{N} \times {\mathbb{M}}_{N}$ such that for any k, $1\leq k\leq K$,

$$ \begin{aligned} &\forall \mathbf{v}_{N} \in { \mathbb{D}}_{N}, \quad \bigl(\mathbf{u}_{N}^{k}, \mathbf{v}_{N}\bigr)_{N} +h_{k}a_{N} \bigl({{\tau }}^{k}_{N},\mathbf{u}^{k}_{N}; \mathbf{v}_{N}\bigr) + h_{k}b_{N}\bigl( \mathbf{v}_{N},p^{k}_{N}\bigr) \\ &\hphantom{\forall \mathbf{v}_{N} \in { \mathbb{D}}_{N},}\quad \quad = \bigl(\mathbf{u}_{N}^{k-1},\mathbf{v}_{N}\bigr)_{N}+h_{k} \bigl({\mathfrak{I}}_{N}\bigl( \mathbf{f}^{k}\bigr),\mathbf{v}_{N}\bigr)_{N} , \\ &\forall q_{N} \in {\mathbb{M}}_{N},\quad b_{N}\bigl(\mathbf{u}^{k}_{N},q_{N} \bigr) = 0, \\ &\forall \boldsymbol {\vartheta }_{N} \in {\mathbb{C}}_{N},\quad c_{N}\bigl({{\tau }}^{k}_{N},\mathbf{u}^{k}_{N}; \boldsymbol {\vartheta }_{N}\bigr) = 0, \end{aligned} $$

(18)

where the bilinear forms $a_{N}(\cdot ,\cdot ;\cdot )$, $b_{N}(\cdot ,\cdot )$, and $c_{N}(\cdot ,\cdot ;\cdot )$ are defined by

$$ \begin{aligned} & a_{N}\bigl({{\tau }}^{k}_{N},\mathbf{u}^{k}_{N}; \mathbf{v}_{N}\bigr) = \nu \bigl({ \mathbf{curl}}\, {{\tau }}^{k}_{N} ,\mathbf{v}_{N}\bigr)_{N},\qquad b_{N}(\mathbf{v}_{N},q_{N}) = - ({ \mathrm{div}}\, \mathbf{v}_{N}, q_{N})_{N} , \\ & c_{N}\bigl({{\tau }}^{k}_{N},\mathbf{u}^{k}_{N}; \boldsymbol {\varphi }_{N} \bigr) = \bigl({{\tau }}^{k}_{N},\boldsymbol {\varphi }_{N}\bigr)_{N} - \bigl(\mathbf{u}^{k}_{N}, {\mathbf{curl}}\, \boldsymbol {\varphi }_{N}\bigr)_{N}. \end{aligned} $$

(19)

Given that

$$\begin{aligned}& \hat{a}_{N}\bigl({{\tau }}^{k}_{N},\mathbf{u}^{k}_{N}; \mathbf{v}_{N}\bigr)=\bigl(\mathbf{u}_{N}^{k},\mathbf{v}_{N}\bigr)_{N} + h_{k}a_{N}\bigl({{\tau }}^{k}_{N}, \mathbf{u}^{k}_{N}; \mathbf{v}_{N}\bigr) \quad \text{and} \\& L_{N}( \mathbf{v}_{N})=\bigl(\mathbf{u}_{N}^{k-1},\mathbf{v}_{N}\bigr)_{N}+h_{k} \bigl({\mathfrak{I}}_{N}\bigl(\mathbf{f}^{k}\bigr), \mathbf{v}_{N}\bigr)_{N}, \end{aligned}$$

using (16) combined with Cauchy–Schwarz inequality, we prove that the bilinear forms $\hat{a}_{N}(\cdot ,\cdot ;\cdot )$, $b_{N}(\cdot ,\cdot )$, and $c_{N}(\cdot ,\cdot ;\cdot )$ are continuous on $({\mathbb{C}}_{N} \times {\mathbb{D}}_{N} ) \times {\mathbb{D}}_{N}$, ${\mathbb{D}}_{N} \times {\mathbb{M}}_{N}$, and $({\mathbb{C}}_{N} \times {\mathbb{D}}_{N} ) \times {\mathbb{C}}_{N}$, respectively, with norms bounded independently of N. Moreover, as a consequence of the exactness property (15), the forms $b(\cdot ,\cdot )$ and $b_{N}(\cdot ,\cdot )$ coincide on ${\mathbb{D}}_{N} \times {\mathbb{M}}_{N}$. We remark also that the functional $L_{N}$ is linear and continuous on ${\mathbb{D}}_{N}$.

Let

$$ V_{N} = \bigl\{ \mathbf{v}_{N} \in {\mathbb{D}}_{N}: \forall q_{N} \in { \mathbb{M}}_{N}, b_{N}(\mathbf{v}_{N},q_{N}) = 0 \bigr\} $$

be the kernel of the bilinear form $b_{N}(\cdot ,\cdot )$. It corresponds to the space of divergence-free polynomials in ${\mathbb{D}}_{N}$.

We consider also

$$ U_{N} = \bigl\{ (\boldsymbol {\vartheta }_{N},\mathbf{v}_{N}) \in { \mathbb{C}}_{N} \times V_{N}: \forall \boldsymbol {\varphi }_{N} \in {\mathbb{C}}_{N}, c_{N}(\boldsymbol {\vartheta }_{N},\mathbf{v}_{N}; \boldsymbol {\varphi }_{N}) = 0 \bigr\} , $$

the kernel of the bilinear form $c_{N}(\cdot ,\cdot ;\cdot )$.

We consider now the following reduced discrete problem:

If $\mathbf{u}_{N}^{0}={\mathfrak{I}}_{N}(\mathbf{u}_{0})$, knowing $\mathbf{u}^{k-1}$, we seek $({{\tau }}^{k}_{N},\mathbf{u}^{k}_{N})$ in $U_{N}$ such that for all k, $1\leq k\leq K$,

$$ \forall \mathbf{v}_{N} \in V_{N},\quad \hat{a}_{N}\bigl({{\tau }}^{k}_{N},\mathbf{u}^{k}_{N}; \mathbf{v}_{N}\bigr) = L_{N}(\mathbf{v}_{N}). $$

(20)

The well-posedness of the reduced discrete problem (20) is based on the positivity condition and the inf–sup condition satisfied by the bilinear form $\hat{a}_{N}(\cdot ,\cdot ;\cdot )$ proved in [8, Lemmas 1 and 2] (see [18, Theorem 2.1] and [19] for similar results in finite element methods). We can presently write the following proposition.

Proposition 3

At each time step and knowing$\mathbf{u}^{k-1}_{N}$, problem (20) admits a unique solution$({{\tau }}^{k}_{N},\mathbf{u}^{k}_{N})$with values in the space$U_{N}$. The discrete velocity satisfies, for anyk, $1\leq k\leq K$,

$$ \bigl\Vert \mathbf{u}^{k}_{N} \bigr\Vert ^{2}_{L^{2}(\varOmega )^{d}}\leq c \Biggl( \bigl\Vert \mathbf{u}_{N}^{0} \bigr\Vert ^{2}_{L^{2}(\varOmega )^{d}} + \sum_{j=1}^{k} h_{j} \bigl\Vert { \mathfrak{I}}_{N}\bigl(\mathbf{f}^{j}\bigr) \bigr\Vert ^{2}_{L^{2}(\varOmega )^{d}} \Biggr), $$

wherecis a positive constant independent ofNandk.

We recall the inf–sup condition on the bilinear form $b_{N}(\cdot ,\cdot )$. We refer to [20, Lemma 3.1] for details of its proof.

Lemma 1

There exists a positive constant β independent of N such that the bilinear form $b_{N}(\cdot ,\cdot )$ satisfies

$$ \forall q_{N} \in {\mathbb{M}}_{N}, \quad \sup _{\mathbf{v}_{N} \in { \mathbb{D}}_{N}} \frac{b_{N}(\mathbf{v}_{N}, q_{N})}{ \Vert \mathbf{v}_{N} \Vert _{H({ \mathrm{div}},\varOmega )}} \ge \beta \Vert q_{N} \Vert _{L^{2}(\varOmega )}. $$

Using above results, we obtain the following theorem.

Theorem 4

Suppose functionfis continuous on$\bar{\varOmega }\times [0,T]$and$\mathbf{u}_{0}$belongs to$H_{0}(\operatorname{div},\varOmega )$that satisfies condition (3). For anyk, $1\leq k\leq K$, problem (18) has a unique solution$({{\tau }}^{k}_{N},\mathbf{u}^{k}_{N},p^{k}_{N})$in${\mathbb{C}}_{N} \times {\mathbb{D}}_{N} \times {\mathbb{M}}_{N}$.

3.3 The error estimate

For the estimation of the error between the solution $({{\tau }}^{k},\mathbf{u}^{k},p^{k})$ of problem (8)–(9) and the solution $({{\tau }}^{k}_{N},\mathbf{u}^{k}_{N},p^{k}_{N})$ of problem (18), we need to define the following space for $r \ge 0$:

$$ H^{r}({\mathbf{curl}},\varOmega ) = \bigl\{ \boldsymbol {\varphi }\in H^{r}( \varOmega )^{\frac{d(d-1) }{2}}: {\mathbf{curl}}\, \boldsymbol {\varphi }\in H^{r}(\varOmega )^{d} \bigr\} . $$

We notice that in dimension $d = 2$ this space coincides with $H^{r+1}(\varOmega )$.

We have the following error estimate (see [8, Sect. 5] for the proof).

Theorem 5

Suppose that function$\mathbf{f}\in L^{2}(0,T;H^{\sigma }(\varOmega )^{d})$for a real number$\sigma > \frac{d}{2}$and$\mathbf{u}_{0}\in H^{\mu }(\varOmega )^{d}$for a real number$\mu > \frac{d}{2}$. For anyk, $1\leq k\leq K$, if the solution$({{\tau }}^{k},\mathbf{u}^{k},p^{k})$of problem (8)–(9) belongs to$H^{r}({\mathbf{curl}},\varOmega ) \times H^{r}(\varOmega )^{d} \times H^{r}( \varOmega )$for a real number$r \ge d-1$, then the following error estimate holds between this solution and the solution$({{\tau }}^{k}_{N}, \mathbf{u}^{k}_{N},p^{k}_{N})$of problem (18):

$$ \begin{aligned}[b] &\bigl\Vert {{\tau }}^{k}- {{\tau }}^{k}_{N} \bigr\Vert _{H({\mathbf{curl}}, \varOmega )} + \bigl\Vert \mathbf{u}^{k}-\mathbf{u}^{k}_{N} \bigr\Vert _{L^{2}\varOmega )^{d}} + \bigl\Vert p^{k}-p^{k}_{N} \bigr\Vert _{L^{2}(\varOmega )} \\ &\quad \le c N^{-\mu } \Vert \mathbf{u}_{0} \Vert _{H^{\mu }(\varOmega )^{d}}+c \vert h \vert \bigl(N^{-r} \bigl( \bigl\Vert {{\tau }}^{k} \bigr\Vert _{H^{r}({\mathbf{curl}},\varOmega )} + \bigl\Vert \mathbf{u}^{k} \bigr\Vert _{H^{r}(\varOmega )^{d}} \\ &\qquad {}+ \bigl\Vert p^{k} \bigr\Vert _{H^{r}(\varOmega )} \bigr) + N^{-\sigma } \Vert \mathbf{f}\Vert _{L^{2}(0,T;H^{\sigma }(\varOmega )^{d})} \bigr). \end{aligned} $$

(21)

Estimate (21) is optimal for the three unknowns, while the estimate of the pressure is not obtained for most spectral discretizations of the Stokes problem.

Corollary 1

Assume that the datafbelongs to$L^{2}(0,T;H^{\sigma }(\varOmega )^{d})$for a real number$\sigma > \frac{d}{2}$and$\mathbf{u}_{0}\in H^{\mu }(\varOmega )^{d}$for a real number$\mu > \frac{d}{2}$. Then for anyk, $1\leq k\leq K$, the following error estimate holds between the solution$({{\tau }}^{k}, \mathbf{u}^{k},p^{k})$of problem (8)–(9) and the solution$({{\tau }}^{k}_{N}, \mathbf{u}^{k}_{N},p^{k}_{N})$of problem (18):

$$ \begin{aligned} & \bigl\Vert {{\tau }}^{k}-{{\tau }}^{k}_{N} \bigr\Vert _{H({\mathbf{curl}}, \varOmega )} + \bigl\Vert \mathbf{u}^{k}-\mathbf{u}^{k}_{N} \bigr\Vert _{L^{2}\varOmega )^{d}} + \bigl\Vert p^{k}-p^{k}_{N} \bigr\Vert _{L^{2}(\varOmega )} \\ &\quad \le c \vert h \vert N^{-\min \{\sigma , \mu ,\sigma _{\varOmega }\}} \bigl( \Vert \mathbf{f}\Vert _{L^{2}(0,T;H^{\sigma }(\varOmega )^{d})}+ \Vert \mathbf{u}_{0} \Vert _{H^{\mu }(\varOmega )^{d}} \bigr), \end{aligned} $$

where$\sigma _{\varOmega }$is a real number ≥1 depending only onΩ.

4 The nonhomogeneous boundary conditions

In this section, we will focus on the nonhomogeneous boundary condition on the velocity and generalize the results of the previous sections to the following problem:

$$ \textstyle\begin{cases} \frac{{\partial \mathbf{u}}}{{\partial t}} + \nu \, {\mathbf{curl}}\, {{\tau }}+ {\mathbf{grad}}\, p = \mathbf{f}&\text{in $ \varOmega \times [0,T]$}, \\ {\mathrm{div}}\, \mathbf{u}= 0 & \text{in $\varOmega \times [0,T]$}, \\ {{\tau }}= {\mathbf{curl}}\, \mathbf{u}& \text{in $\varOmega \times [0,T]$}, \\ \mathbf{u}\cdot \mathbf{n}= g &\text{on $\varGamma \times [0,T]$}, \\ \zeta ({{\tau }}) = {\mathbf{0}} & \text{on $\varGamma \times [0,T]$}, \\ \mathbf{u}(\mathbf{x},0)=\mathbf{u}_{0} & \text{in } \varOmega , \end{cases} $$

(22)

where the function g belongs to $L^{2}(0,T;H^{-\frac{1}{2}}(\varGamma ))$ and satisfies the compatibility condition (here $\langle \cdot ,\cdot \rangle _{\varGamma }$ obviously denotes the duality pairing between $H^{-\frac{1}{2}}(\varGamma )$ and $H^{\frac{1}{2}}(\varGamma )$). For each $0\leq t\leq T$,

$$ \bigl\langle g(\cdot ,t), 1 \bigr\rangle _{\varGamma } = 0. $$

(23)

We consider the following variational problem:

Find $({{\tau }},\mathbf{u},p)$ in $L^{2}(0,T;H_{0}({\mathbf{curl}},\varOmega )) \times L^{2}(0,T;H({\mathrm{div}}, \varOmega )) \times L^{2}(0,T;L^{2}_{0}(\varOmega ))$ such that

$$ \mathbf{u}\cdot \mathbf{n}= g \quad \mbox{on } \varGamma \times [0,T] $$

(24)

and

$$ \begin{aligned} &\forall \mathbf{v}\in H_{0}({ \mathrm{div}},\varOmega ),\quad \biggl(\frac{{\partial \mathbf{u}}}{{\partial t}},\mathbf{v}\biggr) + a({{\tau }},\mathbf{u};\mathbf{v}) + b(\mathbf{v},p) = \prec \mathbf{f},\mathbf{v}\succ , \\ &\forall q \in L^{2}_{0}(\varOmega ),\quad b(\mathbf{u},q) = 0, \\ &\forall \boldsymbol {\vartheta }\in H_{0}({\mathbf{curl}},\varOmega ),\quad c({{\tau }},\mathbf{u};\boldsymbol {\vartheta }) = 0. \end{aligned} $$

(25)

Thanks to the arguments given in [8, Sect. 2], problem (22) is equivalent to (24)–(25). To prove the well-posedness of problem (24)–(25), we need a lifting of the boundary condition (24) given in the following lemma (see [2, Chap. I] for the proof).

Lemma 2

Assumegis in$L^{2}(0,T;H^{-\frac{1}{2}}(\varGamma ))$and satisfies (23). Then there exists a divergence-free and curl-free function$\mathbf{u}_{b}$in$L^{2}(0,T;L^{2}(\varOmega )^{d})$such that$\mathbf{u}_{b} \cdot \mathbf{n}$is equal togon$\varGamma \times [0,T]$. Moreover, this function satisfies

$$ \Vert \mathbf{u}_{b} \Vert _{L^{2}(0,T;H({\mathrm{div}},\varOmega ))} \le c \Vert g \Vert _{L^{2}(0,T;H^{-\frac{1}{2}}(\varGamma ))}, $$

(26)

wherecis a positive constant.

We obtain the following result.

Theorem 6

We assume that the datafis in the dual space of$L^{2}(0,T;H_{0}({\mathrm{div}},\varOmega ))$andgis in$L^{2}(0,T;H^{-\frac{1}{2}}(\varGamma ))$and satisfies (23). The problem (24)–(25) has a unique solution$({{\tau }},\mathbf{u},p)$in$L^{2}(0,T;H_{0}({\mathbf{curl}},\varOmega )) \times L^{2}(0,T;H({\mathrm{div}}, \varOmega )) \times L^{2}(0,T;L^{2}_{0}(\varOmega ))$. Moreover, this solution satisfies

$$ \begin{aligned}[b] &\Vert {{\tau }}\Vert ^{2}_{L^{\infty }(0,t;L^{2}( \varOmega )^{\frac{{d(d-1)}}{{2}}})} + \Vert \mathbf{u}\Vert ^{2}_{L^{\infty }(0,t;L^{2}(\varOmega )^{d})}+ \Vert p \Vert ^{2}_{L^{2}(0,t;L_{0}^{2}( \varOmega ))} \\ &\quad \leq c \bigl( \Vert \mathbf{u}_{0} \Vert ^{2}_{L^{2}(\varOmega )^{d}} + \Vert \mathbf{f}\Vert ^{2}_{L^{2}(0,t;(H_{0}(\operatorname{div},\varOmega ))^{\prime })} + \Vert g \Vert ^{2}_{L^{2}(0,t;H^{-\frac{1}{2}}(\varGamma ))} \bigr), \end{aligned} $$

(27)

wherecis a positive constant.

Proof

Let $\mathbf{u}_{h} = \mathbf{u}-\mathbf{u}_{b}$, where $\mathbf{u}_{b}$ is the function introduced in Lemma 2. Then $({{\tau }},\mathbf{u},p)$ is a solution of problem (24)–(25) if and only if $({{\tau }}, \mathbf{u}_{h},p)$ is a solution of problem (4). Theorem 1 provides the existence and uniqueness of $({{\tau }},\mathbf{u},p)$. Combining (7) and (26), we derive the estimate (27). □

The time nonhomogeneous semidiscrete problem is written as follows:

For any function $\mathbf{f}\in {L^{2}(0,T;(H_{0}(\operatorname{div},\varOmega ))^{\prime })}$ and $u_{0}\in H_{0}(\operatorname{div},\varOmega )$ satisfying condition (3), find $({{\tau }}^{k})_{0\leq k\leq K}\in (H_{0}(\mathbf{curl},\varOmega ))^{K+1}$, $(\mathbf{u}^{k})_{0\leq k\leq K}\in (H(\operatorname{div},\varOmega ))^{K+1}$ and $(p^{k})_{1\leq k\leq K}\in (L^{2}_{0}(\varOmega ))^{K}$ such that

$$ \mathbf{u}^{0}=\mathbf{u}_{0}\quad \text{and}\quad {{\tau }}^{0}=\mathbf{curl}\, \mathbf{u}^{0}\quad \text{in } {\varOmega }, $$

for all $k, 1 \leq k \leq K$,

$$ \mathbf{u}^{k} \cdot \mathbf{n}= g^{k}\quad \mbox{on } \varGamma $$

(28)

and

$$\begin{aligned} \textstyle\begin{cases} \forall \mathbf{v}\in H_{0}(\operatorname{div},\varOmega ),\quad ( \mathbf{u}^{k} ,\mathbf{v}) + h_{k} a({{\tau }}^{k},\mathbf{u}^{k};\mathbf{v}) + h_{k} b(\mathbf{v},p^{k})=( \mathbf{u}^{k-1},\mathbf{v})+h_{k}\prec \mathbf{f}^{k},\mathbf{v}\succ , \\ \forall q\in L_{0}^{2}(\varOmega ),\quad b(\mathbf{u}^{k},q)=0, \\ \forall \boldsymbol {\vartheta }\in H_{0}(\mathbf{curl},\varOmega ),\quad c({{\tau }}^{k},\mathbf{u}^{k};\boldsymbol {\vartheta })=0, \end{cases}\displaystyle \end{aligned}$$

(29)

where $\mathbf{f}^{k}=\mathbf{f}(\cdot ,t_{k})$ and $g^{k}=g(\cdot ,t_{k})$. Let $\mathbf{u}^{k}_{b}=\mathbf{u}_{b}(\cdot ,t_{k})$ be such that $\mathbf{u}^{k}_{b} \cdot \mathbf{n}$ is equal to $g^{k}$ on Γ and

$$ \bigl\Vert \mathbf{u}^{k}_{b} \bigr\Vert _{H({\mathrm{div}},\varOmega )} \le c \bigl\Vert g^{k} \bigr\Vert _{H^{-\frac{1}{2}}(\varGamma )}. $$

(30)

Theorem 7

Assume that for any$k, 1 \leq k \leq K$, $\mathbf{f}^{k}\in (H_{0}(\operatorname{div},\varOmega )){'}$, $g^{k} \in H^{-\frac{1}{2}}(\varGamma )$and that the initial velocity$\mathbf{u}_{0} \in H_{0}(\operatorname{div},\varOmega )$satisfies condition (3). Then problem (28)–(29) has a unique solution$({{\tau }}^{k},\mathbf{u}^{k},p^{k})$in$H_{0}(\mathbf{curl},\varOmega )\times H(\operatorname{div},\varOmega )\times L^{2}(\varOmega )$such that

$$ \bigl\Vert \mathbf{u}^{k} \bigr\Vert ^{2}_{L^{2}(\varOmega )^{d}}\leq c \Biggl( \Vert \mathbf{u}_{0} \Vert ^{2}_{L^{2}(\varOmega )^{d}} + \sum_{j=1}^{k} h_{j} \bigl\Vert \mathbf{f}^{j} \bigr\Vert ^{2}_{(H_{0}(\operatorname{div},\varOmega ))^{\prime }}+ \bigl\Vert g^{j} \bigr\Vert ^{2}_{H^{-\frac{1}{2}}(\varGamma )} \Biggr), $$

(31)

wherecis a positive constant independent ofk.

Proof

For any $k, 1 \leq k \leq K$, if $\mathbf{u}^{k}_{h} = \mathbf{u}^{k}-\mathbf{u}^{k}_{b}$, then $({{\tau }}^{k},\mathbf{u}^{k},p^{k})$ is a solution of problem (28)–(29) if and only if $({{\tau }}^{k}, \mathbf{u}^{k}_{h},p^{k})$ is a solution of problem (8)–(9). Theorem 3 permits us to deduce the well-posedness of problem (28)–(29), and the estimate (31) is derived by combining (11) and (30). □

Since we handle the nonhomogeneous boundary condition on the velocity, we define the following new discrete space for the discrete velocity:

$$ \overline{\mathbb{D}}_{N} = \textstyle\begin{cases} \mathbb{P}_{N,N-1}(\varOmega ) \times \mathbb{P}_{N-1,N}(\varOmega ) &\text{if $d = 2$}, \\ \mathbb{P}_{N,N-1,N-1}(\varOmega ) \times \mathbb{P}_{N-1,N,N-1}(\varOmega ) \times \mathbb{P}_{N-1,N-1,N}( \varOmega ) &\text{if $d = 3$}. \end{cases} $$

Suppose that for any $k, 1 \leq k \leq K$, the function $g^{k}\in L^{2}(\varGamma )$ and $g^{k}_{N}$ is the approximation of $g^{k}$ defined as follows: On each edge ($d = 2$) or face ($d = 3$) $\varGamma _{r}$ of Ω, $1 \le r \le 2d$, $g^{k}_{N\vert \varGamma _{r}}$ is equal to the image of $g^{k}_{\vert \varGamma _{r}}$ by the orthogonal projection operator from $L^{2}(\varGamma _{r})$ onto $\mathbb{P}_{N-1}(\varGamma _{r})$. Thus we write the following spectral discrete problem:

If $\mathbf{u}_{N}^{0}=\mathfrak{I}_{N}(\mathbf{u}_{0})$, knowing $\mathbf{u}_{N}^{k-1}$, find $({{\tau }}^{k}_{N},\mathbf{u}^{k}_{N},p^{k}_{N})$ in ${\mathbb{C}_{N}}\times \overline{\mathbb{D}}_{N} \times {\mathbb{M}}_{N}$ such that for any $k, 1 \leq k \leq K$,

$$ \mathbf{u}^{k}_{N} \cdot \mathbf{n}= g^{k}_{N}\quad \mbox{on } \varGamma , $$

(32)

and

$$ \begin{aligned} &\forall \mathbf{v}_{N} \in { \mathbb{D}}_{N},\quad \bigl(\mathbf{u}^{k}_{N},v_{N} \bigr)_{N} + h_{k}a_{N}\bigl({{\tau }}^{k}_{N},\mathbf{u}^{k}_{N}; v_{N}\bigr) + h_{k}b_{N}\bigl(\mathbf{v}_{N},p^{k}_{N}\bigr) \\ &\hphantom{\forall \mathbf{v}_{N} \in { \mathbb{D}}_{N},}\quad \quad = \bigl(\mathbf{u}^{k-1}_{N},v_{N}\bigr)_{N}+h_{k} \bigl( \mathbf{f}^{k},\mathbf{v}_{N}\bigr)_{N}, \\ &\forall q \in {\mathbb{M}}_{N},\quad b_{N}\bigl(\mathbf{u}^{k}_{N},q_{N}\bigr) = 0, \\ &\forall \boldsymbol {\vartheta }_{N} \in {\mathbb{C}}_{N},\quad c_{N}\bigl({{\tau }}^{k}_{N},\mathbf{u}^{k}_{N}; \boldsymbol {\vartheta }_{N}\bigr) = 0. \end{aligned} $$

(33)

Theorem 8

If for any$k, 1 \leq k \leq K$, the data function$\mathbf{f}^{k}$is continuous onΩ̅and$g^{k}$is in$L^{2}(\varGamma )$satisfying (23), the problem (32)–(33) has a unique solution$({{\tau }}^{k}_{N},\mathbf{u}^{k}_{N},p^{k}_{N})$in${\mathbb{C}}_{N} \times \overline{\mathbb{D}}_{N} \times {\mathbb{M}}_{N}$.

Proof

Writing problem (32)–(33) as a square linear system, we deduce from Theorem 4 that the unique solution of this problem when the data $\mathbf{f}^{k}$ and $g^{k}_{N}$ are zero is $({ 0}, { 0}, 0)$. This yields the existence and uniqueness property. □

We derive the following error estimate from those proved in the homogeneous case (see [8, Sect. 5]) using slight modifications.

Theorem 9

If the assumptions of Theorem 5hold and for any$k, 1 \leq k \leq K$, the data$g^{k}$satisfies condition (23) and is such that each$g^{k}_{\vert \varGamma _{r}}$, $1 \le r \le 2d$, belongs to$H^{\tau }(\varGamma _{r})$for a nonnegative real numberτ, then the following error estimate holds between the solution$({{\tau }}^{k},\mathbf{u}^{k},p^{k})$of problem (24)–(25) and the solution$({{\tau }}^{k}_{N},\mathbf{u}^{k}_{N},p^{k}_{N})$of problem (29)–(30):

$$ \begin{aligned}[b] &\bigl\Vert {{\tau }}^{k}- {{\tau }}^{k}_{N} \bigr\Vert _{H({\mathbf{curl}}, \varOmega )} + \bigl\Vert \mathbf{u}^{k}-\mathbf{u}^{k}_{N} \bigr\Vert _{H({\mathrm{div}},\varOmega )} + \bigl\Vert p^{k}-p^{k}_{N} \bigr\Vert _{L^{2}(\varOmega )} \\ &\quad \le c \Biggl(N^{-s} \bigl( \bigl\Vert {{\tau }}^{k} \bigr\Vert _{H^{s}({\mathbf{curl}}, \varOmega )} + \bigl\Vert \mathbf{u}^{k} \bigr\Vert _{H^{s}(\varOmega )^{d}} + \bigl\Vert p^{k} \bigr\Vert _{H^{s}(\varOmega )} \bigr) \\ &\qquad {} + N^{-\sigma } \bigl\Vert \mathbf{f}^{k} \bigr\Vert _{H^{\sigma }(\varOmega )^{d}} + N^{- \tau -\frac{1}{2}} \sum_{r = 1}^{2d} \bigl\Vert g^{k} \bigr\Vert _{H^{\tau }( \varGamma _{r})} \Biggr). \end{aligned} $$

(34)

Corollary 2

Assume that for any$k, 1 \leq k \leq K$, the data$(\mathbf{f}^{k}, g^{k})$belongs to$H^{\sigma }(\varOmega )^{d} \times H^{\sigma -\frac{1}{2}}(\varGamma )$for a real number$\sigma > \frac{d}{2}$and that condition (23) is satisfied. Then the following error estimate holds between the solution$({{\tau }}^{k},\mathbf{u}^{k},p^{k})$of problem (24)–(25) and the solution$({{\tau }}^{k}_{N},\mathbf{u}^{k}_{N},p^{k}_{N})$of problem (29)–(30):

$$ \begin{aligned}[b] & \bigl\Vert {{\tau }}^{k}- {{\tau }}^{k}_{N} \bigr\Vert _{H({\mathbf{curl}}, \varOmega )} + \bigl\Vert \mathbf{u}^{k}-\mathbf{u}^{k}_{N} \bigr\Vert _{H({\mathrm{div}},\varOmega )} + \bigl\Vert p^{k}-p^{k}_{N} \bigr\Vert _{L^{2}(\varOmega )} \\ &\quad \le c N^{-\min \{\sigma , \sigma _{\varOmega }\}} \bigl( \bigl\Vert \mathbf{f}^{k} \bigr\Vert _{H^{\sigma }(\varOmega )^{d}} + \bigl\Vert g^{k} \bigr\Vert _{H^{\sigma -\frac{1 }{2}}(\varGamma )} \bigr), \end{aligned} $$

(35)

where the real number$\sigma _{\varOmega }$is same as in Corollary 1.

5 The implementation of the discrete problem

In this section, we propose a global method for the resolution of the discrete problem (18). This method was used to solve the stationary Stokes problem for the same formulation (see [11]). This new algorithm enhanced the performance of the previous resolution of Stokes problem (2D and 3D bounded domain) (see [9]) and optimized the execution time. In the following, we start by presenting the linear system. Afterwards, we describe the resolution algorithm.

For the discrete problem (18), as a linear system, we have to choose a basis of the discrete spaces ${\mathbb{C}}_{N}$, ${\mathbb{D}}_{N}$, and ${\mathbb{M}}_{N}$.

The Lagrange polynomials in $\mathbb{P}_{N}(-1,1)$ linked with the nodes $\xi _{j}$ are denoted as $\varphi _{j}$, $0 \le j \le N$. We define

$$ \varphi _{j}^{*}(\zeta ) = \varphi _{j}(\zeta ) \frac{\xi _{j}-\xi _{j^{*}} }{\zeta - \xi _{j^{*}} }, \quad j \in J^{*}, $$

(36)

where $J^{*}$ is the set $\{0,\ldots ,N\} \setminus \{j^{*}\}$ and $j_{*}$ is the integer part of $\frac{N}{2}$. For any k, $1\leq k\leq K$, the unknowns ${{\tau }}^{k}_{N}$, $\mathbf{u}^{k}_{N}=(u^{k}_{Nx},u^{k}_{Ny})$, and a pseudopressure $\tilde{p}^{k}_{N}$ admit the expansions, in dimension $d = 2$ for simplicity,

$$\begin{aligned}& {{\tau }}^{k}_{N}(x,y) = \sum _{i = 1}^{N-1} \sum _{j = 1}^{N-1} {{\tau }}^{k}_{ij} \varphi _{i}(x)\varphi _{j}(y), \\& u^{k}_{Nx}(x,y) = \sum _{i = 1}^{N-1} \sum_{j \in J^{*}} u^{xk}_{ij} \varphi _{i}(x)\varphi _{j}^{*}(y), \qquad u^{k}_{Ny}(x,y) = \sum_{i \in J^{*}} \sum_{j = 1}^{N-1} u^{yk}_{ij} \varphi _{i}^{*}(x) \varphi _{j}(y), \\& \tilde{p}_{N}(x,y) = \sum_{i\in J^{*}, j\in J^{*}, (i,j) \neq (0,0) } p^{k}_{ij} \varphi _{i}^{*}(x) \varphi _{j}^{*}(y). \end{aligned}$$

The function $\tilde{p}^{k}_{N}$ vanishes in $(-1,-1)$ but no longer belongs to $L^{2}_{0}(\varOmega )$; however, the real pressure $p^{k}_{N}$ can easily be recovered in a postprocessing step, thanks to the formula

$$ p^{k}_{N}(x,y) = \tilde{p}^{k}_{N}(x,y) - \frac{1}{2^{d}} \bigl( \tilde{p}^{k}_{N}, 1\bigr)_{N}. $$

(37)

For any k, $1\leq k\leq K$, we denote by $\varPhi ^{k}$, $U^{k}$, and $P^{k}$ the vectors made of these coefficients. Their respective dimensions are equal to $\frac{d(d-1)}{2} N^{d-2}(N-1)^{2}$, $d N^{d-1}(N-1)$, and $N^{d}-1$. Hereinafter, by supposing the viscosity $\nu =1$, problem (18) is equivalent to the following square linear system:

If $U^{0}=(U^{0}_{1},U^{0}_{2})$, the components of the vectors $U^{0}_{1}$ and $U^{0}_{2}$ are respectively $u^{1}_{0}(\xi _{i},\xi _{j})$ and $u^{2}_{0}(\xi _{i},\xi _{j})$ where $\mathbf{u}_{0}=(u^{1}_{0},u^{2}_{0})$, then for any k, $1\leq k\leq K$,

$$ \begin{pmatrix} D^{k} &-{A^{k}}^{T}&0 \\ -A^{k} &I^{k}&{B^{k}}^{T} \\ 0&B^{k} & 0 \end{pmatrix} \begin{pmatrix} \varPhi ^{k} \\ U^{k} \\ P^{k} \end{pmatrix} = \begin{pmatrix} 0 \\ F^{k} \\ 0 \end{pmatrix} , $$

(38)

where ${A^{k}}^{T}$ and ${B^{k}}^{T}$ denote the transposed matrices of $A^{k}$ and $B^{k}$, respectively.

Matrix$A^{k}$:

Matrix $A^{k}$ is written as

$$ A^{k}= \begin{pmatrix} A^{k}_{1} 0 \\ 0 A^{k}_{2} \end{pmatrix} . $$

For any k, $1\leq k\leq K$, ${\mathbf{curl}} (\tau ^{k}_{N})=(\partial _{y}\tau ^{k}_{N}, -\partial _{x} \tau ^{k}_{N})$, and then the coefficients of the matrices $A^{k}_{1}$ and $A^{k}_{2}$ are deduced respectively from the two terms $(\partial _{y}\tau ^{k}_{N},u^{k}_{Nx})_{N}$ and $(\partial _{x}\tau ^{k}_{N},u^{k}_{Ny})_{N}$. Thus,

$$ \bigl(\varphi _{i} \varphi _{j}^{\prime }, \varphi _{r}\varphi ^{*}_{s} \bigr)_{N}= \alpha (j,s)\delta _{ir}\rho _{r};\quad 1\leq i,j,r \leq N-1; s\in J^{*} $$

and

$$ \bigl(\varphi _{i}^{\prime } \varphi _{j} , \varphi ^{*}_{r}\varphi _{s} \bigr)_{N}= \alpha (i,r)\delta _{js}\rho _{s}; \quad 1\leq i,j,s \leq N-1; r\in J^{*}, $$

where

$$ \alpha (j,s)= \varphi _{j}^{\prime }(\xi _{s}) \rho _{s} + (\xi _{s} - \xi _{i^{*}}) \varphi _{j}^{\prime }(\xi _{i^{*}})\varphi _{s}^{\prime }(\xi _{i^{*}})\rho _{i^{*}} $$

and $\delta _{ij}$ is the Kronecker symbol.

We notice that the matrices $A^{k}_{1}$ and $A^{k}_{2}$ are not square. They have $N^{d-2}(N-1)^{2}$ lines and $\frac{d(d-1)}{2} N^{d-2}(N-1)^{2}$ columns.

Matrix$B^{k}$:

For any k, $1\leq k\leq K$, matrix $B^{k}$ is defined as

$$ B^{k}=\bigl[B^{k}_{1}, B^{k}_{2} \bigr]. $$

The coefficients of the matrices $B^{k}_{1}$ and $B^{k}_{2}$ are found respectively from the terms $(\partial _{x} u^{k}_{ Nx}, p^{k}_{N})_{N}$ and $(\partial _{y} u^{k}_{Ny}, p^{k}_{N})_{N}$. Let then

$$ \bigl(\varphi _{i}^{\prime } \varphi _{j}^{*}, \varphi _{r}^{*}\varphi ^{*}_{s} \bigr)_{N}= \alpha (i,r)\beta (j,s),\quad 1\leq i \leq N-1; j,r,s\in J^{*} $$

and

$$ \bigl(\varphi _{i}^{*} \varphi _{j}^{\prime } ,\varphi ^{*}_{r}\varphi _{s}^{*} \bigr)_{N}= \alpha (j,s)\beta (i,r),\quad 1\leq j\leq N-1; i,r,s \in J^{*}, $$

where

$$ \beta (r,s)=\delta _{rs}\rho _{r} + (\xi _{r} - \xi _{i^{*}}) (\xi _{s} - \xi _{i^{*}})\varphi _{r}^{\prime }(\xi _{i^{*}})\varphi _{s}^{\prime }(\xi _{i^{*}}) \rho _{i^{*}}. $$

We also note that matrices $B^{k}_{1}$ and $B^{k}_{2}$ are not square, having $N^{d}-1$ lines and $d N^{d-1}(N-1)$ columns.

Matrix$D^{k}$:

Matrix $D^{k}$ is a diagonal matrix with $\frac{d(d-1)}{2} N^{d-2}(N-1)^{2}$ lines and columns. Its coefficients are equal to

$$ (\varphi _{i} \varphi _{j} ,\varphi _{r} \varphi _{s})_{N}=\delta _{ir} \delta _{js}\rho _{r}\rho _{s},\quad 1\leq i,j,r,s \leq N-1. $$

(39)

Matrix$I^{k}$:

Matrix $I^{k}$ is written as

$$ I^{k}= \begin{pmatrix} I^{k}_{1} 0 \\ 0 I^{k}_{2} \end{pmatrix} . $$

For any k, $1\leq k\leq K$, the coefficients of the matrices $I^{k}_{1}$ and $I^{k}_{2}$ are respectively equal to

$$ \delta _{ir}\rho _{r}\beta (j,s),\quad 1\leq i,r\leq N-1; j,s\in J^{*} $$

and

$$ \delta _{js}\rho _{s}\beta (i,r),\quad 1\leq j,s \leq N-1; i,r\in J^{*}. $$

Lastly, for any k, $1\leq k\leq K$, we denote

$$ F^{k}= \begin{pmatrix} F^{k}_{1} \\ F^{k}_{2} \\ 0 \end{pmatrix} . $$

The components of the two entries $F^{k}_{1}$ and $F^{k}_{2}$ are respectively $(u^{k-1}_{Nx},\varphi _{r} \varphi _{s}^{*})_{N}+h^{k}(f^{k}_{1}, \varphi _{r} \varphi _{s}^{*})_{N}, 1\leq r\leq N-1; s\in I^{*}$ and $(u^{k-1}_{Ny},\varphi _{r}^{*}\varphi _{s})_{N}+h^{k}(f^{k}_{2}, \varphi _{r}^{*}\varphi _{s})_{N}, 1\leq s\leq N-1; r\in I^{*}$, where the data function is $\mathbf{f}=(f_{1},f_{2})$.

We solve the linear system (38) using the gradient conjugate method preconditioned by LU incomplete factorization.

6 Experimental results

6.1 Two-dimensional experiments

In this section, at first, we focus on the time convergence. We consider the square $\varOmega = \,]-1,1[^{2}$. We look at a given solution obtained from the formulas $\mathbf{u}= {\mathbf{curl}}\, \psi$ and ${{\tau }}= {\mathbf{curl}}\, \mathbf{u}$ in the two cases:

Case (i). Regular functions ψ and p, defined by

$$ \psi (x,y ) = t \sin (\pi x)\sin (\pi y) ,\qquad p(x,y) = e^{-t} xy. $$

(40)

Case (ii). Less regular functions ψ and p, defined by

$$ \psi (x,y ) = t \bigl(1-x^{2}\bigr)^{3} \bigl(1-y^{2}\bigr)^{\frac{7}{2}},\qquad p(x,y) =e^{-t} x \bigl(1-x^{2}\bigr)^{\frac{3}{2}} \bigl(1+y^{2} \bigr)^{-\frac{1}{2}}. $$

(41)

The velocity is a Gaussian which is null for $t = 0$. We consider the spectral discrete parameter $N=30$, $T=1$ and the time steps $h\in \{0.1,0.001,0.0001\}$.

Figure 1 presents the curves of convergence for the three terms $\log \|{{\tau }}-{{\tau }}^{n}_{N} \|_{H(\mathbf{curl},\varOmega )}$ (in red), $\log \|\mathbf{u}-\mathbf{u}^{n}_{N}\|_{H(\operatorname{div},\varOmega )}$ (in blue), and $\log \|p-p^{n}_{N}\|_{L^{2}(\varOmega )}$ (in green) as a function of $\log (h)$. Figures 1(a) and 1(b) correspond respectively to the resolution for the continuous solutions defined in (40) and (41). We notice that the time convergence order is almost equal to 1.

Figure 1

The time error curves

Hereinafter, we consider $h=0.0001$ and $T=1$,

Figure 2 presents the spectral error curves on the vorticity, velocity and the pressure ($\log (\vert \mathrm{error} \vert )$ as a function of $\log (N)$ for N varying from 5 to 30).

Figure 2

Error curves for the spectral discretization

In Fig. 2(a), the discrete solution is computed from (40). We obtain a very good (exponential) convergence due to spectral discretization.

In Fig. 2(b), the discrete solution is calculated from (41). We get lower convergence due to the irregularity of the solution.

We notice that the convergence of the pressure is of the same order as for the vorticity and velocity.

Figure 3 corresponds from top to bottom and left to right to the discrete vorticity, the two components of the discrete velocity, and the discrete pressure for the data

$$ \mathbf{f}= (f_{x}, f_{y})= \bigl(tx^{2}y, 0 \bigr), \qquad \mathbf{u}_{0}=(0,xy), $$

(42)

homogeneous boundary conditions $g=0$ and $N=30$.

Figure 3

The solution $({{\tau }}, u_{x}, u_{y},p)$ for the data f defined in (42) and $g = 0$

Figure 4 corresponds from top to bottom and left to right to the discrete vorticity, the two components of the discrete velocity, and the discrete pressure for the data f and $\mathbf{u}_{0}$ defined in (42), with g given by

$$ g(x,-1) = - t\bigl(1-x^{2}\bigr)^{\frac{3}{2}},\qquad g(x,1) = t\bigl(1-x^{2}\bigr)^{\frac{3}{2}},\qquad g(\pm 1,y) = 0, $$

(43)

and $N = 30$.

Figure 4

The solution $({{\tau }}, u_{x}, u_{y},p)$ for the data $(\mathbf{f},g)$ defined in (42)–(43)

We mention that the discrete vorticity ${{\tau }}_{N}$ and the discrete pressure $p_{N}$ are roughly the same in Figs. 3 and 4.

6.2 Three-dimensional experiments

We now work in the cube $\varOmega =\, ]-1,1[^{3}$ with $g = 0$. We consider a less regular solution constructed using on the formulas $\mathbf{u}= {\mathbf{curl}}\, \phi$ and ${{\tau }}= {\mathbf{curl}}\, \mathbf{u}$ with $\phi = (\phi _{x},\phi _{y},\phi _{z})$ and p defined by

$$ \begin{aligned} &\phi _{x}(x,y,z) = t \bigl(1-y^{2}\bigr)^{3}\bigl(1-z^{2} \bigr)^{\frac{7}{2}},\qquad \phi _{y}(x,y,z) = t\bigl(1-x^{2} \bigr)^{\frac{7}{2}}\bigl(1-z^{2}\bigr)^{3}, \\ &\phi _{z}(x,y,z) = t\bigl(1-x^{2} \bigr)^{3}\bigl(1-y^{2}\bigr)^{\frac{7}{2}},\qquad p(x,y,z) = e^{-t}\frac{x(1-x^{2})^{\frac{3}{2}} }{(1+y^{2})^{\frac{1}{2}}(1+z^{2})^{\frac{1 }{2}}}. \end{aligned} $$

(44)

In Fig. 5, we deal with the spectral convergence curves for the discrete solution computed from (44). The error is for the vorticity, velocity and pressure ($\log (\vert \mathrm{error} \vert )$ with respect to $\log (N)$ for N varying from 5 to 18), when $h=0.0001$. We notice that the convergence slopes of the error are similar to those in Fig. 2.

Figure 5

Error curves for the solution defined by (44)

7 Conclusion

This work concerns the numerical implementation of the implicit Euler scheme in time and the spectral discretization in space of the nonstationary vorticity–velocity–pressure formulation of the Stokes problem. We present clearly the details of the linear matrix system and the algorithm of its resolution. Some numerical tests are presented which confirm the optimality error estimates for the three unknowns (vorticity, velocity, and pressure). This estimation depends only on the regularity of the solution. Our forthcoming work concerns the nonlinear nonstationary vorticity–velocity–pressure formulation of the Navier–Stokes problem.