ANNALES MATHÉMATIQUES BLAISE PASCAL

O. BESSON

Finite element solution of Navier-Stokes equations
in shallow domains

Annales mathématiques Blaise Pascal, tome 9, n° 2 (2002), p. 161-180

<http://www.numdam.org/item?id=AMBP_2002__9_2_161_0>

© Annales mathématiques Blaise Pascal, 2002, tous droits réservés.

L'accès aux archives de la revue « Annales mathématiques Blaise Pascal » (http://math.univ-bpclermont.fr/ambp/) implique l'accord avec les conditions générales d'utilisation (http://www.numdam.org/legal.php). Toute utilisation commerciale ou impression systématique est constitutive d'une infraction pénale. Toute copie ou impression de ce fichier doit contenir la présente mention de copyright.

NUMDAM

Article numérisé dans le cadre du programme
Numérisation de documents anciens mathématiques
http://www.numdam.org/
# Finite element solution of Navier-Stokes equations in shallow domains.

O. Besson

## Abstract

A finite element method for the numerical solution of the anisotropic Navier-Stokes equations in shallow domain is presented. This method take into account the low regularity of the vertical component of the velocity in the hydrostatic approximation of the Navier-Stokes equations [2, 3, 5]. A projection method [8] is used for the time discretization. The linear systems are solved via a some preconditioned conjugate algorithm, well adapted to massively parallel computers [4]. Some results are presented for the wind driven water circulation in lakes Geneva and Neuchâtel.

## 1 A common feature of geophysical fluids.

Water flows in oceanography and limnology are governed by the Navier-Stokes equations. In numerical simulations, asymptotic models are in current use (see [10], [12]). These are all based on the following remark:

*The horizontal dimensions are much larger than the vertical one.*

Table (1) illustrates this fact.

<table><caption>Table 1: The puddle law</caption><thead><tr><td></td><td>horizontal width d</td><td>depth h</td><td>ε = h/d</td></tr></thead><tbody><tr><td>Northern Atlantic</td><td>5000 km</td><td>5 km</td><td>0.001</td></tr><tr><td>Lake Neuchâtel (Switzerland)</td><td>38 km</td><td>150 m</td><td>0.004</td></tr><tr><td>Lake Geneva (Switzerland)</td><td>65 km</td><td>300 m</td><td>0.005</td></tr><tr><td>Puddle</td><td>2 m</td><td>1 cm</td><td>0.005</td></tr></tbody></table>
O. BESSON

The simplest model using the fact that

$$ \epsilon = \frac{h}{d} $$

is very small is the hydrostatic model. In this model, we take care of turbulence effects by setting an anisotropic viscosity, much smaller in the vertical direction than in the horizontal one (see [10]).

## 2 Anisotropic Navier-Stokes equations.

Let $W_{\epsilon} \subset \mathbb{R}^3$ be the sufficiently regular domain defined by

$$ W_{\epsilon} = \{ \xi = (\xi_1, \xi_2, \xi_3) \in \mathbb{R}^3, (\xi_1, \xi_2) \in G_s, -h(\xi_1, \xi_2) < \xi_3 < 0 \} $$

where $G_s$ is the surface of the domain and $h: G_s \to \mathbb{R}$ is its depth. Let $G_b = \partial W_{\epsilon} \setminus G_s$ be the bottom of the domain (Fig. 1). It is assumed that the water motion is generated by horizontal tractions, induced by some wind on the surface $G_s$. This motion is driven by the anisotropic Navier-Stokes equations and is influenced by the Coriolis force.

<figure><img src="image_1.png" /><figcaption>Figure 1: Illustration of the shallow domain</figcaption></figure>

$$ v_t + (v \cdot \nabla) \cdot v = \Delta_{\nu} v - 2\omega \wedge v - \nabla p \quad \text{in } W_{\epsilon} $$
NAVIER-STOKES EQUATIONS IN SHALLOW DOMAINS

$$ \operatorname{div} \mathbf{v} = 0 \quad \text{in } W_{\epsilon} $$

