# Résolution de l'équation des plaques minces de Kirchhoff par une approche sans maillage construite à l'aide d'une approximation moindres carrés mobiles

Pierre Joyot, Francisco Chinesta, Pierre Villon

► To cite this version:

Pierre Joyot, Francisco Chinesta, Pierre Villon. Résolution de l'équation des plaques minces de Kirchhoff par une approche sans maillage construite à l'aide d'une approximation moindres carrés mobiles. 9e Colloque national en calcul des structures, CSMA, May 2009, Giens, France. hal-01408415

HAL Id: hal-01408415

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

Submitted on 4 Dec 2016

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
# Résolution de l'équation des plaques minces de Kirchoff par une approche sans maillage construite à l'aide d'une approximation moindres carrés mobiles

P. Joyot¹, F. Chinesta², P. Villon³

¹ ESTIA-Recherche
technopole izarbel, 64210 Bidart
pjoyot@estia.fr

² EADS Corporate Foundation International Chair
GEM CNRS-ECN
1 rue de la Noë BP 92101, 44321 Nantes cedex 3
Francisco.Chinesta@ec-nantes.fr

³ UTC-Roberval UMR 6253
60200 Compiègne
pierre.villon@utc.fr

**Résumé** — Nous présentons une méthode sans maillage de type EFG pour résoudre un problème de plaque mince de Kirchoff. Nous utilisons une approximation Hermite MLS, un schéma d'intégration HSCNI, et une méthode de pénalisation modifiée pour appliquer les conditions aux limites. De cette manière, nous obtenons des taux de convergence optimaux pour des bases d'approximations quadratique ou cubique.

**Mots clés** — Méthode sans maillage, MLS Hermite, plaque mince, Kirchoff

## 1 Introduction

Nous souhaitons résoudre l'équation biharmonique caractérisant le comportement des plaques minces suivant l'hypothèse de Kirchoff, à l'aide d'une technique sans maillage de type AD[1], EFG, RKPM. Rappelons que cette approche consiste à discrétiser la forme faible obtenue par la méthode de Galerkin, en utilisant une approximation moindres carrés mobiles (MLS). L'approximation MLS présente l'intérêt d'être facilement dérivable pour un ordre élevé. Cette propriété permet d'envisager de traiter l'équation biharmonique, qui nous intéresse ici, de manière directe, avec la même facilité que l'équation de Poisson, tout en gardant des propriétés de convergence optimales. Ce dernier point semble être un avantage en comparaison de l'approche éléments finis. Cette intéressante propriété a conduit à un nombre important de publications. Citons par ordre chronologique [2][3][4], ou plus récemment [5]. Malgré tout, les différentes approches proposées ne permettent pas d'obtenir les vitesses de convergence optimales escomptées. Ceci est vrai pour une approximation quadratique, la situation étant pire pour l'approximation cubique qui donne généralement des résultats moins précis que le cas quadratique.

L'objectif de notre travail est de comprendre l'origine de la perte de précision et de proposer les solutions adéquates.
Dans la première partie, nous présenterons la formulation Hermite MLS, ensuite nous préciserons le schéma d'intégration et la méthode utilisée pour appliquer les conditions aux limites. Pour finir, nous donnerons des résultats de convergence pour les bases polynomiales quadratique et cubique.

## 2 L'approximation Moindres Carrés Mobiles Hermite (HMLS)

L'approximation Hermite MLS a été introduite par Pierre Villon [6] avec pour finalité le remaillage de surface. En parallèle Atluri [3] a proposé de résoudre le problème des poutres de Bernoulli avec l'approximation GMLS, dont la formulation est très proche de l'approximation HMLS. Très récemment J.S. Chen a proposé l'approximation HRKA [5] qui est identique à l'approximation GMLS. Il l'a appliqué à la résolution de problème des plaques de Kirchhoff.

Le but de cette partie est de présenter l'approximation HMLS, et de montrer l'importance du paramètre d'adimensionnement dont le rôle prépondérant n'a pas été évoqué dans les publications précédentes.

