3Approche quasistatique push-over. Code_Aster réservoir métallique

Code_Aster

Titre :

Analyse de la tenue sismique des grands réservoirs[...]

Responsable :

Nicolas GREFFET

Version default

Date :

21/04/2009

Page :

7/20

Clé :

U2.06.11

Révision :

931

3 Approche quasistatique push-over

3.1

Chargements imposés

Le chargement imposé est d’origine réglementaire (EC-8), quasistatique [bib1]. On représente les effets d’un séisme sur la structure par un champ de pression variable imposé sur la face interne des viroles. La valeur en chaque point est fonction de la coordonnée courante et croît linéairement avec le temps. Cette évolution monotone en temps est caractéristique des méthodes dites push-over, au sens

EC-8, dont l’objectif est de simuler par un calcul quasistatique la réponse à une sollicitation sismique, et donc de nature physique dynamique transitoire. Les effets dynamiques, comme l’inertie et les efforts générés par le fluide en ballottement sont remplacés par cette distribution de pression imposée. Le problème à traiter ne fait intervenir qu’une modélisation de la structure, sans modélisation du domaine fluide.

Les méthodes push-over ont bien évidemment été construites et justifiées en faisant des hypothèses fortes de linéarisation du problème (petits déplacements, comportement élastique, flambement d’Euler…). Ci-dessous sont données les définitions des champs de pressions imposés [bib2] :

# ----------------------------------------------------------------------

# DEFINITION DES CHAMPS DE PRESSION (PA)

#

# PH(Z) : PRESSION HYDROSTATIQUE

# PIF(Z,TETA) : PRESSION IMPULSIVE FLEXIBLE

# PIR(Z,TETA) : PRESSION IMPULSIVE RIGIDE

# PV(Z) : PRESSION VERTICALE

#

# PIF(Z,TETA) = PIF(Z)*COS(TETA)

# PIR(Z,TETA) = PIR(Z)*COS(TETA)

#

#P(Z,TETA)=PH(Z)+INSTANT*( PIF(Z,TETA) + PIR(Z,TETA) +/- 0.4*PV(Z) )

# ----------------------------------------------------------------------

PHZ = DEFI_FONCTION ( NOM_PARA = 'Z' ,

…

VERIF = 'CROISSANT' , )

PIFZ = DEFI_FONCTION ( NOM_PARA = 'Z' ,

…

VERIF = 'CROISSANT' , )

PIRZ = DEFI_FONCTION ( NOM_PARA = 'Z' ,

…

VERIF = 'CROISSANT' , )

PVZ = DEFI_FONCTION ( NOM_PARA = 'Z' ,

…

VERIF = 'CROISSANT' , )

# ----------------------------------------------------------------------

# P(X,Y,Z) : PRESSION TOTALE

#P(Z,TETA)=PH(Z)+INSTANT*( PIF(Z,TETA) + PIR(Z,TETA) +/- 0.4*PV(Z) )

# ----------------------------------------------------------------------

COSTE=FORMULE(NOM_PARA=('X','Y'),VALE='X/SQRT((X*X)+(Y*Y))')

PH=FORMULE(NOM_PARA=('X','Y','Z'),VALE='PHZ(Z)')

PIF=FORMULE(NOM_PARA=('X','Y','Z'),VALE='PIFZ(Z)*COSTE(X,Y)')

PIR=FORMULE(NOM_PARA=('X','Y','Z'),VALE='PIRZ(Z)*COSTE(X,Y)')

PV=FORMULE(NOM_PARA=('X','Y','Z'),VALE='PVZ(Z)')

PS1=FORMULE(NOM_PARA=('X','Y','Z'),

VALE='PIF(X,Y,Z)+PIR(X,Y,Z)+(0.4*PV(X,Y,Z))')

PS2=FORMULE(NOM_PARA=('X','Y','Z')

VALE='PIF(X,Y,Z)+PIR(X,Y,Z)-(0.4*PV(X,Y,Z))')

On peut alors définir les chargement mécaniques correspondants (pesanteur et pressions suiveuses) :