$$ \mathbf{v} = 0 \quad \text{on } G_b $$

$$ \nu_3 \frac{\partial v_1}{\partial \xi_3} = \tau_1, \quad \nu_3 \frac{\partial v_2}{\partial \xi_3} = \tau_2, \quad v_3 = 0 \quad \text{on } G_s $$

where $\mathbf{v} = (v_1, v_2, v_3)$ is the fluid velocity, $\omega = (0, 0, \omega_3)$ is the angular velocity of the Earth (projected onto the vertical in local coordinates), $p$ is the pressure, $\nu = (\nu_1, \nu_2, \nu_3)$ is the turbulent viscosity tensor, $\theta_1$ and $\theta_2$ are the tractions induced by the wind and

$$ \Delta_{\nu} \varphi = \sum_{j=1}^{3} \nu_j \frac{\partial^2 \varphi}{\partial \xi_j^2}. $$

Let us do the following change of variables and functions

$$ x_1 = \xi_1, \quad x_2 = \xi_2, \quad x_3 = \xi_3/\epsilon $$

$$ u_1 = v_1, \quad u_2 = v_2, \quad u_3 = v_3/\epsilon $$

and set

$$ \Omega = \left\{ \mathbf{x} = (x_1, x_2, x_3); (x_1, x_2) \in \Gamma_s, -\frac{1}{\epsilon} \cdot h(x_1, x_2) < x_3 < 0 \right\}. $$

$$ \Gamma_s = G_s $$

$$ \Gamma_b = \partial\Omega \setminus \Gamma_s, $$

With this scale change we get

$$ \frac{\partial u_1}{\partial t} + \mathbf{u} \cdot \nabla u_1 - \Delta_{\nu} u_1 - f u_2 + \frac{\partial p}{\partial x_1} = 0 \quad \text{in } \Omega $$

$$ \frac{\partial u_2}{\partial t} + \mathbf{u} \cdot \nabla u_2 - \Delta_{\nu} u_2 + f u_1 + \frac{\partial p}{\partial x_2} = 0 \quad \text{in } \Omega $$

$$ \epsilon^2 \left\{ \frac{\partial u_3}{\partial t} + \mathbf{u} \cdot \nabla u_3 - \Delta_{\nu} u_3 \right\} + \frac{\partial p}{\partial x_3} = 0 \quad \text{in } \Omega $$

$$ \operatorname{div} \mathbf{u} = 0 \quad \text{in } \Omega $$

$$ \mathbf{u} = 0 \quad \text{on } \Gamma_b $$
O. BESSON

$$ \nu_3 \frac{\partial u_1}{\partial x_3} = \epsilon \tau_1, \quad \nu_3 \frac{\partial u_2}{\partial x_3} = \epsilon \tau_2, \quad u_3 = 0 \quad \text{on } \Gamma_s $$

$$ \mathbf{u}(\cdot, t=0) = 0 \quad \text{in } \Omega. $$

with $\nu^\epsilon = (\nu_1, \nu_2, \nu_3/\epsilon^2)$, and $f = 2\omega_3$.

Set

$$ \nu_1 = \lambda_1, \quad \nu_2 = \lambda_2, \quad \nu_3 = \epsilon^2 \lambda_3, $$

$$ \tau_i = \epsilon \theta_i \quad i = 1, 2, $$

$$ \lambda = (\lambda_1, \lambda_2, \lambda_3). $$

When $\epsilon \to 0$, this problem becomes the hydrostatic approximation of Navier-Stokes equations.

$$ \frac{\partial u_1}{\partial t} + \mathbf{u} \cdot \nabla u_1 - \Delta_{\lambda} u_1 - f u_2 + \frac{\partial p}{\partial x_1} = 0 \quad \text{in } \Omega $$

$$ \frac{\partial u_2}{\partial t} + \mathbf{u} \cdot \nabla u_2 - \Delta_{\lambda} u_2 + f u_1 + \frac{\partial p}{\partial x_2} = 0 \quad \text{in } \Omega $$

