Application de la méthode de surface de
réponse stochastique à l'analyse de stabilité
d'un tunnel pressurisé

Guilhem Mollon<fnref n="1" />, Daniel Dias<fnref n="2" />, Abdul-Hamid Soubra<fnref n="3" />

RÉSUMÉ. La CSRSM (méthode de la surface de réponse stochastique par collocation) est explicitée et détaillée dans le cas d'un modèle d'analyse de la stabilité d'un front de taille pressurisé. La validation est effectuée par comparaison entre les distributions de probabilité de la variable de sortie du modèle analytique et de son approximation par CSRSM. Une étude paramétrique est ensuite conduite sur les données probabilistes d'entrée du problème.

ABSTRACT. CSRSM (Collocation-based Stochastic Response Surface Methodology) is described and detailed in the case of an analytical model for the stability of a pressurized tunnel face. The validation is carried out by comparison between the output distributions of the analytical model and of its approximation by CSRSM. A parametric study is then driven on the probabilistic input data of the problem

MOTS-CLÉS : Méthode non intrusive, CSRSM, étude paramétrique, Simulation de Monte-Carlo, Propagation d'incertitude

KEYWORDS: Non-intrusive methodology, CSRSM, parametric study, Monte-Carlo
Simulation, Uncertainty propagation

<footnote n="1">¹Doctorant, Laboratoire de Génie Civil et d'Ingénierie Environnementale, INSA Lyon, Université de Lyon, Domaine Scientifique de la Doua, 20 Avenue Einstein, 69621 Villeurbanne. E-mail : Guilhem.Mollon@insa-lyon.fr</footnote>

<footnote n="2">²Maître de Conférence HDR, Laboratoire de Génie Civil et d'Ingénierie Environnementale, INSA Lyon, Université de Lyon, Domaine Scientifique de la Doua, 20 Avenue Einstein, 69621 Villeurbanne. E-mail : Daniel.Dias@insa-lyon.fr</footnote>

<footnote n="3">³Professeur, GeM, Université de Nantes, Boulevard de l'Université, BP152, 44603 Saint-Nazaire cedex. E-Mail : Abed.Soubra@univ-nantes.fr</footnote>
# 1. Introduction

La stabilité du front de taille d'un tunnel a fait l'objet de nombreuses études par le passé, qu'il s'agisse d'approches expérimentales (Chambon et Corté 1994), analytiques (Leca et Dormieux 1990), ou numériques (Eisenstein et Ezzeldine 1994). Néanmoins, ces études considèrent toujours une situation purement déterministe dans laquelle les variables d'entrée sont prises en compte par l'intermédiaire de leur valeur moyenne. Un facteur de sécurité est considéré pour prendre en compte les incertitudes de manière approchée. La valeur cible de ce facteur est généralement laissée au jugement de l'ingénieur. Les outils de l'analyse probabiliste permettent de prendre en compte de manière plus rationnelle les incertitudes sur les variables d'entrées des modèles et d'évaluer l'impact de ces incertitudes sur la réponse du système étudié. Les deux obstacles les plus couramment rencontrés à ces méthodes sont (i) la difficulté à évaluer les distributions statistiques des variables d'entrée du modèle (par exemple les paramètres de cisaillement du sol), et (ii) le temps de calcul, puisqu'il est nécessaire d'appeler un grand nombre de fois le modèle mécanique pour conduire une étude probabiliste.

Une des voies pour lever le premier obstacle est d'effectuer une étude de sensibilité sur les distributions statistiques des variables d'entrée. Parmi les outils permettant de lever le deuxième obstacle, la méthode de surface de réponse stochastique par collocation (abrégée par CSRSM pour *Collocation-Based Stochastic Response Surface Method*) est une méthode non intrusive particulièrement efficace (Sudret 2007). Elle permet d'étudier un modèle numérique complexe de manière probabiliste et fiabiliste en le remplaçant par un méta-modèle analytique approché. On s'intéresse dans cette communication à la CSRSM par régression. Cette méthode consiste à approximer la réponse d'un modèle déterministe par une fonction analytique dont les coefficients inconnus sont déterminés par régression. Ce « méta-modèle » approché a un temps de calcul faible et peut donc être substitué au modèle d'origine dans les méthodes probabilistes classiques (Monte-Carlo etc.). Le modèle déterministe utilisé dans cette étude et le formalisme de la CSRSM (Isukapalli 1999, Sudret et al. 2006, Sudret 2007, Phoon et Huang 2007, Mollon et al. 2010) sont d'abord présentés, puis la CSRSM est validée par une approche de type Monte-Carlo. Finalement, une étude de sensibilité des données probabilistes d'entrée (distributions statistiques et corrélation) est effectuée.

