Stokes

Christophe Prud'homme <christophe.prudhomme@cemosis.fr>, Laurent Navoret

Problème de Stokes

Instationnaire
\[\begin{aligned} \partial_{t} u-\Delta u+\nabla p &=f & & \text { sur } \Omega \times \mathbb{R}^{+} \\ \nabla \cdot u &=0 & \text { sur } \Omega \times \mathbb{R}^{+} \\ u &=0 & \text { sur } \partial \Omega \end{aligned}\]

avec \(u: \mathbb{R}^{} \times \Omega \rightarrow \mathbb{R}^{d}] et stem:[p: \mathbb{R}^{} \times \Omega \rightarrow \mathbb{R}\) pression

Stationnaire
\[\begin{aligned} -\Delta u+\nabla p&=f, \quad \mbox{ sur } \Omega\\ \nabla \cdot u &=0 & & \text { sur } \Omega \\ u &=0 & & \text { sur } \partial \Omega \end{aligned}\]

\(u: \Omega \rightarrow \mathbb{R}^{d}\) et \(p: \Omega \rightarrow \mathbb{R}\) pression

Incompressibilité

Proposition (Théorème de transport)

Soit \(\Omega(t)\) un domaine évoluant avec \(u(t, x)\)

\[\begin{aligned} \frac{d}{d t} \int_{\Omega(t)} F(t, x) d x &=\int_{\Omega(t)} \partial_{t} F(t, x) d x+\int_{\partial \Omega(t)} F(t, x)(u(t, x) \cdot n(x)) d s(x) \\ \frac{d}{d t} \int_{\Omega(t)} d x &=0 \quad+\int_{\partial \Omega(t)}(u(t, x) \cdot n(x)) d s(x) \\ &=0 \quad+\int_{\Omega(t)} \nabla \cdot u(t, x) d x \end{aligned}\]

Incompressibilité (Bis)

Proposition (Divergence et incompressibilité)

  • \(\nabla \cdot u:\) taux d’expansion de volume

  • \(\nabla \cdot u=0:\) flot incompressible

Divergence

Exemple: \(u(x, y)=(y, 0)\) \(\nabla \cdot u=0\)

Champs à divergence nulle

Divergence

Definition (en coordonnées polaires(2D)

Pour \(u(r, \theta)=u_{r}(r, \theta) e_{r}+u_{\theta}(r, \theta) e_{\theta},\) nous avons:

\[\nabla \cdot u=\frac{1}{r} \partial_{r}\left(r u_{r}\right)+\frac{1}{r} \partial_{\theta} u_{\theta}, \quad \nabla \times u=\frac{1}{r} \partial_{r}\left(r u_{\theta}\right)-\frac{1}{r} \partial_{\theta} u_{r}\]

  • Exemple: \(u(r, \theta)=r e_{r}\) \(\operatorname{div} u=2\) and rot \(u=0\)

  • Exemple: \(u(r, \theta)=r e_{\theta}\) \(\operatorname{div} u=0\) and rot \(u=2\)

div2 div3

Contrainte d’incompressibilité

Multiplication par une fonction test \(+\) intégration par partie conduit à

\[\int_{\Omega} \nabla u: \nabla v+\int_{\Omega} p \nabla \cdot v=\int_{\Omega} f \cdot v\]

Formulation variationnelle 1:

\[\text { Trouver } u \in V, \quad \int_{\Omega} \nabla u: \nabla v=\int_{\Omega} f \cdot v, \quad \forall v \in V\]

avec \(V=\left\{u \in\left(H_{0}^{1}(\Omega)\right)^{d}, \nabla \cdot u=0\right\}\) Proposition:: \(u\) solution ssi \(u\) est un minimum de la fonctionnelle \(J(v)=\frac{1}{2} \int_{\Omega}|\nabla v|^{2}-\int f \cdot v\) sur \(V\)

On interprete la pression \(p:\) multiplicateur de Lagrange associé à la contrainte d’incompressibilité
Difficulté : construire sous espace discret \(V_{h} \subset V\)

1) Formulation variationnelle

Supposons \(f \in\left(L^{2}(\Omega)\right)^{d}\) Formulation variationnelle 2:

trouver \((u, p) \in\left(H_{0}^{1}(\Omega)\right)^{d} \times L_{*}^{2}(\Omega)\) tels que

\[\begin{array}{ll} \int_{\Omega} \nabla u: \nabla v-\int_{\Omega} p \nabla \cdot v=\int_{\Omega} f \cdot v, & \forall v \in\left(H_{0}^{1}(\Omega)\right)^{d} \\ \int_{\Omega} q \nabla \cdot u=0 & \forall q \in L_{*}^{2}(\Omega) \end{array}\]

avec \(L_{*}^{2}=\left\{q \in L^{2}(\Omega), \int_{\Omega} q=0\right\}\)

