r4.03.05 Modèles probabilistes paramétriques et non-paramétriques en dynamique

Résumé:

Ce document décrit deux approches probabilistes, l’une paramétrique et l’autre non paramétrique afin de prendre en compte les incertitudes de modèle et de modélisation pour les systèmes dynamiques en mécanique des structures. L’approche non-paramétrique est spécifique à la résolution sur base modale des systèmes linéaires (en transitoire ou en harmonique) et à la résolution sur base modale des systèmes avec des non‑linéarités localisées (en transitoire). L’approche paramétrique est a priori universelle, mais elle est plus particulièrement présentée dans le cas où elle est combinée avec l’approche non paramétrique.

Modélisations du système dynamique

Modèle éléments finis moyen

Résolution transitoire en coordonnées absolues

Dans le repère absolu, le système mécanique est modélisé par la méthode aux éléments finis. Ce modèle de base (en général celui qui aurait été utilisé dans l’étude déterministe) est désigné sous le nom de «modèle éléments finis moyen». Toutes les grandeurs relatives aux modèles moyens sont soulignées.

Soit \(t\mapsto \underline{\mathrm{y}}(t)\) la réponse transitoire dans le repère absolu du «modèle éléments finis moyen» définie sur l’intervalle d’étude \([0,T]\) et à valeur dans \({ℝ}^{k}\) où \(k\) est le nombre de degrés de liberté. Les matrices de masse, d’amortissement et de rigiditésont respectivement notées \(\left[\underline{M}\right]\) , \(\left[\underline{D}\right]\) et \(\left[\underline{K}\right]\) .

La réponse transitoire \(\underline{y}(t)\) du «modèle éléments finis moyen» vérifie l’équation différentielle non linéaire discrétisée suivante:

\(\left[\underline{M}\right]\underline{\ddot{y}}(t)+\left[\underline{D}\right]\underline{\dot{y}}+{f}_{c}(t,\underline{\dot{y}(t)},\underline{y(t)},\underline{w})=f(t)\) , t  , éq 2.1.1-1

avec les conditions initiales,

\(\underline{y}(0)=\dot{\underline{y}}(0)=0\) , éq 2.1.1-2

• \(f(t)\in {ℝ}^{m}\) représente la discrétisation par éléments finis des forces extérieures.

• \({f}_{c}(t,\underline{y}(t),\underline{\dot{y}}(t),\underline{w})\in {ℝ}^{m}\) correspond aux non linéarités localisées (par exemple dues à des butées élastiques de choc). Les éléments \({\underline{w}}_{1,}\mathrm{...},{\underline{w}}_{v}\) du vecteur \(\underline{w}\in {ℝ}^{\nu}\) représentent un jeu de paramètres définissant ces non linéarités ( par exemple jeu, raideur de choc, amortissement de choc, etc.).

Résolution transitoire en coordonnées relatives (séisme)

Comme dans le cas transitoire en coordonnées absolues, le système mécanique est modélisé par un modèle de base, le «modèle éléments finis moyen».

On note \(t\to \underline{z}(t)\) la réponse transitoire en coordonnées absolues de ce modèle sur l’intervalle d’étude \([0,T]\) à valeur dans \({ℝ}^{k}\) ( attention: remarquez le changement de notation par rapport au paragraphe précédent.).

La réponse transitoire du «modèle éléments finis moyen» vérifie l’équation différentielle non linéaire discrétisée suivante:

\(\left[\begin{array}{cc}\left[\underline{M}\right]& \left[{\underline{M}}_{\mathrm{ls}}\right]\\ {\left[{\underline{M}}_{\mathrm{ls}}\right]}^{T}& \left[{\underline{M}}_{s}\right]\end{array}\right]\left[\begin{array}{c}\underline{\ddot{z}}(t)\\ {\underline{\ddot{z}}}_{s}(t)\end{array}\right]+\left[\begin{array}{cc}\left[\underline{D}\right]& \left[{\underline{D}}_{\mathrm{ls}}\right]\\ {\left[{\underline{D}}_{\mathrm{ls}}\right]}^{T}& \left[{\underline{D}}_{s}\right]\end{array}\right]\left[\begin{array}{c}\underline{\dot{z}}(t)\\ {\underline{\dot{z}}}_{s}(t)\end{array}\right]+\left[\begin{array}{cc}\left[\underline{K}\right]& \left[{\underline{K}}_{\mathrm{ls}}\right]\\ {\left[{\underline{K}}_{\mathrm{ls}}\right]}^{T}& \left[{\underline{K}}_{s}\right]\end{array}\right]\left[\begin{array}{c}\underline{z}(t)\\ {\underline{z}}_{s}(t)\end{array}\right]+\left[\begin{array}{c}{F}_{c}(t,\underline{z}(t),\underline{\dot{z}}(t);\underline{w})\\ {0}_{d}\end{array}\right]=\left[\begin{array}{c}\underline{g}(t)\\ {\underline{g}}_{s}(t)\end{array}\right]\)

\(t\in [0,T]\) . éq 2.1.2-1

avec les conditions initiales,

\(\underline{z}(0)=\underline{\dot{z}}(0),{\underline{z}}_{s}(0)={\underline{\dot{z}}}_{s}(0)\) éq 2.1.2-2

• \(\underline{g}(t)\in {ℝ}^{m}\) représente la discrétisation par éléments finis des forces extérieures et \(\underline{{g}_{s}}(t)\in {ℝ}^{d}\) correspond à la discrétisation des forces de réaction dues aux \(d\) conditions de Dirichlet .

• \({F}_{c}(t,\underline{z}(t),\underline{\dot{z}}(t),\underline{w})\in {ℝ}^{m}\) correspond aux non linéarités localisées avec comme précédemment \(\underline{w}\in {ℝ}^{\mathrm{\nu }}\) représentant un jeu de paramètres définissant ces non linéarités.

Après relèvement statique, les équations matricielles [éq 2.1.2-1] et [éq 2.1.2-2] dans le repère absolu sont réécrites en coordonnées « relatives » :

\(\left[\underline{M}\right]\underline{\ddot{y}}(t)+\left[\underline{D}\right]\underline{\dot{y}}+\left[\underline{K}\right]\underline{y(t)}+{f}_{c}(t,\underline{\dot{y}(t)},\underline{y(t)},\underline{w})=f(t)\) , \(t\in [0,T]\) . éq 2.1.2-3

\(\underline{y}(0)=\dot{\underline{y}}(0)=0\) , éq 2.1.2-4

• \(\underline{y}(t)\in {ℝ}^{m}\) est le vecteur des degrés de liberté libres dans le système de coordonnées « relatives » tel que

\(\underline{z}(t)=\underline{y}(t)+[\underline{R}]{\underline{z}}_{s}(t)\) avec \([\underline{R}]=-{[\underline{K}]}^{-1}[\underline{{K}_{\mathrm{ls}}}]\)

• La fonction \(t\to \underline{\mathrm{f}}(t)\) définie sur \([0,T]\) et à valeur dans \({ℝ}^{m}\) et l’application non linéaire \((x,y)\mapsto {f}_{c}(t,x,y;\underline{w})\) de \({ℝ}^{m}\times {ℝ}^{m}\) dans \({ℝ}^{m}\) sont telles que :

\(f(t)=\underline{g}(t)-([\underline{M}][\underline{R}]+[\underline{{M}_{\mathrm{ls}}}])\underline{\ddot{z}}(t)-([\underline{D}][\underline{R}]+[\underline{{D}_{\mathrm{ls}}}])\underline{\dot{z}}(t)\) éq 2.1.2-5

\({f}_{c}(t,x,y;\underline{w})={F}_{c}(t,x+\left[\underline{R}\right]{\underline{z}}_{s}(t),y+\left[\underline{R}\right]{\underline{\dot{z}}}_{s}(t);\underline{w}).\) éq 2.1.2-6

Remarques :