Soit $u_{ex}(\mathbf{x})$ un champ que l'on souhaite reconstruire à partir de valeurs définies en un nombre fini de noeuds. Nous noterons $\chi = \{ \mathbf{x}_1 \cdots \mathbf{x}_{n_d} \}$ cet ensemble de noeuds, avec $n_d$ le nombre de noeuds.

Dans le cas d'une approximation 2D nous aurons 3 valeurs pour chaque noeud. La première valeur definit le champ considéré, les deux autres sont liées à la dérivée de ce champ suivant les variables d'espaces. Nous préciserons par la suite la signification exacte de ces valeurs. Nous représenterons les données nodales en $\bar{\mathbf{x}}_i$ par le vecteur $\mathbf{u}_i^\top = [u_i^{00} \ u_i^{10} \ u_i^{01}]$.

Les valeurs nodales pour l'ensemble des noeuds $\chi$ sont représentées par le vecteur $\mathcal{U}$ défini par $\mathcal{U}^\top = [\mathcal{U}^{00} \ \mathcal{U}^{10} \ \mathcal{U}^{01}]$ avec $\mathcal{U}^\star = [u_1^\star \ \cdots \ u_{n_d}^\star]$ pour $\star \in \{00, 10, 01\}$.

La fonction reconstruite à partir des valeurs modales est notée $u_{app}(\mathbf{x})$. Elle est définie par

$$u_{app}(\mathbf{x}) = \sum_{i=1}^{n_d} N_i^{00}(\mathbf{x}) \mathbf{u}_i^{00} + N_i^{10}(\mathbf{x}) \mathbf{u}_i^{10} + N_i^{01}(\mathbf{x}) \mathbf{u}_i^{01} \quad (1)$$

où $N_i^{00}(\mathbf{x})$, $N_i^{10}(\mathbf{x})$, $N_i^{01}(\mathbf{x})$, sont les fonctions de forme du noeud $i$ que nous allons définir par la suite.

Nous définissons le vecteur des fonctions de forme calculés pour l'ensemble des noeuds par $\mathcal{N}^\top(\mathbf{x}) = [\mathcal{N}^{00\top}(\mathbf{x}) \ \mathcal{N}^{10\top}(\mathbf{x}) \ \mathcal{N}^{01\top}(\mathbf{x})]$ avec $\mathcal{N}^\star(\mathbf{x}) = [N_1^\star(\mathbf{x}) \ \cdots \ N_{n_d}^\star(\mathbf{x})]$ pour $\star \in \{00, 10, 01\}$

Ainsi, $u_{app}(\mathbf{x})$ est également défini par

$$u_{app}(\mathbf{x}) = \mathcal{N}^\top(\mathbf{x}) \mathcal{U} \quad (2)$$

Afin de décrire la méthode permettant de déterminer les fonctions de forme nous introduisons la notion d'approximation locale.

Considérons un point d'évaluation $\mathbf{x}(x, y)$ et soit $\bar{\mathbf{x}}(\bar{x}, \bar{y})$, un autre point situé dans le voisinage de $\mathbf{x}$.

Nous définissons l'approximation locale $v^{00}(\mathbf{x}, \bar{\mathbf{x}})$ centrée en $\mathbf{x}$ et évaluée en $\bar{\mathbf{x}}$ par :

$$v^{00}(\mathbf{x}, \bar{\mathbf{x}}) = \mathbf{P}^{00\top} \left( \frac{\bar{\mathbf{x}} - \mathbf{x}}{a} \right) \theta(\mathbf{x}) \quad (3)$$

ou $\mathbf{P}^{00}(\mathbf{x})$ est une base polynomiale, $a$ est un paramètre d'adimensionnement, $\theta(\mathbf{x})$ est un vecteur à déterminer.

La fonction reconstruite $u_{app}(\mathbf{x})$ correspond à l'évaluation de l'approximation locale en $\mathbf{x}$

$$u_{app}(\mathbf{x}) = v^{00}(\mathbf{x}, \bar{\mathbf{x}}) = \mathbf{P}^{00\top}(0) \theta(\mathbf{x}) \quad (4)$$
$u_{app}(\mathbf{x})$ sera définie si $\theta(\mathbf{x})$ est connu.

Précisons les dérivées de l'approximation locale

