# Éléments finis et hydraulique numérique: 50 ans de la méthode des éléments finis 1955-2005

Gouri Dhatt, Philippe Sergent

► To cite this version:

Gouri Dhatt, Philippe Sergent. Éléments finis et hydraulique numérique: 50 ans de la méthode des éléments finis 1955-2005. 7e colloque national en calcul des structures, CSMA, May 2005, Giens, France. hal-01812998

HAL Id: hal-01812998

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

Submitted on 12 Jun 2018

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 CC0 1.0 - Universal - International License
# Éléments finis et hydraulique numérique

## 50 ans de la méthode des éléments finis 1955-2005

Gouri Dhatt*, Philippe Sergent**

*Laboratoire Roberval, Université de Technologie de Compiègne, BP 20529, 60205 Compiègne cedex

**Centre d'Etudes Techniques Maritimes Et Fluviales, 2, bd Gambetta BP 60039, 60321 Compiègne cedex

Les outils de modélisation numérique simulent le comportement d'un système physique en utilisant l'ordinateur. Ce domaine, issu de l'évolution informatique, est assez jeune. La méthode des éléments finis est un outil efficace utilisé en simulation numérique des problèmes industriels. Aujourd'hui nous trouvons plusieurs codes commerciaux qui intègrent à la fois des modules de Conception Assistée par Ordinateur (CAO) et des modules de simulation par la Méthode des Eléments Finis (MEF).

Dans les sections suivantes, nous présentons l'évolution de la méthode et son application aux problèmes d'hydraulique à surface libre.

### 1. Évolution de la MEF [1, 2, 3]

L'évolution de la MEF est associée étroitement au développement de l'informatique électronique. Nous distinguons trois périodes qui caractérisent les développements significatifs et conjoints des deux technologies.

#### 1.1 Période 1940-1955 : calculateur avec tube électronique

Le premier calculateur avec tube électronique ENIAC est apparu en 1942 occupant 100 m² et effectuant 300 opérations/seconde. En même temps est apparu l'ordinateur Z3 en Allemagne qui est développé par Konrad Zuse. Le calculateur électronique remplace ainsi l'homme calculateur. Le premier ordinateur commercial UNIVAC (universel automatic computer) est construit en 1951, suivi par les machines IBM 650 et 701. A partir de 1954, l'industrie aéronautique devient un client important pour ces nouveaux calculateurs.

Durant la même période, nous trouvons les développements importants en modélisation numérique. Notons les travaux de Courant, Von Neumann, Richtmeyer, Argyris, Southwall etc. en méthodes numériques associées aux problèmes nucléaires et mécaniques.

L'évolution rapide des calculateurs a créé une nouvelle génération de modélisateurs numériques qui sont très pragmatiques, empiriques et intuitifs. Il devient très difficile de suivre cette évolution pour les numériens classiques à cause de la grande complexité des problèmes résolus sur l'ordinateur. L'approche de
l'ingénieur a occupé une place importante dans le développement des nouvelles méthodes numériques. Il n'est plus possible de retarder la mise au point de nouvelles méthodes pour résoudre les problèmes industriels même si on maîtrise encore mal les critères de convergence et la présence des erreurs dans un sens mathématique classique.

## 1.2 Période 1955-1965 : les transistors

Le tube électronique est remplacé par les transistors qui sont plus petits, plus fiables et consomment moins d'énergie. Nous trouvons plusieurs ordinateurs sur la marché : UNIVAC (1956), FORTRAN en 1957 sur IBM 710, DEC-PDP (1960), IBM 7030,1410 et 1620. Le système de stockage sur disque est déjà en place.

Les travaux de chercheurs chez Boeing présentent les premières notions de la méthode des éléments finis. L'article de Turner, Martin, Clough et Topp [4] en 1956 caractérise la rigidité élémentaire d'un élément triangulaire d'élasticité plane. En effet ce travail a été effectué chez Boeing en 1952-1953. L'article de Turner [5] en 1959 dans un congrès AGARD en Allemagne présente la notion de la rigidité élémentaire et d'assemblage. Clough [6] introduit le terme éléments finis dans le cadre d'une formulation discrète du problème d'élasticité plane. En 1963, Wilson présente sa thèse de doctorat à Berkeley sur la formulation éléments finis en élasticité bi- dimensionnelle. Le premier congrès [7] sur la méthode des éléments finis est organisé en 1965 à Wright Patterson Air Force Base Dayton, Ohio.

