r7.05.02 Analyse de stabilité statique des pentes en remblai

Résumé:

La commande CALC_STAB_PENTE permet d’analyser la stabilité statique des pentes en remblai modélisées en 2D dans les ouvrages hydrauliques tels que les barrages, les talus et les digues autour des réservoirs, en prenant compte la pression interstitielle dans le milieu poreux et éventuellement d’autres types de chargement mécanique.

En général, pour analyser la stabilité statique des pentes selon les diverses réglementations du métier dans le monde, on utilise la méthode de réduction progressive des propriétés matériaux (méthode SRM – strengh reduction method ), ou bien la méthode basée sur la théorie d’équilibre limite (méthode LEM). Les deux sont disponibles dans CALC_STAB_PENTE, de sorte que la commande s’adapte à toutes les circonstances d’analyse en condition statique.

Dans ce document on présente en détails les algorithmes implémentés dans [CALC_STAB_PENTE].

Méthode SRM

Généralité de la méthode SRM

La méthode SRM met à profit la précision de la modélisation par la méthode aux éléments finis pour analyser la stabilité des pentes sans avoir besoin de supposer par avance les cercles de glissement. En effet, en utilisant un critère de rupture de type Mohr-Coulomb ou Drucker-Prager, la méthode SRM définit le facteur de sécurité comme étant le rapport entre les paramètres matériaux initiaux et les paramètres réduits conduisant à la divergence du calcul non-linéaire sous les chargements statiques appliqués :

(4616)\[\mathit{FS} = \frac{\cohesion}{{\cohesion}_{i}} = \frac{\tan(\anglFric )}{\tan({\anglFric }_{i})} = \frac{\tan(\anglDila )}{\tan({\anglDila }_{i})}\]

où

• \(\anglDila\) est l’angle de dilatance du matériau

• \(\cohesion\) et \(\anglFric\) représentent respectivement la cohésion et l’angle de frottement du matériau composant le remblai.

Griffiths et al. [bib1] montrent que l’expression (5000) a le même sens physique que la définition traditionnelle par le rapport entre la force de résistance et la force tangentielle le long de la surface de rupture.

Comme l’illustre la Fig. 457 , le traitement informatique dans la commande [CALC_STAB_PENTE] de la méthode SRM se fait en trois étapes:

• Étape 1: Traitement de la zone SRM. On identifie les matériaux existant dans la zone SRM indiquée, et on crée les groupes des mailles dans la zone SRM afin de réaffecter le champ des matériaux dégradés.

• Étape 2: Création du champ des matériaux dégradé. On teste ensuite la convergence du calcul non-linéaire et on raffine l’incrément du facteur de sécurité jusqu’à ce que le calcul diverge et que la précision souhaitée soit atteinte.

• Étape 3: Post-traitement. Si demandé, on obtient en sortie le champ de déformation plastique cumulée (V1 du champ VARI_NOEU) afin de visualiser la surface de rupture.

Fig. 457 Diagramme de la méthode SRM

Traitement de la zone SRM

La commande [CALC_STAB_PENTE] permet de définir la zone où l’algorithme SRM s’applique. Cette fonctionnalité est particulièrement adaptée lorsque le modèle comprend plusieurs pentes et qu’on souhaite vérifier la stabilité d’une pente qui ne subit pas obligatoirement le risque de rupture le plus important. Il est à noter que la zone SRM peut être définie par les groupes des mailles distincts de ceux qui sont affectés dans le champ des matériaux en entrée. Dans ce cas-là on identifie les matériaux dans la zone SRM et on réaffecte les groupes des mailles par le principe de surcharge.

Comme l’illustre la Fig. 458 , pour chaque groupe des mailles dans la zone SRM, on examine s’il possède des mailles en commun avec les groupes des mailles affectés par les matériaux présents dans l’objet aster CHAM_MATER.

Fig. 458 Diagramme de la génération du cham_mater dégradé dans la zone SRM

Réduction des propriétés de résistance des matériaux

Les critères de rupture Mohr-Coulomb et Drucker-Prager sont autorisées dans [CALC_STAB_PENTE]. En même temps que la génération du champ des matériaux dégradés, on vérifie qu’il existe exactement une loi de comportement parmi les deux affectée à la zone SRM par le mot-clé facteur COMPORTEMENT.

D’après l’expression (5000), les propriétés dégradées sont calculées par:

(4617)\[\tan({\anglFric}^{\prime}) = \frac{\tan(\anglFric)}{\mathit{FS}} \text{ et } {\cohesion^{\prime}} = \frac{\cohesion}{\mathit{FS}}\]

On fait l’hypothèse de la loi d’écoulement plastique associée du comportement Mohr-Coulomb, soit \({\anglFric}^{\prime}={\anglDila}^{\prime}\). Dans le cas contraire, l’angle de dilatance prendra automatiquement la valeur de l’angle de frottement lors de l’exécution de l’algorithme SRM.

A part le critère de Mohr-Coulomb, il est également possible d’employer le critère de Drucker-Prager qui améliore le problème de singularité causé par les sommets du cône de la surface de charge de Mohr-Coulomb dans le plan \(\pi\) . Au lieu de réduire directement les paramètres de Drucker-Prager (\({\alpha}\) et \(\yieldStress\) ), la commande privilégie la stratégie phi-c pour éviter le conservatisme du résultat [bib3] :

(4618)\[\tan {\anglFric}^{\prime} = \frac{3{\alpha}}{2\mathit{FS}\sqrt{(2{\alpha}+1)(1-{\alpha})}} \text{ et } {\cohesion}^{\prime} = \frac{\yieldStress}{2\mathit{FS}\sqrt{(2{\alpha}+1)(1-{\alpha})}} \text{ et } \tan{ {\anglDila }^{\prime}} = \frac{\tan{\anglDila }}{\mathit{FS}}\]

Contrôle de la convergence du FS

Afin d’accélérer la recherche du facteur de sécurité (FS), il est possible de contrôler la résolution itérative en configurant les paramètres liés au raffinement de l’incrément \({dP}\) (Fig. 457) selon la méthode indiquée. L’incrément diminue à chaque fois que le calcul non-linéaire diverge jusqu’à ce que le résidu maximal soit atteint.

Deux lois de variation de l’incrément \({dP}\) sont disponibles :

1. “EXPONENTIELLE” : Variation exponentielle

L’incrément du facteur de sécurité approche le résidu maximal de manière exponentielle. Au n-ième raffinement, l’incrément est calculé par:

(4619)\[{p}_{n} = {p}_{0} \, {2}^{-n}\]

jusqu’à ce que \({p}_{N}<{p}_{\mathit{fin}}\) , où

• \({p}_{n}\) est l’incrément courant

• \({p}_{0}\) est l’incrément initial du facteur de sécurité.

1. “LINEAIRE” : Variation linéaire

L’incrément du facteur de sécurité approche le résidu maximal de manière linéaire. Dans ce cas-là un opérande supplémentaire ITER_RAFF_LINE est nécessaire. Au n^e raffinement, l’incrément est calculé par :

(4620)\[{p}_{n} = {p}_{0}-n \, \frac{{p}_{0}-{p}_{\mathit{fin}}}{N}\]

jusqu’à ce que \({p}_{N}={p}_{\mathit{fin}}\).

Note

L’avantage de la méthode exponentielle est que le calcul converge plus rapidement, alors que la méthode linéaire contrôle mieux la précision du résultat. La variation exponentielle de l’incrément s’adapte surtout au calcul avec grand INCR_INIT.

La méthode suivante simplifierait le choix entre les deux loi de variation de l’incrément. Soit \(N\) le nombre de raffinements souhaité, on choisit la loi linéaire si :

(4621)\[\frac{{p}_{0}-{p}_{\mathit{fin}}}{{{log}}_{2}({p}_{0})-{{log}}_{2}({p}_{\mathit{fin}})}<N\]

Dans ce cas, c’est la loi exponentielle qui est choisie.

Méthode LEM

Généralité de la méthode LEM

Toutes les méthodes type LEM se basent sur la définition du facteur de sécurité comme étant le rapport entre la force de résistance au cisaillement et la force tangentielle le long de la surface de rupture.

