# Une application de la méthode LATIN au calcul multirésolution de structures non linéaires

Pierre-Alain Boucard, Pierre Ladevèze

► To cite this version:

Pierre-Alain Boucard, Pierre Ladevèze. Une application de la méthode LATIN au calcul multirésolution de structures non linéaires. Revue Européenne des Éléments Finis, 1999, 8 (8), pp.903-920. 10.1080/12506559.1999.10511417. hal-01698036

HAL Id: hal-01698036

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

Submitted on 7 Jan 2020

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 CC BY 4.0 - Attribution - International License
# Une application de la méthode LATIN au calcul multirésolution de structures non linéaires

Pierre-Alain Boucard — Pierre Ladevèze

LMT - Cachan

(E.N.S. de Cachan / C.N.R.S. / Université Paris VI)

61, Avenue du Président Wilson, F-94235 Cachan

boucard@lmt.ens-cachan.fr, ladeveze@lmt.ens-cachan.fr

**RÉSUMÉ.** Le traitement numérique des problèmes non linéaires complexes (emboutissage, viscoplasticité cyclique, crash...) nécessite, même sur les ordinateurs les plus performants, des temps de calcul prohibitifs. De plus, lors d'études de structures nécessitant de multiples résolutions (études paramétriques, structures avec défauts, ...), il faut effectuer pour chaque jeu de données un calcul complet, même si les problèmes à résoudre sont similaires. Le choix d'une méthode de calcul efficace et adaptée à des résolutions multiples est essentiel. La plupart des algorithmes employés reposent sur une démarche incrémentale et sont de ce fait mal adaptés. La méthode LATIN est une méthode non incrémentale, en rupture avec ces techniques classiques et semble une réponse prometteuse à ce type de problèmes.

**ABSTRACT.** The numerical solution process for complex, non-linear problems requires a computation time that winds up becoming prohibitive. Moreover, with respect to those structural analyses involving multiple solution sequences (parametric studies, structures with defect, etc.), each series of data necessitates performing a full calculation, even if the problems to be solved do resemble one another. The goal of the work presented herein is to develop a strategy that is well-suited to multiple-solution problems. Such a strategy is to be based on the LATIN method.

**MOTS-CLÉS :** mécanique des structures, calculs non linéaires, résolutions multiples.

**KEY WORDS :** Structural mechanics, nonlinear computations, multiple solution.
# 1. Introduction

Lorsque l'étude d'une structure nécessite de multiples résolutions (études paramétriques, analyses probabilistes, structures avec défauts, ...) chaque jeu de données requiert un calcul complet. On résout alors de nombreux problèmes très similaires. Le choix d'une méthode de calcul efficace et adaptée est alors primordial. La méthode LATIN [LAD 91] est une approche non incrémentale et semble à ce titre prometteuse, car plus adaptée de par ses principes que la plupart des algorithmes incrémentaux classiques (éléments finis stochastiques et méthodes de perturbation [BEN 98] [MAC 97], tirages de Monte-Carlo [SHI 72]). Le but de ce travail est de montrer les possibilités de la méthode LATIN en matière de multirésolution.

La méthode LATIN a deux caractéristiques majeures : elle est globale et itérative sur l'espace-temps Ω x [0, T], l'intervalle d'étude pouvant être aussi grand qu'on le souhaite. Elle peut donc réutiliser la solution obtenue pour le problème de base pour initialiser le processus itératif conduisant à la solution du problème perturbé. Cette propriété est l'un des atouts de la méthode LATIN dans le calcul multirésolution de structures non linéaires.

Un autre point également important dans l'efficacité de la méthode concerne la résolution à chaque itération du problème linéaire global que l'on obtient ; on utilise le troisième principe de la méthode qui, pour des chargements quasi-statiques, introduit des approximations « temps-espace » à fort contenu mécanique. Ces approximations « temps-espace » sont des sommes de produits de fonctions du temps par des fonctions de la variable d'espace, approximations qui font exploser le problème global sur l'espace-temps en une série de problèmes en espace et une série de problèmes en temps. Le coût de construction des fonctions de la variable d'espace étant relativement élevé, ces fonctions sont stockées et réutilisées. Pour traiter le problème perturbé, plutôt que d'initialiser avec la solution du problème de base, on exploite les fonctions de la variable d'espace obtenues lors de sa construction. Dans un certain nombre de cas, cette approximation est suffisante. Sans quoi, la procédure itérative est poursuivie, ce qui conduit à compléter la base de fonctions de la variable d'espace.

