Contrib:KeesWouters/plate/thickness

From CAELinuxWiki
Revision as of 20:56, 21 November 2010 by Keeswouters (Talk | contribs)

(diff) ← Older revision | Latest revision (diff) | Newer revision → (diff)
Jump to: navigation, search

Using coque_3d shell elements

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.28 [-] and density RHO=7850 kg/m3.

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

Versions of the software that has been used (april 2009):

  • Salome 4.1.4 (32 bit)
  • Code Aster 9.4 (32 bit)
  • ASTK 1.7.1
  • OS: Ubuntu 9.04 Jaunty (64 bit)

Update Code Aster version 10.1

All files can be downloaded at the end of this contribution.
Changes between CA-9 and CA-10 are discussed also at the end of this contribution.
The result file is now written to the standard file unit 80, type rmed.

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
  • 2D algorithm: quadrangle (Mapping), hypothesis: quadrangle preference
  • 1D algorithm: wire discretisation1D: hypothesis automatic length (set to 1, or very fine)
  • apply and close

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:

Kw astk shelldyn.png

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

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)

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:
  1. Salome: set quadratic mesh directly
  2. Salome: convert to quadratic mesh in Salome after linear mesh has been created
  3. 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 21, corresponding to the LU column in the res file of ASTK (see again Media:ASTKvs171shell33.png).


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

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:

Kw first mod shell.jpg

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:

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)

Retrieved from "http://caelinux.org/wiki/index.php?title=Contrib:KeesWouters/plate/thickness&oldid=5972"