• Dans la suite, selon que l’on a effectué ou non un relèvement statique, \(\underline{y}(t)\) correspond soit à la réponse transitoire en coordonnées absolues définie par le [§ 2.1.1], soit la réponse transitoire en coordonnées « relatives » définie par le [§ 2.1.2].

• On suppose que si les \(d\) conditions de Dirichlet étaient homogènes aucun mouvement de corps rigide ne pourrait se produire. Par conséquent, \(\left[\underline{K}\right]\) est symétrique définie positive et son inverse \({\left[\underline{K}\right]}^{-1}\) est définie, ce qui permet d’introduire la matrice réelle \([\underline{R}]=-{[\underline{K}]}^{-1}[\underline{{K}_{\mathrm{ls}}}]\) de dimension \((m\times d)\) .

• Dans Code_Aster le terme d’amortissement dans [éq 2.1.2-5] est négligé.

Résolution harmonique

Comme dans le cas transitoire, le système mécanique est modélisé par un modèle de base, le «modèle éléments finis moyen». Sur une bande fréquentielle \(\left[{\omega}_{1},{\omega}_{2}\right]\) , la réponse harmonique \(q(\omega )\) du «modèle éléments finis moyen» linéaire vérifie l’équation suivante:

\((-{\omega}^{2}\left[\underline{M}\right]+i\omega \left[\underline{D}\right]+\left[\underline{K}\right])q(\omega )=F(\omega )\) , \(\omega \in \left[{\omega}_{1},{\omega}_{2}\right]\) éq 2.1.3-1

avec \(f(\omega )\) représentant la discrétisation par éléments finis des forces extérieures.

Modèle matriciel réduit moyen

On suppose que l’énergie de vibration de la réponse dynamique est principalement localisée dans le domaine des basses fréquences. On peut donc construire le modèle matriciel réduit moyen en projetant \(\underline{y}(t)\) ou \(\underline{y}(\omega )\) sur le sous espace propre engendré par les n premiers modes du système dynamique linéaire (jeux infinis) conservatif homogène (supports bloqués) associé qui s’écrit,

\(\left[\underline{K}\right]\underline{\varphi}=\lambda \left[\underline{M}\right]\varphi\) . éq 2.2-1

Les matrices \(\left[\underline{M}\right]\) et \(\left[\underline{K}\right]\) étant définies positives (pour \(\left[\underline{K}\right]\) cf. remarque 2 [§2.1.2]), les valeurs propres sont réelles et positives,

\(0\le \underline{{\lambda}_{1}}\le \underline{{\lambda}_{2}}\le \mathrm{...}\le {\lambda}_{n}\) éq 2.2-2

Les modes propres de vibration associés :math:`lbrace {phi }_{1,}{phi }_{2,}dots rbrace ` vérifient les propriétés d’orthogonalité,

\(\text{<}\left[\underline{M}\right]{\varphi}_{\alpha},{\varphi}_{\beta}\text{>}=\underline{{\mu}_{\alpha}}{\delta}_{\alpha \beta }\) éq 2.2-3

\(\text{<}\left[\underline{K}\right]{\varphi}_{\alpha},{\varphi}_{\beta}\text{>}=\underline{{\mu}_{\alpha}}\underline{{\omega}_{\alpha}^{2}}{\delta}_{\alpha \beta }\) éq 2.2-4

avec

\({\omega}_{\alpha}=\sqrt{({\lambda}_{\alpha})}\) éq 2.2-5

On note respectivement la matrice de masse généralisée, la matrice de raideur généralisée et la matrice de amortissement généralisée par:

\(\left[\underline{{M}_{n}}\right]={\left[\underline{{\Phi}_{n}}\right]}^{T}\left[\underline{M}\right]\left[\underline{{\Phi}_{n}}\right]\) éq 2.2-6

\(\left[{\underline{K}}_{n}\right]={\left[{\underline{\Phi}}_{n}\right]}^{T}\left[\underline{K}\right]\left[{\underline{\Phi}}_{n}\right]\) éq 2.2-7

\(\left[{\underline{D}}_{n}\right]={\left[{\underline{\Phi}}_{n}\right]}^{T}\left[\underline{D}\right]\left[{\underline{\Phi}}_{n}\right]\) éq 2.2-8

Résolution en transitoire

La projection \(\underline{{y}^{n}}(t)\) de \(\underline{y}(t)\) sur le sous espace engendré par les n premiers modes du système dynamique linéaire conservatif homogène associé s’écrit:

\(\underline{{y}^{n}}(t)=\left[\underline{{\Phi}_{n}}\right]\underline{{q}^{n}}(t)=\sum_{\alpha =1}^{n}{\underline{q}}_{\alpha}^{n}(t){\underline{\varphi}}_{\alpha}\) , éq 2.2.1-1

Les déplacements généralisés sont solutions du modèle matriciel réduit moyen (système dynamique non linaire),

\(\left[{M}_{n}\right]{\ddot{q}}^{n}(t)+\left[{D}_{n}\right]{\dot{q}}^{n}(t)+\left[{K}_{n}\right]{q}^{n}(t)+{F}_{c}^{n}(t,{q}^{n}(t),{\dot{q}}^{n}(t);W)={F}^{n}(t),\) éq 2.2.1-2

\(\dot{\underline{{q}^{n}}}(0)=\underline{{q}^{n}}(0)=0\) , éq 2.2.1-3

avec

\({F}^{n}(t)={\left[\underline{{\Phi}_{n}}\right]}^{T}f(t)\) , éq 2.2.1-4

\({F}_{c}^{n}(t,q,p;\underline{w})={[{\Phi}_{n}]}^{T}{f}_{c}(t,[{\Phi}_{n}]q,[{\Phi}_{n}]p,\underline{w})\) éq 2.2.1-5

Résolution en harmonique

La projection \({\underline{y}}^{n}(\omega )\) de \(\underline{y}(\omega )\) sur le sous espace engendré par les n premiers modes du système dynamique linéaire conservatif homogène associé s’écrit

\({q}^{n}(\omega )\) = \({[{\Phi}_{n}]}^{T}f(\omega ),\) éq 2.2.2-1

Les déplacements généralisés \({q}^{n}(\omega )\) sont solutions du modèle matriciel réduit moyen

\((-{\omega}^{2}\left[{\underline{M}}_{n}\right]+i\omega \left[{\underline{D}}_{n}\right]+\left[{\underline{K}}_{n}\right]){q}^{n}(\omega )={F}^{n}(\omega )\) éq 2.2.2-2

avec

\({F}^{n}(\omega )={\left[\underline{{\Phi}_{n}}\right]}^{T}f(\omega )\) , éq 2.2.2-3

Modèle probabiliste

Introduction du modèle probabiliste dans le problème dynamique

Afin de prendre en compte les incertitudes de modélisation et les incertitudes sur les données, une formulation probabiliste mixte non paramétrique – paramétrique est utilisée. Pour cela, le vecteur des \(n\) degrés de liberté généralisés \({\underline{q}}^{n}(t)\) (resp. \({\underline{q}}^{n}(\omega )\) ) est remplacé par une variable aléatoire \({Q}^{n}(t)\) (resp. \({Q}^{n}(\omega )\) ).

En transitoire, le processus stochastique \(t\to {Q}^{n}(t)\) indexé par \([0,T]\) et à valeur dans \({R}^{n}\) est solution du système dynamique non linéaire,

\(\left[{M}_{n}\right]{\ddot{Q}}^{n}(t)+\left[{D}_{n}\right]{\dot{Q}}^{n}(t)+\left[{K}_{n}\right]{Q}^{n}(t)+{F}_{c}^{n}(t,{Q}^{n}(t),{\dot{Q}}^{n}(t);W)={F}^{n}(t),\) éq 3.1-1

\(\dot{{Q}^{n}}(0)={Q}^{n}(0)=0\) , éq 3.1-2

et en harmonique, le processus stochastique \(t\mapsto {Q}^{n}(\omega )\) indexé sur \(\left[{\omega}_{1},{\omega}_{2}\right]\) et à valeur dans \({ℝ}^{n}\) est solution du système:

