# Discrete harmonics for the Stokes problem

François Dubois, Michel Salaün, Stéphanie Salmon

► To cite this version:

François Dubois, Michel Salaün, Stéphanie Salmon. Discrete harmonics for the Stokes problem. Comptes Rendus. Mathématique, 2000, 331 (10), pp.827-832. 10.1016/s0764-4442(00)01726-2. hal-03181386

HAL Id: hal-03181386

https://hal.science/hal-03181386

Submitted on 27 Feb 2024

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
Discrete harmonics for the Stokes problem

François Dubois, Michel Salaün, Stéphanie Salmon

► To cite this version:

François Dubois, Michel Salaün, Stéphanie Salmon. Discrete harmonics for the Stokes problem. Comptes Rendus. Mathématique, 2000, 331 (10), pp.827-832. 10.1016/s0764-4442(00)01726-2. hal-03181386

HAL Id: hal-03181386

https://hal.science/hal-03181386v1

Submitted on 27 Feb 2024

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# Discrete harmonics for the Stokes problem

François DUBOIS, Michel SALAÜN, Stéphanie SALMON

Conservatoire national des arts et métiers, UPRES JE n° 2140, 15, rue Marat, 78210 Saint-Cyr-l'École, France
E-mail: fdubois@cnam.fr, salmons@cnam.fr, salaun@iat.cnam.fr

(Reçu le 4 avril 2000, accepté après révision le 20 septembre 2000)

**Abstract.** We consider the bidimensional Stokes problem for incompressible fluids in stream function-vorticity. For this problem, the classical finite elements method of degree one converges as $O(\sqrt{h_T})$ in the $L^2$-norm of the vorticity. We propose to use harmonic functions to approach the vorticity along the boundary. Discrete harmonics are functions that are used in practice to derive a new numerical method. We obtain with this numerical scheme an error of order $O(h_T^{1-e})$, $e > 0$, for the $L^2$-norm of the vorticity and present a first numerical test. © 2000 Académie des sciences/Éditions scientifiques et médicales Elsevier SAS

## Harmoniques discrètes pour le problème de Stokes

**Résumé.** On considère le problème de Stokes bidimensionnel pour les fluides incompressibles en formulation fonction courant-tourbillon. La méthode des éléments finis classiques de degré un donne une convergence très lente pour le tourbillon avec de nombreuses oscillations. Nous proposons d'utiliser des fonctions harmoniques pour approcher le tourbillon au bord du domaine et les harmoniques discrètes sont les fonctions utilisées en pratique pour construire une nouvelle méthode numérique. Nous obtenons avec ce schéma numérique une erreur en moyenne quadratique sur le tourbillon d'ordre arbitrairement proche de un par rapport au pas du maillage et nous présentons un premier test numérique. © 2000 Académie des sciences/Éditions scientifiques et médicales Elsevier SAS

## Version française abrégée

La formulation variationnelle du problème de Stokes en fonction courant-tourbillon [3] consiste à chercher $\omega$ dans $M(\Omega) = \{\varphi \in L^2(\Omega); \Delta\varphi \in H^{-1}(\Omega)\}$ et $\psi$ dans $H_0^1(\Omega)$ satisfaisant le problème (1). Comme la norme $L^2$ est équivalente à la norme dans l'espace $M$ sur $\mathcal{H}(\Omega) = \{\varphi \in M(\Omega); \Delta\varphi = 0\}$, les résultats classiques sur les problèmes mixtes permettent de montrer que ce problème est bien posé [3]. Pour le schéma par éléments finis de degré un (proposé dans [7]), les estimations d'erreur ne sont pas optimales (voir [6]).

La méthode numérique développée ici consiste à discrétiser directement les deux espaces intervenant dans la décomposition suivante de $M(\Omega) : M(\Omega) = H_0^1(\Omega) \oplus \mathcal{H}(\Omega)$. Nous proposons une formulation variationnelle discrète qui consiste à chercher $\psi_T \in H_{0,T}^1$ et $\omega_T \in H_T^{1,\infty} = H_{0,T}^1 \oplus \mathcal{H}_{T,\infty}$ satisfaisant le

<footnote>Note présentée par Olivier PIRONNEAU.</footnote>
# Discrete harmonics for the Stokes problem

François DUBOIS, Michel SALAÜN, Stéphanie SALMON

Conservatoire national des arts et métiers, UPRES JE n° 2140, 15, rue Marat, 78210 Saint-Cyr-l'École, France
E-mail: fdubois@cnam.fr, salmons@cnam.fr, salaun@iat.cnam.fr

(Reçu le 4 avril 2000, accepté après révision le 20 septembre 2000)

**Abstract.** We consider the bidimensional Stokes problem for incompressible fluids in stream function-vorticity. For this problem, the classical finite elements method of degree one converges as $O(\sqrt{h_T})$ in the $L^2$-norm of the vorticity. We propose to use harmonic functions to approach the vorticity along the boundary. Discrete harmonics are functions that are used in practice to derive a new numerical method. We obtain with this numerical scheme an error of order $O(h_T^{1-e})$, $e > 0$, for the $L^2$-norm of the vorticity and present a first numerical test. © 2000 Académie des sciences/Éditions scientifiques et médicales Elsevier SAS

## Harmoniques discrètes pour le problème de Stokes

**Résumé.** On considère le problème de Stokes bidimensionnel pour les fluides incompressibles en formulation fonction courant-tourbillon. La méthode des éléments finis classiques de degré un donne une convergence très lente pour le tourbillon avec de nombreuses oscillations. Nous proposons d'utiliser des fonctions harmoniques pour approcher le tourbillon au bord du domaine et les harmoniques discrètes sont les fonctions utilisées en pratique pour construire une nouvelle méthode numérique. Nous obtenons avec ce schéma numérique une erreur en moyenne quadratique sur le tourbillon d'ordre arbitrairement proche de un par rapport au pas du maillage et nous présentons un premier test numérique. © 2000 Académie des sciences/Éditions scientifiques et médicales Elsevier SAS

## Version française abrégée

La formulation variationnelle du problème de Stokes en fonction courant-tourbillon [3] consiste à chercher $\omega$ dans $M(\Omega) = \{\varphi \in L^2(\Omega); \Delta\varphi \in H^{-1}(\Omega)\}$ et $\psi$ dans $H_0^1(\Omega)$ satisfaisant le problème (1). Comme la norme $L^2$ est équivalente à la norme dans l'espace $M$ sur $\mathcal{H}(\Omega) = \{\varphi \in M(\Omega); \Delta\varphi = 0\}$, les résultats classiques sur les problèmes mixtes permettent de montrer que ce problème est bien posé [3]. Pour le schéma par éléments finis de degré un (proposé dans [7]), les estimations d'erreur ne sont pas optimales (voir [6]).

La méthode numérique développée ici consiste à discrétiser directement les deux espaces intervenant dans la décomposition suivante de $M(\Omega) : M(\Omega) = H_0^1(\Omega) \oplus \mathcal{H}(\Omega)$. Nous proposons une formulation variationnelle discrète qui consiste à chercher $\psi_T \in H_{0,T}^1$ et $\omega_T \in H_T^{1,\infty} = H_{0,T}^1 \oplus \mathcal{H}_{T,\infty}$ satisfaisant le

<footnote>Note présentée par Olivier PIRONNEAU.</footnote>
système (3), où $H_{0,\mathcal{T}}^1$ est l'espace des fonctions continues sur $\bar{\Omega}$, nulles sur le bord, polynomiales de degré un dans chaque triangle, et $\mathcal{H}_{\mathcal{T},\infty}$ défini en (2) est un espace de fonctions harmoniques à trace polynomiale de degré un sur la frontière. La méthode proposée est une discrétisation naturelle du problème précédent puisque $H_{\mathcal{T}}^{1,\infty}$ est un sous-espace de $M(\Omega)$. On obtient les résultats suivants d'existence, d'unicité et de convergence.

