r4.05.07 Interaction sol-structure non-linéaire avec la méthode Laplace-Temps

Résumé:

Ce document est une notice théorique décrivant la méthode Laplace-Temps ( cf . [U2.06.05]) développée dans l’opérateur CALC_MISS ( cf . [U7.03.12]) pour le type de résultat “FICHIER_TEMPS”. Elle s’inscrit dans le cadre des calculs d’interaction sol-structure (ISS) non-linéaires s’appuyant sur le chaînage Code_Aster -MISS3D.

Le principe de résolution d’un tel calcul non-linéaire repose sur une technique de sous-structuration où deux sous-domaines sont généralement considérés: d’une part, un domaine borné non-linéaire correspondant essentiellement à la partie structure (et éventuellement, le sol entourant la fondation) et de l’autre, un domaine non-bornée de sol à comportement linéaire élastique. Cette décomposition de domaines permet d’appliquer la méthode d’éléments finis ( Code_Aster ) pour la modélisation transitoire du problème non-linéaire et une méthode d’éléments de frontière (MISS3D) pour le calcul de la matrice d’impédance représentant le comportement dynamique en temps du domaine de sol semi-infini.

Dans ce cadre, la méthode Laplace-Temps permet de calculer l’opérateur d’impédance en temps à partir d’un échantillonnage particulier de l’espace de fréquences complexes (domaine de Laplace). En particulier, cet échantillonnage est basé sur le polynôme caractéristique d’une méthode multi-pas linéaire de deuxième ordre. De même, la méthode Laplace-Temps fait possible l’évaluation à chaque instant du calcul sismique, de l’intégrale de convolution traduisant les efforts d’interaction sol-structure du système couplé.

La méthode Laplace-Temps

Principe général

Considérons le cas particulier de deux sous-domaines séparés par une interface \(\Gamma ` : l’un noté :math:`{\Omega}_{1}\) , qui est non-borné et qui présente un comportement linéaire et l’autre, \({\Omega}_{2}\) , borné et éventuellement non-linéaire. Les effets d’interaction entre les deux sous-domaines peuvent être représentés sur l’interface par une impédance que nous supposons définie dans le domaine de Fourier ou de Laplace. Dans ce cadre, on résout le problème global en temps dans \({\Omega}_{2}\) en prenant en compte les effets du sous-domaine \({\Omega}_{1}\) au moyen d’une force externe appliquée sur la frontière. Cette force d’interaction, qui met en jeu la fonction impédance dans le domaine temporel \(Z(t)\) correspond à une intégrale de convolution avec le champ inconnu \({u}_{\Gamma}(t)\) , que l’on notera \((Z\ast {u}_{\Gamma})(t)\) .

Ainsi, de manière générale, l’évaluation des efforts d’interaction entre les sous-domaines passe par le calcul d’une intégrale de convolution entre deux fonctions causales définie comme:

(1265)\[(Z\ast {u}_{\Gamma})(t)={\int}_{0}^{t}Z(t-\tau ){u}_{\Gamma}(\tau )\text{}d\tau ,\text{}0\le t\le T\]

Soit \(\widehat{Z}(s)\) la transformée de Laplace de \(Z(t)\) , on supposera qu’elle est analytique dans le demi-plan complexe \(\mathit{\Re e}(s)>{\sigma}_{0}\) et à croissance lente pour \(\mid s\text{|}\) grand:

:math:`

widehat{Z}(s)

le C({sigma}_{0})

s{

}^{mu}text{avec}C({sigma}_{0})text{et}mu in ℝ`

Dans ce cadre, soit \(\widehat{Z}(s)={\widehat{Z}}_{m}(s)\widehat{P}(s)\) avec \(\widehat{P}(s)\) une fonction polynomiale en \(s\) à valeurs matricielles [1]_ de degré \(m\ge \mu\) :

(1266)\[\widehat{P}(s)=\sum_{p=0}^{m}{\Lambda}_{p}{s}^{p}\]

Si cette décomposition polynomiale est considérée, le noyau de convolution \(Z(t)\) défini au sens des distributions peut s’exprimer sous la forme d’un opérateur de différentiation d’ordre \(p\) comme:

(1267)\[Z(t)=\left({Z}_{m}\ast \sum_{p=0}^{m}{\Lambda}_{p}\frac{{d}^{p}\delta }{d{t}^{p}}\right)(t)\]

Par conséquent, le produit de convolution peut finalement s’écrire comme :

