# Contrib:KeesWouters/shell/platedynamics

## Using coque_3d shell elements

A lot of changes between Code Aster versions 9.4, 10, 11.1 and 11.2&3 have taken place. The original contribution has been written for CA9.4. So please note that in some cases names of files and results may have changed. If you find discrepancies please let me know, or even better, apply the changes right here.

### Dynamics of a rectangular plate

This is just a simple example of the use of shells (coque_3d) for dynamical calculations. It is mainly based on the work of Paul Carrico but here the pre- and post processing is carried out in Salome for creating the geometry and the mesh and displaying the results.

The construction is a rectangle of 1.000 * 1.000 m2 and thickness 0.001 m. We only use the 2D information in Salome for meshing, the thickness will be defined in Code-Aster by the shell elements. All four edges of the rectangle are hinging to the fixed world. The material is steel, its mechanical properties are defined by the following parameters: Young's modulus E = 2.1e11 Pa, poisson ratio NU=0.30 [-] and density RHO=7850 kg/m3.

Note that these are ISO values for all quantities. If dimensions are [m], as they are in the input Python script, then this is no problem. If the input dimension is [mm], the following combinations yields frequencies in [Hz]: Young's modulus E, pressure p [MPa] and density rho [ton/mm3]. See also the very nice table of Paul Carrico.

Versions of the software that has been used april 2009-[january 2013]:

• Salome 4.1.4 (32 bit) - [Salome6.6.0 january 2013]
• Code Aster 9.4 (32 bit) - [Code Aster 11.3 january 2013]
• ASTK 1.7.1 - [ASTK 1.10.09]
• OS: Ubuntu 9.04 Jaunty (64 bit) - [Mint 14 Nadia (Ubuntu 12.10 Precise)]

### Update Code Aster version 10.X/11.1/11.2&3

Changes between CA-9, CA-10, CA-11 are discussed at the end of this contribution.
The result file is now written to the standard file unit 80, type rmed.

Note:
In version 11.3 (and 11.2) the commands CALC_ELEM and CALC_NO have been replaced by CREA_CHAMP, CALC_CHAMP and/or POST_CHAMP to extract the field quantities.

December 2012
The command file has been updated to Code Aster version 11.1.XX and slight reformulation of meshes and models.

## The geometry of the plate

We start by defining the four corner nodes, in the Geometry module of Salome by clicking New Entity --> Basic --> Point (or click the point button in the button bar):

• p00 --> [0.0, 0.0, 0.0]
• p10 --> [1.0, 0.0, 0.0]
• p11 --> [1.0, 1.0, 0.0]
• p01 --> [0.0, 1.0, 0.0]

Defining the four connecting lines by clicking New Entity --> Basic --> Line:

• Ly0 --> p00 and p10 (enter name Ly0 and select the two points)
• Lx1 --> p10 and p11
• Ly1 --> p11 and p01
• Lx0 --> p01 and p00

and create the face defined by the four lines by clicking New Entity --> Build --> Face

• Face1 --> Ly0, Lx1, Ly1, Lx0 (shift select the four lines and click Apply).

For Code Aster we have to define groups for boundary conditions and the shell elements coupled to Face1. For the boundary conditions we choose all four lines. We select the lines by clicking New Entity --> Group --> Create:

• keep default choices except:
• shape type: select line
• main shape: Face1
• name: Linex0, select line Lx0, add, apply and repeat
• name: Liney0, select line Ly0, add, apply
• name: Linex1, select line Lx1, add, apply
• name: Liney1, select line Ly1, add, apply and close.

For the shell elements we need the plane Face1. We select the Face1 by clicking New Entity --> Group --> Create:

• keep default choices except:
• shape type: select plane
• main shape: Face1,
• name: shell, select all, and apply and close.

So far for the definition of the geometry. The geometry is also available by the following python script Media:Geom_plate1by1.zip. The script can be run by File --> Load script --> <<path>>/geometry_plate.py in the geometry module of Salome after dearchiving the file. Maybe you need to right click refresh in the object browser window to view the geometry.

## The meshing of the plate

Switch to the meshing module in Salome. We start with a quadrangle mesh, select Mesh --> Create Mesh. Apply the following settings:

