# Modélisation du post-flambement de structures minces par la méthode SPH

Jun Lin, Hakim Naceur, Abdelhoudaded Laksimi, Daniel Coutellier

► To cite this version:

Jun Lin, Hakim Naceur, Abdelhoudaded Laksimi, Daniel Coutellier. Modélisation du post-flambement de structures minces par la méthode SPH. 11e colloque national en calcul des structures, CSMA, May 2013, Giens, France. http://hal-01722090

HAL Id: hal-01722090

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

Submitted on 2 Mar 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
# Modélisation du post-flambement de structures minces par la méthode SPH

J. Lin<fnref n="1" />,†, H. Naceur<fnref n="2" />, A. Laksimi<fnref n="1" />, D. Coutellier<fnref n="2" />

**Résumé** — L'objectif de cette communication est de présenter de récents résultats concernant la modélisation du post-flambement de structures minces par la méthode Smoothed Particle Dynamics (SPH) avec une formulation de type coque. En se basant sur la théorie de Reissner-Mindlin pour la prise en compte du cisaillement transversal dans le cas de coques épaisses, la méthode SPH permet de représenter une coque à l'aide d'une seule couche de particules positionnées dans la surface moyenne de la coque. La prise en compte des non linéarités géométriques est faite en adoptant une formulation lagrangienne totale et la résolution du problème non linéaire qui en résulte, est faite par un schéma de type dynamique explicite. La méthode SPH développée a été appliquée avec succès dans la détermination des charges critiques de flambage ainsi que le trajet de chargement non linéaire et les résultats ont été confrontés à ceux obtenus classiquement par la méthode des éléments finis coques en utilisant le logiciel ABAQUS®. Les premiers résultats obtenus montrent que la méthode SPH peut être une alternative à la méthode EF.

**Mots clés** — Méthode SPH, coques épaisses, non linéarités, dynamique explicite,

## 1 Introduction

De nos jours, la méthode des Eléments Finis (EF) est très répondue et largement implémentée dans les logiciels de simulation industriels tels que : ABAQUS, ANSYS, LS-DYNA, etc. Une spécificité marquante de la méthode EF, consiste à fragmenter un solide en un nombre fini d'éléments, ce qui permet d'approximer la solution localement sur un élément et de reconstruire la solution globale par assemblage de ces éléments, ce qui rend la méthode particulièrement efficace. Cependant, elle présente certains inconvénients lors de la résolution des problèmes fortement non linéaires, tels que les procédés de mise en forme ou les problèmes de crash. En effet lorsque les déformations deviennent grandes, les éléments se trouvent alors trop distordus et la qualité des résultats se dégrade très rapidement ou tout simplement le calcul diverge. La nécessité de refaire le maillage devient obligatoire pour éviter la distorsion des éléments, ce qui entraîne forcément des coûts de calculs prohibitifs.

Dans le but de se libérer des problèmes entraînant de grandes distorsions de maillage, ou de la propagation de fissures, les méthodes dites "sans maillage" ont été développées à partir des années 1970. Une des méthodes sans maillage les plus efficaces, est la méthode baptisée Smoothed Particle Dynamics (SPH) développée initialement pour l'astrophysique par Lucy [10], Gingold et Monaghan [5, 12]. Contrairement à la méthode des EF, la méthode SPH est basée sur le principe de collocation, directement en utilisant l'écriture forte des équations d'équilibre. Toutefois, la méthode SPH classique développée pour les écoulements de fluides, emploie des fonctions noyaux eulériens et par conséquent souffre de plusieurs défauts, tels que : la non consistance, l'instabilité dans le cas de chargements de type traction (membrane), ainsi que la présence de modes parasites (où à énergie nulle) [1, 4].