(1268)\[(Z\ast {u}_{\Gamma})(t)=\sum_{p=0}^{m}({\int}_{0}^{t}{Z}_{m}(\tau ){\Lambda}_{p}{u}_{\Gamma}^{(p)}(t-\tau )d\tau )\]

avec \({u}_{\Gamma}(t)\) une fonction causale suffisamment différentiable.

En se servant de la méthode de quadratures de Lubich ,on peut exprimer l’équation précédente sous forme d’une convolution discrète. En effet, si l’on note \(\Delta t\) le pas de temps de discrétisation, celle-ci peut s’écrire aux instants de temps \(n\Delta t\) (\(0\le n\Delta t\le t\) ) comme suit :

(1269)\[(Z\ast {u}_{\Gamma})(n\Delta t)=\sum_{k=1}^{n}({\Psi}_{1}^{n-k+1}{u}_{\Gamma ,k}+...+{\Psi}_{p}^{n-k+1}{u}_{\Gamma ,k}^{(p)}+...+{\Psi}_{m}^{n-k+1}{u}_{\Gamma ,k}^{(m)})\]

où les coefficients \({\Psi}_{k}^{j}\) contiennent la contribution des matrices \({\Lambda}_{k}\) .

Application aux opérateurs d’impédance de sol

La matrice d’impédance de sol, c’est-à-dire la version discrétisée en espace de l’opérateur d’impédance, relie un vecteur déplacement défini sur les degrés de liberté de l’interface sol-structure (ou de façon généralisée dans le cadre du chaînage Aster-MISS3D, sur la frontière FEM-BEM) au vecteur force défini sur la même interface, notée :math:`Gamma ` dans ce qui suit:

(1270)\[\widehat{Z}(s){\widehat{u}}_{\Gamma}(s)={\widehat{F}}_{\Gamma}(s)\]

De manière générale, la matrice d’impédance de sol, qui est supposée analytique dans le demi-espace \(\mathit{\Re e}(s)>{\sigma}_{0}\) pour des raisons de causalité, peut s’exprimer sous la forme suivante :

(1271)\[\widehat{Z}(s)={\widehat{Z}}_{\mathit{sing}}(s)+{\widehat{Z}}_{\mathit{nsing}}(s)\]

où \({\widehat{Z}}_{\mathit{nsing}}(s)\) correspond à la partie régulière de l’impédance dont la transformée inverse de Laplace \({Z}_{\mathit{nsing}}(t)\) , qui existe au sens classique, est définie comme une fonction exponentiellement bornée et continue en temps qui s’annule pour \(t<0\) :

(1272)\[{Z}_{\mathit{nsing}}(t)=\frac{1}{2\pi i}{\int}_{{\sigma}_{0}+iR}{e}^{st}{\widehat{Z}}_{\mathit{nsing}}(s)\mathit{ds}\]

En revanche, \({\widehat{Z}}_{\mathit{sing}}(s)\) désigne la partie singulière de l’impédance, entendue ici comme une fonction à croissance lente et non-bornée en haute fréquence. En particulier, par analogie à une formulation FEM , cette partie singulière de l’impédance peut s’écrire comme suit :

(1273)\[{\widehat{Z}}_{\mathit{sing}}(s)={M}_{\Gamma}{s}^{2}+{C}_{\Gamma}s+{K}_{\Gamma}\]

où \({M}_{\Gamma}\) , \({C}_{\Gamma}\) et \({K}_{\Gamma}\) modélisent respectivement les contributions du sol en termes d’inertie, amortissement et rigidité au niveau de l’interface \(\Gamma ` . On met en évidence deux points intéressants à mentionner. Le premier vise à souligner que cette formulation permet de montrer facilement que l'impédance vérifie explicitement la condition explicité dans l'équation () pour :math:`m=2\) . Le deuxième porte plutôt sur le caractère non-singulier de la transformée inverse de Laplace qui cette fois-ci n’existe qu’au sens des distributions et qui par conséquent, s’écrit au moyen du delta de Dirac et de ses dérivées de premier et deuxième ordre:

(1274)\[{Z}_{\mathit{sing}}(t)={M}_{\Gamma}\ddot{\delta}(t)+{C}_{\Gamma}\dot{\delta}(t)+{K}_{\Gamma}\delta (t)\]