\(\rightarrow p\) est définie à une constante près : \(p \in L_{*}^{2}(\Omega)\) permet d’assurer l’unicité

2) Résolution

Problème variationnel

Soit \(X, M\) deux espaces de Hilbert

\[\begin{array}{ll} \text { Trouver }(u, p) \in X \times M, & a(u, v)+b(v, q)=\ell(v), \quad \forall v \in X \\ & b(u, q)=g(q), \quad \forall q \in M \end{array}\]

avec \(a: X \times X \rightarrow \mathbb{R}, b: X \times M \rightarrow \mathbb{R}\) bilinéaires, \(\ell: V \rightarrow \mathbb{R}\) et \(g: M \rightarrow \mathbb{R}\) linéaires.

Condition inf-sup

Proposition

\(V\) espace de Hilbert. a: forme bilinéaire, continue, coercive sur \(X\) b: forme bilinéaire, continue sur \(X \times M\) \(\ell, g:\) forme linéaire continue. Le problème variationnel admet une unique solution \((u, p) \in X \times M\) ssi \(\exists \beta>0\)

\[\inf _{q \in M_{v \in X}} \sup _{\|v\|_{X}\|q\|_{M}} \frac{b(v, q)}{\|v\|_{X}\|q\|_{M}} \geqslant \beta\]

Si a est symmétrique, la solution est un point selle du Lagrangien:

\[L(v, q)=\frac{1}{2} a(v, v)-b(v, q)-f(v)-g(q)\]

Méthode de Galerkin

Definition ( Problème variationnel discret)

Soit \(X_{h} \subset X, M_{h} \subset M\)

\[\begin{array}{ll} \text { Trouver }\left(u_{h}, p_{h}\right) \in X_{h} \times M_{h}, & a\left(u_{h}, v_{h}\right)+b\left(v_{h}, p_{h}\right)=\ell\left(v_{h}\right), \quad \forall v_{h} \in X_{h} \\ & b\left(u_{h}, q_{h}\right)=g\left(q_{h}\right), \quad \forall q_{h} \in M_{h} \end{array}\]
Proposition (Condition inf-sup discrète)

Le problème variationnel approché admet une unique solution \(\left(u_{h}, p_{h}\right) \in X_{h} \times M_{h}\) ssi \(\exists \beta>0\)

\[\inf _{q_{h} \in M_{h} v_{h} \in X_{h}} \sup _{\left\|v_{h}\right\|_{X}\left\|q_{h}\right\|_{M}} \frac{b\left(v_{h}, q_{h}\right)}{\left\|v_{h}\right\|_{X}\left\|q_{h}\right\|_{M}} \geqslant \beta\]

Méthode de Galerkin

\(X_{h} \subset X\) s.e.v de dimension \(N_{u},\left(\varphi_{i}\right)\) base de \(X_{h}\) \(M_{h} \subset M\) s.e.v de dimension \(N_{p},\left(\psi_{i}\right)\) base de \(M_{h}\)

\[\begin{aligned} u_{h}(x)=\sum_{j=1}^{N_{u}} u_{j} \varphi_{j}(x), & \quad p_{h}(x)=\sum_{j=1}^{N_{p}} p_{j} \psi_{j}(x) & & \\ & & & \\ & a\left(u_{h}, v_{h}\right)+b\left(v_{h}, p_{h}\right)=\ell\left(v_{h}\right), & & \forall v_{h} \in X_{h} \\ & b\left(u_{h}, q_{h}\right)=g\left(q_{h}\right), & & \forall q_{h} \in M_{h} \\ \Longleftrightarrow & \sum_{j=1}^{N_{u}} u_{j} a\left(\varphi_{j}, v_{h}\right)+\sum_{j=1}^{N_{p}} p_{j} b\left(v_{h}, \psi_{j}\right)=\ell\left(v_{h}\right) & \forall v_{h} \in X_{h} \\ & \sum_{j=1}^{N_{u}} u_{j} b\left(\varphi_{j}, q_{h}\right)=g\left(q_{h}\right), & \forall q_{h} \in M_{h} \end{aligned}\]

Formulation matricielle

\[\Longleftrightarrow \quad \sum_{j=1}^{N_{u}} u_{j} a\left(\varphi_{j}, \varphi_{i}\right)+\sum_{j=1}^{N_{p}} p_{j} b\left(\varphi_{i}, \psi_{j}\right)=\ell\left(\varphi_{i}\right) \quad \forall i \in \mathbb{I}, N_{u} \mathbb{J}\]
\[\Longleftrightarrow \quad\left[\begin{array}{ll} A & B^{T} \\ B & 0 \end{array}\right]\left[\begin{array}{l} U \\ P \end{array}\right]=\left[\begin{array}{l} F \\ G \end{array}\right]\]

