3. Quasi-static push-over approach

3.1. Imposed loads

The imposed load is of regulatory origin (EC-8), almost static [bib1]. The effects of an earthquake on the structure are represented by a variable pressure field imposed on the inner face of the shells. The value at each point is a function of the current coordinate and increases linearly with time. This monotonic evolution in time is characteristic of so-called push-over methods, in the EC-8 sense, whose objective is to simulate the response to seismic stress by means of a quasistatic calculation, and therefore of a transitory dynamic physical nature. Dynamic effects, such as inertia and the forces generated by the sloshing fluid, are replaced by this imposed pressure distribution. The problem to be treated only involves modeling the structure, without modeling the fluid domain.

Push-over methods have obviously been built and justified by making strong assumptions about the linearization of the problem (small displacements, elastic behavior, Euler buckling, etc.). Below are the definitions of the imposed pressure fields [bib2]:

# ——————————————————

# DEFINITION DES CHAMPS FROM 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”, “Y”, “Z”), VALE =” PIFZ (Z) * COSTE (X, Y) “)

PIR = FORMULE (NOM_PARA =( “X”, “Y”, “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)) “)”)

It is then possible to define the corresponding mechanical loads (gravity and subsequent pressures):

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

These pressure fields are multiplied by a linear increasing function of time:

FONCMUL = DEFI_FONCTION (NOM_PARA = “INST”, = (0., 0., 3., 3.),)

The connection with the ground is considered to be complete here (group of nodes” BASE “). We also define the symmetry conditions (group of nodes” SYMETRIE “) since we only mesh half a reservoir:

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. Quasi-static monotonic nonlinear resolution method

We want to solve a nonlinear, almost static evolution problem. We will therefore use the operator STAT_NON_LINE ([U4.51.03], [R5.03.01]). The imposed load will be built with pressure PRESPS1 for example:

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

)

3.3. Notes on calculations and post-treatments

In the example presented above, a Newton algorithm is used: we solve with the global tangent operator (stiffness) updated at each iteration. If the problem is well posed (sufficiently « regular »), it is known that this type of algorithm offers the best convergence. So, in our case, as long as the global tangent operator is « well » defined (far from being singular), the calculation will proceed well with rapid convergence. When the imposed load level will approach the ultimate load, the structure will then become unstable in the sense of buckling [R7.05.01], which results in the tendency of the global tangent operator to become singular. The loss of stability per limit point is in fact the loss of uniqueness of the solution, i.e. therefore the singularity of the resolution operator. In the vicinity of the limit point, Newton’s algorithm will converge less well, which is why it is necessary to impose smaller time increments and an increase in the number of iterations on the residue in equilibrium.

In general, the closer you get to the ultimate load, the smaller the time step should be. Despite this, the risks of stopping the calculation on non-convergence are significant, hence the obligation to carry out the calculation following several successive pursuits.

However, it is possible to improve this convergence by changing the algorithm during the calculation and to switch to a quasi-Newton. To do this, all you have to do is solve on the tangent operator that you only update at each step (between two iterations you keep the same), and if this is insufficient, you can then solve with the elastic operator, which cannot become singular. This choice reinforces the robustness of the algorithm in terms of convergence, but it increases, sometimes considerably, the number of iterations (and/or steps) required to obtain the solution.

For our type of study, we can distinguish two types of quantities of interest for post-treatment. On the one hand, a quantity capable of representing the buckling of the structure and therefore of revealing the ultimate load. To do this, we can draw the pressure level (multiplier coefficient) as a function of the movement from a point located at the top of the tank (at the end of the generator most subjected to compression). On the other hand, a more local indicator of the onset of plasticity: the iso-value of the cumulative plastic deformation, at each calculation moment.

These two post-treatments do not present any particular difficulty in*Code_Aster*.

The finite element models implemented consisted of between 55000 and 110000 degrees of freedom. These models were based on a simplified representation of the anchors: the geometry of the reinforcements and gussets is finely meshed, but the bolting is not present and it is replaced by an embedment condition on all nodes at altitude \(0m\). A more realistic modeling of the anchorages, with detachment, would therefore have resulted in a significantly larger overall problem size.

