# Contrib:KeesWouters/prestressmodal

## Modal analysis of an axial loaded cylinder

To start with, this contribution mainly focuses on the use of Salome and Code Aster, not on the results and the mechanical justifications of the code that has been used. So no guarantee that the results will be correct up to the fifth decimal place, which they are not. I do hope though that this information is useful. For me it has been, because I had to think about some commands and look through the documentation and learn from that. In case of mistakes, errors or remarks, please notify me, or better, you are invited to correct or edit them yourself. Enjoy.

As usual, the files that define the geometry, mesh and command files can be found at the end of this contribution.

## Defining the mechanical properties of the cylinder

### Geometry of the cylinder The model is constructed by extrusion cq piping from a circle with diameter 4 mm. The extrusion is defined along a path along the z axis, defined by six points *), total length of the curve is 5 mm. The main commands of the python script are:

define points:
null = geompy.MakeVertex(0.00, 0, 0)
pz1 = geompy.MakeVertex(0.00, 0, 1)
....
pz5 = geompy.MakeVertex(0.00, 0, 5)
define circle: Circle = geompy.MakeCircle(centre, vector, radius)
define path: Curve = geompy.MakePolyline([null, pz1, pz2, pz3, pz4, pz5 ])
define cylinder: Pipe = geompy.MakePipe(Circle, Curve)

*) I used six points to do some experiments with the Pipe command. Of course, here is not necessary to use six points.

For applying boundary conditions and loads and for defining certain mechanical characteristics to the shell elements in the command file of Code Aster, we define on the cylinder the following entities: top and bottom edges and the cylinder wall:

Surfin = geompy.CreateGroup(Pipe, geompy.ShapeType["FACE"])
geompy.UnionIDs(Surfin, [25, 20, 15, 10, 3])
Pipe = geompy.GetMainShape(Surfin)
Ctop = geompy.CreateGroup(Pipe, geompy.ShapeType["EDGE"])
geompy.UnionIDs(Ctop, )
Pipe = geompy.GetMainShape(Ctop)
Cbot = geompy.CreateGroup(Pipe, geompy.ShapeType["EDGE"])
geompy.UnionIDs(Cbot, )
Pipe = geompy.GetMainShape(Cbot)

### The mesh of the cylinder

The mesh is defined in a Python script. It is generated by NetGen 1D_2D. By allowing quadrangles, a regular mesh with only this type of elements occur. Further options are: maxsize is 0.075, optimise, fineness is 3 and use second order, that is needed for coque_3d elements. Note the these quadratic quadrangles (quad8) need to be adapted to quad9 elements to be compatible with coque_3d elements in the command file.

The names of the axial ends and the cylinder itself are denoted by Cbot, Ctop and Surfin.

### The load and boundary conditions of the cylinder

The load and boundary conditions are defined in the Code Aster command file. The boundary conditions are at Cbot and Ctop:

