# Sur une méthode directe de détermination de réponse transitoire: application à un problème thermomécanique

François Comte, Pierre Burry, Habibou Maitournam, T Mac Lan

Nguyen-Tajan

► To cite this version:

François Comte, Pierre Burry, Habibou Maitournam, T Mac Lan Nguyen-Tajan. Sur une méthode directe de détermination de réponse transitoire: application à un problème thermomécanique. 7e colloque national en calcul des structures, CSMA, May 2005, Giens, France. hal-01812983

HAL Id: hal-01812983

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

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>
# Sur une méthode directe de détermination de réponse transitoire : application à un problème thermomécanique

François Comte*, **, *** — Pierre Burry** — Habibou Maitournam*
— T. Mac Lan Nguyen-Tajan**

* Laboratoire de Mécanique des Solides, Ecole Polytechnique,
91128 Palaiseau

{comte,habibou}@lms.polytechnique.fr

** PSA Peugeot Citroën, Direction de la Recherche et de l'Innovation,
route de Gisy, 78943 Vélizy Villacoublay
{francois.comte,pierre.burry,maclan.nguyen}@mpsa.com

*** UME/ENSTA, Chemin de la Hunière, 91761 Palaiseau

**RÉSUMÉ.** Cet article présente une méthode numérique alternative de résolution par éléments finis de problèmes transitoires. Comme dans le cas de la méthode à Grand Incrément de Temps la solution est recherchée sur la totalité du trajet de chargement en résolvant alternativement un problème dit global, où l'on recherche une solution vérifiant l'équilibre de la structure et respectant les conditions limites en efforts et en déplacements, et un problème local, où on cherche à vérifier la loi de comportement. L'étape globale étant la plus coûteuse en temps de calcul sa résolution se fait en décomposant le chargement sur une base d'ondelettes et en résolvant l'équilibre pour les coefficients d'ondelette les plus importants. La méthode est illustrée sur un exemple de plaque trouée.

**ABSTRACT.** An alternative numerical method to solve time dependant problem is presented. As in the Large Time Increment Method the problem is solved on the whole interval separating a global stage where the solution is searched SA and KA and a local stage where the constitutive law has to be verified. To solve the global stage, consuming the most CPU time, the loading is resolved into a wavelet basis and, for the most important coefficients, the equilibrium is solved.

**MOTS-CLÉS :** solides et structures, méthode numérique, analyse transitoire, ondelettes.

**KEYWORDS:** solids and structures, numerical method, time dependant problem, wavelets.
# 1. Introduction

La simulation numérique est devenue aujourd'hui incontournable dans le processus de conception et de validation des composants automobiles. Or le grand nombre de ddl de ces structures (plus d'un million dans le cas d'une culasse), les chargements complexes appliqués et l'emploi de lois de comportement non linéaires rendent les temps de calcul prohibitifs. Le calcul ne peut ainsi pas être utilisé pour rechercher la solution optimale.

L'objectif du travail présenté est de proposer une nouvelle méthode de simulation permettant de diminuer les durées de simulation, méthode fondée sur l'utilisation de la méthode à Grand Incrément de Temps et la décomposition du chargement sur une base d'ondelettes.

# 2. Formulation du problème et principe de l'algorithme

On considère une structure Ω soumise, sur l'intervalle de temps [0,T], à un chargement volumique $f(x,t)$, à des efforts surfaciques $T^d$ sur $\partial\Omega_t$ et à des déplacement imposés $u(x,t)$ sur la surface $\partial\Omega_u$, complémentaire à $\partial\Omega_t$. Un champ de température $T(x,t)$, peut également lui être appliqué. L'évolution est supposée quasi-statique et vérifie les hypothèses de petites perturbations. Le comportement du matériau est thermoélastoplastique avec, par souci de simplicité, $\varepsilon^p$ comme seule variable interne. La résolution du problème doit permettre d'obtenir sur $[0,T]\times\Omega$ l'évolution de $\{u, \underline{\underline{\varepsilon}}^p\}$. $\underline{\underline{\varepsilon}}$ est alors un champ statiquement admissible (SA), $u$ cinématiquement admissible (CA) et la loi de comportement relie les grandeurs $\underline{\underline{\varepsilon}}$, $\underline{\underline{\varepsilon}}^p$ et $\sigma$.

