RIA
laboria

laboratoire de recherche
en informatique
et automatique

Analyse numérique
du refroidissement
d'un cylindre radioactif

J.F. Bourgat,
R. Trémolières

Institut de Recherche
d'Informatique
et d'Automatique

Domaine de Voluceau
Rocquencourt
78150 - Le Chesnay
Tél. 954.90.20

Rapport de Recherche n°17

avril 1973
ANALYSE NUMÉRIQUE DU REFROIDISSEMENT
D'UN CYLINDRE RADIOACTIF

J.F. Bourgat - R. Tremolières

(présenté par R. Glowinski)

On étudie dans cet article la détérioration de déchets radio-
actifs conservés sous forme de cylindres. Ces cylindres sont constitués
par un mélange de verre et de matière radioactive. Afin d'empêcher
qu'ils ne fondent sous l'effet de l'échauffement naturel dû à la radio-
activité, on place les cylindres dans un courant d'air. Il s'agit de
déterminer les caractéristiques du courant d'air (température initiale,
vitesse) afin de limiter l'échauffement. Le problème est de la catégorie
des problèmes variationnels non-symétriques.
# 1 - HYPOTHESES, NOTATIONS, FORMULATION DU PROBLEME.

On considère un cylindre de révolution $\Omega$, d'axe $Oz$, de rayon $R$ et de longueur $a$. On désigne par $\Gamma_o$ et $\Gamma_a$ les sections droites pour $z=0$ et $z=a$, et par $\Gamma_R$ la frontière latérale du cylindre.

On note $(r, \theta, z)$ les coordonnées cylindriques et $(x, y, z)$ les coordonnées cartésiennes.

Soient $u(x,y,z)$ ou $u(r,\theta,z)$ la température dans le cylindre, et $w(z)$ la température à la cote $z$ d'un courant d'air longeant latéralement le cylindre et destiné à le refroidir. Par symétrie on suppose que $w$ est uniquement fonction de $z$ et pour simplifier on suppose que la température de l'air à l'extérieur du cylindre est, pour $z$ fixé, indépendante de $r > R$.

<figure><img src="image_1.png" /><figcaption>Fig. 1.1</figcaption></figure>
On appelle $w_o$ la température initiale du courant d'air à la cote z=0. L'équilibre du système cylindre-courant d'air se traduit alors par les conditions :

$$ (1.1) \qquad - \Delta u = f \quad \text{dans } \Omega $$

$$ (1.2) \qquad \frac{\partial u}{\partial n} = 0 \quad \text{sur } \Gamma_o \text{ et } \Gamma_a $$

$$ (1.3) \qquad \frac{\partial u}{\partial n} = \frac{h}{\lambda} (w-u) \quad \text{sur } \Gamma_R $$

$$ (1.4) \qquad \frac{\partial w}{\partial z} = k(u-w) \quad \text{sur } \Gamma_R $$

$$ (1.5) \qquad w(0) = w_o $$

avec

$$ f = f(r,z) = \frac{Q}{\lambda} $$

où

Q = quantité de chaleur produite par unité de volume et par seconde dans le cylindre,

λ = conductivité calorifique du cylindre

et

h = coefficient d'arrachement à la surface

et

$$ k = \frac{2 \pi R h}{MC} $$

où

M = débit massique total de l'air par seconde,

C = chaleur spécifique de l'air.
On peut remarquer qu'en intégrant (1.4) (1.5), on obtient

$$
(1.6) \qquad w(u)(z) = k \left( \int_0^z u e^{ky} dv + w_0 \right) e^{-kz}
$$

et le système (1.1)...(1.5) se réduit à un ensemble d'équations en u.

Pour simplifier, nous prendrons, ce qui ne restreint pas la
généralité,

$w_0 = 0$

$$
\frac{h}{\lambda} = 1.
$$

Dans ce qui suit, nous nous proposons d'étudier l'existence,
l'unicité et la résolution du problème ainsi défini :

$$
(1.7) \qquad - \Delta u = f \quad \text{dans } \Omega,
$$

$$
(1.8) \qquad \frac{\partial u}{\partial n} = 0 \quad \text{sur } \Gamma_0 \quad \text{et} \quad \Gamma_a,
$$

$$
(1.9) \qquad \frac{\partial u}{\partial n} = w - u \quad \text{sur } \Gamma_R,
$$

$$
(1.10) \qquad \frac{\partial w}{\partial z} = u - w \quad \text{sur } \Gamma_R,
$$

$$
(1.11) \qquad w(0) = 0.
$$

Remarque 1.1

On pourrait généraliser ce problème sans difficulté en le formulant comme un problème de contrôle, la variable de décision étant $w(0)$. Le critère serait alors de minimiser le coût dû à $w(0)$ tout en respectant une contrainte du genre $u(x) \le T, \forall x \in \Omega$, ou $T$ est une température limite au-delà de laquelle le verre fondrait sous l'échauffement.
## 2 - FORMULATION VARIATIONNELLE DU PROBLEME

On introduit un sous-espace $V$ de $H^1(\Omega)$ défini par

$$V = \{ v \mid v \in H^1(\Omega), v \text{ indépendant de } \theta \}^{(1)}$$

C'est un espace de Hilbert pour le produit scalaire induit par celui de $H^1(\Omega)$.

Soit $v \in V$, nous déduisons de (1.7) :

$$- (\Lambda u, v)_{L^2(\Omega)} = (f, v)_{L^2(\Omega)}$$

puis en utilisant formellement la formule de Green il vient :