Cbot
restrict displacements and rotations dx=dy=dz=drx=dry=drz=0 -->
DDL_IMPO=(_F(GROUP_MA=('Cbot'),DX=0.0,DY=0.0,DZ=0.0,DRX=0.0,...),...
Ctop
restrict displacements and rotations dx=dy=drx=dry=drz=0 except dz -->
DDL_IMPO=(...._F(GROUP_MA=('Ctop'),DX=0.0,DY=0.0,DRX=0.0,...),...

The load on the construction is an axial line force, total force is 7399.515 N -->

dcyl = 4.00 # diameter of cylinder
pie = math.pi
Lcyl = pie*dcyl # circumference of the cylinder
Fax = 7399.515
Fclamp=AFFE_CHAR_MECA(...,FORCE_ARETE=_F(GROUP_MA='Ctop',FZ=+Fax/Lcyl,),);

### Material properties of the beam

The beam consists of steel.

steel=DEFI_MATERIAU(ELAS=_F(E=2.100e5,NU=0.28,RHO=0.007850),);

The dimensions of the quantities are: E [MPa], NU [-] and RHO [g/mm3].???

## General procedure during the solution process

We will compare the modal behaviour of the unloaded and loaded cylinder. For the unloaded cylinder the standard modal analysis will be performed:

determine stiffness matrix Smech and
determine mass matrix Mmass
carry out a modal analysis (Smech + lambda*Mmass)x = 0

For the pre stressed cylinder a modified modal analysis will be performed. The stiffness matrix is now given by the sum of the element stiffness matrix and geometrical stiffness matrix: Stiff = [Smech + Sgeom]:

determine stiffness matrix Smech
determine geometrical stiffness matrix Sgeom
determine mass matrix Mmass
perform modal analysis ([Smech + Sgeom] + lambda*Mmass)x = 0

Determining the geometrical stiffness matrix is exactly the same as in a buckling analysis. Details for the calculating the geometrical stiffness matrix can be found here: buckling analysis.

## Modal analysis of the unloaded structure

### Commands for modal analysis of the unloaded structure

The element stiffness matrix can be calculated by the following commands

optional: perform static analysis [MECA_STATIQUE] *)
determine stiffness matrix Smech [CALC_MATR_ELEM]
determine mass matrix Mmass [CALC_MATR_ELEM]

*) Maybe used for testing the results, eg the boundary conditions and the load, the static analysis can be carried out. Though it is not really needed here.

Then, for the eigenvector calculation, the stiffness and mass matrices need to be correctly ordered according to the dofs (ddl) of the system, so

determine ordering nddl [NUME_DDL]
determine ordered SAmech [ASSE_MATRICE] and
determine ordered SAmass [ASSE_MATRICE]
perform eigenvalue analysis on SAmech and SAmass [MODE_ITER_SIMULT]
project the results [PROJ_CHAMP] suitable for post processing in Salome (coque_3d ie quad9 --> quad8)
```# for testing only
result=MECA_STATIQUE(MODELE=modelc,CHAM_MATER=matprops,
CARA_ELEM=shellch,EXCIT=_F(CHARGE=clamp,),);
```
```# determine stiffness matrix of the structure and DDL vector
Smech = CALC_MATR_ELEM(OPTION = 'RIGI_MECA',
MODELE = modelc,
CARA_ELEM=shellch,
CHARGE = clamp,
CHAM_MATER = matprops,);
```
```Mmass = CALC_MATR_ELEM(OPTION = 'MASS_MECA',
MODELE = modelc,
CHAM_MATER = matprops,
CARA_ELEM = shellch,
CHARGE = clamp,);
```
```nddl   = NUME_DDL(MATR_RIGI=Smech);
SAmech = ASSE_MATRICE(NUME_DDL=nddl,MATR_ELEM=Smech,);
MAmass = ASSE_MATRICE(NUME_DDL=nddl,MATR_ELEM=Mmass,);
```
```Rmodes=MODE_ITER_SIMULT(MATR_A=SAmech,
MATR_B=MAmass,
TYPE_RESU='DYNAMIQUE',
METHODE='SORENSEN',
CALC_FREQ=_F(OPTION='PLUS_PETITE',NMAX_FREQ=6,),);
```
```SALmodes=PROJ_CHAMP(RESULTAT=Rmodes,
MODELE_1=modelc,
MODELE_2=modinit,);
```
```IMPR_RESU(FORMAT='MED',UNITE=21,
RESU=_F(MAILLAGE=meshinit,
RESULTAT=SALmodes,
NOM_CHAM='DEPL',),);
```

### Numerical values of the eigenvalues of the unloaded structure

The eigenvalues can be obtained from the mess file. For the unloaded case, the six smallest eigenvalues are:

```   NUMERO    FREQUENCE (HZ)    NORME D'ERREUR
1      9.57306E+01       3.17056E-11
2      9.57306E+01       3.14144E-11
3      1.13379E+02       2.33041E-11
4      1.13379E+02       2.30014E-11
5      1.22476E+02       2.36309E-11
6      1.22476E+02       1.96973E-11
NORME D'ERREUR MOYENNE:  0.25459E-10
```