La méthode SPH a été développée et utilisée dans la mécanique de structures, ainsi que dans la modélisation de mise en forme des métaux en utilisant l'approche solide tridimensionnel [8, 9]. Cependant la modélisation de structures minces de type coques en utilisant l'approche solide 3D reste coûteuse, car plusieurs particules doivent être placées dans la direction de la fibre épaisseur ce qui rend la densité de particules très importante. Les formulations de coques par méthode sans maillage, utilisant une seule couche de particules dans la surface moyenne ont suscité l'intérêt de nombreux chercheurs. Par

<footnote n="1">¹ Laboratoire ROBERVAL, UMR 7337 CNRS, Université de Technologie de Compiègne</footnote>

<footnote n="2">² Laboratoire LAMIH, UMR 8201 CNRS, Université de Valenciennes</footnote>

<footnote>† Correspondance: jun.lin@utc.fr, Tél: 03 27 51 14 12</footnote>
exemple, Krysl et Belytschko [7] ont combiné le modèle de plaques de Kirchhoff avec la méthode EFG pour résoudre les problèmes de flexion de plaques minces.

Dans ce travail, nous proposons une extension de la méthode SPH proposée par Chen [3], afin de développer une formulation robuste pour l'analyse non linéaire géométrique de structures minces de type coque. La modélisation est basée sur la théorie de Reissner-Mindlin pour prendre en compte l'influence du cisaillement transversal dans le cas de coques épaisses. La formulation lagrangienne totale (FLT) est également adoptée, qui va permettre de gagner du temps CPU en évitant la recherche des particules les plus proches voisines pour la construction de la fonction noyau à chaque pas de temps. La FLT permettra également d'éliminer certains problèmes d'instabilité [2], [6]. Le problème non linéaire résultant, est résolu en utilisant un schéma d'intégration de type dynamique explicite basé sur l'algorithme de Newmark, permettant ainsi de calculer tout le trajet de chargement au cours du temps.

Plusieurs applications numériques de structures minces flexibles en comportement non linéaire géométrique sont étudiées et les résultats obtenus sont confrontés aux solutions classiques obtenues par les éléments finis de coques par utilisation du code ABAQUS®. Les premiers résultats obtenus montrent que la méthode SPH peut être une alternative à la méthode EF.

Le plan de présentation de cet article est comme suit. Tout d'abord une présentation de la méthode SPH, utilisant la formulation lagrangienne totale sera développée dans la section 2. Dans la section 3, sera présentée la méthode SPH adaptée pour les coques en se basant sur la théorie de Reissner-Mindlin pour la modélisation de structures minces. Différentes applications numériques seront ensuite présentées dans la quatrième section. Enfin, nous terminons par des critiques et des conclusions quant à l'efficacité de la méthode SPH par comparaison à la méthode des EF.

## 2 Méthode SPH en Formulation Lagrangienne Totale

Dans la méthode SPH, le domaine étudié est représenté par un ensemble de nombre fini de points dits "particules SPH". On exprime la composante de champ de déplacement $u(x)$ de la particule $a$ définie par son vecteur position $x_a$ en fonction des paramètres calculés sur les particules $b$ se trouvant dans son voisinage le plus proche de $a$, par l'approximation SPH

$$u(x_a) = \sum_{b=1}^{n} u(x_b) W(x_a - x_b, h) V_b = \sum_{b=1}^{n} u_b \phi_{ab} \quad (1)$$

où $n$ est le nombre de particules $b$ de vecteur position $x_b$ se trouvant dans le voisinage proche de la particule $a$ et $V_b$ est la portion de volume de la particule $b$. $\phi_{ab} = W(x_a - x_b, h)$ $V_b$ est la fonction de forme SPH. $W(x_a - x_b, h)$ est la fonction Gaussienne noyau dépendante de deux paramètres : $h$ qui est la longueur de lissage définissant la taille du domaine d'intérêt et $x_{ab} = \|x_{ab}\| = \|x_a - x_b\|$, la distance entre les particules $a$ et $b$.

Une des fonctions noyau les plus utilisées est la fonction spline cubique normalisée, qui peut être définie en deux dimensions [12], en adoptant $q = x_{ab}/h$