Cette technicité spécifique aux calculs multirésolutions est ici détaillée ; elle s'appuie sur les principes de la méthode LATIN dont seules les grandes lignes sont rappelées. Des exemples d'application de cette stratégie dans le cadre, d'une part, du flambage élastique de poutres et d'arcs plans et, d'autre part, de l'élasto-visco-plasticité sont présentés. Enfin, on donne les premières conclusions relatives à cette étude.
## 2. Rappels sur la méthode LATIN

### 2.1. Notations

La structure étudiée occupe, dans la configuration initiale, un domaine de l'espace euclidien à trois dimensions, dont l'espace vectoriel associé est noté $E^3$. Les vecteurs de $E^3$, le déplacement par exemple, sont notés $\mathbf{U}$. Un point courant de $\Omega$ est noté $\mathbf{M}$. Les opérateurs linéaires de $E^3$, déformation et contrainte, sont respectivement notés $\varepsilon$ et $\sigma$. Ces quantités sont représentées dans une base par une matrice $(3,3)$ ; cette notation est aussi utilisée pour les variables internes. Les opérateurs plus complexes, comme, par exemple, le tenseur d'élasticité de Hooke, sont notés $\mathbf{K}$. Pour définir la puissance, on utilise l'opérateur trace noté « Tr », tel que :

$$ \mathrm{Tr}[\sigma \varepsilon] = \sum_{i=1,2,3} [\sigma]_{ij} [\varepsilon]_{ij} $$

où $[\sigma]_{ij}$ et $[\varepsilon]_{ij}$ sont respectivement les composantes des opérateurs $\sigma$ et $\varepsilon$.

### 2.2. Problème de référence

Le milieu étudié occupe dans sa configuration initiale le domaine $\Omega$, de frontière $\Omega$. Il est soumis à chaque instant $t \in [0, T]$ à des forces volumiques $f_d$ dans $\Omega$, à des forces surfaciques $F_d$ sur la partie $2\Omega$ de la frontière $\Omega$ et à des déplacements $\mathbf{U}_d$ sur $1\Omega$, partie complémentaire à $2\Omega$ de la frontière $\Omega$. La méthode est développée en petites perturbations pour simplifier la présentation des idées. Elle est donnée dans [LAD 99] pour les grandes transformations et est développée dans [BOU 96], [BOU 97] pour l'analyse du flambage et du post-flambage.

On se place donc dans le cadre des petites perturbations avec un chargement quasi-statique. $\mathcal{A}[0,T]$, $\mathcal{S}[0,T]$ sont les espaces où le déplacement $\mathbf{U}$ et la contrainte $\sigma$, définis sur $\Omega \times [0,T]$, sont recherchés. Les sous-espaces des déplacements cinématiquement admissibles à zéro sont définis par :

$$ \mathcal{A}_0 = \{ \mathbf{U}(\mathbf{M}) ; \mathbf{U}(\mathbf{M}) \mid 1\Omega = 0 \} $$

$$ \mathcal{A}_0^{[0,T]} = \{ \mathbf{U}(\underline{\mathbf{M}},t) \in \mathcal{A}^{[0,T]} ; \forall t \in [0,T] \quad \underline{\mathbf{U}}(\underline{\mathbf{M}},t) \in \mathcal{A}_0 \} $$

L'état du matériau est complètement défini par la donnée des valeurs de la déformation $\varepsilon$, de la déformation inélastique $\varepsilon^P$ et d'autres variables internes notées $\underline{\mathbf{X}}$. En général, $\underline{\mathbf{X}}$ est un élément de $\mathbf{R}^m$. La quantité conjuguée de $\underline{\mathbf{X}}$ est notée $\underline{\mathbf{Y}}$ ; c'est aussi un élément de $\mathbf{R}^m$.

On peut alors définir le problème de référence pour la structure $\Omega$ :

Trouver $\mathbf{U}(\underline{\mathbf{M}},t)$ et $\sigma(\underline{\mathbf{M}},t)$ $\underline{\mathbf{M}} \in \Omega, t \in [0,T]$ tels que :
- équations de liaison :

$$
\underline{U}(\underline{M},t) \in \mathcal{U}^{[0,T]}, \quad \underline{U} \mid_{1\Omega} = \underline{U}_d
$$
$$
\varepsilon(\underline{M},t) = \frac{1}{2} \{ \mathrm{Grad} \, \underline{U} + \mathrm{Grad}^T \, \underline{U} \}
$$