PESA = AFFE_CHAR_MECA( MODELE = MODELE ,

Manuel d'utilisation

Document diffusé sous licence GNU FDL (http://www.gnu.org/copyleft/fdl.html)

Fascicule u2.06 : Dynamique

Code_Aster

Titre :

Analyse de la tenue sismique des grands réservoirs[...]

Responsable :

Nicolas GREFFET

Version default

Date :

21/04/2009

Page :

8/20

Clé :

U2.06.11

Révision :

931

PESANTEUR = ( 9.81 , 0. , 0. , -1. ) , )

#

PRESPH = AFFE_CHAR_MECA_F( MODELE = MODELE ,

FORCE_COQUE = _F ( GROUP_MA = 'VIROLE',),PRES = PH, PLAN = 'INF',),)

#

PRESPS1 = AFFE_CHAR_MECA_F( MODELE = MODELE ,

FORCE_COQUE = _F ( GROUP_MA = 'VIROLE' , ) ,

PRES = PS1 ,

PLAN = 'INF' , ) , )

#

PRESPS2 = AFFE_CHAR_MECA_F( MODELE = MODELE ,

FORCE_COQUE = _F ( GROUP_MA =

( 'VIROLE' , ) ,

PRES = PS2 ,

PLAN = 'INF' , ) , )

Ces champs de pressions sont multipliés par une fonction croissante linéaire du temps :

FONCMUL = DEFI_FONCTION( NOM_PARA = 'INST' , = ( 0. , 0. , 3. , 3. ) , )

La liaison avec le sol est considérée ici comme étant complète (groupe de noeuds 'BASE'). On définit aussi les conditions de symétrie (groupe de nœuds 'SYMETRIE') puisqu’on ne maille qu’un demi réservoir :

CONDLIM = AFFE_CHAR_MECA( MODELE = MODELE ,

DDL_IMPO = ( _F ( GROUP_NO = 'BASE' ,

DX = 0., DY = 0., DZ = 0.,

DRX = 0., DRY = 0., DRZ = 0.,) ,

_F ( GROUP_NO = 'SYMETRIE' ,

DY = 0., DRX = 0., DRZ = 0.,),) , )

3.2

Méthode de résolution quasistatique monotone non linéaire

On veut résoudre un problème d’évolution quasistatique non linéaire. On va donc utiliser l’opérateur

STAT_NON_LINE ([U4.51.03], [R5.03.01]). Le chargement imposé sera construit avec la pression

PRESPS1 par exemple :

