r4.05.06 Méthode linéaire équivalent pour la propagation des ondes en 1D

Résumé

Ce document décrit la méthode linéaire équivalent et son implémentation dans Code_Aster dans la commande DEFI_SOL_EQUI [U4.81.31].

La méthode linéaire équivalent permet de calculer la réponse d’une colonne de sol unidimensionnelle (1D – une composante) ou multidimensionnelle (1D – trois composantes), stratifiée horizontalement à la propagation verticale d’ondes de cisaillement (SV). Elle rend compte de façon approchée du comportement des sols sous chargement cyclique. La méthode de résolution du linéaire équivalent est une procédure itérative, où on évalue à chaque itération pour chaque couche de sol, des caractéristiques viscoélastiques linéaires équivalentes à partir des courbes de dégradation du module de cisaillement \(G\) et de l’augmentation de l’amortissement hystérétique mesurées expérimentalement.

La simplicité de la mise en œuvre, la qualité des résultats et la rapidité des calculs en font une méthode très utilisée dans l’ingénierie. La méthode linéaire équivalent monodimensionnelle est ainsi implantée dans de nombreux codes de calcul (SHAKE, FLUSH, CYBERQUAKE, DEEPSOIL, EERA,…).

L’application de la méthode linéaire équivalent reste valable tant que la déformation volumique du sol peut être négligée. Dans le cas contraire, il est nécessaire de se tourner vers des modèles non-linéaires intégrant le couplage cisaillement-pression moyenne (ou déformations déviatoriques-dilatance) tel que le modèle élasto-plastique cyclique de Hujeux [R7.01.23] (Loi de comportement cyclique de Hujeux pour les sols).

Traditionnellement le modèle linéaire équivalent est admis jusqu’à des déformations de distorsion en champ libre de l’ordre de \({10}^{-3}\) (Guide ASN/2/01). Ceci reste un ordre de grandeur de validité, qu’il est nécessaire d’adapter à la situation et aux phénomènes que l’on souhaite modéliser.

Relation contrainte-déformation viscoélastique en une dimension

Le modèle linéaire équivalent est basé sur le modèle rhéologique viscoélastique de Kelvin-Voigt. Ce dernier, représenté en , est constitué en parallèle d’un ressort linéaire et d’un amortisseur linéaire représentant une résistance interne à la déformation.

Figure 2-1: Représentation du modèle rhéologique de Kelvin-Voigt.

La relation constitutive de ce modèle rhéologique est la suivante:

(1235)\[{\sigma}_{xz}=G\gamma +H\frac{\partial \gamma }{\partial t}=G\gamma +H\dot{\gamma}\]

avec \({\sigma}_{xz}\) contrainte en cisaillement, \(\gamma\) la distorsion(égale à \(2{\epsilon}_{xz}\) ), \(G\) module de cisaillement et \(H\) le module de viscosité. Le point représente la dérivée temporelle.

Dans le cas d’une colonne de sol en une dimension, en notant \(u(z,t)\) le déplacement horizontal à une profondeur \(z\) et à un instant \(t\) donnés, la distorsion et sa dérive temporelle sont définies à partir du champ de déplacement \(u(z,t)\) :

(1236)\[\gamma =\frac{\partial u(z,t)}{\partial z}\]

Et sa dérivée temporelle:

(1237)\[\dot{\gamma}=\frac{\partial \gamma (z,t)}{\partial t}=\frac{{\partial}^{2}u(z,t)}{\partial z\partial t}\]

Dans le cas d’un chargement harmonique \(u(z,t)=U(z){e}^{i\omega t}\) , la distorsion s’écrit:

(1238)\[\gamma (z,t)=\frac{\partial u(z,t)}{\partial z}=\frac{\mathit{dU}(z)}{\mathit{dz}}{e}^{i\omega t}=\Gamma (z){e}^{i\omega t}\]

Et sa dérivée temporelle:

(1239)\[\dot{\gamma}(z,t)=i\omega \gamma (z,t)\]

avec \(U(z)\) et \(\Gamma (z)\) respectivement l’amplitude du déplacement et de la distorsion.

La relation contrainte-déformation constitutive du modèle de Kelvin-Voigt est définie alors:

(1240)\[{\sigma}_{xz}(z,t)=G\gamma +H\dot{\gamma}=G(1+i\eta \omega )\gamma (z,t)={G}^{\ast }\gamma (z,t)\]

avec \({G}^{\text{*}}\) module de cisaillement complexe \({G}^{\ast }=(G+i\eta \omega )\) ; le facteur de perte est égal ici à \(\eta \omega\) .

Lors d’un chargement cyclique symétrique, la réponse du sol présentent des cycles ou des boucles d’hystérésis dans le plan contrainte-déformation (voir ).

Figure2-2. Réponse d’un sol à un chargement cyclique symétrique, \({\sigma}_{xz}\) la contrainte en cisaillement, \(\gamma\) la distorsion, \(D\) l’amortissement [4].

Ces boucles représentent la quantité d’énergie de déformation dissipée par le sol lors du chargement. Cette énergie peut être quantifiée en définissant le coefficient d’amortissement réduit \(D\) (voir la ) du sol par la relation suivante:

(1241)\[D=\frac{{W}_{d}}{4\pi {W}_{s}}\]