La manipulation numérique de \({\widehat{Z}}_{\mathit{sing}}(t)\) est très délicate et en même temps indispensable, car cette évolution doit être calculée pour l’obtention des efforts transitoires d’interaction sol-structure. Par conséquent, on peut affirmer que le caractère particulier de l’impédance temporelle fait nécessaire l’utilisation d’approches de discrétisation en temps plus efficaces et parmi celles-ci on trouve la méthode Laplace-Temps reposant sur une méthode de quadratures de convolution.

Dans ce cadre, la méthode de quadrature de convolution présentée par Lubich peut être appliquée pour l’évaluation numérique d’une intégrale de convolution telle que donnée dans l’équation (), où on considèrepour \(Z(t)\) la transformée de Laplace inverse de l’impédance de sol \(\widehat{Z}(s)\) et \({u}_{\Gamma}(t)\) le champ déplacement sur l’interface. Cette méthode permet d’évaluer de manière approximative l’équation () par une convolution discrète (avec un pas de temps \(\Delta t>0\) ):

(1275)\[(Z\ast {u}_{\Gamma})(n\Delta t)=\sum_{0\le n\Delta t\le t}{\Phi}_{k}{u}_{\Gamma}(t-n\Delta t)\]

où les coefficients \({\Phi}_{k}\) correspondent aux poids de la série de puissances suivante :

(1276)\[\sum_{k=0}^{+\infty }{\Phi}_{k}{\zeta}^{k}=\widehat{Z}({s}_{\Delta t})\]

Les points d’échantillonnage \({s}_{\Delta t}\) de l’impédance dynamique de sol sont données dans la section suivante.

Arrivés à ce point et en tenant compte que l’expression () est homogène à une force, il paraît judicieux d’exprimer cette convolution non seulement en termes de déplacement, mais aussi en fonction des accélérations et des vitesses. L’approche proposée s’inscrit dans cette optique et repose sur la factorisation de la partie polynomiale \(\widehat{P}(s)\) de l’impédance:

(1277)\[\widehat{Z}(s)={\widehat{Z}}_{m}(s)\widehat{P}(s)={\widehat{Z}}_{m}(s)({\stackrel{~}{M}}_{\Gamma}{s}^{2}+{\stackrel{~}{C}}_{\Gamma}s+{\stackrel{~}{K}}_{\Gamma})\]

Où \({\stackrel{~}{K}}_{\Gamma}\) , \({\stackrel{~}{C}}_{\Gamma}\) et \({\stackrel{~}{M}}_{\Gamma}\) sont respectivement des estimateurs des matrices de rigidité, d’amortissement et d’inertie du sol. Ainsi, la convolution peut s’écrire au moyen de la transformée de Laplace inverse comme suit :

(1278)\[(Z\ast {u}_{\Gamma})(t)=\frac{1}{2\pi i}{\int}_{{\sigma}_{0}+iR}{\widehat{Z}}_{m}(s)\widehat{P}(s){\widehat{u}}_{\Gamma}(s){e}^{st}\mathit{ds}\]

Par conséquent, la fonction polynomiale \(\widehat{P}(s)\) étant vu comme un opérateur de différentiation qui agit sur le champ déplacement, l’équation () devient:

(1279)\[(Z\ast {u}_{\Gamma})(t)=({Z}_{m}\ast {\stackrel{~}{M}}_{\Gamma}{\ddot{u}}_{\Gamma})(t)+({Z}_{m}\ast {\stackrel{~}{C}}_{\Gamma}{\dot{u}}_{\Gamma})(t)+({Z}_{m}\ast {\stackrel{~}{K}}_{\Gamma}{u}_{\Gamma})(t)\]

où les efforts d’interaction (notés dans la suite par \({R}_{\Gamma}(t)\) ) sont calculés par des convolutions avec les accélérations, vitesses et déplacements de l’interface.

En prenant un pas de temps de discrétisation \(\Delta t>0\) , l’intégrale de convolution devient un cas particulier de l’équation () pour \(m=2\) :

(1280)\[{R}_{n}=(Z\ast {u}_{\Gamma})(n\Delta t)=\sum_{k=1}^{n}({\Psi}_{2}^{n-k+1}{\ddot{u}}_{\Gamma ,k}+{\Psi}_{1}^{n-k+1}{\dot{u}}_{\Gamma ,k}+{\Psi}_{0}^{n-k+1}{u}_{\Gamma ,k})\]

où les matrices multipliant les vecteurs de déplacement, vitesse et accélération sont donnés par :