• name: shell
• geometry: Face1
• 1D algorithm: wire discretisation1D: hypothesis automatic length (set to 1, or very fine) 1*)
• apply and close

1*) Set to 0 (very course) for debugging a small mesh.

Now the shell mesh appears in the object browser, right select and hit Compute. A mesh with about 100 edges and 700 quadrangles appear.

Now we have to add the groups already defined in the geometry to the mesh. Using the same strategy as before, add the plane shell and the four edges Linex0, Linex1, Liney0 and Liney1 to the mesh by Mesh --> Create Group.

### Export the mesh to a MED2.3 file

Save the study and export the mesh with filename shell33.med to a MED2.3 file by right select on the name of the mesh in the object browser of Salome.

The definitions of the geometry and mesh are available in the Python script file Media:Mesh_plate1by1.zip. You may need to right click and select refresh to update all the entities in the object browser

## ASTK

Now the mesh is available we can nearly start to use Code Aster. We use ASTK to set up the required files, as defined by the *.astk file. In the resulting *.export file we have the following settings:

• ...
• F mail kees@localhost:/home/kees/shell/shell33.med D 20 (input mesh file)
• F comm kees@localhost:/home/kees/shell/shell33dyn.comm D 1 (input command file)
• F resu kees@localhost:/home/kees/shell/shell33.resu.med R 8 (result file)
• F mess kees@localhost:/home/kees/shell/shell33.mess R 6 (messages file)
• F erre kees@localhost:/home/kees/shell/shell33.erre R 9 (error file)
• R base kees@localhost:/home/kees/shell/shell33.base R 0 (base directory for results)
• F rmed kees@localhost:/home/kees/shell/res33.med R 80 (result file, used in Post Pro of Salome)
• ...

The ASTK window: In this case the result file is a MED file for further processing (post pro) by Salome. I normally copy an old *.astk file into the current working directory and adapt the names of the above mentioned files for the current session.

The settings on the right hand side of ASTK, that do the job for me, are:

• total memory: 2000 Mb (still using 32 bit)
• including Aster: 2000 Mb
• time: 100 (more than enough for this small task)
• interactive (if I choose interactive follow-up the calculation halts)
• nodebug mode
• select Run to start the calculation by Code Aster, but first we need to discuss the command file.

## The command file for Code Aster - overview

The COQUE_3D ekements can be used in combination with quadratic triangular or quadrangular elements. Triangular elements (tria6/7) have 7 nodes, quadrangular elements (quad8/9) have 9 nodes. Every nodes except the centre nodes have 6 dofs: 3 displacements and 3 rotations. The centre nodes only carry rotations, no displacements. Because of this odity a few additional commands are required to handle this.
In the command file the process flow for the calculation is given. In this case the main flow is:

• read the mesh generated by Salome (original mesh initMesh, tria6, quad8 mesh)
• convert the mesh to be suitable for the shell (coque) elements (cocMesh, tria7, quad9 mesh)
• define the model on the modified mesh (cocModel)
• define the characteristics for the shell elements, including thickness of the shells
• set the material properties of the plate, in this case steel
• apply the boundary conditions (and forces) to the construction
• main calculation
• extract the results (cocMesh)
• convert the result to the original mesh to be read by Salome (initMesh)

C'est ca. Now a few remarks.

• The mesh generated by Salome can be either linear or quadratic, depending on your choice in the meshing module of Salome. Since we use quad elements, there are 4 nodes per element (linear mesh, four corner nodes, quad4) or 8 nodes (quadratic mesh, four corner nodes and four midside nodes, quad8). The coque_3d elements however has 9 nodes (four corner, four midside and one centre node) so we need to convert the quad elements into suitable coque_3d elements. Later on, for the post processing, the results have to be converted back to Salome quad elements.
• The corner and midside nodes of the coque_3d element have 6 dofs (dx, dy, dz, drx, dry and drz) whereas the centre node has only 3 dofs (drx, dry and drz).
• the first and last statements are carried out on the initial (Salome mesh). All the other statements are on the modified, say Code Aster suitable mesh.
• For the conversion to coque_3d elements you need quad8 elements or the conversion will fail. So there are a few option to reach this:
• Salome: set quadratic mesh directly
• Salome: convert to quadratic mesh in Salome after linear mesh has been created
• Code Aster: import linear mesh from Salome and convert to quadratic mesh in the command file
• The thickness of the elements is defined in the command file (epais=thickness)
• In this calculation the results are stored on unite 80, corresponding to the LU column in the res file of ASTK.