$$ \frac{\partial}{\partial \bar{x}} v^{00}(\mathbf{x}, \bar{\mathbf{x}}) = \frac{\partial}{\partial \bar{x}} \mathbf{P}^{00\top} \left( \frac{\bar{\mathbf{x}} - \mathbf{x}}{a} \right) \theta(\mathbf{x}) = \frac{1}{a} \mathbf{P}^{10\top} \left( \frac{\bar{\mathbf{x}} - \mathbf{x}}{a} \right) \theta(\mathbf{x}) \quad (5) $$

$$ \frac{\partial}{\partial \bar{y}} v^{00}(\mathbf{x}, \bar{\mathbf{x}}) = \frac{\partial}{\partial \bar{y}} \mathbf{P}^{00\top} \left( \frac{\bar{\mathbf{x}} - \mathbf{x}}{a} \right) \theta(\mathbf{x}) = \frac{1}{a} \mathbf{P}^{01\top} \left( \frac{\bar{\mathbf{x}} - \mathbf{x}}{a} \right) \theta(\mathbf{x}) \quad (6) $$

Notons que $\mathbf{P}^{10\top}(\frac{\bar{\mathbf{x}}-\mathbf{x}}{a})$ est la dérivée de $\mathbf{P}^{00\top}(\frac{\bar{\mathbf{x}}-\mathbf{x}}{a})$ multiplié par le paramètre $a$. Nous avons volontairement mis le terme $\frac{1}{a}$ en facteur afin que le vecteur $\mathbf{P}^{10\top}(\frac{\bar{\mathbf{x}}-\mathbf{x}}{a})$ ne contienne que des puissances de $\frac{\bar{\mathbf{x}}_i-\mathbf{x}}{a}$ qui resteront insensibles à la taille du domaine, à condition que le paramètre $a$ soit correctement choisi. Ceci est essentiel pour s'assurer de la convergence de l'approximation.

Nous pouvons faire la même remarque en ce qui concerne le vecteur $\mathbf{P}^{01\top}(\frac{\bar{\mathbf{x}}-\mathbf{x}}{a})$. Nous notons $v^{10}(\mathbf{x}, \bar{\mathbf{x}}) = a \frac{\partial}{\partial \bar{x}} v^{00}(\mathbf{x}, \bar{\mathbf{x}})$ et $v^{01}(\mathbf{x}, \bar{\mathbf{x}}) = a \frac{\partial}{\partial \bar{x}} v^{10}(\mathbf{x}, \bar{\mathbf{x}})$ les dérivées diffuses adimensionnées.

Nous pouvons maintenant définir la fonctionnelle, notée $J(\theta(\mathbf{x}))$, permettant de calculer $\theta(\mathbf{x})$ telle que la distance entre la fonction reconstruite, ces dérivées diffuses adimensionnées et les valeurs nodales, soit minimisées. Nous avons

$$ \begin{aligned} J(\theta(\mathbf{x})) = & \frac{1}{2} \sum_{i=1}^{n_d} (v^{00}(\mathbf{x}, \bar{\mathbf{x}}_i) - \mathbf{u}_i^{00})^2 w_{r_i}^{00}(\mathbf{x}, \bar{\mathbf{x}}_i) \\ & + \frac{1}{2} \sum_{i=1}^{n_d} (v^{10}(\mathbf{x}, \bar{\mathbf{x}}_i) - \mathbf{u}_i^{10})^2 w_{r_i}^{10}(\mathbf{x}, \bar{\mathbf{x}}_i) \\ & + \frac{1}{2} \sum_{i=1}^{n_d} (v^{01}(\mathbf{x}, \bar{\mathbf{x}}_i) - \mathbf{u}_i^{01})^2 w_{r_i}^{01}(\mathbf{x}, \bar{\mathbf{x}}_i) \end{aligned} \quad (7) $$

Définissons $\underline{\mathbf{P}}^\top = [ \underline{\mathbf{P}}^{00\top} \underline{\mathbf{P}}^{10\top} \underline{\mathbf{P}}^{01\top} ]$ avec $\underline{\mathbf{P}}^\top = [ \underline{\mathbf{P}}^\star(\frac{\bar{\mathbf{x}}_1-\mathbf{x}}{a}) \dots \underline{\mathbf{P}}^\star(\frac{\bar{\mathbf{x}}_{n_d}-\mathbf{x}}{a}) ]$ pour $\star \in \{00, 10, 01\}$. La fonctionnelle $J(\theta(\mathbf{x}))$ se réécrit