### The first modal shape of the unloaded structure

The picture below denote the first modal shape at 96 Hz at different view points. * ## Modal analysis of the structure under prestress

### Commands for modal analysis of the structure under prestress

The element stiffness matrix can be calculated by the following commands

perform static analysis [MECA_STATIQUE], now needed to determine the stresses in the construction,
determine the stresses [CREA_CHAMP] calculated in the static analysis,
determine stiffness matrix Smech [CALC_MATR_ELEM],
determine geometrical stiffness matrix Smech [CALC_MATR_ELEM], based on the stresses from the static analysis and
determine mass matrix Mmass [CALC_MATR_ELEM]

Then, for the eigenvector calculation, the stiffness and mass matrices need to be correctly sorted according to the dofs (ddl) of the system, so

determine ordering nddl [NUME_DDL],
determine ordered SAmech [ASSE_MATRICE],
determine ordered SAgeom [ASSE_MATRICE] and
determine ordered SAmass [ASSE_MATRICE]
add the SAmech & SAgeom [COMB_MATR_ASSE] (COEF the multiplier for each matrix is 1). The total stiffness matrix is Stotal
perform eigenvalue analysis on Stotal and SAmass [MODE_ITER_SIMULT]
project the results [PROJ_CHAMP] suitable for post processing in Salome (coque_3d --> quad8)

```result=MECA_STATIQUE(MODELE=modelc,
CHAM_MATER=matprops,
CARA_ELEM=shellch,
EXCIT=_F(CHARGE=clamp,),);
```
```Fsigma = CREA_CHAMP(OPERATION='EXTR',
TYPE_CHAM='ELGA_SIEF_R',
RESULTAT=result,
NUME_ORDRE=1,
NOM_CHAM='SIEF_ELGA_DEPL',);
```
```# determine stiffness matrix of the structure and DDL vector
Smech = CALC_MATR_ELEM(OPTION = 'RIGI_MECA',
MODELE = modelc,
CARA_ELEM=shellch,
CHARGE = clamp,
CHAM_MATER = matprops,);
```
```Sgeom = CALC_MATR_ELEM(OPTION='RIGI_GEOM',
MODELE=modelc,
CARA_ELEM=shellch,
SIEF_ELGA=Fsigma,);
```

```Mmass = CALC_MATR_ELEM(OPTION = 'MASS_MECA',
MODELE = modelc,
CHAM_MATER = matprops,
CARA_ELEM = shellch,
CHARGE = clamp,);
```
```nddl   = NUME_DDL(MATR_RIGI=Smech);
SAmech = ASSE_MATRICE(NUME_DDL=nddl,MATR_ELEM=Smech,);
SAgeom = ASSE_MATRICE(NUME_DDL=nddl,MATR_ELEM=Sgeom,);
MAmass = ASSE_MATRICE(NUME_DDL=nddl,MATR_ELEM=Mmass,);
```
```Stotal = COMB_MATR_ASSE(COMB_R=(_F(MATR_ASSE=SAmech,COEF_R=1.0,),
_F(MATR_ASSE=SAgeom,COEF_R=1.0,),),);
```
```Rmodes=MODE_ITER_SIMULT(MATR_A=Stotal,
MATR_B=MAmass,
TYPE_RESU='DYNAMIQUE',
METHODE='SORENSEN',
CALC_FREQ=_F(OPTION='PLUS_PETITE',NMAX_FREQ=6,),);
```
```SALmodes=PROJ_CHAMP(RESULTAT=Rmodes,
MODELE_1=modelc,
MODELE_2=modinit,);
```
```IMPR_RESU(FORMAT='MED',UNITE=21,
RESU=_F(MAILLAGE=meshinit,
RESULTAT=SALmodes,
NOM_CHAM='DEPL',),);
```

### Calculating matrices by MACRO_MATR_ASSE

