Approche quasi-statique push-over ================================= .. _Ref121224245: Chargements imposés ------------------- Le chargement imposé est d’origine réglementaire (EC-8), quasi-statique [:ref:`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 quasi-statique 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 [:ref:`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 , PESANTEUR=_F(GRAVITE=9.81, DIRECTION=(0.0,0.0,-1.0),),); # 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 nœuds '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.,),) , ) Méthode de résolution quasi-statique monotone non linéaire ---------------------------------------------------------- On veut résoudre un problème d’évolution quasi-statique non linéaire. On va donc utiliser l’opérateur STAT_NON_LINE ([:external:ref:`U4.51.03 `], [:ref:`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' , ) , ) , COMPORTEMENT = ( _F ( RELATION = 'ELAS' , COQUE_NCOU = 1 , DEFORMATION = 'GROT_GDEP' , GROUP_MA = ( 'SURF2' , 'ANNEAU' , 'TFC2' ) , ) , _F ( RELATION = 'VMIS_ISOT_TRAC' , COQUE_NCOU = 1, DEFORMATION = 'GROT_GDEP' , GROUP_MA = ( 'SURF0' , 'SURF1' , ) , ) , ) , INCREMENT = _F ( LIST_INST = L_INST1 , ) , ) 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 [:ref:`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’iso-valeur de la déformation plastique cumulée, à chaque instant de calcul. Ces deux post-traitements ne présentent aucune difficulté particulière dans *Code_Aster* . Les modèles éléments finis mis en œuvre comportaient entre 55000 et 110000 degrés de liberté. 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 :math:`0m` . 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 quasi-statique en flambage [:ref:`R7.05.01 `] d’un réservoir est très proche du cadre de la notice de calcul au flambage [:ref:`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. .. _Ref121226034: Modélisation fine des ancrages : soulèvement -------------------------------------------- Le type de réservoirs étudiés ici sont boulonnés au sol [:ref:`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 *Code_Aster* , plusieurs modélisations du contact sont disponibles. 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. La méthode la plus simple, dans cette configuration d'étude, est celle de la pénalisation, qui dans ce cas 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 quasi-statique, 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_3Dpour 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 *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 : +--------------------------------------------------------------------------------------------------------------+ | | | .. image:: images/10138D2E00004AD400002FA5C74001194C7E309E.svg | | :width: 403 | | :height: 252 | | | | | | .. _RefImage_10138D2E00004AD400002FA5C74001194C7E309E.svg: | | | | | +--------------------------------------------------------------------------------------------------------------+ | **Figure** :ref:`3.4 ` **-a: Disposition des éléments discrets pour la liaison avec la dalle** | +--------------------------------------------------------------------------------------------------------------+ 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 ; kpenal = kressort * 100. MATRES = DEFI_MATERIAU( DIS_CONTACT = _F( RIGI_NOR = kpenal , 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_MARESS) 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([:external:ref:`U4.51.03 `], [:ref:`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 , … COMP_INCR = ( _F( RELATION = 'ELAS' , DEFORMATION = 'GROT_GDEP' , GROUP_MA = 'ANNEAU' , ) , _F( RELATION = 'ELAS' , GROUP_MA = 'RESS' , ) , _F( RELATION = 'DIS_CHOC' , GROUP_MA = 'RESSC' , ) , _F( RELATION = 'VMIS_ISOT_TRAC' , DEFORMATION = 'GROT_GDEP' , GROUP_MA = ( 'SURF0' , 'ANNEA' , 'SURF1' , 'SURF2' ) , ) , ) , … ); Le reste du fichier de commande est inchangé. 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, soit 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) [:ref:`R7.05.01 `]: :math:`({K}_{T}+\lambda {I}_{d})x=0\mathrm{\iff }{K}_{T}x=\lambda {I}_{d}x` avec : .. csv-table:: ":math:`{K}_{T}=K+{K}^{L}(\underline{u})+{K}^{Q}(\underline{u})+K(\pi )`", "Matrice de rigidité tangente" ":math:`{K}^{L}(\underline{u})`", "Partie linéaire en :math:`\underline{u}` de la matrice :math:`{K}_{T}`" ":math:`{K}^{Q}(\underline{u})`", "Partie quadratique en :math:`\underline{u}` de la matrice :math:`{K}_{T}`" ":math:`K(\pi )`", "Matrice de raideur géométrique" ":math:`\pi`", "Tenseur de Piola-Kirchhoff :math:`\mathrm{II}`" ":math:`{I}_{d}`", "Matrice identité" ":math:`\lambda`", "Valeur propre" La documentation [:ref:`R7.05.01 `] présente ces analyses de stabilité plus en détail. La documentation [:ref:`U2.08.04 `] en présente l'utilisation. **Remarques** *Lorsque les déplacements sont petits, on a simplement* :math:`K(\pi )=K(\sigma )={K}_{G}` *et les matrices* :math:`{K}^{L}` *et* :math:`{K}^{Q}` *peuvent être négligées.* *L'analyse linéaire au sens d'Euler dans Code_Aster ne permet pas de tenir compte de l'aspect suiveur du champ de pression : il faut obligatoirement passer la l'approche présentée ici et basé* e sur *une résolution non-linéaire.* 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([:external:ref:`U4.51.03 `], [:ref:`R5.03.01 `]): CRIT_STAB=_F(NMAX_CHAR_CRIT=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é NMAX_CHAR_CRIT (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 (coefficient de charge critique) 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 [:external:ref:`V3.01.105 `] propose un exemple d’utilisation de ce critère de stabilité pour un cas linéaire et le cas-test SSNL126 [:external:ref:`V6.02.126 `] pour un cas non linéaire (poutre élastoplastique). Pilotage du chargement ---------------------- Afin de faciliter la convergence du calcul incrémental lorsque l’on est proche du niveau de charge ultime, ou afin de pouvoir dépasser ce point critique, il peut être judicieux de ne plus se placer en chargement imposé pour privilégier un pilotage en déplacement ou un pilotage par longueur d’arc (son usage pour le post-critique est brièvement rappelé dans la documentation [:ref:`U2.08.04 `]). Le pilotage ne peut être utilisé avec le contact [:external:ref:`U4.51.03 `].