The calculation methodology presented here, whose objective is to study the quasi-static buckling response [R7.05.01] of a reservoir, is very close to the framework of the buckling calculation instructions [U2.08.04]. This documentation presents linear stability analyses (in the Euler sense) and the development of a nonlinear push-over calculation to obtain the ultimate load.

3.4. Fine modeling of anchorages: lifting

The type of tanks studied here are bolted to the ground [bib5]. These bolts (or tie rods) pass through reinforced flanges welded to the base of the shell). For more speed of calculation, we then present the method implemented on a simplified mesh, but with realistic anchorages: the tank is fixed to the ground by 18 tie rods.

In Code_Aster, several contact models are available.

The modeling of the geometry of the contact zone can be surface (3D problem), linear (2D problem) or composed of discrete elements (keywords DIS_CHOC in STAT_NON_LINE and DIS_CONTACT for the material law).

The contact itself can be treated, either nodally by penalization or by the Lagrangian method, with active or non-active constraints, or continuously by the augmented Lagrangian method.

The simplest method, in this study configuration, is that of penalization, which in this case also has the advantage of being able to take into account the role of the seal without having to mesh it separately. For a semi-static calculation, penalization is not numerically problematic, as in dynamics where this technique can generate high-frequency disturbances (linked to the stiffness of the penalty). In addition, for structural elements, it is possible to adjust the penalization stiffness in order to best approximate the contact stiffness of the part, which is in reality massive.

However, the other methods are more rigorous because they do not lead to interpenetration.

To study the detachment of the tarpaulin, we will therefore use a penalization method, which presents the best modeling quality/calculation cost ratio in our specific case.

Given the option taken to force all the nodes to be able to move only vertically, we can take advantage of this total pairing to content ourselves with using discrete elements for contact. Indeed, it is not necessary to do re-pairing, and therefore, the general master-slave surface models of the contact zones are unnecessarily cumbersome. We will therefore place a carpet of discrete contact elements under the low fret. We will then have a discrete element (DIS_T on mesh SEG2) under each node of altitude 0.

If the calculated solution proves to be too far from the experimental reference solution, it will be necessary to quantify the influence of the vertical displacement condition only. It will suffice to restart calculations without this kinematic hypothesis, with, if the horizontal displacements are large, with a contact management method with re-matching of the nodes of the surfaces concerned.

Remarks

Rigorously, all the nodes of the base shell should be provided with a discrete contact element. However, if we use elements of type COQUE_3Dpour the frets, we cannot place a contact element at the nodes in the center of the elements. In fact, this knot has the particularity of only carrying rotation-type degrees of freedom. The contact condition, which relates to the normal displacement of the face, cannot therefore be expressed in these nodes. The contact elements, as well as the kinematic condition of vertical displacement of the base nodes, must therefore only be carried by the vertices or midpoints of the edges.

In Code_Aster, contact modeling by discrete elements makes it possible to directly define the preload of these elements. It is thus possible to represent the prestress of the tie rods, without having to pre-deform these elements, as would have been the case if another contact model based on the surfaces facing each other had been chosen.

Below, we present the position of the discrete penalized contact elements (in light gray) and of the 9 springs modeling the tie rods (in black). The reinforcement at the base of the shell is indicated in dotted lines:

Figure -a: Arrangement of the discrete elements for the connection with the slab

These discrete elements for anchors are grouped into 2 mesh groups:

• RESSC: contact elements in light grey (one SEG2 element per node of the lower bridle, i.e. 57 SEG2),

• RESS: black tie rods (9 springs per SEG2 elements).

An additional material (MATRES) is therefore introduced which corresponds to the contact elements and which makes it possible to take into account prestress by tightening the nuts on the tie rods (the distance associated with the height between flanges is also defined with DIST_1 and DIST_2):

Bolt = 2.E11;

Devot = 0.026;

Bolt = 3.14159 Bolt Bolt/4;

spring = Bolt * Bolt;

kpenal = spring * 100.

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

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

We complete the material assignment with these discrete elements on the group of elements 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,),),)

