# Une méthode asymptotique numérique et une méthode sans maillage pour le calcul des points de bifurcation

Abdeljalil Tri, Hamid Zahrouni, Michel Potier-Ferry

► To cite this version:

Abdeljalil Tri, Hamid Zahrouni, Michel Potier-Ferry. Une méthode asymptotique numérique et une méthode sans maillage pour le calcul des points de bifurcation. 11e colloque national en calcul des structures, CSMA, May 2013, Giens, France. hal-01722063

HAL Id: hal-01722063

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

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
# UNE METHODE ASYMPTOTIQUE NUMERIQUE ET UNE METHODE SANS MAILLAGE POUR LE CALCUL DES POINTS DE BIFURCATION

A. Tri¹, H. Zahrouni², M. Potier-Ferry²

1 Laboratoire de Mécanique, Faculté des Sciences Aïn Chok Université Hassan II Casablanca, Casablanca, Maroc
jaliltri@yahoo.fr

2 Laboratoire d'Étude des Microstructures et de Mécanique des Matériaux, LEM3, CNRS, UMR 7239,
Université de lorraine – Metz, Ile du Saulcy, 57045 Metz, France
{hamid.zahrouni, michel.potier-ferry}@univ-lorraine.fr

**Résumé** — On propose un nouvel algorithme basé sur la combinaison de la Méthode Asymptotique Numérique (MAN) et d'une méthode sans maillage pour détecter des points de bifurcation sur une branche de solution régie par les équations de Poisson non linéaires. La méthode sans maillage est basée sur les méthodes des solutions fondamentales combinées aux méthodes des solutions particulières (MFS-MPS).

**Mots clés** — Méthode asymptotique numérique- méthode sans maillage - Méthode de solution fondamentale (MFS)- Méthode de solution particulière (MPS)-Indicateur de bifurcation

## 1. Introduction

Trois méthodes, basées principalement sur la méthode asymptotique numérique (MAN) sont mises au point pour détecter les points de bifurcations sur une branche de solution d'un problème de Poisson non linéaire. La MAN consiste en l'association des séries de Taylor et d'une méthode de discrétisation spatiale. Associée à la méthode de discrétisation par éléments finis, l'efficacité de la MAN n'est plus à démontrer pour la résolution des problèmes non linéaires et le suivi des branches de solution aussi bien en mécanique des fluides qu'en mécanique des solides [1-8]. L'un des avantages de cette méthode est d'obtenir une branche de solution analytique dont le domaine de validité est limité par le rayon de convergence des séries entières. La méthode de continuation permet de calculer entièrement la branche de solution pas à pas [2]. La longueur de chaque pas, déterminée a posteriori, varie selon la non linéarité locale de la branche. En effet, on constate une diminution de la longueur des pas au voisinage d'une forte nonlinéarité et une accumulation des pas au voisinage d'une bifurcation. Il a été établi que le rayon de convergence des séries est fortement lié à la distance entre le point de départ de la série et le point de bifurcation [4,6]. La première méthode proposée consiste alors à détecter le point de bifurcation géométriquement en localisant ces zones d'accumulation des pas.

La transformation des séries entières en fractions rationnelles, dites Approximants de Padé, améliore considérablement le rayon de convergence. Il a été montré que la première racine réelle du dénominateur de ces approximants donne exactement le point de bifurcation [4]. C'est notre deuxième technique développée pour la détection des points de bifurcation.