# 2. Présentation du modèle déterministe

Le modèle déterministe utilisé dans cette étude est analytique et repose sur le théorème cinématique de l'analyse limite. Celui-ci permet la détermination de la pression nécessaire à appliquer au front de taille d'un tunnel percé en milieu frottant et/ou cohérent pour en garantir la stabilité (Mollon et al. 2009). Il s'agit d'un mécanisme tridimensionnel de rupture du front par plusieurs blocs de sol en translation. Les différents blocs sont de forme conique et subissent une translation selon leur axe. Les conditions de normalité imposées par l'analyse limite sont donc
respectées du fait que les surfaces de discontinuité de vitesse (i.e. les surfaces externes des blocs en mouvement) forment toujours un angle $φ$ avec le vecteur vitesse, $φ$ étant l'angle de frottement interne du sol. Cette conicité induit un défaut géométrique dans le modèle, puisque l'intersection du premier bloc avec le front de taille est nécessairement une ellipse verticale. Ceci va à l'encontre des observations obtenues à partir de modèles physiques ou numériques, pour lesquels l'ensemble du front circulaire est concerné par la rupture. Néanmoins, le modèle proposé donne des valeurs de la pression de rupture cohérentes avec celles issues des simulations numériques ou physiques (Mollon et al. 2009).

<figure><img src="image_1.png" /><figcaption>Figure 1. Allure du mécanisme de ruine proposé par le modèle déterministe</figcaption></figure>

Ce modèle relativement simple, implémenté sous Excel en Visual Basic, a été choisi par commodité pour son temps de calcul quasi-instantané (0.1s). La variable de sortie du modèle est donc la pression limite d'effondrement du front de taille (nommée $σ_c$ et exprimée en kPa), tandis que les variables aléatoires d'entrée sont au nombre de deux : l'angle de frottement interne du sol ($φ$), et la cohésion (c). Dans une approche probabiliste, ces deux variables peuvent être munies de différentes lois de distribution, et être éventuellement corrélées.

### 3. Présentation de la CSRSM