De manière générale l'approche développée repose sur la méthode à Grand Incrément de Temps [LAD 96] proposée par P. Ladevèze et la méthode cyclique directe [MAI 02] développée par H. Maitournam et B. Pommier. De la première méthode on retient l'idée de traiter le problème sur la totalité de l'intervalle de temps. Les équations y sont séparées en équations globales (on parle de problème global) linéaires et équations locales (on parle de problème local) éventuellement non linéaires. Les deux types d'équations sont résolus alternativement jusqu'à convergence pour tous les instants de chargement. De la seconde méthode, utilisée pour déterminer directement l'état asymptotique d'une structure sous chargement cyclique, on retient l'idée de projection sur une base adaptée des champs durant l'étape globale, étape la plus coûteuse en temps de calcul, afin de réduire, à chaque itération, sa durée de simulation. Ces principes sont ici utilisés pour la recherche, non pas d'un état stabilisé sous chargement cyclique, mais de la solution d'un problème transitoire sous chargement variable sur un intervalle $[0,T]$.

Le processus itératif, résolu sous forme résiduelle, se décompose de la manière suivante :

– étape globale : une solution à la fois SA et CA est recherchée sur tout l'intervalle $[0,T]$. Pour cela la forme résiduelle classique $K\delta u(t) = R(t)$ est projetée sur une base
de fonctions. La base est choisie pour approcher au mieux le résidu avec un nombre fixe N de vecteurs. A cette étape, l'utilisation d'une base de fonctions à support compact permet de localiser la projection du résidu là où l'écart avec la solution obtenue à l'itération précédente est le plus important.

– étape locale : à partir de la solution de l'étape globale on détermine une solution plastiquement admissible (PA) sur tout le trajet. Celle-ci est obtenue de manière incrémentale, aux différents instants, par projection sur le convexe d'élasticité en utilisant l'algorithme de retour radial [NGU 77] [SIM 85] [SIM 98].

– actualisation du résidu R(t) : une fois la solution σ<sup>PA</sup>(t) déterminée le résidu est recalculé.

A chaque itération le nouveau résidu vaudra

$$
R(t) = \int_{\Omega} N^T f(t) dv + \int_{\partial \Omega_t} N^T T^d(t) ds - \int_{\Omega} B^T \sigma^{PA} dv. \quad [1]
$$

L'initialisation de l'algorithme est réalisée en résolvant le problème purement élastique. Le critère d'arrêt porte lui sur la condition d'admissibilité plastique de la solution SA et CA déterminée lors de l'étape globale.

Dans la suite les ondelettes ont été choisies comme base de décomposition lors de
l'étape globale.

3. Décomposition d'un signal en ondelettes

L'analyse en ondelettes [MAL 00] d'un signal $s$ est basée sur sa décomposition
sur un «sous-espace d'approximation »et des «sous-espaces de détails »à différentes
échelles. De manière itérative, à partir d'un sous-espace $V_j$, on construit l'approxima-
tion $V_{j+1}$ et les détails $W_{j+1}$. L'ondelette mère $\psi$, de moyenne nulle, est utilisée pour
définir les détails et $\varphi$ sert à la définition des approximations.

$$
s(t) = \sum_{i=j_0}^{j-1} \sum_k d_k^i \psi_{i,k}(t) + \sum_k a_k^{j_0} \varphi_{j_0,k}(t) \quad [2]
$$

Dans de nombreux cas $\psi$ et $\varphi$ sont réelles et présentent des supports compacts. Aux
différentes échelles d'analyse possibles on utilise les fonctions $\psi_{j,n}(t) = 2^{-\frac{j}{2}}\psi(2^{-j}t-n)$ et $\varphi_{j,n}(t) = 2^{-\frac{j}{2}}\varphi(2^{-j}t-n)$ avec $j, n \in \mathbb{Z}$ pour calculer les coefficients d'approximation et de détails :