## 1.3 Période 1966-... : VLSI, ULSI, microprocesseurs

La technologie des circuits intégrés à grande densité et la fabrication de microprocesseurs ont introduit une véritable révolution : des ordinateurs de plus en plus petits ; une vitesse de calcul de plus en plus rapide. Notons quelques réalisations importantes : IBM système 360 en 1964 de 3ème génération, CDC 7600 en 1968 avec 40 megaflops par Cray, Microprocesseur INTEL 8080 en 1972, CRAY-1 en 1976, APPLE-1 en 1976, VAX 11/780 en 1978, IBM-PC avec architecture ouverte en 1980. Nous connaissons bien la suite des développements informatiques depuis 1980.

Ainsi, très rapidement la MEF est devenue une méthode industrielle en industrie aéronautique, industrie automobile, mécanique, fluide, thermique etc. En 1965, la NASA prend l'initiative de développer un code de calcul aux éléments finis en mécanique des structures utilisant la formulation déplacement. Zienkiewicz et Cheung [8] publie un premier livre sur la méthode des éléments finis en 1967. Notons également les livres de Robinson [9], Przemieniecki [10] et Strang et Fix [11]. Nous laissons la suite aux lecteurs.

Aujourd'hui la MEF est une méthode industrielle et l'on trouve plusieurs codes de qualité remarquable. À titre d'exemples, voici quelques codes utilisés pour le calcul de problèmes linéaires ou non linéaires en mécanique, thermique et fluide:

- ABAQUS, ADINA, ANSYS, CASTEM, CASTOR, CATIA, FIDAP, LS-DYNA, MARC, NE/NASTRAN, RADIOS, SAMCEF, REFLUX (hydraulique), TELEMAC (hydraulique),...
## 2. Démarche éléments finis [12, 13]

Après 50 ans d'évolution, la méthode des éléments finis est utilisée couramment dans tous les milieux industriels. Elle permet aux concepteurs de bien comprendre le comportement des systèmes physiques, d'étudier l'influence des différents paramètres et éventuellement d'optimiser la conception du système. Nous identifions différentes composantes de la démarche des éléments finis pour une meilleure compréhension et utilisation de la méthode par l'ingénieur.

### 2.1 Écriture physique

L'objectif de la méthode des éléments finis est de simuler le comportement d'un système physique en utilisant l'ordinateur. La description correcte d'un système physique constitue une première étape essentielle de la modélisation numérique. L'ingénieur doit utiliser tout son savoir-faire pour identifier avec soin et précision les différentes composantes caractérisant le comportement d'un système physique. Pour bien utiliser un code commercial, il est important de décrire correctement le problème physique à étudier. Pour cela nous présentons les points suivants qui peuvent aider à la préparation du modèle :

- Choix de phénomènes : Pour chaque type de problème, il faut identifier si le problème est linéaire ou non-linéaire, stationnaire ou non, turbulent ou laminaire etc.

- Choix de variables de description $u(x,y,z,t)$ : Il faut identifier les variables qui caractérisent le comportement du système. À priori il semble que ce choix est trivial. Par exemple pour un problème d'élasticité, on choisit les déplacements et les contraintes, pour les fluides on choisit les vitesses, la pression, les contraintes visqueuses etc. Le choix devient complexe si une structure est mince ou épaisse, il faut choisir si le comportement est de type poutre, plaque, coque ou solide. Pour un problème couplé en hydro-sédimentologie, peut-on supposer que le changement du fond ne modifie pas l'hydrodynamique ? Pour un problème thermo-élasto-plastique, comment choisir le couplage élastique et thermique ?