$$ \frac{\partial p}{\partial x_3} = 0 \quad \text{in } \Omega \qquad (2.1) $$

$$ \operatorname{div} \mathbf{u} = 0 \quad \text{in } \Omega $$

$$ \mathbf{u} = 0 \quad \text{on } \Gamma_b $$

$$ \lambda_3 \frac{\partial u_1}{\partial x_3} = \theta_1, \quad \lambda_3 \frac{\partial u_2}{\partial x_3} = \theta_2, \quad u_3 = 0 \quad \text{on } \Gamma_s $$

$$ u_1(\cdot, t=0) = u_2(\cdot, t=0) = 0 \quad \text{in } \Omega. $$

A week formulation of this problem is the following. Define

$$ V = \{ \varphi \in H^1(\Omega) ; \varphi = 0 \text{ sur } \Gamma_b \} $$

$$ H(\partial_i, \Omega) = \left\{ \varphi \in L^2(\Omega) ; \frac{\partial \varphi}{\partial x_i} \in L^2(\Omega) \right\} $$

$$ H_0(\partial_i, \Omega) = \{ \varphi \in H(\partial_i, \Omega) ; \varphi n_i = 0 \text{ on } \partial \Omega \} $$

$$ L_0^2(\Omega) = \left\{ \varphi \in L^2(\Omega) ; \int_{\Omega} \varphi(x) dx = 0 \right\}. $$
NAVIER-STOKES EQUATIONS IN SHALLOW DOMAINS

If $T > 0$ and $K = L^2(0, T; V \times V \times H_0(\partial_3, \Omega)) \times L^2(0, T; L_0^2(\Omega))$, for $\theta_1, \theta_2 \in H^{-1/2}(\Gamma_s)$, we seek for $(\mathbf{u}, p) \in K$, such that

$$
\begin{align}
& \int_{\Omega} \frac{\partial u_1}{\partial t} v_1 \, dx + \int_{\Omega} (\mathbf{u} \cdot \nabla u_1) v_1 \, dx - f \cdot \int_{\Omega} u_2 v_1 \, dx + \sum_{i=1}^{3} \lambda_i \int_{\Omega} \frac{\partial u_1}{\partial x_i} \frac{\partial v_1}{\partial x_i} \, dx \nonumber \\
& + \int_{\Omega} \frac{\partial u_2}{\partial t} v_2 \, dx + \int_{\Omega} (\mathbf{u} \cdot \nabla u_2) v_2 \, dx + f \cdot \int_{\Omega} u_1 v_2 \, dx + \sum_{i=1}^{3} \lambda_i \int_{\Omega} \frac{\partial u_2}{\partial x_i} \frac{\partial v_2}{\partial x_i} \, dx \nonumber \\
& - \int_{\Omega} p \cdot \operatorname{div} \mathbf{v} \, dx - \int_{\Omega} \operatorname{div} \mathbf{u} \cdot q \, dx \nonumber \\
& = \int_{\Gamma_s} \theta_1 v_1 \, ds + \int_{\Gamma_s} \theta_2 v_2 \, ds \tag{2.2}
\end{align}
$$

for all $\mathbf{v}$ in $V \times V \times H_0(\partial_3, \Omega)$ and all $q$ in $L_0^2(\Omega)$.

The above problem (2.2) has been studied in [2] and [3]. Some finite elements adapted to this problem was established in [1], see also [2]. The stationary case was studied in [6, 5].

3 Some finite element discretization.

Let $\tau_h$ be a finite element mesh of the domain $\Omega$ into hexaedric (or prisatic)
elements and set

$$
V_{1,h} = V_{2,h} = \{ \varphi \in C(\bar{\Omega}) \cap V; \hat{\varphi}|_{\hat{T}} \in Q_2(\hat{x}_1, \hat{x}_2, \hat{x}_3) \ \forall T \in \tau_h \}
$$

$$
V_{3,h} = \{ \varphi \in C(\bar{\Omega}) \cap H_0(\partial_3, \Omega); \hat{\varphi}|_{\hat{T}} \in Q_1(\hat{x}_1, \hat{x}_2, \hat{x}_3) \ \forall T \in \tau_h \}
$$