**PROPOSITION 1** (Existence et unicité d'une solution discrète). – Si $f \in (L^2(\Omega))^2$, le problème de Stokes discret en formulation fonction courant-tourbillon (3) a une unique solution $(\psi_{\mathcal{T}}, \omega_{\mathcal{T}}) \in H_{0,\mathcal{T}}^1 \times H_{\mathcal{T}}^{1,\infty}$ qui dépend continûment de la donnée $f$.

**THÉORÈME 1** (Convergence). – Si le maillage $\mathcal{T}$ appartient à une famille régulière de triangulations et $\text{si } (\psi, \omega) \in (H^2(\Omega) \cap H_0^1(\Omega)) \times H^2(\Omega)$, il existe une constante $C > 0$ indépendante de $\mathcal{T}$ telle que l'unique solution $(\psi_{\mathcal{T}}, \omega_{\mathcal{T}})$ du problème (3) vérifie : $\|\omega - \omega_{\mathcal{T}}\|_M + \|\psi - \psi_{\mathcal{T}}\|_{1,\Omega} \le C h_{\mathcal{T}} (\|\omega\|_{2,\Omega} + \|\psi\|_{2,\Omega})$.

Le problème est ensuite d'approcher numériquement l'espace $\mathcal{H}_{\mathcal{T},\infty}$, composé de fonctions harmoniques. On utilise pour cela des éléments finis défini sur la famille de sous-maillages $(\mathcal{T}_k)_{k \in \mathbb{N}}$ d'un maillage donné $\mathcal{T}$ et définie ainsi : (i) $\mathcal{T}_0 = \mathcal{T}$, (ii) $\mathcal{T}_{k+1} \subset \mathcal{T}_k$ est obtenu à partir de $\mathcal{T}_k$ en divisant chaque triangle en quatre triangles homothétiques. On note $H_{\mathcal{T}_k}^1$, l'espace des fonctions continues sur $\bar{\Omega}$, polynomiales de degré 1 dans chaque triangle de $\mathcal{T}_k$ et $H_{0,\mathcal{T}_k}^1 = H_{\mathcal{T}_k}^1 \cap H_0^1(\Omega)$. L'espace des (fonctions) harmoniques discrètes $\mathcal{H}_{\mathcal{T},k}$ est défini de la façon suivante : d'une part les valeurs au bord de $\varphi \in \mathcal{H}_{\mathcal{T},k}$ définissent une fonction linéaire par morceaux sur le maillage initial $\mathcal{T}_0 = \mathcal{T}$ et d'autre part les valeurs aux noeuds internes d'une harmonique discrète sont obtenues par résolution d'un problème de Laplace discret dans l'espace $H_{0,\mathcal{T}_k}^1$ : $(\nabla\varphi, \nabla\xi) = 0$ pour tout $\xi \in H_{0,\mathcal{T}_k}^1$. L'espace $\mathcal{H}_{\mathcal{T},k}$ est de dimension finie égale au nombre de sommets du maillage initial $\mathcal{T}$ sur la frontière $\Gamma$. La formulation harmonique discrétisée du problème de Stokes consiste alors à trouver $\psi_{\mathcal{T},k} \in H_{0,\mathcal{T}}^1$, $\omega_{\mathcal{T},k} \in H_{\mathcal{T}}^{1,k} = H_{0,\mathcal{T}}^1 \oplus \mathcal{H}_{\mathcal{T},k}$ solutions du problème (4). Le théorème suivant traduit le fait que l'utilisation des harmoniques discrètes entraîne la stabilité du schéma discret si on a suffisamment découpé le maillage initial.

**THÉORÈME 2** (Stabilité des harmoniques discrètes). – Pour tout maillage régulier $\mathcal{T}$, il existe $K(\mathcal{T}) \in \mathbb{N}$, tel que :

$$ \forall k \ge K(\mathcal{T}), \quad \forall \varphi \in \mathcal{H}_{\mathcal{T},k}, \quad \|\varphi\|_{0,\Omega} \ge \|\Delta\varphi\|_{-1,\Omega}. $$

**THÉORÈME 3** (Convergence). – Le problème (4) a une unique solution $(\psi_{\mathcal{T},k}, \omega_{\mathcal{T},k}) \in H_{0,\mathcal{T}}^1 \times H_{\mathcal{T}}^{1,k}$. Si $(\psi, \omega) \in (H^2(\Omega) \cap H_0^1(\Omega)) \times H^2(\Omega)$ est solution du problème continu, on a :

$$ \forall \theta, 0 < \theta < 1, \exists C(\theta) > 0, \forall \mathcal{T} \text{ régulière, } \exists K(\mathcal{T}), \forall k \ge K(\mathcal{T}), \\ \|\omega - \omega_{\mathcal{T},k}\|_M + \|\psi - \psi_{\mathcal{T},k}\|_{1,\Omega} \le C h_{\mathcal{T}}^{1-\theta} (\|\omega\|_{2,\Omega} + \|\psi\|_{2,\Omega}). $$

Enfin, les résultats obtenus sur un cas test proposé par Bercovier et Engelman [2] montrent que la difficulté numérique consiste à contrôler les oscillations du tourbillon discret sur la frontière. La figure 1 présente les oscillations initiales [7] qui disparaissent au fur et à mesure que le niveau $k$ de raffinement augmente. Sur trois maillages différents (notés A, B, C dans le tableau 1), nous constatons que la valeur $K(\mathcal{T}) = 4$ conduit en pratique à un calcul stable. Avec ce choix de $K(\mathcal{T})$, la fonction courant et le tourbillon convergent (courbe en figure 2) avec un ordre en accord avec le théorème 3 (et même meilleur !).

# 1. Introduction

Let $\Omega$ be an open bounded domain of $\mathbb{R}^2$, convex with boundary $\Gamma$ assumed to be polygonal. Let also $f$ be a field of given forces, belonging to $(L^2(\Omega))^2$, whose standard inner product is denoted by $(\cdot, \cdot)$. The well-posed weak formulation [3] of the stationary and homogeneous Stokes problem needs the
système (3), où $H_{0,\mathcal{T}}^1$ est l'espace des fonctions continues sur $\bar{\Omega}$, nulles sur le bord, polynomiales de degré un dans chaque triangle, et $\mathcal{H}_{\mathcal{T},\infty}$ défini en (2) est un espace de fonctions harmoniques à trace polynomiale de degré un sur la frontière. La méthode proposée est une discrétisation naturelle du problème précédent puisque $H_{\mathcal{T}}^{1,\infty}$ est un sous-espace de $M(\Omega)$. On obtient les résultats suivants d'existence, d'unicité et de convergence.

**PROPOSITION 1** (Existence et unicité d'une solution discrète). – Si $f \in (L^2(\Omega))^2$, le problème de Stokes discret en formulation fonction courant-tourbillon (3) a une unique solution $(\psi_{\mathcal{T}}, \omega_{\mathcal{T}}) \in H_{0,\mathcal{T}}^1 \times H_{\mathcal{T}}^{1,\infty}$ qui dépend continûment de la donnée $f$.

**THÉORÈME 1** (Convergence). – Si le maillage $\mathcal{T}$ appartient à une famille régulière de triangulations et $\text{si } (\psi, \omega) \in (H^2(\Omega) \cap H_0^1(\Omega)) \times H^2(\Omega)$, il existe une constante $C > 0$ indépendante de $\mathcal{T}$ telle que l'unique solution $(\psi_{\mathcal{T}}, \omega_{\mathcal{T}})$ du problème (3) vérifie : $\|\omega - \omega_{\mathcal{T}}\|_M + \|\psi - \psi_{\mathcal{T}}\|_{1,\Omega} \le C h_{\mathcal{T}} (\|\omega\|_{2,\Omega} + \|\psi\|_{2,\Omega})$.

Le problème est ensuite d'approcher numériquement l'espace $\mathcal{H}_{\mathcal{T},\infty}$, composé de fonctions harmoniques. On utilise pour cela des éléments finis défini sur la famille de sous-maillages $(\mathcal{T}_k)_{k \in \mathbb{N}}$ d'un maillage donné $\mathcal{T}$ et définie ainsi : (i) $\mathcal{T}_0 = \mathcal{T}$, (ii) $\mathcal{T}_{k+1} \subset \mathcal{T}_k$ est obtenu à partir de $\mathcal{T}_k$ en divisant chaque triangle en quatre triangles homothétiques. On note $H_{\mathcal{T}_k}^1$, l'espace des fonctions continues sur $\bar{\Omega}$, polynomiales de degré 1 dans chaque triangle de $\mathcal{T}_k$ et $H_{0,\mathcal{T}_k}^1 = H_{\mathcal{T}_k}^1 \cap H_0^1(\Omega)$. L'espace des (fonctions) harmoniques discrètes $\mathcal{H}_{\mathcal{T},k}$ est défini de la façon suivante : d'une part les valeurs au bord de $\varphi \in \mathcal{H}_{\mathcal{T},k}$ définissent une fonction linéaire par morceaux sur le maillage initial $\mathcal{T}_0 = \mathcal{T}$ et d'autre part les valeurs aux noeuds internes d'une harmonique discrète sont obtenues par résolution d'un problème de Laplace discret dans l'espace $H_{0,\mathcal{T}_k}^1$ : $(\nabla\varphi, \nabla\xi) = 0$ pour tout $\xi \in H_{0,\mathcal{T}_k}^1$. L'espace $\mathcal{H}_{\mathcal{T},k}$ est de dimension finie égale au nombre de sommets du maillage initial $\mathcal{T}$ sur la frontière $\Gamma$. La formulation harmonique discrétisée du problème de Stokes consiste alors à trouver $\psi_{\mathcal{T},k} \in H_{0,\mathcal{T}}^1$, $\omega_{\mathcal{T},k} \in H_{\mathcal{T}}^{1,k} = H_{0,\mathcal{T}}^1 \oplus \mathcal{H}_{\mathcal{T},k}$ solutions du problème (4). Le théorème suivant traduit le fait que l'utilisation des harmoniques discrètes entraîne la stabilité du schéma discret si on a suffisamment découpé le maillage initial.

**THÉORÈME 2** (Stabilité des harmoniques discrètes). – Pour tout maillage régulier $\mathcal{T}$, il existe $K(\mathcal{T}) \in \mathbb{N}$, tel que :

$$ \forall k \ge K(\mathcal{T}), \forall \varphi \in \mathcal{H}_{\mathcal{T},k}, \quad \|\varphi\|_{0,\Omega} \ge \|\Delta\varphi\|_{-1,\Omega}. $$

**THÉORÈME 3** (Convergence). – Le problème (4) a une unique solution $(\psi_{\mathcal{T},k}, \omega_{\mathcal{T},k}) \in H_{0,\mathcal{T}}^1 \times H_{\mathcal{T}}^{1,k}$. Si $(\psi, \omega) \in (H^2(\Omega) \cap H_0^1(\Omega)) \times H^2(\Omega)$ est solution du problème continu, on a :

$$ \forall \theta, 0 < \theta < 1, \exists C(\theta) > 0, \forall \mathcal{T} \text{ régulière, } \exists K(\mathcal{T}), \forall k \ge K(\mathcal{T}), \\ \|\omega - \omega_{\mathcal{T},k}\|_M + \|\psi - \psi_{\mathcal{T},k}\|_{1,\Omega} \le C h_{\mathcal{T}}^{1-\theta} (\|\omega\|_{2,\Omega} + \|\psi\|_{2,\Omega}). $$

Enfin, les résultats obtenus sur un cas test proposé par Bercovier et Engelman [2] montrent que la difficulté numérique consiste à contrôler les oscillations du tourbillon discret sur la frontière. La figure 1 présente les oscillations initiales [7] qui disparaissent au fur et à mesure que le niveau $k$ de raffinement augmente. Sur trois maillages différents (notés A, B, C dans le tableau 1), nous constatons que la valeur $K(\mathcal{T}) = 4$ conduit en pratique à un calcul stable. Avec ce choix de $K(\mathcal{T})$, la fonction courant et le tourbillon convergent (courbe en figure 2) avec un ordre en accord avec le théorème 3 (et même meilleur !).

# 1. Introduction

Let $\Omega$ be an open bounded domain of $\mathbb{R}^2$, convex with boundary $\Gamma$ assumed to be polygonal. Let also $f$ be a field of given forces, belonging to $(L^2(\Omega))^2$, whose standard inner product is denoted by $(\cdot, \cdot)$. The well-posed weak formulation [3] of the stationary and homogeneous Stokes problem needs the
introduction of the classical functional space $H_0^1(\Omega)$ (defined in, e.g., [6]), its dual $H^{-1}(\Omega)$ ($\langle \cdot, \cdot \rangle_{-1,1}$ shall denote the duality product between $H^{-1}(\Omega)$ and $H_0^1(\Omega)$) and the space $M(\Omega)$ defined by: $M(\Omega) = \{\varphi \in L^2(\Omega); \Delta\varphi \in H^{-1}(\Omega)\}$. Let us remark with [3] that space $M(\Omega)$ is a Hilbert space for the norm $\|\varphi\|_M = \sqrt{\|\varphi\|_{0,\Omega}^2 + \|\Delta\varphi\|_{-1,\Omega}^2}$ and that $H^1(\Omega)$ is contained in $M(\Omega)$ with continuous imbedding. We introduce the left kernel $\mathcal{H}(\Omega)$ of the bilinear form $\langle \Delta\cdot, \cdot \rangle_{-1,1}: \mathcal{H}(\Omega) = \{\varphi \in M(\Omega), \Delta\varphi = 0\}$. It is classical that there exists a trace operator from $\mathcal{H}(\Omega)$ into $H^{-1/2}(\Gamma)$ (denoted by $\gamma$ like the trace operator from $H^1(\Omega)$). Moreover, we recall the decomposition: $M(\Omega) = H_0^1(\Omega) \oplus \mathcal{H}(\Omega)$. This decomposition is established by splitting $\varphi \in M(\Omega)$ into two parts: $M(\Omega) \ni \varphi \mapsto (\varphi^0, \varphi^\Delta) \in H_0^1(\Omega) \times \mathcal{H}(\Omega)$. The component $\varphi^0$ is the unique variational solution in $H_0^1(\Omega)$ of the Dirichlet problem $\Delta\varphi^0 = \Delta\varphi$ in $\Omega$, $\gamma\varphi^0 = 0$ on $\Gamma$ and $\varphi^\Delta = \varphi - \varphi^0$ is a harmonic function.

The variational formulation of the Stokes problem in stream function-vorticity consists in finding $\psi$ in $H_0^1(\Omega)$ and $\omega$ in $M(\Omega)$ such that

$$ \left\{ \begin{array}{ll} (\omega, \varphi) + \langle \Delta\varphi, \psi \rangle_{-1,1} = 0, & \forall \varphi \in M(\Omega), \\ \langle -\Delta\omega, \xi \rangle_{-1,1} = (\mathbf{f}, \mathrm{rot}\,\xi), & \forall \xi \in H_0^1(\Omega). \end{array} \right. \qquad (1) $$

As the $L^2$-norm is equivalent to the $M$-norm on $\mathcal{H}(\Omega)$, classical results on mixed problems show that problem (1) is well posed [3]. Error estimates derived for the linear finite element scheme used to approximate (1) are not optimal. If $h_T$ denotes the maximum diameter of the elements of a triangulation, the bound is in $h_T^{1/2}$ for the $L^2$-norm of the vorticity [6]. It is possible to stabilize the numerical scheme by adding jumps at interfaces and prove convergence [1], but the numerical method, developed further, consists in a direct discretization of the above-mentioned decomposition of $M(\Omega)$.

## 2. Harmonic discretization of the Stokes problem

We introduce a mesh $\mathcal{T}$ that belongs to a regular family of triangulations in the sense of [4]. The set $H_T^1$ denotes the space of continuous functions defined on $\bar{\Omega}$, polynomial of degree 1 in each triangle of $\mathcal{T}$ and $H_{0,\mathcal{T}}^1 = H_T^1 \cap H_0^1(\Omega)$. We define the space $\Gamma_T = \{\lambda$ continuous on $\Gamma$, linear on each edge of the boundary mesh $\mathcal{T} \cap \Gamma\}$ of all traces of functions in $H_T^1$. Then dimension of $\Gamma_T$ is equal to the number of vertices of $\mathcal{T}$ on the boundary $\Gamma$. We introduce the space $\mathcal{H}_{\mathcal{T},\infty}$ of harmonic functions whose trace on $\Gamma$ belongs to discrete space $\Gamma_T$:

$$ \mathcal{H}_{\mathcal{T},\infty} = \{\varphi \in H^1(\Omega), \exists \lambda \in \Gamma_T, \gamma\varphi = \lambda, \Delta\varphi = 0 \text{ in } H^{-1}(\Omega)\}. \qquad (2) $$

Space $\mathcal{H}_{\mathcal{T},\infty}$ is finite-dimensional and $\dim \mathcal{H}_{\mathcal{T},\infty} = \dim \Gamma_T$. We set: $H_T^{1,\infty} = H_{0,\mathcal{T}}^1 \oplus \mathcal{H}_{\mathcal{T},\infty}$. Dimension of $H_T^{1,\infty}$ is the same as $H_T^1$ but near the boundary, we use harmonic functions instead of piecewise linear continuous functions. Values of the harmonic functions decrease very quickly. Nevertheless, support is spread over the entire domain. The splitting $\varphi = \varphi^0 + \varphi^\Delta$ proposed in the introduction remains valid for functions in $H_T^{1,\infty}$. Indeed, as $\varphi \in H_T^{1,\infty}$, $\gamma\varphi \in \Gamma_T$, then $\varphi^\Delta \in \mathcal{H}_{\mathcal{T},\infty}$. By construction, $\varphi^0 = \varphi - \varphi^\Delta$ belongs to space $H_{0,\mathcal{T}}^1$.

We propose a natural discrete variational formulation of problem (1). Space $H_T^{1,\infty}$ is a subset of $M(\Omega)$ and is contained in $H^1(\Omega)$, so the following algebra holds: $\langle \Delta\varphi, \xi \rangle_{-1,1} = (\nabla\varphi, \nabla\xi)$ when $\varphi \in H_T^{1,\infty}$ and $\xi \in H_{0,\mathcal{T}}^1$. The discrete problem consists in finding $\psi_T \in H_{0,\mathcal{T}}^1$ and $\omega_T \in H_T^{1,\infty}$ such that:

$$ \left\{ \begin{array}{ll} (\omega_T, \varphi) - (\nabla\psi_T, \nabla\varphi) = 0, & \forall \varphi \in H_T^{1,\infty}, \\ (\nabla\omega_T, \nabla\xi) = (\mathbf{f}, \mathrm{rot}\,\xi), & \forall \xi \in H_{0,\mathcal{T}}^1. \end{array} \right. \qquad (3) $$
introduction of the classical functional space $H_0^1(\Omega)$ (defined in, e.g., [6]), its dual $H^{-1}(\Omega)$ ($\langle \cdot, \cdot \rangle_{-1,1}$ shall denote the duality product between $H^{-1}(\Omega)$ and $H_0^1(\Omega)$) and the space $M(\Omega)$ defined by: $M(\Omega) = \{\varphi \in L^2(\Omega); \Delta\varphi \in H^{-1}(\Omega)\}$. Let us remark with [3] that space $M(\Omega)$ is a Hilbert space for the norm $\|\varphi\|_M = \sqrt{\|\varphi\|_{0,\Omega}^2 + \|\Delta\varphi\|_{-1,\Omega}^2}$ and that $H^1(\Omega)$ is contained in $M(\Omega)$ with continuous imbedding. We introduce the left kernel $\mathcal{H}(\Omega)$ of the bilinear form $\langle \Delta\cdot, \cdot \rangle_{-1,1}: \mathcal{H}(\Omega) = \{\varphi \in M(\Omega), \Delta\varphi = 0\}$. It is classical that there exists a trace operator from $\mathcal{H}(\Omega)$ into $H^{-1/2}(\Gamma)$ (denoted by $\gamma$ like the trace operator from $H^1(\Omega)$). Moreover, we recall the decomposition: $M(\Omega) = H_0^1(\Omega) \oplus \mathcal{H}(\Omega)$. This decomposition is established by splitting $\varphi \in M(\Omega)$ into two parts: $M(\Omega) \ni \varphi \mapsto (\varphi^0, \varphi^\Delta) \in H_0^1(\Omega) \times \mathcal{H}(\Omega)$. The component $\varphi^0$ is the unique variational solution in $H_0^1(\Omega)$ of the Dirichlet problem $\Delta\varphi^0 = \Delta\varphi$ in $\Omega$, $\gamma\varphi^0 = 0$ on $\Gamma$ and $\varphi^\Delta = \varphi - \varphi^0$ is a harmonic function.

The variational formulation of the Stokes problem in stream function-vorticity consists in finding $\psi$ in $H_0^1(\Omega)$ and $\omega$ in $M(\Omega)$ such that

$$ \left\{ \begin{array}{ll} (\omega, \varphi) + \langle \Delta\varphi, \psi \rangle_{-1,1} = 0, & \forall \varphi \in M(\Omega), \\ \langle -\Delta\omega, \xi \rangle_{-1,1} = (\mathbf{f}, \mathrm{rot}\,\xi), & \forall \xi \in H_0^1(\Omega). \end{array} \right. \qquad (1) $$

As the $L^2$-norm is equivalent to the $M$-norm on $\mathcal{H}(\Omega)$, classical results on mixed problems show that problem (1) is well posed [3]. Error estimates derived for the linear finite element scheme used to approximate (1) are not optimal. If $h_T$ denotes the maximum diameter of the elements of a triangulation, the bound is in $h_T^{1/2}$ for the $L^2$-norm of the vorticity [6]. It is possible to stabilize the numerical scheme by adding jumps at interfaces and prove convergence [1], but the numerical method, developed further, consists in a direct discretization of the above-mentioned decomposition of $M(\Omega)$.

## 2. Harmonic discretization of the Stokes problem

We introduce a mesh $\mathcal{T}$ that belongs to a regular family of triangulations in the sense of [4]. The set $H_T^1$ denotes the space of continuous functions defined on $\bar{\Omega}$, polynomial of degree 1 in each triangle of $\mathcal{T}$ and $H_{0,\mathcal{T}}^1 = H_T^1 \cap H_0^1(\Omega)$. We define the space $\Gamma_T = \{\lambda$ continuous on $\Gamma$, linear on each edge of the boundary mesh $\mathcal{T} \cap \Gamma\}$ of all traces of functions in $H_T^1$. Then dimension of $\Gamma_T$ is equal to the number of vertices of $\mathcal{T}$ on the boundary $\Gamma$. We introduce the space $\mathcal{H}_{\mathcal{T},\infty}$ of harmonic functions whose trace on $\Gamma$ belongs to discrete space $\Gamma_T$:

$$ \mathcal{H}_{\mathcal{T},\infty} = \{\varphi \in H^1(\Omega), \exists \lambda \in \Gamma_T, \gamma\varphi = \lambda, \Delta\varphi = 0 \text{ in } H^{-1}(\Omega)\}. \qquad (2) $$

Space $\mathcal{H}_{\mathcal{T},\infty}$ is finite-dimensional and $\dim \mathcal{H}_{\mathcal{T},\infty} = \dim \Gamma_T$. We set: $H_T^{1,\infty} = H_{0,\mathcal{T}}^1 \oplus \mathcal{H}_{\mathcal{T},\infty}$. Dimension of $H_T^{1,\infty}$ is the same as $H_T^1$ but near the boundary, we use harmonic functions instead of piecewise linear continuous functions. Values of the harmonic functions decrease very quickly. Nevertheless, support is spread over the entire domain. The splitting $\varphi = \varphi^0 + \varphi^\Delta$ proposed in the introduction remains valid for functions in $H_T^{1,\infty}$. Indeed, as $\varphi \in H_T^{1,\infty}$, $\gamma\varphi \in \Gamma_T$, then $\varphi^\Delta \in \mathcal{H}_{\mathcal{T},\infty}$. By construction, $\varphi^0 = \varphi - \varphi^\Delta$ belongs to space $H_{0,\mathcal{T}}^1$.

We propose a natural discrete variational formulation of problem (1). Space $H_T^{1,\infty}$ is a subset of $M(\Omega)$ and is contained in $H^1(\Omega)$, so the following algebra holds: $\langle \Delta\varphi, \xi \rangle_{-1,1} = (\nabla\varphi, \nabla\xi)$ when $\varphi \in H_T^{1,\infty}$ and $\xi \in H_{0,\mathcal{T}}^1$. The discrete problem consists in finding $\psi_T \in H_{0,\mathcal{T}}^1$ and $\omega_T \in H_T^{1,\infty}$ such that:

$$ \left\{ \begin{array}{ll} (\omega_T, \varphi) - (\nabla\psi_T, \nabla\varphi) = 0, & \forall \varphi \in H_T^{1,\infty}, \\ (\nabla\omega_T, \nabla\xi) = (\mathbf{f}, \mathrm{rot}\,\xi), & \forall \xi \in H_{0,\mathcal{T}}^1. \end{array} \right. \qquad (3) $$
**PROPOSITION 1 (Existence and uniqueness of a discrete solution).** – If $\mathbf{f} \in (\mathbb{L}^2(\Omega))^2$, problem (3) has a unique solution $(\psi_T, \omega_T) \in H_{0,T}^1 \times H_T^{1,\infty}$ which depends continuously on the datum $\mathbf{f}$: $\exists C > 0$ independent on $T$, $\|\omega_T\|_M + \|\nabla \psi_T\|_{0,\Omega} \le C\|\mathbf{f}\|_{0,\Omega}$.

Although $\mathcal{H}_{T,\infty}$ is generated by a finite number of harmonic functions, we proceed as usual for studying finite element discretization. We first denote by $\Pi_T^\Gamma: H^1(\Gamma) \to \Gamma_T$ the Lagrange interpolation operator associated with mesh $T$ on the boundary. We define the interpolation operator $\phi_T^\infty: \mathcal{H}(\Omega) \cap H^2(\Omega) \to \mathcal{H}_{T,\infty}$ by $\phi_T^\infty \varphi^\Delta = \zeta$, where $\zeta \in H^1(\Omega)$ is the solution of the Dirichlet problem $\Delta\zeta = 0$ in $\Omega$ and $\gamma\zeta = \Pi_T^\Gamma(\gamma\varphi^\Delta)$ on $\Gamma$. We consider also the usual Lagrange interpolation operator associated with mesh $T$ in the domain $\Omega$, denoted by $\Pi_T^\Omega: H^2(\Omega) \to H_T^1$. The interpolation operator $\mathcal{P}_T^\infty$ is now defined by the relations:

$$ \mathcal{P}_T^\infty: M(\Omega) \cap H^2(\Omega) = H^2(\Omega) \ni \varphi = \varphi^0 + \varphi^\Delta \mapsto \mathcal{P}_T^\infty \varphi = \Pi_T^\Omega \varphi^0 + \phi_T^\infty \varphi^\Delta \in H_{0,T}^1 \oplus \mathcal{H}_{T,\infty} = H_T^{1,\infty}. $$

We have (see [5]) the following estimates for every $\varphi \in H^2(\Omega) \cap M(\Omega): \|\varphi^0 - \Pi_T^\Omega \varphi^0\|_M \le C h_T \|\varphi^0\|_{2,\Omega}$ and $\|\varphi^\Delta - \phi_T^\infty \varphi^\Delta\|_M \le C h_T \|\varphi\|_{2,\Omega}$ which lead to the:

**THEOREM 1 (Convergence result).** – Let $(\psi, \omega)$ be the solution of the Stokes problem (1). We suppose that it satisfies the following hypothesis of regularity: $\psi \in (H^2(\Omega) \cap H_0^1(\Omega))$, $\omega \in H^2(\Omega)$. Consider $(\psi_T, \omega_T) \in H_{0,T}^1 \times H_T^{1,\infty}$ the solution of problem (3) associated with mesh $T$. If $T$ belongs to a regular family of triangulations, there exists a constant $C > 0$ independent on $T$ such that:

$$ \begin{aligned} \|\omega - \omega_T\|_M + \|\psi - \psi_T\|_{1,\Omega} &\le C(\|\omega^0 - \Pi_T^\Omega \omega^0\|_M + \|\omega^\Delta - \phi_T^\infty \omega^\Delta\|_M + \|\psi - \Pi_T^\Omega \psi\|_{1,\Omega}) \\ &\le C h_T (\|\omega\|_{2,\Omega} + \|\psi\|_{2,\Omega}). \end{aligned} $$

### 3. Discrete harmonics

We focus now on the fact that discrete space $\mathcal{H}_{T,\infty}$ is composed of rigorously harmonic functions and has to be numerically approached. We use in the following finite elements on sub-grids defined from a given mesh $T$. The family of grids $T_k (k \in \mathbb{N})$ is defined recursively as follows: (i) $\mathcal{T}_0 = T$ and (ii) $\mathcal{T}_{k+1} \subset \mathcal{T}_k$ is obtained from $\mathcal{T}_k$ by dividing each triangle into four homothetic ones. We denote $H_{T_k}^1$ the space of all continuous functions on $\bar{\Omega}$, polynomial of degree 1 in each triangle of $\mathcal{T}_k$ and $H_{0,\mathcal{T}_k}^1 = H_{T_k}^1 \cap H_0^1(\Omega)$. We define $\mathcal{H}_{T,k}$ the space of discrete harmonic functions by: $\mathcal{H}_{T,k} = \{\varphi \in H_{T_k}^1, \exists \lambda \in \Gamma_T, \gamma\varphi = \lambda \text{ on } \Gamma, (\nabla\varphi, \nabla\xi) = 0, \forall \xi \in H_{0,\mathcal{T}_k}^1\}$. Internal values of a discrete harmonic are obtained by solving a discrete Laplace equation in space $H_{0,\mathcal{T}_k}^1$ while trace on the boundary $\mathcal{T}_k \cap \Gamma$ is a piecewise linear function on the initial mesh $\mathcal{T}_0 = T$. Note that discrete space $\mathcal{H}_{T,0}$ is the one used initially by Glowinski and Pironneau [7]. Space $\mathcal{H}_{T,k}$ is finite-dimensional, and dim $\mathcal{H}_{T,k} = \dim \Gamma_T$ is equal to the number of vertices on $\Gamma$. If we refine enough the initial mesh, discrete harmonics lead to a stable numerical scheme.

**THEOREM 2 (Stability of discrete harmonics).** – For any mesh $T$ belonging to a regular family of triangulations, there exists $K(T) \in \mathbb{N}$, such that:

$$ \forall k \ge K(T), \forall \varphi \in \mathcal{H}_{T,k}, \quad \|\varphi\|_{0,\Omega} \ge \|\Delta\varphi\|_{-1,\Omega}. $$

For every regular mesh $T$ and every integer $k \ge 0$, we set $H_{T,k}^{1,k} = H_{0,T}^1 \oplus \mathcal{H}_{T,k}$. The discretized harmonic formulation of the Stokes problem consists now in finding $\psi_{T,k} \in H_{0,T}^1$, $\omega_{T,k} \in H_T^{1,k}$ such that:

$$ \left\{ \begin{array}{ll} (\omega_{T,k}, \varphi) - (\nabla \psi_{T,k}, \nabla \varphi) = 0, & \forall \varphi \in H_T^{1,k}, \\ (\nabla \omega_{T,k}, \nabla \xi) = (\mathbf{f}, \mathbf{rot} \xi), & \forall \xi \in H_{0,T}^1. \end{array} \right. \qquad (4) $$
**PROPOSITION 1** (Existence and uniqueness of a discrete solution). – If $\mathbf{f} \in (\mathbb{L}^2(\Omega))^2$, problem (3) has a unique solution $(\psi_{\mathcal{T}}, \omega_{\mathcal{T}}) \in H_{0, \mathcal{T}}^1 \times H_{\mathcal{T}}^{1, \infty}$ which depends continuously on the datum $\mathbf{f}$: $\exists C > 0$ independent on $\mathcal{T}$, $\|\omega_{\mathcal{T}}\|_{\mathcal{M}} + \|\nabla \psi_{\mathcal{T}}\|_{0, \Omega} \le C \|\mathbf{f}\|_{0, \Omega}$.

Although $\mathcal{H}_{\mathcal{T}, \infty}$ is generated by a finite number of harmonic functions, we proceed as usual for studying finite element discretization. We first denote by $\Pi_{\mathcal{T}}^{\Gamma}: \mathbb{H}^1(\Gamma) \to \Gamma_{\mathcal{T}}$ the Lagrange interpolation operator associated with mesh $\mathcal{T}$ on the boundary. We define the interpolation operator $\phi_{\mathcal{T}}^{\infty}: \mathcal{H}(\Omega) \cap \mathbb{H}^2(\Omega) \to \mathcal{H}_{\mathcal{T}, \infty}$ by $\phi_{\mathcal{T}}^{\infty} \varphi^{\Delta} = \zeta$, where $\zeta \in \mathbb{H}^1(\Omega)$ is the solution of the Dirichlet problem $\Delta\zeta = 0$ in $\Omega$ and $\gamma\zeta = \Pi_{\mathcal{T}}^{\Gamma}(\gamma\varphi^{\Delta})$ on $\Gamma$. We consider also the usual Lagrange interpolation operator associated with mesh $\mathcal{T}$ in the domain $\Omega$, denoted by $\Pi_{\mathcal{T}}^{\Omega}: \mathbb{H}^2(\Omega) \to \mathbb{H}_{\mathcal{T}}^{1}$. The interpolation operator $\mathcal{P}_{\mathcal{T}}^{\infty}$ is now defined by the relations:

$$ \mathcal{P}_{\mathcal{T}}^{\infty}: \mathcal{M}(\Omega) \cap \mathbb{H}^2(\Omega) = \mathbb{H}^2(\Omega) \ni \varphi = \varphi^0 + \varphi^{\Delta} \mapsto \mathcal{P}_{\mathcal{T}}^{\infty} \varphi = \Pi_{\mathcal{T}}^{\Omega} \varphi^0 + \phi_{\mathcal{T}}^{\infty} \varphi^{\Delta} \in \mathbb{H}_{0, \mathcal{T}}^{1} \oplus \mathcal{H}_{\mathcal{T}, \infty} = \mathbb{H}_{\mathcal{T}}^{1, \infty}. $$

We have (see [5]) the following estimates for every $\varphi \in \mathbb{H}^2(\Omega) \cap \mathcal{M}(\Omega): \|\varphi^0 - \Pi_{\mathcal{T}}^{\Omega} \varphi^0\|_{\mathcal{M}} \le C h_{\mathcal{T}} \|\varphi^0\|_{2, \Omega}$ and $\|\varphi^{\Delta} - \phi_{\mathcal{T}}^{\infty} \varphi^{\Delta}\|_{\mathcal{M}} \le C h_{\mathcal{T}} \|\varphi\|_{2, \Omega}$ which lead to the:

**THEOREM 1** (Convergence result). – Let $(\psi, \omega)$ be the solution of the Stokes problem (1). We suppose that it satisfies the following hypothesis of regularity: $\psi \in (\mathbb{H}^2(\Omega) \cap H_0^1(\Omega))$, $\omega \in \mathbb{H}^2(\Omega)$. Consider $(\psi_{\mathcal{T}}, \omega_{\mathcal{T}}) \in \mathbb{H}_{0, \mathcal{T}}^1 \times \mathbb{H}_{\mathcal{T}}^{1, \infty}$ the solution of problem (3) associated with mesh $\mathcal{T}$. If $\mathcal{T}$ belongs to a regular family of triangulations, there exists a constant $C > 0$ independent on $\mathcal{T}$ such that:

$$ \begin{aligned} \|\omega - \omega_{\mathcal{T}}\|_{\mathcal{M}} + \|\psi - \psi_{\mathcal{T}}\|_{1, \Omega} &\le C (\|\omega^0 - \Pi_{\mathcal{T}}^{\Omega} \omega^0\|_{\mathcal{M}} + \|\omega^{\Delta} - \phi_{\mathcal{T}}^{\infty} \omega^{\Delta}\|_{\mathcal{M}} + \|\psi - \Pi_{\mathcal{T}}^{\Omega} \psi\|_{1, \Omega}) \\ &\le C h_{\mathcal{T}} (\|\omega\|_{2, \Omega} + \|\psi\|_{2, \Omega}). \end{aligned} $$

### 3. Discrete harmonics

We focus now on the fact that discrete space $\mathcal{H}_{\mathcal{T}, \infty}$ is composed of rigorously harmonic functions and has to be numerically approached. We use in the following finite elements on sub-grids defined from a given mesh $\mathcal{T}$. The family of grids $\mathcal{T}_k$ ($k \in \mathbb{N}$) is defined recursively as follows: (i) $\mathcal{T}_0 = \mathcal{T}$ and (ii) $\mathcal{T}_{k+1} \subset \mathcal{T}_k$ is obtained from $\mathcal{T}_k$ by dividing each triangle into four homothetic ones. We denote $\mathbb{H}_{\mathcal{T}_k}^1$ the space of all continuous functions on $\bar{\Omega}$, polynomial of degree 1 in each triangle of $\mathcal{T}_k$ and $\mathbb{H}_{0, \mathcal{T}_k}^1 = \mathbb{H}_{\mathcal{T}_k}^1 \cap \mathbb{H}_0^1(\Omega)$. We define $\mathcal{H}_{\mathcal{T}, k}$ the space of discrete harmonic functions by: $\mathcal{H}_{\mathcal{T}, k} = \{\varphi \in \mathbb{H}_{\mathcal{T}_k}^1, \exists \lambda \in \Gamma_{\mathcal{T}}, \gamma\varphi = \lambda \text{ on } \Gamma, (\nabla\varphi, \nabla\xi) = 0, \forall \xi \in \mathbb{H}_{0, \mathcal{T}_k}^1\}$. Internal values of a discrete harmonic are obtained by solving a discrete Laplace equation in space $\mathbb{H}_{0, \mathcal{T}_k}^1$ while trace on the boundary $\mathcal{T}_k \cap \Gamma$ is a piecewise linear function on the initial mesh $\mathcal{T}_0 = \mathcal{T}$. Note that discrete space $\mathcal{H}_{\mathcal{T}, 0}$ is the one used initially by Glowinski and Pironneau [7]. Space $\mathcal{H}_{\mathcal{T}, k}$ is finite-dimensional, and dim $\mathcal{H}_{\mathcal{T}, k} = \text{dim } \Gamma_{\mathcal{T}}$ is equal to the number of vertices on $\Gamma$. If we refine enough the initial mesh, discrete harmonics lead to a stable numerical scheme.

**THEOREM 2** (Stability of discrete harmonics). – For any mesh $\mathcal{T}$ belonging to a regular family of triangulations, there exists $K(\mathcal{T}) \in \mathbb{N}$, such that:

$$ \forall k \ge K(\mathcal{T}), \forall \varphi \in \mathcal{H}_{\mathcal{T}, k}, \quad \|\varphi\|_{0, \Omega} \ge \|\Delta\varphi\|_{-1, \Omega}. $$

For every regular mesh $\mathcal{T}$ and every integer $k \ge 0$, we set $\mathbb{H}_{\mathcal{T}}^{1, k} = \mathbb{H}_{0, \mathcal{T}}^1 \oplus \mathcal{H}_{\mathcal{T}, k}$. The discretized harmonic formulation of the Stokes problem consists now in finding $\psi_{\mathcal{T}, k} \in \mathbb{H}_{0, \mathcal{T}}^1$, $\omega_{\mathcal{T}, k} \in \mathbb{H}_{\mathcal{T}}^{1, k}$ such that:

$$ \left\{ \begin{array}{ll} (\omega_{\mathcal{T},k}, \varphi) - (\nabla \psi_{\mathcal{T},k}, \nabla \varphi) = 0, & \forall \varphi \in \mathbb{H}_{\mathcal{T}}^{1,k}, \\ (\nabla \omega_{\mathcal{T},k}, \nabla \xi) = (\mathbf{f}, \mathbf{rot} \, \xi), & \forall \xi \in \mathbb{H}_{0,\mathcal{T}}^1. \end{array} \right. \qquad (4) $$
PROPOSITION 2 (Solution of the discretized harmonic formulation of the Stokes problem). – The problem (4) has a unique solution $(\psi_{T,k}, \omega_{T,k}) \in H_{0,T}^1 \times H_T^{1,k}$. If $f \in (L^2(\Omega))^2$, if $T$ belongs to a family of regular meshes and if $k \ge K(T)$, with $K(T)$ defined in Theorem 2, the solution depends continuously on the datum $f$.

With classical arguments founded on uniform stability and interpolation theory, taking into account difficulties due to the fact that piecewise linear boundary functions in $\Gamma_T$ (“chinese hat”) belong only to the space $H^{3/2-\epsilon}(\Gamma)$, $(\epsilon > 0)$, we have proved [5] the:

THEOREM 3 (Convergence result). – If $(\psi, \omega) \in (H^2(\Omega) \cap H_0^1(\Omega)) \times H^2(\Omega)$ is the solution of problem (1) we have:

$$ \begin{gathered} \forall \theta, 0 < \theta < 1, \exists C(\theta) > 0, \forall T \text{ regular, } \exists K(T), \forall k \ge K(T), \\ \|\omega - \omega_{T,k}\|_M + \|\psi - \psi_{T,k}\|_{1,\Omega} \le C h_T^{1-\theta} (\|\omega\|_{2,\Omega} + \|\psi\|_{2,\Omega}). \end{gathered} $$

4. Numerical results

Numerical experiments have been performed on a unit square with an analytical test proposed by Bercovier and Engelman [2]. As the integer $K(T)$ a priori depends on the mesh (see Theorem 2), we have worked with different unstructured meshes (A, B, C) [5]. The uncoupled algorithm is the following:

$$ \begin{align} \text{find } \omega_T^0 &\in H_{0,T}^1 \text{ such that } (\nabla \omega_T^0, \nabla \xi) = (f, \operatorname{rot} \xi), \forall \xi \in H_{0,T}^1, \nonumber \\ \text{find } \omega_{T,k}^\Delta &\in \mathcal{H}_{T,k} \text{ such that } (\omega_{T,k}^\Delta, \chi) = -(\omega_T^0, \chi), \forall \chi \in \mathcal{H}_{T,k}, \tag{5} \\ \text{find } \psi_{T,k} &\in H_{0,T}^1 \text{ such that } (\nabla \psi_{T,k}, \nabla \xi) = (\omega_T^0 + \omega_{T,k}^\Delta, \xi), \forall \xi \in H_{0,T}^1. \nonumber \end{align} $$

<figure><img src="image_2.png" /><figcaption>Figure 1. – Value of the vorticity on the boundary (mesh A).</figcaption></figure>

<caption>Figure 1. – Valeur du tourbillon sur le bord (maillage A).</caption>
**PROPOSITION 2 (Solution of the discretized harmonic formulation of the Stokes problem).** – The problem (4) has a unique solution $(\psi_{T,k}, \omega_{T,k}) \in H_{0,T}^1 \times H_T^{1,k}$. If $f \in (L^2(\Omega))^2$, if $T$ belongs to a family of regular meshes and if $k \ge K(T)$, with $K(T)$ defined in Theorem 2, the solution depends continuously on the datum $f$.

With classical arguments founded on uniform stability and interpolation theory, taking into account difficulties due to the fact that piecewise linear boundary functions in $\Gamma_T$ ("chinese hat") belong only to the space $H^{3/2-\epsilon}(\Gamma)$, $(\epsilon > 0)$, we have proved [5] the:

**THEOREM 3 (Convergence result).** – If $(\psi, \omega) \in (H^2(\Omega) \cap H_0^1(\Omega)) \times H^2(\Omega)$ is the solution of problem (1) we have:

$$ \begin{gathered} \forall \theta, 0 < \theta < 1, \exists C(\theta) > 0, \forall T \text{ regular, } \exists K(T), \forall k \ge K(T), \\ \|\omega - \omega_{T,k}\|_M + \|\psi - \psi_{T,k}\|_{1,\Omega} \le C h_T^{1-\theta} (\|\omega\|_{2,\Omega} + \|\psi\|_{2,\Omega}). \end{gathered} $$

## 4. Numerical results

Numerical experiments have been performed on a unit square with an analytical test proposed by Bercovier and Engelman [2]. As the integer $K(T)$ a priori depends on the mesh (see Theorem 2), we have worked with different unstructured meshes (A, B, C) [5]. The uncoupled algorithm is the following:

$$ \begin{align} \text{find } \omega_T^0 &\in H_{0,T}^1 \text{ such that } (\nabla \omega_T^0, \nabla \xi) = (f, \operatorname{rot} \xi), \forall \xi \in H_{0,T}^1, \notag \\ \text{find } \omega_{T,k}^\Delta &\in \mathcal{H}_{T,k} \text{ such that } (\omega_{T,k}^\Delta, \chi) = -(\omega_T^0, \chi), \forall \chi \in \mathcal{H}_{T,k}, \tag{5} \\ \text{find } \psi_{T,k} &\in H_{0,T}^1 \text{ such that } (\nabla \psi_{T,k}, \nabla \xi) = (\omega_T^0 + \omega_{T,k}^\Delta, \xi), \forall \xi \in H_{0,T}^1. \notag \end{align} $$

<figure><img src="image_2.png" /><figcaption>Figure 1. – Value of the vorticity on the boundary (mesh A).</figcaption></figure>

<caption>Figure 1. – Valeur du tourbillon sur le bord (maillage A).</caption>
<caption>Table 1. – Extremum of the vorticity after *k* refinements (*N* = number of discrete harmonics).</caption>

<table><caption>Tableau 1. – Extremum du tourbillon après *k* raffinements (*N* = nombre d’harmoniques discrètes).</caption><thead><tr><th>Mesh</th><th>N</th><th>k = 0</th><th>1</th><th>2</th><th>3</th><th>4</th><th>Reference</th></tr></thead><tbody><tr><td>A : 482 elements</td><td>56</td><td>26.52</td><td>20.37</td><td>17.34</td><td>16.37</td><td>16.10</td><td>16.00</td></tr><tr><td>B : 1562 elements</td><td>100</td><td>27.08</td><td>20.78</td><td>17.41</td><td>16.34</td><td>16.05</td><td>16.00</td></tr><tr><td>C : 3188 elements</td><td>144</td><td>31.77</td><td>22.94</td><td>18.20</td><td>16.60</td><td>16.15</td><td>16.00</td></tr></tbody></table>

<figure><img src="image_3.png" /><figcaption>Figure 2. – Convergence order after 4 refinements.</figcaption></figure>

<caption>Figure 2. – Ordre de convergence après 4 raffinements.</caption>

Notice that $\psi$ and $\omega^0$ are discretized on the initial mesh $\mathcal{T}$ and not on the refined one $\mathcal{T}_k$. The step (5) is the most expensive because we have to calculate all the discrete harmonics on the refined meshes. With the original scheme [7], the vorticity presents oscillations on the boundary. Figure 1 shows that these oscillations decrease when the level $k$ of the sub-grid $\mathcal{T}_k$ increases. Table 1 gives extrema of the vorticity obtained after 0 (corresponding to the non-stabilized scheme), 1, 2, 3 and 4 refinements of each mesh (the reference value is +16.00). Notice that the 4 times refined C-mesh is composed of 816128 elements. Finally, we obtain better results than expected: almost $O(h_T^2)$ for the L$^2$-norm of the vorticity and of the stream function (see Figure 2)!

References

[1] Amara M., Bernardi C., Convergence of a finite element discretization of the Navier-Stokes equations in vorticity and stream function formulation, Modél. Math. Anal. Numér. 33 (5) (1999) 1033–1056.

[2] Bercovier M., Engelman M., A finite element for the numerical solution of viscous incompressible flows, J. Comput. Phys. 30 (1979) 181–201.

[3] Bernardi C., Girault V., Maday Y., Mixed spectral element approximation of the Navier-Stokes equations in stream function and vorticity formulation, IMA J. Numer. Anal. 12 (1992) 565–608.

[4] Ciarlet P.G., The Finite Element Methods for Elliptic Problems, North-Holland, 1987.

[5] Dubois F., Salaün M., Salmon S., Discrete harmonics for the stream function-vorticity Stokes problem, Note interne IAT 325/99, 1999.

[6] Girault V., Raviart P.-A., Finite Element Method for Navier-Stokes Equations, Springer-Verlag, 1986.

[7] Glowinski R., Pironneau O., Numerical methods for the first biharmonic equation and for the bidimensional Stokes problem, SIAM 21 (1979) 167–211.
<caption>Table 1. – Extremum of the vorticity after *k* refinements (*N* = number of discrete harmonics).</caption>

<table><caption>Tableau 1. – Extremum du tourbillon après *k* raffinements (*N* = nombre d’harmoniques discrètes).</caption><thead><tr><th>Mesh</th><th>N</th><th>k = 0</th><th>1</th><th>2</th><th>3</th><th>4</th><th>Reference</th></tr></thead><tbody><tr><td>A : 482 elements</td><td>56</td><td>26.52</td><td>20.37</td><td>17.34</td><td>16.37</td><td>16.10</td><td>16.00</td></tr><tr><td>B : 1562 elements</td><td>100</td><td>27.08</td><td>20.78</td><td>17.41</td><td>16.34</td><td>16.05</td><td>16.00</td></tr><tr><td>C : 3188 elements</td><td>144</td><td>31.77</td><td>22.94</td><td>18.20</td><td>16.60</td><td>16.15</td><td>16.00</td></tr></tbody></table>

<figure><img src="image_3.png" /><figcaption>Figure 2. – Convergence order after 4 refinements.</figcaption></figure>

<caption>Figure 2. – Ordre de convergence après 4 raffinements.</caption>

Notice that $\psi$ and $\omega^0$ are discretized on the initial mesh $\mathcal{T}$ and not on the refined one $\mathcal{T}_k$. The step (5) is the most expensive because we have to calculate all the discrete harmonics on the refined meshes. With the original scheme [7], the vorticity presents oscillations on the boundary. Figure 1 shows that these oscillations decrease when the level $k$ of the sub-grid $\mathcal{T}_k$ increases. Table 1 gives extrema of the vorticity obtained after 0 (corresponding to the non-stabilized scheme), 1, 2, 3 and 4 refinements of each mesh (the reference value is +16.00). Notice that the 4 times refined C-mesh is composed of 816128 elements. Finally, we obtain better results than expected: almost $O(h_T^2)$ for the L$^2$-norm of the vorticity and of the stream function (see Figure 2)!

References

[1] Amara M., Bernardi C., Convergence of a finite element discretization of the Navier-Stokes equations in vorticity and stream function formulation, Modél. Math. Anal. Numér. 33 (5) (1999) 1033–1056.

[2] Bercovier M., Engelman M., A finite element for the numerical solution of viscous incompressible flows, J. Comput. Phys. 30 (1979) 181–201.

[3] Bernardi C., Girault V., Maday Y., Mixed spectral element approximation of the Navier-Stokes equations in stream function and vorticity formulation, IMA J. Numer. Anal. 12 (1992) 565–608.

[4] Ciarlet P.G., The Finite Element Methods for Elliptic Problems, North-Holland, 1987.

[5] Dubois F., Salaün M., Salmon S., Discrete harmonics for the stream function-vorticity Stokes problem, Note interne IAT 325/99, 1999.

[6] Girault V., Raviart P.-A., Finite Element Method for Navier-Stokes Equations, Springer-Verlag, 1986.

[7] Glowinski R., Pironneau O., Numerical methods for the first biharmonic equation and for the bidimensional Stokes problem, SIAM 21 (1979) 167–211.