avec

\[\begin{aligned} A_{i, j} &=\left(a\left(\varphi_{i}, \varphi_{j}\right)\right)_{i, j} \in M_{N_{u}}(\mathbb{R}) & U &=\left(u_{j}\right), \quad F=\left(f\left(\varphi_{j}\right)\right)_{j} \in \mathbb{R}^{N_{u}} \\ B &=\left(b\left(\varphi_{j}, \psi_{i}\right)\right)_{i, j} \in M_{N_{u}, N_{p}}(\mathbb{R}) & P & =\left(p_{j}\right), \quad G=\left(g\left(\psi_{i}\right)\right)_{j} \in \mathbb{R}^{N_{p}} \end{aligned}\]

Formulation matricielle

\[\left[\begin{array}{cc} A & B^{T} \\ B & 0 \end{array}\right]\left[\begin{array}{l} U \\ P \end{array}\right]=\left[\begin{array}{l} F \\ G \end{array}\right]\]

  • \(\rightarrow A\) est symétrique définie positive

  • \(\rightarrow\) La matrice est symétrique mais n’est pas définie positive

  • \(\rightarrow\) Condition inf-sup discrète ssi \(B\) est de rang \(N_{p}\) ssi \(B^{T}\) est de rang \(N_{p}\)

\[\text { ssi } \operatorname{ker}\left(B^{T}\right)=\{0\}\]

Dans ce cas, la matrice est inversible.

Formulation matricielle

La condition inf-sup discrète n’est pas satisfaite dans les cas suivants :

  • \(N_{p}>N_{u}:\) l’espace vectoriel pour la pression est trop grand comparée à celui de la vitesse pour que \(B\) soit de rang \(N_{p}\)

    • \(\rightarrow\) Locking

    • \(\rightarrow B\) est alors en général de rang. \(N_{u}\) et \(\operatorname{ker}(B)=\{0\}:\) la seule solution possible est \(U=0\)

  • il existe \(Q^{*} \in \operatorname{ker}\left(B^{T}\right)\) non nul. Alors \(q_{h}^{*}=\sum_{j=1}^{N_{p}} Q_{k}^{*} \psi_{k}\) vérifie \(b\left(v_{h}, q_{h}^{*}\right)=0\) pour tout \(v_{h} \in V_{h}\)

    • \(\rightarrow\) mode parasite (spurious mode)

Espaces d’approximation

  • 3) Choix du maillage (généralement affine)

  • 4) Choix des espaces d’approximation \(X_{h}, M_{h}\) \(X_{h}=\left\{v_{h} \in(C(\bar{\Omega}))^{d}, \quad \forall K \in \mathcal{T}_{h}, v_{h} \circ T_{K} \in\left(\hat{P}_{u}\right)^{d}\right\} \subset\left(H^{1}(\Omega)\right)^{d}\) \(M_{h}=\tilde{M}_{h} \cap L_{*}^{2}(\Omega)\) avec \(\tilde{M}_{h}=\left\{v_{h} \in C(\bar{\Omega}), \quad \forall K \in \mathcal{T}_{h}, v_{h} \circ T_{K} \in \hat{P}_{p}\right\} \subset L^{2}(\Omega)\)

    Choix de \(\hat{P}_{u}\) et \(\hat{P}_{p}:\)

  • \(\mathbb{P}_{1} / \mathbb{P}_{0} \rightarrow\) locking

  • \(\mathbb{P}_{1} / \mathbb{P}_{1} \rightarrow\) mode parasite

  • \(\cdot \mathbb{P}_{1}\) -bulle/P \(_{1}\) (mini-element) \(\rightarrow\) ok

  • \(\mathbb{P}_{k} / \mathbb{P}_{k-1}(\text { Taylor-Hood }) \rightarrow\) ok

\(\mathbb{P}_{1}\) -bulle \(/ \mathbb{P}_1\)

\[\hat{P}_{u}=\left(\mathbb{P}_{1}+\operatorname{vect}(\hat{b})\right)^{d} \quad \hat{P}_{p}=\mathbb{P}_{1}\]

Fonction bulle: \(\hat{b} \in H_{0}^{1}(\hat{K})\) telle que \(0 \leqslant \hat{b} \leqslant 1\) et \(\hat{b}(\text { barycentre})=1\)

Exemple

\(\hat{b}=3^{3} \lambda_{1} \lambda_{2} \lambda_{3}\)

Proposition