The sorted stiffness and mass matrices have been determined by the two commands CALC_MATR_ELEM and ASSE_MATRICE. It is also possible to use a single command: MACRO_MATR_ASSE.

```Smech = CALC_MATR_ELEM(OPTION = 'RIGI_MECA',MODELE = modelc,CARA_ELEM = shellch,CHAM_MATER = matprops,CHARGE = clamp,);
Mmass = CALC_MATR_ELEM(OPTION = 'MASS_MECA',MODELE = modelc,CARA_ELEM = shellch,CHAM_MATER = matprops,CHARGE = clamp,);
Sgeom = CALC_MATR_ELEM(OPTION = 'RIGI_GEOM',MODELE = modelc,CARA_ELEM = shellch,SIEF_ELGA = Fsigma,);

nddl   = NUME_DDL(MATR_RIGI=Smech);
SAmech = ASSE_MATRICE(NUME_DDL=nddl,MATR_ELEM=Smech,);
SAgeom = ASSE_MATRICE(NUME_DDL=nddl,MATR_ELEM=Sgeom,);
MAmass = ASSE_MATRICE(NUME_DDL=nddl,MATR_ELEM=Mmass,);
```

Instead you can use the MACRO_MATR_ASSE command:

```MACRO_MATR_ASSE(MODELE=modelc,CHAM_MATER=matprops,CARA_ELEM=shellch,CHARGE=clamp,
NUME_DDL=CO('nddl'),  # define 'nddl' here
MATR_ASSE=(_F(MATRICE=CO('SAmech'),OPTION='RIGI_MECA',),
_F(MATRICE=CO('MAmass'),OPTION='MASS_MECA',),),);
```
```MACRO_MATR_ASSE(MODELE=modelc,
#CHAM_MATER=matprops,
CARA_ELEM=shellch,
CHARGE=clamp,
NUME_DDL=nddl,     # nddl already defined prevously
MATR_ASSE=(_F(MATRICE=CO('SAgeom'),OPTION='RIGI_GEOM',SIEF_ELGA=Fsigma,),),);
```

Remark:

• because the geometrical stiffness matrix does not require and does not accept the CHAM_MATER option, the calculation of the geometrical stiffness matrix has been separated from the stiffness and mass matrices.
• the definition of nddl can only occur once, so in the first MACRO_MATR_ASSE command use NUME_DDL=CO('nddl') to define and calculated nddl. Then in the second command it can be used directly as NUME_DDL=nddl.

### Numerical values of the eigenvalues of structure under prestress

The eigenvalues can be obtained from the mess file. For the unloaded case, the six smallest eigenvalues are:

```   NUMERO    FREQUENCE (HZ)    NORME D'ERREUR
1      1.30828E+02       1.63671E-11
2      1.30828E+02       1.71866E-11
3      1.44320E+02       1.49786E-11
4      1.44320E+02       1.45535E-11
5      1.51135E+02       1.51974E-11
6      1.51135E+02       1.48339E-11
NORME D'ERREUR MOYENNE:  0.15520E-10
```

And compared to the non loaded structure, we have

```       loaded  no load difference
[Hz]	 [Hz]	    [%]
1	130.8	  95.7	    37
2	130.8	  95.7	    37
3	144.3	 113.4	    27
4	144.3	 113.4	    27
5	151.1	 122.5	    23
6	151.1	 122.5	    23
```

### The first modal shape of the loaded structure

The picture below denote the first modal shape at 131 Hz from different view points. * ### Write the total mass of the construction

Using the MASS_INER command, the total mass of the construction can be written to a file (here on UNITe 26, defined in ASTK as npcyl.txt, download at the end)

```masshell=POST_ELEM(MASS_INER=_F(GROUP_MA='Surfin',ORIG_INER=(0.0, 0.0, 0.0),),
MODELE=modelc,
CHAM_MATER=matprops,
CARA_ELEM=shellch,);
```

IMPR_TABLE(TABLE=masshell,UNITE=26,NOM_PARA='MASSE',);