$$W(q, h) = \begin{cases} \frac{10}{7\pi h^2} \left(1 - \frac{3}{2}q^2 + \frac{3}{4}q^3\right) & \text{si } 0 \le q < 1 \\ \frac{5}{14\pi h^2} (2-q)^3 & \text{si } 1 \le q < 2 \\ 0 & \text{si } q \ge 2 \end{cases} \quad (2)$$

Il en découle que le gradient du champ de déplacement au point $x_a$ peut être approché par

$$\nabla u(x_a) = \sum_{b=1}^{n} u(x_b) \nabla W(x_a - x_b, h) V_b = \sum_{b=1}^{n} u_b \nabla \phi_{ab} \quad \text{ou} \quad \nabla u_a = \sum_{b=1}^{n} (u_b - u_a) \nabla \phi_{ab} \quad (3)$$

La dernière équation est généralement préférable car elle peut garantir que la dérivée d'un champ constant peut être estimée correctement, c'est-à-dire que l'approximation est consistante à l'ordre 0. Toutefois, la consistance linéaire n'est pas assurée sur les particules se trouvant sur les frontières du domaine d'étude.
Afin de résoudre ce problème, nous avons adopté la méthode dite "SPH corrigée" (SPHC) [3] en normalisant les dérivées de la fonction de forme. Parmi les méthodes les plus efficaces, nous retenons celle de Randles et Libersky [3] qui consiste à introduire une matrice $\mathbb{B}$ telle que

$$ \nabla u_a = \mathbb{B} \sum_{b=1}^{n} (u_b - u_a) \nabla \phi_{ab} = \sum_{b=1}^{n} (u_b - u_a) \nabla \Phi_{ab} \quad (4) $$

avec $\mathbb{B} = \mathbb{A}^{-1}$, défini par

$$ \mathbb{A} = \sum_{b=1}^{n} (\mathbf{x}_b - \mathbf{x}_a) \otimes \nabla \phi_{ab} \quad (5) $$

En mécanique des solides, l'équation de continuité et l'équation de conservation de la quantité de mouvement pour une particule $a$ ayant pour vecteur position $\mathbf{x}_a$, sont données par

$$ \frac{d\rho_a}{dt} = -\rho_a \operatorname{div}(\mathbf{v}) = -\rho_a \sum_{b=1}^{n} (\mathbf{v}_b - \mathbf{v}_a) \cdot \nabla \Phi_{ab} \quad (6) $$

$$ \frac{d\mathbf{v}_a}{dt} = \operatorname{div}(\boldsymbol{\sigma}) / \rho_a = \sum_{b=1}^{n} \rho_b \left( \frac{\boldsymbol{\sigma}_a}{\rho_a^2} + \frac{\boldsymbol{\sigma}_b}{\rho_b^2} + \Pi_{ab} \mathbb{I} \right) \nabla \Phi_{ab} \quad (7) $$

où $\rho$ est la densité du matériau, $\mathbf{v}_a$ le champ de vitesse de la particule $a$ et $\boldsymbol{\sigma}_a, \boldsymbol{\sigma}_b$ sont respectivement les tenseurs de contraintes de Cauchy des particules $a$ et $b$.

En s'inspirant des travaux de Gingold et Monaghan [5, 12], nous avons introduit un terme de viscosité artificielle $\Pi_{ab}$ permettant à l'algorithme d'être capable de modéliser des problèmes dynamiques avec présence de sollicitations de type chocs/impact. Il existe de nombreux modèles de viscosité de dissipation d'énergie de cisaillement et volumétrique, la plus utilisée est donnée par [12]

$$ \Pi_{ab} = \begin{cases} \frac{\mu}{\bar{\rho}_{ab}} (-\bar{c}_{ab} \vartheta_{ab} + 2 \vartheta_{ab}^2) & \text{si } \mathbf{v}_{ab} \cdot \mathbf{x}_{ab} < 0 \\ 0 & \text{si } \mathbf{v}_{ab} \cdot \mathbf{x}_{ab} \ge 0 \end{cases} \quad (8) $$