- équation d'équilibre :

$$
\sigma \in \mathcal{S}^{[0,T]}, \forall \underline{U}^* \in \mathcal{U}^{[0,T]} \\
\int_{\Omega} \mathrm{Tr}[\sigma \underline{\varepsilon}(\underline{U}^*)] d\Omega = \int_{\Omega} \underline{f}_d \cdot \underline{U}^* d\Omega + \int_{2\Omega} \underline{F}_d \cdot \underline{U}^* dS
$$

- lois d'état

$$
\dot{\underline{\varepsilon}} = \dot{\underline{\varepsilon}}^P + \dot{\underline{\varepsilon}}^e \qquad \underline{\sigma} = \underline{K} \underline{\varepsilon}^e \qquad \underline{Y} = \underline{G}(\underline{X})
$$

- lois d'évolution

$$
\frac{d}{dt} \begin{bmatrix} \underline{\varepsilon}^P \\ \underline{\varepsilon} \end{bmatrix}_{\underline{X}} = \mathbf{B} \begin{bmatrix} \underline{\sigma} \\ \underline{Y} \end{bmatrix} \qquad \underline{\varepsilon}^P(t=0) = 0 \text{ et } \underline{Y}(t=0) = 0
$$

K, G sont des opérateurs caractéristiques du matériau. B est un opérateur positif
qui dépend du matériau.

2.3. Principes de la méthode

On introduit $S = (\dot{\underline{\varepsilon}}^P, \dot{\underline{X}}$ ; $\underline{\sigma}, \underline{Y}$) défini sur $\Omega x[0, T]$. L'espace associé est noté
$\mathbf{S}^{[0,T]}$. $S$ est l'élément recherché et est défini par le groupe d'équations précédent.
La méthode LATIN consiste à construire une suite d'éléments $\mathbf{s}_n$ définis sur
$\Omega x[0,T]$ et est basée sur 3 principes.

- **P1** : les équations du problème sont séparées en deux groupes définissant deux ensembles $\Gamma$ et $\mathbf{A}_d$

- $\mathbf{A}_d$ : équations linéaires éventuellement globales,

- Γ : équations locales en variable d'espace éventuellement non linéaires

- P2 : La méthode suit un schéma itératif à deux étapes.

Une itération est constituée d'une étape locale où l'on construit une solution appartenant à $\Gamma$ et d'une étape globale où l'on cherche une solution appartenant à $\mathbf{A}_d$. La solution du problème posé est l'intersection de ces deux ensembles. Pour résoudre ces deux étapes, deux directions de recherche $E^+$ et $E^-$ sont utilisées, qui sont les paramètres de la méthode.