If somehow you forget to apply the quadratic elements you get the following error in Code Aster's mess file, directly after the CREA_MAILLAGE command:

```meshmod=CREA_MAILLAGE(MAILLAGE=meshinit,.....);
....
<EXCEPTION> JEVEUX_26>
....
Objet JEVEUX inexistant dans les bases ouvertes .....
```

### The command file - detailed commands - CA version 11.1

```DEBUT();
```
```# thickness of shells
th = 0.001
```
```# reading med mesh (code aster, salome, default)
initMesh=LIRE_MAILLAGE(FORMAT='MED',INFO=1,);
```
```#define shell (coque_3D) model;
# we have to add centre nodes with CREA_MAILLAGE
cocMesh=CREA_MAILLAGE(MAILLAGE=initMesh,
_F(TOUT='OUI',OPTION='TRIA6_7',),),); # doesnot work in ASTER10.X, warning in Aster11.X
```
```intModel=AFFE_MODELE(MAILLAGE=initMesh,
AFFE=_F(TOUT='OUI',
PHENOMENE='MECANIQUE',
MODELISATION='3D',),);

cocModel=AFFE_MODELE(MAILLAGE=cocMesh,
AFFE=_F(TOUT='OUI',
PHENOMENE='MECANIQUE',
MODELISATION='COQUE_3D',),);

shellch=AFFE_CARA_ELEM(MODELE=cocModel,
COQUE=_F(GROUP_MA='shell',
EPAIS=th,
##ANGL_REP=(0.1,0.3,),
VECTEUR=(1.0,0.0,0.0),  #defines local x-axis equal to (1,0,0)
COQUE_NCOU=1,),);
```
```steel=DEFI_MATERIAU(ELAS=_F(E=2.100e11,NU=0.30,RHO=7850.0,),);
```
```material=AFFE_MATERIAU(MAILLAGE=cocMesh,AFFE=_F(TOUT='OUI',MATER=steel,),);
```
```#define boundary conditions - hinging or clamped edges
false = (0==1)
true = ~false
clamped=false
if clamped:
DDL_IMPO=(_F(GROUP_MA=('Linex0','Linex1'),
DX=0.0,DY=0.0,DZ=0.0,DRY=0.0),
_F(GROUP_MA=('Liney0','Liney1'),
DX=0.0,DY=0.0,DZ=0.0,DRX=0.0),),);
else:
DDL_IMPO=(_F(GROUP_MA=('Linex0','Linex1'),
DY=0.0,DZ=0.0,),
_F(GROUP_MA=('Liney0','Liney1'),
DX=0.0,DZ=0.0),),);
```
```#For Codes Aster version 11.3 see next chapter from here on:
```
```# define stiffness- and mass matrix and dofs
MACRO_MATR_ASSE(MODELE=cocModel,
CHAM_MATER=material,
CARA_ELEM=shellch,
NUME_DDL=CO('NUMEDDL'),
MATR_ASSE=(_F(MATRICE=CO('RIGIDITE'),
OPTION='RIGI_MECA',),
_F(MATRICE=CO('MASSE'),
OPTION='MASS_MECA',),),);
```
```# frequency/mode shapes determination within band of 1 and 30 Hz
modes=MODE_ITER_SIMULT(MATR_A=RIGIDITE,
MATR_B=MASSE,
CALC_FREQ=_F(OPTION='BANDE',
FREQ=(1.0,30.0,),),);
```
```# no need for this
stress=CALC_ELEM(RESULTAT=modes,
OPTION='SIEF_ELNO',);

# project results fo modes from CA-mesh (quad9)
salModes=PROJ_CHAMP(RESULTAT=modes,
MODELE_1=cocModel,
MODELE_2=intModel,);
```
```# write results of converted/projected results
# both salModes (projected results) and modes (original results) are written to file
# modes has zero displacements at centre nodes
# both results (salModes and modes) contain the mesh information as well so no need to define thes in the RESU-part.
IMPR_RESU(FORMAT='MED',
UNITE=80,
RESU=(_F(RESULTAT=salModes,
INFO_MAILLAGE='NON',
NOM_CHAM='DEPL',),
_F(INFO_MAILLAGE='NON',
RESULTAT=modes,
NOM_CHAM='DEPL'),),);
```
```# determine total mass of shell and print to file
masshell=POST_ELEM(MASS_INER=_F(GROUP_MA='shell',ORIG_INER=(0.5, 0.5, 0.0),),
MODELE=cocModel,
CHAM_MATER=material,
CARA_ELEM=shellch,);
```
```IMPR_TABLE(TABLE=masshell,UNITE=22,NOM_PARA='MASSE',);
```
```FIN();
```