- Choix de la géométrie : Pour certains problèmes, la géométrie est bien définie à priori. Mais pour des problèmes d'écoulement naturel, la géométrie n'est pas bien définie. Il est possible aussi que la géométrie soit variable à différents instants. En plus, il faut choisir si le comportement est uni-, bi- ou tri- dimensionnel, soit : $u(x,y,z,t)=u(x,t)$, $u(x,y,t)$ ou $u(x,y,z,t)$. De plus, il faut choisir si le comportement est stationnaire ou non.

- Loi de conservation : Le comportement est régi par les lois universelles de conservation : soit la conservation des efforts (équilibre), la conservation de la quantité de mouvement, la conservation de la masse, la conservation du flux thermique, etc.

- Loi de comportement : Il est important de choisir correctement la loi de comportement : soit la loi contraintes - déformations, la loi contraintes visqueuses - taux de déformations. Pour les grandes déformations, la loi de comportement implique plusieurs approximations. Pour les fluides turbulents, il faut modéliser les contraintes de Reynolds représentant la turbulence par une loi de comportement, c'est un choix assez difficile.
- Conditions aux limites et conditions initiales : Ici encore, l'ingénieur doit bien comprendre la physique de son problème et faire un choix correct des conditions aux limites et des conditions initiales.

- Sollicitations volumiques et surfaciques.

## 2.2 Écriture mathématique

On associe un modèle mathématique à la description physique qui est retenue par l'ingénieur. Ce modèle est écrit sous forme de relations aux dérivées partielles (RDP, appelées écriture forte) ou sous une forme intégrale en utilisant des fonctions test (appelée écriture faible). L'écriture faible est aussi appelée forme intégrale, écriture de type Galerkine, écriture variationnelle, écriture de type résidus pondérés etc. La MEF transforme l'écriture faible en écriture algébrique.

L'erreur de modélisation est associée à la différence entre la solution du modèle mathématique et le comportement du système réel. Remarquons que le modèle mathématique est associé à la description physique retenue après plusieurs simplifications discutées dans le paragraphe précédent. En effet, l'erreur de modélisation représente la différence du comportement entre le modèle retenu et la physique réelle. Dans tous les cas réels, on ne connaît ni le comportement réel, ni la solution du modèle mathématique.

## 2.3 Écriture algébrique

Comme on ne peut pas trouver la solution du modèle mathématique, on associe un modèle algébrique à l'écriture faible en utilisant l'approximation par éléments finis, soit :

$$V = \sum_{e} V^{e} \qquad W = \sum_{e} W^{e} = 0$$

Sur chaque élément, une approximation admissible de type élément fini est :

$$u = \sum_{i} N_{i} u_{i} \quad \text{et} \quad \delta u = \sum_{i} N_{i} \delta u_{i} \qquad u : \text{ensemble des variables cherchées}$$

L'écriture discrétisée W s'écrit : $W = \langle \delta U \rangle = \{ R(u_n) \} = \sum \langle \delta u \rangle = 0$

L'erreur de discrétisation est associée à la différence entre la solution du modèle mathématique et la solution $u_n$ du système discret $R=0$.

La construction du système R et la résolution du problème $R=0$ sont obtenues sur ordinateur en utilisation les techniques numériques de résolution directe ou itérative.

## 2.4 Écriture informatique

Pour obtenir la solution du problème physique, il faut effectuer une mise en forme informatique de toutes les opérations :

- la saisie des données : géométrie, propriétés physiques, conditions aux limites, sollicitations, conditions initiales.

- le maillage : maillage du domaine selon les types d'éléments choisis.
- l'opération de calcul aux éléments finis : matrices élémentaires, assemblage, schémas d'intégration temporelle, schémas de résolution directe ou itérative et la visualisation des résultats.

Un code commercial aux éléments finis effectue toutes ces opérations d'une manière automatique. Souvent un code de calcul aux éléments finis est couplé avec un système de CAO pour obtenir la description de la géométrie du problème.

L'erreur de convergence (soit erreur numérique) est associée à la différence entre la solution du modèle discret et la solution obtenue par le code. Cette erreur inclut l'erreur de convergence d'un schéma et l'erreur d'arrondi.

Nous pouvons résumer la démarche de la méthode des éléments finis par la figure 1.

<figure><img src="image_2.png" /><figcaption>Figure 1. Démarche de la méthode des éléments finis.</figcaption></figure>