La troisième méthode consiste à introduire un indicateur de bifurcation, le long de la branche de solution non linéaire, qui s'annule aux points singuliers. Cette technique est introduite initialement par Boutyour [3] en mécanique des solides et par Tri en mécanique des fluides [5].
L'originalité de ce travail est d'associer ces indicateurs de bifurcation à une méthode de discrétisation sans maillage. Ces dernières méthodes ont été développées dans le but de pallier certaines difficultés liées à la distorsion du maillage par éléments finis. Parmi ces méthodes, on peut citer la méthode des solutions fondamentales (MFS) [9] qui a prouvé son efficacité pour la résolution des équations différentielles homogènes et qui a été étendue aux équations aux dérivées partielles non homogènes grâce à l'utilisation des fonctions de base radiale (RBF) pour la détermination de la solution particulière (MPS). Sachant que la solution fondamentale n'est pas toujours disponible et en particulier pour les opérateurs à coefficients variables, une approximation associant ces deux méthodes (MFS-MPS) est mise au point pour suivre des branches de solutions non linéaires complexes [12]. Dans ce travail, la méthode MFS-MPS est adoptée pour la discrétisation spatiale des problèmes linéaires, résultants de la linéarisation du problème non linéaire par les méthodes asymptotiques, pour sa simplicité et sa robustesse pour la résolution des équations aux dérivées partielles à coefficients variables [8, 11,13]. Dans cet article, l'analyse est consacrée principalement à la méthode basée sur l'indicateur de bifurcation, calculé analytiquement en fonction d'un paramètre de chemin, le long d'une branche solution régie par une équation de Poisson non linéaire.

## 2. Position du problème

Pour montrer l'efficacité et faciliter la compréhension de la méthode proposée, on considère le problème non linéaire (1) dans le domaine $\Omega = [0,1] \times [0,1]$ où les points de bifurcation et les modes propres sont connus analytiquement :

