Difference between revisions of "Contrib:KeesWouters/plate/thickness"

Revision as of 22:19, 21 November 2010

Varying thickness of shells - by Python

Static displacement of a shell under pressure

november 2010 - SalomeMeca 2010 (Code Aster 10.2) - Ubuntu Meerkat 10.10 64bits.

This is just a simple example of the use of shells (coque_3d) for static calculations. The main focus is on varying thickness of the shell by a Python list. The thickness of the shell depends on the average Ycog coordinate of the shell element.

General way of determination of the shell thickness

First the geometry and mesh is generated in Salome in the standard way. In Code Aster the mesh is read, including nodes number, node coordinates, element types and element connectivity. We use quadratic quad8 elements that are converted to quad9 (coque_3d) elements suitable for C Aster. Of each element the nodes are determined by the connectivity matrix, the coordinates of all the nodes and finally the centre of gravity. This includes the y coordinate Ycog. The element thickness th is calculated by th = Ycog*0.1+0.01 and stored in a Python list. This list is passed to the AFFE_CARA_ELEM command. Further processing C Aster is from then on as usual.

Geometry, mesh and material properties

The construction is a rectangle in the xyz plane. Minimum x is 0, maximum x is 1.0 [m]. For y and z these values are: [ymin, ymax] = [0.0, 2.5] and [zmin, zmax] = [1.25, 2.3456]. The varying thickness of the plate depends on the average y coordinate as follows: th = Ycog*0.1+0.01.
Boundary conditions are as follows: the two edges parallel to the y axis are restricted in all six directions. The material is steel, its mechanical properties are defined by the following parameters: Young's modulus E = 2.1e11 Pa, poisson ratio NU=0.28 [-] and density -though not needed- here RHO=7850 kg/m3. For more details see e.g. here where the details of the general calculation procedure are given. In the next part we concentrate on the extraction of the coordinates of the elements and the definition of the thickness of each element.

The command file for Code Aster - overview

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 meshinit)
• convert the mesh to be suitable for the shell (coque) elements (meshmod)
• define the model on the modified mesh (modmod)
• 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 (meshmod)
• convert the result to the original mesh to be read by Salome (meshinit)

The command file - detailed commands

DEBUT();
# reading med mesh (code aster, salome, default)
meshinit=LIRE_MAILLAGE(FORMAT='MED',INFO=2,);

modinit=AFFE_MODELE(MAILLAGE=meshinit,
                   AFFE=_F(TOUT='OUI',
                           PHENOMENE='MECANIQUE',
                           MODELISATION='3D',),);

# define shell (coque_3D) model
# we have to add some centre nodes with CREA_MAILLAGE
# this command converts the Salome mesh into a mesh suitable for the coque_3D model
#  in this case only quad (quadrangle) elements are used, no tria (triangular) elements
meshmod=CREA_MAILLAGE(MAILLAGE=meshinit,
                     MODI_MAILLE=(_F(TOUT='OUI',OPTION='QUAD8_9',),),);
                                  ####_F(TOUT='OUI',OPTION='TRIA6_7',),),); not compatable with C-Aster version 10

# define a model suitable for shell coque_3D elements
modmod=AFFE_MODELE(MAILLAGE=meshmod,
                  AFFE=_F(TOUT='OUI',
                          PHENOMENE='MECANIQUE',
                          MODELISATION='COQUE_3D',),);

# the shell characteristics
# the shells are applied on the group 'shell' defined in the Salome mesh
# the thickness of the shell is epais=0.001
# the angle_rep is a direction angle, use either angle(a,b) or vecteur(x,y,z)
# coque_ncou is the number of gauss nodes along the thickness, for linear analysis one node is sufficient.
# the parameter excentrement can give the offset of the shell wrt to the meshing plane, but this does not apply for coque_3d and is ignored here
shellch=AFFE_CARA_ELEM(MODELE=modmod,
                      COQUE=_F(GROUP_MA='shell',
                               EPAIS=0.001,
                               ####ANGL_REP=(0.0,0.0,), defines local x axis
                               ####VECTEUR=(1.0,0.0,0.0),
                               COQUE_NCOU=1,
                               EXCENTREMENT=0.000),);

# define material properties of steel (ISO values)
steel=DEFI_MATERIAU(ELAS=_F(E=2.100e11,NU=0.28,RHO=7850.0,),);