Par division de la pente en tranches, qui sont considérées comme des solides indéformables, on résout le bilan des forces et des moments de chaque tranche pour en déduire le facteur de sécurité. En respectant ce principe de calcul, les diverses méthodes LEM se distinguent par le nombre des bilans satisfaits et les hypothèses sur la force d’interaction entre les tranches.

Un panorama général des méthodes de résolution par équilibre limite est donné dans la Fig. 459. 4 méthodes LEM sont disponibles dans [CALC_STAB_PENTE], adaptées à l’analyse des pentes avec surface de rupture tant circulaire que non-circulaire:

• Méthode Bishop simplifiée et méthode Fellenius : Ce sont les méthodes les plus classiques qui supposent que la surface de rupture est circulaire.

  Grâce au glissement rotationnel, on réduit le nombre d’équations d’équilibre en moment, qui est égal au nombre des tranches, en une équation décrivant l’équilibre en moment de la partie glissante entière. Ainsi, on calcule le facteur de sécurité de manière simple par une expression explicite (Fellenius) ou par une équation scalaire implicite (Bishop).

  Avertissement

  Les hypothèses du modèle étant trop simplifiées, surtout pour la méthode Fellenius, le résultat peut être éloigné des observations physiques.

• Méthode Spencer et Morgenstern-Price : Ce sont les méthodes les plus classiques supposant que la surface de rupture est non-circulaire.

  On résout la totalité des équations d’équilibre, dont la quantité est égale à 3 fois le nombre des tranches, afin d’obtenir le résultat de facteur de sécurité le plus précis possible. La procédure de résolution devient plus sophistiquée, ce qui permet d’obtenir un résultat plus juste si la surface de rupture est visuellement non-circulaire.

  Avertissement

  Cependant, il est à noter que le résultat ne différera pas significativement de celui obtenu par la méthode de Bishop simplifiée si le champ des matériaux est homogène et isotrope.

---

En outre, les méthodes LEM ne s’appuyant pas sur un calcul aux éléments finis, les algorithmes d’optimisation implémentés dans [CALC_STAB_PENTE] sont indispensables pour repérer la surface de rupture qui minimise globalement le facteur de sécurité :

• Méthode d’exhaustion pour la surface circulaire : On parcourt toutes les combinaisons des points d’extrémité des potentielles surfaces de rupture. Il s’agit de la méthode implémentée dans la plupart des logiciels géotechniques commerciaux tels que SLOPE/W de GeoStudio [bib8].

• Algorithme de feu d’artifice amélioré (EFWA) pour la surface non-circulaire [bib4] : Il s’agit d’un algorithme d’essaims intelligents qui s’inspire du processus de la génération des étincelles lors des explosions des feux d’artifice. Il s’avère plus performant que d’autres algorithmes d’optimisation pour la résolution du problème de stabilité des pentes avec surface de rupture non-circulaire [bib5].

Fig. 459 Panorama des méthodes LEM

Une des particularités de l’implémentation des méthodes LEM est que l’on calcule le poids des tranches à l’aide du même maillage utilisé pour le calcul aux éléments finis. Bien que cela impose des optimisations de maillage, développés dans la commande (de manière à éviter au mieux de ralentir la résolution du facteur de sécurité), l’avantage se traduit par le fait qu’on n’aura plus besoin de basculer entre les solveurs/modèles pour effectuer les différents types d’analyse d’un ouvrage en remblai.

Par exemple il est possible d’enchaîner directement le calcul du champ de pression interstitielle due à l’infiltration d’eau (voir le modèle THM [u2.04.05]) et l’analyse de stabilité des pentes. Ce dernier peut être suivi par la calcul du déplacement irréversible muni d’un résultat d’analyse transitoire dynamique (voir la commande [POST_NEWMARK] ).

La Fig. 460 montre l’algorithme global de la méthode LEM implémentée.

Fig. 460 Diagramme d’implémentation des méthodes LEM dans [CALC_STAB_PENTE]

Calcul des forces et de la géométrie des tranches

Afin de résoudre le bilan des forces et des moments des tranches, il faut dans un premier temps calculer les paramètres géométriques définissant la partie glissante du sol, intégrer les chargements mécaniques extérieurs appliqués sur les tranches, calculer leur poids propre et éventuellement la pression interstitielle à la base des tranches. La Fig. 461 montre le flux de travail global qui évalue les métriques nécessaires au calcul du facteur de sécurité par la suite.

Fig. 461 Flux de travail de la récupération des informations des tranches pour calculer FS

Calcul du centre du cercle de glissement

En cas de surface de rupture circulaire, afin de créer le groupe des mailles de la partie glissante, on retrouve les coordonnées du centre du cercle de glissement de deux manières :

1. Avec comme entrées les points d’extrémités et le rayon :

   En injectant les coordonnées des points d’extrémités (\(({x}_{1},{y}_{1})\) et \(({x}_{2},{y}_{2})\) ) dans l’équation du cercle, on obtient les équations exprimant la relation entre les coordonnées du centre \(({x}_{0},{y}_{0})\) et une équation d’ordre deux de \({y}_{0}\) :

   (4622)\[\begin{split}\left\lbrace \begin{array}{l} x_{0} = c_{1} - c_{2} \, y_{0} \\ (c_{2}^{2} + 1) y_{0}^{2} + (2 \, x_{1} c_{2} - 2 \, c_{1} c_{2} - 2 \, y_{1}) y_{0} + (x_{1} - c_{1})^{2} + y_{1}^{2} - R^{2} = 0 \end{array} \right.\end{split}\]

   où les constantes \({c}_{1}\) et \({c}_{2}\) s’écrivent :

   (4623)\[\begin{split}\begin{array}{c} {c}_{1} = \frac{{x}_{2}^{2}-{x}_{1}^{2}+{y}_{2}^{2}-{y}_{1}^{2}}{2({x}_{2}-{x}_{1})}\\ {c}_{2} = \frac{{y}_{2}-{y}_{1}}{{x}_{2}-{x}_{1}} \end{array}\end{split}\]

   Pour résoudre l’équation (5006), on exprime les coefficients dans la deuxième équation comme :

   (4624)\[\begin{split}\begin{array}{c} A = {c}_{2}^{2}+1\\ B = 2 \, ({x}_{1}-{c}_{1}){c}_{2}-2{y}_{1}\\ C = ({x}_{1}-{c}_{1}^{2}+{y}_{1}^{2}-{R}^{2}) \end{array}\end{split}\]

   Ainsi on obtient l’ordonnée du centre. La surface de rupture étant toujours convexe, on prend la solution la plus grande :

   (4625)\[{y}_{0} = \frac{-B+\sqrt{{B}^{2}-4 {AC}}}{2A}\]

2. Avec comme entrées les points d’extrémités et une droite horizontale / verticale tangentielle :

   Dans le cas où on contrôle manuellement la recherche du rayon critique par les opérandes Y_MINI ou Y_MAXI de [CALC_STAB_PENTE], on en déduit le centre du cercle à partir de l’ordonnée \({y}_{b}\) de la droite horizontale tangentielle. Il faut donc remplacer \(R={y}_{0}-{y}_{b}\) dans l’équation (5006) et reprendre la même procédure de résolution.

   De même, on remplace \(R^2\) par \(R^2 = (x_0 - x_b)^2\) dans l’équation (5006) lors de l’évaluation de la limite inférieure du rayon, soit le rayon du cercle tangentiel à la droite verticale \(x = x_b\) où \(x_b\) est l’abscisse de l’extrémité active.

Création des tranches et calcul des poids propres

On crée le groupe des mailles des tranches à l’aide de l’opérateur [DEFI_GROUP] en suivant les étapes suivantes :

1. Étape 1: Création du groupe des mailles de la partie glissante avec les mailles dont les noeuds sont strictement inclus dans la masse glissante. On crée un maillage réduit à partir de ce groupe des mailles en vue d’économiser du temps dans les étapes suivantes.

2. Étape 2: Création du groupe des mailles de la tranche \(n\) par l’option INTERSEC qui permet de récupérer les mailles en commun de la bande d’abscisse de la tranche et la partie glissante.