- **P3** : Les inconnues du problème sont représentées par une somme de produits de fonctions du temps et de l'espace.
Pour une formulation normale du matériau (lois d'état linéaires, $\bar{Y} = \Lambda \bar{X}$, ce qui est quasiment toujours possible par changement de variables [LAD 99]), les seules relations non linéaires sont les lois d'évolution, elles sont de plus locales en variable d'espace.

Il existe donc un découpage naturel définissant les deux ensembles $A_d$ et $\Gamma$ à partir des variables : $s = (\epsilon P, \bar{X}; \sigma, \bar{Y})$ définies sur $x [0,T]$. Le découpage utilisé est tel que :

équations de liaison
équation d'équilibre
lois d'état, $\sigma = K (\epsilon - \epsilon P)$, $\bar{Y} = \Lambda \bar{X}$

$$ \Gamma \quad \left| \quad \text{lois d'évolutions} \frac{d}{dt} \left[ \epsilon P \right]_{-\bar{X}} = B \left[ \sigma \right]_{\bar{Y}} \quad \epsilon P(t=0) = 0 \text{ et } \bar{X}(t=0) = 0 $$

Le premier groupe d'équation est associé à l'énergie libre, le second à la dissipation. $A_d$ est le sous-espace linéaire des éléments $s$ vérifiant le premier groupe d'équations. Les éléments $s$ vérifiant le second groupe d'équations constituent une variété non linéaire notée $\Gamma$. Le problème est donc de trouver l'intersection de $A_d$ et $\Gamma$.

Le second point de la méthode consiste à définir un algorithme itératif à deux étapes (figure 1). Le principe de la méthode consiste à partir d'un élément initial $s_0$ appartenant à $A_d$ (obtenu par un calcul élastique), à construire des éléments appartenant successivement à $\Gamma$ et à $A_d$ jusqu'à converger vers la solution $s_{ex}$, intersection de ces deux ensembles. Une itération, c'est-à-dire, la construction d'un nouvel élément $s_{n+1}$ de l'ensemble $A_d$, connaissant $s_n$ (élément de $A_d$), comporte deux étapes : l'étape locale et l'étape linéaire globale.

<figure><img src="image_2.png" /><figcaption>Figure 1. Schéma itératif</figcaption></figure>
L'étape locale correspond à la recherche d'un élément $\hat{s}$ appartenant à $\Gamma$, il est déterminé à partir d'un élément $s_n$ de $A_d$ et d'une direction de recherche $E^+$, soit :

$$ \hat{s} = \Gamma \cap (s_n + E^+) $$

$E^+$ est une variété linéaire à caractère local en espace de forme connue de telle sorte que le problème définissant $\hat{s}$ soit local. $E^+$ peut a priori dépendre implicitement de $\hat{s}$. En pratique, à un instant t de l'intervalle d'étude $[0, T]$ le problème à résoudre est local en variable d'espace mais non linéaire.

L'étape linéaire globale correspond à la recherche d'un élément $s_{n+1}$ appartenant à $A_d$. Cet élément est déterminé à partir de $\hat{s}$ (élément de $\Gamma$) et d'une direction de recherche $E^-$, soit :

$$ s_{n+1} = A_d \cap (\hat{s} + E^-) $$

$E^-$ est une variété linéaire qui a priori peut dépendre de $\hat{s}$ calculé précédemment. Ainsi la construction de $s_{n+1}$ est un problème linéaire, sur $x \in [0,T]$.

En pratique, les variétés $E^+$ et $E^-$ sont définis par deux opérateurs $H^+$ et $H^-$ tels que :

$$ E^+ = \{ \Delta s / dt \begin{bmatrix} \Delta \epsilon P \\ -\Delta \underline{X} \end{bmatrix} + H^+ \begin{bmatrix} \Delta \sigma \\ \Delta \underline{Y} \end{bmatrix} = 0 ; \Delta s \in S^{[0,T]} \} $$

$$ E^- = \{ \Delta s / dt \begin{bmatrix} \Delta \epsilon P \\ -\Delta \underline{X} \end{bmatrix} - H^- \begin{bmatrix} \Delta \sigma \\ \Delta \underline{Y} \end{bmatrix} = 0 ; \Delta s \in S^{[0,T]} \} $$

$H^+$ et $H^-$ sont deux opérateurs qui peuvent dépendre de $\hat{s}$ : $H^+$ est un opérateur local en espace, éventuellement non linéaire afin de conserver le caractère local en espace de l'étape locale. $H^-$ est un opérateur linéaire, éventuellement global, pour conserver le caractère linéaire de l'étape globale.

Dans les différents exemples traités par la suite, les directions de recherche sont choisies de la manière suivante :

$H^+$: direction « verticale »

$$ E^+ = \{ \Delta s / \begin{bmatrix} \Delta \sigma \\ \Delta \underline{Y} \end{bmatrix} = 0 ; \Delta s \in S^{[0,T]} \} $$

$H^-$: direction tangente à $\Gamma$ au point $\hat{s}$.

En élasto-visco-plasticité par exemple, l'opérateur B est différentiable, ainsi l'opérateur tangent L à l'ensemble $(\Gamma)$ au point $\hat{s}$ peut être défini et :

$$ H^- = L $$
L'étape locale conduit à un problème d'évolution local en variable d'espace et de temps, en général non linéaire. La résolution de cette étape ne pose pas de difficulté particulière. Par contre, l'étape globale nécessite la résolution d'un problème linéaire global sur l'espace temps pour lequel des techniques de résolution spécifiques sont nécessaires afin de limiter son coût.

## 2.4. Traitement de l'étape globale – Représentation des inconnues

L'étape globale consiste à rechercher un nouvel élément $s_{n+1}$ meilleur que $s_n$ ; la différence entre ces deux éléments de $A_d$ peut être choisie comme l'inconnue de cette étape. En travaillant sur les corrections, deux nouvelles quantités s'introduisent naturellement :

$$ 
\begin{aligned}
\Delta \hat{S} &= \hat{S} - s_n = (\Delta \hat{\epsilon} P, \Delta \hat{X}; \Delta \hat{\sigma}, \Delta \hat{Y}), \\
\Delta s_n &= s_{n+1} - s_n = (\Delta \hat{\epsilon}_n^P, \Delta \hat{X}_n, \Delta \hat{\sigma}_n, \Delta \hat{Y}_n)
\end{aligned}
$$

Il est important de noter que l'élément $s_{n+1}$ doit appartenir à $A_d$ pour rester dans le cadre de la méthode LATIN.

La résolution directe des problèmes globaux linéaires paramétrés par le temps, intervenant dans l'étape globale peut conduire à des temps de calculs importants. Pour abaisser ces temps, on introduit des représentations des inconnues en temps et en espace pour résoudre à faible coût numérique l'étape globale. La séparation de la dépendance fonctionnelle dans l'espace et le temps à donné de bons résultats pour les chargements quasi-statiques [BOU 96], [LAD 91], [LAD 99]. Les corrections à chaque itération sont recherchées comme superposition de solutions de type chargement simple. Ces solutions sont connues pour être de bonnes approximations des problèmes non linéaires en quasi-statique.

Pour chaque variable, la représentation des incréments de variables $\Delta \hat{\epsilon}_n^P, \Delta \hat{X}_n, \Delta \hat{\sigma}_n, \Delta \hat{Y}_n$ est :

$$ 
\begin{aligned}
\Delta \hat{\epsilon}_n^P (\underline{M}, t) &= \dot{\hat{\epsilon}}_{n+1}^P (\underline{M}, t) - \hat{\epsilon}_n^P (\underline{M}, t) = \sum_{i=1}^{m} a_i(t) A_i (\underline{M}) \\
\Delta \hat{X}_n (\underline{M}, t) &= \dot{\hat{X}}_{n+1} (\underline{M}, t) - \hat{X}_n (\underline{M}, t) = \sum_{i=1}^{m} c_i(t) C_i (\underline{M}) \\
\Delta \hat{\sigma}_n (\underline{M}, t) &= \hat{\sigma}_{n+1} (\underline{M}, t) - \hat{\sigma}_n (\underline{M}, t) = \sum_{i=1}^{m} b_i(t) B_i (\underline{M}) \\
\Delta \hat{Y}_n (\underline{M}, t) &= \dot{\hat{Y}}_{n+1} (\underline{M}, t) - \hat{Y}_n (\underline{M}, t) = \sum_{i=1}^{m} d_i(t) D_i (\underline{M})
\end{aligned}
$$

En pratique, m est limité à 1 ou 2. On ajoute donc, à chaque itération de la méthode un ou deux couples fonction du temps, fonction d'espace. Ainsi, après n itérations, et en supposant qu'on ajoute un couple à chaque itération (m=1), chaque
champ est représenté par n couples de fonctions auquel vient s'ajouter la solution de l'initialisation. La contrainte, par exemple, s'écrit :

$$ \sigma_n(\underline{M}, t) = b_0(t) \mathbb{E}_0(\underline{M}) + \sum_{i=1}^{n} b_i(t) \mathbb{E}_i(\underline{M}) $$

où le produit $b_0(t) \mathbb{E}_0(\underline{M})$ est issu de l'initialisation : typiquement $\mathbb{E}_0(\underline{M})$ est obtenu à partir d'une résolution élastique et $b_0(t)$ dépend du chargement.

Les fonctions du temps sont des fonctions scalaires entièrement déterminées par l'algorithmme ; les représentations de type éléments finis sont utilisées pour la discrétisation des fonctions de l'espace. Chaque fonction d'espace générée contient alors les inconnues nodales associées à la discrétisation éléments finis. Cette représentation est intéressante pour analyser les résultats des simulations non linéaires, en particulier si la solution en temps et en espace peut être représentée avec un faible nombre de couples de fonctions. Ainsi il est important de rechercher les meilleurs couples de fonctions pour représenter les corrections dans le but de réduire leur nombre. Il est à noter que pour les problèmes de flambage, la démarche suivie dans [COC 93] bien que différente va dans le même sens ; peu de fonctions de bases sont nécessaires si elles sont bien choisies.

## 2.5. Étape préliminaire

C'est une étape que l'on effectue dans la résolution du problème global, où l'on cherche à construire une approximation en exploitant les fonctions de la variable d'espace calculées lors des itérations précédentes. Cette phase préliminaire fait partie des techniques d'amélioration de la convergence qui consistent à exploiter la base réduite des fonctions qui se construit au fur et à mesure des itérations. Un exemple typique est la technique BFGS associée à la méthode de Newton-Raphson [GER 81].

L'étape préliminaire est insérée dans l'algorithmme comme indiquée figure 2. Elle consiste à utiliser au mieux la base de fonctions d'espace déjà existante pour représenter les corrections. A partir de la base de fonctions d'espace existante, on recherche les corrections sur les fonctions du temps associées à l'étape globale. Un résidu est calculé à partir de la vérification de la direction de recherche $E^-$. Il permet d'estimer si la base de fonctions d'espace existante est suffisante pour résoudre l'étape globale ou si, au contraire, une nouvelle fonction d'espace est nécessaire à cette itération.
<figure><img src="image_3.png" /><figcaption>Figure 2. Présentation de l’algorithme</figcaption></figure>

En pratique, on a donc deux types de situations :

— l'étape préliminaire ne permet pas d'améliorer suffisamment l'approximation de la solution, et il faut résoudre l'étape globale pour obtenir de nouvelles fonctions d'espace. Cette dernière étape nécessite un calcul coûteux car il exige la résolution de l'étape linéaire globale de la méthode LATIN.

— l'étape préliminaire permet, par la détermination de corrections des fonctions du temps existantes, de construire une approximation suffisante : c'est un calcul peu coûteux car il conduit à résoudre un système d'équations différentielles dont la taille est le nombre de fonctions d'espace existantes et est donc très petite. Les coefficients de cette équation sont les intégrales sur le domaine de quantités connues.
### 3. Méthode multirésolution

La méthode LATIN permet d'obtenir une représentation de la solution du problème sous la forme d'une somme de produits de fonctions du temps et d'espace. Ainsi, on peut considérer que la méthode construit automatiquement une base optimale pour représenter la solution. L'idée est donc d'utiliser cette base de représentation pour déterminer la solution d'un problème similaire à celui pour lequel cette base a été construite.

Par similaire, on entend un problème pour lequel, dans un premier temps, on considère que la géométrie et les conditions aux limites restent inchangées (ceci pour respecter la condition d'admissibilité à « 0 » des fonctions d'espace) mais pour lequel les caractéristiques matériaux peuvent varier sur la structure.

On exploite ici le fait que l'algorithme peut être initialisé par une solution quelconque appartenant à l'ensemble Ad (usuellement une solution élastique). L'idée est donc d'initialiser le processus associé au problème similaire (structure avec « défaut ») par le résultat du calcul sur la structure sans « défaut ». On dispose ainsi dès la première itération du calcul sur la structure avec défaut, d'une base de fonctions d'espace à fort contenu mécanique permettant de représenter la solution.

L'étape préliminaire joue ici un rôle fondamental : elle permet d'une part, de vérifier que la base de fonctions d'espace est bien adaptée au problème à résoudre et, d'autre part, de rechercher de nouvelles fonctions du temps permettant d'obtenir la solution du problème avec défaut. Dans le meilleur des cas, si la base est suffisante, on ne générera aucune nouvelle fonction d'espace, et on obtiendra ainsi la solution du problème à faible coût. L'alternative est le calcul automatique de nouvelles fonctions d'espaces pour enrichir la base initiale : le nombre de ces fonctions doit rester faible pour que la méthode soit efficace.

Si les solutions des problèmes avec et sans défaut sont suffisamment proches, on peut alors générer la solution du deuxième problème avec un coût beaucoup plus faible qu'un calcul complet.

### 4. Exemples

La méthode multirésolution est testée sur deux cas d'applications :

- le calcul du flambage de poutres et d'arcs élastiques plans

- le calcul de structures élasto-visco-plastiques.

Le second cas utilise la méthode LATIN dans la situation où elle est rappelée ici. Par contre, le calcul du flambage de poutres et d'arcs élastiques plans est mené avec l'extension de la LATIN proposée dans [LAD 99] et développée dans [BOU 96], [BOU 97].
Pour tous les exemples traités, le critère d’arrêt des itérations de la méthode LATIN est basé sur le fait qu’à chaque itération, après l’étape locale, on a :

- $s_n$ qui vérifie les équations d’état et d’équilibre,

- $\hat{S}$ qui vérifie les lois d’évolution

Si $s_n$ et $\hat{S}$ sont confondus, ils le sont nécessairement avec la solution exacte. La distance utilisée pour caractériser la qualité de la solution est une norme en énergie libre :

$$ \eta^2 = \frac{\| \hat{S} - s_n \|^2}{\frac{1}{2} ( \| \hat{S} \|^2 + \| s_n \|^2 )} $$

où le calcul de la norme $\| \hat{S} \|$ est faite à partir de l'énergie libre associée à chacune des quantités $\hat{\epsilon} P, \hat{\dot{X}}, \hat{\sigma}, \hat{Y}$ [LAD 99]. En pratique, on considère que la solution du problème est obtenue lorsque cette indicateur est inférieur à 5 %.

<figure><img src="image_4.png" /><figcaption>Figure 3. Influence du défaut sur la charge critique</figcaption></figure>
<figure><img src="image_5.png" /><figcaption>**Figure 4.** Nombre de fonctions d'espace supplémentaires générées (base initiale composée de 6 fonctions)</figcaption></figure>

<figure><img src="image_6.png" /><figcaption>**Figure 5.** Influence du défaut</figcaption></figure>

Le premier exemple traité est celui d'une poutre droite encastrée à une extrémité et soumise à un déplacement imposé à l'autre. Les défauts introduits sur la structure sont des variations de - 60 % à +100 % par pas de 10 % du module d'Young sur différents éléments. On étudie ainsi l'influence d'un défaut sur la valeur de la charge
critique de flambement (rapport F_c/F_c0). Les résultats sont présentés sur les figures 3 et 4.

<figure><img src="image_7.png" /><figcaption>**Figure 6.** Nombre de fonctions d'espace supplémentaires générées (base initiale composée de 6 fonctions)</figcaption></figure>

Un calcul similaire a été effectué sur une barre bi-encastrée. Les résultats sont présentés figures 5 et 6. Par souci de clarté, on ne représente sur ces différentes figures que quelques valeurs de défaut : de -60 % à +80 % par pas de 20 % sur les figures 3 et 4, de -40 % à +100% par pas de 20 % sur les figures 5 et 6.

Pour les deux cas traités, les figures 4 et 6 présentent le nombre de fonctions d'espace ajoutées lors du calcul à la base initiale (6 couples de fonctions temps-espace). C'est un indicateur de coût du calcul car cette phase est l'étape la plus coûteuse de l'algorithme. On constate que pour la majorité des cas traités, on ajoute au plus une fonction d'espace. La base initiale de fonctions permet donc de mener de nombreux calculs avec des défauts à un coût beaucoup plus faible qu'un calcul complet : dans les cas présentés, le temps de calcul nécessaire à l'obtention de la solution sur la barre avec défaut est de 10 à 20 % du temps associé à un calcul complet.

L'exemple suivant concerne le flambage d'un arc circulaire non symétrique (rayon 1m, angle d'ouverture 70°, module d'Young E = 200 000 Mpa, section circulaire de rayon 5 10⁻² m). On impose un déplacement strictement croissant au point central de l'arc. Le défaut structurel introduit consiste en une variation du module d'Young compris entre -60% et +100% sur les 7 éléments du maillage. La réponse typique (courbe Effort/Déplacement du point milieu de l'arc) est présentée figure 7. L'influence d'un défaut sur la force critique (point de tangence nulle
indiqué figure 7) peut ainsi être étudié. La figure 8 présente l'évolution du rapport force critique avec défaut sur force critique sans défaut ($F_c/F_{c0}$) en fonction des différents paramètres.

<figure><img src="image_8.png" /><figcaption>Figure 7. Réponse typique</figcaption></figure>

Dans ce cas, le nombre de fonctions d'espace ajouté durant les différents calculs ne dépasse jamais 3 (la base initiale est composée de 7 couples de fonctions temps-espace).

Le dernier exemple concerne une structure élasto-visco-plastique axisymétrique (loi de comportement de Chaboche) dont le maillage (éléments triangulaires à 3 nœuds) comporte 416 degrés de liberté, encastrée sur son rayon extérieur, et soumise à un chargement de type « pression-flexion » (incliné de 10° par rapport à l'horizontale) sur son rayon intérieur. L'évolution temporelle du chargement suit un quart de période de sinusoïde, dont l'amplitude maximale de 30 Mpa est atteinte en 10 secondes. Le défaut introduit est une variation de la limite élastique de – 80% à + 150 % par rapport à sa valeur de référence de 10 Mpa (12 valeurs discrètes). La répartition spatiale est faite avec un découpage géométrique de la structure en 20 zones de l'intérieur vers l'extérieur.

Le problème est de calculer une charge qualifiée de « limite » pour laquelle la raideur de la structure est très faible. Les figures 10 et 11 montrent respectivement le nombre de fonctions d'espace ajoutées à la base initiale (base composée de 7 fonctions d'espace), ainsi que la variation du rapport $F_{lim}/F_0$. L'effort $F_0$ est obtenu à
partir d'une courbe de comportement global décrivant l'évolution de la résultante de l'effort sur le rayon intérieur en fonction du déplacement sur ce même bord. La charge limite $F_0$ (resp. $F_{lim}$) est la valeur de la charge pour laquelle la pente de cette courbe de comportement globale pour la structure sans défaut (resp. avec défaut) est égale au dixième de la pente initiale.

<figure><img src="image_9.png" /><figcaption>Figure 8. Influence d'un défaut sur la charge critique</figcaption></figure>

<figure><img src="image_10.png" /><figcaption>Figure 9. Géométrie de la structure</figcaption></figure>
Comme sur les exemples précédents, on peut observer que pour la majorité des cas traités, au plus deux fonctions d'espace ont été ajoutées à la base initiale pour représenter la solution. Le temps de calcul nécessaire pour obtenir les 240 solutions sur les structures avec défauts est de 2 h 30 mn avec la méthode LATIN sur station de travail. Le temps de calcul associé à 240 résolutions complètes est de l'ordre de 26 h.

<figure><img src="image_11.png" /><figcaption>**Figure 10.** Nombre de fonctions d'espace supplémentaires générées (base initiale composée de 7 fonctions)</figcaption></figure>

**REMARQUE.** — Dans les exemples, on s'est ici intéressé au calcul d'une charge particulière pour une famille de problèmes perturbés. Une autre question, plus large, qui peut être conduite par la même démarche, est de calculer la réponse complète sur un intervalle de temps donné pour une famille de problèmes perturbés. Évidemment, tous les calculs menés le sont avec une précision qui est donnée à l'avance.
<figure><img src="image_12.png" /><figcaption>Figure 11. Influence du défaut sur la réponse ($F_{tim}/F_0$)</figcaption></figure>

## 5. Conclusion

Les premiers exemples numériques montrent un très bon comportement de l'algorithmes dans le cas de résolutions multiples tant pour l'analyse du flambage que pour l'élasto-visco-plasticité.

La base de fonctions d'espace générée par une première résolution sur la structure parfaite est très bien adaptée aux calculs réalisés sur les structures avec défauts, tant que ceux-ci ne perturbent pas de façon excessive la réponse. D'autre part, l'étape préliminaire permet d'exploiter toute la richesse de cette base, assurant ainsi des coûts de calcul très réduits.

Cette première étude donne les éléments permettant, pour une charge et pour une distribution de défauts données, de déterminer la probabilité de ruine de la structure par flambage [BEN 98], [MAC 97].

Enfin, la démarche est très générale et devrait pouvoir s'appliquer à bien d'autres problèmes non linéaires qui nécessitent de multiples résolutions.
## 6. Bibliographie

[BEN 98] BENAROYA H., REHAK M., « Finite element methods in probabilistic structural analysis: A selective review », *Applied Mechanics Reviews*, Vol. 41, No. 5, 1998, p. 201-213.

[BOU 97] BOUCARD P.A., LADEVEZE P., POSS M., ROUGEEP., « A non-incremental approach for large displacement problems », *Computers and Structures* Vol. 64, 1-4, 1997, p. 499-508.

[COC 93] COCHELIN B., AZRAR L., DAMILN., POTTIER-FERRY M., « An asymptotic-numerical method to compute the post-buckling beahviour of elastic plates and shells », *International Journal for Numerical Methods in Engineering*, Vol. 36, 1993, p. 1251-1277.

[GER 81] GERADIN M., IDELSOHN S., HOGGE M., « Computational strategies for the solution of large non linear problems via quasi Newton methods», *Computers and Structures*, Vol. 13, 1981, p. 73-81.

[LAD 91] LADEVEZE P., « New advances in the LArge Time INcrement method », *New Advances in Computational Structural Mechanics*, Ladevèze P., Zienkiewicz O.C. Eds., 1991, p. 3-21.

[LAD 99] LADEVEZE P., *Non Linear Mechanics of Structures - New approaches and non-incremental methods of calculation*, Editions Springer, 1999, existe en français, Editions Hermès, 1996.

[MAC 97] MACIAS O.F. - LEMAIRE M., « *Eléments Finis stochastiques et Fiabilité Application en mécanique de la rupture* », Revue française de génie civil, Vol. 1 - No. 2, 1997.

[SHI 72] SHINBOZUKA M., ASTILL J., « Random eigenvalue problems in structural mechanics », *AIAA Journal*, Vol. 10, No. 4, 1972, p. 456-462.