où \({W}_{d}\) représente l’énergie dissipée dans un cycle complet de chargement (égale à l’aire intérieure formée par la boucle contrainte-déformation, en gris sur la ), et \({W}_{s}\) l’énergie élastique emmagasinée par le sol lors du chargement pour le niveau de distorsion \({\gamma}_{c}\) (triangle inférieur, en hachuré sur la ).

Ces deux énergies sont calculées dans le cas d’un chargement harmonique \(\gamma (t)={\gamma}_{c}{e}^{i\omega t}\) avec le modèle de Kelvin-Voigt:

(1242)\[{W}_{d}=\underset{{{\sigma}_{xz}}_{c}}{\oint}{\sigma}_{xz}d\gamma =\underset{t}{\overset{t+\frac{2\pi }{\omega}}{\oint}}\mathrm{\Re }[{\sigma}_{xz}]\mathrm{\Re }[\frac{d\gamma }{\mathit{dt}}]\mathit{dt}=\pi \eta \omega {\gamma}_{c}^{2}\]

Et:

(1243)\[{W}_{s}=\frac{1}{2}{\tau}_{c}{\gamma}_{c}=\frac{1}{2}G{\gamma}_{c}^{2}\]

Pour le modèle viscoélastique de Kelvin-Voigt, l’énergie dissipée au cours d’un cycle de chargement dépend donc de la fréquence de sollicitation. Or l’expérience montre que pour les sols, l’énergie dissipée pendant le cisaillement est quasiment indépendante de la vitesse de déformation. L’amortissement résulte essentiellement de déformations plastiques irréversibles à l’échelle du grain. L’amortissement dans le sol est de nature plutôt hystérétique que visqueuse [5]. Le modèle de Kelvin-Voigt est ainsi modifié pour que les propriétés dissipatives soient équivalentes au matériau réel.

Relation contrainte-déformation à trois dimensions

On s’intéresse dans cette partie à définir une quantité équivalente scalaire en partant d’un état de contraintes et de déformations tensorielles. On part de la décomposition tensorielle du tenseur de déformation :math:`epsilon ` en partie volumique et déviatoriques:

(1244)\[\epsilon ={\epsilon}_{d}+\frac{1}{3}\mathit{tr}(\epsilon )I\]

où \({\epsilon}_{d}\) est le tenseur déviatorique des déformations. La modélisation linéaire équivalent s’intéressant au comportement en cisaillement de la colonne de sol, la norme du tenseur de déformation déviatorique \({\epsilon}_{d}\) est le candidat le plus approprié car on peut le relier à la déformation de cisaillement :math:`gamma ` pour un chargement monodimensionnel:

(1245)\[\gamma =\sqrt{3}{\epsilon}_{d}\]

Cette expression permet, à partir du calcul de \({\epsilon}_{d}\) obtenu par l’approche 1D – 3 composantes, de remonter à la valeur de :math:`gamma ` à utiliser pour faire évoluer le module élastique.

Pour le cas de la propagation d’ondes planes P, SV et SH perpendiculairement aux plans de stratigraphie (i.e. incidence verticale pour une stratigraphie horizontale), le calcul de la déformation déviatorique se résume à l’expression suivante:

(1246)\[ \begin{align}\begin{aligned}{\epsilon}_{d}=\sqrt{\frac{2}{3}\left({\epsilon}_{d}:{\epsilon}_{d}\right)}\\{\epsilon}_{d}=\frac{2}{3}\sqrt{{\epsilon}_{yy}^{2}+3{\epsilon}_{xy}^{2}+3{\epsilon}_{yz}^{2}}\end{aligned}\end{align} \]

La contrainte équivalente choisie est la contrainte déviatorique \(q\) , calculée en post-traitement via par l’expression suivante:

(1247)\[ \begin{align}\begin{aligned}q=\sqrt{\frac{3}{2}S:S}\\\textrm{et}\\q=3G{\epsilon}_{d}\end{aligned}\end{align} \]

où \(S\) est le tenseur déviatorique, défini par l’expression suivante:

(1248)\[ \begin{align}\begin{aligned}S=\sigma -\frac{1}{3}\mathit{tr}(\sigma )I\\S=\sigma -\frac{1}{3}\mathit{tr}(\sigma )I\end{aligned}\end{align} \]

Formulations du comportement linéaire équivalent

Formulation de Schnabel

Schnabel et al. (1972) [2] relient la viscosité à l’amortissement du matériau \(D\) , par \(H\omega =2\mathrm{GD}\) . L’amortissement critique du matériau \(D\) est la valeur mesurée en laboratoire et reportée sur les courbes \(D\) -gamma.

Le module de cisaillement complexe s’écrit alors:

(1249)\[{G}^{\ast }=(G+i\eta \omega )=(G+2i.\mathit{GD})=G(1+2i.D)\]

La raideur \(G\) est prise égale au module sécant \({G}_{s}\) des courbes (\({\sigma}_{xz}-\gamma\) ) (voir ). La valeur du module de cisaillement complexe est égale à:

(1250)\[∣{G}^{\ast }∣=G\sqrt{1+4{D}^{2}}\]

Le module du module de cisaillement complexe \(∣{G}^{\ast }∣\) et \(D\) augmentent donc conjointement. La représente la valeur de \(\frac{∣{G}^{\ast }∣}{G}\) en fonction de \(D\) . Par exemple, pour un coefficient d’amortissement réduit \(D\) de 20%, le module du module de cisaillement complexe est supérieur de 8% à la valeur de \(G\) .