\(\Omega\) polyédrique et \(\left(\mathcal{T}_{h}\right)\) une famille régulière de triangulations. Alors la condition inf-sup discrète est satisfaite uniformément pour tout \(h\) Soit \(f \in L^{2}(\Omega)\) et \(u \in V\) la solution du problème. Soit \(\left(u_{h}, p_{h}\right) \in X_{h} \times M_{h}\) la solution approchée du problème de Stokes. Alors si \(u \in\left(H^{2}(\Omega) \cap H_{0}^{1}(\Omega)\right)^{d}\) et \(p \in H^{1}(\Omega) \cap L_{*}^{2}(\Omega),\) alors

\[\left\|u-u_{h}\right\|_{\left(H^{1}\right)^{d}}+\left\|p-p_{h}\right\|_{L^{2}} \leqslant C h\left(\|u\|_{H^{2}}+\|p\|_{H^{1}}\right)\]

\(\mathbb{P}_{1}\) -bulle \(/ \mathbb{P}_1\)

Si de plus le problème est régularisant (toute solution du problème (adjoint) est régulière ), alors

\[\left\|u-u_{h}\right\|_{\left(L^{2}\right)^{d}} \leqslant C h^{2}\left(\|u\|_{H^{2}}+\|p\|_{H^{1}}\right)\]

  • \(\rightarrow\) pas de convergence en \(h^{2}\) pour la pression \(\hat{P}_{p}=\mathbb{P}_{1}\)

  • \(\rightarrow\) pas de convergence en \(h^{3}\) pour la vitesse alors que \(\hat{P}_{u}\) contient strictement \(\mathbb{P}_{1}\)

\(\hat{P}_{u}=\mathbb{P}_{2} \quad \hat{P}_{p}=\mathbb{P}_{1} \)

Proposition

\(\Omega\) polyédrique et \(\left(\mathcal{T}_{h}\right)\) une famille régulière de triangulations.

Alors la condition inf-sup discrète est satisfaite uniformément pour tout \(h\) en prenant \(X_{h}\)

Soit \(f \in L^{2}(\Omega)\) et \(u \in V\) la solution du problème. Soit \(\left(u_{h}, p_{h}\right) \in X_{h} \times M_{h}\) la solution approchée du problème de Stokes.

Alors si \(u \in\left(H^{3}(\Omega) \cap H_{0}^{1}(\Omega)\right)^{d}\) et \(p \in H^{2}(\Omega) \cap L_{*}^{2}(\Omega),\) alors

\[\left\|u-u_{h}\right\|_{\left(H^{1}\right)^{d}}+\left\|p-p_{h}\right\|_{L^{2}} \leqslant C h^{2}\left(\|u\|_{\left(H^{3}\right)^{d}}+\|p\|_{H^{2}}\right)\]

\(\hat{P}_{u}=\mathbb{P}_{2} \quad \hat{P}_{p}=\mathbb{P}_{1} \)

Si de plus le problème est régularisant (toute solution du problème (adjoint) est régulière ), alors

\[\left\|u-u_{h}\right\|_{\left(L^{2}\right)^{d}} \leqslant C h^{3}\left(\|u\|_{\left(H^{3}\right)^{d}}+\|p\|_{H^{2}}\right)\]

\(\rightarrow\) convergence en \(h^{2}\) pour la pression \(\hat{P}_{p}=\mathbb{P}_{1}\) \(\rightarrow\) convergence en \(h^{3}\) pour la vitesse avec que \(\hat{P}_{u}=\mathbb{P}_{2}\)

Résolution

Matrice symétrique mais non définie positive : les méthodes itératives classiques fonctionnent mal
Deux méthodes

  • Equation sur la pression: \(\rightarrow\left(B A^{-1} B^{T}\right)\) sym. def. pos. : complément de Schur \(\rightarrow\) Calcul de \(P\) puis \(U\) \(\rightarrow\) méthode itérative d’Uzawa

  • Méthode de compressibilité artificielle : on considère la matrice

\[\begin{bmatrix} A & B^{T} \\ B & \varepsilon M \end{bmatrix}\]

avec \(\varepsilon>0,\) qui est définie positive et \(M=\left(\int \psi_{i} \psi_{j}\right)\)

En pratique

\[M_{h}=\tilde{M}_{h} \cap L_{*}^{2}(\Omega)\]

avec

\[\tilde{M}_{h}=\left\{v_{h} \in C(\bar{\Omega}), \quad \forall K \in \mathcal{T}_{h}, v_{h} \circ T_{K} \in \hat{P}_{p}\right\} \subset L^{2}(\Omega)\]

  • \(\rightarrow\) numériquement, on ne considère pas \(M_{h}\) mais \(\tilde{M}_{h}\)

  • \(\rightarrow\) le système linéaire a alors plusieurs solutions puisque la pression n’est définie qu’à une constante près

    Remèdes

    1. si les méthodes itératives convergent tout de même vers une des solutions,

    2. il faut ensuite retrancher à la pression sa moyenne. la méthode de compressibilité artificielle sélectionne directement la pression à moyenne nulle.