$$ \text{et} \quad \vartheta_{ab} = \frac{\mathbf{v}_{ab} \cdot \mathbf{x}_{ab}}{x_{ab}} $$

Ici $\vartheta_{ab}$ représente la vitesse relative des deux particules, projetée sur le vecteur $\mathbf{x}_{ab}$, $\mathbf{v}_{ab} = \mathbf{v}_a - \mathbf{v}_b$, $\bar{\rho}_{ab} = (\rho_a + \rho_b)/2$ est la densité moyenne du matériau et $\bar{c}_{ab} = (c_a + c_b)/2$ la célérité du son dans le matériau. Le coefficient $\mu$ est un facteur scalaire positif représentant l'effet de viscosité dynamique, généralement fixé entre $1 \le \mu \le 4$.

Les équations Eq.(6) et Eq.(7) sont classiquement utilisées en SPH pour les fluides car écrites en formulation eulérienne exprimées donc dans le repère spatial (configuration finale déformée). Pour les solides et les structures, cette écriture n'est pas convenable car elle est coûteuse notamment par la recherche des plus proches voisins à chaque instant d'une part, et d'autre part elle peut conduire au problème d'instabilité en tension (tensile instability). Ce dernier problème a été signalé initialement par Belytschko [2], il peut être atténué par l'utilisation d'une fonction noyau de lagrangienne.

Dans la formulation lagrangienne totale (FLT), la configuration initiale est considérée comme la configuration de référence. Nous utilisons alors, la densité initiale du matériau et le volume initial des particules, de cette façon il devient pas nécessaire de les actualiser. La seule équation intéressante dans ce cas est Eq.(7), dans laquelle la contrainte de Cauchy doit être remplacée par la contrainte de Piola-Kirchhoff de première espèce (PK1). On obtient ainsi pour la particule $a$

$$ \frac{d\mathbf{v}_a}{dt} = \operatorname{div}(\mathbb{P}) / \rho_0 = \sum_{b=1}^{n} \rho_{0b} \left( \frac{\mathbb{P}_a}{\rho_{0a}^2} + \frac{\mathbb{P}_b}{\rho_{0b}^2} + \Pi_{0ab} \mathbb{I} \right) \nabla_0 \Phi_{0ab} \quad (9) $$

avec $\mathbb{P}$ le tenseur de contrainte PK1 et $\rho_0$ la densité du matériau. Si ce dernier est supposé homogène dans la configuration initiale de référence, l'Eq.(7) s'écrit alors
$$ \frac{dv_a}{dt} = \sum_{b=1}^{n} \left( \frac{P_a}{\rho_0} + \frac{P_b}{\rho_0} + \rho_0 \Pi_{0ab} \mathbb{I} \right) \nabla_0 \Phi_{0ab} \quad (10) $$

$$ \text{avec } \Pi_{0ab} = \frac{\mu}{\rho_0} (-\bar{c}_{ab} \vartheta_{ab} + 2 \vartheta_{ab}^2) $$

### 3 Formulation SPH pour les coques épaisses

La formulation SPHC adaptée pour les structures minces de type coque, est construite en se basant sur la théorie de Reissner-Mindlin pour la prise en compte du cisaillement transversal important dans le cas de coques épaisses. Dans ce cas, seule la surface moyenne de la coque est modélisée par une couche unique de particules pour représenter le comportement global de la structure. Chaque particule possèdera trois degrés de liberté (ddl) de type translation ($u_\xi, u_\eta, u_\zeta$) et deux ddl de type rotation ($\theta_\xi, \theta_\eta$) dans le plan moyen tangent à la coque.

Soit une coque d'épaisseur $t$, dans laquelle on repère une particle $p$ sur la surface moyenne de la coque (Fig. 1), par son vecteur position $X_p$. Le vecteur position $X_q = X$ d'un point $q$ situé à une distance $\zeta$ de la particle $p$, peut être donné par