The comm file (see end of this contribution) may be edited by Eficas [under Tools --> Command file editor (Eficas)] in the ASTK menu bar. You get a number of possible entries for each command.

### The command file - detailed commands - CA version 11.3

First part of the command file copy from CA version 11.1 until where indicated.

```# see also SDLL123b.comm
ASSEMBLAGE(MODELE=cocModel,
CHAM_MATER=material,
CARA_ELEM=shellch,
NUME_DDL=CO('NUMEDDL'),
MATR_ASSE=(_F(MATRICE=CO('Mstiff'),OPTION='RIGI_MECA',),
_F(MATRICE=CO('Mmass'),OPTION='MASS_MECA',),),);
```
```MATK=Mstiff
MATM=Mmass
```
```# CALCUL VIA QZ
## Very expensive in term of in CPU and memory.
## Reserved for small cases (<1000dofs)
## Range into asymmetric and with A symmetric complex
### doesnot run on my pc (need to add memory)
if (0==1):
MG_QZ=MODE_ITER_SIMULT( INFO=1,
MATR_RIGI=MATK,
MATR_MASS=MATM,
METHODE='QZ',
CALC_FREQ=_F(OPTION='CENTRE',FREQ=15.,NMAX_FREQ=15),);
```
```# CALCUL VIA SORENSEN
## default choice
modes=MODE_ITER_SIMULT( INFO=1,
MATR_RIGI=MATK,
MATR_MASS=MATM,
TYPE_RESU=‘DYNAMIQUE’
METHODE='SORENSEN',
CALC_FREQ=_F(OPTION='BANDE',FREQ=(1.0,50.0),),);
#CALC_FREQ=_F(OPTION='CENTRE',FREQ=15.,NMAX_FREQ=15),);

# project results fo modes from CA-mesh (quad9 with centre nodes)
# to Salome mesh (quad8 without centre nodes)
salModes=PROJ_CHAMP(RESULTAT=modes,
MODELE_1=cocModel,
MODELE_2=intModel,);
```
```IMPR_RESU(FORMAT='MED',
UNITE=80,
RESU=(_F(RESULTAT=salModes,INFO_MAILLAGE='NON',NOM_CHAM='DEPL',),
#_F(RESULTAT=modes,INFO_MAILLAGE='NON',NOM_CHAM='DEPL'),
),);
```

Notes:

• ASSEMBLAGE

Here the stiffness and mass matrices are formed. The boundary conditions are taken into account by the NUME_DDL=CO('NUMEDDL'), part.

• MODE_ITER_SIMULT

The eigenvalues and modes shapes are determined based on the stiffness and mass matrices. The most general and default methode is SORENSEN. NOte that the method QZ is both cpu and memory intensive and even for this simple case does not run on my computer due to lack of memory. The options for determining the frequency range are:

