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 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 Code_Aster par une formulation symétrique \((u,p,\varphi )\) ([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

4.2. Modélisation couplée fluide-structure dans Code_Aster

Afin de pouvoir mener un calcul couplé par formulation \((u,p,\varphi )\) [bib7] dans 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 différentes mais nœuds communs) pour le domaine fluide et la structure. Pour des raisons pratiques, il peut être plus simple de générer les mailles distinctes séparément, donc avec des nœuds dédoublés et ensuite d’éliminer ces nœuds doubles pour ne plus avoir que des nœuds simples.

Ensuite, dans Code_Aster , on va générer le maillage support de l’interface comme suit, à partir des groupes de mailles SURF0, SURF1et SURF2qui 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 degrés de liberté 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 degré de liberté de translation, contrairement au nœud correspondant qui vient du domaine fluide et qui, lui, ne porte que des degrés de liberté 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 degrés de liberté de translation: les nœuds sommets et les nœuds milieux des arêtes des éléments.

Il faut donc avoir dans 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 degrés de liberté 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.

Il est impératif que la normale extérieure au domaine fluide soit toujours orientée dans le même sens. Il est fortement conseillée de garder la convention d’orientation de la structure vers le fluide pour toutes les modélisations d’interface fluide-structure, en particulier FLUI_STRU.

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”,

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 FLUID0définissant le fond) et PLANCENT(mailles de bord de FLUID0dé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 IFLUSTRUet à 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=_F( GRAVITE=9.81,

DIRECTION=(0.0,0.0,-1.0),),);

Cependant cette commande est insuffisante car dans la modélisation \((u,p,\varphi )\) [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):

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:

ASSEMBLAGE( 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}_{\mathrm{cor}}=K+\varepsilon M\) \(\varepsilon \ll 1\) qui n’est plus singulière (on la nomme RIGICOMB).

On peut, par exemple, prendre \(\varepsilon =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 );

On peut alors résoudre le problème de statique linéaire \({K}_{\mathrm{cor}}\underline{U}=\underline{{F}_{0}}\) , en utilisant, par exemple, une factorisation de type \({\mathrm{LDL}}^{T}\) :

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,),),

COMPORTEMENT=(_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,),

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),);

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.

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

On peut également utiliser un critère de stabilité tout comme en quasi-statique. La présence du fluide demande cependant quelques options spécifiques. En effet, la matrice de raideur assemblée globale du problème couplé fluide-structure est intrinsèquement singulière (cf. documentation [R4.02.02]) au niveau des degrés de liberté fluides. Il convient donc d’exclure ces degrés de liberté de l’analyse de stabilité, mais aussi de modifier la matrice de raideur (ainsi que la matrice de raideur géométrique quand elle est utilisée). Pour cela, il faut renseigner les mot-clés suivants, sous CRIT_STAB:

• MODI_RIGI = “OUI”,

• DDL_EXCLUS=(“PHI”,”PRES”,”DH”,).

La liste des degrés de liberté exclus doit comporter tous les types de degrés de liberté liés au modèle fluide: dans l’exemple du cas-test FDNV100 on a donc le potentiel PHI, la pression PRES et le déplacement vertical au niveau de la surface libre DH. Si l’on ne fait pas ce traitement, alors l’appel à CRIT_STAB va planter pour cause de matrice singulière et aucune stratégie de décalage ne saurait surmonter cela.

En quasi-statique, ce problème ne se pose pas car la modélisation fluide-structure couplée n’a alors pas de sens.