Code_Aster
Titre :
Analyse de la tenue sismique des grands réservoirs[...]
Responsable :
Nicolas GREFFET
Version default
Date :
21/04/2009
Page :
15/20
Clé :
U2.06.11
Révision :
931
3.6
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. Le pilotage ne peut être utilisé avec le contact [U4.51.03].
4 Approche transitoire couplée
4.1
Problème de référence
On sort ici du cadre réglementaire et l’on va exploiter toutes les possibilités de modélisation offertes par le
Code_Aster. Le modèle du réservoir lui-même reste inchangé (coques volumiques élastoplastiques). En revanche, on va représenter le domaine fluide par un maillage massif. De plus, la résolution se fera en dynamique transitoire avec l’opérateur DYNA_NON_LINE ([U4.53.01] , [R5.05.05]), la sollicitation externe étant du type sismique.
Le domaine fluide est modélisé en acoustique linéaire (barotrope, compressible, non visqueux et avec surface libre). Le problème couplé fluide-structure est résolu dans le Code_Aster par une formulation symétrique (u, p,
φ
) ([R4.02.02], [bib7]), en écriture lagrangienne réactualisée. Le chargement est du type accélérogramme imposé en base de bâche.
Le problème discrétisé se présente alors ainsi [bib6] :
Figure 4.1-a : Maillage complet avec le domaine fluide et la structure
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 :
16/20
Clé :
U2.06.11
Révision :
931
4.2
Modélisation couplée fluide-structure dans le Code_Aster
Afin de pouvoir mener un calcul couplé par formulation (u, p,
φ
) [bib7] dans le Code_Aster , avec surfaces libres et glissement à l’interface (non adhérence entre le fluide non visqueux et la paroi interne du réservoir), il faut respecter une certaine construction du maillage et des modèles correspondants.
On doit donc définir le maillage fluide, le maillage structure et l’interface fluide-structure.
Dans le mailleur on génère donc deux maillages distincts (mailles et nœuds différents pour permettre le glissement à l’interface fluide-structure) pour le domaine fluide et la structure.
Ensuite, dans le Code_Aster, on va générer le maillage support de l’interface comme suit, à partir des groupes de mailles SURF0 , SURF1 et SURF2 qui sont les viroles du réservoir :
MAILLA01=CREA_MAILLAGE(MAILLAGE=MAILLA1,
CREA_GROUP_MA=_F(NOM='IFLUSTRU',
GROUP_MA=('SURF0','SURF1','SURF2',),
PREF_MAILLE='I',),);
On crée donc un nouveau groupe de mailles IFLUSTRU . Les nœuds de ce maillage sont confondus avec ceux de la structure, mais les mailles sont différentes.
Remarques
La structure est maillée en éléments COQUE_3D , dont le support géométrique est un quadrilatère
à 9 nœuds. Le nœud milieu présente la particularité de ne porter que des ddl de rotation.
Si l’on veut coupler un tel élément avec un élément massif de fluide, on ne peut donc écrire de condition cinématique de couplage pour ce nœud milieu puisqu’il ne comporte pas de ddl de translation, contrairement au nœud correspondant qui vient du domaine fluide et qui, lui, ne porte que des ddl de translation.
Pour contourner ce problème, on ne va écrire le couplage fluide-structure que sur les nœuds du maillage coque comportant les ddl de translation : les nœuds sommets et les nœuds milieux des arêtes des éléments.
Il faut donc avoir dans le Code_Aster le maillage structure ne comportant que des éléments quadrilatères à 8 nœuds (ce sont bien les nœuds sommets et les nœuds milieux des arêtes), sur
lequel on définit l’interface. Le maillage structure pour les COQUE_3D étant défini seulement
après, à partir de ce maillage, en rajoutant les nœuds milieux.
Le maillage fluide, pour être conforme avec l’interface, sera lui composé de parallélépipèdes à 20 nœuds.
En outre, même si l’on utilisait un maillage structure avec des éléments de coque dont le support géométrique serait un quadrilatère à 9 nœuds, et où tous les nœuds, même le central portaient des ddl de translation, le couplage fluide-structure pourrait poser un problème. En effet, le maillage fluide massif conforme devrait être réalisé avec des éléments massifs parallélépipédiques à 27 nœuds. Or, certains mailleurs n’offrent pas ce type d’éléments complets qui sont assez peu employés en calcul de structure, contrairement aux éléments quadratiques classiques que sont les parallélépipèdes à 20 nœuds.
Une fois la définition de l’interface, qui est donc maillée en éléments quadrilatères à 8 nœuds (qui sont
équivalents aux faces des éléments massifs employés pour le domaine fluide : éléments parallélépipède à 20 nœuds), on peut donc faire la modification du maillage structure (groupe de mailles RESERVOI) en éléments à 9 nœuds, support géométrique des COQUE_3D :
MAILLA2=CREA_MAILLAGE(MAILLAGE=MAILLA01,
MODI_MAILLE=_F(GROUP_MA='RESERVOIR',
OPTION='QUAD8_9',
PREF_NOEUD='NSQ',
PREF_NUME=1,),);
On peut alors définir les modèles nécessaires au calcul couplé :
MODELE=AFFE_MODELE(MAILLAGE=MAILLA2, INFO=1, VERIF='MAILLE',
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 :
17/20
Clé :
U2.06.11
Révision :
931
AFFE=(_F(GROUP_MA='SURFLIBR',
PHENOMENE='MECANIQUE',
MODELISATION='2D_FLUI_PESA',),
_F(GROUP_MA=('FLUID0','FOND','PLANCENT',),
PHENOMENE='MECANIQUE',
MODELISATION='3D_FLUIDE',),
_F(GROUP_MA='IFLUSTRU',
PHENOMENE='MECANIQUE',
MODELISATION='FLUI_STRU',),
_F(GROUP_MA='RESERVOI',
PHENOMENE='MECANIQUE',
MODELISATION='COQUE_3D',),),);
Le groupe de mailles SURFLIBR (mailles de bord du domaine fluide, situées sur sa face supérieure) porte un modèle de surface libre [R4.02.04] : 2D_FLUID_PESA .
Les groupes de mailles FLUID0 (domaine fluide massif), FOND (mailles de bord de FLUID0 définissant le fond) et PLANCENT (mailles de bord de FLUID0 définissant le plan de symétrie) définissent le domaine fluide total :
3D_FLUIDE .
L’interface fluide-structure ( FLUI_STRU ) est portée par le groupe de mailles IFLUSTRU .
Enfin, le maillage du réservoir ( RESERVOI ) est le support du modèle structure en coques volumiques (
COQUE_3D ).
On définit aussi le matériau compressible fluide :
EAU=DEFI_MATERIAU(FLUIDE=_F(RHO=1000.0, CELE_R=1500.0,),); que l’on affectera au domaine fluide (modèle massif et ses bords) ainsi qu’à l’interface IFLUSTRU et à la surface libre SURFLIBR :
CHMAT=AFFE_MATERIAU(MAILLAGE=MAILLA2,
AFFE=(_F(GROUP_MA=('FLUID0','FOND','PLANCENT','IFLUSTRU','SURFLIBR',),
MATER=EAU,),
…
) ;
4.3
Conditions aux limites
Les conditions aux limites cinématiques portent sur l’encastrement en base de réservoir ( SYMETRI2 ), sur les génératrices dans le plan vertical de symétrie ( SYMETRIE ) et sur la non pénétration (vitesse normale nulle) au fond du domaine fluide ( FOND ), ainsi que dans son plan vertical de symétrie (PLANCENT).
CONDLIM=AFFE_CHAR_MECA(MODELE=MODELE,
DDL_IMPO=(_F(GROUP_NO='SYMETRI2',
DX=0.0, DY=0.0, DZ=0.0,
DRX=0.0, DRY=0.0, DRZ =0.0,),
_F(GROUP_NO='SYMETRIE',
DY=0.0, DRX=0.0, DRZ=0.0,),),
VITE_FACE=_F(GROUP_MA=('FOND','PLANCENT',), VNOR=0.0,),);
Les effets de pesanteur sont pris en compte :
PESA=AFFE_CHAR_MECA(MODELE=MODELE, PESANTEUR=(9.81,0.0,0.0,-1.0,),);
Cependant cette commande est insuffisante car dans la modélisation (u, p,
φ
) [bib7] les effets de pesanteur dans le domaine fluide ne peuvent être pris en compte. Si on ne faisait rien de plus, la pesanteur ne serait donc véritablement imposée que sur la structure.
Pour approcher les effets de pesanteur du fluide sur la paroi, on va imposer une pression hydrostatique équivalente (mais qui ne peut tenir compte de la variation de hauteur du fluide au cours du calcul quand le ballottement va commencer) :
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 :
18/20
Clé :
U2.06.11
Révision :
931
PHZ=DEFI_FONCTION(NOM_PARA='Z',
…
VERIF='CROISSANT',);
#
PH = FORMULE(NOM_PARA=('X','Y','Z'),VALE='PHZ(Z)',),
#
PRESSHYD=AFFE_CHAR_MECA_F(MODELE=MODELE,
FORCE_COQUE=_F(GROUP_MA='VIROLE', PRES= PH ,
PLAN='INF',),);
Le séisme est imposé comme étant un accélérogramme (fonction GASDM_X1) imposé en base suivant la direction X. C’est donc une sollicitation de type mono-appui classique :
ACCELERX=CALC_FONCTION(COMB=(_F(FONCTION=GASDM_X1, COEF=0.5,),),);
#
MULT_X=CALC_CHAR_SEISME(MATR_MASS=MATMAS,
DIRECTION=(1.0,0.0,0.0,), _APPUI='OUI',);
#
CHARG_SE=AFFE_CHAR_MECA(MODELE=MODELE, VECT_ASSE=MULT_X,);
La matrice de masse utilisée ( MATMAS ) est la matrice de masse du système couplé total.
4.4
Conditions initiales
L’état initial du calcul transitoire doit correspondre à l’équilibre du système total lorsqu’il n’est pas soumis au séisme. Cet état d’équilibre statique correspond donc au chargement de pesanteur et aux effets hydrostatiques.
Si l’on commençait le calcul dynamique avec un état initial ne respectant pas cet équilibre, alors cela génèrerait des oscillations de la solution transitoire puisqu’elle ne serait pas initialement à l’équilibre (le niveau sismique
étant cependant alors nul). On peut atténuer ces oscillations en ajoutant un amortissement structural « grand » et en attendant que la solution se stabilise avant d’imposer le séisme, mais cette technique est assez peu
élégante…
Pour calculer cet état initial statiquement équilibré, on peut donc résoudre un problème statique (que l’on suppose en plus linéaire) d’équilibre sous l’action des forces de pesanteur et hydrostatiques.
Pour ce faire, on calcule et assemble préalablement les matrices globales K et M avec le chargement hydrostatique et la pesanteur :
MACRO_MATR_ASSE( MODELE = MODELE, CHAM_MATER=CHMAT,
CARA_ELEM=CARAELEM, CHARGE = (CONDLIM,PESA,PRESSHYD,),
NUME_DDL = CO('NUMSTA'),
SOLVEUR=_F(METHODE='MULT_FRONT',RENUM='METIS'),
MATR_ASSE=(
_F( MATRICE= CO('RIGSTA'), OPTION= 'RIGI_MECA'),
_F( MATRICE= CO('MASSTA'), OPTION= 'MASS_MECA'), ),);
La matrice de raideur assemblée étant singulière à cause du domaine fluide (la formulation (u, p,
φ
) rend cette matrice singulière à fréquence nulle [bib8]), on modifie légèrement le problème en considérant la matrice de raideur
K
cor
=
K
+
ε
M
,
ε
<<
1 qui n’est plus singulière (on la nomme RIGICOMB).
On peut, par exemple, prendre
ε
= 0,001 comme ci-dessous :
RIGICOMB = COMB_MATR_ASSE(COMB_R=(_F( MATR_ASSE= RIGSTA, COEF_R= 1.),
_F( MATR_ASSE= MASSTA, COEF_R= -0.001),);
On assemble ensuite le vecteur chargement F_0 (second membre) :
E_0 = CALC_VECT_ELEM(CARA_ELEM=CARAELEM, CHAM_MATER=CHMAT,
OPTION='CHAR_MECA',CHARGE=(CONDLIM, PESA,PRESSHYD,),);
F_0 = ASSE_VECTEUR(VECT_ELEM= E_0, NUME_DDL= NUMSTA );
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 :
19/20
Clé :
U2.06.11
Révision :
931
On peut alors résoudre le problème de statique linéaire
K
cor
factorisation de type LDL T :
U
=
F
0
, en utilisant, par exemple, une
RIGICOMB = FACTORISER(reuse=RIGICOMB, MATR_ASSE= RIGICOMB,
STOP_SINGULIER= 'NON');
DEP0 = RESOUDRE(MATR_FACT= RIGICOMB, CHAM_NO=F_0);
Le champ solution DEP0 sera donc l’état initial du calcul dynamique transitoire qui suit.
4.5
Résolution transitoire
On utilise l’opérateur DYNA_NON_LINE ([U4.53.01] , [R5.05.05]) comme suit :
RESU=DYNA_NON_LINE(MODELE=MODELE,
CHAM_MATER=CHMAT,
CARA_ELEM=CARAELEM,
EXCIT=(_F(CHARGE=CONDLIM,),
_F(CHARGE=PESA,),
_F(CHARGE=PRESSHYD,
FONC_MULT=FONCMUL0,
TYPE_CHARGE='SUIV',),
_F(CHARGE=CHARG_SE,
FONC_MULT=ACCELERX,),),
COMP_INCR=(_F(RELATION='ELAS',
DEFORMATION='PETIT_REAC',
GROUP_MA=('FLUID0','FOND','PLANCENT',
'IFLUSTRU','SURFLIBR',),),
_F(RELATION='ELAS', DEFORMATION='PETIT_REAC',
GROUP_MA='ANNEAU',),
_F(RELATION='VMIS_ISOT_TRAC',
DEFORMATION='PETIT_REAC',
GROUP_MA=('SURF0','SURF1',
'SURF2','SURF3',),),),
ETAT_INIT=_F(INST_ETAT_INIT=0.0, DEPL=DEP0,),
INCREMENT=_F(LIST_INST=LINST,
SUBD_PAS=4,
SUBD_PAS_MINI=1e-06,
SUBD_METHODE='UNIFORME', ),
SCHEMA_TEMPS=_F(SCHEMA='HHT', ALPHA=-0.1,
FORMULATION='DEPLACEMENT',),
NEWTON=_F(REAC_INCR=1, MATRICE='TANGENTE',REAC_ITER=1,),
SOLVEUR=_F(STOP_SINGULIER='NON',),
CONVERGENCE=_F(RESI_GLOB_RELA=1.e-05, ITER_GLOB_MAXI=20,
ARRET='OUI',),
ARCHIVAGE=_F(LIST_INST=LARCH, ARCH_ETAT_INIT='OUI',),);
La résolution se fait en lagrangien réactualisé (option DEFORMATION='PETIT_REAC') car le domaine fluide est en petites perturbations sur chaque pas. Il faut donc vérifier que le pas de temps est suffisamment petit pour que cette hypothèse soit vérifiée.
On utilise un schéma d’intégration en temps de type accélération moyenne modifiée (SCHEMA='HHT',
ALPHA=-0.1) avec amortissement numérique afin de stabiliser la solution et de faciliter la convergence.
Le cas-test FDNV100 [V8.03.100] présente le calcul d’une cuve rectangulaire pleine d’eau avec une paroi souple. La modélisation mise en œuvre est très proche de celle utilisée ici pour les grands réservoirs.
Manuel d'utilisation
Document diffusé sous licence GNU FDL (http://www.gnu.org/copyleft/fdl.html)
Fascicule u2.06 : Dynamique
