Sciences-industrielles.com

Cours, exercices et corrections en SI

Python scientifique pour les éléments finis

Logiciels 26 décembre 2013, par Steven Masfaraud

Ce cours présente les outils de python utilisables pour les éléments finis (matrices, vecteurs) pour aboutir à un code élément fini pour une poutre en traction. Sont présentés également des méthodes pour visualiser et vérifier les résultats produits.

 1- Installation de Python

1.1- Généralités sur Python

Python est un langage open-source multi-plateformes. Son intérêt réside notemment dans les nombreuses bibliothèques qu’il propose. Dans le domaine du calcul scientifique nous montrerons qu’il remplace avantageusement matlab.
Nous utiliserons Spyder,un éditeur de texte libre lui aussi, qui nous permettra d’avoir une interface interractive qui ne dépaysera pas les habitués de Matlab.
Il existe bien sur d’autres éditeurs de texte vous pouvez utiliser si vous préférez. Les installations proposées après installent Python, numpy et scipy. Si vous voulez utiliser un autre éditeur, assurez vous d’avoir ces paquets installés.

1.2- Sous Linux

Pour les utilisateurs de Debian ou d’Ubuntu, cette ligne de commande suffit à installer. Connectez vous en superutilisateur pour l’executer.

aptitude install spyder

Pour utilisateurs d’autres distributions linux, se référer à votre gestionnaire de paquets qui doit proposer Spyder.

1.3- Sous Windows

Pour installer Spyder, Python et les bibliothèques que nous allons utiliser, le plus simple et d’installer une distribution scientifique qui comporte tous cela.
Je vous propose winPython :
http://sourceforge.net/projects/winpython/files/WinPython_2.7/2.7.6.0/
Sur la page de téléchargements, choisir votre fichier d’installation en fonction de l’architecture de votre ordinateur. En général la version 64 bits convient. La procédure d’installation est classique.
++++

 2- Premiers pas sur Python, utilisation de Spyder

Spyder est un environnement de développement sous Python qui permet sur une même fenêtre d’éditer un fichier Python et de l’executer.
Deux sous fenêtres sont donc absolument indispensables : l’éditeur et la console.

Je vous conseille de maximiser en largeur la taille de l’éditeur car c’est cette fenetre dans laquelle vous passerez le plus de temps, et ensuite de maximiser la taille de la console notamment en hauteur.
Les autres fenetres ne sont pas indispensables, vous pouvez les fermer pour l’instant pour économiser de la place pour l’éditeur et la console.

L’interêt de Spyder et d’éditer un fichier, puis en appuyant sur F5 il est executé dans la console. On peut alors voir le résultat de nos modifications en direct.
Créez un nouveau fichier. Enregistrez le où vous voulez. Spyder vous à généré un code minimaliste qui ressemble à celui-ci :

  1. # -*- coding: utf-8 -*-
  2. """
  3. Created on Tue Dec 24 17:07:37 2013
  4.  
  5. @author: steven
  6. """

Télécharger

