4. Coupled transitional approach#
4.1. Reference problem#
Here we are going beyond the regulatory framework and we are going to exploit all the modeling possibilities offered by*Code_Aster*. The tank model itself remains unchanged (elastoplastic volume shells). On the other hand, we will represent the fluid domain by a massive mesh. In addition, the resolution will be carried out in a transitory dynamic with the operator DYNA_NON_LINE ([U4.53.01], [R5.05.05]), the external stress being of the seismic type.
The fluid domain is modelled in linear acoustics (barotropic, compressible, non-viscous and with free surface). The fluid-structure problem is solved in Code_Aster by a symmetric formulation \((u,p,\varphi )\) ([R4.02.02], [bib7]), in updated Lagrangian script. The load is of the accelerogram type imposed on the base of a tarpaulin.
The discretized problem then looks like this [bib6]:
|
Figure -a: Full mesh with fluid domain and structure |
4.2. Coupled fluid-structure modeling in Code_Aster#
In order to be able to carry out a coupled calculation by formulation \((u,p,\varphi )\) [bib7] in Code_Aster, with free surfaces and sliding at the interface (non-adherence between the non-viscous fluid and the inner wall of the tank), it is necessary to respect a certain construction of the mesh and the corresponding models.
We must therefore define the fluid mesh, the structure mesh and the fluid-structure interface.
In the mesh, two distinct meshes are therefore generated (different meshes but common nodes) for the fluid domain and the structure. For practical reasons, it may be easier to generate the distinct meshes separately, i.e. with split nodes, and then to eliminate these double nodes to have only single nodes.
Then, in Code_Aster, we will generate the support mesh for the interface as follows, from the mesh groups SURF0, SURF1et SURF2qui are the shells of the reservoir:
MAILLA01 = CREA_MAILLAGE (MAILLAGE = MAILLA1,
CREA_GROUP_MA =_F (NOM =” IFLUSTRU “,
GROUP_MA =( “SURF0”, “SURF1”, “SURF2”,),
PREF_MAILLE =”I”,),);
We therefore create a new group of elements IFLUSTRU. The nodes in this mesh are confused with those in the structure, but the meshes are different.
Remarks
The structure is meshed into elements COQUE_3D, whose geometric support is a 9-node quadrilateral. The middle knot has the particularity of carrying only degrees of freedom of rotation.
If one wants to couple such an element with a massive fluid element, one cannot therefore write a kinematic coupling condition for this middle node since it does not have a degree of translational freedom, unlike the corresponding node which comes from the fluid domain and which, in turn, only carries degrees of translational freedom.
To get around this problem, we will only write fluid-structure coupling on the nodes of the shell mesh including the degrees of freedom of translation: the vertex nodes and the middle nodes of the edges of the elements.
In Code_Aster, you must therefore have the mesh structure comprising only quadrilateral elements with 8 nodes (these are the vertex nodes and the middle nodes of the edges), on which the interface is defined. The mesh structures for COQUE_3Détant defined only after, from this mesh, by adding the middle nodes.
The fluid mesh, to be consistent with the interface, will be composed of 20-node parallelepipeds.
In addition, even if a structural mesh was used with shell elements whose geometric support would be a 9-node quadrilateral, and where all the nodes, even the central one, carried translational degrees of freedom, fluid-structure coupling could be a problem. Indeed, the conformal massive fluid mesh should be achieved with massive parallelepipedic elements at 27 knots. However, some manufacturers do not offer this type of complete elements which are fairly rarely used in structural calculation, unlike the classical quadratic elements that are 20-node parallelepipeds.
It is imperative that the normal outside the fluid domain is always oriented in the same direction. It is strongly recommended to keep the convention of orienting the structure towards the fluid for all fluid-structure interface models, in particular FLUI_STRU.
Once the definition of the interface, which is therefore meshed in quadrilateral elements with 8 nodes (which are equivalent to the faces of the massive elements used for the fluid domain: parallelepiped elements with 20 nodes), we can therefore modify the mesh structure (group of elements RESERVOI) in elements with 9 nodes (which are equivalent to the faces of the massive elements used for the fluid domain: parallelepiped elements with 20 nodes), we can therefore modify the mesh structure (group of elements) into elements with 9 nodes, geometric support for COQUE_3D:
MAILLA2 = CREA_MAILLAGE (MAILLAGE = MAILLA01,
MODI_MAILLE =_F (GROUP_MA =” RESERVOIR “,
OPTION =” QUAD8_9 “,
PREF_NOEUD =” NSQ “,
PREF_NUME =1,),);
We can then define the models necessary for the coupled calculation:
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 “,),),);
The group of elements SURFLIBR (edge cells of the fluid domain, located on its upper face) carries a free surface model [R4.02.04]: 2D_ FLUID_PESA.
The cell groups FLUID0 (massive fluid domain), FOND (edge cells of FLUID0définissant the bottom), and PLANCENT (edge cells of FLUID0définissant the plane of symmetry) define the total fluid domain: 3D_ FLUIDE.
The fluid-structure interface (FLUI_STRU) is carried by the group of elements IFLUSTRU.
Finally, the reservoir mesh (RESERVOI) is the support for the solid shell structure model (COQUE_3D).
The fluid compressible material is also defined:
EAU = DEFI_MATERIAU (FLUIDE =_F (RHO =1000.0, CELE_R =1500.0,),);
that will be assigned to the fluid domain (massive model and its edges) as well as to the interface IFLUSTRUet to the free surface SURFLIBR:
CHMAT = AFFE_MATERIAU (MAILLAGE = MAILLA2,
AFFE =( _F (GROUP_MA =( “FLUID0”, “FOND”, “”, “”, “PLANCENT”, “IFLUSTRU”, “SURFLIBR”,),
MATER = EAU,),
…
);
4.3. Boundary conditions#
The kinematic boundary conditions relate to embedding at the base of the tank (SYMETRI2), to the generators in the vertical plane of symmetry (SYMETRIE) and to non-penetration (zero normal speed) at the bottom of the fluid domain (FOND), as well as in its vertical plane of symmetry (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,),);
The effects of gravity are taken into account:
PESA = AFFE_CHAR_MECA (MODELE = MODELE,
PESANTEUR =_F (GRAVITE =9.81,
DIRECTION =( 0.0,0.0, -1.0),),);
However, this command is insufficient because in modeling \((u,p,\varphi )\) [bib7] the effects of gravity in the fluid domain cannot be taken into account. If nothing more was done, gravity would therefore only really be imposed on the structure.
To approximate the effects of gravity of the fluid on the wall, we will impose an equivalent hydrostatic pressure (but which cannot take into account the variation in the height of the fluid during the calculation when the sloshing will start):
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 “,),);
The earthquake is imposed as being an accelerogram (function GASDM_X1) imposed on the base in the direction \(X\). It is therefore a conventional mono-support type of solicitation:
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,);
The mass matrix used (MATMAS) is the mass matrix of the total coupled system.
4.4. Initial conditions#
The initial state of the transitory calculation must correspond to the balance of the total system when it is not subject to the earthquake. This state of static equilibrium therefore corresponds to gravity loading and to hydrostatic effects.
If we started the dynamic calculation with an initial state that did not respect this balance, then this would generate oscillations in the transitory solution since it would not initially be at equilibrium (the seismic level however then being zero). We can mitigate these oscillations by adding « large » structural damping and waiting for the solution to stabilize before imposing the earthquake, but this technique is not very elegant…
To calculate this initial statically balanced state, it is therefore possible to solve a static problem (which is also assumed to be linear) of equilibrium under the action of gravity and hydrostatic forces.
To do this, the global matrices \(K\) and \(M\) are calculated and assembled beforehand with hydrostatic loading and gravity:
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”),),);
Since the assembled stiffness matrix is singular because of the fluid domain (the formulation (u, p,) makes this matrix singular at zero frequency [bib8]), the problem is slightly modified by considering the stiffness matrix \({K}_{\mathrm{cor}}=K+\varepsilon M\) \(\varepsilon \ll 1\) which is no longer singular (it is called RIGICOMB).
For example, you can take \(\varepsilon =\mathrm{0,001}\) as below:
RIGICOMB = COMB_MATR_ASSE (COMB_R =( _F (MATR_ASSE = RIGSTA, COEF_R = 1. ),
_F (MATR_ASSE = MASSTA, COEF_R = -0.001),);
The loading vector F_0 (second member) is then assembled:
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);
We can then solve the linear static problem \({K}_{\mathrm{cor}}\underline{U}=\underline{{F}_{0}}\), using, for example, a factorization of the \({\mathrm{LDL}}^{T}\) type:
RIGICOMB = FACTORISER (reuse= RIGICOMB, MATR_ASSE = RIGICOMB,
STOP_SINGULIER = “NON”);
DEP0 = RESOUDRE (MATR_FACT = RIGICOMB, CHAM_NO =F_0);
The solution field DEP0 will therefore be the initial state of the following transient dynamic calculation.
4.5. Transitional resolution#
We use the DYNA_NON_LINE operator ([U4.53.01], [R5.05.05]) as follows:
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),);
The resolution is done in updated Lagrangian (option DEFORMATION =” PETIT_REAC “) because the fluid domain is in small disturbances on each step. It is therefore necessary to check that the time step is small enough for this hypothesis to be verified.
A modified mean acceleration time integration scheme (SCHEMA =” HHT “, ALPHA =-0.1) with numerical damping is used in order to stabilize the solution and to facilitate convergence.
Test case FDNV100 [V8.03.100] presents the calculation of a rectangular tank full of water with a flexible wall. The modeling used is very similar to that used here for large reservoirs.
4.6. Use of the nonlinear stability criterion#
It is also possible to use a stability criterion just as in quasistatics. However, the presence of the fluid requires some specific options. Indeed, the overall assembled stiffness matrix of the coupled fluid-structure problem is intrinsically singular (cf. documentation [R4.02.02]) in terms of fluid degrees of freedom. It is therefore appropriate to exclude these degrees of freedom from the stability analysis, but also to modify the stiffness matrix (as well as the geometric stiffness matrix when it is used). To do this, you must enter the following keywords, under CRIT_STAB:
MODI_RIGI = “OUI”,
DDL_EXCLUS =( “PHI”, “PRES”, “DH”,).
The list of excluded degrees of freedom must include all the types of degrees of freedom linked to the fluid model: in the example of test case FDNV100, we therefore have the potential PHI, the pressure PRES and the vertical displacement at the level of the free surface DH. If we don’t do this treatment, then the call to CRIT_STAB will crash because of a singular matrix and no lag strategy can overcome this.
In quasistatics, this problem does not arise because coupled fluid-structure modeling then makes no sense.