Pour obtenir une solution convergente du modèle aux éléments finis, il faut vérifier les conditions suivantes :

- **Consistance** : L'approximation par éléments finis doit satisfaire les conditions nécessaires de continuité inter éléments. De plus la précision de l'approximation doit assurer une représentation constante de tous les termes de W. Cette condition est nécessaire mais pas suffisante.

- **Stabilité** : Le modèle discrétisé doit être stable conduisant à une solution sans oscillations (absence de modes parasites) et sans blocage. Par exemple, pour une poutre ou une plaque épaisse, il faut éviter le blocage par cisaillement. On utilise une formulation qui réduit le rang de la matrice de cisaillement, soit par une approche mixte, par une projection de type collocation par ligne DKT ou par intégration sélective. Pour des problèmes en fluide incompressible par formulation pénalité, il faut éviter le blocage d'incompressibilité. Le rang du terme d'incompressibilité doit être petit par rapport au rang de la matrice élémentaire. L'intégration numérique des quantités élémentaires doit assurer un rang correct. Pour des schémas temporels, nous utilisons la technique spectrale de Neumann pour vérifier la stabilité.

- **Positivité** : Il est utile de vérifier que la discrétisation éléments finis respecte la condition de positivité. Cette condition est couramment utilisée en volumes finis pour construire des schémas stables. La condition est basée sur le raisonnement physique, soit pour une discrétisation donnée : *Il faut s'assurer que la valeur u du nœud i est positive si les valeurs de u de tous les noeuds autour du noeud i sont positives.* Il est facile de vérifier si une discrétisation spatiale ou spatio-temporelle respecte la positivité. Le problème peut être uni-, bi- ou tri dimensionnel. Si un
schéma ne vérifie pas la positivité, on peut rencontrer des oscillations parasites. Il est possible qu’un schéma stable ne soit pas positif.

- Exemples stationnaires :

Soit le problème de convection diffusion : $a u_x + k u_{xx} = 0$

Le schéma discret au nœud i est : $c_{i-1} u_{i-1} - c_i u_i + c_{i+1} u_{i+1} = 0$

$$ \text{avec : } c_{i-1} = \frac{(k+a) \frac{dx}{dx}}{2}, \quad c_{i+1} = \frac{(k-a) \frac{dx}{dx}}{2}, \quad c_i = \frac{2k}{dx^2}. \text{ Pour un } c_i > 0, \text{ il faut que } c_{i+1} \text{ et } c_{i-1} > 0. $$