$$ J(\theta(\mathbf{x})) = \frac{1}{2} (\underline{\mathbf{P}}(\mathbf{x}) \theta(\mathbf{x}) - \underline{\mathbf{U}})^{\top} \underline{\mathbf{W}}(\mathbf{x}) (\underline{\mathbf{P}}(\mathbf{x}) \theta(\mathbf{x}) - \underline{\mathbf{U}}) \quad (8) $$

Le vecteur $\theta(\mathbf{x})$ qui minimise $J(\theta(\mathbf{x}))$ est obtenu en résolvant le système

$$ \mathbf{M}(\mathbf{x}) \theta(\mathbf{x}) = \underline{\mathbf{P}}^{\top}(\mathbf{x}) \underline{\mathbf{W}}(\mathbf{x}) \underline{\mathbf{U}} \quad (9) $$

où, $\mathbf{M}(\mathbf{x})$ est la matrice moment définie par

$$ \begin{aligned} \mathbf{M}(\mathbf{x}) &= \underline{\mathbf{P}}^{\top}(\mathbf{x}) \underline{\mathbf{W}}(\mathbf{x}) \underline{\mathbf{P}}(\mathbf{x}) \\ &= \underline{\mathbf{P}}^{00\top}(\mathbf{x}) \underline{\mathbf{W}}^{00}(\mathbf{x}) \underline{\mathbf{P}}^{00}(\mathbf{x}) + \underline{\mathbf{P}}^{10\top}(\mathbf{x}) \underline{\mathbf{W}}^{10}(\mathbf{x}) \underline{\mathbf{P}}^{10}(\mathbf{x}) + \underline{\mathbf{P}}^{01\top}(\mathbf{x}) \underline{\mathbf{W}}^{01}(\mathbf{x}) \underline{\mathbf{P}}^{01}(\mathbf{x}) \end{aligned} \quad (10) $$

Les équations 2, 4, 9, permettent finalement de définir les fonctions de forme Hermite MLS

$$ \underline{\mathbf{N}}^{\star}(\mathbf{x}) = \underline{\mathbf{W}}^{\star}(\mathbf{x}) \underline{\mathbf{P}}^{\star}(\mathbf{x}) \mathbf{M}^{-1}(\mathbf{x}) \mathbf{P}^{00}(0) \quad \text{pour } \star \in \{00, 10, 01\} \quad (11) $$

*Remarque :* Le paramètre d'adimensionnement est essentiel pour s'assurer de la convergence de l'approximation. Les fonctions de forme MLS ne dépendent pas de ce paramètre, ce n'est pas le cas de l'approximation Hermite MLS. Ce dernier point n'a jamais été précisé dans la littérature.

*Remarque :* Les valeurs nodales de l'approximation sont reliées à la valeur de la fonction et à ces deux dérivées diffuses multipliées par le paramètre d'adimensionnement. Les valeurs nodales sont donc liées au paramètre d'adimensionnement ce qui n'est pas le cas des fonctions de forme MLS classiques.
### 3 Équations des plaques minces

Soit $\Omega \subset \mathbb{R}^2$ un domaine représentant le plan moyen d'une plaque d'épaisseur $h$. Soit $\Gamma$ la frontière de $\Omega$. Cette frontière est coupée en plusieurs parties telles que $\Gamma = \Gamma_w \cup \Gamma_Q = \Gamma_{\theta_n} \cup \Gamma_{M_n}$ et $\Gamma_w \cap \Gamma_Q = \Gamma_{\theta_n} \cap \Gamma_{M_n} = \emptyset$. On définit un système de coordonnées cartésien $x, y, z$ telque $z$ soit normal à la surface moyenne de la plaque. Dans ce qui suit, $\alpha, \beta$ seront des indices prenant les valeurs $x$ et $y$. $\mathbf{n}$ (resp. $\mathbf{s}$) représente la normale sortante (resp. tangente) unitaire de $\Gamma$.

