SImulation par ÉLéments finiX Code SILEX

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

Introduction

• SILEX est un code de calcul par éléments finis écrit en Python avec éventuellement une partie en Fortran pour accélérer les calculs.
  
  
• Le langage Python est utilisé pour définir les paramètres du calcul, lire le maillage, résoudre le système, écrire les résultats.
  
  
• Le langage Fortran est éventuellement utilisé pour effectuer les calculs au niveau local et construire la matrice de rigidité.
  
  
• Le logiciel libre Gmsh est utilisé pour faire les maillages et visualiser les résultats.
  
  
• Seules les routines pour le tétraèdre à 4 nœuds et pour le triangle à 3 nœuds dans le cas de statique linéaire sont disponibles pour l'instant en téléchargement libre.
  
  
• Les applications proposées dans la suite permettent de prendre en main le code afin de pouvoir faire d'autre calculs dans un premier temps, puis de développer d'autres routines dans un second temps.
  
  
• Le document de cours présent sur ce site web ( ici ) propose une initiation à la programmation. Les programmes du documents sont téléchargeables ici .

Téléchargement des routines du code SILEX et librairies nécessaires:

• ROUTINES SILEX

  
  
  

• Python 3.4 et un ide, Idle-python par exemple
  
• Éventuellement Fortran avec gfortran
  
• gmsh