RESU = STAT_NON_LINE( MODELE = MODELE ,

CHAM_MATER = CHMAT ,

CARA_ELEM = CARAELEM ,

EXCIT = ( _F ( CHARGE = CONDLIM ,

TYPE_CHARGE = 'FIXE_CSTE' , ) ,

_F ( CHARGE = PESA ,

TYPE_CHARGE = 'FIXE_CSTE' , ) ,

_F ( CHARGE = PRESPH ,

FONC_MULT = FONCMUL ,

TYPE_CHARGE = 'SUIV' , ) ,

_F ( CHARGE = PRESPS1 ,

FONC_MULT = FONCMUL ,

TYPE_CHARGE = 'SUIV' , ) , ) ,

COMP_INCR = ( _F ( RELATION = 'ELAS' ,

COQUE_NCOU = 1 ,

DEFORMATION = 'GREEN_GR' ,

GROUP_MA = ( 'SURF2' ,

'ANNEAU' , 'TFC2' ) , ) ,

_F ( RELATION = 'VMIS_ISOT_TRAC' ,

COQUE_NCOU = 1,

DEFORMATION = 'GREEN_GR' ,

GROUP_MA =

( 'SURF0' , 'SURF1' , ) , ) , ) ,

INCREMENT = _F ( LIST_INST = L_INST1 ,

Manuel d'utilisation Fascicule u2.06 : Dynamique

Document diffusé sous licence GNU FDL (http://www.gnu.org/copyleft/fdl.html)

Code_Aster

Titre :

Analyse de la tenue sismique des grands réservoirs[...]

Responsable :

Nicolas GREFFET

Version default

Date :

21/04/2009

Page :

9/20

Clé :

U2.06.11

Révision :

931

SUBD_PAS = 4 ,

SUBD_PAS_MINI = 1.E-4 ,

SUBD_METHODE='UNIFORME' ) ,

…

)

3.3

Remarques sur les calculs et post-traitements

Dans l’exemple présenté ci-dessus on utilise un algorithme de Newton : on résout avec l’opérateur tangent global (raideur) réactualisé à chaque itération. Dans le cas où le problème est bien posé

(suffisamment « régulier »), il est connu que ce type d’algorithme offre la meilleure convergence. Donc, dans notre cas, tant que l’opérateur tangent global est « bien » défini (loin d’être singulier), le calcul va bien se dérouler avec une convergence rapide. Quand le niveau de charge imposé va se rapprocher de la charge ultime, la structure va alors devenir instable au sens du flambage [R7.05.01], ce qui se traduit par la tendance de l’opérateur tangent global à devenir singulier. La perte de stabilité par point limite est en fait la perte d’unicité de la solution, soit donc la singularité de l’opérateur de résolution. Au voisinage du point limite l’algorithme de Newton convergera moins bien, d’où la nécessité d’imposer des incréments de temps plus petits et une augmentation du nombre d’itérations sur le résidu en

équilibre.

D’une manière générale, plus l’on va s’approcher de la charge ultime, plus le pas de temps devra être petit. Malgré cela les risques d’arrêt du calcul sur non convergence sont importants, d’où l’obligation de mener le calcul suivant plusieurs poursuites successives.

Il est néanmoins possible d’améliorer cette convergence en changeant en cours de calcul d’algorithme et de basculer sur un quasi-Newton. Pour cela, il suffit de résoudre sur l’opérateur tangent que l’on ne réactualise qu’à chaque pas (entre deux itération on garde le même), et si cela est insuffisant, on peut alors résoudre avec l’opérateur élastique qui, lui, ne peut devenir singulier. Ce choix renforce la robustesse de l’algorithme en termes de convergence, mais il augmente, parfois considérablement, le nombre d’itérations (et/ou de pas) nécessaires pour obtenir la solution.

Pour notre type d’étude on peut distinguer deux types de quantités d’intérêt pour le post-traitement.

D’une part une grandeur apte à traduire le flambement de la structure et donc à faire apparaître la charge ultime. Pour cela on peut tracer le niveau de pression (coefficient multiplicateur) fonction du déplacement d’un point situé au sommet du réservoir (à l’extrémité de la génératrice la plus soumise à de la compression). D’autre part un indicateur plus local de l’apparition de la plasticité : l’isovaleur de la déformation plastique cumulée, à chaque instant de calcul.

Ces deux post-traitements ne présentent aucune difficulté particulière dans le Code_Aster.

Les modèles éléments finis mis en œuvre comportaient entre 55000 et 110000 ddl, menant à un temps de calcul total de l’ordre de 5 h à 10 h sur Compaq AlphaServer. Ces modèles se basaient sur une représentation simplifiée des ancrages : la géométrie des renforts et goussets est maillée finement, mais le boulonnement n’est pas présent et il est remplacé par une condition d’encastrement sur tous les nœuds d’altitude 0 m. Une modélisation plus réaliste des ancrages, avec décollement aurait donc entraîné une taille de problème global sensiblement plus grande.

La méthodologie de calcul présentée ici, dont l’objectif est d’étudier la réponse quasistatique en flambage [R7.05.01] d’un réservoir est très proche du cadre de la notice de calcul au flambage

[U2.08.04]. Cette documentation présente les analyses de stabilité linéaire (au sens d’Euler) et la mise au point d’un calcul non linéaire de type push-over pour obtenir la charge ultime.

3.4

Modélisation fine des ancrages : soulèvement

Le type de réservoirs étudiés ici sont boulonnés au sol [bib5]. Ces boulons (ou tirants) traversent des brides renforcées soudées en base de virole). Pour plus de rapidité de calcul nous présentons par la suite la méthode mise en œuvre sur un maillage simplifié, mais avec des ancrages réalistes : le réservoir est fixé au sol par 18 tirants.

Dans le Code_Aster, plusieurs modélisations du contact sont disponibles.