(1281)\[{\Psi}_{0}^{k}={Z}_{m}^{k}{\stackrel{~}{K}}_{\Gamma} {\Psi}_{1}^{k}={Z}_{m}^{k}{\stackrel{~}{C}}_{\Gamma}\]

En effet, de manière analogue au raisonnement suivi pour les équations () et (), les coefficients des matrices \({\Psi}_{0}^{k}\) , \({\Psi}_{1}^{k}\) et \({\Psi}_{2}^{k}\) correspondent aux poids des séries de puissances suivantes:

(1282)\[\sum_{k=0}^{+\infty }{\Psi}_{j}^{k}{z}^{k}={\widehat{Z}}_{m}({s}_{\Delta t}){\Lambda}_{j}\]

Avec \({\Lambda}_{0}={\stackrel{~}{K}}_{\Gamma}\) , \({\Lambda}_{1}={\stackrel{~}{C}}_{\Gamma}\) et \({\Lambda}_{2}={\stackrel{~}{M}}_{\Gamma}\) .

[1]

Pour simplifier la démarche théorique, nous considérons dans la suite toutes les grandeurs déjà discrétisées dans l’espace (par exemple, avec une méthode éléments finis).

Calcul de l’impédance en temps

Comme mentionné précédemment, la méthode Laplace-Temps est une approche numérique permettant de discrétiser en temps une intégrale de convolution.

Implantation dans code_aster

L’équation () montre que l’évaluation numérique des poids \({\Psi}_{k}^{j}\) passe d’abord par celle de \({\widehat{Z}}_{m}^{k}\) . A son tour, celle-ci peut être obtenue à partir d’une intégrale de Cauchy appliquée sur un contour :math:`|z|=rho ` :

:math:`{Z}_{m}(kDelta t)=frac{1}{2pi i}{int}_{

z

=rho }{widehat{Z}}_{m}left(frac{delta (z)}{Delta t}right){z}^{-k-1}dz`

Si on exprime cette intégrale en coordonnées polaires \(z=\rho {e}^{i\theta }\) , puis on applique la règle de trapèzes à la phase pour la discrétiser en \(L\) pas identiques de valeur \(\Delta \theta =\frac{2\pi }{L}\) , on obtient:

(1283)\[{Z}_{m}(k∆t)\approx {Z}_{m}^{k}=\frac{{\rho}^{-k}}{L}\sum_{l=0}^{L-1}{\widehat{Z}}_{m}({s}_{l}){e}^{-i\frac{2\pi l}{L}k}\]

où l’opérateur défini dans le domaine de Laplace est échantillonné sur \({s}_{l}=\frac{\delta (\rho {e}^{2\pi il∕L})}{\Delta t}\) avec \(\delta (z)=\frac{3}{2}-2z+\frac{1}{2}{z}^{2}\) le polynôme caractéristique d’une méthode linéaire multi-pas d’ordre deux. Notons que \(N\) et \(L\) ne sont pas forcement égaux, \(N\) étant le nombre de pas de temps de la fenêtre de temps d’intérêt (par exemple, en fonction de la durée du signal sismique) et \(L\) étant le nombre de pas de temps de calcul de l’impédance en temps.

Si on considère que \({\widehat{Z}}_{m}({s}_{l})\) est évalué avec une variabilité [2]_ \({ϵ}_{\mathit{CQM}}\) , les coefficients \({\widehat{Z}}_{m}^{k}\) seront calculés avec une précision [3]

\(O(\sqrt{{ϵ}_{\mathit{CQM}}})\) s i \(L=N\) et \({\rho}^{N}=\sqrt{{ϵ}_{\mathit{CQM}}}\) . Avec les valeurs de \(L\) et \(\rho\) fixées, les poids de la série peuvent être obtenus à partir d’une FFT classique, ce qui ramène la complexité de l’algorithme de calcul à \(O(L\logL)\) au lieu de \(O({L}^{2})\) . Néanmoins, la valeur de \({ϵ}_{\mathit{CQM}}\) n’est pas évidente à calibrer et dépend en partie de la nature du problème à résoudre. En effet, des valeurs trop faibles de \({ϵ}_{\mathit{CQM}}\) , telles que \({10}^{-20}\) ou \({10}^{-30}\) peuvent entraîner des rayons d’intégration \(\mid z\mid =\rho\) trop petits. Vu que l’intégrande de Cauchy contient une singularité en zéro, des valeurs trop petites du rayon peuvent se traduire en une mauvaise évaluation de l’intégrale. A ce titre, des études paramétriques ont été réalisés dans pour conclure qu’un bon compromis est atteint lorsque \({ϵ}_{\mathit{CQM}}={10}^{-10}\) [4]