$$ X(\xi, \eta, \zeta) = X_p(\xi, \eta) + \zeta V(\xi, \eta) \quad \text{avec} \quad \zeta \in \left\{ -\frac{t}{2}, +\frac{t}{2} \right\} \quad (11) $$

où $V(\xi, \eta)$ est le vecteur fibre épaisseur au niveau de la particule $p$. Le champ de déplacement $U_q = U$ du point $q$ peut être calculé de façon similaire

<figure><img src="image_2.png" /><figcaption>Fig. 1 – Cinématique SPH de coques épaisses</figcaption></figure>

$$ U(\xi, \eta, \zeta) = U_p(\xi, \eta) + \zeta U_n(\xi, \eta) \quad (12) $$

Pour la formulation de coques, il est nécessaire d'exprimer les déformations et les contraintes dans le repère tangent local orthonormé {$t_1, t_2, n$}, pour cela on définit la matrice $Q$ en tout point $q$ par

$$ Q = [t_1; t_2; n] \quad (13) $$

En utilisant Eq. (11), le vecteur gradient $dX$ peut être exprimé en terme de $d\xi = \{d\xi, d\eta, d\zeta\}$

$$ dX = \mathbb{F} \, d\xi \quad (14) $$

avec $\mathbb{F} = [a_{1\zeta}; a_{2\zeta}; a_{3\zeta}]$ le tenseur gradient de déformations exprimé en fonction des vecteurs de la base covariante ou naturelle

$$ \begin{aligned} a_{1\zeta} &= X_{p,\xi} + \zeta V_{,\xi} \\ a_{2\zeta} &= X_{p,\eta} + \zeta V_{,\eta} \\ a_{3\zeta} &= V \end{aligned} \quad (15) $$
Le vecteur gradient de déplacement dans la base curviligne peut être exprimé par

$$dU = L \, ds \quad \text{et} \quad L = F^T L_\zeta \qquad (16)$$

avec $dU = L_\zeta d\xi$. Enfin nous pouvons définir le tenseur de déformations dans la base covariante par

$$E_c = \frac{1}{2} (L^T + L + L^T L) \qquad (17)$$

Dans le but d'exprimer le tenseur de déformations directement dans la base curviligne, nous définissons le tenseur $C = F^{-1} Q$. On obtient finalement le tenseur de déformation de Green-Lagrange dans la base curviligne

$$E = C^T E_c C \qquad (18)$$

Ce dernier tenseur peut être décomposé en deux parties : une partie constante dans l'épaisseur $E_m$ associée aux effets de membrane et de cisaillement transversal et une partie linéaire dans l'épaisseur $E_f$ associée aux effets de flexion.