• CALC_FREQ=_F(OPTION='BANDE',FREQ=(1.0,50.0),),)
• search resonance frequencies between 1.0 and 50.0 [Hz]
• CALC_FREQ=_F(OPTION='CENTRE',FREQ=15.0,NMAX_FREQ=15),)
• search maximum 15 resonance frequencies around 15.0 [Hz]
• CALC_FREQ=_F(OPTION='PLUS_PETITE',NMAX_FREQ=3)
• this is default choice, search lowest 3 resonance frequencies
• CALC_FREQ=_F(OPTION='TOUT')
• only interesting/allowed? for METHODE='QZ': On cherche tous les modes associés à des ddls physiques. Option utilisable qu’avec la méthode QZ.

The keyword TYPE_RESU is used to indicate a frequency ‘DYNAMIQUE’ or bluckling analysis:

• TYPE_RESU=‘DYNAMIQUE’ (default value, frequency) or
• TYPE_RESU=‘MODE_FLAMB’ (buckling factor with respect to applied load, probably multiplied by -1)
• PROJ_CHAMP

Due to limitations of the postprocessor of Salome we have to project the displacements of the TRIA7/QUAD9 to TRIA6/QUAD8 meshes and finally

• IMPR_RESU

print the result to file, in this case the displacements DEPL in med-format on unite 80, see ASTK file.

• RESU=_F(RESULTAT=salModes,INFO_MAILLAGE='NON',NOM_CHAM='DEPL',)
• projected displacements to tria6/quad8 mesh, usefull --> model initMode, mesh initMesh
• RESU=_F(RESULTAT=modes,INFO_MAILLAGE='NON',NOM_CHAM='DEPL')
• displacements of tria7/quad9 mesh, no use anymore --> model cocModel, mesh comMesh

## Post processing with Salome - Based on CA11.1

Note1 - update:
These commands and results are based on Code Aster version 11.1. In version 11.2&3 commands have been greatly changed.

Note2 - update:
Three result files have been written on units 80, 81 and 82, the names defined by ASTK. Select eg. the result file from unit 81.
Return back to Salome, select Post Pro and select by File --> Import --> res33.med ie. the result file defined by ASTK (res33.med). Now the nice part of the work starts by making some colourful pictures. In the left window click on Post Pro and the + res33 appears. Click onthe plus (+) sign and the two meshes appear: meshinit and meshmod. Each have some subfields. The most interesting is the field Fields in meshmod. Click again on the + in front of the Fields and the frequencies of the MODES___DEPL___ appear. Now right click on the + sign in front of one of the frequencies and play around with the settings of the template window and modes shapes appear in the VTK viewer window. Have fun with different setting, by right clicking on the Def.Shape field below the frequency.

A view of the first mode shape is given here: : Note:
The picture on the left shows the displacement of the shell for the projected salome modes (salmodes) whereas the picture on the right show the displacement with zero displacement at the centre nodes of the elements.

## Post processing with Salome - Based on CA11.3

Just a few pictures to show the results. Nothing fancy. The egg-holder is a gift.
Note that the no results are displayed for the mesh with centre nodes, see right picture for CA11.1. So the projection is needed to see the displacements. : frequencies: 4.91 (1x1 shape) and 44.24 [Hz] (3x3 shape)

### Comparison to analytic results

[update --> need to recalculate with nu=0.3]
All four boundaries of the plate are hinged. Comparing the first five analytic (Fth) to the numerical (Fnum) results [Hz] and relative differences, gives the following table:

``` shape       Fth        Fnum         dif %
1*1       7.429       7.425        -0.054
1*2      14.757      14.751        -0.041
2*1      22.389      22.381        -0.036
1*3      26.974      26.963        -0.041
2*2      29.715      29.702        -0.044
```

In total 8 eigenvalues are present in the range below 50 Hz, from 7 to is 47 Hz. The analytical results are according to G.B. Warburton, the vibration of rectangular plates, 1954. [??]

```        1      8.90783E+00       8.61615E-09
2      1.81686E+01       1.98617E-09
3      1.81686E+01       2.09953E-09
4      2.67922E+01       9.84001E-10
```

## Total mass of the plate

Finally we add a command to determine the mass of the shell. At the end of the command file add the following lines:

```masshell=POST_ELEM(MASS_INER=_F(GROUP_MA='shell',),
MODELE=modmod,
CHAM_MATER=material,
CARA_ELEM=shellch,);
```
```IMPR_TABLE(TABLE=masshell,NOM_PARA='MASSE',);
FIN();
```