\((-{\omega}^{2}\left[{\underline{M}}_{n}\right]+i\omega \left[{\underline{D}}_{n}\right]+\left[{\underline{K}}_{n}\right]){Q}^{n}(\omega )={F}^{n}(\omega )\) éq 3.1-3

où, dans les deux cas transitoire et harmonique, \(\left[{\underline{M}}_{n}\right]\) , \(\left[{\underline{D}}_{n}\right]\) , et \(\left[{\underline{K}}_{n}\right]\) sont des matrices réelles symétriques définies positives pleines aléatoires et où \(W\) est une variable aléatoire à valeur dans \({ℝ}^{\nu}\) . L’introduction de matrices aléatoires dans les équations [éq 3.1-1] et [éq 3.1-3] permet de modéliser les incertitudes aléatoires associées à la partie linéaire du système dynamique. La variable aléatoire \(W\) à valeur vectorielle introduite dans l’équation [éq 3.1-1] permet de modéliser les incertitudes aléatoires concernant les paramètres des non linéarités de choc.

L’approche probabiliste paramétrique et l’approche probabiliste non paramétrique introduisent des matrices aléatoires (, et ) et une variable aléatoire \(W\) dont les lois de probabilité sont a priori non connues. Le choix d’un modèle probabiliste plutôt qu’un autre doit reposer uniquement sur l’information disponible (propriétés algébriques des matrices généralisées, valeurs moyennes des paramètres et des matrices généralisées, etc.). Afin de construire objectivement les lois de probabilité du modèle probabiliste des incertitudes, ([bib20] à [bib24]), le principe du maximum d’entropie ([bib14], [bib15], [bib17]) est utilisé avec un système de contraintes définies par cette information disponible. L’information disponible et le modèle probabiliste qui en découle sont présentés dans le prochain paragraphe.

Modèle probabiliste pour les matrices du système dynamique (incertitudes non paramétriques)

Information disponible sur les matrices du système dynamique

Le modèle probabiliste non paramétrique est construit en substituant les matrices \(\left[{\underline{M}}_{n}\right]\) , \(\left[{\underline{D}}_{n}\right]\) , et \(\left[{\underline{K}}_{n}\right]\) par des matrices aléatoires respectivement notées \(\left[{M}_{n}\right]\) , \(\left[{D}_{n}\right]\) , et \(\left[{K}_{n}\right]\) . Afin que le système dynamique probabiliste ainsi construit soit mécaniquement et statistiquement correct, la construction des matrices aléatoires \(\left[{M}_{n}\right]\) , \(\left[{K}_{n}\right]\) et \(\left[{D}_{n}\right]\) doit être telle que:

\(\left[{M}_{n}\right]\) , \(\left[{K}_{n}\right]\) et \(\left[{D}_{n}\right]\) sont des variables aléatoires du second ordre à valeurs dans l’ensemble des matrices symétriques réelles définies positives et de dimension \((n\times n)\) .

\(\left[{M}_{n}\right]\) , \(\left[{K}_{n}\right]\) et \(\left[{D}_{n}\right]\) \(\in {M}_{n}^{+}\) p.s. (presque sûrement), éq 3.2.1-1

où \({M}_{n}^{+}\) est l’ensemble des matrices réelles symétriques définies positives de dimension \((n\times n)\) .

Cette propriété algébrique est absolument nécessaire pour avoir un modèle d’équation aléatoire qui correspond à celle d’un système dynamique du second ordre amorti.

Les valeurs moyennes des matrices aléatoires \(\left[{M}_{n}\right]\) , \(\left[{K}_{n}\right]\) et \(\left[{D}_{n}\right]\) sont respectivement \(\left[{\underline{M}}_{n}\right]\) \(\left[{\underline{K}}_{n}\right]\) \(\left[{\underline{D}}_{n}\right]\) :

\(E\left\lbrace \left[{M}_{n}\right]\right\rbrace =\left[{\underline{M}}_{n}\right]\) , \(E\left\lbrace \left[{K}_{n}\right]\right\rbrace =\left[{\underline{K}}_{n}\right]\) et \(E\left\lbrace \left[{D}_{n}\right]\right\rbrace =\left[{\underline{D}}_{n}\right]\) , éq 3.2.1-2

où \(E\) désigne l’espérance mathématique.

Afin que la solution du système dynamique probabiliste soit aussi une variable du second ordre, on impose aux moments du second ordre des normes de Frobenius des matrices inverses \({\left[{M}_{n}\right]}^{\text{-1}}\) , \({\left[{K}_{n}\right]}^{\text{-1}}\) et \({\left[{D}_{n}\right]}^{\text{-1}}\) d’être finis:

\(E\lbrace {∥{\left[{M}_{n}\right]}^{\text{-1}}∥}_{F}^{2}\rbrace <+\infty\) \(E\lbrace {∥{\left[{K}_{n}\right]}^{\text{-1}}∥}_{F}^{2}\rbrace <+\infty\) \(E\lbrace {∥{\left[{D}_{n}\right]}^{\text{-1}}∥}_{F}^{2}\rbrace <+\infty\) éq 3.2.1-3

avec \({\parallel \left[A\right]\parallel }_{F}={(\text{tr}\left\lbrace \left[A\right]{\left[A\right]}^{T}\right\rbrace )}^{1/2}\) .

Remarque :