# apply material properties to the whole mesh
material=AFFE_MATERIAU(MAILLAGE=meshmod,AFFE=_F(TOUT='OUI',MATER=steel,),);

# define BC and loads
# for all the edges of the plate the z displacement is fixed [DZ=0.0]
# one edge [Linex0] is kept fixed in y direction [DY=0.0] and
# another edge [Liney0] is kept fixed in x direction [DX=0.0]
Loaddyn=AFFE_CHAR_MECA(MODELE=modmod,
                      DDL_IMPO=(_F(GROUP_MA='Linex0',
                                   DY=0.0,
                                   DZ=0.0,),
                                _F(GROUP_MA='Liney0',
                                   DX=0.0,
                                   DZ=0.0,),
                                _F(GROUP_MA=('Liney1','Linex1',),
                                   DZ=0.0,),),);

# define the calculation
# determine stiffness and mass matrix depending on material, mesh and applied  model (coque_3D) for dynamic calculation (mode shapes and frequencies)
MACRO_MATR_ASSE(MODELE=modmod,
               CHAM_MATER=material,
               CARA_ELEM=shellch,
               CHARGE=Loaddyn,
               NUME_DDL=CO('NUMEDDL'),
               MATR_ASSE=(_F(MATRICE=CO('RIGIDITE'),
                             OPTION='RIGI_MECA',),
                          _F(MATRICE=CO('MASSE'),
                             OPTION='MASS_MECA',),),);

# determine the frequencies within the band 1 and 50 Hz
MODES=MODE_ITER_SIMULT(MATR_A=RIGIDITE,
                      MATR_B=MASSE,
                      CALC_FREQ=_F(OPTION='BANDE',
                                   FREQ=(1.0,50.0,),),);

# determine / project the result fields (champ) from modified model modmod to the initial model modinit
salomres=PROJ_CHAMP(RESULTAT=MODES,
                   MODELE_1=modmod,      #project modes/results of model_1
                   MODELE_2=modinit,);   #to model_2

# store the results in a MED file
# only the displacements 'DEPL' will be stored
# note that UNITe 21 corresponds to the LU column of ASTK.
IMPR_RESU(FORMAT='MED',
         UNITE=80,
         RESU=_F(MAILLAGE=meshinit,
                 RESULTAT=MODES,
                 NOM_CHAM='DEPL',),);
#fin
FIN();

The comm file (see end of this contributionn) 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.

Post processing with Salome

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:

Comparison to analytic results

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.

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)

download the files here:

• >>>Media:Fea_plate1by1.zip: error on TRIA6_7 in Code Aster10, please do not use anymore
• >>>Media: kw_platedyn_ca10.zip: suitable for Aster10

CODE ASTER version 10.X - TRIA6_7 and QUAD8_9

Update in the second version:

shellch=AFFE_CARA_ELEM(MODELE=modmod,
                      COQUE=_F(GROUP_MA='shell',
                               EPAIS=th,
                               COQUE_NCOU=1,),);

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

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_MAILLAGE 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,
                      MODI_MAILLE=(_F(TOUT='OUI',OPTION='QUAD8_9',),
                                   _F(TOUT='OUI',OPTION='TRIA6_7',),),);

If only QUAD8 elements (quadrangles) are available use:

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

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
 MODIFICATION DE    729 MAILLES    QUAD8 EN    QUAD9
! <EXCEPTION> <MODELISA3_32>                               !
!  Erreur utilisateur dans CREA_MAILLAGE/MODI_MAILLE :     !
!   Vous avez demandé la transformation de type QUAD8_9.   !
!   Mais il n'y a aucune maille quadratique à transformer. !

References

From Code Aster references:
[u2.02.01] Notice d’utilisation des éléments plaques, coques et coques volumiques SHB
[u4.81.22] Opérateur POST_ELEM
[u5.02.00] Manuel d’Utilisation - Fascicule U5.0- : Structure de données resultat Document: U5.02.00 - Présentation des tables
[u3.12.03] Modélisation COQUE_3D
[R3.07.04] Eléments finis de coques volumiques
[R3.07.05] Éléments de coques volumiques en non linéaire géométrique

Overview of units - Paul Carrico

--keeswouters 17:11, 24 April 2009 (CEST)