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

From CAELinuxWiki
Jump to: navigation, search
m (→‎'''Using coque_3d shell elements''')
m (→‎General - extraction of coordinates)
 
(126 intermediate revisions by the same user not shown)
Line 1: Line 1:
 
=='''Varying thickness of shells - by Python'''==
 
=='''Varying thickness of shells - by Python'''==
−
==='''Static displacement of a shell under pressure'''===
+
If you are already familiar to Code Aster, enjoy this. If you are new to Code Aster have a look [http://www.caelinux.org/wiki/index.php/Contrib:KeesWouters/platedynamics here] first please.<br/>
−
 
+
 
november 2010 - SalomeMeca 2010 (Code Aster 10.2) - Ubuntu Meerkat 10.10 64bits.
 
november 2010 - SalomeMeca 2010 (Code Aster 10.2) - Ubuntu Meerkat 10.10 64bits.
 +
<br/>
 +
January 2013 - update to Code Aster version 11.3 - Salome 6.6.0 - Linux Mint 14 Nadia<br/>
  
 +
For a similar contribution with variable pressue see: [[Contrib:KeesWouters/plate/variable_pressure]]
 +
 +
==='''Static displacement of a shell under pressure'''===
 
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.
 
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.
  
 +
Two meshes are used throughout this contribution: a course mesh (3 by 4 elements) and a fine mesh. In the geometry and mesh file (download as at the end) the meshsize can be set by the variable ''meshsize'' (number between 0.0 and 1.0). The course mesh I have been using for experimenting. The fine mesh to show more reasonable results.
  
 
==='''General way of determination of the shell thickness'''===
 
==='''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.
+
The geometry and mesh are 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.<br/>
−
 
+
Of each element  
 +
* determine the nodes are defined by the connectivity matrix,  
 +
* determine the coordinates of all the nodes and finally
 +
* determine the centre of gravity and
 +
* generate the thickness of the shell
 +
by Python commands.
 +
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'''===
 
==='''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.<br/>
+
The construction is a rectangle in the xyz plane. The minimum and maximum x, y and z values are:
−
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.
+
* x: [xmin, xmax] = [0.00, 1.00]
 +
* y: [ymin, ymax] = [0.00, 2.50]
 +
* z: [zmin, zmax] = [1.25, 2.35]
 +
The thickness of the plate we make dependent on the y coordinate of the centre of the element as follows: th = Ycog*0.1+0.01 (edit this for your own definition here of course).<br/>
 +
The  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 of the general calculation procedure see e.g. [http://www.caelinux.org/wiki/index.php/Contrib:KeesWouters/platedynamics here.] 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 geometry of the plate'''==
+
: [[Image:kw_varshell_geom.jpg]]
  
−
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):
+
=='''Extraction of the coordinates'''==
−
* p00 --> [0.0, 0.0, 0.0]
+
Among others, in [http://www.code-aster.org/forum2/viewtopic.php?id=13692 Code Aster forum] the procedure is discussed.<br/>  
−
* p10 --> [1.0, 0.0, 0.0]
+
In Code Aster, after reading the initial mesh, the quad8 elements are converted to quad9 (coque_3d compatible) elements. This mesh is called meshmod (''modified mesh'') and is used as the basis for the further operations.  
−
* 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:
+
First we adapt the DEBUT() command since we need ''real'' Python commands:
−
*Ly0 --> p00 and p10 (enter name Ly0 and select the two points)
+
DEBUT(PAR_LOT='NON');
−
*Lx1 --> p10 and p11
+
We import the Python '''partition''' module in the command file to read mesh quantities. The definition of the thickness of the shells is stored in a separate Python script/module/file (construct_shell_thickness) in the path defined below (change for your own needs [sys.path.append('/your/path/shellplate_thlist')]):
−
*Ly1 --> p11 and p01
+
from Utilitai.partition import *
−
*Lx0 --> p01 and p00
+
import sys
 +
sys.path.append('/cae/caexample/shell/shellplate_thlist')
 +
from construct_shell_thickness import *
  
−
and create the face defined by the four lines by clicking New Entity --> Build --> Face
+
After reading the initial mesh '''meshinit''' with quad8 elements
−
*Face1 --> Ly0, Lx1, Ly1, Lx0 (shift select the four lines and click Apply).
+
meshinit=LIRE_MAILLAGE(FORMAT='MED',INFO=1,);
 +
and converting it to quad9, coque_3d elements
 +
meshmod=CREA_MAILLAGE(MAILLAGE=meshinit,
 +
                      MODI_MAILLE=(_F(TOUT='OUI',
 +
                                      OPTION='QUAD8_9',),),);
  
−
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 define the thickness of each element by calling
−
We select the lines by clicking New Entity --> Group --> Create:
+
info = 1  ## true/false - print additional information
−
* keep default choices except:
+
ThicknessGroupShell = construct_shell_thickness(meshmod,info)
−
* 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.  
+
The Python script '''construct_shell_thickness''' does the work.<br/>
−
We select the Face1 by clicking New Entity --> Group --> Create:
+
We have a look at this in detail. Standard, the characteristics of the shell are defined by the AFFE_CARA_ELEM command with, among others, the thickness parameter EPAIS:
−
* keep default choices except:
+
shellch=AFFE_CARA_ELEM(MODELE=modmod,
−
* shape type: select plane
+
                      COQUE=_F(GROUP_MA='shell',
−
* main shape: Face1,  
+
                                EPAIS=0.001,),);
−
* name: shell, select all, and apply and close.
+
  
−
So far for the definition of the geometry.
+
Now remember that the parameters in the _F(...) operator can be replaced by a Python list, in this case: ellist = {'GROUP_MA': 'shell','EPAIS':0.001}. Or more general:
−
The geometry is also available by the following python script
+
  ellist = [{'MAILLE': 'Mxxx','EPAIS':thxxx},
−
[[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.
+
            {'MAILLE': 'Myyy','EPAIS':thyyy},                     
 +
              ::::
 +
            {'MAILLE': 'Mzzz','EPAIS':thzzz}]
  
−
=='''The meshing of the plate'''==
+
where 'Mxxx' is the element number and thxxx is the corresponding element thickness. Eg in this example the following list ''thicknessGroupShell'' is generated:
 +
  [{'EPAIS': 0.16625, 'MAILLE': 'M15'},
 +
  {'EPAIS': 0.22875, 'MAILLE': 'M16'},
 +
    ::::
 +
  {'EPAIS': 0.22875, 'MAILLE': 'M25'},
 +
  {'EPAIS': 0.10375, 'MAILLE': 'M26'}]
  
−
Switch to the meshing module in Salome. We start with a quadrangle mesh, select Mesh --> Create Mesh. Apply the following settings:
+
==='''General - extraction of coordinates'''===
−
* name: shell
+
In this example, the actual thickness of the shell depends on the y coordinate at the centre of the element. This is extracted as follows. First read the mesh:
−
* geometry: Face1
+
meshCA = MAIL_PY()
−
* 2D algorithm: quadrangle (Mapping), hypothesis: quadrangle preference
+
meshCA.FromAster(mesh)
−
* 1D algorithm: wire discretisation1D: hypothesis automatic length (set to 1, or very fine)
+
In '''meshCA''' the available mesh quantities can be extracted (nodes, elements, connectivity, coordinates):
−
* apply and close
+
  nonu    = meshCA.dime_maillage[0]    ## number of nodes
 +
  elnu    = meshCA.dime_maillage[2]    ## number of elements
 +
  ncoord  = meshCA.cn                  ## coordinates of the nodes
 +
  connect  = meshCA.co                  ## connectivity of the element
 +
  NodeList = list(meshCA.correspondance_noeuds)   ## says it all
 +
  ElemList = list(meshCA.correspondance_mailles)
 +
  ElemType = list(meshCA.tm)
  
−
Now the shell mesh appears in the object browser, right select and hit Compute. A mesh with about 100 edges and 700 quadrangles appear.
+
The procedure is then as follows:
 +
* step through element list
 +
** check if element type is quad9/coque_3d/number 16
 +
** read connectivity of the element, ie all attached node numbers: connectivity = meshCA.co
 +
** step through all node of the element
 +
*** determine nodes and its coordinates
 +
*** compute average (y) coordinate
 +
*** define thickness depending on (y) coordinate
 +
*** add thickness and element to list {...}
  
−
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.
+
In 'real' Python code:
 +
for ii in xrange(len(ElemList)):
 +
  print ElemType[ii]
 +
  if ElemType[ii]==16:  # select shell (coque_3d) elements
 +
    #print connect[ii]
 +
    shellcount +=1
  
−
==='''Export the mesh to a MED2.3 file'''===
+
ThicknessGroupShell=[]  # initialise variables
−
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.
+
thickness = [None]*shellcount
  
−
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
+
shellcount = 0  # reset shellcount
 +
for ii in xrange(len(ElemList)):
 +
  print 'elementtype: ',ii,ElemType[ii]
 +
  if ElemType[ii]==16:           # select shell elements coque_3d (14:quad8 but only 16:quad9 available)
 +
    sumxyz = [0.0, 0.0, 0.0]
 +
    ncount = len(connect[ii])    # should be 7 for tria7 and/or 9 quad9 for coque_3d (elementType 16)
 +
    for jj in xrange(0,ncount):
 +
      Gnode = connect[ii][jj]    # globalnode = connect[elemii,internal_nodejj]
 +
      sumxyz += ncoord[Gnode]    # this might be written more compact --> for later
  
−
=='''ASTK'''==
+
    sumxyz /=ncount
 +
    # define your own code for thickness here
 +
    # thickness now depends on average y coordinate only
 +
    thickness[shellcount] = sumxyz[1]*0.01+0.1  # make thickness dependending on y coordinate
 +
    tempth = thickness[shellcount]
 +
    shellcount +=1
 +
    #print 'average xyz - centre point: ',sumxyz,ncoord[Gnode]  #centre point is same average nodes
 +
    ThicknessGroupShell.append({'EPAIS':tempth,'MAILLE':'M%d'%(ii+1)});  ## add 1 for element number Mxx
  
−
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:
+
If needed, vecteur, angle or NCOU can be defined as well.
  
−
* ...
+
==='''Things that might go wrong'''===
−
* F mail kees@localhost:/home/kees/shell/shell33.med D 20 (input mesh file)
+
Of course, nothing will ever go wrong. Except code written in Python and used in conjunction with Code Aster. Then things definitely will go wrong. At least the first time by a unstructered programmer (like me?).
−
* F comm kees@localhost:/home/kees/shell/shell33dyn.comm D 1 (input command file)
+
First all kind of syntax errors occur. These are relatively easy to address and remove, because Python offers help by issuing nice error messages. Things are less trivial if you forget the +1 to inidicate the mesh element in the append command: ''' 'MAILLE':'M%d'%(ii+1)''' because of the 1 offset numbering of Salome to Code Aster mesh. To see how the mesh is built up you can print the whole mesh FromAster:<br/>
−
* F resu kees@localhost:/home/kees/shell/shell33.resu.med R 8 (result file)
+
meshCA = MAIL_PY()
−
* F mess kees@localhost:/home/kees/shell/shell33.mess R 6 (messages file)
+
meshCA.FromAster(mesh)
−
* F erre kees@localhost:/home/kees/shell/shell33.erre R 9 (error file)
+
print meshCA,
−
* 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:  
+
Then you get eg:
−
: [[Image:kw_astk_shelldyn.png]]
+
COOR_3D
−
 
+
N1        1.00000000000000E+00  2.50000000000000E+00  2.34560000000000E+00
−
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.
+
SEG3
−
 
+
  M1        N3      N17      N18   
−
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:
+
−
# 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 ''unit''e 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>
+
  QUAD9
 +
M15      N10    N29    N30    N12    N36    N37    N38    N13    NS1 
 +
M16      N7    N5    N32    N29    N8    N42    N43    N35    NS2   
 
  ....
 
  ....
−
  Objet JEVEUX inexistant dans les bases ouvertes .....
+
  GROUP_MA
 +
Linex0
 +
M5      M7      M10      M11   
 +
FINSF
 +
%
 +
GROUP_MA
 +
Linex1
 +
M3      M9      M12      M14   
 +
FINSF
 +
%
 +
GROUP_MA
 +
shell
 +
M15      M16      M17      M18      M19      M20      M21      M22   
 +
M23      M24      M25      M26   
  
−
===''The command file - detailed commands''===
+
This allows you to check the correct node numbers, elements, connectivity and groups defined in your mesh. Starting with a small mesh is advicable.<br/>
−
DEBUT();
+
Also carefully check the code that generates the thickness of the shell. It is very easy to make a mistake: selecting the wrong index of the coordinate array is an obvious one.<br/>
−
# reading med mesh (code aster, salome, default)
+
For a small mesh (3 by 4 elements) the Python list is as follows:
−
meshinit=LIRE_MAILLAGE(FORMAT='MED',INFO=2,);
+
  
−
  modinit=AFFE_MODELE(MAILLAGE=meshinit,
+
  thicknessGroupShell: 
−
                    AFFE=_F(TOUT='OUI',
+
[{'EPAIS': 0.16625, 'MAILLE': 'M15'},  
−
                            PHENOMENE='MECANIQUE',
+
  {'EPAIS': 0.22875, 'MAILLE': 'M16'},  
−
                            MODELISATION='3D',),);
+
  {'EPAIS': 0.04125, 'MAILLE': 'M17'},  
 +
  {'EPAIS': 0.22875, 'MAILLE': 'M18'},
 +
  {'EPAIS': 0.16625, 'MAILLE': 'M19'},
 +
  {'EPAIS': 0.10375, 'MAILLE': 'M20'},
 +
  {'EPAIS': 0.10375, 'MAILLE': 'M21'},
 +
  {'EPAIS': 0.04125, 'MAILLE': 'M22'},
 +
  {'EPAIS': 0.04125, 'MAILLE': 'M23'},
 +
  {'EPAIS': 0.16625, 'MAILLE': 'M24'},
 +
  {'EPAIS': 0.22875, 'MAILLE': 'M25'},
 +
  {'EPAIS': 0.10375, 'MAILLE': 'M26'}]
  
−
# ''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
+
'''Remarks'''<br/>
−
modmod=AFFE_MODELE(MAILLAGE=meshmod,
+
Note that you could select the midplane node '''NSxxx''' just as well to make the thickness dependending on (NSxxx is added to the quad8 element by the  modi_maillage command). Then there is no need to compute the average values of the coordinates of the 8 nodes on the egde.
−
                  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)
+
: [[Image:kw_varshell_3x4.jpg]]
−
steel=DEFI_MATERIAU(ELAS=_F(E=2.100e11,NU=0.28,RHO=7850.0,),);
+
  
 +
=='''Continue with the commands'''==
 +
Now we follow the standard procedure for calculating the static displacement of the shell under pressure. Normally we define the the thickness of the shell by:
  
−
  # apply material properties to the whole mesh
+
  th = 0.02
−
  material=AFFE_MATERIAU(MAILLAGE=meshmod,AFFE=_F(TOUT='OUI',MATER=steel,),);
+
  shellch=AFFE_CARA_ELEM(MODELE=modelc,
 +
                      COQUE=_F(GROUP_MA='shell',
 +
                                EPAIS=th,),);
 +
                                #VECTEUR=(1.0,0.0,0.0,),
 +
                                #COQUE_NCOU=1,),);
  
−
# ''define BC and loads''
+
This will be replaced by:
−
# for all the edges of the plate the z displacement is fixed [DZ=0.0]
+
  shellch=AFFE_CARA_ELEM(MODELE=modelc,COQUE=(ThicknessGroupShell),);
−
# one edge [Linex0] is kept fixed in y direction [DY=0.0] and
+
Very compact, I would say. Unfortunately the whole process defined in '''construct_shell_thickness''' hides the effort. The MECA_STATIQUE command carries out the static analysis:
−
# another edge [Liney0] is kept fixed in x direction [DX=0.0]
+
result=MECA_STATIQUE(MODELE=modelc,
−
  Loaddyn=AFFE_CHAR_MECA(MODELE=modmod,
+
                    CHAM_MATER=material,
−
                      DDL_IMPO=(_F(GROUP_MA='Linex0',
+
                    CARA_ELEM=shellch,
−
                                    DY=0.0,
+
                    EXCIT=_F(CHARGE=clamped,),);
−
                                    DZ=0.0,),
+
−
                                _F(GROUP_MA='Liney0',
+
−
                                    DX=0.0,
+
−
                                    DZ=0.0,),
+
−
                                _F(GROUP_MA=('Liney1','Linex1',),
+
−
                                    DZ=0.0,),),);
+
  
−
# ''define the calculation''
+
where of course the model modelc, the material properties, the shell characteristics shellch (oeps) and the charge ''clamped'' -being boundary conditions in this case- have been defined in advance. After calculating the element results and printing the results by
−
# determine stiffness and mass matrix depending on material, mesh and applied  model (coque_3D) for dynamic calculation (mode shapes and frequencies)
+
result = CALC_ELEM(...)
−
MACRO_MATR_ASSE(MODELE=modmod,
+
IMPR_RESU(...)  
−
                CHAM_MATER=material,
+
we can display the results in Salome (or any other post processing program).
−
                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
+
But have a look in the *.mess file first.
−
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
+
=='''Results'''==
−
salomres=PROJ_CHAMP(RESULTAT=MODES,
+
All pictures below show the z deformation due to a pressure of 1000 [Pa].<br/>
−
                    MODELE_1=modmod,      #project modes/results of model_1
+
Note that the mesh finess has been set to maximum ie. 1 instead of minimum 0 in the previous part.
−
                    MODELE_2=modinit,);  #to model_2
+
: [[Image:kw_dz_1a.png]] : [[Image:kw_dz_12.png]] <br/>
  
−
# ''store the results in a MED file''
+
left picture: shell thickness increases in positive y direction and undeformed mesh<br/>
−
# only the displacements 'DEPL' will be stored
+
right picture shows both increasing and decreasing shell thickness in positive y direction<br/>
−
# 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.
+
: [[Image:kw_dz_2.png]] <br/>
 +
detail of the maximum deformation
  
−
=='''Post processing with Salome'''==
+
Clearly the solution is asymmetric due to the increasing thickness and hence stiffness of the shells in y direction.
  
−
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.
+
=='''Input files - download'''==
−
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.
+
* [[Media: kw_shell_list.zip]]
 +
* [[Media: kw_variable_shell_thickness_v11p3.zip]] version 11.3
 +
* [[Media: Kw variable thickness pressure v22.zip]] version 11.3
  
−
A view of the first mode shape is given here:
 
−
: [[Image: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 ''unit''e 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:
 
Input files:
−
* Python script for defining the geometry (geom_plate1by1.py)
+
* Python script for defining the geometry and mesh (geom_mesh_shell.py), load this file by File --> load script (cntrl T in the object browser), refresh (F5) after running. Change on line 11 the variable ''meshsize'' between 0.0 (course) and 1.0 (fine)<br/>
−
* Python script for defining the mesh (mesh_plate1by1.py)
+
** Save the mesh by:
−
* python script of geometry and mesh (combination of those above, geom_mesh_plate1by1.py)
+
*** change to Mesh module
−
* ASTK file (shelldyn.astk, you need to edit the path to your requirements ...)
+
*** right clicking in the object browser by right clicking on the mesh name ''Mshell''; select ''Export to MED''
−
* command file (shelldyn.comm)
+
*** change mesh file name to mshell.med and click ''save''
−
output file:
+
* ASTK file (shellist.astk, you need to edit the path to your requirements ...). To start Code Aster calculation press ''Run''. Interactive enabled, Interactive-follow-up disabled, nodebug enabled.
−
* total mass returned by CA (masse.txt)
+
* command file (shell_static_d.comm), change sys.path.append('.....') for next file
−
download the files here:
+
* construct_shell_thickness Python file: place this file on your system and add the Python system path in your command file
−
* >>>[[Media:Fea_plate1by1.zip]]: error on TRIA6_7 in Code Aster10, please do not use anymore
+
** import sys #line 9
−
* >>>[[Media: kw_platedyn_ca10.zip]]: suitable for Aster10
+
** sys.path.append('/<your_path_to_construct_shell_thickness>')
  
−
=='''CODE ASTER version 10.X - TRIA6_7 and QUAD8_9'''==
+
The last file also has a revised function for applying the thickness to the shells. Also a function for applying a variable pressure on the shell elements is provided. It uses the same procedure as for the variable thickness.
−
Update in the second version:<br>
+
−
shellch=AFFE_CARA_ELEM(MODELE=modmod,
+
−
                      COQUE=_F(GROUP_MA='shell',
+
−
                                EPAIS=th,
+
−
                                COQUE_NCOU=1,),);
+
  
−
meshmod=CREA_MAILLAGE(MAILLAGE=meshinit,
+
: [[Image:kw_shell_shellplate_thlist.png]]
−
                      MODI_MAILLE=(_F(TOUT='OUI',OPTION='QUAD8_9',),),);
+
−
 
+
−
 
+
−
'''CREA_MAILLAGE'''<br>
+
−
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'''==
 
=='''References'''==
 
From Code Aster references:<br/>
 
From Code Aster references:<br/>
−
[[http://www.code-aster.org/V2/doc/default/man_u/u2/u2.02.01.pdf u2.02.01]] Notice d’utilisation des éléments plaques, coques et coques volumiques SHB<br/>
+
[[http://www.code-aster.org/V2/doc/default/en/man_u/u2/u2.02.01.pdf u2.02.01]] Notice d’utilisation des éléments plaques, coques et coques volumiques SHB, Note of use of the voluminal elements plates, shells,
−
[[http://www.code-aster.org/V2/doc/default/man_u/u4/u4.81.22.pdf u4.81.22]] Opérateur POST_ELEM<br/>
+
shells SHB, grids and Résumé<br/>
 +
[[http://www.code-aster.org/V2/doc/default/en/man_u/u4/u4.81.22.pdf u4.81.22]] Opérateur POST_ELEM<br/>
 
[[http://www.code-aster.org/DOCASTER/Man_U/U5/U50200b.pdf u5.02.00]] Manuel d’Utilisation - Fascicule U5.0- : Structure de données resultat Document: U5.02.00 - Présentation des tables<br/>
 
[[http://www.code-aster.org/DOCASTER/Man_U/U5/U50200b.pdf u5.02.00]] Manuel d’Utilisation - Fascicule U5.0- : Structure de données resultat Document: U5.02.00 - Présentation des tables<br/>
−
[[http://www.code-aster.org/V2/doc/default/man_u/u3/u3.12.03.pdf u3.12.03]] Modélisation COQUE_3D<br/>
+
[[http://www.code-aster.org/V2/doc/default/en/man_u/u3/u3.12.03.pdf u3.12.03]] Modélisation COQUE_3D<br/>
−
[[http://www.code-aster.org/V2/doc/default/man_r/r3/r3.07.04.pdf R3.07.04]] Eléments finis de coques volumiques<br/>
+
[[http://www.code-aster.org/V2/doc/default/en/man_r/r3/r3.07.04.pdf R3.07.04]] Eléments finis de coques volumiques, Finite elements of voluminal shells<br/>
−
[[http://www.code-aster.org/V2/doc/default/man_r/r3/r3.07.05.pdf R3.07.05]] Éléments de coques volumiques en non linéaire
+
[[http://www.code-aster.org/V2/doc/default/en/man_r/r3/r3.07.05.pdf R3.07.05]] Éléments de coques volumiques en non linéaire géométrique, Voluminal shell elements into nonlinear geometrical ''(i.e. machine translation France2Frenglish)''.
−
géométrique
+
<br/>
−
 
+
[[http://www.code-aster.org/V2/doc/default/en/man_u/u4/u4.42.01.pdf u4.42.01]] Opérateur AFFE_CARA_ELEM --8.3.1 EPAIS/EPAIS_F, Opérateur AFFE_CARA_ELEM<br/>
−
[http://www.caelinux.org/wiki/downloads/docs/PCarrico/CAELINUX_plasticite/CAELINUX_plasticite.html#SECTION000180000000000000000 Overview of units - Paul Carrico]
+
−
 
+
−
 
+
−
 
+
  
 +
Code Aster forum - discussion:<br/>
 +
[[http://www.code-aster.org/forum2/viewtopic.php?id=13692 Code Aster forum]] MAX number of GROUP_MAs & Variable shell thickness<br/>
 +
[[http://www.code-aster.org/forum2/viewtopic.php?pid=25803 CA forum - 2 -this is a circular reference]] coque à épaisseur variable<br/>
 +
[[http://www.code-aster.org/forum2/viewtopic.php?id=14692 More on Python and CA - by Bracchesimo]] Structure Optimization with Code_Aster<br/>
  
−
--[[User:Keeswouters|keeswouters]] 17:11, 24 April 2009 (CEST)
+
--[[User:Keeswouters|keeswouters]] november 2010
 +
-- update januari 2013

Latest revision as of 23:00, 3 December 2017

Varying thickness of shells - by Python

If you are already familiar to Code Aster, enjoy this. If you are new to Code Aster have a look here first please.
november 2010 - SalomeMeca 2010 (Code Aster 10.2) - Ubuntu Meerkat 10.10 64bits.
January 2013 - update to Code Aster version 11.3 - Salome 6.6.0 - Linux Mint 14 Nadia

For a similar contribution with variable pressue see: Contrib:KeesWouters/plate/variable_pressure

Static displacement of a shell under pressure

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.

Two meshes are used throughout this contribution: a course mesh (3 by 4 elements) and a fine mesh. In the geometry and mesh file (download as at the end) the meshsize can be set by the variable meshsize (number between 0.0 and 1.0). The course mesh I have been using for experimenting. The fine mesh to show more reasonable results.

General way of determination of the shell thickness

The geometry and mesh are 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

  • determine the nodes are defined by the connectivity matrix,
  • determine the coordinates of all the nodes and finally
  • determine the centre of gravity and
  • generate the thickness of the shell

by Python commands. 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. The minimum and maximum x, y and z values are:

  • x: [xmin, xmax] = [0.00, 1.00]
  • y: [ymin, ymax] = [0.00, 2.50]
  • z: [zmin, zmax] = [1.25, 2.35]

The thickness of the plate we make dependent on the y coordinate of the centre of the element as follows: th = Ycog*0.1+0.01 (edit this for your own definition here of course).
The 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 of the general calculation procedure see e.g. here. In the next part we concentrate on the extraction of the coordinates of the elements and the definition of the thickness of each element.

Kw varshell geom.jpg

Extraction of the coordinates

Among others, in Code Aster forum the procedure is discussed.
In Code Aster, after reading the initial mesh, the quad8 elements are converted to quad9 (coque_3d compatible) elements. This mesh is called meshmod (modified mesh) and is used as the basis for the further operations.

First we adapt the DEBUT() command since we need real Python commands: 

DEBUT(PAR_LOT='NON');

We import the Python partition module in the command file to read mesh quantities. The definition of the thickness of the shells is stored in a separate Python script/module/file (construct_shell_thickness) in the path defined below (change for your own needs [sys.path.append('/your/path/shellplate_thlist')]): 

from Utilitai.partition import *
import sys
sys.path.append('/cae/caexample/shell/shellplate_thlist')
from construct_shell_thickness import *

After reading the initial mesh meshinit with quad8 elements 

meshinit=LIRE_MAILLAGE(FORMAT='MED',INFO=1,);

and converting it to quad9, coque_3d elements 

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

we define the thickness of each element by calling 

info = 1   ## true/false - print additional information
ThicknessGroupShell = construct_shell_thickness(meshmod,info)

The Python script construct_shell_thickness does the work.
We have a look at this in detail. Standard, the characteristics of the shell are defined by the AFFE_CARA_ELEM command with, among others, the thickness parameter EPAIS: 

shellch=AFFE_CARA_ELEM(MODELE=modmod,
                      COQUE=_F(GROUP_MA='shell',
                               EPAIS=0.001,),);

Now remember that the parameters in the _F(...) operator can be replaced by a Python list, in this case: ellist = {'GROUP_MA': 'shell','EPAIS':0.001}. Or more general: 

 ellist = [{'MAILLE': 'Mxxx','EPAIS':thxxx}, 
           {'MAILLE': 'Myyy','EPAIS':thyyy},                       
             ::::
           {'MAILLE': 'Mzzz','EPAIS':thzzz}]

where 'Mxxx' is the element number and thxxx is the corresponding element thickness. Eg in this example the following list thicknessGroupShell is generated: 

 [{'EPAIS': 0.16625, 'MAILLE': 'M15'}, 
  {'EPAIS': 0.22875, 'MAILLE': 'M16'}, 
   ::::
  {'EPAIS': 0.22875, 'MAILLE': 'M25'}, 
  {'EPAIS': 0.10375, 'MAILLE': 'M26'}]

General - extraction of coordinates

In this example, the actual thickness of the shell depends on the y coordinate at the centre of the element. This is extracted as follows. First read the mesh: 

meshCA = MAIL_PY()
meshCA.FromAster(mesh)

In meshCA the available mesh quantities can be extracted (nodes, elements, connectivity, coordinates): 

 nonu     = meshCA.dime_maillage[0]    ## number of nodes
 elnu     = meshCA.dime_maillage[2]    ## number of elements
 ncoord   = meshCA.cn                  ## coordinates of the nodes
 connect  = meshCA.co                  ## connectivity of the element
 NodeList = list(meshCA.correspondance_noeuds)   ## says it all
 ElemList = list(meshCA.correspondance_mailles)
 ElemType = list(meshCA.tm)

The procedure is then as follows:

  • step through element list
    • check if element type is quad9/coque_3d/number 16
    • read connectivity of the element, ie all attached node numbers: connectivity = meshCA.co
    • step through all node of the element
      • determine nodes and its coordinates
      • compute average (y) coordinate
      • define thickness depending on (y) coordinate
      • add thickness and element to list {...}

In 'real' Python code: 

for ii in xrange(len(ElemList)):
  print ElemType[ii]
  if ElemType[ii]==16:  # select shell (coque_3d) elements
    #print connect[ii]
    shellcount +=1

ThicknessGroupShell=[]  # initialise variables
thickness = [None]*shellcount

shellcount = 0  # reset shellcount
for ii in xrange(len(ElemList)):
  print 'elementtype: ',ii,ElemType[ii]
  if ElemType[ii]==16:           # select shell elements coque_3d (14:quad8 but only 16:quad9 available)
    sumxyz = [0.0, 0.0, 0.0]
    ncount = len(connect[ii])    # should be 7 for tria7 and/or 9 quad9 for coque_3d (elementType 16)
    for jj in xrange(0,ncount):
      Gnode = connect[ii][jj]    # globalnode = connect[elemii,internal_nodejj]
      sumxyz += ncoord[Gnode]    # this might be written more compact --> for later

    sumxyz /=ncount
    # define your own code for thickness here 
    # thickness now depends on average y coordinate only
    thickness[shellcount] = sumxyz[1]*0.01+0.1  # make thickness dependending on y coordinate
    tempth = thickness[shellcount]
    shellcount +=1
    #print 'average xyz - centre point: ',sumxyz,ncoord[Gnode]  #centre point is same average nodes
    ThicknessGroupShell.append({'EPAIS':tempth,'MAILLE':'M%d'%(ii+1)});  ## add 1 for element number Mxx

If needed, vecteur, angle or NCOU can be defined as well.

Things that might go wrong

Of course, nothing will ever go wrong. Except code written in Python and used in conjunction with Code Aster. Then things definitely will go wrong. At least the first time by a unstructered programmer (like me?). First all kind of syntax errors occur. These are relatively easy to address and remove, because Python offers help by issuing nice error messages. Things are less trivial if you forget the +1 to inidicate the mesh element in the append command: 'MAILLE':'M%d'%(ii+1) because of the 1 offset numbering of Salome to Code Aster mesh. To see how the mesh is built up you can print the whole mesh FromAster:

meshCA = MAIL_PY()
meshCA.FromAster(mesh)
print meshCA,

Then you get eg: 

COOR_3D
N1         1.00000000000000E+00   2.50000000000000E+00   2.34560000000000E+00
.... 
SEG3
M1        N3       N17      N18     
....
QUAD9
M15       N10    N29    N30    N12    N36    N37    N38     N13    NS1  
M16       N7     N5     N32    N29    N8     N42    N43     N35    NS2     
....
GROUP_MA
Linex0
M5       M7       M10      M11     
FINSF
%
GROUP_MA
Linex1
M3       M9       M12      M14     
FINSF
%
GROUP_MA
shell
M15      M16      M17      M18      M19      M20      M21      M22     
M23      M24      M25      M26

This allows you to check the correct node numbers, elements, connectivity and groups defined in your mesh. Starting with a small mesh is advicable.
Also carefully check the code that generates the thickness of the shell. It is very easy to make a mistake: selecting the wrong index of the coordinate array is an obvious one.
For a small mesh (3 by 4 elements) the Python list is as follows: 

thicknessGroupShell:  
[{'EPAIS': 0.16625, 'MAILLE': 'M15'}, 
 {'EPAIS': 0.22875, 'MAILLE': 'M16'}, 
 {'EPAIS': 0.04125, 'MAILLE': 'M17'}, 
 {'EPAIS': 0.22875, 'MAILLE': 'M18'},  
 {'EPAIS': 0.16625, 'MAILLE': 'M19'}, 
 {'EPAIS': 0.10375, 'MAILLE': 'M20'}, 
 {'EPAIS': 0.10375, 'MAILLE': 'M21'}, 
 {'EPAIS': 0.04125, 'MAILLE': 'M22'},
 {'EPAIS': 0.04125, 'MAILLE': 'M23'}, 
 {'EPAIS': 0.16625, 'MAILLE': 'M24'}, 
 {'EPAIS': 0.22875, 'MAILLE': 'M25'}, 
 {'EPAIS': 0.10375, 'MAILLE': 'M26'}]


Remarks
Note that you could select the midplane node NSxxx just as well to make the thickness dependending on (NSxxx is added to the quad8 element by the modi_maillage command). Then there is no need to compute the average values of the coordinates of the 8 nodes on the egde.


Kw varshell 3x4.jpg

Continue with the commands

Now we follow the standard procedure for calculating the static displacement of the shell under pressure. Normally we define the the thickness of the shell by: 

th = 0.02
shellch=AFFE_CARA_ELEM(MODELE=modelc,
                      COQUE=_F(GROUP_MA='shell',
                               EPAIS=th,),);
                               #VECTEUR=(1.0,0.0,0.0,),
                               #COQUE_NCOU=1,),);

This will be replaced by: 

shellch=AFFE_CARA_ELEM(MODELE=modelc,COQUE=(ThicknessGroupShell),);

Very compact, I would say. Unfortunately the whole process defined in construct_shell_thickness hides the effort. The MECA_STATIQUE command carries out the static analysis: 

result=MECA_STATIQUE(MODELE=modelc,
                    CHAM_MATER=material,
                    CARA_ELEM=shellch,
                    EXCIT=_F(CHARGE=clamped,),);

where of course the model modelc, the material properties, the shell characteristics shellch (oeps) and the charge clamped -being boundary conditions in this case- have been defined in advance. After calculating the element results and printing the results by 

result = CALC_ELEM(...)
IMPR_RESU(...)

we can display the results in Salome (or any other post processing program).

But have a look in the *.mess file first.

Results

All pictures below show the z deformation due to a pressure of 1000 [Pa].
Note that the mesh finess has been set to maximum ie. 1 instead of minimum 0 in the previous part.

Kw dz 1a.png : Kw dz 12.png

left picture: shell thickness increases in positive y direction and undeformed mesh
right picture shows both increasing and decreasing shell thickness in positive y direction

Kw dz 2.png

detail of the maximum deformation

Clearly the solution is asymmetric due to the increasing thickness and hence stiffness of the shells in y direction.

Input files - download

Input files:

  • Python script for defining the geometry and mesh (geom_mesh_shell.py), load this file by File --> load script (cntrl T in the object browser), refresh (F5) after running. Change on line 11 the variable meshsize between 0.0 (course) and 1.0 (fine)
    • Save the mesh by:
      • change to Mesh module
      • right clicking in the object browser by right clicking on the mesh name Mshell; select Export to MED
      • change mesh file name to mshell.med and click save
  • ASTK file (shellist.astk, you need to edit the path to your requirements ...). To start Code Aster calculation press Run. Interactive enabled, Interactive-follow-up disabled, nodebug enabled.
  • command file (shell_static_d.comm), change sys.path.append('.....') for next file
  • construct_shell_thickness Python file: place this file on your system and add the Python system path in your command file
    • import sys #line 9
    • sys.path.append('/<your_path_to_construct_shell_thickness>')

The last file also has a revised function for applying the thickness to the shells. Also a function for applying a variable pressure on the shell elements is provided. It uses the same procedure as for the variable thickness.

Kw shell shellplate thlist.png

References

From Code Aster references:
[u2.02.01] Notice d’utilisation des éléments plaques, coques et coques volumiques SHB, Note of use of the voluminal elements plates, shells, shells SHB, grids and Résumé
[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, Finite elements of voluminal shells
[R3.07.05] Éléments de coques volumiques en non linéaire géométrique, Voluminal shell elements into nonlinear geometrical (i.e. machine translation France2Frenglish).
[u4.42.01] Opérateur AFFE_CARA_ELEM --8.3.1 EPAIS/EPAIS_F, Opérateur AFFE_CARA_ELEM

Code Aster forum - discussion:
[Code Aster forum] MAX number of GROUP_MAs & Variable shell thickness
[CA forum - 2 -this is a circular reference] coque à épaisseur variable
[More on Python and CA - by Bracchesimo] Structure Optimization with Code_Aster

--keeswouters november 2010 -- update januari 2013

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