Si $k, a > 0$, il faut : $a \frac{dx}{k} = \text{Péclet} < 2$ pour la positivité; $dx = \text{taille d'élément}$

Soit un problème de dynamique forcée de pulsation $\omega$ :

$$ a u + k u_{,xx} = 0 \text{ avec } a = \rho \omega^2 \text{ et } k = \text{module E} $$

Le schéma discret ou nœud i est avec masse consistante :

$$ c_{i-1} u_{i-1} - c_i u_i + c_{i+1} u_{i+1} = 0 $$

$$ \text{avec : } c_{i-1} = \frac{k+a}{6} \frac{dx^2}{dx} \qquad c_{i+1} = \frac{k+a}{6} \frac{dx^2}{dx} \qquad c_i = \frac{2k-4a}{6} \frac{dx^2}{dx} $$

Pour la positivité, il faut : $dx^2 < 3\frac{k}{a}$, soit 6 éléments par longueur d'onde pour assurer une solution sans oscillation.

- Exemples non-stationnaires :

Soit le problème thermique temporel

$$ u_t - k u_{,xx} = 0 $$

Le schéma explicite au nœud i avec masse diagonale est :

$$ u_i^{t+dt} = c_{i-1} u_{i-1}^t + c_i u_i^t + c_{i+1} u_{i+1}^t $$

$$ \text{avec } c_{i-1} = c_{i+1} = k \frac{dt}{dx^2} \quad \text{et} \quad c_i = 1 - 2k \frac{dt}{dx^2}. \text{ Pour la positivité, il faut } k \frac{dt}{dx^2} < 0.5 $$

Pour les problèmes thermiques, il est courant de choisir une matrice masse diagonale pour assurer la positivité d'un schéma implicite.

### 3. Hydraulique numérique [14, 15, 16, 17]

Le domaine de l'hydraulique numérique est issu de l'évolution informatique. L'objectif principal est d'obtenir une simulation numérique d'écoulement à surface libre en mer, dans les fleuves et rivières. Nous pouvons également coupler l'écoulement de surface avec l'écoulement souterrain en milieu poreux. Ce domaine couvre également les problèmes de type hydro sédimentaire et de pollution d'eau. Les outils de simulation numérique en hydraulique sont couramment utilisés par les concepteurs pour comprendre les problèmes d'inondation, d'érosion sédimentaire et de pollution. Dans ce travail nous nous limitons aux problèmes d'inondation rurale et urbaine induite par le débordement des rivières et par ruissellement. Il est possible d'étendre ces modèles à l'étude des marées et tsunamis. Les modèles d'hydraulique numérique sont obtenus à partir des relations de Navier Stokes incompressible et turbulent.

#### 3.1 Modèles de type Saint-Venant
Les relations générales de Navier Stokes sont :

Conservation de la masse div( u ) = 0

Conservation de la quantité de mouvement (forme conservative)

$$
\mathbf{u}_t + \text{div} ( \mathbf{u} \mathbf{u} + \mathbf{I} \cdot \mathbf{p} / \rho - \tau / \rho ) - \mathbf{g} = 0
$$

avec **u**=(u, v, w) composantes cartésiennes de la vitesse, verticale selon z

$$
\mathbf{u} \cdot \mathbf{u} = (u.u + u.v + u.w, u.v + v.v + v.w, u.w + v.w + w.w)
$$

p : pression et τ : contraintes visqueuses turbulentes

La surface libre est définie par : h<sub>t</sub> + div (H U) = 0 avec h : surface libre, H :

profondeur, **U**: (U, V) vitesses moyennes sur la profondeur

<figure><img src="image_3.png" /><figcaption>**Figure 2.** *Variables utilisées pour les écoulements à surface libre.*</figcaption></figure>

En intégrant les équations de Navier Stokes selon la verticale, nous obtenons le
modèle de Saint-Venant avec les approximations suivantes :

- Longueur d'onde h/L << 1 et pression hydrostatique. L : longueur d'onde.

- Fond avec faible pente.

L'écriture conservative est (x,y):

$$
(H,U)_t + \text{div}(H U \cdot \mathbf{U} - H \frac{\tau}{\rho}) + g \text{ grad } h - \frac{\tau_f}{\rho} - \frac{\tau_s}{\rho} - f_c + \frac{\text{grad } P_a}{\rho} = 0
$$

avec **U** = (U, V) vitesse moyenne , τ : contraintes visqueuses moyennes,

τ<sub>f</sub> = contrainte au fond de Chézy :

$$
\text{Chézy : } \rho g \text{ abs(U)} \mathbf{U} / C^2, \text{ Manning-strickler : } \rho g n^2 \text{ abs(U)} \mathbf{U} / H^{1/3}
$$

La conservation de la masse est : h<sub>t</sub> + div (H U) = 0

Ces relations sont résolues par la méthode des éléments finis en utilisant un élément T6 composé de 4 sous-triangles. Pour assurer la stabilité, on utilise une approximation de la vitesse sur chaque sous triangle et une approximation pour h sur le triangle de base. La solution numérique est obtenue à l'aide d'un schéma implicite de type Euler avec une linéarisation de Newton. Pour des problèmes avec front mobile, nous utilisons un schéma explicite de type Kawahara ou Lax-Wendroff. Pour étudier l'inondation urbaine, nous avons développé des modèles uni- et bidimensionnels représentant le réseau des rues. Le modèle uni dimensionnel est obtenu en intégrant les relations de Navier Stokes sur la section de chaque rue. On utilise le schéma explicite de Lax-Wendroff avec une capture de choc pour obtenir l'évolution des fronts mobiles.
# 4. Résultats Numériques

## 4.1 Inondations rurales à Remiremont (Département des Vosges)

De nombreuses zones urbanisées sont protégées des crues par des digues. Même si elles ne sont plus touchées par les inondations directes des rivières, elles peuvent être le sujet de remontées de nappe. Poligot-Pitsch et al. [18] ont étudié les remontées de la nappe alluviale dans le val d'Avaray derrière les digues de la Loire. Ils ont montré que le ruissellement de surface pouvait influencer les écoulements souterrains et augmenter l'exfiltration. Le même phénomène peut exister pour des remblais routiers et a été étudié pour le compte de la Direction Départementale de l'Equipement des Vosges (DDE 88) dans le cadre d'une collaboration entre le CETMEF, l'UTC et le BRGM sur le site de Remiremont situé le long de la Moselle à l'aval du confluent avec la Moselotte. La crue de février 1990 a occasionné des inondations et des dommages importants dans les vallées de la Moselle et de la Moselotte. La crue a traversé Remiremont après la construction de la déviation de la RN 66 (appelée rocade) et de sa jonction avec la CD 417 (appelée pénétrante). Ces deux voies routières sont implantées dans le lit majeur de la Moselle à l'aval du confluent avec la Moselotte. Les terrains compris entre la RN 66 et la voie SNCF actuellement démantelée se composent de trois casiers principaux (voir figure 3). Dans le premier casier, situé en amont, sont implantées diverses usines. Le second casier en aval du remblai de la pénétrante comprend un plan d'eau. A l'arrière de la voie SNCF qui forme un écran, se trouve le troisième casier où sont implantées diverses constructions, notamment une filature qui a subi de graves dommages lors du passage de la crue de 1990. La filature n'a pas été directement inondée par les eaux du lit actif de la Moselle, mais plutôt par un effet de remplissage en série de ces casiers. L'eau qui a rempli le casier de la filature s'est écoulée en résurgence par les réseaux d'assainissement. L'étude a consisté d'abord à évaluer l'impact hydraulique successif des constructions des deux voies routières, la RN66 et la CD 417, implantées dans le lit majeur du confluent Moselle-Moselotte. Cette étude hydraulique est en régime permanent avec des calculs réalisés sur la base des débits au maximum de la crue du 15 Février 1990. Les scénarii étudiés sont donc les suivants : (1) Rivière sans aménagements routiers ; (2) Rivière avec rocade ; (3) Rivière avec rocade et pénétrante.

<figure><img src="image_4.png" /><figcaption>**Figure 3.** Plans de situation des voies routières et des casiers.</figcaption></figure>
Divers aménagements ont été réalisés par la Direction Départementale de l'Equipement des Vosges après la crue de 1990 (voir figure 4). Le premier casier est notamment protégé désormais d'un écoulement direct des eaux de la Moselle. Un calcul a donc été ensuite mené en régime transitoire pour donner la dynamique de la crue de 1990 ainsi que l'impact hydraulique des travaux d'aménagements contre les crues qui ont suivi l'épisode d'inondations. (4) Rivière avec aménagements routiers et aménagements actuels de protection. Comme les casiers à enjeux sont désormais protégés d'un écoulement direct de la crue par un remblai routier, il a enfin été décidé d'estimer aussi les inondations de surface par remontées de nappe.

<figure><img src="image_5.png" /><figcaption>**Figure 4.** Aménagements réalisés après la crue de 1990.</figcaption></figure>

Les modèles hydrauliques sont des modèles emboîtés. Une première modélisation filaire (1D) repose sur trois biefs, deux biefs couvrant le lit principal de la Moselle en amont et aval de la confluence avec la Moselotte et le troisième bief décrivant la Moselotte. Le MNT est composé de 270 profils en travers, dont 18 profils de ponts et 19 profils de barrages en lit mineur. Le secteur d'étude du modèle bidimensionnel (2D) couvre la confluence de la Moselle et de la Moselotte. Le MNT de cette configuration possède 29 752 nœuds et 14 139 éléments. Le modèle 2D est le modèle REFLUX® qui résout l'équation de Saint-Venant par la méthode des éléments finis.

Le modèle numérique a été calé sur de nombreuses laisses de crue à la fois pour les modélisations 1D et 2D. En régime transitoire, le limnigramme à l'échelle de crue de Remiremont a également été validé (voir figure 5).

<figure><img src="image_6.png" /><figcaption>**Figure 5.** Calage du modèle hydraulique à l'échelle de Remiremont.</figcaption></figure>
Les calculs ont d'abord montré que la pénétrante était responsable d'une surélévation de 50 cm du lit majeur de la Moselle avec une progression de la crue vers l'aval dans le casier 1 qui va buter contre le remblai de la pénétrante. Cela provoque alors une surélévation de 80 cm à cet endroit. Les aménagements de protection contre les inondations dans le casier 1 ont ensuite fait perdre 58 000 m³ de volume de stockage et provoquent une surélévation inférieure à trois centimètres en lit majeur de la Moselle. Le casier 1 étant maintenant isolé du lit majeur, il reste à examiner les remontées de nappe éventuelles qui pourraient se produire lors d'une crue exceptionnelle. C'est le sujet des paragraphes suivants.

Le calcul de l'écoulement souterrain est réalisé à l'aide du modèle MARTHE® du BRGM (Thiéry [19]). Ce modèle utilise une discrétisation par volumes finis. Il a été utilisé ici en système monocouche avec une perméabilité de $5 \cdot 10^{-4} \text{ ms}^{-1}$ soit une transmissivité de $1 \cdot 10^{-2} \text{ m}^2\text{s}^{-1}$ pour une épaisseur de nappe de 20 m. Le substratum est en effet supposé situé 20 m en dessous du terrain naturel. Le maillage comprend 6 098 mailles carrées de 10 m de coté. Pour les calculs en régime transitoire, le coefficient d'emmagasinement libre SL est pris égal à 10 %.

Le protocole de couplage est le suivant à la suite d'un calcul de REFLUX®. Dans les mailles pour lesquelles REFLUX® a calculé une hauteur nulle, un index de suintement est introduit dans MARTHE®. Si le niveau piézométrique tend à passer au-dessus de la topographie, le niveau est maintenu à la topographie et un débit de débordement est calculé par MARTHE®. Le débit qui déborde est injecté au pas de temps suivant dans la maille correspondante de REFLUX®. Dans les mailles pour lesquelles REFLUX® a calculé une hauteur positive, la charge est imposée dans MARTHE® qui calcule alors un débit positif (exfiltration) ou négatif (infiltration). Ce débit sera lui aussi injecté dans REFLUX® au pas de temps suivant.

<figure><img src="image_7.png" /><figcaption>**Figure 6.** Hauteur d'eau au-dessus du terrain naturel après 50 heures.</figcaption></figure>

Les calculs sont réalisés sur cinquante heures entre le 14 et le 16 Février sur la configuration (4) avec un pas de temps de couplage de 30 minutes. Sur la figure 6 sont synthétisés les différents résultats de cette étude : d'abord l'inondation de 1990 qui a touché globalement les deux casiers et principalement le casier 1 à enjeux, ensuite la situation actuelle où le casier 2 joue seul le rôle de zone de stockage, enfin l'effet des remontées de nappe. Un effet bénéfique des remontées de nappe est un stockage plus fort dans le casier 2. Le stockage passe en effet de 7 000 m3 à 22 000 m3. En revanche de légères remontées sont aussi notées juste derrière les remblais dans le casier 1 qui est occupé par des usines. On peut estimer que ce constat
pourrait être aggravé pour une simulation plus longue ou pour une crue centennale du type de celle de 1947.

## 4.2 Inondations urbaines à Marseille et Bordeaux

Les tragédies récentes dans les zones urbanisées comme Bordeaux (1982), Nîmes (1988) (Al Mikdad [17]) ou Marseille (2000) ont amené les gestionnaires et techniciens à se pencher sur le problème. Les systèmes de collecte des eaux de pluie sont souvent dimensionnés pour des événements décennaux alors que des événements plus graves génèrent du ruissellement dans les rues et endommagent malheureusement les biens et provoquent parfois des morts. Le travail réalisé a pour objectif d'améliorer les modèles d'écoulement en ville à travers un projet rassemblant différents organismes du Ministère de l'Equipement ainsi que des scientifiques. Des applications ont été réalisées avec les modèles sur le bassin-versant du quartier du Vieux-Port à Marseille (figure 7) et sur le quartier de La Grenouillère à Bordeaux (figure 8).

<figure><img src="image_8.png" /><figcaption>**Fig 7.** Submersion des rues au début de la crue (gauche) et après 45 min. (droite).</figcaption></figure>

<figure><img src="image_9.png" /><figcaption>**Fig. 8.** Submersion des rues au début de la crue (gauche) et après 150 min. (droite).</figcaption></figure>

Ce travail vise en particulier à tester un modèle bidimensionnel résolvant les équations de Saint-Venant et qui simule les fortes vitesses de courant dans les rues avec zone d'urbanisation dense durant des crues éclair. Le modèle prend en compte les changements de régime avec capture du ressaut hydraulique à l'aide d'une viscosité artificielle. Il inclut également les apports des précipitations, les écoulements sur front sec ainsi que différents types de frontière ouverte.
Ce modèle est testé sur plusieurs cas-tests : l'écoulement sur un fond en forme de cloche, le ressaut hydraulique oblique et des carrefours à angle-droit.

Des outils de pré-traitement sont utilisés afin de créer un maillage adapté à l'écoulement le long des rues et des carrefours. Ces outils extraient les frontières des rues d'une base de données urbaine, déterminent l'axe de la rue à l'aide d'un diagramme de Voronoi et produisent un maillage régulier des rues. Un mailleur de Delaunay est utilisé pour les carrefours complexes tandis qu'un maillage régulier est conservé pour les carrefours à angle-droit.

Références

[1] Carlson B., Burgess A., Miller C., Evolution of computer and computer
technology, Timeline research, 1995

[2] Jones international Computers : History and development,
www.digitalcetury.com/encyclo/update/comp_hd.html, 1999

[3] Felippa C., A Historical outline of matrix structural analysis, CU-CAS-00_13,
Univ of Colorado, Boulder, 2000

[4] Turner M.J., Clough R.W., Martin H.C., Topp L.J. Stiffness and deflection analysis complex structures, J. Aero. Sc. 23, 805-823,1956

[5] Turner A.J., Direct stiffness method of structural analysis, structural and material panel paper, AGARD Aachen, 1959

[6] Clough R.W., The finite element in plane stress analysis Proc. 2nd ASCE conf. on electronic computation, Pittsburgh Pa, 1960

[7] Proceedings of conference on Matrix methodes in structural analysis, Wright Patterson Air Force Base, Dayton, Ohio, 1965, 68, 71

[8] Zienkiewicz O.C., Cheung Y.K., The finite element method in structural an solid mechanics, Mcgraw Hill 1967

[9] Robinson J., Structural matrix analysis for engineer, Wiley 1966

[10] Przemieniecki J.S., Theory of matrix structural analysis, Mcgraw Hill 1968 (Dover 86)

[11] Strang G., Fox O.J., Analysis of finite element method, Prentice hall 1973

[12] Dhatt G., Touzot G, Lefrançois E., Méthode des éléments finis, Édition Hermès-Lavoisier, 2005

[13] Batoz JL., Dhatt G., Modélisation des structures par éléments finis, Vols 1,2,3, Hermès 1990

[14] Tanguy J.M., De la goutte de pluie jusqu'à la mer, thèse HDR, UTC, 2003

[15] Sergent P., De l'acoustique à l'hydraulique numérique, thèse HDR, UTC, 2004

[16] Zhang B., Modélisation d'écoulements à surface libre avec fronts mobiles par éléments finis, thèse PhD, UTC, 1992

[17] Al Mikdad 0., Modélisation des écoulements en réseau maillé et ramifié, thèse
PhD, UTC, 2000

[18] Poligot-Pitsch.S., Sergent.P., Dhatt G., Modélisation tridimensionnelle des remontées de nappe et inondations au val d'Avaray. Revue Européenne des Eléments Finis, 12 (2/3), 235-255, 2003

[19] Thiéry D., Modélisation des aquifères complexes - Prise en compte de la zone non saturée et de la salinité. Calcul des intervalles de confiance. Revue Hydrogéologie 4, 325-336, 1993