.

De même, on remarquera que la dépendance en \({\rho}^{-k}\) mène, pour des valeurs de \(k\) grandes, à des imprécisions sur les coefficients \({\widehat{Z}}_{m}^{k}\) calculés. Sachant que l’utilisation des méthodes d’intégration plus précises, de type Simpson par exemple, ne corrigent pas ce problème, l’idée est d’élargir la plage de temps (i.e. \(L=\mathit{mN},m\in \mathrm{ℝ}\) ) sur laquelle on calcule l’impédance afin de réduire les perturbations sur la solution d’intérêt. A ce titre, on a observé que, en général, les imprécisions de \({\widehat{Z}}_{m}^{k}\) ne se propagent sur la réponse qu’au-delà de \(t\approx 0,7{T}_{\mathit{FIN}}\) , où \({T}_{\mathit{FIN}}\) est l’instant final de calcul de l’impédance [5] .

Considérations numériques principales

Nous avons vu précédemment qu’en générale, la partie singulière de l’impédance de sol peut être approchée par un polynôme d’ordre deux. Si les coefficients de ce polynôme sont réels, il est facile à démontrer que l’impédance de sol vérifie \(\widehat{Z}(\sigma -i\omega )=\mathit{conj}[\widehat{Z}(\sigma +i\omega )]\) dans le plan analytique du domaine de Laplace [6] . On remarquera que lorsque un amortissement hystérétique est utilisé dans le sol (cas de MISS3D), la partie imaginaire de sa raideur statique est non-nulle et l’hypothèse de coefficients réels n’est plus vérifiée.

Une deuxième propriété intéressante venant directement de la définition de transformé unilatérale de Laplace est la suivante:

(1284)\[\widehat{Z}(s=0)={\int}_{0}^{ꝏ}Z(t)\mathit{dt}\approx \sum_{k}^{N}Z(k\Delta t)\]

pour \(N\) choisi suffisamment grand (~ une centaine de pas de temps). Notons que l’impédance de sol évaluée à \(s=0\) correspond à la raideur statique de sol (partie réelle) qui elle, est très différente de la raideur instantanée du sol. Ceci doit être bien pris en compte lorsque un calcul dynamique (\(Z(t=0)]\) est assemblé au premier membre) est réalisé à la suite d’un calcul statique (\(\mathrm{\Re }[\widehat{Z}(s=0)]\) est assemblé au premier membre) [7]_ .

[1]

Faisant référence à la notion de reproductibilité (precision en anglais).

[2]

Faisant référence à la notion d’exactitude (accuracy en anglais).

[3]

Cette valeur initialise par défaut le mot-clé PRECISION de l’opérateur CALC_MISS (cf. [U7.03.12]) lorsque le type de résultat attendu est “FICHIER_TEMPS”.

[4]

Pour ce faire, CALC_MISS (TYPE_RESU = “FICHIER_TEMPS”) prévoit le mot-clé COEF_SURECH (cf. [U7.03.12]) initialisé à 1.35, soit l’équivalent de faire:math:`L=mathrm{1.35N}`dans l’équation 23.

[5]

Cette propriété est exploitée dans code_aster : l’impédance est échantillonnée sur le demi-plan supérieur, puis complétée numériquement avec le conjugué. Par conséquent, il faut faire attention avec les sols présentant des amortissements hystérétiques très forts (supérieur à 25%).

[6]

Il existe dans code_aster une macro-commande (PRE_SEISME_NONL, cf. [U4.63.02])) permettant de réaliser la transition statique-dynamique d’un calcul d’interaction sol-structure non-linéaire.

Analyse sismique

Le calcul sismique d’interaction sol-structure non-linéaire doit être résolu dans le domaine temporel. Quelques éléments sur le calcul de la force sismique et la discrétisation en temps des équations du système couplé sont fournis ci-dessous.

Calcul de la force sismique en temps