The results are written to a table in the resu file, in this case the shell33.resu.med file.

```#ASTER  9.04.06 CONCEPT masshell CALCULE LE 14/05/2009 A 11:54:43 DE TYPE
#TABLE_SDASTER
MASSE
7.85000E+00
```

Standard the tables are written to unite 8, the result file. To write to a seperate file, add a line in ASTK with type 'resu', filename ./masse.txt (in Linux) and unit number (LU) e.g. 22, change the command line into:

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

and the results are nicely written to the file masse.txt seperately from all messages in the result file.

[You can also get the centre of gravity and the inertia of the construction, but I did succeed in doing that.]

## Input and output files for the FE Analysis

Input files:

• Python script for defining the geometry (geom_plate1by1.py)
• Python script for defining the mesh (mesh_plate1by1.py)
• python script of geometry and mesh (combination of those above, geom_mesh_plate1by1.py)
• ASTK file (shelldyn.astk, you need to edit the path to your requirements ...)
• command file (shelldyn.comm)

output file:

• total mass returned by CA (masse.txt)

Update:
In the latest zip a slightly different naming convention has been used. The meshes and models are renamed to init.... and coc... for initial mesh/model and converted to coque_3d mesh/model. The results are called modes for coque-3d-modes and salmodes for quad8/tria6 modes.

## Differences between version MODI_MAILLE

### Code Aster version 11.X - TRIA6_7 and QUAD8_9

In Version 11.X you can once again use OPTION='QUAD8_9' combined with OPTION='TRIA6_7' even though only one of quad8 or tria6 is available. A warning is issued. In this case only quad8 elements are present in the mesh:

```cocMesh=CREA_MAILLAGE(MAILLAGE=initMesh,
INFO=1,
MODI_MAILLE=(_F(PREF_NOEUD='NS',
PREF_NUME=1,
TOUT='OUI',
_F(PREF_NOEUD='NS',
PREF_NUME=1,
TOUT='OUI',
OPTION='TRIA6_7')),);
```

MOT CLE FACTEUR "MODI_MAILLE", OCCURRENCE 1

``` MODIFICATION DE    729 MAILLES    QUAD8 EN    QUAD9

!----------------------------------------------------------------!
! <A> <MODELISA3_32>                                             !
!                                                                !
!  Alarme utilisateur dans CREA_MAILLAGE/MODI_MAILLE :           !
!   Occurrence du mot clé facteur MODI_MAILLE : 2.               !
!   Vous avez demandé la transformation de type TRIA6_7.         !
!   Mais il n'y a aucune maille à transformer.                   !
!                                                                !
!                                                                !
! Ceci est une alarme. Si vous ne comprenez pas le sens de cette !
! alarme, vous pouvez obtenir des résultats inattendus !         !
!----------------------------------------------------------------!
```

### Code Aster version 10.X - TRIA6_7 and QUAD8_9

```shellch=AFFE_CARA_ELEM(MODELE=modmod,
COQUE=_F(GROUP_MA='shell',
EPAIS=th,
COQUE_NCOU=1,),);
```
```meshmod=CREA_MAILLAGE(MAILLAGE=meshinit,
```

CREA_MAILLAGE
The problem is that no TRIA6 elements are present in the mesh. So version Fea_plate1by1.zip works on Code-Aster 9.xx but not on version 10.xx. In this example only Quad8 elements are present. The MODI_MAILLE keyword in the CREA_MAILLAGE comnmand in C-Aster 10 does not accept this part: _F(TOUT='OUI',OPTION='TRIA6_7',) anymore.

```....MODI_MAILLE=(_F(TOUT='OUI',OPTION='QUAD8_9',),
_F(TOUT='OUI',OPTION='TRIA6_7',),) ....;
```

If both TRIA6 and QUAD8 elements are present use:

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

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

If only TRIA6 elements (triangles) are available use:

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

The message file reports QUAD8 elements not available, but I guess it it the absence of TRIA6 elements that causes a F_COPY_ERROR:

```MOT CLE FACTEUR "MODI_MAILLE", OCCURRENCE    1