$$
V_h = V_{1,h} \times V_{2,h} \times V_{3,h}
$$

$$
W_h = \{ \varphi : \Omega \to \mathbf{R}, \varphi |_{\hat{T}} \in Q_1(\hat{x}_1, \hat{x}_2, \hat{x}_3) \ \forall T \in \tau_h \}
$$

where $\hat{T}$ is as usual a reference element. These spaces take into account the degenerated characteristics of the hydrostatic approximation of the Navier-Stokes equations (2.1) (see also [1, 2, 5]).

In this section, some results in the case of the full Navier-Stokes equations
are presented, using the above spaces. For this the following projection
method like in [8] is used, see also [11].
For all $n \ge 0$, assume $\mathbf{u}^n$ and $p^n$ given, find $\mathbf{w}^n \in \mathbf{V}_h$, $\mathbf{w}^n = (w_1^n, w_2^n, w_3^n)$ such that for all $\varphi \in \mathbf{V}_h$

$$ 
\begin{aligned}
\frac{1}{k} \int_{\Omega} w_1^n \varphi_1 + \sum_{i=1}^{3} \lambda_i \int_{\Omega} \frac{\partial w_1^n}{\partial x_i} \frac{\partial \varphi_1}{\partial x_i} &= \int_{\Gamma_s} \theta_1 \varphi_1 \, \mathrm{d}s + \int_{\Omega} p^n \frac{\partial \varphi_1}{\partial x_1} \\
&\quad + f \cdot \int_{\Omega} u_2^n \varphi_1 - \int_{\Omega} (\mathbf{u}^n \cdot \nabla u_1^n) \varphi_1 \\
&\quad + \frac{1}{k} \int_{\Omega} u_1^n \varphi_1
\end{aligned}
\quad (3.3) $$

$$ 
\begin{aligned}
\frac{1}{k} \int_{\Omega} w_2^n \varphi_2 + \sum_{i=1}^{3} \lambda_i \int_{\Omega} \frac{\partial w_2^n}{\partial x_i} \frac{\partial \varphi_2}{\partial x_i} &= \int_{\Gamma_s} \theta_2 \varphi_2 \, \mathrm{d}s + \int_{\Omega} p^n \frac{\partial \varphi_2}{\partial x_2} \\
&\quad - f \cdot \int_{\Omega} u_1^n \varphi_2 - \int_{\Omega} (\mathbf{u}^n \cdot \nabla u_2^n) \varphi_2 \\
&\quad + \frac{1}{k} \int_{\Omega} u_2^n \varphi_2
\end{aligned}
\quad (3.4) $$

$$ 
\begin{aligned}
\frac{1}{k} \int_{\Omega} w_3^n \varphi_3 + \sum_{i=1}^{3} \lambda_i \int_{\Omega} \frac{\partial w_3^n}{\partial x_i} \frac{\partial \varphi_3}{\partial x_i} &= \int_{\Omega} p^n \frac{\partial \varphi_3}{\partial x_3} - \int_{\Omega} (\mathbf{u}^n \cdot \nabla u_3^n) \varphi_3 \\
&\quad + \frac{1}{k} \int_{\Omega} u_3^n \varphi_3
\end{aligned}
\quad (3.5) $$

Then find $q_n \in V_{3,h}$ such that

$$ \sum_{i=1}^{3} \int_{\Omega} \frac{\partial q_n}{\partial x_i} \frac{\partial \psi}{\partial x_i} = \int_{\Omega} \mathrm{div} \, \mathbf{w}^n \, \psi \quad (3.6) $$

for all $\psi \in V_{3,h}$, and set

$$ p^{n+1} = p^n + q^n \quad (3.7) $$

Finally $\mathbf{u}^{n+1}$ is the solution of

$$ \int_{\Omega} \mathbf{u}^{n+1} \cdot \varphi = - \int_{\Omega} q^n \, \mathrm{div} \, \varphi \quad (3.8) $$

for all $\varphi \in \mathbf{V}_h$.