3. Étape 3 : Correction de la tranche \(n\) par l’option DIFFE qui permet d’exclure les mailles appartenant à la tranche précédente. Cette étape assure que la masse totale des tranches est égale à la masse de la partie glissante.

Un exemple des groupes de mailles est montré dans la Fig. 462. Les poids propres des tranches sont ensuite calculés via l’opérateur [POST_ELEM].

On calcule en même temps l’angle d’inclinaison \(\alpha\) de la base des tranches telle que :

(4626)\[\tan \alpha = \frac{{y}_{n+1}-{y}_{n}}{{x}_{n+1}-{x}_{n}}\]

Fig. 462 Exemple de groupe des mailles des tranches et les centres à la base

Intégration des chargements mécaniques extérieurs

Intégration de la pression répartie

La pression répartie appliquée sur la surface de la pente est introduite via le mot-clef FONC_PRES sous la forme de fonction ou formule avec comme variables les coordonnées X et/ou Y. Comme le montre la Fig. 463, une série d’opérations s’effectuent de manière à évaluer la force nodale induite par la pression répartie. On profite de ce champ nodal pour déterminer le vecteur résultant des efforts et de moments, par tranche, en intégrant la pression le long des segments 1D :

(4627)\[\vector{P}_i = \int_{l_{i-1}}^{l_i} { - \pres(l) \, \normal(l) \measLine }\]

(4628)\[\vector{M}_{Q, i} = \int_{l_{i-1}}^{l_i} { [\vector{r}_c - \vector{r}(l)] \times \pres(l) \, \normal(l) \measLine }\]

où \({\normal}\) représente le vecteur normal unitaire des éléments 1D, \(\vector{r}\) est le vecteur coordonnées des noeuds sur la pente, \(\vector{r}_c\) sont les coordonnées du point par rapport auquel le moment est défini. Pour les méthodes Bishop et Fellenius, \(\vector{r}_c\) est le centre du rayon. Pour les méthodes Morgenstern et Spencer, il s’agit du point au milieu de la base des tranches. On adopte l’intégration trapézoïdale pour calculer numériquement (5009) et (5010).

Fig. 463 Calcul de la force nodale sur la pente induite par la pression répartie.

Intégration de la force inertielle sismique

Il est possible d’appliquer la force inertielle volumique due au séisme via le mot-clé facteur ACCE. Sachant le poids propre des tranches, on détermine la force et le moment inertielle par tranche telle que :

(4629)\[\vector{F}_{i} = k_a W_i \frac{\vector{r}}{\Vert \vector{r} \Vert_2}\]

(4630)\[\vector{M}_{S, i} = (\vector{r}_{b, i} - \vector{r}_c) \times \vector{F}_{i}\]

où \(\vector{r}_{b, i}\) représente les coordonnées du barycentre des tranches, \(k_a\) est le coefficient d’accélération, \(W_i\) est le poids des tranches, \(\vector{r}\) est la direction d’excitation sismique.

Propriétés des matériaux et pression interstitielle

Récupération des propriétés de résistance

Le champ des matériaux en entrée de la commande [CALC_STAB_PENTE] permet de retrouver la cohésion et l’angle de frottement le long de la surface de rupture. Les paramètres du comportement Drucker-Prager sont convertis par la même formule que dans la méthode SRM avec \(\mathit{FS}=1\).

Sans l’opérande PHI_C_EQUI, on suppose que les propriétés sont homogènes sur la base de chaque tranche et égales à ceux du matériau affectant la maille qui contient le centre de la base. Cette maille «centre» est définie par le critère de distance minimum entre le centre de la maille et le point milieu de la base des tranches.

Cependant, pour les structures non homogènes ou stratifiées, il existe probablement des tranches dont la base traverse plusieurs matériaux. Dans ce cas-là, en activant PHI_C_EQUI, le programme calcule la cohésion et l’angle de frottement équivalentes sous l’hypothèse que la contrainte normale distribue uniformément sur la base. Pour ce faire, on identifie d’abord les points d’intersection entre la base des tranches et les interfaces des matériaux. Par une approche dichotomique, on divise le segment de la base en deux de manière récursive jusqu’à ce que les matériaux aux extrémités soient identiques ou que la longueur du segment soit inférieur à la valeur définie dans l’opérande TOLE_DIST. Les propriétés équivalentes sont déterminées telles que :

(4631)\[\cohesion_{eq} = \frac{\sum_{k=1}^{N_c} \cohesion_k \, l_k}{L}, \quad \tan \anglFric_{eq} = \frac{\sum_{k=1}^{N_c} \tan \anglFric_k \, l_k} {L}\]

où \(N_c\) est le nombre des matériaux sur la base des tranches, \(l_k\) est la longueur de base dans le matériau \(k\), \(L\) est la longueur totale de la base.

Le calcul des propriétés équivalentes a pour effet d’éliminer la discontinuité du FS/KC dans la zone de recherche, ce qui accélère la localisation du minima global et rend le résultat plus pertinent dans certains cas spécifiques. Comme le montre la Fig. 464, avec PHI_C_EQUI on obtient une courbe de variation du coefficient d’accélération critique plus lisse sans oscillation locale lors d’une analyse de stabilité par la méthode Bishop.

Fig. 464 Comparaison des résultats avec et sans PHI_C_EQUI (cf. [V6.03.011]).

Interpolation du champ de pression hydraulique

Il est possible de définir le champ de pression hydraulique soit par un simple coefficient de pression interstitielle \(R_u\), soit par la ligne phréatique, soit par un résultat evol_noli contenant les composants PRE1 et PRE2. Les deux premières approches permettent de reconstruire le champ de pression interstitielle de manière analytique [CALC_STAB_PENTE]. Dans le mot-clé CHAM_PRES, deux méthodes d’interpolation sont disponibles au choix de l’utilisateur, qui servent à évaluer la pression hydraulique au centre de la base des tranches.

• Interpolation pondérée par l’inverse de la distance (IDW)

  Soit \(u_i\) la pression capillaire ou la pression d’air sur les noeuds de la maille qui contient le centre de la base des tranches, l’interpolation se fait par :

  (4632)\[u = \frac{\sum_{i=1}^{N}{(\frac{1}{{d}_{i}})}^{2}{u}_{i}}{\sum_{i=1}^{N}{(\frac{1}{{d}_{i}})}^{2}}\]

  où \({d}_{i}\) désigne la distance entre les nœuds de la maille et le centre de la base.

  Il s’agit d’une méthode simple et efficace. Cependant, la pondération des noeuds distants diminuant rapidement en fonction de la distance, seuls les noeuds au voisinage immédiat influencent le résultat d’interpolation. Donc cette méthode est plutôt adapté sur un maillage fin.

• Interpolation à l’aide de la fonction de forme

  On utilise la commande [PROJ_CHAMP] pour projeter le champ nodal des pressions hydrauliques sur un maillage constitué des éléments discrets seuls dont les coordonnées sont les centres de base des tranches. Comparée avec l’IDW, cette méthode est plus fiable quand le maillage est grossier.

Prise en compte de la succion matricielle

Il est bien connu que la résistance au cisaillement des sols augmente avec la succion matricielle, qui se traduit par la pression capillaire positive (pression d’eau négative) dans le sol non saturé au-dessus de la ligne phréatique. Il existe deux options pour calculer la résistance avec succion matricielle :

Option 1

On définit un angle constant \(\anglFric_b\) traduisant l’augmentation de la résistance en fonction de la succion. On détermine la résistance telle que :

(4633)\[S = {\cohesion}^{\prime} + (\stressCmp_n - \presAir) \tan \anglFric^{\prime} + \presCapi \tan \anglFric_b\]

où \({\cohesion}^{\prime}\) et \({\anglFric}^{\prime}\) représentent la cohésion et l’angle de frottement effectifs, \(\stressCmp_n\) est la contrainte normale à la base des tranches, \(\presAir\) est la pression d’air, \(\presCapi = \presAir - \presEau\) est la pression capillaire.

Option 2

On calcule la résistance selon le comportement hydraulique renseigné par l’utilisateur dans DEFI_MATERIAU/THM_DIFFU. Comme proposée par Vanapalli et al[9], la résistance s’exprime par :