La méthode Laplace-Temps pourrait être appliquée aussi pour le calcul de la force sismique, néanmoins la nature du champ incident développé dans MISS3D invalide les hypothèses sur les champs solution recherchées dans le domaine de Laplace. Dans ce cadre, on ne peut pas garantir que le calcul de la force sismique à fréquence complexe avec MISS3D soit correct pour tout type de fondation (superficielle et enterrée) et stratigraphie. L’approche préconisée dans les études passe ainsi par la transformée de Fourier inverse de la force sismique évaluée par MISS3D dans le domaine fréquentiel.

Schéma d’intégration en temps

Cette partie cherche à illustrer le traitement numérique de l’intégrale de convolution dans le cadre d’une résolution avec un schéma d’intégration en temps de la famille de Newmark. Néanmoins, le raisonnement s’applique aussi au schéma d’intégration en temps :math:`alpha ` -HHT, disponible également dans code_aster.

On cherchera ainsi à résoudre en temps le problème dynamique linéaire suivant, supposé discrétisé en espace avec la méthode classique des éléments finis:

(1285)\[\begin{split}\left[\begin{array}{cc}{L}_{bb}(\cdot )& {L}_{b\Gamma }(\cdot )\\ (sym.)& {L}_{\Gamma \Gamma }(\cdot )\end{array}\right]\left[\begin{array}{c}{u}_{b}(t)\\ {u}_{\Gamma}(t)\end{array}\right]+\left[\begin{array}{c}0\\ {R}_{\Gamma}(t)\end{array}\right]=\left[\begin{array}{c}{F}_{b}(t)\\ {F}_{\Gamma}(t)\end{array}\right]\end{split}\]

où un opérateur de différentiation \({L}_{\mathit{\alpha \beta }}(\cdot )\) en temps est introduit pour \(\alpha ,\beta \in \lbrace b,\Gamma \rbrace\) , \(\Gamma\) pour les degrés de liberté de l’interface en contraposition au reste des degrés de liberté du système notés \(b\) et où le vecteur d’efforts d’interaction sol-structure noté \({R}_{\Gamma}(t)\) , correspond à l’intégrale de convolution qu’on discrétisera avec la méthode Laplace-Temps.

La résolution du système couplé sol-structure passe par la discrétisation de toutes les grandeurs en temps et parisoler les inconnues à l’instant \(t=n\Delta t\) , où le pas de temps est noté \(\Delta t\) . En particulier, pour les efforts d’interaction sol-structure, il est intéressant de partir de l’équation () et de regrouper les termes non-inconnus correspondant aux instants précédents à \(t=n\Delta t\) dans un seul terme afin de les séparer des inconnues à résoudre au pas de temps \(n\) . Le vecteur \({R}_{n}\approx {R}_{\Gamma}(n\Delta t)\) correspondant s’écrit ainsi comme suit :

(1286)\[{R}_{n}={\Psi}_{2}^{1}{\ddot{u}}_{\Gamma ,n}+{\Psi}_{1}^{1}{\dot{u}}_{\Gamma ,n}+{\Psi}_{0}^{1}{u}_{\Gamma ,n}+{R}_{\Sigma (n-1)}\]

Avec \({\Psi}_{0}^{1}\) , \({\Psi}_{1}^{1}\) et \({\Psi}_{2}^{1}\) faisant référence respectivement à la raideur, l’amortissement et l’inertie instantanés [8]_ et \({R}_{\Sigma (n-1)}\) dépendant uniquement des quantités calculées aux instants précédents [9]_ :

(1287)\[{R}_{\Sigma (n-1)}=\sum_{k=1}^{n-1}({\Psi}_{2}^{n-k+1}{\ddot{u}}_{\Gamma ,k}+{\Psi}_{1}^{n-k+1}{\dot{u}}_{\Gamma ,k}+{\Psi}_{0}^{n-k+1}{u}_{\Gamma ,k})\]

Ce regroupement de termes permet d’introduire au second membre la partie de la convolution venant des instants précédents et dans le premier membre les termes portant sur des inconnues. En effet, ceci peut facilement mis en évidence en partant d’une formulation en déplacement d’un schéma d’intégration Newmark inconditionnellement stable (\(\beta =0.25\) , \(\gamma =0.5\) ) appliquée à l’équation () :

(1288)\[\begin{split}\left[\begin{array}{cc}{\stackrel{̃}{K}}_{11}& {\stackrel{̃}{K}}_{12}\\ {\stackrel{̃}{K}}_{21}& {\stackrel{̃}{K}}_{22}\end{array}\right]\left[\begin{array}{c}{u}_{b,n}\\ {u}_{\Gamma ,n}\end{array}\right]=\left[\begin{array}{c}{F}_{b,n}\\ {F}_{\Gamma ,n}\end{array}\right]-\left[\begin{array}{c}0\\ {R}_{\Sigma (n-1)}\end{array}\right]\end{split}\]