Figure4.1-1. Évolution de \(\frac{∣{G}^{\ast }∣}{G}\) en fonction de \(D\) pour le modèle de Schnabel, comparé à celle pour le modèle de Lysmer.

Pour résoudre ce problème, Lysmer a formulé un second modèle linéaire équivalent.

Formulation de Lysmer

La formulation de Lysmer consiste à ré-écrire le module de cisaillement complexe \({G}^{\ast }\) pour s’assurer que le rapport \(\frac{∣{G}^{\ast }∣}{G}\) soit toujours égal à un. Le module de cisaillement complexe s’écrit alors:

(1251)\[{G}^{\ast }=G((1-2{D}^{2})+2\mathit{iD}\sqrt{1-{D}^{2}})\]

L’énergie dissipée dans un cycle est alors:

(1252)\[{W}_{d}=2\pi G\sqrt{1-{D}^{2}}{\gamma}_{c}^{2}\]

Comparaison des énergies dissipées par les deux formulations

La figure montre que le module complexe de cisaillement \(∣{G}^{\ast }∣\) est surévalué avec la formulation de Schnabel, et rigoureusement égal à \(G\) pour la formulation de Lysmer. L’énergie dissipée durant un cycle de chargement est également comparée entre les deux formulations en fonction de \(D\) .

Figure 4.3-1: Comparaison de l’énergie dissipée durant un cycle de chargement pour les deux formulations du linéaire équivalent en fonction de l’amortissement \(D\) [5].

L’énergie dissipée entre les deux formulations est la même jusqu’à environ \(D=0,2\) , puis la formulation de Lysmer dissipe moins d’énergie. Le domaine de validité du modèle linéaire équivalent fait que le domaine d’intérêt de ces courbes est sur l’intervalle \([0,0.3]\) . On remarque alors que même si \({D}_{\mathit{Schnabel}}<{D}_{\mathit{Lysmer}}\) , le modèle de Schnabel peut dissiper davantage d’énergie que le modèle de Lysmer.

Prise en compte simplifiée de la pression d’eau via le modèle de Byrne

La modélisation simplifiée de l’influence de la pression d’eau sur la réponse en cisaillement d’une colonne de sol avec DEFI_SOL_EQUI repose sur les hypothèses suivantes:

• les déformations volumiques plastiques développées pendant un chargement sismique sont essentiellement dues à la déformation de cisaillement y subie par le sol,

• incompressibilité de l’eau interstitielle, les variations de pression d’eau dans un sable saturé lors d’un chargement en cisaillement cyclique peuvent être estimées à partir des déformations volumiques obtenues lorsd’un essai de cisaillement cyclique réalisé en conditions drainées à la même contrainte de confinement.

Avec ces hypothèses, l’évaluation de l’augmentation de pression d’eau pendant un séisme (chargement en cisaillement cyclique) requiert évaluer les déformations volumiques plastiques générées par ce séisme en conditions saturées. Ces déformations sont fonction de la déformation volumique plastique cumulée et de l’amplitude de la déformation de cisaillement appliquée (Martin et al., 1975; Byrne1991). Martin et al. (1975) considèrent que le glissement entre grains générant un incrément de déformation volumique plastique \(\Delta {\epsilon}_{v}^{p}\) lors d’un essai en cisaillement drainé ou non drainé est le même. Toute contrainte hydrostatique transférée par le glissement des grains à l’eau interstitielle génère une augmentation de la pression de l’eau, selon la relation suivante:

(1253)\[\Delta u={E}_{r}\Delta {\epsilon}_{v}^{p}\]

où \(\Delta {\epsilon}_{v}^{p}\) représente l’incrément de déformation volumique plastique accumulé pendant une période de l’histoire de déformation de cisaillement, \(\Delta u\) est l’incrément de pression d’eau et \({E}_{r}\) est le module de décharge/recharge du squelette pour le niveau de contrainte effective actuel (Martin et al. 1975).

• Évaluation de \(\Delta {\epsilon}_{v}^{p}\)

L’évaluation de \(\Delta {\epsilon}_{v}^{p}\) est effectué à l’aide du modèle de Byrne (1991):

(1254)\[\frac{\Delta {\epsilon}_{v}^{p}}{\gamma}={C}_{1}\exp(-{C}_{2}\frac{{\epsilon}_{v}^{p}}{\gamma})\]

Où \({\epsilon}_{v}^{p}\) est le cumul de déformation volumique (somme des \(\Delta {\epsilon}_{v}^{p}\) ) depuis le début du calcul et \({C}_{1}\) et \({C}_{2}\) sont les paramètres du modèle, qui dépendent de la densité relative \({D}_{r}\) ou de la valeur de l’essai SPT \({({N}_{1})}_{60}\) . Byrne (1991) propose une relation empirique afin d’identifier les paramètres \({C}_{1}\) et \({C}_{2}\) à partir de la valeur \({({N}_{1})}_{60}\) :

(1255)\[{C}_{1}=\frac{8,7}{{({N}_{1})}_{60}^{1,25}} {C}_{2}=\frac{0,4}{{C}_{1}}\]

Ces expressions permettent d’obtenir l’incrément de déformation volumique plastique pour un chargement harmonique. Pour un chargement sismique (i.e. non uniforme), l’évaluation des déformations plastiques est faite par demi-cycle successifs de façon incrémentale, via une méthode de comptage de cycles par passage à zéro de la déformation de cisaillement. Cette méthode de comptage sur la déformation restreint l’usage de cette approche au cas 1D – 1 composantes uniquement.