La forme forte du problème de la plaque en flexion suivant les hypothèses de Kirchhoff est définie de la manière suivante. Connaissant $q$ une force distribuée, $\bar{w}$ le déplacement transverse sur $\Gamma_w$, $\bar{\theta}_n$ la rotation normale sur $\Gamma_{\theta_n}$, $\bar{Q}$ la force transverse sur $\Gamma_Q$, $\bar{M}_n$ le moment normal sur $\Gamma_{M_n}$, trouver $w$ tel que

$$m_{\alpha\beta,\alpha\beta} = q \quad \text{sur } \Omega \qquad (12)$$

$$w = \bar{w} \quad \text{sur } \Gamma_w \qquad (13)$$

$$\theta_n = w, \alpha n_\alpha = \bar{\theta}_n \quad \text{sur } \Gamma_{\theta_n} \qquad (14)$$

$$M_n = m_{\alpha\beta} n_\alpha n_\beta = \bar{M}_n \quad \text{sur } \Gamma_{M_n} \qquad (15)$$

$$Q = -m_{\alpha\beta,\beta} n_\alpha - (m_{\alpha\beta} n_\beta s_\alpha)_{,s} = \bar{Q} \quad \text{sur } \Gamma_Q \qquad (16)$$

avec $w \in H^4(\Omega)$ quand $q \in H^2(\Omega)$.

La fonction $w$, solution de la forme forte, est également le minimum de l'énergie potentielle définie par

$$\pi(w) = \pi_{\text{int}}(w) - (\pi_{\text{load}}(w) + \pi_{\text{bend}}(w)) \qquad (17)$$

où, $\pi_{\text{int}}(w)$ est l'énergie interne

$$\pi_{\text{int}}(w) = \frac{1}{2} \int_{\Omega} \mathbf{b}_{\nabla}^{\top}(w) \mathbf{D} \mathbf{b}_{\nabla}(w) \, d\Omega \qquad (18)$$

$$\mathbf{b}_{\nabla}^{\top}(w) = [ w_{,xx} \quad w_{,yy} \quad 2w_{,xy} ] \qquad (19)$$

$$\mathbf{D} = \begin{bmatrix} 1 & \mathbf{v} & 0 \\ \mathbf{v} & 1 & 0 \\ 0 & 0 & \frac{1-\mathbf{v}}{2} \end{bmatrix} \qquad (20)$$

$\pi_{\text{load}}(w)$ est l'énergie créée par la force distribuée

$$\pi_{\text{load}}(w) = \int_{\Omega} qw \, d\Omega \qquad (21)$$

et $\pi_{\text{bend}}(w)$ est l'énergie de flexion créée par les moments et les efforts tranchants

$$\pi_{\text{bend}}(w) = \int_{\Gamma_Q} \bar{Q}w \, d\Gamma - \int_{\Gamma_{M_n}} \bar{M}_n w_{,n} \, d\Gamma \qquad (22)$$

L'application des conditions aux limites par pénalisation nous amène à introduire la forme augmentée de l'énergie potentielle. $\lambda_{\bar{w}}$ et $\lambda_{\bar{\theta}_n}$ sont les coefficients de pénalisation. Afin de résoudre numériquement le problème, nous introduisons la forme discrétisée du déplacement transverse $w = \mathbf{N}^{\top}\mathbf{U}$, où $\mathbf{U}$ représente les inconnues nodales et $\mathbf{N}$ les fonctions de forme Hermite MLS décrites précédemment. L'énergie potentielle augmentée discrétisée s'écrit

$$\pi(\mathbf{U}) = \frac{1}{2} \mathbf{U}^{\top} \mathbf{K} \mathbf{U} - \mathbf{U}^{\top} \mathbf{f} + \frac{1}{2} \lambda_{\bar{w}} \int_{\Gamma_w} (\mathbf{N}^{\top} \mathbf{U} - \bar{w})^{\top} (\mathbf{N}^{\top} \mathbf{U} - \bar{w}) \, d\Gamma \qquad (23) \\ + \frac{1}{2} \lambda_{\bar{\theta}_n} \int_{\Gamma_{\theta_n}} (\mathbf{N}_{,n}^{\top} \mathbf{U} - \bar{\theta}_n)^{\top} (\mathbf{N}_{,n}^{\top} \mathbf{U} - \bar{\theta}_n) \, d\Gamma$$
avec