The solution of the linear systems (3.3) to (3.8) are performed via a preconditioned conjugate gradient. Various band matrix preconditioners, well adapted to massively parallel computers can be used. These results are presented in [4].
# 4 Some numerical results: lakes Neuchâtel and Geneva.

Lakes Neuchâtel and Geneva are located in the North-West and West part of Switzerland. Lake Neuchâtel has 38 km long, 7 km width and 150 m depth and lake Geneva has 65 km long, 13 km width and 300 m depth. The level curves are presented in figures (2) and (3). On figure (2), the level curves show *la Motte*, this is a submarine hill culminating at 8 m below the lake level. This hill has an important influence on the water circulation.

<figure><img src="image_2.png" /><figcaption>Figure 2: Lake Neuchâtel, level curves: 30 m</figcaption></figure>
O. BESSON

<figure><img src="image_3.png" /><figcaption>Figure 3: Lake Geneva, level curves: 50 m</figcaption></figure>

The dominant winds in this region are the SW and NW wind. They are exactly along the main axe of lake Neuchâtel.

The currents are presented in lakes Neuchâtel and Geneva. They are induced by a Beaufort 3 SW wind (7.5 m/s) blowing on the lake during 12 hours. The following physical constants are used in the equations (see [7, 9]):

* $\nu = (100, 10, 2.5 \cdot 10^{-3}) [m^2/s]$,

* $f = 10^{-4}$ (latitude 47° N),

* $\theta_i = \alpha \cdot (U_1^2 + U_2^2)^{1/2} U_i$, where $U$ is the wind velocity and $\alpha$ is a traction coefficient ($\alpha = 2.5 \cdot 10^{-6}$).

In the following figures, some numerical results for the currents in the mentioned lakes at different depth and times are presented. First we give some pictures of the currents after 12 hours of wind tractions (i.e. at the time when the wind just stops). Then the residual currents after a 24 hours simulation (i.e. 12 hours after the wind stops) are given. The influence of the Coriolis force on the currents in both lakes can be seen as well as the evident effect of *la Motte* on the water circulation in lake Neuchâtel.
NAVIER-STOKES EQUATIONS IN SHALLOW DOMAINS

<figure><img src="image_4.png" /><figcaption>Figure 4: Lake Neuchâtel, flow after 12h. at the surface, max. max. : 45 cm/s</figcaption></figure>

<figure><img src="image_5.png" /><figcaption>Figure 5: Lake Neuchâtel, flow after 12h. depth 50 m, max. velocity : 6.8 cm/s</figcaption></figure>

<figure><img src="image_6.png" /><figcaption>Figure 6: Lake Neuchâtel, flow after 12h. depth 100 m, max. velocity : 3.2 cm/s</figcaption></figure>
O. BESSON

<figure><img src="image_7.png" /><figcaption>Figure 7: Lake Neuchâtel, flow after 24h. at the surface, max. velocity : 3.5 cm/s</figcaption></figure>

<figure><img src="image_8.png" /><figcaption>Figure 8: Lake Neuchâtel, flow after 24h. depth 50 m, max. velocity : 1.8 cm/s</figcaption></figure>

<figure><img src="image_9.png" /><figcaption>Figure 9: Lake Neuchâtel, flow after 24h. depth 100 m, max. velocity : 0.4 cm/s</figcaption></figure>
NAVIER-STOKES EQUATIONS IN SHALLOW DOMAINS

<figure><img src="image_10.png" /><figcaption>Figure 10: Streamlines in lake Neuchâtel, flow after 12h.</figcaption></figure>
<figure><img src="image_11.png" /><figcaption>Figure 11: Lake Neuchâtel, flow after 12h, vertical velocities in the y-z plane.</figcaption></figure>
NAVIER-STOKES EQUATIONS IN SHALLOW DOMAINS

<figure><img src="image_12.png" /><figcaption>Figure 12: Lake Neuchâtel, flow after 12h, vertical velocities in the y-z plane.</figcaption></figure>
<figure><img src="image_13.png" /><figcaption>Figure 13: Lake Geneva, flow after 12h. at the surface, max. velocity : 50 cm/s</figcaption></figure>