(4634)\[S = {\cohesion}^{\prime} + (\stressCmp_n - \presAir) \anglFric^{\prime} + \presCapi \left(\frac{\theta_w - \theta_r}{\theta_s - \theta_r}\right) \tan {\anglFric}^{\prime}\]

où \(\theta_w\), \(\theta_s\), \(\theta_r\) représente respectivement la saturation en eau, la saturation à l’état saturé et la saturation résiduelle. Selon les propriétés renseignées dans le champ des matériaux en entrée, on détermine la saturation effective qui s’exprime par :

(4635)\[S_{\textrm{we}} = \frac{\theta_w - \theta_r}{\theta_s - \theta_r}\]

soit par les paramètres du modèle de Van-Genuchten, soit par la fonction \(\theta_{w}(\presCapi)\) définie via le mot-clef SATU_PRES dans DEFI_MATERIAU/THM_DIFFU. Dans le deuxième cas, \(\theta_s\) et \(\theta_r\) prennent la valeur du premier et du dernier élément dans la liste définissant la fonction.

Si les paramètres de Van-Genuchten sont disponibles dans le champ des matériaux, le programme évalue analytiquement la saturation effective, au lieu de tenter de chercher la saturation en eau dans le champ VARI_NOEU. Cela évite l’usage de [CALC_CHAMP] avec le résultat du calcul hydraulique avant de passer à [CALC_STAB_PENTE]. La Fig. 465 met en évidence l’erreur négligeable de l’évaluation analytique par rapport au résultat numérique et au résultat de référence. Cette figure prouve également la cohérence entre les résultats issus du champ de pression définie par CHAM_PRES et LIGN_PHREA.

Fig. 465 Evaluation analytique de la résistance due à la succion par le modèle Van-Genuchten (cf. [V6.03.009]).

Surface de rupture circulaire

Méthode de Fellenius

La méthode de Fellenius est appelée aussi «méthode suédoise» ou « ordinary method of slices » [bib2]. Cette méthode néglige complètement la force d’interaction entre les tranches. Ainsi la formule du coefficient de sécurité \(\mathit{FS}\) devient explicite. Le schéma des forces est montré dans la Fig. 466.

Fig. 466 Schéma des forces de la méthode Fellenius [bib2]

Étant donnés le poids des tranches, la pression répartie et la force inertielle sismique, la force normale appliquée sur la base des tranches s’écrit:

(4636)\[N = W\left[k_h \sin\alpha + (1 - k_v)\cos\alpha\right] + P \cos(\alpha - \beta)\]

où

• \(W\) désigne le poids,

• \(k_h\) et \(k_v\) désignent les coefficients d’accélération normalisés dans la direction horizontale et verticale respectivement.

• \(P\) est force totale de la pression répartie appliquée sur la surface des tranches.

• \(\alpha\) et \(\beta\) désignent l’angle d’inclinaison de la base des tranches et l’angle de la force de pression par rapport à l’horizontale.

En négligeant l’interaction entre les tranches, par résolution du bilan des forces dans le plan de la base des tranches, on en déduit la force de cisaillement subie par les tranches :

(4637)\[S = -W \left[ k_h \cos\alpha - (1-k_v) \sin\alpha\right] + P \sin(\alpha - \beta)\]

Enfin, avec la présence du champ de pression hydraulique, un terme supplémentaire lié à l’effet hydraulique s’ajoute à l’expression de la résistance :

(4638)\[\begin{split}\Omega_{\textrm{hydr}} = \left\lbrace \begin{array}{l} \presCapi \Delta l \cos^2 \alpha \tan {\anglFric}^{\prime} & \text{ si } \quad \presCapi < 0 \\ -\presAir \Delta l \cos^2 \alpha \tan {\anglFric}^{\prime} + \presCapi \Delta l \tan \anglFric_{\textrm{hydr}} & \text{ sinon} \end{array} \right.\end{split}\]

où \(\Delta l\) est la longueur de base des tranches, on a \(\anglFric_{\textrm{hydr}} = \anglFric_b\) si les valeurs de l’angle de succion matricielle sont fournies, sinon

(4639)\[\tan \anglFric_{\textrm{hydr}} = S_{\textrm{we}}(\presCapi)\tan {\anglFric}^{\prime}\]

avec \(S_{\textrm{we}}\) la saturation effective évaluée par le comportement hydraulique du sol.

Par définition l’expression explicite du facteur de sécurité s’écrit alors :

(4640)\[\mathit{FS} = \frac{\sum\cohesion \, {\Delta}l + N \tan{\anglFric}^{\prime} + \Omega_{\textrm{hydr}}}{\lvert \sumS \rvert}\]

Note

A cause de l’hypothèse très forte, la méthode Fellenius ne vérifie pas les bilans des forces et des moments en même temps, mais tous les chargements s’alignent dans la direction verticale. On a choisi arbitrairement la formule qui satisfait les bilans des forces. L’autre possibilité est d’exprimer la force de cisaillement par \(S = W \sin\alpha - \frac{M_n}{R}\) où \(M_n\) désigne le moment des chargements mécaniques à part de la gravité par rapport au centre de rotation et \(R\) est le rayon de rotation. Ainsi l’expression (4849) satisfait les bilans des moments.

Selon les tests, il s’avère que les deux formules du \(S\) conduisent aux valeurs de \(\mathit{FS}\): différentes, mais l’écart par rapport à la solution exacte dépendant du cas étudié. Il est donc difficile de juger laquelle est plus pertinente. De ce fait, la méthode Fellenius demeure un moyen de calcul simple mais peu précis. Il vaut mieux de passer en Bishop simplifié pour un meilleur équilibre entre l’efficacité et la fidélité.

En mettant \(\mathit{FS} = 1\) dans l’équation (4849) et supposant l’accélération sismique horizontale, on en déduit l’expression du coefficient d’accélération sismique :

(4641)\[k_c = \frac{ \sum\limits_i \pm \left[ W_i \sin\alpha_i + P_i \sin(\alpha_i - \beta_i) \right] - \left[ {\cohesion}^{\prime} \Delta l + \left[ W_i \cos\alpha_i + P_i \cos(\alpha_i - \beta_i)\right] \tan {\anglFric_i}^{\prime} + \Omega_{\textrm{hydr}, i} \right] }{ \sum\limits_i W_i \left[\sin\alpha_i \tan {\anglFric_i}^{\prime} \pm \cos\alpha_i \right] }\]

où la signe \(\pm\) est positif pour les pentes en amont et négative pour les pentes en aval.

Méthode de Bishop simplifiée

La méthode de Bishop simplifiée suppose que la force d’interaction entre les tranches est horizontale, comme le montre la Fig. 467.

Fig. 467 Schéma des forces de la méthode Bishop [bib2]

Le bilan des forces suivant la direction verticale s’écrit:

(4642)\[N\cos\alpha +S\sin\alpha -W(1 - k_v) - P \cos \alpha = 0\]

La contribution du champ de pression hydraulique en prenant en compte la succion matricielle s’écrit :

(4643)\[\Omega_{\textrm{hydr}} = \Delta l \cos\alpha \, \left[ -\presAir \tan {\anglFric}^{\prime} + \presCapi \tan\anglFric_{\textrm{hydr}} \right]\]

Selon le critère de Mohr-Coulomb, la force tangentielle de résistance s’écrit :

(4644)\[S = \frac{1}{\mathit{FS}} \left[ {\cohesion}^{\prime} \Delta l + N \tan {\anglFric}^{\prime} + \frac{\Omega_{\textrm{hydr}}}{\cos\alpha} \right]\]

On en déduit donc l’expression implicite du facteur de sécurité:

(4645)\[\mathit{FS} = \frac{ \sum\limits_i \frac{\deri{\cohesion} \Delta l \cos\alpha_i + \left[W_i(1 - k_v) + P_i \cos\beta_i\right] \deri{\anglFric} + \Omega_{\textrm{hydr}, i}} {\cos\alpha_i \pm \frac { \sin\alpha_i \tan\deri{\anglFric} }{\mathit{FS}}} }{ \lvert \sum\limits_i W_i \sin\alpha_i - \frac{ M_{n,i}}{ R} \rvert }\]