Manuel d'utilisation

Document diffusé sous licence GNU FDL (http://www.gnu.org/copyleft/fdl.html)

Fascicule u2.06 : Dynamique

Code_Aster

Titre :

Analyse de la tenue sismique des grands réservoirs[...]

Responsable :

Nicolas GREFFET

Version default

Date :

21/04/2009

Page :

10/20

Clé :

U2.06.11

Révision :

931

La modélisation de la géométrie de la zone de contact peut être surfacique (problème 3D), linéique

(problème 2D) ou constituée d’éléments discrets (mot-clés DIS_CHOC dans STAT_NON_LINE et

DIS_CONTACT pour la loi matériau).

Le contact lui-même peut être traité, soit de façon nodale par pénalisation ou par la méthode des lagrangiens, avec contraintes actives ou non, soit de façon continue par la méthode des lagrangiens augmentés.

Manuel d'utilisation

Document diffusé sous licence GNU FDL (http://www.gnu.org/copyleft/fdl.html)

Fascicule u2.06 : Dynamique

Code_Aster

Titre :

Analyse de la tenue sismique des grands réservoirs[...]

Responsable :

Nicolas GREFFET

Version default

Date :

21/04/2009

Page :

11/20

Clé :

U2.06.11

Révision :

931

La méthode la plus simple est celle de la pénalisation, qui, pour notre problème, présente aussi l’avantage de pouvoir tenir compte du rôle du joint d’étanchéité sans avoir à le mailler séparément.

Pour un calcul quasistatique, la pénalisation n’est pas numériquement problématique, comme en dynamique où cette technique peut engendrer des perturbations hautes fréquences (liées à la raideur de pénalisation). De plus, pour des éléments de structure, on peut caler la raideur de pénalisation de manière à approcher au mieux la raideur de contact de la pièce qui est en réalité massive.

Les autres méthodes sont cependant plus rigoureuses car elles n’engendrent pas d’interpénétration.

Pour l’étude du décollement de la bâche, nous allons donc utiliser une méthode de pénalisation, qui présente le meilleur rapport qualité de modélisation / coût de calcul sur notre cas précis.

Étant donnée l’option prise d’obliger tous les nœuds à ne pouvoir se déplacer que verticalement, on peut profiter de cet appariement total pour se contenter d’utiliser des éléments discrets pour le contact.

En effet, il n’est pas nécessaire de faire du réappariement, et donc, les modélisations générales surfaciques maître-esclave des zones de contact sont inutilement lourdes. On va donc placer un tapis d’éléments discrets de contact sous la frette basse. On aura alors un élément discret (DIS_T sur maille SEG2) sous chaque nœud d’altitude 0.

Si la solution calculée s’avère trop éloignée de la solution de référence expérimentale, il sera nécessaire de quantifier l’influence de la condition de déplacement uniquement vertical. Il suffira de relancer des calculs sans cette hypothèse cinématique, avec, si les déplacements horizontaux sont grands, une méthode de gestion du contact avec réappariement des nœuds des surfaces concernées.

Remarques

Rigoureusement, il faudrait doter tous les nœuds de la virole de base d’un élément discret de

contact. Or, si l’on utilise des éléments de type COQUE_3D pour les frettes, on ne peut placer

d’élément de contact au niveau des nœuds au centre des éléments. En effet, ce nœud présente la particularité de ne porter que des degrés de liberté de type rotation. La condition de contact, qui porte sur le déplacement normal à la face, ne peut donc être exprimée en ces nœuds. Les

éléments de contact, ainsi que la condition cinématique de déplacement vertical des nœuds de

base ne doivent donc être portées que par les nœuds sommets ou milieux des arêtes.

La modélisation du contact par éléments discrets permet, dans le Code_Aster, de définir directement la précharge de ces éléments. On peut ainsi représenter la précontrainte des tirants, sans avoir à prédéformer ces éléments, comme on devrait le faire si l’on avait choisi une autre modélisation du contact s’appuyant sur les surfaces en vis-à-vis.

Ci-dessous, nous présentons la position des éléments discrets de contact pénalisé (en gris clair) et des 9 ressorts modélisant les tirants (en noir). Le renforcement en base de virole est indiqué en pointillés :

Figure 3.4-a : Disposition des éléments discrets pour la liaison avec la dalle

Manuel d'utilisation

Document diffusé sous licence GNU FDL (http://www.gnu.org/copyleft/fdl.html)

Fascicule u2.06 : Dynamique

Code_Aster

Titre :

Analyse de la tenue sismique des grands réservoirs[...]

Responsable :

Nicolas GREFFET

Version default

Date :

21/04/2009

Page :

12/20

Clé :

U2.06.11

Révision :

931

Ces éléments discrets pour les ancrages sont groupés dans 2 groupes de mailles :

• RESSC : éléments de contact en gris clair (un élément SEG2 par nœud de la bride basse, soit 57 SEG2),

•

RESS : tirants en noir (9 ressorts par éléments SEG2).

On introduit donc un matériau supplémentaire (MATRES) qui correspond aux éléments de contact et qui permet de prendre en compte la précontrainte par le serrage des écrous sur les tirants (on définit aussi la distance associée à la hauteur entre brides avec DIST_1 et DIST_2) :

Eboulon = 2.E11 ;

Dboulon = 0.026 ;

Sboulon = 3.14159 * Dboulon * Dboulon / 4 ; kressort = Eboulon * Sboulon ; prechar = DEFI_CONSTANTE ( VALE= -1000.,) kpenal = kressort * 100.

MATRES = DEFI_MATERIAU( DIS_CONTACT = _F( RIGI_NOR = kpenal ,

EFFO_N_INIT = prechar ,

DIST_1 = 0.05, DIST_2 = 0.05,) , );

On complète l’affectation matériau avec ces éléments discrets sur le groupe de mailles RESSC :

CHMAT = AFFE_MATERIAU( MAILLAGE = MAILLA2 ,

AFFE = ( _F ( GROUP_MA = 'VIROLE' , ) ,

MATER = MAA240 , ) ,

_F ( GROUP_MA = ( 'ANNET' , ) ,

MATER = MABID , ) ,

_F ( GROUP_MA = ( 'ANNEA' , ) ,

MATER = MAA240 , ) ,

_F ( GROUP_MA = ( 'RESSC' , ) ,

MATER = MATRES , ) , ) , )

Les caractéristiques élémentaires des nouveaux éléments discrets sont définies ainsi :

CARAELEM = AFFE_CARA_ELEM( MODELE = MODELE ,

DISCRET = (_F( GROUP_MA = 'RESSC' ,

REPERE = 'LOCAL', CARA = 'K_T_D_L',

VALE = ( 10. , 0.0 , 0.0 ,) ,) ,

_F( GROUP_MA = 'RESS' ,

REPERE = 'LOCAL', CARA = 'K_T_D_L',

VALE = ( kressort , 0., 0.,),),),)

Les ressorts qui modélisent les tirants (GROUP_MA RESS) ont bien la raideur équivalente kressort.

Il reste à modifier les conditions aux limites en déplacement pour les ancrages : on n’autorise que les déplacements suivant la verticale des nœuds de la face inférieure de la bride basse (nœuds initialement en contact avec le sol).

L’opérateur STAT_NON_LINE ([U4.51.03], [R5.03.01]) voit aussi ses arguments impactés par l’introduction du soulèvement en base (mot-clé DIS_CHOC pour le groupe de maille RESSC) :

RESU = STAT_NON_LINE( MODELE = MODELE , _MATER = CHMAT ,

…

C OMP_INCR = ( _F( RELATION = 'ELAS' ,

DEFORMATION = 'GREEN_GR' ,

GROUP_MA = 'ANNEAU' , ) ,

_F( RELATION = 'ELAS' ,

GROUP_MA = 'RESS' , ) ,

_F( RELATION = 'DIS_CHOC' ,

GROUP_MA = 'RESSC' , ) ,

_F( RELATION = 'VMIS_ISOT_TRAC' ,

DEFORMATION = 'GREEN_GR' ,

GROUP_MA = ( 'SURF0' , 'ANNEA' ,

'SURF1' , 'SURF2' ) , ) , ) ,

…

);

Manuel d'utilisation Fascicule u2.06 : Dynamique

Document diffusé sous licence GNU FDL (http://www.gnu.org/copyleft/fdl.html)

Code_Aster

Titre :

Analyse de la tenue sismique des grands réservoirs[...]

Responsable :

Nicolas GREFFET

Le reste du fichier de commande est inchangé.

Version default

Date :

21/04/2009

Page :

13/20

Clé :

U2.06.11

Révision :

931

Manuel d'utilisation

Document diffusé sous licence GNU FDL (http://www.gnu.org/copyleft/fdl.html)

Fascicule u2.06 : Dynamique

Code_Aster

Titre :

Analyse de la tenue sismique des grands réservoirs[...]

Responsable :

Nicolas GREFFET

Version default

Date :

21/04/2009

Page :

14/20

Clé :

U2.06.11

Révision :

931

3.5

Utilisation du critère de stabilité non linéaire

On peut également utiliser un critère de stabilité basé sur la matrice tangente : on a une instabilité si la matrice de raideur tangente devient singulière, i.e. si au moins une de ses valeurs propres s’annule.

On résout alors le problème aux valeurs propres suivant, écrit en grands déplacements (écriture en lagrangien avec le tenseur de déformation de Green-Lagrange) [R7.05.01] :

(

K

T

+

λ

I d

)

x

=

0

⇔

K

T

x

=

λ

I d

x

Avec





K





K

K

T

L

(

=

u

)

K

+

K

L

(

u

)

+

K

Q

(

: partie linéaire en

u

)

+

K

(

π

) : matrice de rigidité tangente,

u

de la matrice

K

T

,





K

π





I



λ

d

Q

π

:

(

u

)

: partie quadratiqu matrice identité,

: valeur propre.

e en

u

: matrice de raideur géométriqu e,

: le tenseur de Piola Kirchhoff de la matrice

II,

K

T

,

La documentation [R7.05.01] présente ces analyses de stabilité plus en détail.

Remarque

Lorsque les déplacements sont petits, on a simplement

K

(

π

)

=

K

(

σ

)

=

K

G

et les matrices

K

L

et

K

Q

peuvent être négligées.

Si l’on veut se servir du critère de stabilité, il suffit de rajouter la ligne suivante parmi les arguments de

STAT_NON_LINE ([U4.51.03], [R5.03.01]) :

CRIT_FLAMB=_F(NB_FREQ=1,),

On cherche alors la première valeur propre de l’opérateur global tangent de notre système.

Si au cours du calcul on observe que cette valeur propre diminue, voire change de signe, cela signifie que l’on s’est approché de la première charge critique et qu’ensuite on l’a même dépassée.

Le nombre de valeurs propres à déterminer peut être imposé par le mot-clé NB_FREQ (3 par défaut).

Il est également possible, en utilisant la commande CHAR_CRIT, de choisir la bande dans laquelle il faut chercher ces valeurs propres (de –10 à 10 par défaut).

Remarque

L’indication d’une bande de fréquence est utile surtout pour des calculs en petites perturbations où un test de Sturm est effectué pour la bande de fréquences fournie. On peut ainsi gagner du temps en ne calculant les valeurs propres que s’il y en a dans la bande indiquée. Le test de

Sturm n’est pas fait en grandes déformations et les valeurs propres sont calculées à chaque pas de temps.

Le mode de flambement ainsi que les valeurs propres déterminés peuvent être récupérés en utilisant la commande IMPR_RESU :

IMPR_RESU( MODELE = MODELE ,FORMAT = 'RESULTAT' ,RESU=(_F(RESULTAT=RESU,

NOM_PARA='CHAR_CRIT',’MODE_FLAMB’,),),)

Le cas-test SSLL105 [V3.01.105] propose un exemple d’utilisation de ce critère de stabilité pour un cas linéaire et le cas-test SSNL126 [V6.02.126] pour un cas non linéaire (poutre élastoplastique).

Manuel d'utilisation

Document diffusé sous licence GNU FDL (http://www.gnu.org/copyleft/fdl.html)

Fascicule u2.06 : Dynamique