$$ \mathbf{K} = \int_{\Omega} \mathbf{B}_{\nabla}^{\top} (\mathbf{N}) \mathbf{D} \mathbf{B}_{\nabla} (\mathbf{N}) \, d\Omega \qquad (24) $$

$$ \mathbf{B}_{\nabla}^{\top}(\mathbf{N}) = \left[ \begin{array}{ccc} \mathbf{N}_{,xx} & \mathbf{N}_{,yy} & 2\mathbf{N}_{,xy} \end{array} \right] \qquad (25) $$

$$ \mathbf{f} = \int_{\Omega} \mathbf{N}q \, d\Omega + \int_{\Gamma_Q} \mathbf{N}\bar{Q} \, d\Gamma - \int_{\Gamma_{M_n}} \mathbf{N}_{,n}\bar{M}_n \, d\Gamma \qquad (26) $$

La minimisation, par rapport à $\mathcal{U}$, de l'énergie potentielle conduit au système linéaire suivant

$$ (\mathbf{K} + \lambda_{\tilde{w}} \mathbf{K}_{\tilde{w}} + \lambda_{\tilde{\theta}_n} \mathbf{K}_{\tilde{\theta}_n}) \mathcal{U} = \mathbf{f} + \lambda_{\tilde{w}} \mathbf{f}_{\tilde{w}} + \lambda_{\tilde{\theta}_n} \mathbf{f}_{\tilde{\theta}_n} \qquad (27) $$

# 4 Résultats numériques

La qualité des résultats de la méthode présentée est sensible au choix du schéma d'intégration et à la manière d'appliquer les conditions aux limites. En effet, il n'existe pas de méthode permettant d'intégrer exactement les matrices du système linéaire 27 quand des fonctions de forme MLS sont utilisées. Une solution à ce problème est d'imposer au schéma d'intégration d'être exact dans le cas où le champ solution est une combinaison linéaire des éléments contenus dans la base d'approximation MLS. Ainsi la satisfaction des contraintes obtenues permet d'atteindre le taux de convergence optimal dans le cas général. Dans [5], Wang donne l'expression de ces contraintes dans le cas du problème des plaques de Kirchhoff.

Nous avons comparé les performances de trois schémas d'intégration. Nous définissons des cellules d'intégration dont la taille correspond à la distance entre deux noeuds. Ces cellules peuvent être divisées en sous cellules. Nous obtenons les résultats suivants : Le schéma d'intégration de Gauss ne satisfait pas les contraintes d'intégration. Malgré tout, il est possible d'obtenir un taux de convergence optimal à condition que le nombre de points de Gauss soit très grand (5x5 sous cellules d'intégration avec 5x5 points d'intégration par sous cellules), le temps de calcul obtenu est prohibitif.

Le schéma d'intégration SCNI [5] satisfait les contraintes d'intégration. Dans ce cas, l'intégration est réalisée uniquement sur le bord de la cellule d'intégration, le gradient des fonctions de forme est considéré constant. Pour cette raison, le taux de convergence en norme L2 est limité à 2 même si les fonctions de forme sont quadratiques ou cubiques. De plus, dans le cas de l'opérateur biharmonique nous devons utiliser au moins 3x3 sous-divisions, pour avoir une précision acceptable.

Le schéma d'intégration HSCNI [7] satisfait également les contraintes d'intégration. Dans ce cas, l'intégration est réalisée sur le bord et à l'intérieur de la cellule d'intégration. Ce schéma semble donc plus coûteux. Il permet d'atteindre la convergence optimale. Notons que cette approche n'a jamais été testée avec l'opérateur biharmonique. Notre expérience montre qu'il est possible d'obtenir un schéma d'intégration précis avec seulement 2x2 sous-divisions et 3x3 points de Gauss par sous-divisions.