$$\int_{\Omega} \text{grad } u \cdot \text{grad } v \, d\Omega - \int_{\Gamma} \frac{\partial u}{\partial n} \cdot v \, d\Gamma = (f, v)$$

et compte tenu de (1.8) et (1.9) il vient :

$$\int_{\Omega} \text{grad } u \cdot \text{grad } v \, d\Omega - \int_{\Gamma_R} (w(u)-u)vd\Gamma = (f,v)$$

où $w(u)$ est donné par (1.6).

On considère alors sur $V$ la forme bilinéaire définie par :

$$(2.0) \qquad a(u,v) = \int_{\Omega} \text{grad } u \cdot \text{grad } v \, d\Omega - \int_{\Gamma_R} (w(u)-u)v \, d\Gamma.$$

<footnote>(1) Pour les espaces $L^2(\Omega)$, $H^1(\Omega)$, $H^{\frac{1}{2}}(\Gamma)$, etc...et divers théorèmes que nous utiliserons ici cf. LIONS-MAGENES [1] et DUVAUT-LIONS [1]. Tous les espaces considérés ici seront définis sur le corps des réels.</footnote>
Sous forme variationnelle le problème (1.7)...(1.11) s'écrit alors :

$$
(2.1) \quad \left[
\begin{array}{l}
\underline{\text{trouver}} \ u \in \mathcal{V} \text{ tel que} \\
\quad a(u,v) = (f,v), \ \forall v \in \mathcal{V}.
\end{array}
\right.
$$

Nous allons étudier l'existence et l'unicité de u solution
de (2.1) en nous appuyant sur le lemme de Lax-Milgram suivant lequel
le problème (2.1) a une solution unique si $a(u,v)$ est un opérateur
linéaire **coercif** (i.e. $\rho a(v,v) \ge b(v,v)$, $\forall v \in \mathcal{V}$, où $b(u,v)$ est un
produit scalaire faisant de $\mathcal{V}$ un espace de Hilbert) et **continu**
(i.e. $|a(u,v)| \le C b(u)b(v)$, où l'on note $b(u) = \sqrt{b(u,u)}$, $b(v) = \sqrt{b(v,v)}$).

2.1 - *Définition d'une norme équivalente sur V et coercité de a(u,v).*

Nous allons introduire un nouveau produit scalaire $\langle\cdot\rangle_b^{(1)}(u,v)$ sur
$V$ tel que

$$
(2.2) \qquad b(v,v) = a(v,v),
$$

$$
(2.3) \qquad \sqrt{b(v,v)} \text{ soit une norme équivalente à la norme usuelle sur } H^1(\Omega).
$$

On a

$$
a(v,v) = \int_{\Omega} |\operatorname{grad} v|^2 d\Omega + \int_{\Gamma_R} (v^2 - vw(v)) d\Gamma,
$$

or

$$
v^2 - vw(v) = (v - w(v))^2 + w(v) (v - w(v))
$$

<footnote>(1) une forme bilinéaire symétrique strictement définie positive.</footnote>
ce qui compte tenu de (1.10) donne

$$v^2 - vw(v) = \left(\frac{\partial w}{\partial z}\right)^2 + \frac{1}{2} \frac{\partial^2 w}{\partial z^2}.$$

On a alors

$$a(v,v) = \int_{\Omega} |\text{grad } v|^2 d\Omega + \int_{\Gamma_R} \left(\frac{\partial w(v)}{\partial z}\right)^2 d\Gamma + \frac{1}{2} \int_{\Gamma_R} \frac{\partial^2 w(v)}{\partial z^2} d\Gamma,$$

et en intégrant le dernier terme

$$(2.4) \qquad a(v,v) = \int_{\Omega} |\text{grad } v|^2 d\Omega + \int_{\Gamma_R} \left(\frac{\partial w(v)}{\partial z}\right)^2 d\Gamma + \pi_R w^2(v) (a).$$

On introduit alors

$$(2.5) \qquad b(u,v) = \int_{\Omega} \text{grad } u \text{ grad } v \, d\Omega + \int_{\Gamma_R} \frac{\partial w(u)}{\partial z} \frac{\partial w(v)}{\partial z} \, d\Gamma_R + \pi_R w(u) (a) w(v) (a).$$

Il est immédiat que $h(u,v)$ définit un produit scalaire et que (2.2) est assuré.

D'autre part $\sqrt{h(v,v)}$ est une norme car

$$\begin{align*}
\frac{\partial w(v)}{\partial z} &= 0 \text{ sur } \Gamma_R \text{ et } w(v)(a) = 0 \Rightarrow w(v) = 0 \Rightarrow \\
&\Rightarrow v|_{\Gamma_R} = 0 \text{ d'après (1.10).}
\end{align*}$$
Il reste à montrer (2.3).

On sait (cf. NECAS [1]) que la norme usuelle sur $H^1(\Omega)$ :

$$ (2.6) \qquad \| v \| = \left( \int_{\Omega} |\operatorname{grad} v|^2 d\Omega + \int_{\Omega} |v|^2 d\Omega \right)^{\frac{1}{2}} $$

est équivalente à la norme

$$ (2.7) \qquad \| |v| \| = \left( \int_{\Omega} |\operatorname{grad} v|^2 d\Omega + \int_{\Gamma_R} |v|^2 d\Gamma \right)^{\frac{1}{2}}, $$

et compte tenu de (2.4) il suffira de montrer qu'il existe deux constantes strictement positives $C_1, C_2$, telles que :

$$ (2.8) \qquad \begin{aligned} & C_1 \left[ \int_{\Gamma_R} \left( \frac{\partial w(v)}{\partial z} \right)^2 d\Gamma + \pi R w^2(v)(a) \right] \le \int_{\Gamma_R} v^2 d\Gamma \le \\ & \le C_2 \left[ \int_{\Gamma_R} \left( \frac{\partial w(v)}{\partial z} \right)^2 d\Gamma + \pi R w^2(v)(a) \right]. \end{aligned} $$

Or d'après (1.10)

$$ \frac{\partial w(v)}{\partial z} = v - w(v), $$

ce qui entraîne

$$ \int_{\Gamma_R} v^2 d\Gamma = \int_{\Gamma_R} w^2(v) d\Gamma + \int_{\Gamma_R} \left( \frac{\partial w(v)}{\partial z} \right)^2 d\Gamma + \int_{\Gamma_R} \frac{\partial w^2(v)}{\partial z} d\Gamma . $$

Mais sur l'espace

$$ v(\Gamma_R) = \{ v \in H^1(\Gamma_R), v \text{ indépendant de } \theta \} $$
les normes

$$ \left( \left( \int_{\Gamma_R} \left( \frac{\partial \varphi}{\partial z} \right)^2 d\Gamma + \int_{\Gamma_R} \varphi^2 d\Gamma \right)^{\frac{1}{2}} \right) $$

et

$$ \left( \left( \int_{\Gamma_R} \left( \frac{\partial \varphi}{\partial z} \right)^2 d\Gamma + \varphi^2 (a) \right)^{\frac{1}{2}} \right) $$

sont équivalentes, d'où l'on déduit :

$$ \begin{aligned} & c_1' \left[ \int_{\Gamma_R} \left( \frac{\partial w(v)}{\partial z} \right)^2 d\Gamma + \int_{\Gamma_R} \frac{\partial w^2(v)}{\partial z} d\Gamma + w^2(v) \quad (a) \right] \le \int_{\Gamma_R} v^2 d\Gamma \le \\ & \le c_2' \left[ \int_{\Gamma_R} \left( \frac{\partial w(v)}{\partial z} \right)^2 d\Gamma + \int_{\Gamma_R} \frac{\partial w^2(v)}{\partial z} d\Gamma + w^2(v) \quad (a) \right], \end{aligned} $$

or

$$ \int_{\Gamma_R} \frac{\partial w^2(v)}{\partial z} d\Gamma = 2\pi R_w^2(v) (a) $$

d'où (2.8).

Ainsi l'espace V muni de la norme $\sqrt{b(v,v)}$ est un espace de Hilbert et $a(u,v)$ est une forme coercive sur $V.(^1)$

<footnote>(^1) On donne une autre démonstration en Annexe.</footnote>
## 2.2 - Continuité de $a(u,v)$

Posant $b(v)=\sqrt{b(v,v)}$, il faut montrer que :

$$ (2.10) \qquad |a(u,v)| \le C b(u)b(v), \forall u,v \in V. \text{ avec } C > 0. $$

En utilisant pour (2.0) la majoration

$$ (\sum a_i b_i)^2 \le (\sum a_i^2)(\sum b_i^2) $$

on obtient

$$ a(u,v) \le \left( \int_{\Omega} |\text{grad } u|^2 d\Omega + \int_{\Gamma_R} (w(u)-u)^2 d\Gamma \right) \left( \int_{\Omega} |\text{grad } v|^2 d\Omega + \int_{\Gamma_R} v^2 d\Gamma \right), $$

ce qui compte tenu de (1.10) et (2.8) entraîne (2.10).

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

Les résultats des paragraphes 2.1 et 2.2 assurent alors l'existence et l'unicité d'une solution de (2.1).

<figure><img src="image_3.png" /></figure>
### 3 - REDUCTION DU PROBLEME DANS $\mathbb{R}^3$ A UN PROBLEME DANS $\mathbb{R}^2$.

Nous allons tenir compte de la symétrie pour simplifier le problème (2.1) et le ramener à un problème dans $\mathbb{R}^2$.

En utilisant les relations

$$r = \sqrt{x^2 + y^2},$$

$$ \frac{\partial r}{\partial x} = \frac{x}{r}, \quad \frac{\partial r}{\partial y} = \frac{y}{r}, \quad \frac{\partial u}{\partial x_i} \frac{\partial v}{\partial x_i} = \frac{\partial u}{\partial r} \frac{\partial v}{\partial r} \frac{x_i^2}{r^2}, \quad i=1,2,3, $$

la forme bilinéaire

$$ a(u,v) = \sum_{i=1}^{3} \int_{\Omega} \left( \frac{\partial u}{\partial x_i} \frac{\partial v}{\partial x_i} - \frac{\partial v}{\partial x_i} \frac{\partial u}{\partial x_i} \right) d\Omega - \int_{\Gamma_R} (w(u)-u) v d\Gamma $$

s'écrit

$$ a(u,v) = \int_{\Omega} \left( \frac{\partial u}{\partial r} \frac{\partial v}{\partial r} + \frac{\partial u}{\partial z} \frac{\partial v}{\partial z} \right) rdrd\theta dz - \int_{\Gamma_R} (w(u)-u) v^R d\theta dz. $$

Par symétrie, en intégrant en $\theta$, il vient :

$$ a(u,v) = 2\pi \int_{0}^{a} \int_{0}^{R} \left( \frac{\partial u}{\partial r} \frac{\partial v}{\partial r} + \frac{\partial u}{\partial z} \frac{\partial v}{\partial z} \right) rdrdz - 2\pi R \int_{0}^{a} (w(u)-u) v dz. $$

Dans l'espace $(r,z)$ définissons

$$ \tilde{\Omega} = ]0,R[x, 0,a[. $$

<figure><img src="image_4.png" /><figcaption>Fig. 3.1</figcaption></figure>

Le problème (2.1) est alors équivalent au problème variationnel dans $\mathbb{R}^2$ :

$$ \begin{array}{l} \underline{\text{trouver}} \quad \tilde{u} \in H^1(\tilde{\Omega}) \quad \text{tel que} \\ \tilde{a}(\tilde{u}, \tilde{v}) = (\text{fr}, \tilde{v})_{L^2(\tilde{\Omega})}, \quad \forall \tilde{v} \in H^1(\tilde{\Omega}), \end{array} \qquad (3.1) $$
avec

$$
\tilde{a}(\tilde{u}, \tilde{v}) = \int_{0}^{a} \int_{0}^{P} \left( \frac{\partial \tilde{u}}{\partial r} \frac{\partial \tilde{v}}{\partial r} + \frac{\partial \tilde{u}}{\partial z} \frac{\partial \tilde{v}}{\partial z} \right) rdrdz - \int_{\Gamma_R} (w(\tilde{u}) - \tilde{u}) v dz.
$$

On a évidemment

$$
\begin{align*}
& \text{u solution de (2.1) } \Rightarrow \text{u} \Big|_{\partial \Omega} \text{ solution de (3.1),} \\
& \tilde{u} \text{ solution de (3.1) } \Rightarrow \text{u}(r, \theta, z) = \tilde{u}(r, z), \forall z, \text{ solution de} \\
& \qquad (2.1).
\end{align*}
$$

On peut démontrer sans difficulté l'existence et l'unicité de
\tilde{u} solution de (3.1). Pour cela on introduit

$$
\tilde{b}(\tilde{u}, \tilde{v}) = \int_0^a \int_0^R \left( \frac{\partial \tilde{u}}{\partial r} \frac{\partial \tilde{v}}{\partial r} + \frac{\partial \tilde{u}}{\partial z} \frac{\partial \tilde{v}}{\partial z} \right) rdrdz + \\
+ \int_0^a \frac{\partial \tilde{w}(\tilde{v})}{\partial z} \frac{\partial \tilde{w}(\tilde{v})}{\partial z} + \tilde{w}(\tilde{u})(a) \tilde{w}(\tilde{v})(a),
$$

et, raisonnant de la même façon qu'en 2, on montre que

- $\sqrt{h(\tilde{v}, \tilde{v})}$ est une norme sur $H^1(\Omega)$,

- $H^1(\bar{\Omega})$ est complet pour cette norme,

- $\tilde{a}(\tilde{u}, \tilde{v})$ est coercive pour cette norme, en effet
  $\tilde{a}(\tilde{v}, \tilde{v}) = h(\tilde{v}, \tilde{v})$,

- $\tilde{a}(\tilde{u}, \tilde{v})$ est continu pour cette norme.
# 4 - DISCRETISATION

## 4.1 - Notations, formulation du problème approché

D'après les résultats du paragraphe précédent, le problème (2.1) est entièrement résolu dès que l'on connait la solution $\tilde{u}$ du problème (3.1). C'est ce dernier problème que nous allons résoudre de façon approchée en utilisant la technique de discrétisation par différence finie. (cf. AUBIN [1], BOURGAT [1], CEA [1], GLOWINSKI-LIONS-TREMOLIFRES [1], LIONS [1], RAVIART [1], TREMOLIERES [2]).

On pose $h=(h_r, h_z)=(h_1, h_2)$, $h_r>0$, $h_z>0$, $(N_z=\frac{R}{h_r}, N_z=\frac{a}{h_z}$ étant entiers) et l'on introduit :

- un réseau $\mathbb{R}_h$ :

$$ \mathbb{R}_h = \{ \mathbf{M} \mid \mathbf{M} \in \mathbb{R}^2, \mathbf{M} = (m_1 h_1, m_2 h_2), m_i \in \mathbb{Z}, i=1,2 \}, \\ \text{où } \mathbb{Z} \text{ est l'espace des entiers positifs et négatifs,} $$

- des pavés centrés en chaque noeud du maillage :

$$ \bar{\omega}_h^0(\mathbf{M}) = \prod_{i=1}^{2} (m_i - \frac{1}{2})h_i, \quad (m_i + \frac{1}{2})h_i [ , $$

- des ensembles en forme de croix :

$$ \bar{\omega}_h^1(\mathbf{M}) = \bigcup_{i=1,2} \bar{\omega}_h^0(M \pm \frac{1}{2} i), \\ \text{où l'on note } M \pm \frac{1}{2} i \text{ les points } M \pm \frac{1}{2} e_i, e_i \text{ étant le vecteur unitaire du i-ème axe de coordonnée,} $$

- une discrétisation de $\tilde{\Omega} = ]0, \mathbb{R}[x] ]0, a[$ :

$$ \tilde{\Omega}_h = \{ \mathbf{M} \mid \bar{\omega}_h^1(\mathbf{M}) \cap \Omega \neq \emptyset \}, $$

- les opérateurs de différences finies :
$$ \hat{\nabla}_r u(x) = \frac{u(x+\frac{1}{2}h_r) - u(x-\frac{1}{2}h_r)}{h_r}, $$

$$ \hat{\nabla}_z u(x) = \frac{u(x+\frac{1}{2}h_z) - u(x-\frac{1}{2}h_z)}{h_z}, $$

- les fonctions caractéristiques :

$\theta_h^M$ : fonction caractéristique du pavé $\bar{\omega}_h^O(M)$,

- un espace $V_h$ de fonctions en escalier : $v_h \in V_h$ si.

$$ v_h = \sum_{M \in \Omega_h} v_h^M \partial_h^M, \quad v_h^M \in \mathbb{R}, $$

- une discrétisation de $a(u,v)$ définie sur $V_h$ :

$$ (4.1) \qquad \left| \begin{array}{@{}l@{}} a_h(u_h, v_h) = \int_0^a \int_0^R (\hat{\nabla}_r u_h \hat{\nabla}_r v_h + \hat{\nabla}_z u_h \hat{\nabla}_z v_h) r dr dz - \\ \qquad \quad - R \int_0^a (w(u_h) - u_h) v_h dz, \end{array} \right. $$

avec $w(u_h)$ donné par (1.6),

- une discrétisation de $(fr,v)$ définie sur $V_h$ par

$$ (4.2) \qquad (\text{fr}, v_h) = \int_0^a \int_0^R \text{fr } v_h \, dr \, dz. $$

Le problème approché s'écrit alors

$$ (4.3) \qquad \left| \begin{array}{@{}l@{}} \text{Trouver } u_h \in V_h \text{ tel que} \\ \qquad a_h(u_h, v_h) = (\text{fr}, v_h), \quad \forall v_h \in V_h. \end{array} \right. $$
## 4.2 - Existence et unicité de la solution du problème approché

Comme pour le problème continu nous définissons sur $V_h$ une norme $\sqrt{b_h(v_h, v_h)}$ où $b_h$ est le produit scalaire défini par :

$$b_h(u_h, v_h) = \int_0^a \int_0^R (\hat{\nabla}_r u_h \hat{\nabla}_r v_h + \hat{\nabla}_z u_h \hat{\nabla}_z v_h) r dr dz + \int_0^a \frac{\partial w(u_h)}{\partial z} \frac{\partial w(v_h)}{\partial z} dz + w(u_h) (a) w(v_h) (a).$$

$\sqrt{b_h(v_h, v_h)}$ est une norme car :

$$\left. \begin{array}{l} \frac{\partial w(v_h)}{\partial z} = 0 \\ w(v_h)(a) = 0 \end{array} \right\} \Rightarrow w(v_h) = 0 \Rightarrow v_h \Big|_{\Gamma_R} = 0 \text{ d'après (1.10),}$$

et que d'autre part on sait que :

$$\int_0^a \int_0^R ((\hat{\nabla}_r v_h)^2 + (\hat{\nabla}_z v_h)^2) r dr dz = 0 \Rightarrow v_h = \text{constante}.$$

On vérifie comme en (2.1) que :

$$a_h(v_h, v_h) = b_h(v_h, v_h),$$

donc $a_h(v_h, v_h)$ est une forme bilinéaire coercive et continue (même démonstration qu'en 2.2) et par suite le problème (4.3) admet une solution unique d'après le lemme de Lax-Milgram.
4.3 - Le système linéaire associé au problème discret

La résolution du problème (4.3) conduit à la résolution du
système linéaire

(4.4)

$$
a_h(u_h, \hat{\theta}_h^M) = (\text{fr}, \hat{\theta}_h^M), \quad \forall M \in \Omega_h.
$$

Nous choisissons $h_r$ et $h_z$ de façon que la frontière de $\bar{\Omega}$
coïncide avec les lignes du réseau $\mathbb{R}$, et nous désignons les segments
de la frontière de $\bar{\Omega}$ par $\Gamma_o$, $\Gamma_a$, $\Gamma_R$, $\Gamma_{axe}$ et les noeuds du
maillage qui sont aux coins de $\Omega$ par A,B,C,D.

<figure><img src="image_5.png" /><figcaption>Fig. 4.1</figcaption></figure>

L'équation (4.4) s'écrit

(.5)

$$
\begin{align*}
& \int_0^a \int_0^R (\hat{\nabla}_r u_h \cdot \hat{\nabla}_r \theta_h^M + \hat{\nabla}_z u_h \cdot \hat{\nabla}_z \theta_h^M) r dr dz - R \int_0^a (w(u_h) - u_h) \theta_h^M dz = \\
& = \int_0^a \int_0^R f r \theta_h^M dr dz.
\end{align*}
$$
Mais pour $i=r$ et $i=z$ nous avons :

$$ \hat{\nabla}_i u_h = \left\{ \begin{array}{ll} \displaystyle \frac{u_h^{(M)} - u_h^{(M-h_i)}}{h_i} & \text{sur } \bar{\omega}_h^{\circ (M-\frac{h_i}{2})}, \\ \\ \displaystyle \frac{u_h^{(M+h_i)} - u_h^{(M)}}{h_i} & \text{sur } \bar{\omega}_h^{\circ (M+\frac{h_i}{2})}, \end{array} \right. $$

et

$$ \nabla_i \theta^M = \left\{ \begin{array}{ll} \displaystyle \frac{1}{h_i} & \text{sur } \bar{\omega}_h^{\circ (M-\frac{h_i}{2})}, \\ \\ \displaystyle \frac{1}{h_i} & \text{sur } \omega_h^{\circ (M+\frac{h_i}{2})}, \\ \\ 0 & \text{ailleurs.} \end{array} \right. $$

Donc (4.5) peut s'écrire :

$$ (4.6) \qquad \sum_{i=1,2} \left\{ \begin{array}{l} \displaystyle \frac{u_h^{(M)} - u_h^{(M-h_i)}}{h_i} \frac{1}{h_i} \int_{\Sigma \cap \bar{\omega}_h^o (M-\frac{h_i}{2})} r dr dz \\ \\ - \frac{u_h^{(M+h_i)} - u_h^{(M)}}{h_i} \frac{1}{h_i} \int_{\Sigma \cap \bar{\omega}_h^o (M+\frac{h_i}{2})} r dr dz \end{array} \right\} = \int_{\Gamma_R \cap \bar{\omega}_h^o (M)} fr dr dz, $$

ce qui, en posant :

$$ I(M \pm \frac{h_i}{2}) = \int_{\Sigma \cap \bar{\omega}} r dr dz, $$
nous donne l'équation au noeud $M \in \Omega_h$

$$ (4.7) \quad \left\{ \begin{aligned} & \frac{1}{h_r^2} \left[ I(M-\frac{h}{2}r) + I(M+\frac{h}{2}r) \right] + \frac{1}{h_z^2} \left[ I(M-\frac{h}{2}z) + I(M+\frac{h}{2}z) \right] \\ & - \frac{1}{h_r^2} I(M-\frac{h}{2}r) u_h(M-h_r) - \frac{1}{h_r^2} I(M+\frac{h}{2}r) u_h(M+h_r) - \\ & - \frac{1}{h_z^2} I(M-\frac{h}{2}z) u_h(M-h_z) - \frac{1}{h_z^2} I(M+\frac{h}{2}z) u_h(M+h_z) - \\ & - \int_{\Gamma_R} (w(u_h) - u_h) dz = \int_{\partial \Omega_h^o(M)} \text{fr } dr \, dz. \end{aligned} \right. $$

## 4.4 - Résolution du système linéaire

Pour résoudre le système linéaire (4.7) nous utilisons la méthode des relaxations successives (cf. VARGA [1], GLOWINSKI-LIONS-TREMOLIERES [1]).

### Descriptions de l'algorithme

Soit à résoudre le système linéaire

$$ Ax = \sum_{j=1}^{p} a_{ij} x_j = b_i, \quad i=1, \dots, p, \quad a_{ii} \neq 0. $$

On se donne $x^o \in \mathbb{R}^p$ et l'on calcule, pour $n=0, \dots, n$ le vecteur $x^{n+1} = x^n + \omega [x^{n+1} - x^n]$

$$ x_i^{n+\frac{1}{2}} = - \frac{1}{a_{ii}} \left( \sum_{j=1}^{i-1} a_{ij} x_j^{n+1} + \sum_{j=i+1}^{p} a_{ij} x_j^n - b_i \right) $$

$$ x_i^{n+1} = x_i^n + \omega \left[ x_i^{n+\frac{1}{2}} - x_i^n \right] $$
où $\omega$ est le paramètre de relaxation.

L'algorithme converge pour $0 < \omega < 2$ si la matrice A est par exemple à diagonale dominante et irréductible (1), ce que l'on vérifie aisément ici.

## Paramètre optimal de relaxation

La rapidité de la méthode est fonction du paramètre $\omega$, et passe par un maximum pour $\omega = \omega_{opt}$, valeur "optimale" du paramètre de relaxation. Cette valeur est donnée asymptotiquement par :

$$ \omega_{opt} = \frac{2}{1 + \sqrt{1 - \rho(\mathcal{L}_1)}} $$

où $\rho(\mathcal{L}_1)$ est le rayon spectral (la plus grande valeur propre) de la matrice dite de Gauss-Seidel :

$$ \mathcal{L}_1 = -(D+L)^{-1} U, $$

avec

$D$ = matrice diagonale extraite de $A$,

$L$ = matrice triangulaire inférieure extraite de $A$,

$U$ = matrice triangulaire supérieure extraite de $A$.

Une façon de calculer $\rho(\mathcal{L}_1)$ est d'appliquer la méthode des relaxations successives à la résolution de

$$ Au = 0, $$

<footnote>(1) A est irréductible s'il n'existe pas de matrice de permutation $P$ telle que $PA^T P$ soit bloc-triangulaire supérieur.</footnote>
avec $\omega = 1$ et $x^0 > 0$, et en calculant à chaque itération

$$ \rho_{p}^{n+1} = \min_{i=1, \dots, p} \frac{x_i^{n+1}}{x_i^n}, \quad \rho_{g}^{n+1} = \max_{i=1, \dots, p} \frac{x_i^{n+1}}{x_i^n}, $$

et, avec

$$ \omega_p^{n+1} = \frac{2}{1 + \sqrt{1 - \omega_p^{n+1}}}, \quad \omega_g^{n+1} = \frac{2}{1 + \sqrt{1 - \omega_g^{n+1}}}, $$

on peut montrer que $\rho_p^n \to \rho(\mathcal{L}_1)$, $\rho_g^n \to \rho(\mathcal{L}_1)$ et que les ? suites $(\omega_p^n)$ et $(\omega_g^n)$ convergent vers $\omega_{opt}$ en vérifiant

$$ \dots \rho_p^n \le \rho_p^{n+1} \le \dots \le \rho(\mathcal{L}_1) \le \dots \le \rho_g^{n+1} \le \rho_g^n \le \dots, $$

$$ \dots \omega_p^n \le \omega_p^{n+1} \le \dots \le \omega_{opt} \le \dots \le \omega_g^{n+1} \le \omega_g^n \le \dots, $$

en particulier

$$ \omega_{opt}^n = \frac{\omega_p^n + \omega_g^n}{2} \longrightarrow \omega_{opt} $$

<figure><img src="image_6.png" /></figure>
# 5 - <u>RESULTATS NUMERIQUES</u>

## 5.1 - <u>Exemples traités</u>

Nous avons pris pour domaine

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

avec

$R=1 \text{ et } a=1,2,4.$

Pour le coefficient de refroidissement on a considéré les valeurs.

$k=\frac{1}{10}, 1,10.$

On a pris $w(0)=0$, et pour pas de discrétisation

$h_r=1/20, h_z=a/20.$

<caption>Fig 5.1 : Nombre d'itérations de relaxation en fonction du paramètre $\omega$ pour $a=4, k=0,1$</caption>

## 5.2 - <u>Recherche du paramètre optimal de relaxation</u>

Afin de marquer l'importance de la détermination de $\omega_{opt}$ pour le problème considéré on a tracé sur les

Figures 5.1 et 5.2 le nombre d'itérations requises pour résoudre le problème en fonction de $\omega$.

<figure><img src="image_8.png" /><figcaption>Fig 5.2 : Nombre d'itérations de relaxation en fonction du paramètre $\omega$ pour $a=1, k=10$.</figcaption></figure>
La convergence des suites $\omega_p^n$ et $\omega_g^n$ vers $\omega_{opt}$ se fait comme le montre la Figure 5.3. La moyenne

$$ \omega_{opt} = (\omega_p^n + \omega_g^n)/2 \quad \text{donne une assez bonne estimation.} $$

<figure><img src="image_9.png" /><figcaption>FIGURE 5.3 : Convergence de l'algorithme de recherche du paramètre optimal de relaxation.</figcaption></figure>
## 5.3 - Résultats de l'algorithme des relaxations successives

Nous avons pris pour test d'arrêt de l'algorithme

$$DRM(n) \le 10^{-5},$$

où DRM(n) désigne la différence relative moyenne entre deux itérés définie par

$$DRM(n) = \frac{\sum_{i=1}^{p} |x_{i}^{n+1} - x_{i}^{n}|}{\sum_{i=1}^{p} |x_{i}^{n+1}|},$$

et nous obtenons pour les différents exemples traités.

<table><thead><tr><th></th><th>a</th><th>k</th><th>ω<sub>opt</sub></th><th>itérations</th><th>temps<sup>(1)</sup></th></tr></thead><tbody><tr><td>a=1</td><td>k=10</td><td>1,97</td><td>338</td><td>40s</td></tr><tr><td>a=1</td><td>k=1</td><td>1,937</td><td>143</td><td>20s</td></tr><tr><td>a=1</td><td>k=0,1</td><td>1,928</td><td>129</td><td>18s</td></tr><tr><td>a=2</td><td>k=10</td><td>1,973</td><td>424</td><td>46s</td></tr><tr><td>a=2</td><td>k=1</td><td>1,927</td><td>143</td><td>20s</td></tr><tr><td>a=2</td><td>k=0,1</td><td>1,908</td><td>109</td><td>16s</td></tr><tr><td>a=4</td><td>k=10</td><td>1,973</td><td>706</td><td>75s</td></tr><tr><td>a=4</td><td>k=1</td><td>1,933</td><td>141</td><td>20s</td></tr><tr><td>a=4</td><td>k=0,1</td><td>1,901</td><td>82</td><td>14s</td></tr></tbody></table>

<footnote>(1) temps d'exécution sur ordinateur CII 10070, recherche du paramètre optimal comprise.</footnote>
## 5.4 - Etude des solutions

On trouvera dans les figures 5.4 à 5.12 le tracé des lignes isothermes de la fonction température u.

<figure><img src="image_10.png" /><figcaption>FIGURE 5.4. Courbes isothermes de la fonction u pour a=1, k=10.</figcaption></figure>
<figure><img src="image_11.png" /><figcaption>FIGURE 5.5 : Courbes isothermes de la fonction u pour a=1, k=1.</figcaption></figure>
<figure><img src="image_12.png" /><figcaption>FIGURE 5.6. : Courbes isothermes de la fonction u pour a=1, k=1/10.</figcaption></figure>
<figure><img src="image_13.png" /><figcaption>FIGURE 5.7 : Courbes isothermes de la fonction u pour
$a=2$ , $k=10$.</figcaption></figure>
<figure><img src="image_14.png" /><figcaption>FIGURE 5.8 : Courbes isothermes de la fonction u pour a=2, k=1.</figcaption></figure>
<figure><img src="image_15.png" /><figcaption>FIGURE 5.9 : Courbes isothermes de la fonction u pour
a=2, k=0.1.</figcaption></figure>
<figure><img src="image_16.png" /><figcaption>FIGURE 5.10 : Courbes isothermes de la fonction u
pour a=4.,k=10.</figcaption></figure>
<figure><img src="image_17.png" /><figcaption>FIGURE 5.11 : Courbes isothermes de la fonction u pour a=4, k=1.</figcaption></figure>
<figure><img src="image_18.png" /><figcaption>FIGURE 5.12 : Courbes isothermes de la fonction u pour a=4, k=1/10.</figcaption></figure>
# ANNEXE

Au lieu de démontrer (2.3) (paragraphe 2.1) on peut établir l'existence en remplaçant (2.3) par

$$v \text{ est complet pour } \sqrt{b(v,v)}.$$

Dans cette annexe nous allons démontrer directement ce résultat.

Pour cela on considère une suite de Cauchy d'éléments $v_n$ de $V$ pour la norme $\sqrt{b(v,v)}$ :

$$ \int_{\Omega} \| \text{grad}(v_n - v_m) \|^2 d\Omega + \int_{\Gamma_p} \left( \frac{\partial [w(v_n) - w(v_m)]}{\partial z} \right)^2 d\Gamma_p + \Pi R \left( w(v_n)^2(a) - w(v_m)^2(a) \right) \to_\Omega, \quad \text{si } m, n \to +\infty $$

et nous allons montrer que la suite $v_n$ est convergente dans $V$. comme

$$ \int_{0}^{a} \left( \frac{d\varphi}{dz} \right)^2 dz + \varphi^2(a) $$

est une norme sur $H^1(\partial O, a)$ équivalente à la norme usuelle $\langle \cdot \rangle^{(1)}$ la suite $w(v_n)$ est une suite de Cauchy dans $H^1(\partial O, a)$ donc il existe $w \in H^1(\partial O, a)$ tel que

$$ w(v_n) \underset{\sim}{\longrightarrow} w \text{ dans } H^1(\partial O, a). $$

<footnote>(1) cf. NECAS [1].</footnote>
On fait alors l'hypothèse suivante (cf. LIONS-MAGENES [ 1] Tome 1, th. 8.3):

$H v \in V$

1) ∃ une application $\gamma$ linéaire continue et surjective :

$\gamma : V \longrightarrow H^2(\text{o}, a[),$

2) ∃ un relèvement R :

$$ \text{(5)} \quad \left[ \begin{array}{l} \text{R linéaire continu,} \\ \text{P : } H^2(\text{o}, a[) \to V \text{ tel que} \\ \gamma R(w) = w, \forall w H^2(\text{o}, a[). \end{array} \right. $$

Rappelons d'autre part (cf. CEA-GOULAQUIC [1], p. 123, Rem. 42) que

$H^1(\text{o}, a[) \subset H^2(\text{o}, a[).$

avec (4) on a alors

$w_n = w(v_n) \underset{w}{\longrightarrow} w \text{ dans } H^1(\text{o}, a[).$

D'après (5) (6) on peut considérer $R(w_n)$, $R(w)$ avec

$R(w_n) \to R(w)$ dans $V$.

Posons :

$C_n = P(w_n) - v_n.$

Nous allons montrer que $C_n \to C$ dans $V$. On aura alors, d'après (8)

$v_n \to C+R(w)$
d'où la convergence de la suite $v_n$, en sorte que $V$ sera complet pour
$\sqrt{b(v,v)}$.

Montrons donc que $C_n \to C$ dans $V$ pour la norme

$$
\int_{\Omega} |\operatorname{grad} v|^2 d\Omega + \int_{0}^{a} (w)^2 dz.
$$

Pour cela il suffit de voir que $C_n$ est une suite de Cauchy dans $V$
pour (9). On a

$$
\begin{align*}
& \int_{\Omega} |\operatorname{grad}(c_n - c_m)|^2 d\Omega + \int_{0}^{a} |\gamma_{c_n} - \gamma_{c_m}|^2 dz = \\
&= \int_{\Omega} |\operatorname{grad}(c_n + v_n - v_n - c_m + v_m - v_m)|^2 d\Omega + \int_{0}^{a} |w_n - v_n - (w_m - v_m)| dz \\
&\le 2 \int_{\Omega} |\operatorname{grad}(p(w_n) - R(w_m))|^2 d\Omega + 2 \int_{0}^{a} |\operatorname{grad}(v_n - v_m)|^2 d\Omega + \\
& \phantom{\le 2 \int_{\Omega} |\operatorname{grad}(p(w_n) - R(w_m))|^2 d\Omega + 2 \int_{0}^{a} |\operatorname{grad}(v_n - v_m)|^2 d\Omega +} \int_{0}^{a} \left| \frac{dw_n}{dz} - \frac{dw_m}{dz} \right|^2 dz.
\end{align*}
$$

Dans cette dernière inégalité,

- le premier terme tend vers 0 d'après (7),

- le second terme tend vers 0 d'après (2),

- le dernier terme tend vers 0 d'après (2),

la suite $c_n$ est donc une suite de Cauchy dans $V$ pour (9) ; elle est
donc convergente vers $c \in V$, d'où

$$
v_n \to c+R(w) \in V.
$$

L'espace $V$ est donc complet pour $\sqrt{b(v,v)}$.
BIBLIOGRAPHIE

[1] J.P. AUBIN : Approximation des espaces de distributions et des opérateurs différentiels. THESE, PARIS, 1966.

[1] J.F. BOURGAT : Analyse numérique du problème de la torsion élastoplastique. THESE 3eme Cycle, PARIS, 1971.

[1] J. CEA : Approximation variationnelle des problèmes aux limites. THESE, PARIS 1965.

[1] J. CEA GOULAOUIC : C4 Analyse fonctionnelle et numérique. RENNES, 1967 - 1968.

[1] G. DUVAUT J.L. LIONS : Les inéquations en mécanique et en physique. DUNOD, PARIS, 1972.

[1] R. GLOWINSKI J.L. LIONS R. TREMOLIERES : Méthodes numériques de résolution des problèmes d'inéquations en Mécanique et en Physique. DUNOD, PARIS, 1972 (à paraître).

[1] J.L. LIONS : Equations différentielles opérationnelles. SPRINGER-VERLAG, BERLIN, 1961.

[1] J.L. LIONS E. MAGENES : Problèmes aux limites non homogènes et applications. DUNOD, PARIS, 1968.

[1] J. NECAS : Les méthodes directes en théorie des équations elliptiques. MASSON, PARIS 1967.

[1] P.A. RAVIART : Sur l'approximation de certaines équations d'évolution linéaires et non-linéaires. THESE, PARIS, 1965.

[1] R. TREMOLIERES : Inéquations non symétriques. Note I.R.I.A. 1971.

[2] R. TREMOLIERES : Inéquations variationnelles : existence, approximation, résolution. THESE, PARIS, 1972.

[1] R.S. VARCA : Matrix Iterative Analysis. Prentice Hall 1962.