The elementary characteristics of the new discrete elements are defined as follows:

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 = (spring, 0., 0.,),),),)

The springs that model the tie rods (GROUP_MARESS) do have the equivalent stiffness kspring.

It remains to modify the moving boundary conditions for the anchorages: only vertical movements of the nodes on the underside of the lower flange (knots initially in contact with the ground) are authorized.

The operator STAT_NON_LINE ([U4.51.03], [R5.03.01]) also sees its arguments impacted by the introduction of the uprising in base (keyword DIS_CHOC for the mesh group 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”),),),

…

);

The rest of the command file is unchanged.

3.5. Use of the nonlinear stability criterion

We can also use a stability criterion based on the tangent matrix: we have instability if the tangent stiffness matrix becomes singular, or if at least one of its eigenvalues is cancelled out. We then solve the problem with the following eigenvalues, written in large displacements (written in Lagrangian with the Green-Lagrange deformation tensor) [R7.05.01]:

\(({K}_{T}+\lambda {I}_{d})x\mathrm{=}0\mathrm{\iff }{K}_{T}x\mathrm{=}\lambda {I}_{d}x\) with:

\({K}_{T}\mathrm{=}K+{K}^{L}(\underline{u})+{K}^{Q}(\underline{u})+K(\pi )\)	Tangent stiffness matrix
\({K}^{L}(\underline{u})\)	Linear part in \(\underline{u}\) of the matrix \({K}_{T}\)
\({K}^{Q}(\underline{u})\)	Quadratic part in \(\underline{u}\) of the matrix \({K}_{T}\)
\(K(\pi )\)	Geometric stiffness matrix
\(\pi\)	Piola-Kirchhoff tensor \(\mathrm{II}\)
\({I}_{d}\)	Identity matrix
\(\lambda\)	Eigenvalue

The [R7.05.01] documentation describes these stability analyses in more detail. The [U2.08.04] documentation shows how to use it.

Remarks

When the trips are small, we just have \(K(\pi )=K(\sigma )={K}_{G}\) and the matrixs \({K}^{L}\) and \({K}^{Q}\) can be overlooked.

Linear analysis in the sense of Euler in Code_Aster does not make it possible to take into account the follower aspect of the pressure field: it is mandatory to pass the approach presented here and based on a non-linear resolution.

If you want to use the stability criterion, simply add the following line among the arguments of STAT_NON_LINE ([U4.51.03], [R5.03.01]):

CRIT_STAB =_F (NMAX_CHAR_CRIT =1,),

We then look for the first eigenvalue of the global tangent operator of our system.

If during the calculation we observe that this eigenvalue decreases, or even changes sign, this means that we have approached the first critical load and that we have then even exceeded it.

The number of eigenvalues to be determined can be imposed by the NMAX_CHAR_CRIT keyword (3 by default).

It is also possible, using the CHAR_CRIT command, to choose the band in which to look for these eigenvalues (from —10 to 10 by default).

Note

The indication of a frequency band is useful especially for calculations in small disturbances where a Sturm test is performed for the frequency band provided. We can thus save time by only calculating the eigenvalues if there are any in the indicated band. The Sturm test is not done in large deformations and the eigenvalues are calculated at each time step.

The buckling mode as well as the determined eigenvalues (critical load coefficient) can be retrieved using the IMPR_RESU command:

IMPR_RESU (MODELE = MODELE, FORMAT = “RESULTAT”, RESU =( _F (RESULTAT = RESU,

NOM_PARA =” CHAR_CRIT “,” MODE_FLAMB “,),)

Test case SSLL105 [V3.01.105] proposes an example of using this stability criterion for a linear case and test case SSNL126 [V6.02.126] for a non-linear case (elastoplastic beam).

3.6. Load control

In order to facilitate the convergence of the incremental calculation when one is close to the ultimate load level, or in order to be able to exceed this critical point, it may be wise to no longer use imposed loading in order to favor piloting while moving or piloting by arc length (its use for post-criticism is briefly mentioned in the documentation [U2.08.04]). Control cannot be used with the [U4.51.03] contact.