où \(M_{n,i} = M_{S,i} + M_{Q,i}\) représente le moment par rapport au centre de rotation des chargements mécaniques extérieurs (ref. (5009) et (5011)).

On résout l’équation (4852) par une méthode de point fixe.

Similaire à la méthode de Fellenius, on détermine le coefficient d’accélération critique de manière explicite tel que :

(4646)\[k_c = \frac{ \sum\limits_i W_i \sin\alpha_i - \frac{M_{n,i}}{R} - \sum\limits_i \frac{\deri{\cohesion} \Delta l \cos\alpha_i + \left[W_i(1 - k_v) + P_i \cos\beta_i\right] \deri{\anglFric} + \Omega_{\textrm{hydr}, i}} {\cos\alpha_i \pm (\sin\alpha_i \tan\deri{\anglFric})} }{ \sum\limits_i W_i \frac{y_c - y_{b,i}}{R} }\]

où \(y_c\) et \(y_b\) sont les ordonnées du centre de rotation et du barycentre des tranches respectivement.

La méthode de Bishop simplifiée satisfait à la fois les équations d’équilibre des moments pour la masse glissante et des forces verticales pour chaque tranche de sol, soit \(n+1\) équations au total.

Localisation de la surface critique

Afin de localiser la surface critique, on parcourt toutes les combinaisons possibles des points d’extrémités des surfaces potentielles de rupture. Pour chaque combinaison, on détermine les bandes d’essai du rayon par les deux approches suivantes :

• Option 1 : On parcourt tous les rayons qui forment une surface de rupture cinématiquement admissible. La limite inférieure est la valeur minimale parmi :

  • Le rayon du cercle tangentiel à la droite verticale \(x=x_b\) où \(x_b\) est l’abscisse de l’extrémité active (située à la crête de la pente) de la surface de rupture.

  • Le rayon du cercle tangentiel à la base du maillage.

  Pour la limite supérieure du rayon, on examine d’abord si tous les points de discontinuité convexes se situent en dehors du segment reliant les points d’entrée et de sortie. Si la condition est satisfaite, la limite supérieure est égale à l’infini, ce qui signifie qu’on parcourt les surfaces circulaires même à courbure très faible si le paramètre de discrétisation NB_RAYON est assez élevé.

  Sinon, on détermine la limite supérieure du rayon par un algorithme adaptatif qui choisit le rayon maximal en fonction du maillage fourni. Par analyse de maillage, on récupère la taille minimale globale des éléments 2D , noté \(s_{min}\), et les tailles verticales des éléments 2D contenant les neouds du groupe des mailles représentant le profil de pente défini par le mot-clef GROUP_MA, noté \(s_p\) qui est donc une fonction des coordonnées \(X\). À partir du rayon maximal des cercles passant par les points de discontinuité convexes, noté \(R_{max}\), on corrige itérativement \(R_{max}\) jusqu’à ce qu’il satisfasse les deux conditions suivantes :

  1. La distance entre le cercle et le profil de pente est supérieure \(s_{min}\) dans la direction verticale.

  2. Tous les points constituant la base des tranches s’éloignent du profil de pente à une distance verticale supérieure à \(s_p\) à la même abscisse.

• Option 2 : On définit manuellement la bande d’essai en définissant les droites horizontales tangentielles aux surfaces de rupture. Cette option de recherche accélère le calcul si l’utilisateur expérimenté peut anticiper la zone réduite où se situe probablement la surface de rupture. On profite également de cette option pour examiner les surfaces dans l’entourage de celle donnée par le résultat. Il est à noter qu’on peut utiliser Y_MINI ou Y_MAXI seul pour imposer la limite inférieure ou supérieure du rayon, le reste est détemriné par la procédure décrite dans l’option 1.

Afin de générer les surfaces potentielles de rupture distribuées uniformément dans l’espace de recherche, on définit la variable de contrôle \(d\) qui mesure la distance entre le segment reliant les extrémités, noté X1-X2, et la droite tangentielle du cercle parallèle à X1-X2, comme l’illustre la Fig. 126. Les limites du rayon étant converties en \(d_{min}\) et \(d_{max}\), on génère les surfaces circulaires avec un incrément constant :

(4647)\[\begin{split}\Delta d = \left\lbrace \begin{array}{l} \frac{d_{max}}{N_r} & \text{ si } \quad d_{min} = +\infty \\ \frac{d_{max} - d_{min}}{N_r - 1} & \text{ sinon } \end{array} \right.\end{split}\]

Fig. 468 Génération des cercles dans les méthodes Bishop et Fellenius.

Modifier le paramètre de discrétisation NB_RAYON est susceptible d’être insuffisant pour localiser la surface critique à l’issue de la recherche globale. Face à ce problème, après le calcul du facteur de sécurité correspondant à tous les rayons testés pour une combinaison des points d’extrémités, si un nouveau minima apparaît, on l’enregistre en faisant une vérification locale autour du rayon critique. Cette vérification se fait en raffinant itérativement \(\Delta d\) de manière exponentielle, et en testant les rayons autour du rayon qui minimise le facteur de sécurité.

A l’itération \(n\) , l’incrément raffiné s’écrit :

(4648)\[\Delta d_{n} = {2}^{-n} \, \Delta d_{0}\]

On teste alors les deux surfaces potentielles de rupture correspondants à :

(4649)\[{d}_{n}={d}_{\mathit{minFS},n-1} \pm \Delta d_{n}\]

Au cours des itérations, on enregistre toujours le rayon minimisant le facteur de sécurité et on continue à tester son voisinage avec l’incrément raffiné, jusqu’à la stabilisation du rayon.

Fig. 469 Visualization des surfaces circulaires parcourues dans le cas-test ssnp008 (cf. [V6.03.008]).

D’après les tests, cet algorithme d’optimisation conçu pour les surfaces circulaires s’avère efficace même pour les pentes inhomogènes. La Fig. 469 montre qu’avec la pente modélisée dans le cas test ssnp008 [V6.03.008], les surfaces testées et colorées en fonction du FS distribuent uniformément dans l’espace de recherche. La concentration des surfaces autour du point optimum reflète l’effet de l’optimisation locale permettant de localiser corretement la surface critique.

Surface de rupture non-circulaire

Calcul de FS par les méthodes Spencer ou Morgenstern-Price

Les méthodes Spencer et Morgenstern-Price sont conçues pour analyser la stabilité des pentes dont la surface de rupture est de forme plus complexe en raison de l’inhomogénéité des matériaux, telle que la stratification de la fondation et la présence d’une couche faible (cas test ssnp007 [V6.03.007]).

Fig. 470 Schéma des forces des méthodes Spencer et Morgenstern-Price [bib2]

La Fig. 470 montre le schéma des forces subies par les tranches quand on choisit les méthodes Spencer ou Morgentern-Price. \({F}_{v}\) désigne la composante verticale des chargements mécaniques, qui s’exprime par (4848) et \({F}_{h}\) désigne leur composante horizontale.

Ces méthodes satisfont toutes les équations d’équilibre, soit \(n\) équations d’équilibre des forces dans la direction verticale, \(n\) dans la direction horizontale et \(n\) équations d’équilibre du moment, où \(n\) est le nombre des tranches. Les deux méthodes se distinguent par l’hypothèse sur la direction des forces d’interaction entre les tranches \({Z}_{i}\) . On représente cette direction par l’angle \({\theta}_{i}\) par rapport à la direction horizontale :

(4650)\[\tan {\theta}_{i} = {f}_{0}({x}_{i}) + \lambda {f}_{1}({x}_{i})\]

où \({f}_{0}\) et \({f}_{1}\) sont les fonctions définissant la variation de l’angle d’inclinaison des forces d’interaction entre les tranches, et \(\lambda\) est un coefficient multiplicateur qui est une des inconnues des équations d’équilibre. Les fonctions \({f}_{0}\) et \({f}_{1}\) dépendent de la méthode choisie:

• Méthode Spencer : \({f}_{0}=0,{f}_{1}=1\) .