La seule propriété de positivité des matrices ne suffit pas et il faut s’assurer que leurs inverses sont du second ordre, d’où (3 (une variable aléatoire du second ordre presque sûrement inversible n’a pas dans le cas général une variable aléatoire inverse du second ordre). Pour plus de détails voir [bib19].

Construction du modèle probabiliste par le principe du maximum d’entropie

L’entropie «mesure» le niveau d’incertitude d’une loi de probabilité. Ainsi, si \({p}_{[A]}\) est la fonction de densité de probabilité correspondant à une matrice aléatoire \([A]\) (représentant les matrices \(\left[{M}_{n}\right]\) , \(\left[{K}_{n}\right]\) ou \(\left[{D}_{n}\right]\) ) de loi donnée, alors l’entropie (ou incertitude probabiliste) \(S({p}_{[A]})\) de \({p}_{[A]}\) est définie par:

\(S({p}_{[A]})=-{\int}_{{M}_{n}^{+}}{p}_{[A]}([A])\ln({p}_{[A]}([A]))\tilde{\mathrm{dA}}\) éq 3.2.2-1

Le principe du maximum d’entropie de Jayne consiste à construire la fonction de densité de probabilité qui maximise l’entropie probabiliste tout en vérifiant un système de contraintes. Dans le cas présent, le système de contraintes est défini par l’information disponible correspondant aux équations [éq 3.2.1-1] à [éq 3.2.1-3]. Pour la matrice aléatoire , ce système de contraintes s’écrit

\([A]\in {M}_{n}^{+}\) p.s. , \(E\lbrace [A]\rbrace =[\underline{A}]\) \(E\lbrace {∥{\left[A\right]}^{\text{-1}}∥}_{F}^{2}\rbrace <+\infty\) . éq 3.2.2-2

On montre alors que la matrice aléatoire \([A]\) est telle que (voir [bib20] à [bib24])

\([A]={[\underline{{L}_{A}}]}^{T}[{G}_{A}][\underline{{L}_{A}}]\) éq 3.2.2-3

où \([\underline{{L}_{A}}]\) est la matrice triangulaire inférieure issue de la factorisation de Cholesky de la matrice moyenne \([A]\) et où la fonction de densité de probabilité de la matrice aléatoire \([{G}_{A}]\) est définie sur l’ensemble \({M}_{n}^{+}\) par rapport à la mesure \(\tilde{d}A\) telle que:

\(\tilde{d}G={2}^{n(n-1)/4}{\Pi}_{1\le i\le j\le n}{\mathrm{dG}}_{ij}\) éq 3.2.2-4

\({p}_{[{G}_{A}]}([G])={1}_{{M}_{n}^{+}}([G])\times {C}_{{G}_{A}}\times {(\det[G])}^{(1-{\delta}_{A}^{2}){(2{\delta}_{A}^{2})}^{\text{-1}}(n+1)}\times {e}^{-(n+1){(2{\delta}_{A}^{2})}^{\text{-1}}\mathrm{tr}[G]}\) éq 3.2.2-5

avec

\({C}_{{G}_{A}}=\frac{{(2\pi )}^{-n(n-1)/4}{(\frac{n+1}{2{\delta}_{A}^{2}})}^{n(n+1){(2{\delta}_{A}^{2})}^{\text{-1}}}}{{\Pi}_{j=1}^{n}\Gamma (\frac{n+1}{2{\delta}_{A}^{2}}+\frac{1-j}{2})}\) éq 3.2.2-6

où \(\Gamma (z)={\int}_{0}^{+\infty }{t}^{z-1}{e}^{-t}\mathrm{dt}\) , éq 3.2.2-7

et où \({1}_{{M}_{n}^{+}}\) est la fonction indicatrice de \({M}_{n}^{+}\) , et où le paramètre \({\delta}_{A}\) contrôlant la dispersion de la matrice aléatoire \([A]\) est défini par:

\({\delta}_{A}={\left\lbrace \frac{E\left\lbrace {\parallel \left[{G}_{A}\right]-\left[{\underline{G}}_{A}\right]\parallel }_{F}^{2}\right\rbrace }{{\parallel \left[{\underline{G}}_{A}\right]\parallel }_{F}^{2}}\right\rbrace }^{1/2}\) , éq 3.2.2-8

La construction théorique du modèle fournit une borne admissible pour le niveau d’incertitude introduit. doit être choisi de sorte que

\(0<{\delta}_{A}<\sqrt{\frac{{n}_{0}+1}{{n}_{0}+5}}\) , éq 3.2.2-9

où \({n}_{0}\in \mathbb{N}\) est une constante du modèle probabiliste choisie de sorte que .

On démontre de plus que, sous les seules contraintes des équations [éq 3.2.1-1] à [éq 3.2.1-3], le principe du maximum d’entropie conduit à ce que les matrices aléatoires , ou soient statistiquement indépendantes dans leur ensemble.

Ce modèle probabiliste pour les matrices aléatoires réelles symétriques définies positives diffère des modèles plus classiques des matrices aléatoires basées sur les Ensembles Gaussiens et les ensembles Circulaires (références dans [bib19]). L’ensemble gaussien orthogonal utilisé par ailleurs pour des domaines hautes fréquences conduirait dans le domaine basses fréquences (dans lequel on se place) à des valeurs propres négatives, ce que l’on ne peut pas admettre pour les systèmes considérés. De plus, une matrice de l’ensemble gaussien orthogonal ne possède pas dans le cas général une matrice inverse du second ordre, ce qui conduirait à une solution du système dynamique de variance infinie, ce que l’on ne peut pas non plus admettre.

Modèle probabiliste pour les variables réelles (incertitudes paramétriques)

Information disponible sur les variables réelles

Dans l’approche probabiliste mixte, la modélisation probabiliste paramétrique consiste à substituer le paramètre \(\underline{w}\) des non-linéarités dans les systèmes dynamiques non linéaires données par [éq 2.1.1-1] ou [éq 2.1.2-3] par une variable aléatoire notée \(W=({W}_{1},\dots ,{W}_{\nu})\) . Dans une approche purement paramétrique (i.e. sans rendre aléatoires les matrices du système dynamique), la modélisation probabiliste paramétrique consiste à substituer certains paramètres \(\underline{w}\) des matrices \(\left[{\underline{M}}_{n}(\underline{w})\right]\) , \(\left[{\underline{K}}_{n}(\underline{w})\right]\) et \(\left[{\underline{D}}_{n}(\underline{w})\right]\) du système dynamique réduit moyen par une variable aléatoire \(W\) . Ces paramètres peuvent être par exemple des paramètres du matériau.

On suppose que les composantes de \(W\) sont des variables aléatoires réelles indépendantes entre elles et indépendantes des matrices aléatoires du système dynamique. Dans la suite, pour alléger l’écriture, on note \(W\) une coordonnée quelconque \({W}_{j}\) . La construction du modèle probabiliste nécessite de définir l’information disponible, laquelle constitue un système de contrainte sous lequel l’entropie de la densité de probabilité de la variable aléatoire \(W\) est maximisée.

L’information disponible est la suivante:

Le support de la variable aléatoire \(W\) est un intervalle \(D\) de \({ℝ}^{\nu}\)

\(W\in D,p.s..\) éq 3.3.1-1

La valeur moyenne de la variable aléatoire \(W\) est \(\underline{w}\) :

\(E\left\lbrace W\right\rbrace =\underline{w}\) . éq 3.3.1-2

Eventuellement, selon l’information effectivement disponible, le moment du second ordre de la variable aléatoire \({\mid W\mid }^{-1}\) est fini:

\(E\left\lbrace {\mid W\mid }^{-2}\right\rbrace \text{<+}\infty .\) éq 3.3.1-3

Construction du modèle probabiliste par le principe du maximum d’entropie

Si \({p}_{W}\) est la fonction de densité de probabilité correspondant à la variable aléatoire \(W\) alors l’entropie (incertitude probabiliste) \(S({p}_{W})\) de \({p}_{W}\) est définie par:

\(S({p}_{W})=-{\int}_{-\infty }^{+\infty }{p}_{W}(w)\ln({p}_{W}(w))\text{dw}\) , éq 3.3.2-1

En utilisant le principe du maximum d’entropie, on obtient trois densités de probabilité suivant la nature du support \(D\) et selon que la contrainte correspondant à l’équation [éq 3.2.2-6] est considérée ou non.

Support fermé borné sans information sur l’inverse

S’il existe deux réels \(a\) et \(b\) tels que \(D=\left[a,b\right]\) et si l’information disponible est donnée par les équations [éq 3.3.1-1] et [éq 3.2.2-5], alors la variable aléatoire \(W\) suit une loi exponentielle tronquée dont la fonction de densité de probabilité est:

\({p}_{W}(w)={1}_{[a,+\infty ]}(w)\frac{k}{\alpha (k)}\exp(-\text{kw})\) éq 3.3.3-1

où \({1}_{[a,b]}\) est la fonction indicatrice de et où \(\alpha (k)\) et \(k\) sont tels que :

\((\underline{{w}_{k}}-1)\alpha (k)-k\beta (k)=0\) éq 3.3.3-2

avec

\(\alpha (k)={e}^{-\mathrm{ak}}-{e}^{-\mathrm{bk}}\) , éq 3.3.3-3

et

\(\beta (k)=a{e}^{-\mathrm{ak}}-b{e}^{-\mathrm{bk}}\) éq 3.3.3-4

Support semi fermé non borné sans information sur l’inverse

S’il existe un réel \(a\) tel que \(D=[a,+\infty [:ref:`\) et si l’information disponible est donnée par les équations [éq 3.3.1-1 <` et si l’information disponible est donnée par les équations [éq 3.3.1-1>`] et [éq 3.3.1-2] alors la variable aléatoire \(W\) suit une loi exponentielle dont la fonction de densité de probabilité est:

\({p}_{W}(w)={1}_{[a,+\infty ]}(w)\frac{1}{\underline{w}-a}\exp(-\frac{w-a}{\underline{w}-a})\) , éq 3.3.4-1

où \({1}_{[a,+\infty ]}\) est la fonction indicatrice de \([a,+\infty [\) .

Support semi fermé non borné avec information sur l’inverse

S’il existe un réel \(a\) tel que \(D=[a,+\infty [:ref:`\) et si l’information disponible est donnée par les équations. [éq 3.3.1-1 <` et si l’information disponible est donnée par les équations. [éq 3.3.1-1>`], [éq 3.3.1-2] et [éq 3.3.1-3], alors la variable aléatoire \(W\) suit une loi gamma dont la fonction de densité de probabilité est,

\({p}_{W}(w)={1}_{[a,+\infty ]}(w)\frac{(\underline{w}{\delta}^{2}-a{\delta}^{2}{)}^{-1/{\delta}^{2}}}{\Gamma (1/{\delta}^{2})}(w-a{)}^{(\text{1-}{\delta}^{2})/{\delta}^{2}}\exp\left\lbrace -\frac{w-a}{(\underline{w}-a){\delta}^{2}}\right\rbrace\) , éq 3.3.5-1

où \(\delta ` est un paramètre contrôlant le niveau d’incertitude de la variable aléatoire qui s’écrit :math:`W\) (de façon analogue au cas non paramétrique [éq 3.2.2-8]) :

\(\delta ={\left\lbrace \frac{E\left\lbrace {(w-\underline{w})}^{2}\right\rbrace }{{\underline{w}}^{2}}\right\rbrace }^{1/2}\) , éq 3.3.5-2

Construction de la réponse stochastique et des statistiques associées

Cas transitoire

Réponse transitoire stochastique

Les excitations du système dynamique sont supposées déterministes, mais au paragraphe [§3.1], des matrices et des paramètres aléatoires ont été introduits dans le modèle matriciel réduit. Donc, la réponse transitoire \(t\to {Q}_{n}(t)\) est un processus stochastique non stationnaire indexé par \([0,T]\) à valeur dans \({ℝ}^{n}\) (en utilisant quelques hypothèses supplémentaires d’existence, d’unicité et de régularité de la solution déterministe, cf. [bib19]).

Par voie de conséquence, au vecteur des \(m\) degrés de liberté libres \(\underline{{y}^{n}}(t)\) correspond le processus stochastique \({Y}^{n}(t)\) indexé par \([0,T]\) et à valeur dans \({ℝ}^{m}\) tel que

\({Y}^{n}(t)=\left[{\underline{\Phi}}_{n}\right]{Q}^{n}(t),\) éq 3.4.1.1-1

Dans le cas du passage en coordonnées relatives, au processus stochastique \(t\to {Y}_{n}(t)\) indexé par \([0,T]\) et à valeur dans \({ℝ}^{m}\) défini par l’équation [éq 3.4.1.1-1] correspond le processus stochastique \(t\to {Z}_{n}(t)\) des d.d.l. libres de la structure en coordonnées absolues indexée par et à valeur dans \({ℝ}^{m}\) telle que

\({Z}_{n}(t)={Y}_{n}(t)+[R]{z}_{s}(t)\) . éq 3.4.1.1-2

Spectre de réponse élastique

On note \({Z}_{j}^{n}(t)\) la \(j\) ème composante du vecteur \({Z}_{n}(t)\) correspondant à une réalisation aléatoire de la réponse stochastique du \(j\) ème degré de liberté libre de la structure. \({Z}_{j}^{n}(t)\) peut-être caractérisée par son spectre de réponse élastique (appelé également spectre d’oscillateur dans la documentation de Code_Aster ) que l’on note \({S}_{j}(\xi ,\omega )\) où \(\xi\) et \(\omega\) sont respectivement le taux d’amortissement et la pulsation associés.

Avec des hypothèses raisonnables, en particulier sur la régularité de l’application non linéaire \((t,q,p;w)\mapsto {F}_{c}^{n}(t,q,p;w)\) , on peut démontrer que est un processus du second ordre dont les trajectoires sont presque sûrement continues. Par conséquent, pour tout \(\xi\) fixé dans un intervalle \({J}_{\xi}\) donné, \(\omega \to {S}_{j}(\xi ,\omega )\) est un processus stochastique indexé sur la bande d’analyse \({J}_{\omega}\) à valeur dans \({ℝ}^{+}\) . On admet que ce processus est du second ordre, c’est-à-dire:

\(E\lbrace {S}_{j}{(\xi ,\omega )}^{2}\rbrace <+\infty ,\forall \omega \in {J}_{\omega}\) . éq 3.4.1.2-1

Cas harmonique

Au paragraphe [§3.1], des matrices et des paramètres aléatoires ont été introduits dans le modèle matriciel réduit. La réponse harmonique \(t\mapsto {Q}^{n}(\omega )\) est donc un processus stochastique indexé sur \(\left[{\omega}_{1},{\omega}_{2}\right]\) à valeur dans \({R}^{n}\) .

Par voie de conséquence, au vecteur des \(m\) degrés de liberté libres \({\underline{y}}^{n}(\omega )\) correspond le processus stochastique \({Y}^{n}(\omega )\) indexé sur \(\left[{\omega}_{1},{\omega}_{2}\right]\) et à valeur dans \({R}^{m}\) tel que

\({Y}^{n}(\omega )=\left[{\underline{\Phi}}_{n}\right]{Q}^{n}(\omega )\) éq 3.4.2-1

\({Y}_{j}^{n}(\omega )\) la \(j\) ème composante du vecteur \({Y}^{n}(\omega )\) est une variable aléatoire que l’on admettra du second ordre.

Construction de la réponse stochastique par la méthode de Monte Carlo

Choix et mise en œuvre de la méthode de Monte Carlo

Les réponses et les spectres de réponse correspondent à des transformations fortement non linéaires des matrices aléatoires et des paramètres aléatoires qui résultent de la modélisation probabiliste des incertitudes. De plus, on ne peut bien sûr construire que des approximations numériques des ces réponses et de ces spectres de réponse. Les statistiques (premiers moments statistiques, probabilité de dépassement d’un seuil, …) s’écrivent formellement comme des intégrales multiples de très grande dimension car le nombre de variables aléatoires du modèle probabiliste est par construction élevé. Enfin, le nombre de grandeurs observées est très grand (plusieurs ddl pour plusieurs fréquences). Pour toutes ces raisons, la méthode la plus adaptée pour calculer la solution probabiliste (réponse stochastique et statistiques associées) est la méthode de simulation numérique de Monte-Carlo.

La méthode de simulation de Monte-Carlo présente l’avantage de donner des résultats dont on peut contrôler la précision (vérification de la convergence, cf. [§4.1]), contrairement à la plupart des méthodes basées sur des approximations. Elle peut être coûteuse en temps de calculs, mais l’utilisation des techniques de réduction de la variance peut permettre de réduire le nombre de simulations nécessaires (cf. [bib8] ou [bib9]).

La mise en œuvre de la méthode de Monte Carlo consiste pour le problème qui nous concerne à générer \({n}_{s}\) réalisations des matrices aléatoires \(\left[{M}_{n}\right]\) , \(\left[{K}_{n}\right]\) et \(\left[{D}_{n}\right]\) du système dynamique et/ou \({n}_{s}\) réalisations de la variable aléatoire vectorielle \(W\) . Les résolutions du système dynamique déterministe pour chacune des \({n}_{s}\) réalisations de ( \(\left[{M}_{n}\right]\) , \(\left[{K}_{n}\right]\) , \(\left[{D}_{n}\right]\) , \(W\) ) produisent \({n}_{s}\) réalisations du processus stochastique solution \(t\mapsto {Q}^{n}(t)\) (resp. \(t\mapsto {Q}^{n}(\omega )\) ) et par suite de \(t\to {Y}_{n}(t)\) , de \(t\to {Z}_{n}(t)\) et de \(\omega \to {S}_{j}(\xi ,\omega )\) (resp. \(\omega \mapsto {Y}_{j}^{n}(\omega )\) ). La génération des matrices aléatoires est traitée dans le paragraphe suivant ; la génération de la variable aléatoire \(W\) est plus classique et n’est pas rappelée.

Génération des matrices pseudo-aléatoires

Afin de générer les réalisations de la matrice aléatoire \(\left[{G}_{A}\right]\) , on utilise la représentation algébrique suivante de la matrice aléatoire \(\left[{G}_{A}\right]\) dont la loi de probabilité est définie par les équations [éq3.2.1‑2], [éq 3.2.2-1]:

\(\left[{G}_{A}\right]={\left[L\right]}^{T}\left[L\right]\) , éq 3.4.3.2-1

la matrice aléatoire triangulaire \(\left[L\right]\) étant telle que:

• Les variables aléatoires \({\left[L\right]}_{ij},i\le j\) sont indépendantes.

• Pour \(i<j\) , les variables aléatoires réelles \({\left[L\right]}_{ij}\) s’écrivent \({\left[L\right]}_{ij}={\sigma}_{n}{U}_{ij}\) où \({\sigma}_{n}={\delta}_{A}{(n+1)}^{-1/2}\) et où \({U}_{ij}\) est une variable aléatoire réelle gaussienne de moyenne 0 et de variance 1.

• Pour \(i=j\) , les variables aléatoires réelles \({\left[L\right]}_{ij}\) s’écrivent \({\left[L\right]}_{ij}={\sigma}_{n}{(2{V}_{j})}^{1/2}\) où \({\sigma}_{n}\) est défini précédemment et où \({V}_{j}\) est variable réelle positive aléatoire de loi gamma dont la fonction de densité de probabilité \({p}_{{v}_{j}}(v)\) par rapport à la mesure \(\mathit{dv}\) s’écrit:

\({p}_{{v}_{j}}(v)=\frac{{1}_{[0,+\infty [}(v)}{\Gamma ((n+1)/(2{\delta}_{A}^{2})+(1-j)/2)}{v}^{(n+1)/(2{\delta}_{A}^{2})+(1-j)/2}{e}^{-v}\) , éq 3.4.3.2-2

où \({1}_{[0,+\infty [}\) est la fonction indicatrice de \([0,+\infty [\) .

Statistiques sur les spectres

Dans ce chapitre, on présente la définition des statistiques des spectres de réponse élastique \({S}_{j}(\xi ,\omega )\) , dans le cas d’une résolution transitoire. Dans le cas harmonique, les statistiques sur les variables aléatoires \({Y}_{j}^{n}(\omega )\) se définissent de la même façon et ne sont donc pas présentées.

Estimation des quantiles

Pour tout \((\xi ,\omega )\in {J}_{\xi}\times {J}_{\omega}\) , \({S}_{j}(\xi ,\omega )\) est une variable aléatoire à valeur dans \({ℝ}^{+}\) . On cherche à estimer le quantile associé à la probabilité \(\alpha\) noté \({S}_{j,\alpha }(\xi ,\omega )\) et défini par:

\({S}_{j,\alpha }(\xi ,\omega )={F}_{\xi ,\omega }^{-1}(1-\alpha )\) éq 3.4.4.1-1

où \({F}_{\xi ,\omega }\) est la fonction de répartition inconnue de \({S}_{j}(\xi ,\omega )\) .

Soit \(({S}_{j}(\xi ,\omega ;{\theta}_{1})\text{, ... ,}{S}_{j}(\xi ,\omega ;{\theta}_{r})\text{, ... ,}{S}_{j}(\xi ,\omega ;{\theta}_{{n}_{s}}))\) l’échantillon composé des \({n}_{s}\)

réalisations de \({S}_{j}(\xi ,\omega )\) et \(({S}_{j}(\xi ,\omega ;{\theta}_{(1)})\text{, ... ,}{S}_{j}(\xi ,\omega ;{\theta}_{(r)})\text{, ... ,}{S}_{j}(\xi ,\omega ;{\theta}_{({n}_{s})}))\) l’échantillon ordonné associé.

Un estimateur naturel du quantile \({S}_{j,\alpha }(\xi ,\omega )\) pour \(\alpha =\frac{r}{{n}_{S}}\) , \(1\le r\le {n}_{S}\) est:

\({\hat{S}}_{j,\frac{r}{{n}_{S}}}(\xi ,\omega )={S}_{j}(\xi ,\omega ;{\theta}_{(r)})\) éq 3.4.4.1-2

Pour obtenir un estimateur plus robuste du quantile, on peut « moyenner » l’estimateur sur plusieurs séries de \({n}_{s}\) réalisations. Si la probabilité souhaitée est telle que \(\alpha <1/{n}_{S}\) , ou si l’on souhaite réduire le nombre de simulations, il est possible d’utiliser des estimateurs plus sophistiqués, par exemple en supposant que la fonction de répartition \({F}_{\xi ,\omega }\) appartient à un domaine d’attraction donné (théorie des valeurs extrêmes) ou par exemple en utilisant une méthode de régularisation bayésienne (cf.[bib12]).

Valeurs extrêmes d’échantillon

Pour un échantillon de \({n}_{s}\) réalisations de \({S}_{j}(\xi ,\omega )\) notées \({S}_{j}(\xi ,\omega ;{\theta}_{1})\text{, ... ,}{S}_{j}(\xi ,\omega ;{\theta}_{{n}_{s}})\) on définit les valeurs extrêmes d’échantillon par:

\(\omega \mapsto {\text{dB}}_{j\text{,min}}(\xi ,\omega ;{n}_{S})={\log}_{10}(\underset{r=1,...,{n}_{S}}{\min}{S}_{j}(\xi ,\omega ;{\theta}_{r}))\) éq 3.4.4.2-1

\(\omega \mapsto {\text{dB}}_{j\text{,max}}(\xi ,\omega ;{n}_{S})={\log}_{10}(\underset{r=1,...,{n}_{S}}{\max}{S}_{j}(\xi ,\omega ;{\theta}_{r}))\) éq 3.4.4.2-2

« Domaine de confiance » établi à partir de l’inégalité de Tchebychev

Pour un échantillon de \({n}_{s}\) réalisations du processus \(\omega \mapsto {S}_{j}(\xi ,\omega )\) notées \(\omega \mapsto {S}_{j}(\xi ,\omega ;{\theta}_{1})\) , …, \(\omega \mapsto {S}_{j}(\xi ,\omega ;{\theta}_{{n}_{S}})\) , on peut construire le « domaine de confiance » de la variable aléatoire pour tout , en utilisant l’inégalité de Tchebychev associé à un niveau de probabilité \({P}_{C}\) :

\(\text{Proba}\left\lbrace {\text{dB}}_{j}^{-}(\xi ,\omega )<{\text{dB}}_{j}(\xi ,\omega )<{\text{dB}}_{j}^{+}(\xi ,\omega )\right\rbrace \ge {P}_{C}\) , éq 3.4.4.3-1

où l’enveloppe inférieure et l’enveloppe supérieure sont définies par:

\({\text{dB}}_{j}^{+}(\xi ,\omega )={\log}_{10}({m}_{\mathrm{1j}}(\xi ,\omega )+\frac{{\sigma}_{j}(\xi ,\omega )}{\sqrt{1-{P}_{C}}})\) , éq 3.4.4.3-2

\({\text{dB}}_{j}^{-}(\xi ,\omega )=2{\log}_{10}({m}_{\mathrm{1j}}(\xi ,\omega ))-{\text{dB}}_{j}^{+}(\xi ,\omega )\) . éq 3.4.4.3-3

avec la moyenne et l’écart type de \({\text{dB}}_{j}(\xi ,\omega )\) :

\({m}_{\mathrm{1j}}(\xi ,\omega )=E\lbrace {S}_{j}(\xi ,\omega )\rbrace\) éq 3.4.4.3-4

\({\sigma}_{j}(\xi ,\omega )=E{\left\lbrace {({S}_{j}(\xi ,\omega )-{m}_{\mathrm{1j}}(\xi ,\omega ))}^{2}\right\rbrace }^{1/2}\) . éq 3.4.4.3-5

Le « domaine de confiance » ainsi construit s’est avéré être une bonne approximation des valeurs extrêmes d’échantillon pour le cas traité dans [bib21]. Or, ce « domaine de confiance » ne fait intervenir que les deux premiers moments dont les estimateurs convergent plus rapidement vis-à-vis du nombre \({n}_{s}\) de simulations que les valeurs extrêmes d’échantillon. Il peut donc être intéressant d’utiliser cette construction du « domaine de confiance » plutôt qu’une construction basée sur l’estimation des quantiles plus coûteuse en nombre de simulations.

Remarque :

Le terme « domaine de confiance », peut être considéré par certains comme un abus de langage. On devrait plutôt utiliser la terminologie moins intuitive « domaine inter-quantiles ». En effet, dans la littérature statistique, un intervalle de confiance est théoriquement l’intervalle dans lequel se trouve la vraie valeur d’un paramètre d’une variable aléatoire (par exemple sa moyenne) avec une probabilité donnée. Cette terminologie est employée dans le cadre très précis de la théorie de l’estimation ensembliste. L’intervalle de confiance n’est pas une caractérisation de la variabilité d’une variable aléatoire, contrairement à un écart-type ou à des quantiles. On utilise quand même »domaine de confiance » avec parcimonie dans la suite, car c’est certainement un peu plus parlant pour les non-spécialistes des statistiques.

Mise en œuvre dans Code_Aster

Etude de la convergence stochastique du modèle numérique

Cas transitoire

La convergence de la solution stochastique doit être étudiée par rapport au nombre \(n\) de modes et au nombre \({n}_{s}\) de simulations de Monte Carlo. Comme la solution stochastique est un processus du second ordre (par hypothèse, cf. [§3.4.1.2]), sa convergence peut être analysée en étudiant les applications tel que :

\({\mid \parallel {\ddot{Z}}_{j}^{n}\parallel \mid }^{2}={\int}_{0}^{T}E{\left\lbrace {\ddot{Z}}_{j}^{n}(t)\right\rbrace }^{2}\text{dt}\) , éq 4.1.1-1

où est un processus stochastique du second ordre indexé par et à valeur dans \(R\) représentant l’accélération du \(j\) ème degré de liberté de la structure.

Dans le cadre des simulations de Monte-Carlo, cette norme est estimée pour \(n\) fixé à partir d’un ensemble de \({n}_{s}\) réalisations aléatoires par l’approximation avec

\({\text{conv}}_{j}(n,{n}_{S}{)}^{2}={\int}_{0}^{T}{(\frac{1}{{n}_{S}}\sum_{i=1}^{{n}_{S}}{\ddot{Z}}_{j}^{n}(t;{\theta}_{i}))}^{2}\text{dt}\) éq 4.1.1-2

La convergence stochastique du modèle est ainsi analysée suivant la dimension du modèle réduit (c’est-à-dire le nombre de mode \(n\) du sous espace propre du modèle éléments finis moyen sur lequel le système dynamique non linéaire stochastique a été projeté au paragraphe [§2.2]) et le nombre \({n}_{s}\) de simulations de Monte-Carlo en étudiant la fonction .

Cas harmonique

La convergence dans le cas d’une résolution transitoire peut se transposer directement dans le cas d’une résolution harmonique, avec la norme :

\({\mid \parallel {Z}_{j}^{n}\parallel \mid }^{2}={\int}_{{\omega}_{1}}^{{\omega}_{2}}E{\left\lbrace {Z}_{j}^{n}(\omega )\right\rbrace }^{2}d\omega\) , éq 4.1.2-1

Choix des paramètres de dispersion

Pour utiliser la méthode, les paramètres de dispersion \(\delta\) doivent être fixés. Deux approches peuvent être a priori utilisées pour fixer la valeur de ces paramètres.

La première approche consiste à identifier la valeur des paramètres \(\delta\) pour une structure donnée ou pour une classe de structure à l’aide de méthodes appropriées. Pour cela, on peut utiliser des résultats expérimentaux des réponses dynamiques de la structure. On peut aussi utiliser des simulations numériques construites en utilisant une approche paramétrique des incertitudes. Dans ce dernier cas, il est à noter que seules les erreurs sur les données du modèle sont prises en compte, puisque les erreurs de modélisation ne peuvent pas être prise en compte par l’approche paramétrique.

La deuxième approche consiste à ne pas fixer a priori une valeur fixe des paramètres :math:`delta ` mais à les faire varier dans une plage donnée (seulement 3 scalaires à faire varier pour les matrices de masse, de raideur et d’amortissement sur la partie non-paramétrique en comparaison avec le très grand nombre de paramètres à faire varier simultanément dans une étude paramétrique classique). Cette approche permet de mener une analyse globale de sensibilité aux incertitudes. Dans le cas d’absence d’information objective sur les paramètres dispersion à choisir, il est préférable d’utiliser une telle approche. La méthode non paramétrique proposée apparaît alors comme une approche robuste et simple d’analyse de sensibilité aux incertitudes.

Principales étapes

La mise en œuvre dans Code_Aster est composée de trois principales étapes: la construction du modèle matriciel réduit moyen, la génération des réalisations de la réponse vue comme un processus stochastique, et enfin le post-traitement statistique de ces réalisations. Les deux dernières étapes constituent en fait la méthode de simulation numérique de Monte Carlo directe.

Étape 1: construction du modèle matriciel réduit moyen

Le modèle matriciel réduit moyen est construit à l’aide d’un enchaînement classique d’opérateurs dépendant de l’analyse précise effectuée dont les principaux peuvent être: ASSE_MATRICE, CALC_MODES, MODE_STATIQUE, CALC_CHAR_SEISME, PROJ_BASE, …

Étape 2: génération des réalisations de la réponse transitoire

Les \({n}_{s}\) réalisations de la réponse transitoire stochastique sont calculées dans une boucle en langage Python composée de:

• Génération des \(p\) ième réalisations des matrices généralisées aléatoires de masse, de raideur et d’amortissement par genE_matR_alea (doc. [U4.36.06]). Ces matrices ne sont pas diagonales et nécessitent donc un stockage plein.

• Génération des \(p\) ième réalisations des variables aléatoires des paramètres des non-linéarités par genE_VARI_alea (doc. [U4.36.07]).

• Calcul de la \(p\) ième réalisation \({Q}^{n}(t;p)\) ou \({Q}^{n}(\omega ;p)\) solution du système matriciel stochastique \(s\) . Cette réalisation est la solution du système matriciel classique dont les matrices et les seconds membres sont les réalisations précédemment générées. Le calcul est donc effectué par dyna_tran_modal ou DYNA_LINE_HARM (avec des matr_asse_GENE_R et vect_asse_GENE en entrée).

• 1-Extraction des observations temporelles des degrés de liberté physiques prédéfinis (par exemple \({\ddot{Z}}_{{}_{i}}^{n}(t;p)\) ou \({Y}_{j}^{n}(\omega ;p)\) , mais aussi éventuellement les champs de déplacement, vitesse, contraintes, etc) via recu_fonction (après un REST_GENE_PHYS pour \({Y}_{j}^{n}(\omega ;p)\) ).

2-Calcul des spectres correspondant (par calc_fonction(SPEC_OSCI) pour \(\omega \to {S}_{j}(\xi ,\omega ;p)\) et calc_fonction(MODULE) pour \({Y}_{j}^{n}(\omega ;p)\) ) .

• Evaluation, via calc_fonction mots-clé COMB ou PUISSANCE ou ENVELOPPE, des contributions aux estimateurs des moyennes, des moments d’ordre deux, des valeurs extrêmes max. et min. d’échantillon pour les spectres normalisés:

\({\stackrel{ˆ}{m}}_{\mathrm{1j}}(x,w;p)={S}_{j}(x,w;p)+{\stackrel{ˆ}{m}}_{\mathrm{1j}}(x,w;p-1)\) , \({\stackrel{ˆ}{m}}_{\mathrm{2j}}(x,w;p)={S}_{j}(x,w;p{)}^{2}+{\stackrel{ˆ}{m}}_{\mathrm{2j}}(x,w;p-1)\) ,

\({\stackrel{ˆ}{S}}_{j\text{,max}}(x,w;p)=\text{Max}\left\lbrace {S}_{j}(x,w;p),{\stackrel{ˆ}{S}}_{j\text{,max}}(x,w;p-1)\right\rbrace\) , \({\stackrel{ˆ}{S}}_{j\text{,min}}(x,w;p)=\text{Min}\left\lbrace {S}_{j}(x,w;p),{\stackrel{ˆ}{S}}_{j\text{,min}}(x,w;p-1)\right\rbrace\) .

Etape 3: post-traitements statistiques

Les moyennes, les écart-types, les valeurs extrêmes max. et min. d’échantillon pour les spectres normaliséspeuvent être évalués via calc_fonction(COMB):

\({m}_{\mathrm{1j}}(x,w)=\frac{1}{{n}_{s}}{\stackrel{ˆ}{m}}_{\mathrm{1j}}(x,w;{n}_{s})\) , \({m}_{\mathrm{2j}}(x,w)=\frac{1}{{n}_{s}}{\stackrel{ˆ}{m}}_{\mathrm{2j}}(x,w;{n}_{s})\) .

Les intervalles de confiance peuvent alors être tracés à partir des valeurs extrêmes d’échantillon ou des bornes obtenues par Tchebychev cf. [§3.4.4].

Dans le cas transitoire, un exemple est donné par un cas test d’une plaque en flexion avec non linéarités de choc, cf. doc. [V5.06.001] [bib1]. D’autres détails sont donnés dans la doc. [U2.08.05] [bib2].

Efficacité numérique de l’approche non paramétrique

L’approche non paramétrique est plus économique en temps de calcul qu’une approche purement paramétrique dans laquelle les paramètres de géométrie, matériaux, etc. sont des variables aléatoires.

Dans l’approche purement paramétrique, le modèle éléments finis dépend des paramètres incertains. Pour chaque simulation de Monte-Carlo, le modèle éléments finis est différent. Il faut donc, pour chaque simulation, calculer les matrices élémentaires, effectuer les assemblages, passer en coordonnées relatives, résoudre le problème aux valeurs propres, projeter sur base modale, résoudre le système réduit et revenir en base physique puis en coordonnées relatives.

Dans l’approche non paramétrique, seul le système réduit est différent à chaque simulation. Il est donc simplement nécessaire, à chaque simulation, de résoudre le système réduit et revenir en base physique puis en coordonnées relatives. En particulier, la résolution du problème aux valeurs propres du modèle éléments finis moyens est effectuée une fois pour toutes, avant les simulations de Monte‑Carlo.

Le gain en temps de calcul qui en résulte est variable, mais il peut être important. En première approximation, ce gain en temps calcul dépend du ratio entre le temps CPU nécessaire à la résolution aux valeurs propres et le temps CPU nécessaire à la résolution du système réduit. Plus ce ratio est grand, plus l’approche non paramétrique est avantageuse par rapport à l’approche purement paramétrique. En particulier, le gain en temps de calcul peut être très important pour des structures avec un très grand nombre de degrés de liberté et une base modale de petite dimension.

Vérification

Les fonctionnalités et les méthodes présentées dans ce document sont vérifiées par les cas tests suivants:

SDNS01	Modèle probabiliste non paramétrique – paramétrique d’une plaque en flexion avec non linéarités localisées de choc	[V5.06.001]
SHLS200	Modèle probabiliste non paramétrique : Réponse harmonique d’une plaque sous-structurée	[V2.06.200]

Bibliographie

1. CAMBIER S., DESCELIERS C.: Modèle probabiliste non paramétrique – paramétrique d’une plaque en flexion avec non linéarités localisées de choc, documentation Code_Aster [V5.06.001].

2. CAMBIER S., DESCELIERS C.: Simulation numérique de Monte Carlo, documentation Code_Aster [U2.08.05].

3. CAMBIER S., DESCELIERS C.: Opérateur GENE_VARI_ALEA, documentation Code_Aster [U4.36.06].

4. CAMBIER S., DESCELIERS C.: Opérateur GENE_MATR_ALEA, documentation Code_Aster [U4.36.07].

5. CAMBIER S., DESCELIERS C., SOIZE S.: Prise en compte probabiliste des incertitudes dans l’estimation du comportement sismique d’un circuit primaire, Note EDF-R&D HT‑62/03/005/A, 2003

6. CAMBIER S., DESCELIERS C., SOIZE S.: Prise en compte probabiliste des incertitudes dans l’estimation du comportement sismique d’un circuit primaire. 7ème Colloque National AFPS, 2003.

7. CAMBIER S.: Application de l’approche probabiliste à l’assemblage combustible dans le cadre des « 3% AMA » en 2002, CR-AMA-03-004, Décembre 2002

8. CAMBIER S.: Approche probabiliste pour la prise en compte de la dispersion de paramètres mécaniques - Application à la fatigue vibratoire, Note EDF-R&D HT-64/00/040/A, Décembre 2000

9. CAMBIER S., GUIHOT P., COFFIGNAL G. : Computational methods for accounting of structural uncertainties, applications to dynamic behaviour prediction of piping systems. Structural Safety, 24: 29-50, 2002.

10. CHEBLI H.: Modélisation des Incertitudes Aléatoires non homogènes en Dynamique des Structures pour le Domaine des Basses Fréquences Thèse CNAM-ONERA, 2002

11. DESCELIERS C., SOIZE C., CAMBIER S. : Nonparametric–parametric model for random uncertainties in nonlinear structural dynamics – Application to earthquake engineering, soumis à Earthquaque Engineering and Structural Dynamics en 2003

12. DURBEC V., TROTTIER C (INRIA), DIEBOLT J. (CNRS): Régularisation de distributions pour une meilleure adéquation extrême: rapport de synthèse du CRECO EDF/INRIA EP 902, Note EDF-R&D HP-26/00/008/A

13. IBRAHIM R.A. : Structural dynamics with parameter uncertainties, Applied Mechanics Reviews, Vol. 40, no. 3, pp.309-328, 1987

14. JAYNES E. T., Information theory and statistical mechanics, Physical Review, 106(4), 620‑630 and 108(2), 171-190, 1957

15. KAPUR J.N, KEYSAVAN H. K. : Entropy Optimization Principles with Applications, Academic Press, San Diego, 1992

16. KATAFYGIOTIS L.S., PAPADIMITRIOU C. : Dynamic Response Variability of Structures with Uncertain Properties, Earthquake Engineering and Structural Dynamics, Vol. 25, pp.775-793, 1996

17. SHANNON C. E. : A mathematical theory of communication, Bell System Tech. J., 27, 379‑423 and 623-659, 1948

18. SOIZE C.: Méthodes d’études des problèmes classiques de dynamiques stochastiques, Techniques de l’Ingénieur, traité Sciences fondamentales, A 1 346, 61p, 1988

19. SOIZE C.: Modèle Probabiliste Mixte Paramétrique – Non Paramétrique des Incertitudes en Dynamique Non Linéaire des Structures, Rapport scientifique CRECO T62/E28858, 2001

20. SOIZE C.: Modélisation des incertitudes aléatoires en élastodynamique transitoire, Comptes Rendus de l’Académie des Sciences - Séries IIB – Mechanics, 329(3), 225-230, 2001.

21. SOIZE C. : Nonlinear dynamical System with nonparametric model of random uncertainties, The Uncertainties in Engineering Mechanics Journal, e-journal form Resonance Publication, http://www.resonance-pub.com/eprints.htm, 2001

22. SOIZE C. : Transient responses of dynamical systems with random uncertainties, Probabilistic Engineering Mechanics, 16(4), 363-372, 2001

23. SOIZE C. : Maximum entropy approach for modeling random uncertainties in transient elastodynamics, J. Acoust. Soc. Am., 109(5), 1979-1996, 2001.

24. SOIZE C. : Nonlinear dynamical System with nonparametric model of random uncertainties, The Uncertainties in Engineering Mechanics Journal, e-journal form Resonance Publication, http://www.resonance-pub.com/eprints.htm, 2001

Description des versions du document

Version Aster	Auteur(s) Organisme(s)	Description des modifications
6.4	S. CAMBIER, C. DESCELIERS EDF-R&D/AMA	Texte initial