La première ligne précise l’encodage. Cela permet notamment d’afficher les caractéres spéciaux proprement. Dans mon cas, sous linux c’est de l’utf8.
Les lignes suivantes, entourées de """ sont une description de votre code, vous pouvez renseigner le but de votre script, l’auteur, etc...

Ajoutez la ligne suivante à la fin de votre fichier :

  1. print "Ca marche!"

Appuyez ensuite sur F5. Si une boîte de dialogue s’ouvre pour vous demander les options d’execution, cochez "Executer dans l’interpréteur Python" actif puis validez.
Vous devriez voir dans la console le message "Ca marche !" être affiché. Vous avez exécuté votre premier script Python !

++++

 3- Résolution éléments finis pour une poutre en traction

3.1- Problème, résolution éléments finis

On s’intéresse au comportement d’une poutre console (encastrée à une extrémitée) en traction de longueur L, de section S, de module de Young E soumise à une force F en bout de poutre et à une force linéique sur l’ensemble de la poutre de valeur f.

3.2- Formulation du problème

La théorie de la résistance des matériaux nous fournit les deux équations suivantes :

\left\{\begin{array}
\frac{dN}{dx}=-f \\
E\times S\times \frac{du}{dx}=N \\
\end{array}\right

avec :

  • N l’effort normal dans la poutre
  • u le déplacement longitudinal dû à la traction.

En combinant ces deux équations on obtient :

E\times S \times \frac{d^2u}{dx^2}=-f

Cette équation différentiell peut être facilement résolue analytiquement (à la main). Dans le cadre général des problèmes d’élasticité 2D ou 3D se n’est pas possible. D’oû la méthode d’approximation que sont les éléments finis.

3.3- Rappels sur la résolution éléments finis

En multipliant l’équation précédente par un champ test v et en l’intégrant sur la longeur de la poutre, on obtient :

\int_0^L E\times S \times \frac{d^2u}{dx^2} \times v =-\int_0^L f \times v

En intégrant par partie on obtient :

 [E\times S \times \frac{du}{dx} v]_0^L-\int_0^L E\times S \times \frac{du}{dx} \times \frac{dv}{dx} =-\int_0^L f \times v

On impose à v d’être nul en L et 0, c’est à dire sur les bords du domaine. D’où :

 \int_0^L E\times S \times \frac{du}{dx} \times \frac{dv}{dx} =\int_0^L f \times v

C’est alors que l’on injecte l’approximation éléments finis :

 u=\sum_{i=1}^{nn} \varphi_i u_i u=\pmatrix{
\varphi_1\cr
\varphi_2\cr
\vdots\cr
\varphi_{nn}\cr} \pmatrix{
u_1 & u_2 & \dots & u_{nn}}

avec nn le nombre de noeuds.

On choisit pour v :

 v=\sum_{i=1}^{nn} \varphi_i v=\pmatrix{
\varphi_1\cr
\varphi_2\cr
\vdots\cr
\varphi_{nn}\cr}

Ce qui conduit à :

 \underbrace{\int_0^L E S \pmatrix{
\frac{d\varphi_1}{dx}\cr
\frac{d\varphi_2}{dx}\cr
\vdots\cr
\frac{d\varphi_{nn}}{dx}\cr}
\pmatrix{
\frac{d\varphi_1}{dx} &
\frac{d\varphi_2}{dx} &
\vdots &
\frac{d\varphi_{nn}}{dx} & }}_{\mathbb{K}}
\pmatrix{
u_1 \cr u_2 \cr \vdots \cr u_{nn}\cr}
 =\underbrace{\int_0^L f \pmatrix{
\varphi_1\cr
\varphi_2\cr
\vdots\cr
\varphi_{nn}\cr}}_{F}

On peut montrer que la construction de ces opérateurs peut se décomposer en l’addition de contribution des éléments à la rigidité du système et au vecteur force.

On a donc affaire à la résolution d’un système matriciel :

 \mathbb{K}u=F

Voilà pour la théorie, voyons comment coder cela en pratique !

+++++

 4- Implémentation en Python

Dans la suite, nous écrirons notre programme dans un fichier Python. Pour cela, ouvrez un nouveau fichier et enregistrez le où vous le souhaitez.

4.1- Chargement des paquets

Comme nous l’avons vu, Python fonctionne avec des bibliothèques d’outils regroupés en paquets (packages). Le langage en lui meme définit des éléments de base comme les boucles, les conditions, la gestion des variables et certaines structures de données comme les listes.
Pour résoudre un problème éléments finis, nous aurons besoin de manipuler des vecteurs et des matrices, et de tracer des graphiques pour visualiser les résultats. Pour cela nous utiliserons des paquets qui sont fournis avec les distributions qui ont été présentées précédemment.

Nous allons utiliser :

  • numpy : un paquet définissant les matrices et vecteurs et permettant de les manipuler
  • linalg, sous-paquet de numpy
  • pyplot, sous-paquet de matplotlib, comme son nom l’indique une bibliothèque de tracé pour les maths.

Certains paquets peuvent en contenir d’autres, ce qui permet dans le cas de paquets volumineux de ne pas charger trop de fonctions dont on ne va pas se servir.

Pour importer un paquet on utilise le mot-clé import.
Pour importer la paquet numpy on doit écrire :

  1. import numpy

Tout simplement. Nous pouvons alors utiliser les outils de numpy par exemple la fonction identity en tapant :

  1. numpy.identity(3)

Ecrivez ces lignes à la fin de votre script et executez le. La console vous affiche le résultat :

array([[ 1.,  0.,  0.],
      [ 0.,  1.,  0.],
      [ 0.,  0.,  1.]])

La fonction numpy.identity crée une matrice identitée (1 sur la diagonale, zéros ailleurs) de la taille passée en paramètre.

Pour appeller la fonction identity inclue dans le paquet numpy il faut taper numpy.identity
Comme nous risquons d’appeler plusieurs fonctions de ce paquet, il faudra à chaque fois les préfixer par "numpy.". Pour ne pas avoir à faire cela, une solution consiste à modifier l’import du paquet en tapant :

  1. from numpy import *

On peut alors taper directement identity(3). Mais cette méthode pose un problème pour les grands projets où sont utilisés beaucoup de paquets. Deux paquets peuvent avoir chacun une fonction avec le meme nom pour des utilisations radicalement différentes. On est alors plus trop sûr de quelle fonction va etre appellée.
Une solution à mi-chemin est d’importer le paquet sous un nom raccourci :

  1. import numpy as np

Il faut alors taper np.identity(3). Cette solution garde l’avantage de savoir quel paquet on utilise sans pour autant alourdir trop le code. Nous garderons cette écriture pour la suite du projet.

Les autres paquets à inclure sont des sous-paquets. Pour cela :

  1. from numpy import linalg as lg
  2. from matplotlib import pyplot as plt

Télécharger

Pour résumer, voici à quoi doit ressembler votre script. Remarquez la ligne de commentaire, précédée d’un # pour commenter le code et le rendre plus sympa.

  1. # -*- coding: utf-8 -*-
  2. """
  3. @author: steven
  4. """
  5.  
  6. # Import des packages
  7. import numpy as np
  8. from numpy import linalg as lg
  9. from matplotlib import pyplot as plt

Télécharger

4.2- Mise en données

L’étape précédente était purement informatique, nous nous attaquons maintenant à la mise en données du problème avant sa résolution.
Dans chacune des étapes de la résolution, nous aurons besoin de connaitre certaines grandeurs comme les paramètres matériau ou les dimensions de la poutre ou encore des paramètres de discrétisation.

  1. # Parametres
  2. E=200e9 #Module de young
  3. S=0.01 # Section de la poutre
  4. L=1 # Longueur de la poutre
  5. Fp=1000 # Intensite de la force ponctuelle en Newtons
  6. ne=5 # Nombre d'elements
  7. fl=500 # Intensite de la force lineique en Newtons/metre

Télécharger

4.3- Maillage

Un mailleur est un logiciel qui pour une géométrie en entrée crée un maillage et génère une table de connectivité et une table de coordonnées du maillage.
L’étape du maillage consiste donc à générer ces deux tables.

La table de coordonnées permet de connaitre les coordonnées du noeud i du maillage.

La table de connectivité permet, pour un élément i, de savoir comment sont appelés dans le système global chacun de ces noeuds


Ces deux tables suffisent à définir complètement le maillage.

Dans notre cas, nous choisissons de numéroter comme suit :

Pour n_e nombre d’éléments, il y a donc n_e+1 nombre de noeuds.
La suite cherche donc à générer de manière automatique ces tables pour notre problème.

4.3.1- Table de coordonnées

Dans notre cas, le problème est à une dimension. Une coordonnée d’un noeud est donc uniquement sa distance algébrique à l’origine.
dans le cas où ne=5 et L=1, la table de coordonnées est :

 coor=\pmatrix{
0.0\cr
0.2\cr
0.4\cr
0.6\cr
0.8\cr
1.0\cr}

La table est de taille (ne+1,1). Le plus simple est de partir d’une table de la bonne dimension remplie de 0 puis de la remplir case par case. Pour l’initialiser, on utilise la fonction np.zeros qui remplit exactement ce but. Pour lui indiquer la taille de la matrice à créer :

  1. # Maillage
  2. # Generation de la table coor
  3. coor=np.zeros((ne+1,1)) # Initialisation

Télécharger

Pour ceux qui ont l’habitude de matlab, la différence est qu’il faut rajouter une paire de parenthèse. En effet, np.zeros ne prend qu’un seul paramètre.

Il faut maintenant peupler le vecteur qui pour l’instant est vide. Pour s’en convaincre, exécuter le script puis taper dans la console :

  1. print coor

Le plus simple pour remplir le vecteur est de parcourir les noeuds, et de remplir leur coordonnée dans le vecteur. Il faut donc parcourir ne+1 éléments, il est tout indiqué d’utiliser une boucle for qui aura ne+1 itérations, une pour chaque noeud. On cherche à faire prendre à une variable i les valeurs 1,2,3,...,ne+1. Nous irons insérer dans la i-ème case du vecteur la valeur \frac{i-1}{ne}L
Pour générer une liste de valeur, on utilise la fonction range. Dans le terminal taper :

  1. range(ne+1)

Python répond :

[0, 1, 2, 3, 4, 5]

Or dans notre programme ne vaut 5 et nous voulions avoir la réponse :

[1, 2, 3, 4, 5, 6]

En fait, Python à le mauvais goût, comme le C++ et d’autres langages, de compter à partir de 0, question de conventions.
De la même manière ,pour accéder à la i-ème valeur du vecteur on appellera donc coor[i,0], 0 signifiant la première colonne.

La syntaxe de la boucle for en Python est la suivante :

  1. for variable in possibilites:
  2. action 1
  3. action 2
  4. ...
  5. action n
  6. action en dehors de la boucle

Télécharger

Notez bien les deux points " :" après avoir déclaré la condition de parcours des valeurs. Pour délimiter les actions qui vont être exécutés à chaque itération, il est obligatoire d’indenter (décaler vers la droite) le code. Pour cela il suffit d’appuyer sur la touche [TAB].

Nous voulions insérer dans la i-ème case la valeur \frac{i-1}{ne}L, si i allait de 1 à ne+1. Or pour Python, i va de 0 à ne. Il faut alors insérer la valeur \frac{i}{ne}L dans le vecteur (on a décalé de +1)

La boucle for qui fait prendre à une variable i ces valeurs s’écrit :

  1. for i in range(ne+1):
  2. coor[i,0]=(float(i)/(ne))*L

Télécharger

Pour calculer la valeur à insérer dans la i-ème case,on doit se servir de la valeur de i qui est un entier pour le diviser. Le problème est que pour Python, la division d’un entier donne un entier (division euclidienne). Pour le forcer à considérer i comme un réel, on lui convertit explicitement i en réel en appellant float(i).

4.3.2- Table de connectivité

La table de connectivité se génère en parcourant les éléments cette fois et en listant le numéro des noeuds de l’élément dans le système global.

  1.  
  2. # Generation de la table connec
  3. connec=np.zeros((ne,2))
  4. for i in range(ne):
  5. connec[i,0]=i+1
  6. connec[i,1]=i+2

Télécharger

4.4- Assemblage

Cette étape consiste à créer la matrice de rigidité et le vecteur force. Comme précedemment, on les initialise de la bonne taille à zéro puis on vient les remplir par les contributions de chaque élément.

  1.  
  2. K=np.zeros((ne+1,ne+1))
  3. F=np.zeros((ne+1,1))

Télécharger

On itère ensuite sur chaque élément grâce à une boucle for. On génère pour chaque élément sa matrice de rigidité élémentaire, dans notre cas de taille (2,2). lors de l’assemblage, on "éclate" la matrice élémentaire entre n_1 et n_2, qui sont les numéros des noeuds de l’élément dans le système global, information que l’on trouve dans la table de connectivité.

La première chose à faire dans la boucle est de récupérer ces nombres n1 et n2 dans la table de connectivité, qui sont rangés à la ligne i.

  1.  
  2. # Assemblage
  3. K=np.zeros((ne+1,ne+1))
  4. F=np.zeros((ne+1,1))
  5. for i in range(ne):
  6. n1=connec[i,0]
  7. n2=connec[i,1]
  8.  

Télécharger

On calcule ensuite la matrice de rigidité. Le seul paramètre qui change d’une matrice à l’autre est la longueur l de l’élément. Il faut la calculer, en allant chercher pour cela les coordonnées des noeuds de l’élément.

Comme Python compte à partir de 0, on décale nos valeurs de -1. Par exemple, au lieu d’assigner K[n1,n1] comme on le ferait sous matlab, on est obligés d’écrire K[n1-1,n1-1].

Pour l’assemblage de la force ponctuelle, chaque noeud reçoit la moitié de la force totale subie par la force surfacique.

  1.  
  2. # Assemblage
  3. K=np.zeros((ne+1,ne+1))
  4. F=np.zeros((ne+1,1))
  5. for i in range(ne):
  6. n1=connec[i,0]
  7. n2=connec[i,1]
  8. l=abs(coor[n1-1,0]-coor[n2-1,0])
  9. Ke=E*S/l*np.array([[1,-1],[-1,1]])
  10. K[n1-1,n1-1]+=Ke[0,0]
  11. K[n2-1,n1-1]+=Ke[1,0]
  12. K[n1-1,n2-1]+=Ke[0,1]
  13. K[n2-1,n2-1]+=Ke[1,1]
  14. F[n1-1]+=0.5*fl
  15. F[n2-1]+=0.5*fl
  16.  

Télécharger

Remarquez que la syntaxe de Python pour incrémenter une variable selon le schéma a=a+b s’écrit a+=b, ce qui permet une ecriture plus compacte.

Il ne faut pas oublier d’aditionner la contribution de la force ponctuelle, qui s’ajoute au dernier noeud.

  1. # Ajout de la force pontcuelle
  2. F[ne]+=Fp
  3.  

Télécharger

4.5- Prise en compte des conditions limites

Il faut maintenant imposer le déplacement à certains noeuds. On sépare, quitte à réorganiser la matrice de rigidité les déplacements connus u_c et inconnus u_i. Le problème s’écrit alors :

 
\pmatrix{
K_{ii} & K_{ic} \cr
K_{ci} & K_{cc} \cr}
\pmatrix{
u_i \cr u_c}=
\pmatrix{
F_i \cr F_c}

Si on ne garde que la première ligne du système :

 K_{ii}u_i+K_{ic} u_c=F_i

D’où :

 u_i=K_{ii}^{-1}(F_i-K_{ic} u_c)

Dans notre exemple, u_c est réduit à un seul élément, u_1, déplacement du premier noeud et égal à 0. Dans ce cas particulier, l’expression devient :

 u_i=K_{ii}^{-1}(F_i)

Pour obtenir K_{ii}, il faut supprimer la première ligne et première colonne de K et pour F_i, supprimer le premier élément de F. Ceci s’écrit grâce à la fonction np.delete.

Elle prend respectivement comme paramètre :

  • la matrice ou le vecteur à partir duquel la fonction va travailler
  • la place de l’élément à supprimer. Ici les premières lignes ou colonnes donc 1, donc 0 en convention Python
  • 0 si c’est une ligne à supprimer, 1 si c’est une colonne

Il est important de remarquer que la fonction ne modifie pas la matrice ou le vecteur passé en paramètre. Elle le copie, effectue les opération demandées et renvoie le résultat. C’est pour cela qu’il faut assigner ce résultat :

  1. # Prise en compte des CL
  2. Kii=np.delete(K,0,0)
  3. Kii=np.delete(Kii,0,1)
  4. Fi=np.delete(F,0,0)
  5.  

Télécharger

Cette série d’instruction réalise :

  • copier K et supprimer dans cette copie la première ligne. Assigner le résultat dans une nouvelle variable K_{ii}. Il faut remarquer que K n’a donc pas été modifié !
  • copier K_{ii}, supprimer dans cette copie la première colonne. Assigner le résultat dans K_{ii} pour le mettre à jour
  • copier F, supprimer dans cette copie la. Assigner le résultat dans F_{i} pour le mettre à jour

4.6- Vérification de la construction des matrices et vecteurs

Nous avons effectué des opérations relativement complexes sur des matrices et vecteurs, et on peut rapidement se tromper surtout quand on débute dans un nouveau langage. Pour vérifier les opérations effectuées, on peut vouloir afficher une matrice ou un vecteur par exemple K_{ii}. La manière simple est de taper dans la console ou d’écrire à la fin du script :

  1. print Kii
  2.  

Télécharger

La réponse en console est du genre :

[[  2.00000000e+10  -1.00000000e+10   0.00000000e+00   0.00000000e+00
   0.00000000e+00]
[ -1.00000000e+10   2.00000000e+10  -1.00000000e+10   0.00000000e+00
   0.00000000e+00]
[  0.00000000e+00  -1.00000000e+10   2.00000000e+10  -1.00000000e+10
   0.00000000e+00]
[  0.00000000e+00   0.00000000e+00  -1.00000000e+10   2.00000000e+10
  -1.00000000e+10]
[  0.00000000e+00   0.00000000e+00   0.00000000e+00  -1.00000000e+10
   1.00000000e+10]]

Ce qui n’est pas très visuel. Si vous changez en haut du script ne=20 par exemple, le résulat sera une bouillie de nombres...

C’est pourquoi je vous propose de taper ces lignes après les opérations sur $Kii$ notamment :

  1. plt.figure()
  2. plt.pcolor(Kii)
  3. plt.gca().invert_yaxis()
  4. plt.colorbar()
  5.  

Télécharger

Vous avez alors une visualisation de la matrice, les nombres étant transformés en couleurs.

4.7- Résolution de système linéaire

Finalement, la résolution du système linéaire tient en une seule ligne :

  1. ui=lg.solve(Kii,Fi)
  2.  

Télécharger

4.8- Post-traitement

Nous avons obtenu les u_i par calcul. Il faut maintenant réintégrer les u_c pour former u :

  1. u=np.insert(ui,0,0)
  2. plt.plot(coor,u)
  3.  

Télécharger

L’instruction plot permet de tracer la courbe, en prenant le premier vecteur comme valeurs sur x, le deuxième sur y.

Nous obtenons alors le déplacement on fonction de la coordonnée des noeuds.
Maintenant que l’on a un code éléments finis qui donne des résultats qui semblent cohérents, on cherche à comparer le déplacement obtenu avec le déplacement analytique. En effet, en programmation, il est obligatoire de vérifier le bon fonctionnement du code en essayant sur un cas test dont on connait à l’avance le résultat.

La première équation donnée par la RDM permet la résolution de l’effort normal :

 N(x)=-fx+c_1


Or en bout de poutre, N(L)=F, d’où

 N(x)=F(l-x)

La relation de comportement nous donne le déplacement :

 u(x)=\frac{1}{ES}(Fx-\frac{f(L-x)^2}{2}+c_2)


Or u(0)=0, d’où :

 u(x)=\frac{1}{ES}(Fx+fLx-f\frac{x^2}{2})

Nous voulons donc tracer sur un même graphique le déplacement éléments finis et le déplacement analytique. Nous n’allons pas tracer cette fonction de manière exacte, nous allons l’évaluer en un certain nombre de points. Pour cela on se donne pour x un nombre fini de valeurs grâce à la fonction linspace, qui retourne un vecteur de valeurs régulièrement espacées en donnant la valeur de début, de fin et le nombre de valeurs.

  1. x=np.linspace(0,L,num=100)
  2.  

Télécharger

Nous pouvons alors écrire la formule de uth en utilisant le x ainsi défini. On le trace en superposition avec la solution éléments finis avec une légende pour distinguer les deux courbes avec le code suivant :

  1. uth=1/(E*S)*((Fp+fl*L)*x-fl*np.square(x)/2)
  2.  
  3. plt.figure()
  4. plt.plot(coor,u,x,uth)
  5. legend(('Deplacement EF','Deplacement analytique'),loc=2)
  6.  

Télécharger

Il n’est pas possible de mettre tous les éléments d’un vecteur au carré en tappant x^2 comme cela se fait pour les réels. On utilise alors la fonction np.square

Les deux courbes se superposent, on valide ainsi le calcul éléments finis. On peut remarquer que le résultat est exact à l’endroit ou sont positionnés les noeuds. Entre ceux-ci, on a postulé avec l’approximation éléments finis que le déplacement était linéaire, or la résolution analytique nous montre qu’elle est quadratique. On ne peut donc pas représenter convenablement entre les noeuds la solution.

Le code ainsi créé est une suite d’instructions simples qui travaillent directement sur les variables pour arriver au résultat. Cette technique de programmation est intuitive mais devient rapidement complexe lorsque l’on veut empiler les fonctionnalités. Nous pourrions par exemple vouloir que l’utilisateur puisse définir de manière plus complexe le chargement ou les conditions aux limites. Nous présenterons dans un autre cours une manière différente d’envisager la manière de coder pour permettre d’avoir un code bien plus maintenable et modulaire.

++++

 5- Version finale du code

Pour récapituler, voici la version complète du code :

  1. # -*- coding: utf-8 -*-
  2. """
  3. @author: steven
  4. """
  5.  
  6. # Import des packages
  7. import numpy as np
  8. from numpy import linalg as lg
  9. from matplotlib import pyplot as plt
  10.  
  11. # Paramètres
  12. E=200e9 #Module de young
  13. S=0.01 # Section de la poutre
  14. L=1 # Longueur de la poutre
  15. Fp=100 # Intensité de la force ponctuelle en Newtons
  16. ne=5 # Nombre d'éléments
  17. fl=500 # Intensité de la force linéique en Newtons/mètre
  18.  
  19. # Maillage
  20. # Génération de la table coor
  21. coor=np.zeros((ne+1,1)) # Initialisation
  22. for i in range(ne+1): # Parcours des noeuds
  23. coor[i,0]=(float(i)/(ne))*L
  24. #print coor
  25.  
  26. # Génération de la table connec
  27. connec=np.zeros((ne,2))
  28. for i in range(ne):
  29. connec[i,0]=i+1
  30. connec[i,1]=i+2
  31.  
  32. # Assemblage
  33. K=np.zeros((ne+1,ne+1))
  34. F=np.zeros((ne+1,1))
  35. for i in range(ne):
  36. n1=connec[i,0]
  37. n2=connec[i,1]
  38. l=abs(coor[n1-1,0]-coor[n2-1,0])
  39. Ke=E*S/l*np.array([[1,-1],[-1,1]])
  40. K[n1-1,n1-1]+=Ke[0,0]
  41. K[n2-1,n1-1]+=Ke[1,0]
  42. K[n1-1,n2-1]+=Ke[0,1]
  43. K[n2-1,n2-1]+=Ke[1,1]
  44. F[n1-1]+=0.5*fl*l
  45. F[n2-1]+=0.5*fl*l
  46.  
  47. # Ajout de la force pontcuelle
  48. F[ne]+=Fp
  49.  
  50. # Prise en compte des CL
  51. Kii=np.delete(K,0,0)
  52. Kii=np.delete(Kii,0,1)
  53. Fi=np.delete(F,0,0)
  54.  
  55. #plt.figure()
  56. #plt.pcolor(Kii)
  57. #plt.gca().invert_yaxis()
  58. #plt.colorbar()
  59.  
  60. # Résolution
  61. ui=lg.solve(Kii,Fi)
  62. # On remet l'élément supprimé
  63. u=np.insert(ui,0,0)
  64.  
  65. x=np.linspace(0,L,num=100)
  66. uth=1/(E*S)*((Fp+fl*L)*x-fl*np.square(x)/2)
  67.  
  68. plt.figure()
  69. plt.plot(coor,u,x,uth)
  70. legend(('Deplacement element fini','Deplacement analytique'),loc=2)

Télécharger

Discussion sur cet article

Une question, une réaction sur ce sujet?

modération a priori

Ce forum est modéré a priori : votre contribution n’apparaîtra qu’après avoir été validée par un administrateur du site.

Qui êtes-vous ?
Votre message

Pour créer des paragraphes, laissez simplement des lignes vides.