The content of the results file is:

```.....
#ASTER  9.04.00 CONCEPT masshell CALCULE LE 03/09/2009 A 17:53:48 DE TYPE
#TABLE_SDASTER
MASSE
4.93230E-02
```

This agrees with the mass calculated by:

mass = rho * Volume = rho * pi*Dcyl*Hcyl*th = 0.00785*pi*4*5*0.1 = 0.0493 kg.

## Buckling analysis

The prestress in the construction is determined such that the buckling load is -1, eg. the both the element stiffness and the geometrical stiffness have roughly the same order of magnitude. The command is as follows:

```# multiply results by -1 or reverse force, since we solve
# A*x = lamdba*B*x iso (A+mu*B)x=0
Rbuckle = MODE_ITER_SIMULT(INFO=1,
MATR_A  = SAmech,
MATR_B  = SAgeom,
METHODE = 'SORENSEN',
PREC_SOREN=2e-16,
NMAX_ITER_SOREN=250,
TYPE_RESU='MODE_FLAMB',
CALC_FREQ=_F(OPTION='PLUS_PETITE',NMAX_FREQ=4,),);
```
```IMPR_RESU(FORMAT='MED',
UNITE=21,
RESU=_F(MAILLAGE=meshinit,
RESULTAT=Rbuckle,
NOM_CHAM='DEPL',),);
```

and the corresponding four lowest buckling loads are, see mess file:

```   NUMERO    CHARGE CRITIQUE    NORME D'ERREUR
1      9.86304E-01       6.80828E-12
2      9.86305E-01       6.84256E-12
3      9.99999E-01       5.58740E-12
4      1.00000E+00       5.56880E-12
NORME D'ERREUR MOYENNE:  0.62018E-11
```

The load at number 4 corresponds to a total axial force of Fax = 7399.515 N. The corresponding bluckling mode is given below.
Note:

• that we applied an axial pulling force, and again, the lowest buckling load at 4 need to be multiplied by -1 because of the implementation of the eigenvalue analysis in Code Aster.
• that the highest value of charge-critique corresponds to the lowest critical load.

### The buckling mode of the loaded structure

The picture below shows the lowest buckling mode for an axial load of 7399.515 N at different view points. * The displacement according to Code Aster at the outer edge Ctop of the cylinder is 0.1387 mm for the given axial load. The analytical displacement is

• dL = F*L/(E*A) = 0.1402:
• F = 7399.515 N
• L = 5.00 mm (height of the cylinder)
• E = 2.1e5 MPa
• A = pi*D*th = 1.2566 mm2, D is diameter and th is thickness of the cylinder.

## ASTK files

In the picture below the various input and output files of the calculation are given. Media:kw_cylinder.zip

This zip contains the following files

• python geometry and mesh files (load in Salome by File --> Load script (cntrl T)) and right select refresh (F5) in the object browser). Export the med file in the mesh module under cylinder.med for further processing by Code-Aster, controlled by ASTK.
• command files for Code-Aster:
•     - cyl_modal.comm, modal analysis of the unloaded structure
•     - cyl_loaded.comm, modal analysis of the structure under pre stress, buckling also included
• astk file for file control by ASTK
• npcyl.txt, text file with nodal displacements and force on plane top and total mass of the structure, not very spectacular.

The result file buck2res.med can be viewed in the post processor module of Salome by File --> Import --> buck2res.med or cntrl I in the object browser.

Note on Code Aster 10:
Code Aster 10+ is stricter then CA9- with respect to the MODI_MAILLAGE keyword in CREA_MAILLAGE. In case no TRIA6 or QUAD8 elements are available in the mesh, change near line 30 in both command files:

```meshmod=CREA_MAILLAGE(MAILLAGE=meshinit,
_F(TOUT='OUI',OPTION='TRIA6_7',),),);
```

```meshmod=CREA_MAILLAGE(MAILLAGE=meshinit,
```meshmod=CREA_MAILLAGE(MAILLAGE=meshinit,