• Méthode Morgenstern-Price (proposée par Chen et Morgenstern en 1983) :

  • \({f}_{0}(x) = \tan{\beta}_{0}+\dfrac{\tan{\beta}_{1}-\tan{\beta}_{0}}{L}x\) , où \({\beta}_{0}\) et \({\beta}_{1}\) sont les angles d’inclinaison du profil de la pente aux deux extrémités de la surface de rupture.

  • \({f}_{1}(x)=\sin \left( \dfrac{\pi}{L}x \right )\)

On adopte la procédure de résolution du facteur de sécurité proposée par Zhu et al. [bib6]. L’avantage de cette procédure, comparée avec les procédures traditionnelles, est qu’elle en déduit le schéma d’itération du point fixe pour calculer le facteur de sécurité et propose une démarche faiblement couplée pour résoudre FS et \(\lambda\). Selon les auteurs, cette procédure de résolution est en même temps plus performante et plus efficace que les procédures traditionnelles.

• Le bilan d’équilibre en force dans la direction normale de la tranche \(i\) s’écrit :

  (4651)\[N_i = {W}_i (1 - k_v) \cos\alpha_i + k_hW_i\sin\alpha_i + Z_i \sin(\theta_i-\alpha_i) - Z_{i-1}\sin(\theta_{i-1}-\alpha_i) + P_i \cos(\beta_i - \alpha_i)\]

  où

  • \(N\): l’effort normal sur la base des tranches.

  • \(Z\): la force d’interaction entres les tranches.

• Le bilan d’équilibre en force dans la direction tangentielle de la tranche \(i\) s’écrit :

  (4652)\[\begin{split}\begin{array}{l} S_i &= \frac{N_i \tan {\anglFric}^{\prime}_i + \Omega_{\textrm{hydr},i} + {\cohesion}^{\prime}_i\Delta l}{F_S} \\ &= W_i(1 - k_v)\sin\alpha_i - P_i\sin(\beta_i - \alpha_i) - k_h W_i\cos\alpha_i + Z_i \cos(\theta_i-\alpha_i) - Z_{i-1}\cos(\theta_{i-1}-\alpha_i) \end{array}\end{split}\]

  où \(S\) désigne la force de cisaillement, \(\Omega_{\textrm{hydr}}\) rassemble les termes liés à l’effet de couplage hydro-mécanique :

  (4653)\[\begin{split}\Omega_{\textrm{hydr}} = \left\lbrace \begin{array}{l} \presCapi\Delta l \tan {\anglFric}^{\prime} & \text{ si } \quad \presCapi < 0 \\ -\presAir \Delta l \tan {\anglFric}^{\prime} + \presCapi \Delta l \tan \anglFric_{\textrm{hydr}} & \text{ sinon } \end{array} \right.\end{split}\]

En éliminant la force normale dans les équations (4794) et (4795), on a le bilan des forces simplifié :

(4654)\[{Z}_{i+1}{\Phi}_{i+1}={\Psi}_{i}{Z}_{i}{\Phi}_{i} \pm {R}_{i+1}-{F}_{s}{T}_{i+1}\]

où

(4655)\[\begin{split}\begin{array}{l} {\Phi}_{i}={F}_{s}\cos({\alpha}_{i}-{\theta}_{i}) \pm \sin({\alpha}_{i}-{\theta}_{i})\tan{\anglFric }_{i}\hfill \\ {\Psi}_{i}=[{F}_{s}\cos({\alpha}_{i+1}-{\theta}_{i}) \pm \sin({\alpha}_{i+1}-{\theta}_{i})\tan{\anglFric }_{i+1}]/{\Phi}_{i}\\ {R}_{i}=\left[{W}_{i}(1 - k_v)\cos{\alpha}_{i} + k_hW_i\sin\alpha_i + P_i\cos(\beta_i - \alpha_i)\right]\tan{\anglFric }_{i} + \Omega_{\textrm{hydr},i}+{c}_{i}\Delta l\hfill \\ {T}_{i}={W}_{i}(1 - k_v)\sin{\alpha}_{i} - P_i\sin(\beta_i - \alpha_i) - k_hW_i\cos\alpha_i\hfill \end{array}\end{split}\]

Les signes \(\pm\) sont positifs pour les pentes en amont et négatives pour les pentes en aval.

En fait, \({R}_{i}\) représente la force de résistance au cisaillement générée par toutes les forces appliquées sur la tranche sauf les forces d’interaction entre les tranches. \({T}_{i}\) représente la force d’entraînement conduisant à la déstabilisation de la tranche. Ainsi l’équation (4796) satisfait les bilans des forces complets dans toutes les directions. Si on rajoute l’hypothèse des forces d’interaction nulles, on retrouvera l’expression explicite de FS (4852) de la méthode Fellenius.

Ainsi, avec les conditions \({Z}_{0}={Z}_{N}=0\) (les forces d’interaction aux deux côtés de la partie glissante sont nulles), on en déduit le schéma de point fixe pour calculer FS:

(4656)\[\mathit{FS} = \frac{ \sum_{i=1}^{N-1} \left( \lvert{R}_{i}\rvert\prod_{j=i}^{N-1}{\Psi}_{j}+\lvert{R}_{N}\rvert \right) } {\sum_{i=1}^{N-1} \left( {T}_{i} \prod_{j=i}^{N-1}{{\Psi}}_{j}+{T}_{N} \right) }\]

Enfin, le bilan d’équilibre en moment par rapport au centre de la base de la tranche \(i\) s’écrit:

(4657)\[{E}_i \left( {z}_{i+1}+\frac{b}{2}\tan{\alpha}_i \right) = {E}_{i-1} \left( {z}_{i-1}-\frac{b}{2}\tan{\alpha}_i \right) + \frac{b}{2} \left( \tan{\theta}_{i}{E}_{i} + \tan{\theta}_{i-1}{E}_{i-1} \right) - M_{Q, i} - M_{S, i}\]

où

• \({E}_{i}\) est la composante horizontale de la force d’interaction qu’on peut calculer à l’aide de l’équation (4851),

• \(b\) est la largeur horizontale de la tranche,

• \({z}_{i}\) désigne la distance verticale entre le point où s’applique la force d’interaction et le centre de la base des tranches.

• \(M_{Q, i} = - P_i \sin\beta_i \left( y_{s,i} - y_{c, i} \right)\) est le moment de la pression répartie par rapport au centre de la base des tranches.

• \(M_{S, i} = - k_hW_i \left( y_{b, i} - y_{c, i} \right)\) est le moment de la force inertielle sismique par rapport au centre de la base des tranches.

Soit le moment de la composante horizontale par rapport au centre de la base \({M}_{i}={E}_{i}{z}_{i}\), on réécrit l’équation (4799) comme :

(4658)\[{M}_{i+1} = {M}_{i}-\frac{b}{2} \tan{\alpha}_{i+1}({E}_{i}+{E}_{i+1}) + \frac{b}{2} \left(\tan{\theta}_{i}{E}_{i}+\tan{\theta}_{i+1}{E}_{i+1}\right) - M_{Q, i+1} - M_{S, i+1}\]

Compte tenu de la condition \({M}_{0}={M}_{N}=0\), on en déduit l’expression de \(\lambda\) dans l’équation (4793):

(4659)\[\lambda =\frac{ \sum\limits_{i=1}^\limits{N} { \left[ \tan{\alpha}_{i} \left( {E}_{i}+{E}_{i+1} \right) - \left( {f}_{0,i}{E}_{i}+{f}_{0,i-1}{E}_{i-1} \right) \right] } + \frac{2}{b}\sum\limits_{i=1}\limits^{N}\left[M_{Q, i} + M_{S, i}\right] } { \sum_{i=1}^{N} \left( {f}_{1,i}{E}_{i}+{f}_{1,i-1}{E}_{i-1} \right) }\]

On résout FS et \(\lambda\) par l’évaluation alternative de ses valeurs à l’aide des équations (4798) et (48) jusqu’à la convergence des deux inconnues, comme l’illustre la Fig. 471.

Fig. 471 Diagramme de la résolution du facteur de sécurité des méthodes Spencer et Mogenstern-Price

D’après le résultat des tests, cette procédure de résolution de FS converge en deux itérations, comme le montre la Fig. 472.

Fig. 472 Convergence du FS et de lambda