$$ 
\left\{ 
\begin{array}{ll}
\Delta u + u^2 - \lambda u = \lambda \xi f & \text{dans } \Omega \\
u = 0 & \text{sur } \partial \Omega
\end{array}
\right. 
\qquad (1) 
$$

La variable $u$ et le paramètre $\lambda$ sont les inconnues du problème. La fonction $f$ et le paramètre $\xi$ sont des données du problème. Pour introduire l'indicateur de bifurcation $\mu$, la solution fondamentale $u$ est perturbée par une fluctuation $\delta u$ qui satisfait le système d'équations suivant :

$$ 
\left\{ 
\begin{array}{ll}
L_t(\delta u) = \Delta(\delta u) + 2u \delta u - \lambda \delta u = \mu f_{\mu} & \text{dans } \Omega \\
\delta u = 0 & \text{sur } \partial \Omega \\
\langle \delta u - \delta u_0, \delta u_0 \rangle = 0 & \text{dans } \Omega
\end{array}
\right. 
\qquad (2) 
$$

où $f_{\mu}$ est une fonction aléatoire donnée, $\langle \cdot, \cdot \rangle$ représente le produit scalaire euclidien et $\delta u_0$ est la réponse à la perturbation $f_{\mu}$ pour $\mu = 1$. L'indicateur $\mu$ peut être calculé directement à partir de l'équation (2) par la formule suivante :

$$ \mu = \frac{\langle \delta u_0, \delta u_0 \rangle}{\langle L_t^{-1} f_{\mu}, \delta u_0 \rangle} \qquad (3) $$

Cette méthode, dite méthode directe, permet de donner les courbes de référence.

## 3. Méthode de perturbation

La méthode de perturbation consiste à chercher les inconnues $(u, \lambda, \delta u, \mu)$ sous forme d'un développement de Taylor à partir des solutions connues $(u_0, \lambda_0, \delta u_0, 1)$. En effet,

$$ (u(a), \lambda(a), \delta u(a), \mu(a)) = (u_0, \lambda_0, \delta u_0, 1) + \sum_{i=1}^{N} a^i (u_i, \lambda_i, \delta u_i, \mu_i) \qquad (4) $$

où N représente l'ordre de troncature des séries entières et le paramètre de chemin $a'$ est défini par le produit scalaire suivant:
$$ a = \langle u - u_0, u_1 \rangle + (\lambda - \lambda_0)\lambda_1 \qquad (5) $$

Les équations (1-5) mènent à une succession de problèmes linéaires pour chaque ordre $p \ge 1$ de la série aussi bien pour la solution fondamentale $u$ que pour la fluctuation $\delta u$ :

$$
\left\{
\begin{array}{ll}
L_t^0(u_p) = \lambda_p(\xi f + u_0) - \sum_{i=1}^{p-1} u_i u_{p-i} + \sum_{i=1}^{p-1} \lambda_i u_{p-i} & \text{dans } \Omega \\
\langle u_p, u_1 \rangle + \lambda_p \lambda_1 = \eta & \text{dans } \Omega \\
u_1 = 0 & \text{sur } \partial\Omega \\
\hline
L_t^0(\delta u_p) = \mu_p f_\mu + \sum_{i=1}^{p} \lambda_i \delta u_{p-i} - \sum_{i=1}^{p} 2u_i \delta u_{p-i} & \text{dans } \Omega \\
\langle \delta u_p, \delta u_0 \rangle = 0 & \text{dans } \Omega \\
\delta u_p = 0 & \text{sur } \partial\Omega
\end{array}
\right.
\qquad (6)
$$

où $\eta = 1$ à l'ordre 1 et nul aux ordres supérieurs. L'opérateur tangent $L_t^0(\cdot) = \Delta(\cdot) + 2u_0(\cdot) - \lambda_0(\cdot)$ est commun à tous les problèmes linéaires. La méthode de continuation permet de calculer toute la branche solution pas à pas [2]. L'indicateur de bifurcation est alors évalué à chaque pas de calcul. Pour améliorer la solution, les séries entières sont remplacées par des fractions rationnelles appelées approximants de Padé [4].

## 4. Méthode MFS-MPS

Dans cette section, nous proposons de résoudre par MFS-MPS les problèmes linéaires (6) pris sous la forme simple suivante à chaque pas,

$$
\left\{
\begin{array}{ll}
\Delta(u) + \gamma(x,y)u = f(x,y) & \text{dans } \Omega \\
u = g(x,y) & \text{sur } \partial\Omega
\end{array}
\right.
\qquad (7)
$$

Ce sont des EDP avec des opérateurs différentiels à coefficients variables et dont la solution fondamentale n'est pas disponible. Pour pallier cette difficulté, la méthode des équations analogues (AEM) [12] est utilisée pour transformer le problème (7) en un problème de Poisson standard. Ceci consiste à maintenir le Laplacien, dont la solution fondamentale est connue, à gauche et les autres termes sont déplacés à droite et sont considérés comme un second membre. L'équation (7) devient alors,

$$
\left\{
\begin{array}{ll}
\Delta(u) = -\gamma(x,y)u + f(x,y) = h(x,y,u) & \text{dans } \Omega \\
u = g(x,y) & \text{sur } \partial\Omega
\end{array}
\right.
\qquad (8)
$$

qui est facilement résolue par la méthode MFS si le terme h est connu. Or h dépend de l'inconnue u, d'où la nécessité de la méthode MFS-MPS.

Ici, nous présentons brièvement la technique MFS-MPS nouvellement développée. L'idée suppose que la solution d'une EDP combine simultanément la solution particulière et la solution homogène. Ceci permet d'écrire cette solution sous la forme suivante:

$$ u(x, y) \approx \tilde{u}(x, y) = \sum_{i=1}^{N_t} \alpha_i v(r_i) + \sum_{i=1}^{N_b} \beta_i G(\rho_i) \qquad (9) $$

où $N_t = N_d + N_b$, $N_d$ et $N_b$ représentent respectivement les points d'interpolation à l'intérieur du domaine $\Omega$ et le nombre de points sur sa frontière. $\{\alpha_j\}_{j=1}^{N_t}$ et $\{\beta_j\}_{j=1}^{N_b}$ sont des coefficients inconnus à déterminer. Les fonctions $v(r_i)$ représentent les solutions particulières de l'opérateur de Laplace défini par,

$$ \Delta(v(r_i)) = \phi(r_i) \qquad (10) $$

où les $\phi(r_i)$ sont des fonctions de base radiale qui permettent d'approcher le second terme de la première équation de (8), ainsi:
$$h(x, y, u) = \sum_{j=1}^{N_t} \alpha_j \phi(r_j) \qquad (11)$$

où $r_j = \|(x, y) - (x_j, y_j)\|$, $(x_j, y_j)_{j=1..N_t}$ sont des points de collocation dans le domaine. On note également que,

$$G(\rho_i) = -\frac{1}{2\pi}\log(\rho_i) \qquad (12)$$

est la solution fondamentale du Laplacien pour le problème 2D caractérisée par,

$$\Delta G(\rho) = 0 \qquad (13)$$

Les $\rho_j = \|(x, y) - (x_s_j, y_s_j)\|$, $(x_s_j, y_s_j)_{j=1..N_b}$ représentent les points sources pris sur une frontière fictive $\Gamma$ contenant le domaine $\Omega$. Au lieu de chercher la solution particulière et la solution homogène séparément, nous calculons les coefficients $\{\alpha\}_{j=1}^{N_t}$ et $\{\beta\}_{j=1}^{N_b}$ simultanément. En utilisant l'approximation (10) et en introduisant (9) dans (8), on a:

$$\sum_{i=1}^{N_t} \alpha_i \Psi(r_i) + \sum_{i=1}^{N_b} \beta_i \Theta(\rho_i) = f(x, y) \quad \text{dans } \Omega \qquad (14)$$

Telle que

$$\Psi(r_i) = \Delta v(r_i) + \gamma(x, y)v(r_i) = \phi(r_i) + \gamma(x, y)v(r_i) \qquad (15)$$

$$\Theta(\rho_i) = \Delta G(\rho_i) + \gamma(x, y)G(\rho_i) = \gamma(x, y)G(\rho_i) \qquad (16)$$

La condition aux limites (12b) devient alors,

$$\sum_{i=1}^{N_t} \alpha_i v(r_i) + \sum_{i=1}^{N_b} \beta_i G(\rho_i) = g(x, y) \quad \text{sur } \partial\Omega \qquad (17)$$

Il est clair que $\Psi(r)$ et $\Theta(\rho)$ sont des fonctions complètement connues.

L'efficacité et la précision de l'interpolation dépend du choix des fonctions de base radiale $\phi$. En effet, les polynômes de Lagrange, les interpolants de Fourier, les fonctions radiales polynomiales, les thin plate Splines (TPS) ainsi que les fonctions multi-quadriques (MQ) peuvent être utilisés [9, 10]. Dans ce papier, le choix s'est restreint aux fonctions de base multi-quadriques où la méthode MFS-MPS peut coïncider avec la méthode de Kansas, une des méthodes les plus connues dans la communauté Meshless. On prend alors,

$$\phi(r) = \sqrt{c^2 + r^2} \qquad (18)$$

d'où

$$v(r) = \frac{1}{9}(4c^2 + r^2)\sqrt{c^2 + r^2} - \frac{c^3}{3}\ln(c + \sqrt{c^2 + r^2}) \qquad (19)$$

en 2D. Le paramètre c est une constante choisie.

Pour la mise en œuvre de cette technique, nous choisissons $N_t$ points intérieurs dans le domaine $\Omega$, $N_b$ points sur la frontière $\partial\Omega$ et $N_b$ points source sur la frontière fictive $\Gamma$. Par collocation, on obtient un système d'ordre $(N_t + N_b) \times (N_t + N_b)$ à résoudre:

$$\begin{bmatrix} \Psi(r_{ij}) & \Theta(\rho_{ij}) \\ v(r_{ij}) & G(\rho_{ij}) \end{bmatrix} \begin{Bmatrix} \alpha_j \\ \beta_j \end{Bmatrix} = \begin{Bmatrix} f \\ g \end{Bmatrix} \qquad (20)$$

Une fois les coefficients $\{\alpha_j\}_{j=1}^{N_t}$ et $\{\beta_j\}_{j=1}^{N_b}$ sont déterminés par (20), la solution approchée est obtenue à partir de (9).
# 5. Applications numériques

## 5.1 Exemple 1

Dans cet exemple, 9 points sont choisis à l'intérieur du domaine, 16 points sur la frontière $\partial\Omega$ et le même nombre de points sources distribués uniformément sur une frontière fictive $\Gamma = [-4,5] \times [-4,5]$, (Fig. 1).

<figure><img src="image_2.png" /><figcaption>Fig. 1: Discrétisation du domaine et de la frontière fictive</figcaption></figure>

Les points de bifurcation $\lambda_{bif}$ et les modes propres $u_{bif}$ sont connus analytiquement et sont donnés par la formule suivante :

$$ \begin{cases} \lambda_{bif} = (m^2 + n^2)\pi^2 \\ u_{bif} = \sin(m\pi x) \sin(n\pi y) \end{cases} \qquad (21) $$

où m et n sont des entiers. Ici, on s'intéresse principalement au premier point de bifurcation avec $m = 1$ et $n = 1$. Le terme $\lambda\xi f$ de l'équation (1) représente un défaut permettant de suivre automatiquement les branches bifurquées. Pour s'approcher du point de bifurcation, le paramètre $\xi$ est choisi égal à $\pm 10^{-6}$. Tenant compte des travaux [10], le paramètre c est pris égal à 1.5

Pour illustrer l'efficacité de notre algorithme, la variable u au point $(x = 0.5, y = 0.5)$ est reportée sur la figure 2 en fonction du paramètre $\lambda$.

<figure><img src="image_3.png" /><figcaption>Fig. 2: Branches de solutions par MAN-MFS-MPS pour $\xi = 10^{-6}$ et $\xi = -10^{-6}$</figcaption></figure>

Sur la figure 2, une accumulation de pas est observée pour $\lambda \approx 19.79$ (chaque pas de continuation est représenté par un symbole circulaire ou un carré sur la courbe). Après cette accumulation de pas, la solution suit la branche bifurquée au lieu de rester sur la branche fondamentale en localisant géométriquement le premier point de bifurcation. Sur la figure 3, on reporte l'évolution de l'indicateur $\mu$ en fonction du paramètre $\lambda$ en utilisant des séries entières tronquées aux ordres 10 et 30, des approximants de Padé à l'ordre 20 et la méthode directe considérée comme courbe de référence. Cet indicateur s'annule pour $\lambda = 19.79$. Des trois méthodes résulte la même valeur critique qui
coïncide exactement avec la valeur critique analytique. Ces résultats confirment la robustesse de l’algorithme proposé pour la détection des points de bifurcation. Des résultats similaires ont été obtenus en appliquant cet algorithme à d’autres exemples et pour différents domaines.

<figure><img src="image_4.png" /><figcaption>Fig. 3: Indicateur μ en fonction de λ par MAN-MFS-MPS</figcaption></figure>

## 5.2 Exemple 2

Pour étendre la méthode à d'autres formes géométriques, on considère le problème suivant dans un domaine circulaire :

$$ \begin{cases} \Delta u + u^2 = \lambda f \\ u = \lambda g \end{cases} \quad \begin{array}{l} \text{dans } \Omega \\ \text{sur } \partial \Omega \end{array} \qquad (22) $$

où $\Omega = \{x^2 + y^2 - 1 \le 0\}$, $f = 4 + (x^2 + y^2)^2$ et $g = x^2 + y^2$. Dans cet exemple, 69 points sont choisis à l'intérieur du domaine, 10 points sur la frontière $\partial\Omega$ et le même nombre de points sources distribués uniformément sur un cercle centré à l'origine et de rayon $R = 8$ (Fig. 4).

<figure><img src="image_5.png" /><figcaption>Fig.4 : Discrétisation du domaine et de la frontière</figcaption></figure>

On reporte, sur la figure 5, la solution $u$ de (22) au point $(x = -0.4286, y = 0.6429)$ en fonction du paramètre $\lambda$ calculée par la méthode de continuation MAN-MFS-MPS. On observe une accumulation de pas au voisinage de $\lambda \approx 2.88$ permettant de localiser géométriquement le point de bifurcation. Nous avons également calculé l'indicateur de bifurcation à l'aide des séries entières et des approximants de Padé. Les résultats sont reportés sur la figure 6 et montrent que l'indicateur s'annule effectivement pour la même valeur $\lambda \approx 2.88$. Ce résultat est confirmé en déterminant la première racine réelle du dénominateur du Padé. Ces résultats numériques témoignent de l'efficacité des algorithmes proposés dans ce travail.
<figure><img src="image_6.png" /><figcaption>Fig. 5: Branches de solutions par MAN-MFS-MPS</figcaption></figure>

<figure><img src="image_7.png" /><figcaption>Fig. 6: Indicateur $\mu$ en fonction de $\lambda$ par MAN-MFS-MPS</figcaption></figure>

# 6. CONCLUSION

Nous avons présenté dans ce travail une nouvelle technique associant une méthode sans maillage et la méthode asymptotique numérique pour la localisation des points de bifurcation sur une branche solution d'un problème de Poisson non linéaire en introduisant des indicateurs de bifurcation bien adaptés aux méthodes asymptotiques. Ces points de bifurcation ont été cherchés par trois méthodes: l'indicateur géométrique, les racines du dénominateur des approximants de Padé et l'indicateur de bifurcation analytique. Les exemples numériques montrent l'efficacité de la méthode proposée. Cette technique peut être étendue à la recherche de bifurcations multiples ou des bifurcations de Hopf.

# REFERENCES

[1] N. Damil, M. Potier-Ferry, A new method to computed bifurcations: application to the buckling of imperfect elastic structures. International Journal of Engineering Sciences, 28 934-957. (1990)

[2] B. Cochelin, A path-following technique via an asymptotic-numerical method, Computers and Structures 53 (5) 1181-1192, (1994)

[3] E.H. Boutyour, Méthode Asymptotique Numérique pour le calcul de bifurcations: application aux structures élastiques, Ph.D. Thesis, Université de Metz, France (1994).

[4] E.H Boutyour., H. Zahrouni, M. Potier-Ferry, Boudi M., Bifurcation points and bifurcated branches by an asymptotic numerical method and Padé approximants, International Journal of Numerical Methods in Engineering (60) 1987-2012, 2004

[5] A. Tri, B. Cochelin, M. Potier-Ferry, Résolution des équations de Navier Stokes et détection des bifurcations stationnaires par une méthode asymptotique numérique, Revue européenne des éléments finis, (5) 415-442, 1996
[6] S. Baguet, B. Cochelin, On the behaviour of the ANM continuation in the presence of bifurcations, Communication in Numerical Methods in Engineering. (19) 459–471, 2003

[7] A. Tri, H. Zahrouni, M. Potier-Ferry, High order continuation algorithm and meshless method to solve nonlinear Poisson problems, Engineering Analysis with Boundary Elements (35) 273–278, 2011

[8] A. Tri, H. Zahrouni, M. Potier-Ferry, Perturbation technique and method of fundamental solution to solve nonlinear Poisson problems, Engineering Analysis with Boundary Elements (36) 1705–1714, 2012

[9] M.A. Golberg, The method of fundamental solutions for Poisson's equation, Engineering Analysis with Boundary Elements, Volume 16 (3), 205-213, 1995

[10] M.A. Golberg, C.S Chen, S.R Karur, Improve multiquadric interpolation for partial differential equations, Engineering Analysis with Boundary Elements, (18) 9-17, 1996

[11] G.Yao, C.H. Tsai, W. Chen, The comparison of three meshless methods us in gradial basis functions for solving fourth-order partial differential equations Engineering Analysis with Boundary Elements (34) 625–631, 2010

[12] H. Wang, Q.H. Qin, A meshless method for generalized Linear or nonlinear Poisson-type problems Engineering Analysis with Boundary Elements (30) 515–521, 2006

[13] C.S. Chen, C.M. Fan, J. Monroe, The method of fundamental solutions for solving elliptic pdes with variable coefficients, In : Chen CS, Karageorghis A, Smyrlis YS, editors. The method of fundamental solutions-meshless method. Atlanta : Dynamic Publishers 75–105, 2008