Le calcul des contraintes PK2 S locales associées aux déformations de Green-Lagrange $E$ est basé sur l'hypothèse d'un matériau élastique homogène et isotrope (module d'Young $E$ et coefficient de Poisson $\nu$) avec prise en compte de l'hypothèse des contraintes planes. L'intégration analytique des contraintes à travers l'épaisseur dans le cas d'un matériau élastique, permet d'aboutir aux 3 expressions des efforts normaux $N_\xi$, $N_\eta$, $N_{\xi\eta}$, 2 efforts tranchants $T_\xi$, $T_\eta$ et 3 moments fléchissants $M_\xi$, $M_\eta$, $M_{\xi\eta}$, qui peuvent être organisés sous forme matricielle

$$N = \begin{bmatrix} N_{\xi} & N_{\xi\eta} & T_{\xi} \\ N_{\eta\xi} & N_{\eta} & T_{\eta} \\ T_{\xi} & T_{\eta} & 0 \end{bmatrix} \quad \text{et} \quad M = \begin{bmatrix} M_{\xi} & M_{\xi\eta} & 0 \\ M_{\eta\xi} & M_{\eta} & 0 \\ 0 & 0 & 0 \end{bmatrix} \qquad (19)$$

Le tenseur des contraintes nominales $P$ dans la base globale est alors obtenu par

$$P = (Q S Q^T) F^T \qquad (20)$$

On en déduit finalement l'équation d'équilibre dynamique en membrane et en cisaillement transverse L'équation d'équilibre dynamique de les contraintes interieures est ensuite évolué dans le repère local initial,

$$\rho_0 \ddot{U} = \nabla_0 \cdot P^T \qquad (21)$$

Les deux équations d'équilibre des moments dans le repère local permettent de calculer les deux accélé- rations angulaires

$$I_{\xi} \ddot{\theta}_{\xi} = M_{\xi\eta,\xi} + M_{\eta,\eta} + T_{\eta}; \quad I_{\eta} \ddot{\theta}_{\eta} = M_{\xi,\xi} + M_{\xi\eta,\eta} - T_{\xi} \qquad (22)$$

où $I_\xi, I_\eta$ sont les inerties de rotation.

Les deux dernières équations peuvent s'écrire sous forme matricielle dans le repère global [11]

$$I \ddot{\theta} = \nabla_0 \cdot M_0^T + T_0 \qquad (23)$$

avec

$$T = \{T_{\eta}, -T_{\xi}, 0\}^T, M_0 = J C M Q^T \text{ et } T_0 = J Q^T T.$$

# 4 Résolution non linéaire en dynamique explicite

Dans cette étude, nous utilisons un schéma d'intégration explicite en temps afin de résoudre les équations d'équilibre Eq.(21) et Eq.(23). Nous pouvons ainsi obtenir à chaque pas de temps $\Delta t$, le vecteur incrément de déplacement incrémental $\Delta U$ et le vecteur incrément de rotations $\Delta \theta$.

Notons cependant que les vitesses angulaires $\dot{\theta}$ à différents pas de temps n'appartiennent pas au même espace tangent et ne sont donc pas additives directement. Il devient nécessaire alors d'adopter la matrice $W$ afin d'actualiser $\dot{\theta}$ et tenir compte de grandes rotations, telle que
$$ \dot{\theta}_{t+\Delta t} = \mathbb{W}(\Delta\theta) (\dot{\theta}_t + \Delta t \ddot{\theta}_t) \qquad (24) $$

avec $\mathbb{W}(\theta) = \frac{\sin \theta}{\theta} I + \frac{1 - \cos \theta}{\theta^2} \Xi + \frac{\theta - \sin \theta}{\theta^3} \theta \times \theta$, $\theta = \|\theta\|$ et $\Xi$ est le tenseur antisymétrique pour lequel $\theta$ est le vecteur axial. La mise à jour du vecteur pseudo normal $n$ est faite par

$$ \mathbf{n}_{t+\Delta t} = \mathbb{R}_{t+\Delta t} \cdot \mathbf{n}^0 \qquad (25) $$

où $\mathbb{R}$ est la matrice de Rodrigues. Nous terminons enfin par signaler que pour éviter un comportement irrégulier des normales moyennes de particules, en particulier sur les bords, nous introduisons le "filtre" d'homogénéisation suivant [1].

$$ \hat{\mathbf{n}}'_{a} = \mathbf{n}_{a} + \alpha \left( \frac{\sum_{b=1}^{n} \mathbf{n}_{b} \phi_{ab}}{\sum_{b=1}^{n} \phi_{ab}} - \mathbf{n}_{a} \right) \qquad (26) $$

où $\hat{\mathbf{n}}'$ est la nouvelle normale moyenne filtrée et $\alpha$ est la valeur du filtre fixée à 0.01.

# 5 Applications numériques

Dans cette section, nous présentons l'application de l'algorithme exposé ci-haut à travers deux exemples numériques en utilisant le schéma dynamique explicite. les résultats obtenus sont comparés avec des solutions de référence et à ceux obtenus par le logiciel ABAQUS® avec ses deux versions explicite et standard.

## 5.1 Post-flambement d'une colonne d'Euler

On considère dans un premier temps un exemple avec une solution analytique de flambement afin de vérifier l'efficacité et la pertinence de la méthode SPH que nous avons développée. La colonne d'Euler est de longueur $L = 200mm$, d'épaisseur $h = 5mm$ et de largeur $b = 10mm$. La structure est encastrée à ses deux bords et soumise à une force axiale de compression. La colonne est fabriquée d'un matériau de module d'Young $E = 210GPa$ et de coefficient de Poisson $\nu = 0.3$.

Cet exemple présente un flambement par point de bifurcation et la charge critique analytique est donnée par $P_c^{ref} = 4\pi^2 E I/L^2 = 21.59kN$.

Nous avons discrétisé la colonne avec 50 × 5 particules, puis nous avons effectué un pilotage en déplacement sur une durée de 30ms pour atteindre un grand déplacement maximal de 60mm. L'analyse non linéaire est effectuée par un algorithme de type dynamique explicite. La Figure 2 montre la déformée obtenue par la méthode SPH, pour un déplacement maximal de 60mm sachant que le flambement a lieu très tôt pour un déplacement inférieur à 2mm.

<figure><img src="image_3.png" /><figcaption>Fig. 2 – Configuration déformée SPH</figcaption></figure>

<caption>Fig. 3 – Courbes charge-déplacement</caption>

Nous pouvons remarquer sur la Figure 3 que l'allure générale de la courbe charge déplacement obtenue par la méthode SPH est proche de la solution analytique d'Euler. Nous voyons également que la valeur
de la charge critique obtenue par méthode SPH ($P_c^{SPH} = 21.75kN$) très proche de la solution analytique d'Euler à savoir $P_c^{ref} = 21.59kN$. Ceci témoigne de l'efficacité de la méthode SPH à prédire correctement le comportement non linéaire en post-flambement de coques.

## 5.2 Post-flambement d'un arc semi-circulaire sous charge

On considère une coque formée d'un arc semi-circulaire de rayon $R = 100mm$, d'épaisseur $h = 5mm$ et de largeur $b = 20mm$. La structure est encastrée à ses deux bords et soumise à une force concentrée en son sommet, comme l'indique la figure 4. Lencastré à deux extrémités et soumis à une force ponctuelle $P$ en son centre. Les données géométriques et celles de matériau sont présentées sur la Figure 4. Le choix du chargement est tel qu'il provoque un flambement par point limite, l'objectif est de trouver la solution non linéaire à ce problème par la méthode SPH.

<figure><img src="image_4.png" /><figcaption>Fig. 4 – Arc circulaire : géométrie et chargement</figcaption></figure>

<figure><img src="image_5.png" /><figcaption>Fig. 5 – Courbe de déplacement imposé en fonction du temps</figcaption></figure>

Dans un premier temps, nous avons résolu ce problème en utilisant ABAQUS® pour avoir une solution de référence. Un maillage de 5x47 éléments de coque SC4R est alors utilisé. Ce même maillage est utilisé pour la génération de la discrétisation en particules SPH en ne gardant que les noeuds précédemment générés.

Un pilotage en longueur d'arc est imposé, cette technique est la plus adaptée pour les instabilités de type flambage. Pour toutes nos simulations, le déplacement imposé augmente graduellement pendant 3ms jusqu'à atteindre 100% de sa valeur ou 110mm tel que montré sur la courbe présentée dans la Figure 5. Les Figures 6 et 7 montrent les deux déformées obtenues respectivement par la méthode SPH et par ABAQUS, au point limite juste avant flambage. Comme nous pouvons l'observer les deux déformées sont identiques.

<figure><img src="image_6.png" /><figcaption>Fig. 6 – Configuration déformée SPH</figcaption></figure>

<caption>Fig. 7 – Configuration déformée ABAQUS</caption>

La Figure 8 montre une comparaison des courbes charge-déplacement au sommet de l'arc. Nous pouvons voir que la solution obtenue par notre modèle SPH ($P_c^{SPH} = 23.56kN$) est en parfait accord avec la solution obtenue par ABAQUS Standard ($P_c^{ABQ-Std} = 23.60kN$). Nous remarquons également que la
solution obtenue par ABAQUS Explicite est en légère baisse par rapport aux deux précédentes solution avec une charge critique $P_c^{ABQ-Exp} = 22.18kN$ soit avec une erreur de moins de 6%. Enfin, nous remarquons que la courbe de solution SPH coincide parfaitement avec celle de ABAQUS Standard, ce qui prouve la robustesse du modèle SPH même en utilisant un algorithme dynamique explicite.

<figure><img src="image_7.png" /><figcaption>Fig. 8 – Comparaison des courbes charge-déplacement au sommet de l'arc</figcaption></figure>

# 6 Conclusion

Nous avons présenté dans cet article, une analyse de post-flambement de structures minces en utilisant une formulation SPH de type coque. Les deux principaux inconvénients de la méthode SPH classique ont été éliminés à l'aide de fonctions noyau corrigées et d'une FLT. En se basanr sur la théorie de Reissner-Mindlin, un modèle SPH coque a été élaboré reproduisant le comportement d'une coque modélisée par une seule couche de particules placées sur sa surface moyenne.

Plusieurs applications numériques ont été traitées et les résultats comparés à ceux obtenus par ABAQUS. La bonne qualité des résultats obtenus, prouve que la méthode SPH peut constituer une alternative à la méthode EF.

# Références

[1] D.S. Balsara, von Neumann stability analysis of smoothed particle hydrodynamics-suggestions for optimal algorithms, J. Comput. Phys. 121 1995, pp. 357-372.

[2] T. Belytschko., Y. Guo, W. K. Liu and S.P. Xiao. A unified stability analysis of meshless particle methods. International Journal for Numerical Methods in Engineering, 2000, 48 pp. 1359–1400.

[3] J.K. Chen, J.E. Beraun, T.C. Carney, A corrective smoothed particle method for boundary value problems in heat conduction, Internat. J. Numer. Meth. Engrg. 1999, 46, pp. 231–252.

[4] C. T. Dyka, P. W. Randies, R. P. Ingel. Stress points for tension instability in smoothed particle hydrodynamics, International Journal for Numerical Methods in Engrgineering, 1997, 40, pp. 2325-2341.

[5] R.A. Gingold, J.J. Monaghan. Smoothed particle hydrodynamics : theory and application to non-spherical stars. MonthlyNotice of the Royal Astronomical Society. 1977(181), pp. 375-389.

[6] J.P. Gray, J.J. Monaghan, R.P. Swift. SPH elastic dynamics. Computer Methods in Applied Mechanics and Engineering. 2001, 190, pp. 6641–6662.

[7] P. Krysl and T. Belytschko. Analysis of thin shells by the element free Galerkin method, Int. J. Numer. Meth. Engng. 1996(33), pp. 3057-3078.

[8] J. Limido, C. Espinosa, M. Salaün and J.L. Lacome, SPH method applied to high speed cutting modelling. Int. J. Mech. Sci. 4 (7), 2007, pp. 898-908.

[9] J. Liu, Meshless study of dynamic failure in shells, J. Eng. Math. 71(2), 2011, pp. 205-222.

[10] L.B. Lucy. A numerical approach to the testing of the fission hypothesis. Journal of Astronomy. 1977(82), pp. 1013-1024.

[11] B. Maurel, A. Combescure. An SPH shell formulation for plasticity and fracture analysis in explicit dynamics. International Journal for Numerical Methods in Engineering, 2008, 76, pp. 949–971.

[12] J.J. Monaghan, Particle methods for hydrodynamics, Comput. Phys. Rep. 3, 1985, pp. 71-124.