• Évaluation de \({E}_{r}\)

L’évaluation de \({E}_{r}\) est faite à l’aide de l’expression proposée par Wu (1995), déterminée à partir des essais expérimentaux:

(1256)\[{E}_{r}=M\sigma {'}_{v}=M(\sigma {'}_{v0}-u)\]

où la valeur de \(M\) dépend de la densité relative \({D}_{r}\) ou de la valeur de l’essai SPT \({({N}_{1})}_{60}\) . Wu (2001) propose la corrélation suivante pour le calcul de \(M\) :

(1257)\[M=10{({N}_{1})}_{60}+160\]

Cette expression est intégrée dans la commande DEFI_SOL_EQUI lors de l’évaluation du module \({E}_{r}\) .

• Évaluation du taux de pression interstitielle \({r}_{u}\)

En combinant les équations et on peut estimer directement le rapport \({r}_{u}\) par l’équation suivante:

(1258)\[{r}_{u}=\frac{\Delta u}{\sigma {'}_{v0}}=1-\exp(-M{\epsilon}_{v}^{p})\]

La valeur de \({r}_{u}\) est estimée pour chaque demi-cycle de façon incrémentale, à partir du calcul de \({\epsilon}_{v}^{p}\) obtenu à l’aide l’expression .

On suppose que l’impact de la pression d’eau dans le module de cisaillement peut être estimé à l’aide de l’expression suivante:

(1259)\[{G}_{max}={G}_{0max}{(1-{r}_{\mathit{ueff}})}^{0,5}\]

où \({G}_{0max}\) est le module de cisaillement maximal à la pression de confinement in-situ et \({r}_{\mathit{ueff}}\) est estimé comme un facteur \({\mathrm{\chi }}_{\mathit{ru}}\) de la valeur maximale de \({r}_{u}\) , \({r}_{\mathit{umax}}\) , de façon similaire au coefficient de pondération de la méthode linéaire équivalent classique § 7 .

(1260)\[{r}_{\mathit{ueff}}={\mathrm{\chi }}_{\mathit{ru}}{r}_{\mathit{umax}}\]

Donc, la méthode présentée ci-dessous et appliquée par Kteich et al. (2018) à l’estimation de liquéfaction et tassements sismo-induits à la ville d’Urayasu au Japon pendant le séisme de 2011 consiste à enrichir le comportement en cisaillement décrit par les courbes \(G/{G}_{\max}\) en fonction de la déformation de cisaillement, en tenant compte de l’influence de la pression d’eau pour le calcul du module de cisaillement dégradé § 6.3 . Néanmoins, elle n’intervient pas pour l’obtention de la valeur d’amortissement hystérique à partir des courbes \(D\) en fonction de la déformation effective.

Implantation du linéaire équivalent dans code_aster

La commande DEFI_SOL_EQUI[U4.84.31] permet de calculer la réponse d’une colonne de sol unidimensionnelle à une sollicitation sismique avec le modèle linéaire équivalent (formulation de Schnabel ou Lysmer) en convolution ou en déconvolution.

La méthode de résolution du modèle linéaire équivalent est une procédure itérative, où sont évaluées à chaque itération pour chaque couche des caractéristiques linéaires équivalentes à partir des courbes de dégradation du module de cisaillement \(G\) et de l’augmentation de l’amortissement hystérétique.

La commande fournit en sortie les signaux temporels, les Spectres de Réponse d’Oscillateur (SRO), les déformations et contraintes temporelles à différentes profondeurs ainsi que les propriétés dégradées de sol en accord avec les niveaux de déformations de cisaillement calculées.

L’utilisateur doit donner en entrée de la commande:

• Une table donnant le rapport \(G/{G}_{max}\) en fonction de Gamma,

• Une table donnant l’amortissement réduit du matériau \(D\) en fonction de Gamma,

• Une table décrivant les caractéristiques matériaux initiaux des différentes couches (module de Young E, masse volumique RHO, coefficient de poisson NU, l’amortissement hystérétique AMOR_HYST, éventuellement la valeur de \({({N}_{1})}_{60}\) N1 pour l’obtention des paramètres du modèle de Byrne). L’amortissement hystérétique est égal à deux fois l’amortissement critique \(\mathit{AH}=2D\) (voir justification au § 4.2 .).

• Le maillage 2D/3D de la colonne respectant une taille de maille \({L}_{max}<\frac{{V}_{s}}{8{f}_{\mathit{coupure}}}\) ,

• Un signal temporel (unecomposante ou troiscomposantes) imposé en condition de champ libre ou en rocher affleurant.

Conditions limites du modèle

Le modèle sous-jacent est une colonne de sol 2D avec les conditions aux limites suivantes: Bords droit et gauche de la colonne : condition de périodicité. Cela signifie que les déplacements des nœuds des faces gauche et droite en vis-à-vis sont rendus égaux. Des couches horizontales infinies sont ainsi modélisées. Bas de la colonne : affectation d’un élément de frontière absorbante (cf. [R4.02.05] Éléments de frontière absorbante). Le substratum est ainsi modélisé comme un milieu infini. Selon le coefficient donné à l’élément, une partie des ondes se réfléchit à l’interface, tandis que l’autre partie se dissipe (par amortissement radiatif). La frontière absorbante adopte la modélisation D_PLAN_ABSO. Les autres éléments adoptent la modélisation D_PLAN.

Formulations du modèle linéaire équivalent

La commande modélise l’amortissement par un amortissement hystérétique (opérateur DEFI_MATERIAU mot clé AMOR_HYST). À partir du module de Young \(E\) et de l’amortissement hystérétique \(\mathrm{AH}\) renseigné dans DEFI_MATERIAU, code_aster construit le module complexe \({E}^{\ast }\) selon l’équation:

(1261)\[{E}^{\ast }=E(1+\mathrm{i.AH})\]

La construction du module complexe correspond donc à la formulation de Schnabel décrite au § 4.1 si \(\mathit{AH}=\mathrm{2D}\) . C’est pourquoi il faut bien préciser \(2{D}_{\min}\) en valeur initiale dans la table de sol.

Cependant si l’utilisateur a choisi la formulation de Lysmer, la commande calcule à partir des données matériaux \(E\) et \(\mathit{AH}\) , un module \({E}_{1}\) et un amortissement \({\mathit{AH}}_{1}\) équivalent tel que la formulation de Lysmer soit retrouvée:

(1262)\[{E}_{1}\left(1+i{\mathit{AH}}_{1}\right)=E\left(\left(1-\frac{{\mathit{AH}}^{2}}{2}\right)+ii\mathit{AH}\sqrt{1-\frac{{\mathit{AH}}^{2}}{4}}\right)\]

En égalisant les parties réelles et imaginaires de deux équations précédentes, les relations suivantes sont obtenues:

(1263)\[ \begin{align}\begin{aligned}{E}_{1}=E\left(1-\frac{{\mathit{AH}}^{2}}{2}\right)\\\textrm{et}\\{\mathit{AH}}_{1}=\mathit{AH}\sqrt{1-\frac{{\mathit{AH}}^{2}}{4}}\end{aligned}\end{align} \]

La méthode linéaire équivalent ne présente pas sous cette forme de dépendance de l’amortissement à la fréquence.

Méthode de résolution

La procédure de résolution du problème de la commande est la procédure itérative suivante:

1. Affectation des modules de cisaillement et d’amortissement initiaux (\(G={G}_{\max}\) , \(D={D}_{\min}\) ) en chaque point de discrétisation de la colonne de sol.

2. Calcul des modules complexes de cisaillement selon la formulation du modèle linéaire équivalent choisie (Schnabel ou Lysmer) à partir de \(G\) et \(D\) renseignés.

3. Calcul de la réponse harmonique de la colonne de sol: excitation par un bruit blanc unitaire à la base de la colonne (bord noté RA): \({A}_{\mathit{RA}}^{\mathit{harm}}=1\) sur la bande de fréquence choisie DYNA_VIBRA, TYPE_CALCUL = ‘HARM’.

4. Les accélérations harmoniques à chaque profondeur sont récupérées et les fonctions de transfert entre le rocher affleurant RA, et le champ libre CL, et les différentes couches \(i\) sont calculées.

En notant \({A}^{\mathit{harm}}{}_{i}\) , l’accélération harmonique de la couche \(i\) , il vient:

\(\mathit{Fonction}\mathit{de}{\mathit{transfert}}_{\mathit{CL}/\mathit{RA}}=\frac{1+{A}^{\mathit{harm}}{}_{\mathit{CL}}}{{A}^{\mathit{harm}}{}_{\mathit{RA}}}=\frac{1+{A}^{\mathit{harm}}{}_{\mathit{CL}}}{1}\) \(\mathit{Fonction}\mathit{de}{\mathit{transfert}}_{i/\mathit{RA}}=\frac{1+{A}^{\mathit{harm}}{}_{i}}{1}\) Le calcul étant en mouvement relatif (CHARGEMENT=’MONO_APPUI’), il convient d’ajouter le mouvement d’entraînement aux accélérations harmoniques de chaque couche. Ici, le mouvement d’entraînement est \(1\) , ce qui explique l’expression du numérateur de la fonction de transfert \(1+{A}^{\mathit{harm}}{}_{\mathit{CL}}\) . Si le calcul est en mouvement absolue (CHARGEMENT=‘ONDE_PLANE’), les fonctions de transfert se réduisent à la réponse harmonique de chaque couche.

1. Calcul de la transformée de Fourier (FFT) du signal d’entrée (passage en fréquentiel).

2. Convolution/déconvolution du signal d’entrée dans la colonne (Multiplication des fonctions de transferts par la FFT du signal d’entrée). Si le signal imposé est en rocher affleurant RA, il est convolué en utilisant les fonctions de transfert définies ci-dessus. Si le signal imposé est en champ libre CL, on calcule le signal rocher affleurant correspondant avec la fonction de transfert RA/CL pour CHARGEMENT=’MONO_APPUI’:

\(\mathit{Fonction}\mathit{de}{\mathit{transfert}}_{\mathit{RA}/\mathit{CL}}=\frac{1}{1+{A}^{\mathit{harm}}{}_{\mathit{CL}}}\)

1. Calcul des déformations \(\gamma (f)\) dans chaque couche.

Pour le cas d’une modélisation 1D – 3 composantes, la valeur des déformations \({\epsilon}_{xy}(f)\) , \({\epsilon}_{yz}(f)\) et \({\epsilon}_{yy}(f)\) sont récupérées afin de calculer la déformation équivalente.

2. Calcul des transformées de Fourier inverse \(\gamma (t)\) (retour en temporel).

Pour le cas d’une modélisation 1D – 3 composantes, la valeur de \({\epsilon}_{d}(t)\) est obtenue à partir des déformations \({\epsilon}_{xy}(t)\) , \({\epsilon}_{yz}(t)\) et \({\epsilon}_{yy}(t)\) par l’expression suivante:

\({\epsilon}_{d}=\frac{2}{3}\sqrt{{\epsilon}_{yy}^{2}+3{\epsilon}_{xy}^{2}+3{\epsilon}_{yz}^{2}}\)

La valeur de \(\gamma (t)\) est ensuite calculée à partir de \({\epsilon}_{d}(t)\) .

• Calcul dans chaque couche de la déformation de cisaillement effective définie comme:

\({\gamma}_{\mathit{effective}}=\mathrm{\chi }.{\gamma}_{max}\)

avec \({\gamma}_{\max}\) déformation de cisaillement maximale calculée dans la couche et \(\chi\) coefficient de pondération de la déformation effective (COEF_GAMMA, par défaut est égal à \(0,65\) , cf. §7).

• Évaluation pour chaque couche du module de cisaillement et de l’amortissement correspondant à \({\gamma}_{\mathit{effective}}\) sur les courbes \(G/{G}_{max}\) -Gamma et \(D\) -Gamma correspondantes.

Dans le cas de la modélisation simplifiée de la pression d’eau via le modèle de Byrne, la formule suivante est utilisée pour l’évaluation du module de cisaillement:

\(G/{G}_{max}({\gamma}_{\mathit{effective}}).{(1-{\mathrm{\chi }}_{\mathit{ru}}{r}_{\mathit{umax}})}^{0,5}\)

où \({r}_{\mathit{umax}}\) représente le maximum de l’historique d’évolution de \({r}_{u}\) et \({\mathrm{\chi }}_{\mathit{ru}}\) est un coefficient de pondération qui permet de moduler l’impact de \({r}_{\mathit{umax}}\) sur la déformation effective (COEF_KSI, par défaut est égal à \(0,667\) ). On multiplie donc la valeur de \(G/{G}_{max}\) évaluée à \({\gamma}_{\mathit{effective}}\) par l’expression prenant compte de la montée de la pression d’eau.

L’amortissement est évaluéavec l’expression de \({\gamma}_{\mathit{effective}}\) sans prise en compte de la pression d’eau.

• Si \(∣\frac{{G}_{k+1}-{G}_{k}}{{G}_{k}}∣>\theta\) (avec \(\theta =5\text{\%}\) par défaut (RESI_RELA)), retour à l’étape 2 avec les nouveaux modules \({G}_{k+1}\) pour chaque couche.

• Si la convergence est atteinte, le post-traitement s’effectue. La commande calcule les accélérations, contraintes, déformations, en temporel (FFT inverse), et les SRO pour chaque couche. Les propriétés dégradées de la colonne sont aussi données.

Choix du coefficient de pondération de la déformation effective

Le coefficient de pondération de la déformation effective \(\chi\) est un paramètre-clé reliant les conditions de chargement durant le séisme aux essais de laboratoire qui sont réalisés pour mesurer les courbes \(G/{G}_{\max}\) et \(D\) -gamma (Kramer, 1996 [7]). Sa valeur est prise inférieure à \(1\) . L’introduction de ce paramètre dans la théorie du modèle linéaire équivalent est pratique car il permet d’ajuster les résultats calculés à une réponse enregistrée. Cependant lorsque l’on souhaite prédire la réponse d’une colonne de sol, il faut faire le choix d’une valeur de \(\chi\) .

Dans la littérature, les pratiques et les recommandations sur la valeur de ce coefficient sont assez différentes:

• Dans le manuel d’utilisation de SHAKE (2012 et versions antérieures) [8], Schnabel et al. (2012) recommande de le prendre égal à une valeur de \(0.65\) . Cette valeur est aujourd’hui restée comme la valeur référence utilisée dans la pratique de l’ingénierie sans plus de considération. Cependant cette recommandation reste empirique.

• Dans le manuel d’utilisation de DEEPSOIL (2012) [9], Y. Hashash recommande de prendre ce paramètre à une valeur de \(0,65\) .

• Dans le manuel d’utilisation de SHAKE91 (1992) [10], Idriss and Sun recommande de faire varier ce coefficient en fonction de la magnitude \({M}_{w}\) du séisme du signal d’entrée, selon la formule:

(1264)\[\chi =\frac{({M}_{w}-1)}{10}\]

• Selon Idriss and Sun [10], les valeurs de \(\chi\) sont ainsi comprises dans l’intervalle \([0.4,0.75]\) , ce qui correspond à des séismes de magnitude \({M}_{w}\) de \(5\) à \(8,5\) .

• Dans le manuel d’utilisation de EERA (2000) [5], ce coefficient est noté \({R}_{y}\) et dépend de la magnitude du séisme selon la formulation de Idriss et Sun [10] ci-dessus.

• Dickenson (1994) [11] montre que la réponse calculée avec \(\chi =0,65\) peut être légèrement améliorée en utilisant une valeur de \(0,35\) à \(0,55\) pour des séismes de magnitude \({M}_{s}=6-7\) , et entre \(0,55\) et \(0,7\) pour des séismes de magnitude \({M}_{s}=7-8\) , \({M}_{s}\) étant la magnitude des ondes de surface.

• Yoshida et al (2002) [12] estiment que la valeur de \(\chi\) doit varier selon le niveau sismique du signal enregistré et qu’il varie en réalité entre \(0,2\) et \(1\) . Il précise qu’il n’y a pas de bonne méthode pour l’évaluation de \(\chi\) . Yoshida et al. montrent que l’utilisation du coefficient de pondération pour l’évaluation du module de cisaillement et de l’amortissement est à l’origine de deux défauts dans la théorie du linéaire équivalent:

1. • La surestimation de la contrainte de cisaillement maximale qui implique une surestimation de l’accélération maximale au pic.

Le mécanisme est décrit sur la . Le trait plein O-A-B représente la courbe contrainte-déformation du sol calculée à partir de la relation \(G\) -Gamma donnée par la courbe d’essai. Prenons l’exemple où une déformation maximale \({\gamma}_{\max}\) est calculée. L’estimation du module de cisaillement (module sécant - ligne droite OA) est faite par rapport à la déformation effective \({\gamma}_{\mathrm{eff}}\) . Or ceci implique qu’avec ce module de cisaillement sécant (OA) à la déformation \({\gamma}_{\max}\) , la contrainte de cisaillement maximale devient \({\tau}_{1}\) et non plus \({\tau}_{2}\) . Il y a donc une surestimation de la contrainte maximale de cisaillement qui implique une surestimation de l’accélération au pic. Cette surestimation reste modérée à faible déformation mais peut être importante pour des séismes forts. Yoshida et al. (2012) montre qu’à forte déformation (sur le palier plastique), la surestimation est en \(\frac{1}{\chi }\) , soit \(\mathrm{50 }\text{\%}\) pour \(\chi =0,65\) . Pour pallier cette difficulté, la première idée est de prendre \(\chi\) égal à \(1\) mais cette solution implique un autre problème expliqué plus bas.

Figure 7-1: Schéma illustrant la surestimation de la contrainte maximale de cisaillement avec la méthode linéaire équivalent [12].

1. • La sous-estimation de l’amplification à haute fréquence

Yoshida et al. [12] illustrent le second défaut de la méthode linéaire équivalent par la où l’amplification (ratio des FFT des accélérations) a été calculée avec SHAKE pour un sondage instrumenté (borehole) dans la baie de Tokyo et comparée au ratio mesuré. L’amplification calculée par SHAKE est inférieure à partir de \(7\mathrm{Hz}\) . Yoshida et al. [12] montrent que le problème vient du fait que le module de cisaillement et l’amortissement sont évalués à partir de la déformation effective, et cela pour toutes les fréquences. Le module de cisaillement est alors plus faible et l’amortissement plus grand que nécessaire pour le domaine des hautes fréquences, où les amplitudes sont plus petites que la déformation effective. Ce point montre la simplification d’utiliser un même coefficient de pondération de la déformation effective pour le comportement à hautes et basses fréquences.

C’est la raison pour laquelle certains auteurs recommandent de faire varier \(\chi\) en fonction de la magnitude. \(\chi\) doit prendre des valeurs basses pour des séismes faibles et des valeurs plus importantes pour des séismes plus forts.

Figure 7-2. Comparaison de l’amplification entre SHAKE et les mesures du borehole [12].

Les points 1 et 2 sont donc contradictoires puisqu’il faudrait augmenter \(\chi\) pour abaisser la surestimation de la contrainte maximale de cisaillement (et donc de l’accélération au pic) et réduire \(\chi\) pour éviter la sous-estimation de l’amplification à haute fréquence.

Le choix de \(\chi\) doit donc se faire selon le problème posé en tenant compte des deux remarques précédentes.

Des modifications du modèle linéaire équivalent ont été proposées pour résoudre ces difficultés (Modèle FDEL, Assimaki et al. [13], Modèle DYNEQ Yoshida et al. [12]). Dans ces modèles, le coefficient de pondération de la déformation effective, le module de cisaillement et l’amortissement dépendent de la fréquence.

Ces modèles restent cependant peu utilisés dans la pratique du fait des résultats nuancés.

Pour le modèle FDEL, Assimaki et al [13] et Furumoto (2000) [14] montrent une bonne amélioration des résultats par rapport au modèle classique tandis que Yoshida et al. [12] et Kwak et al. [15] mettent en avant des résultats moins bons que le modèle classique avec le modèle FDEL.

Pour le modèle DYNEQ, Yoshida et al. [12] montrent une amélioration certaine dans la qualité des résultats obtenus particulièrement en déconvolution par rapport au modèle classique (décrit ici) et au modèle FDEL. Cette amélioration est néanmoins à nuancer par les travaux de Kwak et al. [15].

Choix du facteur de pondération de la montée de pression d’eau effective

Le facteur de pondération de la montée de pression d’eau effective \({\chi }_{\mathit{ru}}\) est un paramètre-clé permettant de considérer l’effet de la montée de pression d’eau sur la dégradation du module de cisaillement. L’introduction de ce paramètre est pratique car il permet d’ajuster l’influence de l’augmentation de la pression d’eau dans le comportement en cisaillement des sols.

Lors de la thèse de Kteich (2018), une étude paramétrique sur le cas réel de la ville d’Urayasu a été menée pour calibrer ce facteur en comparant les observations de liquéfaction dans la ville aux prédictions pour différentes valeurs de \({\chi }_{\mathit{ru}}\) . La valeur permettant d’obtenir toutes les observations d’occurrence de liquéfaction est égale à 0,6667.

Vérification et validation

Le lecteur se reportera à la documentation d’utilisation de la commande DEFI_SOL_EQUI [U4.84.31] ou aux cas-tests ZZZZ412 (Validation de la commande DEFI_SOL_EQUI), SDLS128 (Sollicitation sismique d’une colonne de sol avec le modèle linéaire équivalent, Comparaison avec SHAKE en déconvolution), SDLS141 (Sollicitation sismique d’une colonne de sol – Comparaison Benchmark PRENOLIN – Comparaison avec le modèle de Hujeux en faible déformation), SDLX101e,h (Vérification du chaînage MISS3D‑code_aster dans le cas d’un bâtiment de grande taille), pour avoir des cas d’applications.

Bibliographie

1. Idriss, I.M., Seed, H.B., 1967. «Response of Horizontal Soil Layers During Earthquakes» Soil Mechanics and Bituminous Materials Research Laboratory, University of California, Berkeley, CA.

2. Schnabel, P.B., Lysmer, J., Seed, H.B., 1972. SHAKE: A Computer Program for Earth- quake Ground Response Analysis for Horizontally Layered Sites. Report No: EERC 72-12. Earthquake Engineering Research Center, University of California, Berkeley, CA.

3. Vandeputte, D., Dévesa G., 2014, » Réponse sismique d’une colonne de sol avec Code_Aster Méthode linéaire équivalente «  - CEIDRE TEGG – Présentation à la journée utilisateur Code_Aster 18 Mars 2014 https://www.code-aster.org/V2/UPLOAD/DOC/Presentation/10-vandeputte-expose.pdf.

4. Yoshida, N., Iai, S., « Nonlinear site response and its evaluation and prediction », Proc. 2ndInternational Symposium on the Effect of Surface Geology on Seismic Motion, Yokosuaka, Japan, pp. 71-90, 1998.

5. Bardet, J.P., Ichii, I, « EERA A computer Program for Equivalent-linear Earthquake site Response Analyses of Layered Soil Deposits » University of Southern California, 2000.

6. Clouteau, D., « Rapport CRSA/EDF-DRD Contrat T62/F28242 Couplage ProMiss-Cyberquake-Code_Aster ».

7. Kramer, S.L., 1996. Prentice Hall, Upper Saddle River, New Jersey. Geotechnical Earthquake Engineering.

8. Schnabel, P.B., Lysmer, J., Seed, H.B., 2012. Computer Program SHAKE: A Computer Program for Earthquake Response Analysis of Horizontally Layered Sites. University of California, Berkeley, California.

9. Hashash, Y.M.A, Groholski, D.R., Phillips, C. A., Park, D, Musgrove, M., 2012. “DEEPSOIL 5.1, User Manual and Tutorial.” 107 p.

10. Idriss, I.M., Sun, J.I., 1992. User’s Manual for SHAKE91: A Computer Program for Conducting Equivalent Linear Seismic Response Analyses of Horizontally Layered Soil Deposits.

11. Dickenson, S.E., Ph.D. Dissertation 1994. Dynamic Response of Soft and Deep Cohesive Soils During the Loma Prieta Earthquake of October, 17, 1989. University of California, Berkeley, CA.

12. Yoshida, N., Kobayashi, S., Suetomi, I., Muira, K., 2002. Equivalent linear method considering frequency dependent characteristics of stiffness and damping. Soil Dyn. Earthquake Eng. 22, 205–222.

13. Assimaki, D., Kausel, E., Whittle, A., 2000. Model for dynamic shear modulus and damping for granular soils. ASCE J. Geotech. Geoenviron. Eng. 126, 859–869.

14. Furumoto, Y., Sugito, M., Yashima. A. Frequency dependent equivalent linearized technique for FEM response analysis of ground. The 12th World Conference on Earthquake Engineering, 2000.

15. Kwak, D-Y-K., Jeong, C-G., Comparaison of frequency dependent equivalent linear analysis methods, The 14th World Conference on Earthquake Engineering, 2008.

16. Kteich, Z., Labbé, P., Javelaud, E., Bennabi, A., Semblat, J.-F. Coupled approach for estimating seismically induced liquefaction and settlements in Urayasu City during the 2011 Great East Japan Earthquake, Soils and Foundations, in revision.

17. Byrne, Peter M., 1991. A Cyclic Shear-Volume Coupling and Pore Pressure Model for Sand. International Conferences on Recent Advances in Geotechnical Earthquake Engineering and Soil Dynamics.

18. Martin, G.R, Finn, W.D.L. and Seed, H.B., 1975. Fundamentals of liquefaction under cyclic loading. Journal of the Geotechnical Engineering Division, ASCE, 101(GT5): 423-438.

19. Wu, G., 1996. Volume change and residual pore water pressure of saturated granular soils to blast loads. A research report submitted to the Natural Sciences and Engineering Research Council of Canada. AGRA Earth & Environmental Ltd., Burnaby

20. Wu, G., 2001. Earthquake induced deformation analyses of the Upper San Fernando dam under the 1971 San Fernando earthquake. Canadian Geotechnical Journal, 38: 1-15