SImulation par ÉLéments finiX Code SILEX

It seems you are using an outdated browser that does not support canvas :-(

Introduction

• SILEX is a finite element code written in Python language, eventually with a Fortran part in order to speed up the computations.
  
  
• The Python language is used to define parameters, to read the mesh, to solve the system, to write the results.
  
  
• The Fortran language is eventually used for elemental computations as well as to build the stiffness matrix.
  
  
• Le open source software Gmsh is used to create the meshes as well as to show the results.
  
  
• The only free routines available on-line concern the 4-node tetrahedral element and the 3-node triangle element in the case of linear static analysis. They can be adapt to other elements.
  
  
• The following proposed applications allow to understand the code, in order to perform other computations for other mechanical systems. Later on, the user can develop new elements or new method, and thus extend the possibilities of the code.
  
  
• The following course document available on-line ( here ) gives a programming introduction.

Download the SILEX routines:

• SILEX ROUTINES

  
  
  

• Python 3.4 and an ide, Idle-python for instance
  
• eventually Fortran and gfortran
  
• gmsh

• Anaconda for Windows is the best way to do it: it contains all the necessary Python librairies and text editor (Spider):
  Download Anaconda
  
  
• Gmsh:
  Download Gmsh
  
  
• The text editor Notepad++ is perfect with Windows to edit text files:
  Download Notepad++
  
  

Computation of a piston: first steps with SILEX

Create the mesh with Gmsh

• With the help of the mouse, we can easily turn the solid (left button), move it (right button) or zoom it (central button).
  
  
• Click on Mesh; Define; Element size at points: select all the points with the left click of the mouse and with holding the Ctrl key on: points should appear with a red color. Confirm the command by pressing the e key.
  
  
• The physical groups can be used to give names (numbers) to elements and nodes on which we want to apply conditions: Geometry; physical group; add; surface: select the 2 following: then press the e key
  
  Surface S1:
  
  
  
  
  
  
  
• Do the same for the surface S2:
  
  
  
  
  
  
  Do the same for the surface S3:
  
  
  
  
  
  
  Then, do the same for surface S4 (select the 1/4 of the sphere and the 1/4 of the "horizontal" ring):
  
  
  
  
  
• Define the volume:
  
  geometry/physical group/add/volume
  
  select the small yellow ball, then confirm by pressing the e key
  
  
• In the gedit (or Notepad++) window of the file piston-tet4.geo, click on reload.
  
  
• Modify the numbers of the "physical groups" in this text file in order to have the correct numbers: 1, 2, 3 and 4 for the surfaces, and 5 for the volume
  
  
• Modify at the end of the line Characteristic Length the value 0.1 into h
  
  
• Save the file in the gedit (or Notepad++) window
  
  
• Go back in the gmsh window, and then, click on reload
  
  
• Click on Mesh, 2D, the surface mesh should appear (otherwise go in tools/options/mesh/visibility in order to change the entities you want to plot).
  
  
• Click on 3D, the volume mesh should appear.
  
  
• In the File drop-down menu, choose Save Mesh, the mesh is saved in the file piston-tet4.msh.
  
  

With Ubuntu: use of the Fortran library for Python (it is not necessary for the full Python version)

Main Python program Main-Piston.py

• Part Import libraries:
  
  This first part indicates to the Python the specific used libraries in this program.
  
  
  
• By moving the # comment symbol, choose the elemental library:
  
  
  • For a full Python version (typically for Windows, at home):
    
    #import silex_lib_tet4_fortran as silex_lib_elt
    import silex_lib_tet4_python as silex_lib_elt
    
    
  • For a compiled Fortran library (typically for Ubuntu, at the university):
    
    import silex_lib_tet4_fortran as silex_lib_elt
    #import silex_lib_tet4_python as silex_lib_elt
  
  
  Note that you save you work one a USB key, and change the # comment symbol according to the used computer.
  
  
  
• Part USER PART: Import mesh, boundary conditions and material:
  
  
  • The variable MeshFileName gives the name of the mesh file coming from Gmsh (with no .msh extension).
    
    
  • The variable ResultsFileName gives the name of the result file (Gmsh format) in which the selected results are written (with no .msh extension).
    
    
  • The variable eltype gives the element type used in the program: 4 for 4-node tetrahedral. Other types of elements are possible, see the silex_lib_gmsh.py library to have the full list.
    
    
  • The variable ndim gives the dimension of the geometry: 2 for 2D, 3 for 3D.
    
    
  • The variable flag_write_fields gives the type of results that are saved in the result file:
    • 1 for all fields,
    • 0 for a "short" output.
    • Other values can be defined by the user at the end of this file (Main-Piston-tet4.py).
    
    
    
  • The routine ReadGmshNodes of the library silex_lib_gmsh is used to import node coordinates from the mesh file in the table nodes.
    
    Use:
    
    • nodes=silex_lib_gmsh.ReadGmshNodes('mymesh.msh', 2 ou 3)
    
    
    
  • The routine ReadGmshElements of the library silex_lib_gmsh is used to import the connectivity table from the mesh file in the table elements, the corresponding node numbers are in the table Idnodes.
    
    Use:
    
    • elements,Idnodes=silex_lib_gmsh.ReadGmshElements('mymesh.msh',element_type,physical_group)
    
    
    
  • As soon as the volume mesh is imported, we use the routine ReadGmshElements in order to import the 3-node triangle surface meshes (type of element 2) of the important surfaces of the computation: elementsS1 to elementsS4 and IdnodeS1 to IdnodeS4.
    
    
  • In order to check if the meshes are correctly read and imported, we can write them in Gmsh format files with the help of the WriteResults routine of from library silex_lib_gmsh.
    
    Use:
    
    • silex_lib_gmsh.WriteResults('checkmesh.msh',nodes,elements,element_type)
    
    
    
  • The variable Young gives the Young's modulus, here 200000 MPa for steel.
    
    
  • The variable nu gives the Poisson's coefficient, here 0.3 for steel.
    
    
  • The fixed nodes in x direction are given in the table IdNodesFixed_x (here IdnodeS3)
    
    
  • The fixed nodes in y direction are given in the table IdNodesFixed_y (here IdnodeS1)
    
    
  • The fixed nodes in z direction are given in the table IdNodesFixed_z (here IdnodeS2).
    
    
  • We use in the program the following numbering:
    • Local numbering:
      
      

      Local numbering of the dofs in an element

      Name of the unknown	Ux1	Uy1	Uz1	Ux2	Uy2	Uz2	Ux3	Uy3	Uz3	Ux4	Uy4	Uz4
      Python numbering	0	1	2	3	4	5	6	7	8	9	10	11
      Fortran numbering	1	2	3	4	5	6	7	8	9	10	11	12

      
      
      
    • Global numbering:
      

      Global numbering of the dofs

      Name of the unknown	Ux1	Uy1	Uz1	Ux2	Uy2	Uz2	.............	Uxi	Uyi	Uzi	.............	Uxn	Uyn	Uzn
      Python numbering	0	1	2	3	4	5	.............	(i-1)*3	(i-1)*3+1	(i-1)*3+2	.............	(n-1)*3	(n-1)*3+1	(n-1)*3+2
      Fortran numbering	1	2	3	4	5	6	.............	(i-1)*3+1	(i-1)*3+2	(i-1)*3+3	.............	(n-1)*3+1	(n-1)*3+2	(n-1)*3+3

      
      
      
  • The upper face of the piston is loaded under a pressure of 10MPa (100bar). The triangle elements from this surface are in the table elementsS4.
    The routine forceonsurface of the library silex_lib_tet4 (renamed into silex_lib_elt) allow to compute the forces on the nodes for a given surface force with the norm press in the direction direction:
    
    Use:
    
    • F = silex_lib_elt.forceonsurface(nodes,surface_elements_where_the_load_is_applied,press,direction)
    
    
    If the given direction is [0.0,0.0,0.0], then the local unit normal to the surface is used, it is the case for a surface pressure load. In this case, the normal is computed from the local numbering in the elements, so check if the load is correctly applied. If it is in the wrong direction, just change the sign of the pressure value.
    
    
• EXPERT part
  
  
  • nnodes : number of nodes.
    
    
  • ndof : number of degrees of freedom (dofs).
    
    
  • nelem : number of elements.
    
    
  • Fixed_Dofs : fixed dofs.
    
    
  • SolvedDofs : free dofs.
    
    
  • Q : nodal displacements.
    
    
  • Compute the stiffness matrix and sparse format assembly:
    
    
    • Ik,Jk,Vk=silex_lib_elt.stiffnessmatrix(nodes,elements,[Young,nu])
      K=scipy.sparse.csc_matrix( (Vk,(Ik,Jk)), shape=(ndof,ndof) ,dtype=float)
    
    
    
  • Solve the system (note that we assume the the imposed displacements are zero, otherwise, it has to be modified):
    
    
    • Q[SolvedDofs] = scipy.sparse.linalg.spsolve(K[SolvedDofs,:][:,SolvedDofs],F[SolvedDofs])
    
    
    
  • Compute:
    • the stresses (in the elements as well as the nodal smooth stress)
    • the strain (in the elements as well as the nodal smooth stress)
    • the elemental contributions to the global error
    • the global energy error (ZZ1)
    
    
    
    • SigmaElem, SigmaNodes, EpsilonElem, EpsilonNodes, ErrorElem, ErrorGlobal = silex_lib_elt.compute_stress_strain_error(nodes,elements,[Young,nu],Q)
    
    
    
  • disp : table of the displacements of the nodes written using 3 columns (or 2 for a 2D plane computation) for the Gmsh format output.
    
    
  • load : given and imposed forces on the nodes written using 3 columns (or 2 for a 2D plane computation) for the Gmsh format output.
    
    
  • fields_to_write : Python list of the fields to save in the result file. See the routine WriteResults of the library silex_lib_gmsh for more details.
    
    

Results "Piston": displacements for h=0.8

Results "Piston": Von Mises stress in the elements for h=0.8

Results "Piston": smooth Von Mises stress at the nodes for h=0.8

Results "Piston": elemental contribution to the global error for h=0.8

A new possible analysis : Landing gear fork