Pour trouver le coefficient d’accélération critique qui est égal à la valeur absolue du \(k_h\) quand \(F_S = 1.0\), il suffit de remplacer dans l’équation (4798) la force tangentielle mobilisante \(T\) par la force déstabilisante due à la séisme \(A\) :

(4660)\[A_i = - W_i \left( \cos\alpha_i \pm \sin\alpha_i \tan {\anglFric}^{\prime}_i \right)\]

et la force de résistance \(R\) par la force tangentielle stabilisant la pente \(S\) :

(4661)\[S_i = P_i \sin (\beta_i - \alpha_i) - W_i \sin\alpha_i \pm : \left[ {\cohesion}^{\prime} \Delta l + \left( W_i \cos\alpha_i + P_i \cos( \beta_i - \alpha_i) \right) \tan {\anglFric}^{\prime}_i + \Omega_{\textrm{hydr},i} \right]\]

Algorithme d’optimisation EFWA

Procédure d’optimisation

Due au grand nombre de variables d’état, le problème de localisation de la surface critique non-circulaire s’avère trop difficile pour être résolu par la méthode essai-erreur. Récemment, grâce à l’émergence des algorithmes d’essaims intelligents, de nombreux chercheurs les utilisent pour résoudre le problème de localisation de la surface de rupture non-circulaire. Par exemple, l’algorithme d’essaim particulaire (PSO), l’algorithme génétique et l’algorithme de la colonie de fourmis ont montré leur capacité de trouver l’optimum avec un coût temporel acceptable.

Zheng et al. (2013) [bib4] ont proposé l’algorithme de feu d’artifice amélioré (EFWA). Xiao et al. (2019) [bib5] a pour la première fois appliqué l’EFWA à l’analyse de stabilité des pentes, et d’après ses résultats, l’EFWA aboutit à l’optimum avec moins d’itérations que d’autres algorithmes d’essaim [bib5].

Dans [:ref:`CALC_STAB_PENTE <u4.84.47>`] on adopte donc l’EFWA comme l’algorithme d’optimisation.

L’EFWA simule le processus d’explosion des feux d’artifices. Dans l’espace des variables d’état figuré comme le ciel nocturne, on fait exploser les feux d’artifice à différents endroits. Chaque explosion génère localement des étincelles pour chercher l’optimum dans l’espace voisin. Parmi ces étincelles, certaines qui présente la valeur sélective importante (inversement proportionnelle au FS) sont gardées pour la prochaine explosion.

L’EFWA se déroule en trois étapes:

1. Génération aléatoire des N feux d’artifices initiaux

   Le vecteur des variables d’état comportent \(N+1\) composants, qui sont les abscisses des points d’extrémité de la surface de rupture, et \(N-1\) qui sont les ordonnées des points intermédiaires, où \(N\) est le nombre des tranches :

   (4662)\[X=\lbrace {x}_{0},{y}_{1}, \cdots ,{y}_{N-1},{x}_{N}\rbrace\]

   Si on ne précise pas l’état initial via l’opérande ETAT_INIT de [CALC_STAB_PENTE], on génère aléatoirement \({N}_{\mathit{FA}}\) positions initiales des feux d’artifices (loi de distribution uniforme) :

   (4663)\[{X}_{i} = \rand \left( {x}_{\min,i},{x}_{\max,i} \right )\]

   On note \(\rand([i,j])\) est un générateur de nombres aléatoires dans l’intervalle \([i,j]\) et où les limites des degrés de libertés sont déterminées par l’algorithme décrit dans la section 3.4.2.2 .

   Si on fournit à ETAT_INIT la table en sortie de [CALC_STAB_PENTE], on l’ajoute dans la génération initiale des feux d’artifice.

2. Génération des étincelles ordinaires et gaussiennes

   Lors de la génération des étincelles ordinaires à partir de la position d’explosion \(i\), l’amplitude de l’explosion est calculée par :

   (4664)\[{A}_{i} = A \, \frac{\mathit{FS} \left({X}_{i} \right)-{\mathit{FS}}_{\min}+\varepsilon} {\sum_{i=1}^{N}(\mathit{FS}({X}_{i})-{\mathit{FS}}_{\min})+\varepsilon}\]

   où \(A\) est la valeur du mot-clef A de [CALC_STAB_PENTE] désignant l’amplitude totale des étincelles ordinaires, et \(\varepsilon=1.0 \times 10^-8\) est une constante proche de zéro pour éviter la division par zéro. Par l’équation (51) on cherche dans un intervalle plus large si la position centrale est de qualité insuffisante, et on se concentre sur l’espace proche autour de la valeur minimale de FS.

   Selon l’équation (51), l’amplitude d’explosion sera très faible autour de la surface minimisant le FS dans l’itération courante, ce qui restreint la capacité de recherche locale. Donc, on corrige l’amplitude d’explosion par une limite inférieure qui décroît au cours des itérations :

   (4665)\[A_{\textrm{min}}(N) = A_{\textrm{init}} - \frac{A_{\textrm{init}} - A_{\textrm{fina}}}{N_{\textrm{max}}}\sqrt{(2N_{\textrm{max}} - N)N}\]

   où \(N\) est le numéro de l’itération courante et \(N_{\textrm{max}}\) est le nombre maximum d’itération EFWA défini dans le mot-clef ITER_MAXI. Par défaut les paramètres \(A_{\textrm{init}}\) et \(A_{\textrm{fina}}\) sont fixés à 100% et 10% de la largeur des tranches. Il est également possible de les définir manuellement via les mots-clefs A_INIT et A_FINAL.

   Le nombre d’étincelles ordinaires d’une explosion est calculé par:

   (4666)\[{S}_{i} = M \, \frac{{\mathit{FS}}_{\max}-\mathit{FS}({X}_{i})+\varepsilon} {\sum_{i=1}^{N}({\mathit{FS}}_{\max}-\mathit{FS}({X}_{i}))+\varepsilon}\]

   où \(M\) est la valeur du mot-clef M de [CALC_STAB_PENTE] désignant le nombre total d’étincelles ordinaires. Par l’équation (52) on génère plus d’étincelles si la position centrale est de bonne qualité. Afin d’éviter l’émergence des nombres absurdes, on corrige :

   (4667)\[\begin{split}\widehat{S}_{i} = \left\lbrace \begin{array}{l} \round(aM) & \text{si } S_{i} < aM \\ \round(bM) & \text{si } S_{i} > bM \\ \round(S_{i}) & \text{sinon } \end{array} \right.\end{split}\]

   On génère les étincelles ordinaires en modifiant arbitrairement les degrés de liberté des positions de la génération courante:

   (4668)\[\widehat{X}_{i}^{k} = {X}_{i}^{k} + {\Delta}{X}^{k} \text{ avec } \Delta}{X}^{k}={A}_{i}\rand{-1,1})\]

   Afin de garder toujours la diversité de la recherche de l’optimum, on génère également les étincelles gaussiennes de quantité \({M}_{g}\) . On sélectionne de manière arbitraire les positions de feu d’artifice, et on modifie arbitrairement les degrés de liberté:

   (4669)\[\widehat{X}_{i}^{k} = {X}_{i}^{k} + \left( {X}_{B}^{k}-{X}_{i}^{k} \right) e \text{ avec } e=N(\left{0,1\right})\]

   où \({X}_{B}^{k}\) est la position correspondante au FS minimum de la génération courante.

   Si certains degrés de liberté sont en dehors des limites déterminés par l’algorithme décrit dans la section 3.4.2.2 , on les projette dans la bande admissible:

   (4670)\[\widehat{X}_{i}^{k} = {X}_{\min,i}^{k} + \rand(\left{0,1 \right}) \, \left( {X}_{\max,i}^{k}-{X}_{\min,i}^{k} \right)\]

   Du fait que la modification d’un degré de liberté influence la bande admissible des points intermédiaires à droite d’elle, on les ajuste aussi afin que la surface de rupture générée soit cinématiquement admissible.

   Enfin, au cas où la surface de rupture générée est tellement proche du profil de la pente qu’il est impossible de créer le groupe des mailles des tranches, on évalue la bande admissible des degrés de liberté en prenant en compte la taille des éléments 2D autour du profil de pente ou la MARGE_PENTE imposée, et on ajuste les paramètres des nouvelles surfaces par un algorithme récursif (ref. 3.4.2.2) en revérifiant les bandes admissibles des degrés de liberté qui suivent.