où la matrice \(\stackrel{̃}{K}\) correspond à l’opérateur de Newmark qu’il faudra inverser dans chaque pas de temps et le second membre à la somme du vecteur équivalent de Newmark et la partie connue du produit de convolution discret. En particulier, le terme \(\stackrel{̃}{{K}_{22}}\) qui s’applique aux degrés de liberté de l’interface \(\Gamma\) , s’écrit :

(1289)\[{\stackrel{̃}{K}}_{22}=\frac{1}{\beta \Delta {t}^{2}}({M}_{22}+{\Psi}_{2}^{1})+\frac{\gamma}{\beta \Delta t}({C}_{22}+{\Psi}_{1}^{1})+({K}_{22}+{\Psi}_{0}^{1})\]

On remarquera que les termes \({\Psi}_{0}^{1}\) , \({\Psi}_{1}^{1}\) et \({\Psi}_{2}^{1}\) jouent un rôle sur le conditionnement de l’opérateur de Newmark \(\stackrel{̃}{K}\) et par conséquent, leur valeur ont un impact sur la convergence de la solution. Notons également que \(\stackrel{~}{{K}_{22}}\) , restant constant pendant le calcul, on ne le calcule qu’une fois.

[1]

Dans code_aster, ces termes instantanés sont assemblés au premier membre via l’opérateur MACR_ELEM_DYNA (cf. [U4.65.01]) et éventuellement une liaison LIAISON_INTERF dans l’opérateur AFFE_CHAR_MECA (cf. [U4.44.01]) si une base modale réduite d’interface a été utilisée.

[2]

Le terme:math:`{R}_{Sigma (n-1)}`est modélisé dans code_aster lors de l’analyse dynamique non-linéaire comme une charge au second membre obtenue avec l’option « CHARGE_SOL » de l’opérateur AFFE_CHAR_MECA (cf. [U4.44.01]).

Bibliographie

[1] G.R. Darbre and J.P. Wolf, Criterion of stability and implementation issues of hybrid frequency time domain procedure for non-linear dynamic analysis, Earthquake Engineering and Structural Dynamics, 16, 4, 569-581, 1988.

[2] R. Cottereau, D. Clouteau, C. Soize and S. Cambier, Probabilistic nonparametric models of impedance matrices. Application to the seismic design of a structure., European Journal of Computational Mechanics, 15, 1-3, 131-142, 2006.

[3] J.P. Wolf, Soil-Structure-Interaction Analysis in Time Domain, Prentice-Hall, New Jersey, USA, 1988.

[4] S. François and G. Degrande, A time domain coupled boundary element-finite element method for the dynamic response of structures, Proceedings of the 12th International Congress on Sound and Vibration, 2005.

[5] L. Gaul and M. Schanz, A comparative study of three boundary element approaches to calculate the transient response of viscoelastic solids with unbounded domains, Computer Methods in Applied Mechanics and Engineering, 179, 111-123, 1999.

[6] M. Schanz and H. Antes, Aplication of ’Operational Quadrature Methods’ in Time Domain Boundary Element Methods, Meccanica, 32, 3, 179-186, 2006.

[7] A. Pereira and G. Beer, Interface dynamic stiffness matrix approach for three-dimensional transient multi-region boundary element analysis, International Journal for Numerical Methods in Engineering, 80, 1463-1495, 2009.

[8] C. Lubich, Convolution quadrature and discretized operational calculus I, Numerische Mathematik, 52, 129-145, 413-425, 1988.

[9] C. Lubich, Convolution quadrature and discretized operational calculus II, Numerische Mathematik, 52, 413-425,1988.

[10] Nieto Ferro, A., Nonlinear Soil-Structure Interaction in Earthquake Engineering, PhD Thesis, Ecole Centrale Paris, 2013.

[11] W. Moser and H. Antes and G. Beer, A Duhamel integral based approach to one-dimensional wave propagation analysis in layered media, Computational Mechanics, 35, 115-126, 2005.

[12] W. Moser and H. Antes and G. Beer, Soil-structure interaction and wave propagation problems in 2D by a Duhamel integral based approach and the convolution quadrature method, Computational Mechanics, 36, 431-443, 2005.