• Anaconda pour Windows contient tous les paquets nécessaires (avec l'ide Spider):
  Téléchargement d'Anaconda (Python 3)
  
  
• Gmsh:
  Téléchargement de Gmsh
  
  
• En complément, l'éditeur de texte Notepad++ est parfait sous windows pour éditer les fichiers textes:
  Téléchargement de Notepad++
  
  

Présentation du calcul d'un piston pour la prise en main de SILEX

Maillage avec Gmsh

• A l'aide de la souris, on peut de façon très simple faire tourner le modèle (clic gauche), le déplacer (clic droit) ou faire des zooms (bouton central).
  
  
• Cliquer sur Mesh; Define; Element size at points: sélectionner tous les points avec le bouton gauche de la souris en laissant la touche Ctrl enfoncée: les points doivent apparaitrent en rouge. Valider en tapant e.
  
  
• Les physical group affectent des noms aux éléments et aux nœuds pour pouvoir ensuite leur appliquer des conditions: Geometry; physical group; add; surface: sélectionner les 2 surfaces suivantes: puis taper e
  
  Surface S1:
  
  
  
  
  
  
  
• Faire de même avec pour la surface S2:
  
  
  
  
  
  
  Puis pour la surface S3:
  
  
  
  
  
  
  Enfin pour la surface S4 (cliquer sur le 1/4 de sphère et sur le 1/4 d'anneau "horizontal"):
  
  
  
  
  
• Définir le volume:
  
  geometry/physical group/add/volume
  
  sélectionner la boule jaune, puis taper e
  
  
• Taper reload dans la fenêtre gedit (ou Notepad++) du fichier Piston.geo
  
  
• Modifier les numéros des "physical groups" du fichier texte Piston.geo pour avoir les numéros corrects: 1, 2, 3 et 4 pour les surfaces, et 5 pour le volume
  
  
• Modifier aussi à la fin de Characteristic Length le 0.1 en h
  
  
• Sauver le fichier dans gedit (ou Notepad++)
  
  
• Retourner dans la fenêtre gmsh puis taper reload
  
  
• Cliquer sur Mesh, 2D, le maillage surfacique doit apparaitre (sinon aller dans tools/options/mesh/visibility pour changer les entités affichées).
  
  
• Cliquer sur 3D, le maillage volumique doit apparaitre.
  
  
• Version 2 de Gmsh : Dans le menu déroulant File, choisir Save Mesh, le maillage est ainsi sauvegardé dans le fichier piston-tet4.msh.
  
  
• Version 4 de Gmsh : Dans le menu déroulant File, choisir Export, puis le format Mesh - Gmsh (*.msh), et enfin Version 2 ASCII dans le menu déroulant. Le maillage est ainsi sauvegardé dans le fichier piston-tet4.msh.
  
  

Sous Ubuntu, compilation de la librairie fortran pour Python (non nécessaire pour la version en Python seul)

Explications du programme en Python Main-Piston.py

• Partie Import libraries:
  
  Cette première partie indique à python les librairies spécifiques que l'on va utiliser.
  
  
  
• En déplaçant le # du commentaire, choisir la librairie de calcul au niveau élémentaire:
  
  
  • Pour une version Python (typiquement pour Windows):
    
    #import silex_lib_tet4_fortran as silex_lib_elt
    import silex_lib_tet4_python as silex_lib_elt
    
    
  • Pour une version avec une librairie Fortran compilée:
    
    import silex_lib_tet4_fortran as silex_lib_elt
    #import silex_lib_tet4_python as silex_lib_elt
  
  
  
• Partie USER PART: Import mesh, boundary conditions and material:
  
  
  • La variable MeshFileName indique le nom du fichier maillage de provenant de gmsh (sans l'extension .msh).
    
    
  • La variable ResultsFileName indique le nom du fichier dans lequel les résultats sont écrits (sans l'extension .msh).
    
    
  • La variable eltype indique le type d'éléments que l'on utilise: 4 pour les tétraèdres à 4 nœuds. Les autres types d'éléments sont détaillés dans la librairie silex_lib_gmsh.py.
    
    
  • La variable ndim indique la dimension de la géométrie: 2 pour 2D, 3 pour 3D.
    
    
  • La variable flag_write_fields indique le type de résultats que l'on souhaite écrire dans le fichier résultat: 1 pour tous les champs, 0 pour une sortie courte. D'autres valeurs peuvent être définies à la fin du fichier Main-Piston-tet4.py.
    
    
  • La routine ReadGmshNodes de la librairie silex_lib_gmsh permet d'importer les coordonnées des nœuds depuis le fichier de maillage dans le tableau nodes.
    
    Utilisation:
    
    • nodes=silex_lib_gmsh.ReadGmshNodes('monmaillage.msh', 2 ou 3)
    
    
    
  • La routine ReadGmshElements de la librairie silex_lib_gmsh permet d'importer la table de connectivité des éléments depuis le fichier de maillage dans le tableau elements, les numéros des nœuds correspondants sont stockés dans le tableau Idnodes.
    
    Utilisation:
    
    • elements,Idnodes=silex_lib_gmsh.ReadGmshElements('monmaillage.msh',element_type,physical_group)
    
    
    
  • Une fois le maillage volume importé, on utilise la routine ReadGmshElements pour importer les maillages surfaciques de triangles à 3 nœuds (type d'élément 2) des surfaces fonctionnelles du calcul: elementsS1 à elementsS4 et IdnodeS1 à IdnodeS4.
    
    
  • Afin de vérifier si les maillages sont bien importés, on peut les écrire dans un fichier au format gmsh à l'aide de la routine WriteResults de la librairie silex_lib_gmsh.
    
    Utilisation:
    
    • silex_lib_gmsh.WriteResults('maillage_lu.msh',nodes,elements,element_type)
    
    
    
  • La variable Young indique le module d'Young, ici 200000 MPa pour de l'acier.
    
    
  • La variable nu indique le coefficient de Poisson, ici 0.3 pour de l'acier.
    
    
  • Les nœuds bloqués dans la direction x sont indiqués dans le tableau IdNodesFixed_x (ici IdnodeS3), ceux bloqués dans la direction y dans le tableau IdNodesFixed_y (ici IdnodeS1), enfin ceux bloqués dans la direction z dans le tableau IdNodesFixed_z (ici IdnodeS2).
    
    
  • On utilise dans le programme les numérotations suivantes:
    • Numérotation locale:
      
      

      Numérotation des ddls au niveau local

      Noms des inconnues	Ux1	Uy1	Uz1	Ux2	Uy2	Uz2	Ux3	Uy3	Uz3	Ux4	Uy4	Uz4
      Indexage Fortran	1	2	3	4	5	6	7	8	9	10	11	12

      
      
      
    • Numérotation globale:
      

      Numérotation des ddls au niveau global

      Noms des inconnues	Ux1	Uy1	Uz1	Ux2	Uy2	Uz2	.............	Uxi	Uyi	Uzi	.............	Uxn	Uyn	Uzn
      Indexage Python	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
      Indexage Fortran	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

      
      
      
  • La face supérieure est soumise à une pression de 10MPa (100bar). Les éléments triangles constituants cette surface sont dans le tableau elementsS4.
    La routine forceonsurface de la librairie silex_lib_tet4 (renomée en silex_lib_elt) permet de calculer les forces aux nœuds dues à une force surfacique de norme press dans la direction direction:
    
    Utilisation:
    
    • F = silex_lib_elt.forceonsurface(nodes,surfacic_elements_where_the_load_is_apllied,press,direction)
    
    
    Si la direction vaut [0.0,0.0,0.0], alors c'est la normale locale à la surface qui est utilisée, c'est le cas pour une force surfacique due à une pression. Dans ce cas la direction de la normale étant calculée à partir de la numérotation locale des nœuds, elle peut être en sens opposé, il suffit de changer le signe de press pour obtenir une force dans la bonne direction.
    
    
• Partie EXPERT
  
  
  • nnodes : nombre de nœuds.
    
    
  • ndof : nombre de degrés de liberté.
    
    
  • nelem : nombre d'éléments.
    
    
  • Fixed_Dofs : ddl bloqués.
    
    
  • SolvedDofs : ddl libres à calculer.
    
    
  • Q : déplacements aux nœuds.
    
    
  • Calcul de la matrice de rigidité et assemblage au format sparse :
    
    
    • Ik,Jk,Vk=silex_lib_elt.stiffnessmatrix(nodes,elements,[Young,nu])
      K=scipy.sparse.csc_matrix( (Vk,(Ik,Jk)), shape=(ndof,ndof) ,dtype=float)
    
    
    
  • Résolution du système en déplacement (on suppose que les déplacements imposés sont nuls, à modifier dans le cas contraire):
    
    
    • Q[SolvedDofs] = scipy.sparse.linalg.spsolve(K[SolvedDofs,:][:,SolvedDofs],F[SolvedDofs])
    
    
    
  • Calcul des contraintes (dans les éléments et lissées aux nœuds), des déformations (dans les éléments et lissées aux nœuds), des contributions des éléments à l'erreur, et enfin de l'erreur en énergie ZZ1 :
    
    
    • SigmaElem, SigmaNodes, EpsilonElem, EpsilonNodes, ErrorElem, ErrorGlobal = silex_lib_elt.compute_stress_strain_error(nodes,elements,[Young,nu],Q)
    
    
    
  • disp : déplacements des nœuds écrit sur 3 colonnes (ou 2 en 2D) pour la sortie au format gmsh.
    
    
  • load : efforts imposés aux nœuds écrit sur 3 colonnes (ou 2 en 2D) pour la sortie au format gmsh.
    
    
  • fields_to_write : liste python des champs à écrire dans le fichier résultat. Voir routine WriteResults de la librairie silex_lib_gmsh pour plus de détails.
    
    

Résultats "Piston": déplacements pour h=0.8

Résultats "Piston": Contraintes dans les éléments pour h=0.8

Résultats "Piston": Contraintes lissées pour h=0.8

Résultats "Piston": erreurs dans les éléments pour h=0.8

Calcul d'une autre pièce : fourche de train d'atterrissage