<figure><img src="image_14.png" /><figcaption>Figure 14: Lake Geneva, flow after 12h. depth 50 m, max. velocity : 5.5 cm/s</figcaption></figure>
NAVIER-STOKES EQUATIONS IN SHALLOW DOMAINS

<figure><img src="image_15.png" /><figcaption>Figure 15: Lake Geneva, flow after 12h. depth 100 m, max. velocity : 3.2 cm/s</figcaption></figure>

<figure><img src="image_16.png" /><figcaption>Figure 16: Lake Geneva, flow after 12h. depth 200 m, max. velocity : 2.0 cm/s</figcaption></figure>
<figure><img src="image_17.png" /><figcaption>Figure 17: Lake Geneva, flow after 24h. at the surface, max. velocity : 15 cm/s</figcaption></figure>

<figure><img src="image_18.png" /><figcaption>Figure 18: Lake Geneva, flow after 24h. depth 50 m, max. velocity : 2.7 cm/s</figcaption></figure>
NAVIER-STOKES EQUATIONS IN SHALLOW DOMAINS

<figure><img src="image_19.png" /><figcaption>Figure 19: Lake Geneva, flow after 24h. depth 100 m, max. velocity : 3.4 cm/s</figcaption></figure>

<figure><img src="image_20.png" /><figcaption>Figure 20: Lake Geneva, flow after 24h. depth 200 m, max. velocity : 2.3 cm/s</figcaption></figure>
<figure><img src="image_21.png" /><figcaption>Figure 21: Streamlines in lake Geneva, flow after 12h.</figcaption></figure>
References

[1] P. Azérad. Analyse et approximation du problème de stokes dans un bassin peu profond. *C.R. Acad. Sci. Paris*, 318:53-58, 1994.

[2] P. Azérad. *Equations de Navier-Stokes en bassin peu profond*. PhD thesis, Université de Neuchâtel, 1995.

[3] P. Azérad and F. Guillén-González. Mathematical justification of the hydrostatic approximation in the primitive equations of geophysical fluid dynamics. *SIAM J. Math. Anal.*, 33:847-859, 2001.

[4] O. Besson. Band preconditioners: application to preconditioned conjugate gradient methods on parallel computers. In J. Dongarra and E. J. Kontoghiorghes, editors, *Parallel Numerical Linear Algebra*, pages 111-123. Nova Science Publisher, 2001.

[5] O. Besson and M.R. Laydi. Some estimates for the anisotropic navier-stokes equations and for the hydrostatic approximation. *M2AN - Mod. Math. Ana. Num.*, 26:855-865, 1992.

[6] O. Besson, M.R. Laydi, and R. Touzani. Un modèle asymptotique en océanographie. *C.R. Acad. Sci. Paris*, 310:661-665, 1990.

[7] W. Graf. Hydrodynamics of lake geneva. *Schweiz. Z. Hydrol.*, 45:62-100, 1983.

[8] J. Guermond. Some practical implementations of the projection methods for navier-stokes equations. *M2AN - Mod. Math. Ana. Num.*, 30:637-667, 1996.

[9] K. Hutter. Mathematische vorhersage von barotropen und baroklinen prozessen im zürich und luganersee. *Vierteljahrschrift der Naturforschenden Gesellschaft in Zürich*, 129:51-92, 1984.

[10] J. Pedlosky. *Geophysical fluid dynamics*. Springer-Verlag, 1987.

[11] R. Temam. *Navier-Stokes equations*. North-Holland Elsevier, 1985.

[12] R. Zeytounian. *Modélisation asymptotique en mécanique des fluides newtoniens*. Springer-Verlag, 1994.
O. BESSON

O. BESSON
UNIVERSITÉ DE NEUCHÂTEL
INSTITUT DE MATHÉMATIQUES
11, RUE E. ARGAND
CH-2007 NEUCHÂTEL
SWITZERLAND
olivier.besson@maths.unine.ch