$$
\left\{
\begin{array}{l}
a_j^n = \langle s, \varphi_{j,n} \rangle = 2^{-\frac{j}{2}} \int_{-\infty}^{+\infty} s(t) \bar{\varphi}(2^{-j}t - n) \\
d_j^n = \langle s, \psi_{j,n} \rangle = 2^{-\frac{j}{2}} \int_{-\infty}^{+\infty} s(t) \bar{\psi}(2^{-j}t - n)
\end{array}
\right.
\qquad [3]
$$

$2^{-j}$ est le facteur d'échelle de l'ondelette et n le paramètre de translation (et donc de localisation). $d_j^n$ et $a_j^n$ traduisent donc ce qui se passe autour de l'instant $2^{-j}n$. Un signal décomposé sur une base d'ondelettes peut ainsi être représenté de manière exacte
par son approximation à une échelle grossière complétée par les détails aux différentes échelles intermédiaires.

Les coefficients d'ondelette des échelles les plus basses décrivent le signal dans sa globalité. A l'inverse les échelles les plus fines servent à décrire ce qui se passe de manière locale. Ainsi aux échelles les plus fines les coefficients d'ondelette seront de faibles valeurs lorsqu'aucune variation rapide du signal n'est observée sur le support de l'ondelette et les seuls coefficients significatifs seront ceux correspondant à des zones proches des points singuliers.

Pour un nombre de coefficients N fixé, l'erreur d'approximation sera minimisée en considérant uniquement les N fonctions de base $\psi_{j,n}$ et $\varphi_{j,n}$ pour lesquels les coefficients $a_j^n$ ou $d_j^n$ sont les plus importants.

De la forme exacte développée (2) on passe à l'approximation

$$s(t) \approx \sum_{i=j_0}^{j-1} \sum_{k' \in [1, m_i]} d_{k'}^i \psi_{ik'}(t) + \sum_{k' \in [1, n]} a_{k'}^{j_0} \varphi_{j_0 k'}(t) \text{ avec } N = n + \sum_i m_i \ll N_{\text{instants}} \quad [4]$$

Dans la suite on notera, sans distinction des fonctions d'approximation ou de détails,

$$s(t) \approx \sum_{i=1}^{N} \alpha_i q_i(t) \left( \text{où } \sum_{i=1}^{N} \beta_i q_i(t) \right)$$

## 4. Nouvelle résolution de l'étape globale

La décomposition du résidu $R(t)$, a priori complexe, sur une base d'ondelettes permet de représenter son évolution avec un nombre réduit de coefficients. De plus, la taille et la localisation des supports étant variables, la sélection de certaines fonctions de base particulières permet de ne corriger la solution que sur une partie de l'intervalle $[0,T]$. Ces 2 propriétés valident l'utilisation des ondelettes comme base de décomposition.

A partir du système $K(t)\delta u(t) = R^{eq}(t)$ on obtient, en décomposant $\delta u$ et $R^{eq}$ sur la base d'ondelettes,

$$\sum_{j=1}^{N} K(t) \alpha_j(M) q_j(t) = \sum_{j=1}^{N} \beta_j(M) q_j(t) \text{ et en multipliant par la fonction test } q_i(t)$$

$$\sum_{j=1}^{N} K(t) \alpha_j(M) q_i(t) q_j(t) = \sum_{j=1}^{N} \beta_j(M) q_i(t) q_j(t). \text{ En intégrant sur } [0,T] \text{ il vient}$$

$$\sum_{j=1}^{N} \left[ \int_{0}^{T} K(t) q_i(t) q_j(t) dt \right] \alpha_j(M) = \sum_{j=1}^{N} \left[ \int_{0}^{T} q_i(t) q_j(t) dt \right] \beta_j(M) \text{ qui se ré-écrit}$$
$$
\tilde{K}\tilde{\alpha}(M) = \tilde{\beta}(M) \qquad [5]
$$

$$
\text{où } \tilde{K}_{ij} = \int_{\text{Supp } q_i(t) \cap \text{Supp } q_j(t)} K(t) q_i(t) q_j(t) dt, \tilde{\alpha} = \begin{bmatrix} \alpha_1(M) \\ \vdots \\ \alpha_N(M) \end{bmatrix} \text{ et } \tilde{\beta} = \begin{bmatrix} \beta_1(M) \\ \vdots \\ \beta_N(M) \end{bmatrix}.
$$

Le système (5) est de taille importante puisque la matrice [K̃] est de taille NxN_ddl
sur NxN_ddl. Cependant [K̃] peut présenter un certain nombre de matrices élémen-
taires [K̃_{ij}] nulles. Ainsi, d'après l'approche multi-échelles de la décomposition en
ondelettes, les supports de q_i(t) et q_j(t) peuvent être disjoints. Dans ce cas, l'inter-
section des supports étant nulle, [K̃_{ij}] sera nulle. Par ailleurs la base {q_k(t)}_{k∈{1,N}}
étant orthonormée, si K(t) est constante sur Supp q_i(t) ∩ Supp q_j(t) alors pour i ≠ j
[K̃_{ij}] sera également nulle. Ces deux cas mènent à des matrices élémentaires nulles.
Par ailleurs [K̃] étant symétrique réelle certaines méthodes de résolutions de systèmes
de grande taille telles MINRES peuvent être utilisées afin d'optimiser la durée de ré-
solution de l'équation (5) [GRE 97].

Autrement, dans un cas de faible plasticité ou si le module tangent est relativement
important par rapport au module d'Young, un modèle de type Newton modifié peut
être utilisé. La matrice de rigidité est alors prise constante et égale à sa valeur initiale.

Dans ce cas particulier $[\tilde{K}_{ij}] = K \delta_{ij}$ et $[\tilde{K}] = \begin{bmatrix} K & & & \\ & K & & \\ & & \ddots & \\ & & & K \end{bmatrix}$. Seuls les termes diagonaux sont conservés et pour chaque coefficient d'ondelettes le système à résoudre devient

$$
K \alpha_i = \beta_i \forall i \in \{1; N\} \qquad [6]
$$

5. Résultats numériques

Afin d'illustrer la méthode l'exemple ci-dessous traite du cas d'une plaque trouée soumise à un chargement thermomécanique avec un comportement élastoplastique à écrouissage cinématique linéaire où la limite d'éléasticité et le module d'écrouissage dépendent de la température (figure 1). De part la symétrie du problème seul un quart de plaque est maillé et des conditions aux limites adéquates sont appliquées. Un premier calcul découplé a permis de déterminer le champ de température T(x, t) dans toute la pièce. Dans le calcul mécanique, simulé sous matlab, la plaque subit simultanément un effort de traction sur la face en x = L/2 dans la direction e_y et un gradient thermique dans l'épaisseur. Le résultat du calcul est présenté sur la figure 2 où sont superposés, au point le plus chargé, les résultats de la méthode incrémentale et de la méthode directe développée ici. Une bonne corrélation entre les deux calculs est ob-
<figure><img src="image_2.png" /><figcaption>**Figure 1.** Cas test de la plaque trouée</figcaption></figure>

<figure><img src="image_3.png" /><figcaption>**Figure 2.** Comparaison des résultats numériques</figcaption></figure>

servée. Par ailleurs, par rapport à une approche incrémentale classique, une diminution des temps de calcul est observée.

## 6. Conclusion

La méthode directe présentée dans cet article permet de déterminer l'état d'une structure soumise à un chargement transitoire. Contrairement à l'approche incrémentale, à chaque itération, la résolution est effectuée sur tout le trajet de chargement. Lors de l'étape globale, le choix de la base de décomposition s'est porté sur les ondelettes de part leur capacité d'adaptation au signal traité. La solution peut alors se représenter avec relativement peu de coefficients même dans des cas complexes. L'exemple traité a permis de montrer l'adéquation des résultats entre méthodes directe et incrémentale.

## 7. Bibliographie

[GRE 97] GREENBAUM A., *Iterative methods for solving linear systems*, SIAM, Philadelphia, 1997.

[LAD 96] LADEVÈZE P., *Mécanique non linéaire des structures*, Hermès, 1996.

[MAI 02] MAITOURNAM H., POMMIER B., THOMAS J., « Détermination de la réponse asymptotique d'une structure anélastique sous chargement thermomécanique cyclique », C. R. Mécanique, vol. 330, 2002, p. 703-708.

[MAL 00] MALLAT S., *Une exploration des signaux en ondelettes*, Les éditions de l'Ecole polytechnique, 2000.

[NGU 77] NGUYEN Q., « On the elastic plastic initial-boundary value problem and its numerical integration », *Int. J. Numer. Methods Engrg.*, vol. 11, 1977, p. 817-832.

[SIM 85] SIMO J., TAYLOR R., « Consistent tangent operators for rate-independent elastoplasticity », *Comput. Methods Appl. Mech. Engrg.*, vol. 48, 1985, p. 101-118.

[SIM 98] SIMO J., HUGUES T., *Computational Inelasticity*, Springer-Verlag, New York, 1998.