L'application des conditions aux limites dans le cas d'une approximation MLS est difficile pour deux raisons. Premièrement, cette approximation n'est pas interpolante, deuxièmement l'application des conditions aux limites aux noeuds n'implique pas leur respect pour les points de la frontière situés entre deux noeuds. Nous obtenons de bons résultats avec la méthode de pénalisation, à condition de prendre en compte le phénomène de sur-contraintes, décrit dans [8] pour l'équation de Poisson, lié au caractère non interpolant des fonctions de forme. Nous proposons une solution nouvelle adaptée au cas de l'équation biharmonique. Afin d'éliminer les équations parasites, il est nécessaire que le nombre de valeurs propres de la matrice $\mathbf{K}_{\tilde{w}}$ (resp. $\mathbf{K}_{\tilde{\theta}_n}$) soit égal au nombre de noeuds portant des conditions aux limites en déplacement (resp. rotation). Notre
solution consiste à reconstruire, à partir d'une décomposition SVD, les matrices $K_{\tilde{w}}$ et $K_{\tilde{\theta}_n}$ en ne gardant que les modes corespondant aux valeurs propres les plus grandes.

Nous avons appliqué notre schéma de résolution pour simuler une plaque carrée, avec ses bords encastrés, et soumise à une pression uniforme. Nous avons choisi une distribution de noeuds ré-gulière. Nous avons utilisé une base HMLS quadratique ou cubique. Nous avons obtenu en norme $L^2$ (resp. $H^1$) un taux de convergence de 3 (resp 2) pour une base quadratique et 4 (resp 3) pour une base cubique. Nous avons également constaté que l'approximation HMLS est plus précise que l'approximation MLS d'un ordre de grandeur. L'approximation HMLS permet d'appliquer directement les conditions aux limites en déplacement et en rotation, mais les résultats obtenus ne sont pas optimaux.

# 5 Conclusion

Le schéma de résolution présenté permet pour la première fois de résoudre de manière optimale et efficace le problème de la plaque de Kirchhoff avec des fonctions d'approximations HMLS quadratique ou cubique. La suite de notre travail consistera à tester l'approche de Nitsche [9] pour appliquer les conditions aux limites.

# Références

[1] B. Nayroles, G. Touzot, and P. Villon. Generalizing the finite element method : Diffuse approximation and diffuse elements. *Computational Mechanics*, 10 :307-318, 1992.

[2] P. Krysl and T. Belytschko ER. Analysis of thin plates by the element-free galerkin method. *Computational Mechanics*, 17(1) :26-35, 12 1995.

[3] S. N. Atluri, J. Y. Cho, and H.-G. Kim ER. Analysis of thin beams, using the meshless local petrov-galerkin method, with generalized moving least squares interpolations. *Computational Mechanics*, 24(5) :334-347, 11 1999.

[4] Suleiman BaniHani and Suvranu De. Development of a genetic algorithm-based lookup table approach for efficient numerical integration in the method of finite spheres with application to the solution of thin beam and plate problems. *International Journal for Numerical Methods in Engineering*, 67(12) :1700-1729, 2006.

[5] Dongdong Wang and Jiun-Shyan Chen. A hermite reproducing kernel approximation for thin-plate analysis with sub-domain stabilized conforming integration. *International Journal for Numerical Methods in Engineering*, 74(3) :368-390, 07 2007.

[6] A. Rassineux, P. Villon, J.-M. Savignat, and O. Stab. Surface remeshing by local hermite diffuse interpolation. *International Journal for Numerical Methods in Engineering*, 49(1-2) :31-49, 2000.

[7] J.S. Chen, W. Hu, and M. Puso. Orbital hp-clouds for solving schrödinger equation in quantum mechanics. *Computer Methods in Applied Mechanics and Engineering*, 196(37-40) :3693-3705, 08 2007.

[8] Jin Yeon Cho, You Me Song, and Yun Hyuk Choi. Boundary locking induced by penalty enforcement of essential boundary conditions in mesh-free methods. *Computer Methods in Applied Mechanics and Engineering*, 197(13-16) :1167-1183, 02 2008.

[9] Sonia Fernandez-Mendez and Antonio Huerta. Imposing essential boundary conditions in mesh-free methods. *Computer Methods in Applied Mechanics and Engineering*, 193(12-14) :1257-1275, 03 2004.