3. Sélection des étincelles pour la prochaine génération

   La nouvelle génération des feux d’artifice se produit au travers de l’algorithme de la roulette. La probabilité qu’une étincelle ordinaire ou gaussienne soit gardée dans la prochaine génération est:

   (4671)\[{p}_{i} = \frac{{\mathit{FS}}_{\max}-{\mathit{FS}}_{i}} {{\mathit{FS}}_{\max}-{\mathit{FS}}_{\min}}\]

La diagramme de l’algorithme EFWA est représenté sur la Fig. 473.

Fig. 473 Diagramme de l’algorithme EFWA

Contraintes dynamiques des variables d’état

Fig. 474 Illustration des contraintes des variables d’état

Afin d’assurer que la surface de rupture est toujours convexe et cinématiquement admissible, on adopte la procédure proposée par Cheng et al. pour déterminer les contraintes dynamiques des variables d’état [bib7].

Dans cette procédure, les limites des ordonnées des points intermédiaires dépendent d’un ou deux points précédents (à gauche). Prenons l’exemple dans la Fig. 474 :

• Borne inférieure :

  • La surface de rupture faisant un angle inférieur de 45 degrés avec le profil de la pente au point d’extrémité A, ce qui détermine la borne inférieure du point C.

  • On prolonge la droite AC jusqu’à l’abscisse du point D pour déterminer la borne inférieure du point D.

• Borne supérieure :

  • On relie CB et vérifie si CB intersecte le profil de la pente. Si c’est le cas, on trouve le point d’intersection d’ordonnée minimum J et relie CJ pour déterminer la borne supérieure D2 du point D. Sinon, D1 sera la borne supérieure du point D.

Afin d’assurer que la surface de rupture reste dans le domaine, on ajuste la borne inférieure au fond du domaine si elle se situe en dessous de ce dernier.

Note

1. Dans [CALC_STAB_PENTE] on adopte une stratégie stochastique qui choisit un nombre aléatoire de surfaces à chaque itération EFWA et génère / modifie les points intermédiaires de sens inverse, soit de l’extrémité active (en haut) à l’extrémité passive (en bas). Cela permet d’augmenter la diversité topologique des surfaces potentielles, et donc d’accélérer la convergence de l’EFWA. Il fait croître également la probabilité de tomber sur la solution exacte de surface critique grâce à l’espace de recherche élargi.

2. De manière similaire à la génération des surfaces circulaires, un algorithme récursif s’applique à la vérification de la cohérence des surfaces et à la correction automatique dans la phase d’explosion de l’EFWA. En l’absence du paramètre MARGE_PENTE explicitement défini par l’utilisateur, on récupère la taille des éléments 2D autour du profil de pente par analyser le maillage et on s’assure que la surface s’éloigne de la pente à une distance supérieure à la taille verticale de ces éléments aux bords.

Adaptation du maillage

Fig. 475 Diagramme de l’analyse de stabilité avec raffinement du maillage

L’analyse de stabilité via [CALC_STAB_PENTE] se basant sur le maillage, sa finesse influence significativement le résultat du facteur de sécurité. De ce fait, [CALC_STAB_PENTE] automatise le raffinement de maillage. Cependant, le temps de calcul augmente d’autant plus que le maillage est fin. De ce fait, comme le montre la Fig. 475, on isole le raffinement du maillage au moment de la recherche de la surface critique, étant donné que la finesse du maillage a une influence négligeable sur la localisation de la surface de rupture.

Afin d’accélérer la création des groupes des tranches, à l’aide de la commande MACR_ADAP_MAIL, on raffine le maillage uniquement sur le contour des tranches jusqu’à l’atteinte du nombre maximum de raffinement ou à la stabilisation de FS. Ce processus est piloté par un champ d’indicateur de proximité exprimé sous la forme d’une fonction gaussienne de moyenne nulle et d’écart type \(\sigma\) :

(4672)\[\mathit{val} = \frac{1}{\sqrt{2\pi }\sigma }\exp\left(\frac{-{d}_{\min}^{2}}{2{\sigma}^{2}} \right)\]

où \({d}_{\min}\) est la distance minimale entre un nœud du maillage et le contour des tranches. On raffine les mailles selon un critère donné par la valeur de \(\sigma\) , soit 68% des mailles les plus proches du contour, et on dé-raffine les mailles selon le critère de valeur \(2 \,\sigma\) , soit 5% des mailles les plus lointaines du contour.

On diminue de manière exponentielle l’écart type pour contrôler dynamiquement la zone de raffinement.

Fig. 476 Algorithme d’adaptation du maillage

La Fig. 476 illustre l’algorithme d’adaptation de maillage et le résultat du maillage adapté.

Validation et vérification des fonctionnalités

ZZZZ135	[V1.01.135]	Validation informatique des méthodes SRM dans la commande CALC_STAB_PENTE
ZZZZ146	[V1.01.146]	Validation informatique des méthodes LEM dans la commande CALC_STAB_PENTE
SSNP507	[V6.03.507]	Analyse de stabilité d’une pente argileuse non-drainée sur une fondation faible
SSNP006	[V6.03.006]	Analyse de stabilité d’une pente homogène par la méthode Bishop
SSNP007	[V6.03.007]	Analyse de stabilité d’une pente non-drainée avec couche faible par la méthode Morgenstern-Price
SSNP008	[V6.03.006]	Stabilité d’une pente à deux niveaux avec ligne phréatique et succion matricielle
SSNP009	[V6.03.006]	Stabilité d’une pente homogène avec la succion évaluée par la courbe de rétention d’eau
SSNP010	[V6.03.006]	Stabilité d’une pente subie à la force inertielle sismique
SSNP011	[V6.03.006]	Calcul du coefficient d’accélération critique des pentes homogènes et stratifiées

Bibliographie

[bib1]

Griffiths, D. V., & Lane, P. A. (1999). Slope Stability Analysis by Finite Elements. Geotechnique, 49, 387–403. doi:10.1680/geot.1999.49.3.387

[bib2] (1,2,3,4)

Duncan, J. M., Wright, S. G., & Brandon, T. L. (2014). Soil Strength and Slope Stability (2nd ed.). Wiley.

[bib3]

Wang, D., Chen, X., Qi, J., et al. (2021). Assessment on Strength Reduction Schemes for Geotechnical Stability Analysis Involving the Drucker-Prager Criterion. Journal of Central South University, 28, 3238–3245. doi:10.1007/s11771-021-4814-2

[bib4] (1,2)

Zheng, S., Janecek, A., & Tan, Y. (2013). Enhanced Fireworks Algorithm. In Proceedings of the 2013 IEEE Congress on Evolutionary Computation (CEC) (pp. 2069–2077). IEEE. doi:10.1109/CEC.2013.6557854

[bib5] (1,2,3)

Xiao, Z., Tian, B., & Lu, X. (2019). Locating the Critical Slip Surface in a Slope Stability Analysis by Enhanced Fireworks Algorithm. Cluster Computing, 22(1), 719–729. doi:10.1007/s10586-018-2304-9

[bib6]

Zhu, D. Y., Lee, C. F., & Qian, Q. H. (2005). A Concise Algorithm for Computing the Factor of Safety Using the Morgenstern-Price Method. Canadian Geotechnical Journal, 42(1), 272–278. doi:10.1139/t04-086

[bib7]

Cheng, Y. M. (2003). Location of Critical Failure Surface and Some Further Studies on Slope Stability Analysis. Computers and Geotechnics, 30(3), 255–267. doi:10.1016/S0266-352X(03)00021-X

[bib8]

GeoSlope International Ltd. (2021). Stability Modeling with GeoStudio (2021 Edition).

[bib9]

19. 11. Vanapalli, D. G. Fredlund, D. E. Pufahl, et A. W. Clifton (1996). Model for the Prediction of Shear Strength With Respect To Soil Suction, Can. Geotech. J., 33(3): 379-392