L'objectif de la CSRSM est de définir un méta-modèle approché qui vient se substituer au modèle déterministe dans les études probabilistes, et dont le temps de calcul est quasi-instantané. Ce méta-modèle (ou PCE pour *Polynomial Chaos Expansion*) s'exprime dans une base nommée chaos polynomial. Pour un chaos polynomial d'ordre $n$, les vecteurs de cette base sont des polynômes orthogonaux (Hermite, Legendre, etc.) multidimensionnels de degrés $\le n$. Pour un jeu donné de paramètres d'entrée (i.e. corrélations et distributions statistiques des différentes variables aléatoires d'entrée), on peut déterminer les coefficients du chaos par un
nombre réduit d'appels au modèle déterministe. Dans l'étude proposée, le chaos est formé par des polynômes d'Hermite multidimensionnels, et les coefficients du chaos sont déterminés par régression à partir de la réponse du modèle déterministe en un certain nombre de points de collocation. Cette partie décrit brièvement l'utilisation de la CSRSM dans le cas de deux variables. Pour une présentation plus générale, le lecteur peut se référer à Isukapalli (1999), ou à Sudret (2007).

Les deux variables aléatoires (φ et c) doivent être représentées dans le chaos polynomial par deux variables standards (variables normales centrées et réduites) nommées respectivement $\xi_1$ et $\xi_2$. Pour un PCE d'ordre n donné, les points de collocation disponibles sont déterminés dans l'espace standard ($\xi_1, \xi_2$) : chacune des variables standards peut prendre les valeurs des racines du polynôme d'Hermite de degré n+1, et les points de collocation disponibles résultent de l'ensemble des combinaisons de ces racines. Pour deux variables standards, la variable de sortie U peut être calculée à partir du chaos polynomial d'ordre n comme suit :

$$U = \sum_{i=1}^{p} a_i \cdot \Gamma_i(\xi_1, \xi_2) \quad (1)$$

Dans cette expression, les termes $\Gamma_i$ sont des polynômes d'Hermite multidimensionnels, qui constituent les vecteurs de la base du chaos polynomial. Les coefficients $a_i$ du chaos sont les inconnues du problème. Pour un PCE d'ordre n et un nombre $n_v$ de variables aléatoires d'entrée, le nombre p de coefficients et le nombre M de points de collocation disponibles sont donnés par (Isukapalli 1999) :

$$p = \frac{(n + n_v)!}{n!n_v!} \quad (2)$$

$$M = (n+1)^{n_v} + \begin{cases} 0 & \text{si l'une des racines du polynôme d'Hermite de degré n+1 est nulle} \\ 1 & \text{sinon} \end{cases} \quad (3)$$

Ces expressions permettent de remarquer que, dès lors que l'on considère un nombre important de variables aléatoires, le nombre de points disponibles est bien plus élevé que le nombre de coefficients à déterminer. Il est donc nécessaire d'opérer un choix judicieux des points de collocation pour lesquels le calcul déterministe sera effectué, afin de réduire les temps de calcul sans affecter la qualité de l'approximation. Isukapalli (1999) propose d'utiliser un nombre de points de collocation environ deux fois supérieur au nombre de coefficients du PCE. Il est également recommandé par le même auteur de disposer ces points symétriquement, et aussi près de l'origine que possible. Enfin, il est proposé d'ajouter aux points de collocation l'origine de l'espace des variables standard pour le cas où le polynôme d'Hermite de degré n+1 n'aurait pas de racine nulle, comme explicité dans l'équation (3). Dans la présente étude, le faible nombre de variables implique que $2p>M$ quelque soit l'ordre du PCE, il sera donc nécessaire d'utiliser tous les points de collocation disponibles.
Pour être introduits dans le modèle déterministe, les points de collocation (sous forme de couples $(\xi_{1,m}, \xi_{2,m})$ avec m compris entre 1 et M) doivent d'abord être exprimés dans l'espace des variables physiques (couples $(\varphi_m, c_m)$ correspondants), en utilisant les expressions suivantes :

$$ \begin{bmatrix} \xi_{1,C,m} \\ \xi_{2,C,m} \end{bmatrix} = H \cdot \begin{bmatrix} \xi_{1,m} \\ \xi_{2,m} \end{bmatrix} \qquad (4) $$

$$ \left\{ \begin{array}{l} \varphi_m = F_{\varphi}^{-1}[\Phi(\xi_{1C,m})] \\ c_m = F_{c}^{-1}[\Phi(\xi_{2C,m})] \end{array} \right. \qquad (5) $$

avec H : transformée de Cholesky de la matrice de corrélation de c et φ, Φ : fonction de répartition normale centrée réduite, F_φ : fonction de répartition de φ, et F_c : fonction de répartition de c. Après obtention de la réponse du modèle déterministe aux points de collocation, la détermination des coefficients du PCE s'effectue par régression, en résolvant l'équation matricielle suivante :

$$ \underline{N}^t \cdot \underline{N} \cdot \underline{a} = \underline{N}^t \cdot \underline{f} \qquad (6) $$

avec a : vecteur colonne des coefficients $a_i$, f : vecteur colonne des réponses du modèle aux points de collocation, et N la matrice suivante :

$$ \underline{N} = \begin{bmatrix} \Gamma_{0,0} & \Gamma_{1,0}(\xi_{1,1}) & \Gamma_{0,1}(\xi_{2,1}) & \Gamma_{2,0}(\xi_{1,1}) & \Gamma_{1,1}(\xi_{1,1}, \xi_{2,1}) & \Gamma_{0,2}(\xi_{2,1}) \\ \Gamma_{0,0} & \Gamma_{1,0}(\xi_{1,2}) & \Gamma_{0,1}(\xi_{2,2}) & \Gamma_{2,0}(\xi_{1,2}) & \Gamma_{1,1}(\xi_{1,2}, \xi_{2,2}) & \Gamma_{0,2}(\xi_{2,2}) \\ \cdot & \cdot & \cdot & \cdot & \cdot & \cdot \\ \Gamma_{0,0} & \Gamma_{1,0}(\xi_{1,M}) & \Gamma_{0,1}(\xi_{2,M}) & \Gamma_{2,0}(\xi_{1,M}) & \Gamma_{1,1}(\xi_{1,M}, \xi_{2,M}) & \Gamma_{0,2}(\xi_{2,M}) \end{bmatrix} \qquad (7) $$

La connaissance des coefficients du PCE fournit finalement une approximation analytique du modèle déterministe. Le post-traitement de ce PCE consiste dans la présente étude à déterminer la distribution statistique de la réponse par une méthode de Monte-Carlo. On peut aussi envisager d'appliquer à ce méta-modèle les outils de l'analyse fiabiliste (Monte-Carlo, Tirages d'Importance, FORM) afin d'étudier la probabilité de ruine due à une charge de service.

## 4. Validation de la CSRSM

Dans cette section, le méta-modèle fourni par la CSRSM est comparé au modèle déterministe original afin d'évaluer son exactitude et de choisir l'ordre optimal du PCE. Pour cette comparaison, on utilise des caractéristiques de références pour les deux variables d'entrée, qui sont considérées gaussiennes et non corrélées, avec les
caractéristiques suivantes : $\mu_{\phi}=17^{\circ}$ et $COV(\phi)=10\%$; $\mu_{c}=7$ kPa et $COV(c)=20\%$. Pour ce cas de référence, des tirages de Monte-Carlo sont effectués avec (i) les méta-modèles d'ordres 2, 3, 4, et 5, et (ii) le modèle déterministe original. Le nombre de tirages est de $10^6$. La figure 2a présente la densité de probabilité de la réponse pour le modèle déterministe et les 4 méta-modèles. Seul le PCE d'ordre 2 peut être discerné, tandis que les PCE d'ordres 3 et supérieurs fournissent une courbe de densité de probabilité confondue avec celle du modèle analytique. La figure 2b propose la probabilité de ruine du tunnel en queue de distribution et apporte plus de précision. L'ordre 2 donne des résultats peu satisfaisants, mais les ordres supérieurs se rapprochent de la probabilité de ruine obtenue pour le modèle analytique. Les ordres 4 et 5 fournissent des probabilités de ruines très similaires aux résultats obtenus par FORM (Mollon et al. 2009).

<figure><img src="image_2.png" /><figcaption>Figure 2. (a) Comparaison entre les PDF obtenus pour le modèle analytique et pour les PCE d'ordres 2 à 5 ; (b) Comparaison entre les probabilités de ruine obtenues pour le modèle analytique et pour les PCE, et les valeurs obtenues par FORM</figcaption></figure>

<table><caption>Tableau 1. Moments statistiques de la réponse obtenus pour les PCE d'ordres 2 à 5</caption><thead><tr><th rowspan="2">Ordre du PCE</th><th colspan="4">Moments Statistiques</th></tr><tr><th>Moyenne</th><th>Variance</th><th>Asymétrie</th><th>Aplatissement</th></tr></thead><tbody><tr><td>2</td><td>28.92</td><td>47.90</td><td>0.496</td><td>0.409</td></tr><tr><td>3</td><td>28.93</td><td>46.12</td><td>0.555</td><td>0.703</td></tr><tr><td>4</td><td>28.92</td><td>46.00</td><td>0.552</td><td>0.744</td></tr><tr><td>5</td><td>28.92</td><td>46.02</td><td>0.552</td><td>0.749</td></tr></tbody></table>

Le tableau 1 présente les quatre premiers moments statistiques de la distribution de la variable de sortie, pour les PCE d'ordres 2 à 5. Il apparaît que l'ordre 3 est suffisant pour obtenir une convergence sur la moyenne, la variance, et l'asymétrie, mais qu'il faut aller à l'ordre 4 pour que la convergence soit également satisfaisante en terme de coefficient d'aplatissement.
La figure 3 présente les surfaces de réponse fournies dans le plan (φ,c), en terme de lignes d'égales valeurs de pression d'effondrement (en kPa). Ces surfaces sont tracées pour le modèle analytique, ainsi que pour les PCE d'ordre 2 et 4. La correspondance est très satisfaisante pour les deux méta-modèles au centre de l'espace physique (c'est-à-dire dans la zone de forte probabilité d'occurrence des variables φ et c), mais l'ordre 4 est plus performant sur les zones « périphériques », correspondant aux queues de distribution. Ceci explique clairement les résultats fournis par la figure 2b. L'ordre 4 sera choisi pour conduire la suite de cette étude.

<figure><img src="image_3.png" /><figcaption>Figure 3. Comparaison entre les surfaces de réponses du modèle analytique et des PCE d'ordres 2 et 4 dans l'espace physique</figcaption></figure>

## 5. Etude de sensibilité des données probabilistes d'entrée

Le protocole décrit dans les sections précédentes est utilisé pour étudier la sensibilité de la réponse vis-à-vis des données probabilistes du problème. Il est en effet important de déterminer quelles caractéristiques des distributions de c et φ doivent être connues avec le plus de précision. La figure 4 montre l'impact du coefficient de variation (COV) de φ sur la distribution de la réponse. Le COV de φ varie de 5% à 15%, tandis que le COV de c est maintenu égal à 20%. Dans la figure 5, la démarche inverse est appliquée : le COV de c varie de 10% à 30%, tandis que celui de φ reste égal à la valeur de référence de 10%.
<figure><img src="image_4.png" /><figcaption>Figure 4. Impact du COV de φ sur le PDF (a.) et sur les moments statistiques (b.) de la pression critique d'effondrement</figcaption></figure>

<figure><img src="image_5.png" /><figcaption>Figure 5. Impact du COV de c sur le PDF (a.) et sur les moments statistiques (b.) de la pression critique d'effondrement</figcaption></figure>

Il apparaît que les COV des deux variables ont une influence déterminante sur la variance de la pression limite d'effondrement, et donc sur sa dispersion. Logiquement, un COV important pour c ou pour φ conduit à une plus grande dispersion de la réponse. Il apparaît que, pour un COV donné, l'impact de φ sur la dispersion de la réponse est plus important que celui de c. Néanmoins, cet effet est contrebalancé par la plus grande variabilité naturelle de c. Ceci confère aux deux variables un poids comparable dans une analyse probabiliste. On observe également que le COV de φ influence légèrement la moyenne, l'asymétrie, et l'aplatissement de la distribution de σ_c, ce qui n'est pas le cas du COV de c.
<figure><img src="image_6.png" /><figcaption>Figure 6. Impact des lois de distribution de c et φ sur le PDF de la réponse (a.) et sur la probabilité de ruine en queue de distribution (b.)</figcaption></figure>

La figure 6 présente l'influence de la non-normalité des distributions de c et φ sur le PDF de σ_c et sur la probabilité de ruine en queue de distribution. Le cas « non-normal » correspond aux distributions les plus couramment utilisées pour ces grandeurs : distribution lognormale pour c et distribution bêta pour φ (bornée par φ=0° et φ=45°). Cette comparaison est effectuée pour des COV standards (COV(φ)=10% et COV(c)=20%), et pour des COV dits « élevés » (COV standards majorés de 50%). Les figures 6a et 6b montrent que la non-normalité des variables d'entrée n'a pas une influence déterminante sur l'allure du PDF de la réponse et sur la probabilité de ruine, pour peu que les coefficients de variation de ces variables soient constants. Une différence importante apparaît en revanche si ces COV sont modifiés.

<figure><img src="image_7.png" /><figcaption>Figure 7. Impact de la corrélation entre c et φ sur le PDF (a.) et sur les moments statistiques (b.) de la pression critique d'effondrement</figcaption></figure>
La figure 7 présente l'influence du coefficient de corrélation entre $c$ et $φ$ sur le PDF de la pression limite d'effondrement et sur les quatre premiers moments statistiques de cette variable. Dans la communauté géotechnique peu d'informations sont disponibles sur le signe de cette corrélation, une plage de valeurs de $ρ_{cφ}$ comprise entre -0.6 à +0.6 a été retenue. Il apparaît que la corrélation entre les variables a une influence prépondérante sur la distribution de la pression limite d'effondrement, et en particulier sur sa dispersion. On observe ainsi que la variance de la réponse est multipliée par 4 lorsque la corrélation passe de $ρ_{cφ}=-0.6$ à $ρ_{cφ}=0.6$.

## 6. Conclusion

Cette communication rappelle les principes de la Méthode de la Surface de Réponse Stochastique par Collocation, et valide cette méthode par l'intermédiaire d'un modèle analytique relativement simple d'étude de la stabilité tridimensionnelle d'un front de taille pressurisé. Cette validation est effectuée par comparaison entre les résultats probabilistes obtenus par tirages de Monte-Carlo sur le modèle analytique et sur son approximation par chaos polynomial. Un ordre 4 pour ce chaos est jugé optimal en termes de précision et de temps de calcul. Dans une deuxième partie, une analyse de sensibilité est conduite vis-à-vis des données probabilistes des deux variables d'entrée $c$ et $φ$. Il apparaît que la prise en compte de lois non-normales ne modifie que très peu la distribution statistique de la pression limite d'effondrement, tandis que les coefficients de variation de $c$ et $φ$ et le coefficient de corrélation entre ces variables ont un effet déterminant sur cette distribution, et en particulier sur sa dispersion. Il apparaît donc souhaitable d'investiguer ces paramètres avec attention dans de futurs travaux.

## 7. Références Bibliographiques

[CHA 94] Chambon. P. and Corté, J.F. (1994). "Shallow tunnels in cohesionless soil: Stability of tunnel face.", J. of Geotech. Engrg., ASCE, 120(7), 1148-1165.

[EIS 94] Eisenstein, A.R., and Ezzeldine, O. (1994). "The role of face pressure for shields with positive ground control.", *Tun. and grd. conds.*, Balkema :Rotterdam, 557-571.

[ISU 99] Isukapalli, S.S. (1999). An uncertainty analysis of transport-transformation models, *Ph.D. thesis*, The State University of New Jersey, New Brunswick, New Jersey.

[LEC 90] Leca, E., and Dormieux, L. (1990). "Upper and lower bound solutions for the face stability of shallow circular tunnels in frictional material.", *Géot.*, 40(4), 581-606.

[MOL 09] Mollon, G., Dias, D., and Soubra, A.-H. (2009). "Probabilistic analysis and design of circular tunnels against face stability." *Int. J. of Geomech.*, ASCE, 9(6), 237-249.

[MOL 10] Mollon, G., Dias, D., and Soubra, A.-H. (2010). "Probabilistic Analysis of Pressurized Tunnels Against Face Stability Using Collocation-based Stochastic Response Surface Method" *J. of Geotech. & Geoenv. Engrg.*, Under review.

[PHO 07] Phoon, K.K., and Huang, S.P. (2007). "Geotechnical probabilistic analysis using collocation-based stochastic response surface method", *Applications of Statistics and Probability in Civil Engineering*, Kanda, Takada and Furada (eds.), London.

[SUD 06] Sudret, B., Berveiller, M., and Lemaire, M. (2006). "A stochastic finite element procedure for moment and reliability analysis", *Eur. J. Comput. Mech.*, (15)7-8, 825-866.

[SUD 07] Sudret, B. (2007). "Global sensitivity analysis using polynomial chaos expansion.", *Reliability Engineering and System Safety*, 93, 964-979.