Cours Ateliers Numériques en Mécanique

Contenu

Cours Ateliers Numériques en Mécanique#

par Marc BUFFAT, dpt Mécanique, Université Lyon 1 [1]

_images/programmation.png — Fig. 1 démarche illustré du calcul scientifique#

1. Introduction#

1.1. Objectif du cours#

L’objectif de ce cours est d’apprendre à maîtriser les outils numériques utilisées pour la modélisation de problèmes en Mécanique

En particulier, on va mettre en oeuvre une méthode, que l’on peut qualifier de démarche du calcul scientifique, dont les étapes sont les suivantes:

Analyse physique du problème
Définition d’un modèle mathématique
Recherche d’une solution (approchée) numérique
Mise en oeuvre d’une solution algorithmique
Programmation sur un ordinateur
Validation et Analyse physique du résultat

Important

L’assimilation du cours nécessite une méthode de travail avec:

une démarche rigoureuse pour méthodes de base
un travail régulier pour savoir appliquer ces méthodes
apprendre à vérifier systématiquement sa démarche
savoir l’expliquer (avec des phrases)
savoir prendre des notes manuscrites pour apprendre !

1.1.1. Outils informatiques:#

Pour mettre en oeuvre cette démarche, on utilisera en particulier les outils informatiques suivants:

Linux, qui est un système d’exploitation de type Unix
des outils de programmation numériques sous Python

avec en particulier les bibliothèques suivantes:

numpy, scipy pour le calcul numérique:
sympy pour le calcul formel:
matplotlib pour la visualisation:

des logiciels Éléments Finis: FreeFem, Comsol

1.1.2. Applications#

Cette démarche sera utilisée pour résoudre divers problèmes de mécanique, dont:

la modélisation d’un treillis
la vibration d’une membrane
l’étude d’une turbo-voile
l’écoulement autour d’une pale d’éolienne

2. Histoire et enjeux de l’informatique scientifique#

boulier chinois

2.1. Historique de l’informatique#

Dès l’apparition des mathématiques, on a cherché à traiter automatiquement des données, ce qui a conduit aux premiers algorithmes

2.1.1. Premiers Algorithmes#

Euclide ~ 300 av JC: premier algorithme pour le pavage d’une pièce (PGCD)

2.1.2. Origine du nom Algorithme#

Mathématicien perse Al-Khawarizmi né vers 783

2.1.3. Prémices de l’informatique#

machine mécanique : la pascaline de Blaise Pascal en 1642

2.1.4. Les pionniers: machine de Turing#

machine éléctro-mécanique

On Computable Numbers article de 1936 d’Alan Turing (1912-1954)
The imitation Game film sur Turing et sa machine Enigma

2.1.5. Premiers ordinateurs électronique#

machine à tubes: UNIVAC I en 1951 (Hergé Objectif lune)
mainframes (transistor) IBM, CDC .. à partir de 1960

2.1.6. Naissance de la micro-informatique ~ 1970#

premier ordinateur personnel français Micral en 1975
TRS 80 à base de microprocesseur 8 bits Zilog Z80 en 1977 (à droite)
Apple II à base de microprocesseur 8 bits Mos 6502 en 1977
Tavernier à base de microprocesseur 16 bits 6809 en 1982 (à gauche)

2.1.7. L’ère internet (à partir 1990)#

Arpanet (1969 réseau militaire) devient internet TCPI/IP (1989)
applications World Wide Web: WWW, HTML, mail (1989 CERN)
premier noyau Linux (Linus Towarld + Unix) 1991
- création de la Free Software Fondation

2.1.8. L’ère actuelle: mobilité , données#

1er IPhone Apple 2007
Supercalculateur IDRIS HPC-IA = 14 PetaFlops ($14. 10^{15}$ ops) (à gauche)
micro-ordinateur ARM Rasberry PI4 14 GigaFlops mais consommation de 6W (à droite)
big DATA, IA

2.1.9. Conclusion#

le numérique (informatique) est maintenant présente partout
puissance du matériel augmente sans cesse (mais pble consommation!)
nécessité de développer des outils logiciels (algorithmique)
- IA aucune intelligence dans la machine, mais dans les algorithmes (traitement de données: Big Data)
nécessité de maîtriser et comprendre les outils pour
- modéliser (simulation numérique)
- analyser (les données expérimentales)

2.2. Utilité d’un ordinateur en science#

Un ordinateur exécute très rapidement des instructions élémentaires: $\approx 4*30. 10^9$ FLOPS sur un ordinateur de bureau (Intel Core I7 4 coeurs). Et il les exécute de façon mécanique.

Question: Comment peut-on utiliser un ordinateur pour simuler un problème physique ?

Par exemple calculer le mouvement d’un pendule, puisque l’ordinateur ne connaît ni la mécanique, ni les équations, ni les méthodes de résolutions de ces équations !

2.2.1. Modélisation numérique#

le scientifique imagine un algorithme pour résoudre le problème de façon mécanique, et le traduit ensuite dans un langage de programmation pour être exécuté par un ordinateur. Il peut ensuite faire l’étude paramétrique du problème comme avec une expérience.

2.2.2. Méthode#

Problème physique

Modèle mathématique: mise en équation

\[ m l \frac{d^2 \theta}{d t^2} = m g \sin(\theta) \]

\[ Y = [ \theta, \frac{d \theta}{dt} ] \]

\[ Y_0=[\theta_0, 0]\]

\[ \frac{d Y}{dt} = F(Y,t) \]

Solution Algorithme RK2

    ....
    Y = Y0
    t = 0
    Pour i de 1 a n
        solution a t+dt/2
        Y12 = Y + 0.5 * dt *F(Y,t)
        solution a t+dt
        Y = Y +  dt *F(Y12,t+dt/2)
        t=t+dt
    Fin Pour

Programmation Python: traduction dans un langage de programmation

  .... 
  Y = np.array([theta0, 0.])
  t = 0.0
  for i in range(n):
       # prediction Y(t+dt/2)
       Y12 = Y + 0.5*dt*F(Y,t)
       # solution Y(t+dt)
       Y = Y + dt*F(Y12,t+0.5*dt)
       t = t + dt
  # fin

2.3. Utilisation de l” informatique en science#

2.3.1. modélisation en mécanique de système complexe#

2.3.2. expérience numérique fondamentale#

DNS simulation numérique directe en mécanique des fluides

2.3.3. CFD = Computational Fluid Dynamics#

simulation numérique à l’ONERA pour Airbus Industrie

2.3.4. CFD = Color Fluid Dynamics ?#

2.3.5. Attention au mirage de la simulation numérique#

Les logiciels sont de plus en plus sophistiqués et fournissent toujours une réponse
Mais la réponse est-elle correcte ?
Il est important de comprendre les principes de base
Importance de la validation des simulations !

Un problème majeur est que la plupart des scientifiques n’ont reçu presque aucune formation en informatique, ce qui devient désormais problématique car la plupart d’entre eux passent une grande partie de leur temps à utiliser des ordinateurs, pour de la modélisation, l’acquisition et le traitement de données expérimentales

acculturation l’IA $\rightarrow$ fournit des corrélations mais pas d’explication !

échecs dans la modélisation

lors du 1er tir d’Ariane V en 1996 , un bug a détruit la fusée après le décolage
destruction d’une plateforme pétrolière concue entièrement sur ordinateur (erreur de modélisation)
biais sur les données d’apprentissage en IA (W.A.S.P White Anglo-Saxon Male Protestant)

2.4. FIN de la leçon#

3. Python pour les scientifiques#

Cette partie est une introduction à la programmation sous Python pour les scientifiques.

I learned Python last night! Everything is so simple! I just typed:

   import antigravity

Python COMICS d'après le site XKCD

Python COMICS d’après le site XKCD

Contenu

3.1. Introduction au calcul scientifique avec Python#

Marc BUFFAT, dpt mécanique, Université Lyon 1

(**) voir aussi le cours « Scientific Python Lectures » by Robert Johansson

Mise à disposition selon les termes de la Licence Creative Commons
Attribution - Pas d’Utilisation Commerciale - Partage dans les Mêmes Conditions 2.0 France.

%matplotlib inline
import numpy as np
import scipy as sp
import matplotlib.pyplot as plt
# option de mise en page
from IPython.core.display import HTML,display

/tmp/ipykernel_268606/3841631963.py:6: DeprecationWarning: Importing display from IPython.core.display is deprecated since IPython 7.14, please import from IPython.display
  from IPython.core.display import HTML,display

3.1.1. Objectifs#

3.1.1.1. Apprendre:#

concept fondamentaux en programmation
différentes techniques de base
démarche de programmation scientifique
syntaxe du langage Python

3.1.1.2. À la fin du cours, vous pourrez:#

résoudre des problèmes scientifiques (à l’aide de Python)
appliquer la démarche de programmation scientifique avec d’autres langages (matlab, scilab, ..)

3.1.2. Fonctionnement d’un ordinateur#

ordinateur = matériel (CPU, mémoire, disque,..) + logiciel (OS, programmes,..)

3.1.2.1. composants d’un ordinateur#

OS

Système d’exploitation (OS) = gestion d’un ordinateur + environement

Linux/Unix = vue abstraite et générique indépendante du matériel

programmes + données

programme = processus
données = fichiers (files, devices, streams)
gestion des droits: notion d’utilisateur (user), groupes (group), autres (other)
super-utilisateur: root

Principes

small is beautiful
generality is better than specificity
communication is important (format ouvert pour les données)

3.1.2.2. programmes / processus#

processus = programme executé sur un ordinateur

pid = numéro du processus
programme (utilisateur ou commande système)
E/S: entrée par défaut (input stream=stdin), sortie par défaut (output stream=stdout)
sortie programme = entrée d’un autre programme : pipe ou |
pge 1 | pge 2 | pge 3
controle des processus: kill pid ou ctrl+C, stop (ctrl+Z) , tache de fond (bg) , ps, top
éxecution des commandes à travers une interface
- un terminal avec un shell (interpérteur de commande): bash (avec auto-complétion)
- un GUI (gnome, KDE, ..)
syntaxe (attention majuscule # minuscule, espace = séparateur d’arguemnts)

cde [-options] arguments

%%bash
ps -au buffat | head

    PID TTY          TIME CMD
 245164 ?        00:00:00 systemd
 245165 ?        00:00:00 (sd-pam)
 

245172 ?        00:00:00 pipewire
 245173 ?        00:00:00 pipewire-media-
 245174 ?        00:06:2

0 pulseaudio
 245177 ?        00:00:00 gnome-keyring-d
 245181 ?        00:00:06 dbus-daemon
 245187

 tty2     00:00:00 gdm-x-session
 245189 tty2     00:06:44 Xorg

3.1.2.3. Système de fichiers (Unix)#

vu arborescente abstraite indépendante du matériel à partir de la racine (root) /
notation
- fichier = nom de fichier + répertoire + droits d’accés
- / séparateur de répertoire (directory)
- . répértoire courant
- .. répértoire parent
- pwd : affiche répertoire courant
- lien symbolique
- fichiers cachés: commence par .
principales commandes
- ls : liste des fichiers
- mkdir : création répertoire
- rm : éfface un fichier
- cd : change de répertoire (directory)
- help cde : aide
- cat file : affiche le contenu du fichier

%%bash
ls -al | head

total 21364
drwxrwxr-x 6 buffat buffat     4096 avril  9 09:36 .
drwxrwxr-x 4 buffat buffat     4096

 févr. 10  2025 ..
lrwxrwxrwx 1 buffat buffat       39 févr. 10  2025 Algorithme.ipynb -> MGC1061M

/04_Algorithme/Algorithme.ipynb
-rw-rw-r-- 1 buffat buffat     1854 févr. 10  2025 alunissage.py
-r

w-rw-r-- 1 buffat buffat     2136 févr. 10  2025 anim_pendule.py
-rw-rw-r-- 1 buffat buffat    2543

0 févr. 10  2025 BaseProgrammation.ipynb
lrwxrwxrwx 1 buffat buffat       39 févr. 10  2025 BasePy

thon.ipynb -> MGC1061M/02_BasePython/BasePython.ipynb
-rw-rw-r-- 1 buffat buffat    34861 févr. 10 

 2025 BasePython.md
-rw-rw-r-- 1 buffat buffat    27355 févr. 10  2025 BibliothequeScientifique.ipy

nb

3.1.3. Execution d’un programme#

3.1.3.1. Calcul de la somme d’une série:#

\[ S = x - \frac{x^2}{2} + \frac{x^3}{3} + .. = \sum_{i=1}^n (-1)^{i+1} x^i/i \]

démarche:

algorithme
programmation
execution du code machine sur l’ordinateur

3.1.3.2. Algorithme de base#

Algorithme Serie(x,n)
   somme = 0
   pour i de 1 a n
      somme=somme + (-1)^(i+1)*x^i/i
   retour somme

Ecrire un algorithme plus efficace en notant que:

\[ (-1)^{i+1} \frac{x^i}{i} = - \frac{x}{i} (-1)^{i} x^{i-1} \]

3.1.3.3. Execution sur un ordinateur#

L’ordinateur (le CPU) n’éxécute que des instructions machines codées en binaire.

Problème: comment transformer un algorithme en langage machine ?

3.1.3.4. Langage machine (binaire)#

suite d’instruction binaire executer par le processeur: 01010111 (suite 0 ou 1 : bits)
représentation en hexadécimal (octet ou bytes):
- f8 = 248
instructions et données sont stockées en mémoire (à une adresse donnée)
- @000141 instruction f8
spécifiques au processeur (intel, ARM, power, AMD, ..)
spécifique au système d’exploitation OS

%%bash
od -x -N 100 data/a.out
./data/a.out

0000000 457f 464c 0102 0001 0000 0000 0000 0000
0000020 0002 003e 0001 0000 0870 0040 0000 0000
0000

040 0040 0000 0000 0000 21b0 0000 0000 0000
0000060 0000 0000 0040 0038 0009 0040 001e 001b
0000100 

0006 0000 0005 0000 0040 0000 0000 0000
0000120 0040 0040 0000 0000 0040 0040 0000 0000
0000140 01f8

 0000
0000144

Calcul de la serie pour n=20 et x=0.2
Somme de la serie = 0.182322
Log(1+x)=0.182322

3.1.4. Langages de programmation#

conversion d’une programme en code machine
plusieurs centaine de langages différents

-langage le plus complexe: malbolge

3.1.4.1. Langage assembleur#

L’assembleur traduit le code source en code binaire

notation des instructions machines (move, push, load, add, mul)
programme + assembleur = code machine

%%bash
head -15 data/serie.s

	.file	"serie.C"
	.local	_ZStL8__ioinit
	.comm	_ZStL8__ioinit,1,1
	.text
	.globl	_Z5seriedi
	.type	_

Z5seriedi, @function
_Z5seriedi:
.LFB971:
	.cfi_startproc
	pushq	%rbp
	.cfi_def_cfa_offset 16
	.cfi_

offset 6, -16
	movq	%rsp, %rbp
	.cfi_def_cfa_register 6
	movsd	%xmm0, -40(%rbp)

3.1.4.2. Langage compilé: compilateur#

le compilateur transforme le code source en assembleur (compilation), éffectue l’édition avec des librairies pour obtenir un binaire qui dépend de l’ordinateur cible, mais qui s’exécute sans le compilateur et sans le code source.

C++ Fortran C
programme source + compilateur = programme binaire
programme source = portable
programme binaire = spécifique à l’ordinateur (non portable)

_images/C.jpeg langage C

%%bash
head -17 data/serie.C

#include <stdlib.h>
#include <math.h>
#include <iostream>
//
// calcul de la somme de n termes de la

 serie x-x^2/2+x^3/3-
//
double serie(double x,int n)
{
	double coef=x;
	double somme=0.0;
	for (int

 i=1; i<=n; i++)
	{
		somme = somme + coef/i;
		coef = -coef*x;
	}
	return somme;
}

%%bash
g++ data/serie.C -o serie
./serie

Calcul pour n=20 et x=0.2
Somme   = 0.182322
Log(1+x)= 0.182322

3.1.4.3. Langage interprété:#

l’interpéteur exécute de façon interactive chaque ligne du fichier pour la traduire en code machine et l’execute interactivement. Pour exécuter le code, on a donc toujours besoin de l’interpréteur et du code source.

Python, Matlab
programme + interpréteur
portable = indépendant du système (Linux, MacOS, Windows)

3.1.4.3.1. Matlab (Matrix Laboratory)#

Matlab est un logiciel propriétaire qui est une boite à outils (toolbox) de calcul matriciel et numérique développé par Mathworks Inc. Il fonctionne sous Windows, Mac et Unix et permet d’écrire rapidement des scripts pour faire du calcul scientifique.

De nombreuses alternatives libres et de qualité existent, et notamment les logiciels : Scilab, Octave et surtout Python avec numpy.

%%bash
cat data/serie.m

% programme matlab
% ================
x=0.2; n=20;
% calcul de la serie
coef=x;
somme=0.0;
for i=1:n

;
somme = somme+coef/i;
coef  = -coef*x;
end;
display(sprintf('Calcul de la serie pour x=%g et n=%d'

,x,n))
display(sprintf('somme    = %g',somme))
display(sprintf('log(1+x) = %g',log(1+x)))
exit

3.1.4.3.2. langage python#

Python est un language de programmation moderne de haut niveau, logiciel libre, généraliste, multi-plateformes, multi-architecture, multi-OS (linux, windows, MacOS,..). Un programme python peut s’executer sur un tout petit raspberry Pi à 30€ (ARM), sur des smartphones, des portables, des PC jusqu’aux super-calculateurs HPC avec $10^9$ coeurs de calcul.

%%bash
cat data/serie.py

#! /usr/bin/env python
from math import log

def serie(x,n):
    """ calcul de la somme de n termes 

de la serie x-x^2/2+x^3/3- """
    coef=x
    somme=0.0
    for i in range(1,n+1):
        somme = s

omme + coef/i
        coef  = -coef*x
    return somme
#
x=0.2
n=20
print("Calcul de la serie pour x

=",x," et n=",n)
print("somme    = ",serie(x,n))
print("log(1+x) = ", log(1+x))

%%bash
python data/serie.py

Calcul de la serie pour x= 0.2  et n= 20
somme    =  0.18232155679395456
log(1+x) =  0.1823215567939

3.1.5. Historique des langages de programmation#

Années 1950 (approches expérimentales) :
- FORTRAN, LISP, COBOL…
Années 1960 (langages universels) :
- ALGOL, PL/1, PASCAL…
Années 1970 (génie logiciel) :
- C, MODULA-2, ADA…
Années 1980 (programmation objet) :
- C++, Eiffel…
Années 1980 (boites à outils):
- LabView, Matlab…
Années 1990 (langages interprétés objet) :
- Java, Perl, Python…

3.1.5.1. Langages les plus utilisées par les scientifiques#

Pour le calcul intensif (HPC) (Argonne Lab)

C,C++
Fortran
Python

Matlab est très utilisé pour le traitement des données, mais peut être avantageusement remplacer par Python avec ses librairies numpy et matplotlib et pandas

3.1.6. Exigences du calcul scientifique#

« Le logiciel (software) est une des pierres angulaires de la science moderne. Sans logiciel, la science du vingt et unième siècle serait impossible. Mais sans de meilleurs logiciels, la science ne peut pas progresser » [http://sciencecodemanifesto.org/]

Réplicabilité et de reproductibilité sont des pierres angulaires de la méthode scientifique. En ce qui concerne le travail numérique, se conformer à ces concepts a des implications pratiques suivantes:

Réplicable: L’auteur d’un article scientifique qui implique des calculs numériques doit être en mesure de relancer les simulations et de reproduire les résultats sur demande. D’autres scientifiques devraient également être en mesure d’effectuer les mêmes calculs et d’obtenir les mêmes résultats, compte tenu des informations sur les méthodes utilisées dans une publication.
Reproductible: Les résultats obtenus à partir de simulations numériques doivent être reproductibles avec une mise en oeuvre indépendante du procédé, ou en utilisant un procédé tout à fait différent.

En résumé: un résultat scientifique solide doit être reproductible, et une étude scientifique solide doit être réplicable.

importance de la documentation
importance de la validation
importance de la maîtrise du code et des transformations de données (logiciels libres)

3.1.6.1. Démarche du calcul scientifique#

Modèle mathématique
Discrétisation numérique: solution approchée
Recherche solution algorithmique
Programmation
Validation

3.1.7. Documentation#

Python. The official Python web site.
Python tutorials. The official Python tutorials.
Think Python. “”How to Think Like a Computer Scientist”” by Allen B. Downey (free book).
COURS InProS “”Introduction à la Programmation Scientifique””

3.1.7.1. Installation (logiciel libre)#

Distribution pour Linux, Windows, Mac
Distribution Anaconda
accès aux serveurs Jupyter à l’UCB Lyon 1:

3.1.8. Fin de la leçon#

3.2. Base de programmation en Python#

Marc BUFFAT, dpt mécanique, Université Lyon 1 et [1]

[1] inspiré librement du cours « Engineering Computations » du Pr L. Barba (Washington Univ.) et des Monty Python scriptol.fr

Python & Monty Python

from IPython.display import display, Markdown, clear_output
def printmd(string):
    display(Markdown(string))
from metakernel import register_ipython_magics
register_ipython_magics()

3.2.1. Questions préliminaires#

Cliquez sur Next pour répondre à la question suivante

#%activity /usr/local/commun/ACTIVITY/IntroPython/questionLangage

#%activity /usr/local/commun/ACTIVITY/IntroPython/questionProg

#%activity /usr/local/commun/ACTIVITY/IntroPython/questionPython

#%activity /usr/local/commun/ACTIVITY/IntroPython/questionProgrammation

3.2.2. Historique#

Python a été créé en 1989 par Guido van Rossum en hommage aux Monty Python

3.2.3. Pourquoi Python#

3.2.3.1. Avantages du langage Python#

Langage interprété
Usage général : on peut tout faire
- interfaces graphiques
- calcul scientifique
- applications webs
- base de données
- IA
- robotique
Vaste librairie de modules
Syntaxe cohérente
- langage orienté objet
- langage fonctionnel
Facile à apprendre / agréable à utiliser
Excellent premier langage

3.2.3.2. Défauts du langage Python#

selon les circonstances, les programmes écrits en Python peuvent comporter des problèmes de performance
c’est le cas de tous les langages interprétés

3.2.3.3. Approche générale de développement#

on commence à développer en Python
on identifie éventuellement les sections de code qui posent un problème de performance
- on réécrit les sections de code problématiques dans un autre langage tel que le C/C++
Python est conçu pour s’interfacer facilement aux autres langages
utilisation des bibliothèques

3.2.4. Philosophie de Python (Zen of python)#

import this

The Zen of Python, by Tim Peters

Beautiful is better than ugly.
Explicit is better than implicit.
Simple is better than complex.
Complex is better than complicated.
Flat is better than nested.
Sparse is better than dense.
Readability counts.
Special cases aren't special enough to break the rules.
Although practicality beats purity.
Errors should never pass silently.
Unless explicitly silenced.
In the face of ambiguity, refuse the temptation to guess.
There should be one-- and preferably only one --obvious way to do it.
Although that way may not be obvious at first unless you're Dutch.
Now is better than never.
Although never is often better than *right* now.
If the implementation is hard to explain, it's a bad idea.
If the implementation is easy to explain, it may be a good idea.
Namespaces are one honking great idea -- let's do more of those!

3.2.5. Concepts de base#

3.2.5.1. Algorithme#

Un algorithme est une suite finie et non ambigüe d’opérations ou d’instructions permettant de résoudre un problème. Les algorithmes sont connus depuis l’antiquité (Euclide).

Le mot algorithme vient du nom du mathématicien perse du 9ième siècle (AJC)

Abu Abdullah Muhammad ibn Musa al-Khwarizmi (photo ci dessus).

L’algorithmique correspond à la phase préparatoire avant une quelconque programmation. Elle permet de décrire un problème sous une forme que l’on peut ensuite programmer sur un ordinateur et ceci dans un langage naturel, indépendant d’un langage de programmation.

algorithme numérique suite finie et non ambiguë d’opérations ou d’instructions sur des nombres permettant de résoudre un problème.

Et il n’est pas nécessaire d’avoir un ordinateur pour exécuter un algorithme (ex: machine de Turing inventé en 1936 avant l’ordinateur)!

3.2.5.2. Programme informatique#

ORDINATEUR est une machine qui exécute (très rapidement) des algorithmes sous la forme de programmes informatiques

PROGRAMME est une suite non ambigüe d’instructions permettant de traiter des données stockées dans des variables pour obtenir un résultat.

VARIABLE est le nom (symbol) donné à un espace de stockage dans la mémoire de l’ordinateur pour stocker des données numériques. Attention une variable informatique n’est pas une variable mathématique (pourquoi ?).

INSTRUCTION est une étape élémentaire décrit dans un langage de programmation avec une syntaxe rigoureuse.

FONCTION permet est un regroupement d’instructions qui traitent des données qui sont les arguments de la fonction pour calculer un résultat qui est la valeur de la fonction. Attention une fonction informatique n’est pas une fonction mathématique (pourquoi ?).

CODE MACHINE est une suite d’instructions binaires exécutée par le processeur. La traduction des instructions en code machine se fait automatiquement soit avec un interpréteur pour les langages interprétés (comme Python ou Matlab) ou avec un compilateur pour les langages compilés (comme C ou C++)

3.2.6. Les chaînes de caractères en Python#

Les chaînes de caractères sont un exemple de données textuelles que sait manipuler un ordinateur. Pour manipuler des caractères sur un ordinateur, il faut traduire les caractères en nombre, en utilisant un code. Le premier système de codage date des années 1960: c’est le code ASCII American Standard Code for Information Interchange, qui est encore la norme de codage utilisée pour les caractères latins sans accents. Pour pouvoir coder les accents ou les caractères non latin, on utilise maintenant un codage standard UTF-8 ou UTF-16 Universal Character Set Transformation Format sur 8 ou 16 bits, qui inclut le codage ASCII.

Sous Python, une chaîne est composée d’une séquence de caractères entre deux guillemets, et elle est de type str (string en anglais).

3.2.6.1. exercice#

Créer 2 variables mot_1 et mot_2 contenant respectivement Bonjour et Python et créer une nouvelle variable ma_chaîne contenant la concaténation des 2 variables avec un espace entre les deux.

   mot_1 = "Bonjour"
   mot_2 = "Python"

Comment peut-on ajouter un point d’exclamation !à la fin de cette chaîne.

# écrire ici votre solution
mot_1 = ""
mot_2 = ""
ma_chaine = ""

# solution
mot_1 = "Bonjour"
mot_2 = "Python"
ma_chaine = mot_1 + " " + mot_2 + "!"
print(ma_chaine)

Bonjour Python!

3.2.6.2. Indexation#

Nous pouvons accéder à chacun des caractères d’une chaîne de caractère en utilisant un indice, c.a.d un entier indiquant la position du caractère par rapport au début de la chaîne. Les indices sont placés entre crochets après le nom de la variable de chaîne et commencent à partir de 0. Par exemple, le 1er élément (caractère) de ma_chaine correspond à ma_chaine[0], et le 3ième à ma_chaine[2]. On peut aussi compter à partir du dernier caractère en utilisant des indices négatifs. Ainsi le dernier caractère correspond à ma_chaine[-1] et l’avant dernier à ma_chaine[-2].

Et Oui! en Python, nous commençons à compter à partir de 0 (comme en C/C++).

la fonction len(chaine) renvoie le nombre de caractères de la chaîne.

exercice

afficher dans la cellule suivante, le premier et second caractère de ma_chaine, ainsi que l’avant-dernier et le dernier.

# écrire ici votre solution

# solution
ma_chaine[0],ma_chaine[1],ma_chaine[-2],ma_chaine[-1]

('B', 'o', 'n', '!')

3.2.6.3. sous-chaînes#

Pour sélectionner plusiers caractères consécutifs d’une chaîne, on utilise la notation slicing entre crochets avec : «[debut:fin]», où «debut» est l’index pour commencer la tranche, et «end» est l’index (non inclusif,i.e. non inclus) pour terminer la tranche. Par exemple, pour saisir le mot Bonjour dans notre chaîne, on utilise:

    ma_chaine[0:7]

si debut est omis, il correspond au premier (i.e. 0) et si fin est omis, il correspond au dernier caractère.

exercise

Dans la variable ma_chaine, sélectionnez le mot Python en utilisant le slicing. De même sélectionnez toute la chaine avec un slicing le plus simple possible.
Définissez une variable b contenant banana et afficher le premier et le dernier a en utilisation une indexation de b.
De même, saisissez les 2 slices possibles qui correspondent au mot ana et affichez-les.

# écrire ici votre solution

# solution 1
print(ma_chaine[8:-1])
print(ma_chaine[:])
# solution 2
b = 'banana'
print(b[1],b[-1])
# solution 3
print(b[1:4],b[3:])

Python
Bonjour Python!
a a
ana ana

3.2.6.4. manipulation des chaînes de caractères#

Python a de nombreuses fonctions intégrées utiles pour les chaînes.Un détail technique : en Python, certaines fonctions sont associées à une classe particulière d’objets (par exemple, des chaînes ou string). Le mot méthode est utilisé dans ce cas, et nous avons une nouvelle façon de les appeler avec l’opérateur point. C’est un peu contre-intuitif en ce que le nom de la méthode vient après le point, tandis que le nom de l’objet particulier sur lequel elle agit vient en premier. Comme ceci: ma_chaine.method().

3.2.6.4.1. commentaires#

on peut ajouter des commentaires dans un code python en utilisant le caractère #

3.2.6.4.2. application#

Utilisons une citation d’Albert Einstein comme chaîne et appliquons quelques méthodes de chaîne utiles.

 AE_quote = "Tout le monde est un génie. Mais si vous jugez un poisson sur ses capacités à grimper à un arbre, il passera sa vie à croire qu'il est stupide."

# création d'une chaine contenant la citation
AE_quote = "Tout le monde est un génie. Mais si vous jugez un poisson sur ses capacités à grimper à un arbre, il passera sa vie à croire qu'il est stupide."

3.2.6.4.3. `count()`#

La méthode count() donne le nombre d’occurrences d’une sous-chaîne dans une plage. Les arguments de la plage sont facultatifs.

Syntaxe:

str.count(sous-chaîne, début, fin)

Ici, «debut» et «fin» sont des entiers qui indiquent les indices où commencer et terminer le comptage. Par exemple, pour calculer combien de lettres e il y a dans AE_quote, il suffit d’écrire:

AE_quote.count('e')

3.2.6.4.4. `find()` & `index()`#

La méthode find() nous indique si une chaîne substr apparaît dans la chaîne sur laquelle nous appliquons la méthode. Les arguments de la plage sont facultatifs.

Syntaxe:

str.find(substr, debut, fin)

«debut» et «fin» sont des indices indiquant où commencer et terminer la tranche sur laquelle appliquer la méthode find().

Si la chaîne substr est dans la chaîne d’origine, la méthode find() retournera l’index où la sous-chaîne commence, sinon elle retournera -1.

Par exemple, trouvons le mot «poisson» dans la citation d’Albert Einstein et sélectionnons ensuite le mot dans la phrase

AE_quote.find('poisson')

AE_quote[50:50+len('poisson')]

'poisson'

3.2.6.4.5. `index()`#

Une méthode similaire est index(): elle fonctionne comme la méthode find(), mais renvoie une erreur si la chaîne que nous recherchons n’est pas trouvée.

Syntaxe:

str.index(substr, début, fin)

Essayez les instructions suivantes:

    AE_quote.index('poisson')
    AE_quote.index('albert')

la dernière instructions renvoie une erreur, car albert n’est pas un mot de la pharse

# écrire ici votre solution

# solution
AE_quote.index('poisson')

3.2.6.4.6. exercices:#

Utilisez la méthode count() pour compter le nombre de lettres 'a' dans AE_quote?
En utilisant la même méthode, combien de lettres isolées 'à' se trouvent dans AE_quote?
Utilisez la méthode index() pour trouver la position des mots « génie », « jugez » et « arbre » dans AE_quote.
À l’aide d’un slicing, extrayez les mots de la question précédente dans AE_quote.

#  écrire ici votre solution

# solution
print("nbre de a=",AE_quote.count('a'))
print("nbre de à isolé =",AE_quote.count(' à '))
print(AE_quote.index("génie"),AE_quote.index("jugez"),AE_quote.index("arbre"))
# remplacer la valeur de mot par le mot choisie
mot = "génie"
pos = AE_quote.index(mot)
print(AE_quote[pos:pos+len(mot)])

nbre de a= 7
nbre de à isolé = 3
21 41 91
génie

3.2.6.4.7. `split()`#

Cette méthode split() renvoie une liste de tous les mots d’une chaîne. Nous pouvons également définir un séparateur autre que l’espace et diviser notre chaîne en fonction de ce séparateur, et éventuellement limiter le nombre de divisions à «num».

Syntaxe:

str.split(séparateur, num)

print(AE_quote.split())

['Tout', 'le', 'monde', 'est', 'un', 'génie.', 'Mais', 'si', 'vous', 'jugez', 'un', 'poisson', 'sur', 'ses', 'capacités', 'à', 'grimper', 'à', 'un', 'arbre,', 'il', 'passera', 'sa', 'vie', 'à', 'croire', "qu'il", 'est', 'stupide.']

3.2.7. Les listes en Python#

Les crochets ci-dessus indiquent une liste (list) en Python. Une liste est un type de données constitué d’une séquence de valeurs, par exemple des nombres ou des chaînes. Les listes fonctionnent à bien des égards de la même manière que les chaînes: leurs éléments sont numérotés à partir de zéro, le nombre d’éléments est donné par la fonction len (), ils peuvent être manipulés avec la notation de slicing.

Le moyen le plus simple de créer une liste consiste à placer une séquence de valeurs séparées par des virgules entre crochets.

Les éléments d’une liste ne sont pas forcément du même type, et on peut mélanger des entiers, des réels et des chaînes de caractères.

# une liste d'entiers
entiers = [1, 2, 3, 4, 5]
print(entiers)

[1, 2, 3, 4, 5]

# une liste avec des types d'éléments differents 
ma_liste = [2, 'pomme', 4.5, [5, 10]]
print(ma_liste)

[2, 'pomme', 4.5, [5, 10]]

3.2.7.1. exercice:#

Dans la liste entiers, écrire le slicing contenant «[2, 3, 4]» puis «[4, 5]».

# écrire ici votre solution

entiers[1:4], entiers[3:]

([2, 3, 4], [4, 5])

3.2.7.2. `append()`: ajout d’éléments#

Nous pouvons ajouter des éléments à une liste en utilisant la méthode append(): elle ajoute l’objet que nous passons dans la liste existante. Par exemple, pour ajouter l’élément 6 à notre liste entiers, nous pouvons écrire:

   entiers.append(6)

# écrire ici votre solution

# solution
entiers.append(6)
entiers

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

3.2.7.3. `in` : test si dans la liste#

La vérification de l’appartenance à une liste en Python ressemble assez à un anglais simple!

Syntaxe:

Pour vérifier si un élément est dans une liste:

element in liste

Pour vérifier si un élément n’est pas dans une liste:

element not in liste

3.2.7.3.1. exercice#

Étant donné la liste alist = [1, 2, 3, '4', [5, 'six'], [7]] exécutez ce qui suit dans des cellules séparées et discutez le résultat:

in alist
in alist
in alist 
   [7] in alist

alist = [1, 2, 3, '4', [5, 'six'], [7]]
# écrire ici votre solution

# solution
print(4 in alist)
print(5 in alist)
print(7 in alist)
print([7] in alist)

False
False
False
True

3.2.7.4. Modifier les éléments d’une liste#

Nous pouvons non seulement ajouter des éléments à une liste, mais également modifier un élément spécifique. Réutiliser la liste de l’exercice ci-dessus en remplaçant le caractère '4' par l’entier 4.

     print(4 in alist)
     alist[alist.index('4')] = 4
     print(4 in alist)

Remarque: on peut modifier avec une affectation un élément d’une liste mais pas un caractère d’une chaîne. Une chaîne de caractère est par défaut immuable.

# écrire ici votre solution

# solution
print(alist)
print(4 in alist)
alist[alist.index('4')] = 4
print(alist)
print(4 in alist)

[1, 2, 3, '4', [5, 'six'], [7]]
False
[1, 2, 3, 4, [5, 'six'], [7]]
True

3.2.8. Itérations et tests#

Un ordinateur est très efficace pour répéter une série d’instructions. Dans les langages de programmation, c’est la notation d’itération et de boucle.

3.2.8.1. Itération avec `for`#

L’idée de l’itération est de répéter des instructions plusieurs fois. Si vous connaissez un autre langage de programmation C/C++ ou java, vous savez créer une itération avec des instructions for. Mais l’instruction for est un peu différente en Python.

En Python, l’instruction for permet d’itérer sur les éléments d’une séquence (ou liste). Supposons que l’on a créé une liste Fruits contenant des noms de fruits, on peut écrire

    for  fruit in fruits:
        faire quelque chose avec chaque élément (fruit) de la liste

Ici, pour la première fois, nous rencontrerons une particularité du langage Python: l’indentation. Pour délimiter ce que Python doit faire avec chaque fruit dans la liste des fruits, nous plaçons la ou les déclarations suivantes en retrait à partir de la gauche.

Comment indenter ? C’est une question de style, et chacun a une préférence : deux espaces, quatre espaces, une tabulation… sont tous des styles valables : mais choisissez un style et soyez cohérent !

un conseil : utiliser plutôt des espaces qu’une tabulation. Nous choisirons dans la suite une indentation de 4 espaces.

fruits = ['pomme', 'banane', 'orange', 'cerise', 'mandarine']

for fruit in fruits:
    print("Manger une ",fruit)

Manger une  pomme
Manger une  banane
Manger une  orange
Manger une  cerise
Manger une  mandarine

3.2.8.2. Attention:#

l’instruction for se termine par un deux-points,:
la variable fruit est implicitement définie dans l’instruction for
fruit prend la valeur (chaîne) de chaque élément de la liste fruits, dans l’ordre
l’instruction indentée print() est exécutée pour chaque valeur de fruit
une fois que Python est à court de fruits, il s’arrête
nous n’avons pas besoin de savoir à l’avance le nombre d’articles dans la liste!

3.2.8.3. Exercice:#

A partir de la liste de listes (c.a.d, une liste imbriquée) suivante:

prenom_noms = [['samuel', 'dupont'], ['zoe', 'martin'], ['nael', 'abu ammar'], ['thomas', 'perez']]

Écrire le code qui crée deux listes simples : une avec les prénoms, une autre avec les noms de la liste imbriquée ci-dessus, mais en majuscules.

Pour commencer, vous devez créer deux listes vides (empty) en utilisant les crochets sans rien à l’intérieur et utiliser les méthodes de liste append ().

    noms = []
    prenoms = []

Nous l’avons fait pour vous ci-dessous. (conseil: utilisez la méthode de liste append ()!)

prenom_noms = [['samuel', 'dupont'], ['zoe', 'martin'], ['nael', 'abu ammar'], ['thomas', 'perez']]
prenoms = []
noms    = []
# ecrire votre code ici

for pnom in prenom_noms:
    prenoms.append(pnom[0])
    noms.append(pnom[1])
print(prenoms,noms)

['samuel', 'zoe', 'nael', 'thomas'] ['dupont', 'martin', 'abu ammar', 'perez']

3.2.8.4. Test avec `if`#

Parfois, on a besoin de vérifier une condition et modifier le comportement du programme en fonction de la condition. C’est le but de l’instruction if, qui peut prendre l’une des trois formes suivantes:

if: test simple
```
 si condition alors instruction1
```

if-else : test avec 2 branches

si condition alors instruction1 
sinon instruction2

if-elif-else: test le plus complet

si condition1 alors instruction1 
sinon-si condition 2 alors instruction2
sinon instruction 3

# test simple
a = 8 
b = 3
if a > b:
    print('a est plus grand que b')

a est plus grand que b

# test avec 2 branches
x = 1547
if x % 17 == 0: 
    print('x est un multiple de 17.')
else:
    print('x n"est pas un multiple de 17.')

x est un multiple de 17.

a = 3
b = 5
if a > b:
    print('a est plus grand que b')
elif a < b:
    print('a est plus petit que b')
else:
    print('a est egale a b')

a est plus petit que b

3.2.8.5. Exercice#

En utilisant les instructions if, elif et else, écrire un code où on test un nombre entier N à 4 chiffre:

s’il est divisible par 2 et 3, vous affichez: Votre nombre n'est pas seulement divisible par 2 et 3 mais aussi par 6.
S’il est divisible par 2, vous affichez: Votre nombre est divisible par 2.
S’il est divisible par 3, vous affichez: Votre nombre est divisible par 3.
Dans tous les autres cas, vous affichez: Votre nombre n'est pas divisible par 2, 3 ou 6.

# écrire ici votre solution

N = 1200

if N%2 == 0 :
    if N%3 == 0 :
        print(N," est divisible par 2, 3, 6")
    else:
        print(N," est divisible par 2")
elif N%3 == 0:
    print(N," est divisible par 3")
else:
    print(N," n'est divisible par 2, 3 ou 6")

1200  est divisible par 2, 3, 6

3.2.9. Type de variables#

variable:

case mémoire pour stocker de l’information
doit être déclaré (initialisé) avant d’être utilisé
son type dépend de la valeur d’initialisation et est fonction du type de l’information (entier, réel, chaine)

la fonction python type(var) permet d’avoir le type de la variable

x=1
print(x,type(x))
x=1.0
print(x,type(x))
x='1'
print(x,type(x))

<class 'int'>
0 <class 'float'>
<class 'str'>

# erreur variable non définit
try:
    x = y
except Exception as err:
    print("Exception:",err) 

Exception: name 'y' is not defined

3.2.10. Fonction#

implémentation d’un algorithme qui traite des données pour obtenir un résultat.
définition de la fonction (aucun code n’est exécuté)
utilisation de la fonction (exécution de la fonction avec les données fournies)

3.2.10.1. définition d’une fonction#

  def  MaFonction(donnees)
        '''documentation'''
        ...
        resultat = calculer a partir des donnees
        ...
        return resultat

Attention: dans une fonction les arguments et les variables dans la fonction sont des variables locales.

3.2.10.2. utilisation d’une fonction#

  mon_res = MaFonction(mes_donnees)

Lors de l’appel d’une fonction, les arguments (mes_donnees) peuvent être des expressions, des valeurs ou des variables (si elles sont initialisées)

3.2.10.3. Exemple#

Ecrire une fonction qui calcule $n!$ pour $n\ge 0$

# definition
def factoriel(n):
    '''calcul n!  pour n>=0 '''
    fac = 1

    return fac

# validation 3! -> 6
factoriel(3)

# utilisation calcul de n!! pour n=4
n = 4
fn = factoriel(n)
fnn = factoriel(fn)
print(fnn)
factoriel(factoriel(n))

3.2.10.3.1. exécution du code avec des « print »#

# verification avec des print
a=3
c=-3
def func1(a,b):
    c=0
    print("dans func1 a={} b={} c={}".format(a,b,c))
    def func2(a,b):
        c = a - b
        print("dans func2 a={} b={} c={}".format(a,b,c))
        return c
    c=func2(b,a)
    print("fin func1 a={} b={} c={}".format(a,b,c))
    return c
a=func1(2,1)
try:
    print("a={} c={}".format(a,c))
    print("b={}".format(b))
except Exception as err:
    print("Exception:",err) 

dans func1 a=2 b=1 c=0
dans func2 a=1 b=2 c=-1
fin func1 a=2 b=1 c=-1
a=-1 c=-3
b=5

3.2.11. Liste et tableau#

liste: ensemble ordonné de valeurs
- ajout et suppression d’éléments
- la taille et le type peuvent variés
tableau: ensemble ordonné de valeurs de même type
- vecteurs, matrices
- taille fixée
indice: on compte à partir de 0
- indice à partir de 0
- [ ] pour sélectionner un élèment
- [0] premier element
- [-1] dernier element
- [n0:n1:p] selection des elements de l’indice n0 (defaut 0) à n1 (exclus) avec un pas p (defaut 1)
aliasing attention à la copie
- A = B aliasing (A et B sont identiques)
- A = B[:] copie

3.2.11.1. Test aliasing#

# question 1 uniquement
%activity /usr/local/commun/ACTIVITY/IntroPython/questionAlgorithme1

3.2.11.2. vérification#

# pble aliasing
L=[1,2,3]
L1=L
print(L,L1)
L[0]=2
print(L,L1)

[1, 2, 3] [1, 2, 3]
[2, 2, 3] [2, 2, 3]

#  manipulation de liste
L=[1,2,3]
L1=L[:]
print(L,L1)
L1[0]=2
print(L,L1)

[1, 2, 3] [1, 2, 3]
[1, 2, 3] [2, 2, 3]

3.2.12. Erreur algorithmique#

3.2.12.1. exemple 1#

Pour éffectuer une permutation circulaire à droite d’un tableau X : $X_{i+1}=X_i$ (décalage à droite), on utilise l’un des deux programmes suivants. L’un est algorithmiquement faux. Lequel ?

# question 3 uniquement
%activity /usr/local/commun/ACTIVITY/IntroPython/questionAlgorithme3

X = [1,2,3,4,5]
n = len(X)
# decalage a droite: version 1
a = X[-1]
for i in range(n-1):
    X[i+1]=X[i]
X[0]=a

Y = [1,2,3,4,5]
n = len(Y)
# decalage a droite: version 2
x=Y[-1]
for i in range(1,n):
    Y[n-i]=Y[n-i-1]
Y[0]=x

# test de validation
print(X)
print(Y)

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

3.2.12.2. exemple 2#

Calcul de l’expression suivante: itération de Gauss-Seidel

  pour i de a N

\[ Y_i = Y_i - \frac{\sum_{j=1}^{N} A_{ij} Y_j - B_i}{A_{ii}} \]

Pour calculer cette expression, on utilise l’une des deux fonctions suivantes. L’une est algorithmiquement fausse. Laquelle ?

# version 1
def iteration1(A,B,X):
    n = len(X)
    Y = X[:]
    for i in range(n):
        sum=0.0
        for j in range(n):
            sum = sum + A[i][j]*Y[j]
            Y[i] = Y[i] - (sum - B[i])/A[i][i]
    return Y

# version 2
def iteration2(A,B,X):
    n = len(X)
    Y = X[:]
    for i in range(n):
        sum=0.0
        for j in range(n):
            sum = sum + A[i][j]*Y[j]
        Y[i] = Y[i] - (sum - B[i])/A[i][i]
    return Y

# exemple de validation

X=[1,2,3]
A=[[1,2,3],[4,5,6],[7,8,9]]
B=[1,1,2]
# 
X1=iteration1(A,B,X)
print("X1=",X1)
X2=iteration2(A,B,X)
print("X2=",X2)

X1= [-16.0, 7.4, 6.088888888888887]
X2= [-12.0, 6.2, 4.044444444444444]

3.2.13. programmation récursive#

la fonction s’appelle elle-même
calcul factorielle n!
- n! = n*(n-1)!

def fac(n):
    ''' calcul recursif de n!'''
    if n>1 :
        return n*fac(n-1)
    else :
        return 1

# validation
fac(3)

3.2.13.1. example 2#

Que calcule la fonction récursive suivante ?

def fonc(L):
    print("appel fonc avec ",L)
    if not L:
        return 0
    else:
        res = 1 + fonc(L[1:])
        print("res=",res)
        return res

# resultat
fonc([1,2,3,4])

appel fonc avec  [1, 2, 3, 4]
appel fonc avec  [2, 3, 4]
appel fonc avec  [3, 4]
appel fonc avec  [4]
appel fonc avec  []
res= 1
res= 2
res= 3
res= 4

3.2.14. Erreur sous Python#

retour erreur

   Traceback (most recent call last):
    File "test.py", line 6, in <module>
    test()
    File "test.py", line 3, in test
    print table[4]
    IndexError: list index out of range

** code d’erreurs**

IndentationError :
```
  expected an indented block
```
IndexError:
```
  list index out of range
```

SyntaxError :

  inconsistent use of tabs and spaces in indentation

NameError :
```
  name 'X' is not defined
```
ImportError :
```
  no module named X
```
ZeroDivisionError :
```
  X division or modulo by zero
```

3.2.15. Programmation sctructuré#

principe: « Divide and Conquer »

3.2.15.1. analyse descendante: top-down design#

définition des différentes étapes pour résoudre le problème
on découpe le problème en une série de sous-problèmes plus simples (si possible indépendant)
on spécifie ce qui doit résolu dans chacun des sous-problèmes sans forcément dire comment
puis on itère au niveau des sous problèmes.

analyse top-down

3.2.15.2. programmation ascendante: bottom-up programming#

on programme d’abord les sous-problèmes(sous forme de fonction)
on validation les fonctions
puis on réitère en remontant dans l’arbre jusqu’au programme principal
on crée une bibliothéque
on éffectue l’analyse et la validation globale

règles réutiliser les fonctions déjà écrites et validées (bibliothèques): principe du moindre effort !

3.2.15.3. Méthodologie#

Réflexion algorithmique
Création d’une bibliothèque de fonctions pour résoudre les sous-pbles:
```
 fichier mabib.py
```
Utilisation d’un environnement de programmation (éditeur / jupyter lab/ shell)
Validation de la bibliothèque
Ecriture du programme principale en utilisant la bibliothéque
```
 import mabib.py
```
Utilisation d’un notebook jupyter
Analyse et validation du résultat (tracé de courbes)
Ecriture d’un rapport au format $\LaTeX$

3.2.16. Bibliographie#

Python. The official Python web site.
Python tutorials. Le tutoriel officiel Python.
Think Python. “”How to Think Like a Computer Scientist”” by Allen B. Downey (free book).
COURS InProS Cours en vidéo d“« INtroduction à la Programmation Scientifique””
Scientific Python Le site officiel de SciPy, qui regroupe les bibliothèques scientifiques les plus utilisées en Python: numpy, scipy, matplotlib, sympy et pandas

3.2.17. FIN de la leçon#

3.3. Tableaux avec numpy et courbes avec matplotlib#

Marc BUFFAT, dpt mécanique, Université Lyon 1 et [1]

images/numpy_array.png

from IPython.display import display, Markdown, clear_output
def printmd(string):
    display(Markdown(string))
from metakernel import register_ipython_magics
register_ipython_magics()

3.3.1. Introduction#

Pour beaucoup d’applications en mécanique, on est amené à manipuler des tableaux qui sont des séquences de données du même type. Ils ressemblent beaucoup aux listes, à l’exception de la contrainte sur le type des éléments. Cela apporte un énorme avantage en termes d’efficacité aux tableaux par rapport aux listes. Ainsi, les méthodes sur les tableaux s’exécutent beaucoup plus rapidement que celles sur les listes.

Grâce aux bibliothèques, le langage Python permet de traiter efficacement des applications numériques en sciences et technologie. Pour le calcul scientifique, la bibliothèque la plus importante est NumPy (Numerical Python), qui fournit les structures de données de tableaux (array) à n dimensions ( ndarray), avec toutes les fonctions, opérations et algorithmes pour les calculs d’algèbre linéaire. Nous utiliserons aussi la bibliothèque Matplotlib pour tracer les données en deux dimensions.

3.3.2. Importation de bibliothèques#

En Python, les bibliothèques ne sont pas chargées automatiquement, mais doivent etre importés explicitement en utilisant la commande import. Par exemple, pour importer NumPy, avec toutes ses fonctions d’algèbre linéaire, on utilise

# soit
import numpy
# ou
import numpy as np

La deuxième forme, que nous utiliserons, permet de donner un nom abrégé (ici sp) à la bibliothèque, que l’on utilise à la place du nom de la bibliothèque (i.e. numpy).

Une fois la bibliothéque chargée, on peut utiliser les fonctions de la bibliothèque en ajoutant devant le nom de la bibliothéque ou son nom abrégé. Par exemple, certaines fonctions couramment utilisées avec numpy sont:

Suivez les liens pour voir la documentation de ces fonctions NumPy très utiles, ou utiliser le menu Aide->Numpy Reference!

import numpy as np

3.3.3. Création de tableaux#

Pour créer un tableau NumPy à partir d’une liste existante de valeurs ( homogènes), on utilise numpy.array():

   numpy.array([3, 5, 8, 17])

NumPy propose de nombreuses façons de créer des tableaux (cliquez sur le lien).

Par exemple numpy.ones() et numpy.zeros() créent des tableaux remplis de uns et de zéros, respectivement. On passe en argument le nombre souhaité d’éléments du tableau.

Créer dans la cellule suivante un tableau de 1 de dimension 5 et un tableau de 0 de dimension 3.

# écrire votre code ici

# solution

Une autre fonction numpy.arange() crée un tableau de valeurs régulièrement espacées de pasdans un intervalle définit [debut,fin] mais où finest exclu.

Syntaxe:

numpy.arange(debut, fin, pas)

où debut vaut par défaut zéro,fin est exclus du tableau et pas a une valeur par défaut de 1.

Utiliser les commandes suivantes et analyser le résultat:

np.arange(4)
np.arange(2, 6)
np.arange(2, 6, 2)
np.arange(2, 6, 0.5)

# écrire votre code ici

la fonction numpy.linspace() est similaire à numpy.arange(), mais utilise un nombre d’échantillons N au lieu du pas. Il renvoie un tableau de N valeurs avec des nombres régulièrement espacés sur l’intervalle spécifié [debut,fin] où finest inclu.

Syntaxe:

numpy.linspace(debut,fin,N)

Utiliser les commandes suivantes et analyser les résultats:

  np.linspace(2.0, 3.0)
  len(np.linspace(2.0, 3.0))
  np.linspace(2.0, 3.0, 6)
  np.linspace(-1, 1, 9)

# écrire votre code ici

3.3.4. Opérations sur les tableaux#

En affectant un tableau à une variable, on peut ensuite faire simplement des opérations sur le tableau, qui correspondent à des opérations éffectuées sur chaque élément du tableau (i.e. terme à terme).

Utiliser les commandes suivantes et analyser les résultats

   X = np.linspace(-1, 1, 9)
   Y = X * X
   Z = np.sqrt(Y)
   W = X + Z

# écrire votre code ici

Nous pouvons également diviser des tableaux, mais il faut faire attention à ne pas diviser par zéro. Cette opération se traduira par un nan qui signifie Not a Number. Python effectuera toujours la division, mais nous informera du problème.

Tester le résultat de l’opération

   X / Y

# écrire votre code ici

3.3.5. Tableaux multidimensionnels#

3.3.5.1. tableaux 2D#

NumPy peut créer des tableaux de N dimensions. Par exemple, un tableau 2D ressemble à une matrice et est créé à partir d’une liste imbriquée comme suit:

  A = np.array([[1, 2], [3, 4]])

Les opérations sur matrices sont des opérations terme à terme sur chaque élément.

A partir des 2 matrices A, B suivantes, affichez le résultat de

A + B
A - 2B
A * B

questions

Que se passet-t-il si les matrices ne sont pas de la même dimension ?
Peut on utiliser une notation implicite pour les opérations p.e. 2B ?

A = np.array([[1, 2], [3, 4]])
B = np.array([[1, -1], [0, 1]])
# votre code ici

Attention la dernière opération A*B ne correspond pas à la multiplication matricielle (réfléchissez pourquoi ?)

La multiplication matricielle est obtenue avec l’opérateur @ ou la fonction numpy.dot()

    A @ B
    np.dot(A,B)

Dans la cellule suivante calculer le produit matriciel des 2 matrice A et B.

Si X est le vecteur donné ci dessous, calculer le produit matrice vecteur $A X$

Que fournit l’expression

A*X

X = np.array([1,2])

# votre code ici

3.3.5.2. algèbre linéaire#

En mécanique, on doit manipuler assez souvent des systèmes d’équations linéaires. Numpy possède une sous-bibliothèque d’algèbre linéaire linalg fournissant les fonctions suivantes pour un système linéaire $$ A X = D $$

np.linalg.det(A) calcule le déterminant de A
np.linalg.inv(A) calcule $A^{-1}$ l’inverse de A
np.linalg.solve(A,D) calcule la solution X du système linéaire $A X = D$
np.linalg.eigvals(A)calcul les valeurs propres de A
np.linalg.eig(A)calcul les valeurs et les vecteurs propres de A

remarque importante

pour résoudre un système linéaire numériquement on n’utilise jamais la méthode mathématique $ X = A^{-1} D$ en calculant $A^{-1}$ l’inverse de A, car même si on sait calculer l’inverse de A numériquement avec np.linalg.inv(A). Cette méthode est sujette à des erreurs importantes d’arrondie et est inefficace. On utilise à la place la méthode d’élimination de Gauss (que vous verrez dans le cours de méthodes numériques) qui correspond à la fonction np.linalg.solve(A,D).
pensez à consultez la documentation numpy et numpy.linalg (menu Aide->numpy)

exercise En utilisant le vecteur D suivant et la matrice A précédente, calculer la solution Y de $$ A Y = D $$ et mettre le résultat dans la variable Y

Comment peut-on vérifier le résultat ?

D = np.array([5,11])
# votre code ici
Y = 0

3.3.5.3. sélection des éléments d’un tableau#

Pour accéder aux éléments d’une matrice on utilise la notation classique [ligne,colonne].

Pour sélectionner une ligne entière on utilise : à la place de ligne, et idem pour une colonne.

Pour sélectionner une partie du tableau on utilise la notation debut:fin pour afficher les valeurs de debut à fin (attention fin est exclu)

Exercices

afficher la valeur de la première ligne, première colonne de A,
afficher la valeur de la première ligne, deuxième colonne de A,
afficher la première ligne de A,
afficher la seconde colonne de A.

# votre code ici

3.3.5.4. boucle sur les éléments d’un tableau#

la bibliothèque numpy possède de nombreuses fonctions pour manipuler des tableaux. par contre dans certain cas on doit effectuer des calculs particuliers qui nécessitent de faire des boucles itératives sur les éléments d’un tableau.

Python offre divers formes de boucles itératives à utiliser suivant les besoins. Nous verrons ici uniquement des boucles sur des tableaux 1D ou vecteurs:

3.3.5.4.1. boucle sur les valeurs#

soit un tableau Tab, on s’intéresse à la valeur val des éléments successifs du tableau? C’est la boucle classique en Python:

   for val in Tab:
      print(val)

exercise on veut déterminer le nombre de valeurs d’un tableau qui sont supérieures à une valeur donnée. Par exemple dans le tableau X suivant, on veut connaître le nombre de valeurs supérieure à 0.5

X = np.random.rand(5)
print(X)
# votre code ici

[0.33661121 0.75305157 0.6049483  0.00343809 0.25299808]

3.3.5.4.2. boucle sur les indices et les valeurs#

Si l’on souhaite modifier les valeurs du tableau, la boucle précédente ne le permet. Il faut avoir accés aussi à l’indice i de la valeur val dans le tableau Tab.

En python la boucle s’écrit alors:

   for i,val in enumerate(Tab):
      print(val,"se trouve a l'indice ",i)

exercise en reprenant l’exemple précédent, on souhaite écrêter le tableau X en remplaçant les valeurs supérieures à 0.5 par 0.5.

# votre code ici

3.3.5.4.3. boucle sur les indices uniquement#

enfin on peut vouloir faire une boucle uniquement sur les indices, dans le cas ou manipule plusieurs tableaux.

En python la boucle s’écrit alors:

   for i in range(len(Tab)):
      print("pour l'indice ",i," la valeur vaut ",Tab[i])

On a besoin de la taille du tableau Tab que l’on peut obtenir par la fonction générique len() ou la méthode .size() pour un tableau ndarray

exercise en reprenant l’exemple précédent, on souhaite écrêter le tableau X en remplaçant les valeurs supérieures à 0.5 par la valeur d’un autre tableau Y de même dimension que X

Y = 0.5*np.random.rand(5)
# votre code ici

3.3.5.5. Aliasing#

lors d’un instruction d’affectation

    X = Y

si X,Y sont des scalaires pas d’aliasing
si X,Y sont des listes ou des tableaux , alors aliasing
utilisation d’une copie explicite

print("pas d'aliasing avec les scalaires")
x=2
y=x
x=3
print(x,"#",y)

pas d'aliasing avec les scalaires
3 # 2

print("mais aliasing avec les listes et tableaux")
X=[1,2,3]
Y=X
X[0]=0
print(X,"==",Y)

mais aliasing avec les listes et tableaux
[0, 2, 3] == [0, 2, 3]

print("pas d'aliasing si copie explicite des valeurs")
X=np.array([1,2,3])
Y=X.copy()
X[0]=0
print(X,Y)
Y=np.zeros(X.shape)
Y[:]=X
X[0]=3
print(X,Y)

pas d'aliasing si copie explicite des valeurs
[0 2 3] [1 2 3]
[3 2 3] [0. 2. 3.]

3.3.5.5.1. Question test#

%activity /usr/local/commun/ACTIVITY/IntroPython/questionAlgorithme1

3.3.5.6. Notation matricielle#

Attention par defaut opération sur les éléments
les fonctions mathématiques dans numpy ont pour aguments des tableaux

np.dot produit scalaire
np.outer produit tensorielle
np.inner

X=np.arange(1,4,1)
print("X=",X)
print("X*X=",X*X)
print("X.X=",np.dot(X,X))
print("XoX=",np.outer(X,X))
print("XiX=",np.inner(X,X))

X= [1 2 3]
X*X= [1 4 9]
X.X= 14
XoX= [[1 2 3]
 [2 4 6]
 [3 6 9]]
XiX= 14

A=np.array(np.outer(np.arange(3),np.arange(3)))
print(type(A),A.shape)
print("A*X=\n",A*X)
print("A.X=\n",np.dot(A,X))
print("AoX=\n",np.outer(A,X))
print("AiA=\n",np.inner(A,A))
print("A.At=\n",np.dot(A,A.transpose()))

<class 'numpy.ndarray'> (3, 3)
A*X=
 [[ 0  0  0]
 [ 0  2  6]
 [ 0  4 12]]
A.X=
 [ 0  8 16]
AoX=
 [[ 0  0  0]
 [ 0  0  0]
 [ 0  0  0]
 [ 0  0  0]
 [ 1  2  3]
 [ 2  4  6]
 [ 0  0  0]
 [ 2  4  6]
 [ 4  8 12]]
AiA=
 [[ 0  0  0]
 [ 0  5 10]
 [ 0 10 20]]
A.At=
 [[ 0  0  0]
 [ 0  5 10]
 [ 0 10 20]]

3.3.5.7. attention:#

on évite les boucles explicites pour des questions d’efficacité
mais on commence toujours par écrire la forme explicite avec des boucles , puis ensuite la forme vectorielle

exemple: calcul de la forme quadratique $$ \Phi = \sum_{i=1}^n \sum_{j=1}^n A_{ij} X_i X_j - \sum_{j=1}^n B_j X_j$$

# calcul explicite
Phi=0.0
n=X.size
B=np.ones(X.shape)
### BEGIN SOLUTION
for i in range(n):
    sum=0.0
    for j in range(n):
        sum = sum + A[i,j]*X[j]
    Phi = Phi + (sum - B[i])*X[i]
print(Phi)
### END SOLUTION

58.0

# fonction numpy
print(np.dot(np.dot(A,X)-B,X))

58.0

3.3.5.7.1. Question test#

%activity /usr/local/commun/ACTIVITY/IntroPython/questionAlgorithme2

%activity /usr/local/commun/ACTIVITY/IntroPython/questionAlgorithme3

3.3.6. Tracé avec matplotlib#

Pour tracer des courbes, nous utiliserons la bibliothèque matplotlib et son module pyplot. Sous Jupyter on utilise les instructions suivantes:

  %matplotlib inline
  import matplotlib.pyplot as plt

pour utiliser cette librairie avec le nom raccourci plt. La première instruction permet d’afficher les courbes dans le notebook (inline) plutôt que dans une fenêtre à part.

exemple

pour $x\in[0,2]$, on veut tracer les courbes $y_1=x^2$, $y_2=x^3$ et $y_3=\sqrt{x}$. Pour cela on crée un tableau X de valeurs entre 0 et 2, et on calcule ensuite les valeurs des courbes, puis on trace les 3 courbes avec des commandes du type

  plt.plot(X,Y1)

%matplotlib inline
import matplotlib.pyplot as plt

X  = np.linspace(0,2,41)
Y1 = X*X
Y2 = X**3
Y3 = np.sqrt(X)

3.3.6.1. tracé basique#

faire un tracé basique des 3 courbes dans la cellule suivante.

Attention ce genre de tracé n’est à utiliser que pour des tests et pas pour présenter des résultats, car les courbes n’ont pas de titre, les axes pas de label. Il faut regarder le code pour éventuellement comprendre.

N’utilisez jamais ce tracé basique pour présenter vos résultats !!!

# tracer basique

3.3.6.2. tracé de résultats#

Pour avoir un résultat de qualité, il faut ajouter un titre lisible, des labels avec les fonctions:

plt.title
plt.xlabel
plt.figure
plt.legend

# tracer avec des titres et des labels
# utilisation de taille de caractères plus grande
plt.rc('font', family='serif', size='18')
# tracer d'une figure plus grande
plt.figure(figsize=(10,8))
# plot x^2
plt.plot(X, Y1, color='r', linestyle='-', label='$x^2$')
# plot x^3
plt.plot(X, Y2, color='g', linestyle='--', label='$x^3$')
# plot sqrt(x)
plt.plot(X, Y3, color='b', linestyle=':',  label='$\sqrt{x}$')
# ajoute un titre et des labels 
plt.title('Courbes étudiées')
plt.xlabel('x')
plt.ylabel('y')
# ajoute la legende au meilleur endroit
plt.legend(loc='best');

_images/a241d50206d5af0023ff7888c90751a9a7e52b64bb4b9e03b74d329bcacf7fff.png

On a utiliser une notation $\LaTeX$ pour afficher des expressions mathématiques: $x^2$ pour afficher $x^2$ ou $\sqrt{x}$ pour afficher $\sqrt{x}$.

Consulter le site de matplotlib https://matplotlib.org pour avoir un aperçu des immenses possibilités de tracé de matplotlib.

3.3.6.3. tracé avec des sous-figures#

# tracer de plusieurs sous figure
fig = plt.figure(figsize=(10,8))
# plot x^2
plt.subplot(3,1,1)
plt.plot(X, Y1, color='r', linestyle='-', label='$x^2$')
plt.title('Courbes étudiées')
plt.ylabel('$y_1$ $[m]$')
plt.legend(loc='best')
# plot x^3
plt.subplot(3,1,2)
plt.plot(X, Y2, color='g', linestyle='--', label='$x^3$')
plt.ylabel('$y_2$ $[m/s]$')
plt.legend(loc='best')
# plot sqrt(x)
plt.subplot(3,1,3)
plt.plot(X, Y3, color='b', linestyle=':',  label='$\sqrt{x}$')
# ajoute un titre et des labels 
plt.ylabel('$y_3$ $[m/s^2]$')
plt.xlabel('x [m]')
plt.legend(loc='best');

_images/793cf9bbf871d8bb56f5b3d9c286bccc9af736c246e8eb6a0ac8c1febabff127.png

3.3.6.4. mise sur fichier#

fonction savefig()

attention: sauvegarde la figure dans une variable fig

fig.savefig(nom_fichier)

format pdf (vectoriel) ou image (png)

# sauvegarde de la figure
fig.savefig("mafig.png")
! ls -al mafig.png

-rw-rw-r-- 1 buffat buffat 63616 avril  9 13:57 mafig.png

3.3.6.5. Gallerie matplotlib#

from IPython.display import HTML
HTML('<iframe src=https://matplotlib.org/gallery.html width=800 height=500></iframe>')

/home/buffat/venvs/jupyter/lib/python3.10/site-packages/IPython/core/display.py:475: UserWarning: Consider using IPython.display.IFrame instead
  warnings.warn("Consider using IPython.display.IFrame instead")

3.3.7. Méthodologie#

Réflexion algorithmique
Création d’une bibliothèque de fonctions pour résoudre les sous-pbles:
```
 fichier mabib.py
```
Utilisation d’un environnement de programmation (éditeur / jupyter lab/ shell)
Validation de la bibliothèque
Ecriture du programme principale en utilisant la bibliothéque
```
 import mabib.py
```
Utilisation d’un notebook jupyter
Analyse et validation du résultat (tracé de courbes)
Ecriture d’un rapport au format $\LaTeX$

3.3.8. Bibliographie#

Python. The official Python web site.
Python tutorials. Le tutoriel officiel Python.
Think Python. “”How to Think Like a Computer Scientist”” by Allen B. Downey (free book).
COURS InProS Cours en vidéo d“« INtroduction à la Programmation Scientifique””
Scientific Python Le site officiel de SciPy, qui regroupe les bibliothèques scientifiques les plus utilisées en Python: numpy, scipy, matplotlib, sympy et pandas

3.3.9. FIN de la leçon#

3.4. Vidéo d’introduction à Python#

Il s’agit de la première leçon d’un ensemble de modules d’apprentissage pour les étudiants en sciences et technologie. Ces modules utilisent le langage Python, mais en présupposant aucune expérience de programmation préalable. Le premier objectif sera de vous apprendre à utiliser un environnement numérique Python avec des notebooks Jupyter permettant de résoudre des problèmes et de gérer des données scientifiques.

vidéo du cours introduction à Python

3.5. Vidéo sur les notebooks Jupyter#

Cette cidéo décrit l’utilisation des notebooks Jupyter pour un cours ou TP sur les plateformes jupyter nbgrader du département mécanique.

Une erreur s’est glissée dans la vidéo du cours: seriez vous capable de la trouver?

vidéo du cours sur les notebook jupyter

3.6. Bibliographie#

Engineering Computations cours du Pr L. Barba (Washington Univ.) qui a inspiré ce cours
Python. The official Python web site.
Python tutorials. Le tutoriel officiel Python.
Think Python. “”How to Think Like a Computer Scientist”” by Allen B. Downey (free book).
COURS InProS Cours en vidéo d’INtroduction à la Programmation Scientifique
Scientific Python Le site officiel de SciPy, qui regroupe les bibliothèques scientifiques les plus utilisées en Python: numpy, scipy, matplotlib, sympy et pandas

et pour les plus curieux d’entre vous

Algorithms in Python implementation en Python de nombreux algorithmes
Tensor Flow in Python pour utiliser la bibliothèque d’IA de Google avec Python

1. Introduction à Unix#

1.1. Notions de système d’exploitation: « OS »#

1.1.1. Principes#

types: Unix (Linux), Windows (Win 95/98, Win NT, XP), Mac OS …

BIOS
Chargeur=Boot (LILO)
Noyau=Kernel (vmlinuz)
Interpréteurs de commande=shell (sh, bash, csh)
Interfaces graphiques=GUI (X11, gestionnaires de fenetres KDE, NexStep)

1.2. Système de fichiers#

arborescence logique: répertoires (dossier, directory)

racine ou root « / »
répertoires systèmes « /bin » « /lib » « /usr »
répertoires utilisateurs « /home/nom_utilisateur »

arborescence physique: disque et partitions

partitions locales et distantes\

[buffat@ufrmeca Cours]$ df 
Filesystem           1k-blocks Use% Mounted on
/dev/sda2                54441 55% /                      
/dev/sda6               995115 40% /home
/dev/sda7               792800 87% /home1
/dev/sda5               995115 80% /usr
/dev/sda8               567548 75% /home2
/dev/sdb1              1981000 47% /home3
/dev/sdb2              1981000 49% /home4
/dev/sdb3              1981000  1% /home5  
/dev/sdb4              2687444  0% /home6
mecapar:/home7         2640869 36% /mnt

1.3. Controle d’accès#

login + password

permissions
protection des fichiers

[buffat@ufrmeca cours_unix]$ ls -al base.tex
-rw-r--r--   1 buffat   aero        28044 Jan 29  1995 base.tex                

et des répertoires

[buffat@ufrmeca Latex]$ ls -al cours_unix
total 732
drwxr-xr-x   8 buffat   aero         1024 Sep 10 13:56 .
drwxr-xr-x  24 buffat   aero         1024 Jul 28 16:00 ..
drwxr-xr-x   2 buffat   aero         1024 May 18 13:21 BE
drwxr-xr-x   3 buffat   aero         1024 Apr 20 09:42 Cours
drwxr-xr-x   6 buffat   aero         1024 Feb 17  1997 Demo
-rw-r--r--   1 buffat   aero        28044 Jan 29  1995 base.tex  

Modification des protections d’un fichier: chmod.

chmod {a,u,g,o}{+,-}{r,w,x}

avec all (tous), user (propriétaire), group (groupe) , ou other (les autres).

1.3.1. liens sur les fichiers#

commande ln

   ln {-as}

   ln -s ../Cours/doc.tex doc1.tex

1.4. Commandes du shell#

cd Change le répertoire de travail courant.
Syntaxe: cd \
ls Affiche des informations sur les fichiers ou répertoires.
Syntaxe: ls …
Exemple : ls -lF /usr/bin affichera le contenu du répertoire /usr/bin.\
cp Copie un (des) fichier(s) dans un autre fichier ou répertoire.
Syntaxe: cp …
Exemple : cp ../frog joe copie le fichier ../frog dans le fichier ou le répertoire joe.\
mv Déplace un (des) fichier(s) vers un autre fichier ou répertoire.
Syntaxe: mv …
Exemple : mv ../frog joe déplace le fichier ../frog dans le fichier ou le répertoire joe.\
rm Supprime des fichiers..
Syntaxe: rm …
Options: -i demandera confirmation avant chaque effacement de fichier.
Exemple : rm -i /home/dupont/recettes/cassoulet /home/dupont/recettes/vomitifs/macdo effacera les fichiers cassoulet et macdo si l’opération est bien confirmée par l’opérateur.\
mkdir Crée de nouveaux répertoires.
Syntaxe: mkdir …
Exemple : mkdir /tmp/test créera le répertoire test dans le répertoire /tmp.\
rmdir Cette commande supprime les répertoires vides.
Syntaxe: rmdir …
Où à sont les répertoires à supprimer.
Exemple : rmdir /tmp/test supprime le répertoire test dans /tmp, si il est vide (et s’il existe).\
man Affiche la page de manuel pour la commande ou la ressource donnée.
Syntaxe: man
Exemple : man ls donne la description de la commande ls.\
more Affiche le contenu des fichiers, un écran à la fois.
Syntaxe: more …
Exemple : more /etc/termcap affiche le fichier
/etc/termcap.
cat Normalement destinée à concaténer des fichiers, la commande cat est aussi utilisée pour afficher tout le contenu d’un fichier d’un coup.
Syntaxe: cat …
Exemple : cat /etc/passwd affiche le contenu du fichier /etc/passwd.\
echo Affiche simplement les arguments qu’on lui passe.
Syntaxe: echo …
Example: echo "Bonjour tout le monde" affichera la chaı̂ne « Bonjour tout le monde ».\
grep Affiche toutes les lignes dans le (les) fichier(s) correspondant à l’expression donnée.
Syntaxe: grep …
Où est une expression rationnelle, et à les fichiers dans lesquels la rechercher.
Exemple : grep local /etc/hosts affichera toutes les lignes du fichier /etc/hosts qui contiennent l’expression « local ».\

1.5. Caractères génériques du shell#

* remplace n’importe quelle suite de caractère
? remplace 1 caractère
[ccc] remplace un caractère parmi ccc; (0-9) et (a-b) sont possibles
x protection du caractère x (i.e non interprétation)
“xxx” prend xxx littéralement sans interprétation
"xxx" idem mais avec interprétation de
‘cde‘ remplace par le résultat de la commande
# début de commentaire
; fin de commande

1.6. les entrées sorties#

entrée standard stdin, sortie standard stdout, sortie d’erreur stderr
redirection des E/S

: *prog file* remplace l’entrée standard de prog par le contenu de file.
str : prog `` str l’entrée standard suit jusqu’à l’apparition de « str » dans le flot de données.
: *prog file* écrit la sortie standard de prog dans file.
: *prog file* ajoute la sortie standard de prog dans file.
| : prog1 | prog2 création d’un « pipe » en connectant la sortie de prog1 à l’entrée de prog2

1.7. contrôle de processus#

commande ps

[buffat@ufrmeca COURS_UNIX]$ ps
  PID TTY          TIME CMD
 6289 pts/1    00:00:00 bash
11272 pts/1    00:00:00 ps

processus en avant-plan, arrière-plan
arrêt d’un processus
kill pid

1.8. Editeur de texte#

un éditeur simple, mais puissant vi ou autre kwrite,.. xemacs..

2. Réseau#

2.1. Définition#

nom de machine hostname
nom de domaine domainname\

[buffat@ufrmeca COURS_UNIX]$ hostname -a
ufrmeca
[buffat@ufrmeca COURS_UNIX]$ hostname -d
univ-lyon1.fr  
[buffat@ufrmeca COURS_UNIX]$ hostname
ufrmeca.univ-lyon1.fr  

numéro IP = numéro carte ethernet unique

[buffat@ufrmeca COURS_UNIX]$ /sbin/ifconfig
eth0      Lien encap:Ethernet  HWaddr 00:C0:4F:C9:62:03
          inet adr:134.214.93.120  Bcast:134.214.95.255  Masque:255.255.252.0
          UP BROADCAST RUNNING MULTICAST  MTU:1500  Metric:1
          Paquets Reçus:5443720 erreurs:13 jetés:13 débordements:0 trames:13
          Paquets transmis:697928 erreurs:0 jetés:0 débordements:0 carrier:0
          collisions:6850 lg file transmission:100
          Interruption:19 Adresse de base:0xdce0
 
lo        Lien encap:Boucle locale
          inet adr:127.0.0.1  Masque:255.0.0.0
          UP LOOPBACK RUNNING  MTU:3924  Metric:1
          Paquets Reçus:4636 erreurs:0 jetés:0 débordements:0 trames:0
          Paquets transmis:4636 erreurs:0 jetés:0 débordements:0 carrier:0
          collisions:0 lg file transmission:0                                   

Numero de réseau et masque

[buffat@ufrmeca COURS_UNIX]$ cat ifcfg-eth0
DEVICE=eth0
IPADDR=134.214.93.120
NETMASK=255.255.252.0
NETWORK=134.214.92.0

2.2. Protocol de communication#

transmission de packets sur le réseau avec un protocol de communication (TCP/IP)
accès à un service entre 2 machines

| num. IP destination | num. IP source | service (port) |  message |\\

2.3. Serveur de noms#

transformation $IP <-> hostname$
fichier local $/etc/host$

[buffat@ufrmeca Cours]$ cat /etc/hosts | more
127.0.0.1       localhost
134.214.94.123  alien.univ-lyon1.fr     alien
134.214.93.120  ufrmeca.univ-lyon1.fr   ufrmeca
# UCB
#134.214.92.151 biomeca
134.214.92.164  mecaflu4.univ-lyon1.fr  mecaflu4
134.214.94.124  mecapar.univ-lyon1.fr   mecapar
134.214.92.100  secret-meca.univ-lyon1.fr       secret-meca 

serveur de noms

[buffat@ufrmeca Cours]$ cat /etc/resolv.conf
search univ-lyon1.fr
nameserver 134.214.100.6
nameserver 134.214.100.245
nameserver 156.18.22.3

resolution des noms

[buffat@ufrmeca Cours]$ nslookup
Default Server:  cismsun.univ-lyon1.fr
Address:  134.214.100.6
 
> athena.mecaflu.ec-lyon.fr
Server:  cismsun.univ-lyon1.fr
Address:  134.214.100.6
 
Non-authoritative answer:
Name:    athena.mecaflu.ec-lyon.fr
Address:  156.18.40.105        

2.4. Routage#

Porte=GATEWAY

[buffat@ufrmeca Cours]$ cat /etc/sysconfig/network
NETWORKING=yes
HOSTNAME=ufrmeca.univ-lyon1.fr
DOMAINNAME=univ-lyon1.fr
GATEWAY=134.214.92.1
GATEWAYDEV=eth0
FORWARD_IPV4=no     

Table de routage

[buffat@ufrmeca Cours]$ /sbin/route -n
Table de routage IP du noyau
Destination     Passerelle      Genmask         Indic Metric Ref    Use Iface
134.214.93.120  0.0.0.0         255.255.255.255 UH    0      0        0 eth0
134.214.92.0    0.0.0.0         255.255.252.0   U     0      0        0 eth0
127.0.0.0       0.0.0.0         255.0.0.0       U     0      0        0 lo
0.0.0.0         134.214.92.1    0.0.0.0         UG    0      0        0 eth0   

Test de routage

[root@ufrmeca Cours]# traceroute athena.mecaflu.ec-lyon.fr
traceroute to athena.mecaflu.ec-lyon.fr (156.18.40.105), 30 hops max, 40 byte packets
accel101-1ercycle (134.214.92.1)  3.336 ms  18.720 ms  2.632 ms
cs7505-cism (134.214.200.225)  1.202 ms  1.156 ms  0.964 ms
u-1-cism-villeurbanne (193.48.222.2)  2.021 ms  5.537 ms  2.071 ms
lyon.aramis.ft.net (193.48.66.13)  21.626 ms  9.631 ms  8.879 ms
centrale-lyon.aramis.ft.net (193.48.66.62)  6.949 ms  29.308 ms  60.353 ms
193.54.216.1 (193.54.216.1)  19.696 ms  47.226 ms  22.045 ms
eclgw.servers.ec-lyon.fr (156.18.19.254)  9.213 ms  110.096 ms  36.818 ms
athena.mecaflu.ec-lyon.fr (156.18.40.105)  41.103 ms  39.818 ms  80.468 ms 

Test d’une connection: $ping$

[root@ufrmeca Cours]# ping athena.mecaflu.ec-lyon.fr
PING athena.mecaflu.ec-lyon.fr (156.18.40.105): 56 data bytes
64 bytes from 156.18.40.105: icmp_seq=0 ttl=248 time=51.9 ms
64 bytes from 156.18.40.105: icmp_seq=1 ttl=248 time=53.0 ms
64 bytes from 156.18.40.105: icmp_seq=2 ttl=248 time=79.7 ms
64 bytes from 156.18.40.105: icmp_seq=3 ttl=248 time=16.5 ms
 
--- athena.mecaflu.ec-lyon.fr ping statistics ---
4 packets transmitted, 4 packets received, 0% packet loss
round-trip min/avg/max = 16.5/50.2/79.7 ms            

2.5. Services#

service= port + protocol (TCP,UDP)// messagerie (mail)= port 25 TCP

[buffat@ufrmeca Cours]$ telnet mecalpha 25
Trying 134.214.95.86...
Connected to mecalpha.univ-lyon1.fr.
Escape character is '^]'.
220 mecalpha.univ-lyon1.fr ESMTP Sendmail 8.8.7/8.8.7; Wed, 22 Sep 1999 01:00:39 +0200
helo mecalpha
250 mecalpha.univ-lyon1.fr Hello buffat@ufrmeca [134.214.93.120], pleased to meet you
......  

ssh nom@hostname\

connexion sécurisée
[buffat@p2chpd-visu2]$ ssh buffat@athena.mecaflu.ec-lyon.fr
buffat@athena.mecaflu.ec-lyon.fr's password:xxxxxxx
Last login: Thu Jul 10 09:46:01 from lns-th2-10-82-64
Digital UNIX V4.0D  (Rev. 878); Mon Dec 13 15:29:15 MET 1999
athena.mecaflu.ec-lyon.fr>exit
athena.mecaflu.ec-lyon.fr> logout
Connection to athena.mecaflu.ec-lyon.fr closed.

telnet « hostname »
ancien mode de connexion (non sécurisé)
scp user@host1:fichier1 fichier2

scp fichier1 user@host1:fichier2 transfert de fichier sécurisé

scp lstree.txt buffat@athena.mecaflu.ec-lyon.fr:toto.txt
buffat@athena.mecaflu.ec-lyon.fr's password: xxxxxx
lstree.txt       100%  206   962.5KB/s   00:00

ftp « hostname »
transfert de fichiers (non scurisé)

mail : messagerie
envoi (en général avec un client)

[buffat@ufrmeca COURS_UNIX]$ mail buffat@mecapar.univ-lyon1.fr
Subject: essai mail
Hello,
essai de mail
Marc
^D
Cc: buffat 

reception

You have mail in /var/spool/mail/buffat
[buffat@ufrmeca COURS_UNIX]$ mail
Mail version 8.1 6/6/93.  Type ? for help.
"/var/spool/mail/buffat": 1 message 1 new
>N  1 buffat@ufrmeca.univ-  Wed Sep 15 14:12  16/439   "essai mail"
& q

3. Cours algorithmique#

Marc Buffat dpt mécanique, UCB Lyon 1

%matplotlib inline
import numpy as np
import matplotlib.pyplot as plt
from IPython.display import HTML

3.1. Algorithmes numériques de base#

Un algorithme est une suite finie et non ambigüe d’opérations ou d’instructions permettant de résoudre un problème. Les algorithmes sont connus depuis l’antiquité (Euclide).

Le mot algorithme vient du nom du mathématicien perse du 9ième siècle

Abu Abdullah Muhammad ibn Musa al-Khwarizmi.

L’algorithmique correspond à la phase préparatoire avant une quelconque programmation. Elle permet de décrire un problème sous une forme que l’on peut ensuite programmer sur un ordinateur et ceci dans un langage naturel, indépendant d’un langage de programmation.

algorithme numérique suite finie et non ambiguë d’opérations ou d’instructions sur des nombres permettant de résoudre un problème.

Et il n’est pas nécessaire d’avoir un ordinateur pour exécuter un algorithme!

3.2. Méthode d’analyse algorithmique#

programmation procédurale en informatique, la programmation procédurale est un paradigme qui se fonde sur le concept d’appel procédural.

Une procédure, aussi appelée routine, sous-routine ou fonction contient simplement une série d’étapes à réaliser en fonction de données (arguments) et fournit des résultats (valeur).

top down algorithm / design

on découpe le problème en sous problèmes plus simples

bottom up programming

on programme d’abord les algorithmes des sous-problèmes, avant de passer à la résolution globale:

 - utilisation de librairies
 - définition des fonctions

3.3. Différents types d’algorithme#

algorithme itératif

Un algorithme itératif est basé sur les procédures d’itération que sont le Tant que et le Pour. Il réalise donc une boucle jusqu’à ce que la condition d’arrêt soit respectée (test , nbre d’itérations, ..).

algorithme récursif

récursion méthode algorithmique où la solution du problème s’exprime en fonction de solutions du même problème, mais sur des données plus petites.

3.3.1. Exemple : calcul de n!#

version itérative

\[ n! = \prod_{i=1}^n i \]
version récursive

\[ n! = n * (n-1)! \]

programmation dans le notebook !

3.4. Algorithme d’Euclide#

Soit une pièce rectangulaire de taille $a \times b$, déterminer la taille du plus grand carré permettant un pavage exacte de la pièce:

Analyse

propriétés
1. pavage(a,b) = pavage(a-b,b) si a>b
2. pavage(a,b) = pavage(a,b-a) si a<b
3. pavage(a,b) = a si a=b
algorithmes
1. version récurrente
2. version itérative

def PGCD(a,b):
    if a==b :
        return a
    elif a>b :
        return PGCD(a-b,b)
    else :
        return PGCD(a,b-a)

%timeit PGCD(100,75)

380 ns ± 1.87 ns per loop (mean ± std. dev. of 7 runs, 1,000,000 loops each)

def PGCD1(a,b):
    while a!=b :
        if a>b :
            a = a-b
        else:
            b = b-a
    return a

timeit PGCD1(100,75)

244 ns ± 1.76 ns per loop (mean ± std. dev. of 7 runs, 1,000,000 loops each)

3.5. Algorithme du calcul du déterminant d’une matrice#

soit $A$ une matrice d’ordre $n$, la méthode de Cramer permet le calcul récursif du déterminant par développement par rapport à la 1ere colonne

\[ det(A) = \sum_{k=1}^n (-1)^{k+1} A_{k,1} \times det(A^k) \]

où $A^k$ est la sous-matrice obtenu par élimination de la colonne 1 et de la ligne $k$ de $A$

Pour une matrice de dimension 1 $$ det(a) = a $$

### Version recursive
def determinant(A):
    
    return

3.5.1. Algorithme déterminant#

 Algorithme determinant(A)
     n = dim(A)
     si n==1 alors
         retour A[0,0]
     fin si
     # A1 sous matrice d'ordre n-1
     A1 = tableau(n-1,n-1)
     det = 0
     signe = 1
     # développement par rapport a la 1ere colonne
     pour k de 0 a n-1
         # sous matrice Ak
         A1[0:k,:] = A[0:k ,1:n]
         A1[k:,:]  = A[k+1:,1:n]
         det = det + signe*A[k,0]*determinant(A1)
         signe = -signe
         fin pour
     retour det

3.5.2. Programme Python (solution)#

from numpy.linalg import det

def determinant(A):
    """ calcul du déterminant d'une matrice A """
    n=A.shape[0]
    if n==1 :
        return A[0,0]
    A1 = np.zeros((n-1,n-1))
    det = 0.
    signe = 1
    for k in range(n):
        A1[0:k,:] = A[0:k  ,1:n]
        A1[k: ,:] = A[k+1:n,1:n]
        det = det + signe*A[k,0]*determinant(A1)
        signe = -signe
    return det
# verification
M = np.array([[2.,0.,0],[0.,3.,0],[0,0,4.]])
delta = determinant(M)
print("determinant :",delta, det(M))

determinant : 24.0 23.999999999999993

# comparaison avec scipy
M = np.random.rand(8,8)
print(M)
delta = determinant(M)
print("determinant M:",delta,det(M))

[[0.7393828  0.2230019  0.8562197  0.99433771 0.1565653  0.79930796
  0.90192473 0.3823761 ]
 [0.94994486 0.89859604 0.43061165 0.04686198 0.14548221 0.53549263
  0.49002993 0.41966703]
 [0.26416228 0.37261767 0.18743199 0.11021367 0.04673572 0.14852054
  0.10712103 0.53171711]
 [0.88452029 0.74424527 0.89223583 0.56639971 0.92099731 0.70744181
  0.8429787  0.98633387]
 [0.9449187  0.4625425  0.68669302 0.57388277 0.26252077 0.95048458
  0.18627193 0.81260822]
 [0.16122013 0.81540555 0.39614714 0.26859707 0.55721694 0.1317386
  0.10457908 0.4932543 ]
 [0.44160171 0.34201614 0.59957683 0.43623267 0.13137884 0.34305066
  0.68123181 0.62833921]
 [0.11252061 0.33153549 0.69243258 0.47728407 0.1615201  0.80683352
  0.64448575 0.19659104]]
determinant M: 0.00944264046079138 0.009442640460791378

## Efficacité
%time determinant(M)

CPU times: user 142 ms, sys: 0 ns, total: 142 ms
Wall time: 142 ms

0.00944264046079138

%time det(M)

CPU times: user 42 μs, sys: 0 ns, total: 42 μs
Wall time: 44.1 μs

0.009442640460791378

3.6. Algorithmes classiques#

méthode de Newton
méthode des trapèzes
méthode d’Euler
tour de Hanoi
tri par insertion
bubble sort

3.7. Bibliothèques scientifiques#

numpy
scipy

3.8. FIN#

4. Stucture de données#

Marc BUFFAT, dpt mécanique, Université Claude Bernard Lyon 1

Mise à disposition selon les termes de la Licence Creative Commons
Attribution - Pas d’Utilisation Commerciale - Partage dans les Mêmes Conditions 2.0 France.

Table des matières

1 Structure simple
2 Exemple
3 Structure complexe
- 3.1 Enregistrement (structure)
  - 3.1.1 structure avec définition explicite
  - 3.1.2 structure dynamique
4 TP: système solaire
5 Fin

%matplotlib inline
import numpy as np
import scipy as sp
import matplotlib.pyplot as plt
from matplotlib import rcParams
rcParams['font.family'] = 'serif'
rcParams['font.size'] = 14
from IPython.core.display import HTML
from IPython.display import display
from matplotlib import animation

4.1. Structure simple#

4.1.1. Vecteurs, Matrices#

tableau numpy

accès avec un indice
dimension fixe
type uniforme

4.1.2. Liste d’objets non uniforme#

liste python

accès avec un indice
dimension variable
type différent

4.1.3. Dictionnaire#

liste avec un indicage par mot clés - accès avec une clé (key) - dimension variable - type différent

4.1.4. Boucle sur les élèments#

par indice
par valeur
par indice et valeur

# boucle classique sur les indices
X = np.linspace(0.,1.,6)
for i in range(X.size):
    print("X[{}]={}".format(i,X[i]))

X[0]=0.0
X[1]=0.2
X[2]=0.4
X[3]=0.6000000000000001
X[4]=0.8
X[5]=1.0

# boucle sur les valeurs
for x in X:
    print("X[]=",x)

X[]= 0.0
X[]= 0.2
X[]= 0.4
X[]= 0.6000000000000001
X[]= 0.8
X[]= 1.0

# boucle sur les valeurs et indices
for i,x in enumerate(X):
    print("X[{}]={}".format(i,x))

X[0]=0.0
X[1]=0.2
X[2]=0.4
X[3]=0.6000000000000001
X[4]=0.8
X[5]=1.0

4.1.5. boucle sur un dictionnaire#

clé : valeur
dico = { cle1:valeur1, cle2:valeur2, .. }

# dictionnaire
Dico={'mon':'my','personne':'nobody','nom':'name','est':'is'}
for mot in Dico:
    print("traduction de {} : {}".format(mot,Dico[mot]))

traduction de mon : my
traduction de personne : nobody
traduction de nom : name
traduction de est : is

# utilisation
phrase="mon nom est toto"
traduction=""
for mot in phrase.split():
    if mot in Dico:
        traduction += Dico[mot] + " "
    else:
        traduction += mot + " "
print(traduction)

my name is toto

4.2. Exemple#

on souhaite manipuler une liste d’étudiants avec leur nom (une chaine de caractère) et leur note (un nombre réel)

# version numpy
Noms  = np.array(['toto','bidule','machin'])
Notes = np.array([10.,16.0,13.])
print(Noms,Notes)
# recherche de la note d'un etudiant
etudiant='bidule'
for k in range(Noms.size):
    nom = Noms[k]
    if etudiant == nom :
        print("{} a pour note {}".format(nom,Notes[k]))

['toto' 'bidule' 'machin'] [10. 16. 13.]
bidule a pour note 16.0

# version avec liste
Listes = [['toto',10.],['bidule',16.],['machin',13.]]
print(Listes)
# recherche de la note d'un etudiant
nom='bidule'
for etudiant in Listes:
    if etudiant[0] == nom:
        print("{} a pour note {}".format(nom,etudiant[1]))

[['toto', 10.0], ['bidule', 16.0], ['machin', 13.0]]
bidule a pour note 16.0

# version avec dictionnaire
Etudiants = {'toto':10.,'bidule':16.,'machin':13.}
print(Etudiants)
# recherche de la note d'un etudiant
if nom in Etudiants:
    print("{} a pour note {}".format(nom,Etudiants[nom]))

{'toto': 10.0, 'bidule': 16.0, 'machin': 13.0}
bidule a pour note 16.0

4.3. Structure complexe#

On veut représenter un élève qui est caractérisé son nom, son prénom, son numéro d’étudiant, sa moyenne générale, etc. On voudrait qu’une seule variable conserve et donc donne accès à toutes ces informations.

En algorithmique, on définirait alors un type enregistrement Eleve regroupant ces informations.

Le type Eleve contiens alors une chaîne de caratère (le nom), une deuxième chaîne de caractère (le prénom), un entier (numero étudiant), un réel (moyenne générale)

# representation avec des listes
eleve=['machin','chouette',123456,12.5]
print(eleve)
print("eleve numero:{} moyenne:{}".format(eleve[2],eleve[3]))

['machin', 'chouette', 123456, 12.5]
eleve numero:123456 moyenne:12.5

Problème

accés aux données avec un indice (peu lisible)
complexification

solution - accès avec un mot clé - eleve.nom - eleve.note

4.3.1. Enregistrement (structure)#

Définition algorithmique : Un enregistrement est un type correspondant à un agrégat d’élément de types éventuellement différent auxquels on accède grâce à un nom.

en Python on peut l'implémenter avec la notion de classe

4.3.1.1. structure avec définition explicite#

à préférer !!!

# structure avec définition explicite
class Eleve():
    def __init__(self,name,forname,num,moy):
        self.nom     = name
        self.prenom  = forname
        self.numero  = num
        self.moyenne = moy
        return  

eleve = Eleve('machin','chouette',123456,12.5)
print("Eleve {} numero:{} moyenne:{}".
      format(eleve.nom,eleve.numero,eleve.moyenne))

Eleve machin numero:123456 moyenne:12.5

4.3.1.2. structure dynamique#

uniquement si nécéssaire !!!
préférer une définition explicite

class Eleve():
    pass

eleve = Eleve()
eleve.nom    = 'machin'
eleve.prenom = 'chose'
eleve.numero = 123456
eleve.moyenne= 12.5
print("nom {} numero:{} moyenne:{}".
      format(eleve.nom,eleve.numero,eleve.moyenne))

nom machin numero:123456 moyenne:12.5

4.4. TP: système solaire#

On veut manipuler les planétes du système solaire.

Une planéte est caractérisée par:

son nom
sa masse
sa distance au soleil
sa période de rotation

4.4.1. Définition d’une structure planete#

# unite de masse : terre en kg
masse_terre = 5.9736e24 
class Planete():
    def __init__(self,name,dist,diametre,density,period):
        self.nom   = name
        self.masse = np.pi/6.*(density*1.e12)*\
                    diametre**3/masse_terre
        self.rayon = diametre/2.
        self.distance = dist
        self.periode  = period
        return

4.4.2. Manipulation de la structure planete#

Terre = Planete("terre",150,12756,5.5,365.256)
Mars  = Planete("mars",228,6794,4.0,686.98)

# distance  terre/mars
dmin = -Terre.distance + Mars.distance
dmax =  Terre.distance + Mars.distance
print("distance terre mars: min {} max {} (10^6 km)".format(dmin,dmax))

distance terre mars: min 78 max 378 (10^6 km)

# durée du voyage (en jour) (vitesse du vaiseau en km/h)
V = 15000.0
tmin = (dmin*1e6/V)/24.
tmax = (dmax*1e6/V)/24.
print("duree min={:.2f}  max={:.2f} (jours)".format(tmin,tmax))

duree min=216.67  max=1050.00 (jours)

# liste de planetes
Jupiter=Planete("jupiter",778,143884,1.3,4332.6)
Venus  =Planete("venus",108,12104,5.3,224.701)
Saturne=Planete("saturne",1427,120536,0.7,10759.2)
SystemeSolaire=[Saturne,Terre,Mars,Jupiter,Venus]
#
def affiche(SSolaire):
    for planete in SSolaire:
        print("{:8s} masse {:6.2f} distance {:5.0f}mkm periode {:7.1f}j".format(
          planete.nom, planete.masse, planete.distance, 
          planete.periode))
    return
#
affiche(SystemeSolaire)

saturne  masse 107.45 distance  1427mkm periode 10759.2j
terre    masse   1.00 distance   150mkm periode   365.3j
mars     masse   0.11 distance   228mkm periode   687.0j
jupiter  masse 339.42 distance   778mkm periode  4332.6j
venus    masse   0.82 distance   108mkm periode   224.7j

4.4.3. Tri bulle (bubble sort)#

Le tri bulle est un algorithme classique de tri basé sur la comparaison répétitive d’éléments consécutifs du tableau. Il doit son nom au fait qu’il déplace rapidement les plus grands éléments en fin de tableau, comme des bulles d’air qui remonteraient rapidement à la surface d’un liquide Wikipedia. Cette technique est bien plus efficace que l’algorithme basique consistant à déterminer le maximum du tableau.

L’algorithme de tri bulle parcourt le tableau et compare les éléments consécutifs. Lorsque deux éléments consécutifs ne sont pas dans l’ordre, ils sont échangés. Après un premier parcours complet du tableau du dernier au premier, le plus grand élément est forcément en fin de tableau, à sa position définitive. Le reste du tableau est en revanche encore en désordre. Il faut donc le parcourir à nouveau jusqu’à l’avant-dernier élément. Après ce deuxième parcours, les deux plus grands éléments sont à leur position définitive. Il suffit de répéter cette opération N-1 fois, N étant la dimension du tableau.

Principe de l’algorithme de Tri bulle pour un tableau T de dimension n

Algorithm tri_à_bulles(Tableau T)
       pour n allant de taille(T) - 1 à 1
          "n=nbre d'elements à trier"
          pour j allant de 0 à n - 1
           si T[j+1] < T[j]
               échanger(T[j+1], T[j])

# Implémentation en Python

N = len(SystemeSolaire)
for i in range(N):
    n = N-i  # nbre d'elements à trier    
    for j in range(n-1):
        planete_j  = SystemeSolaire[j]
        planete_j1 = SystemeSolaire[j+1]
        if planete_j1.distance < planete_j.distance:
            # echange des planetes dans le tableau
            SystemeSolaire[j]  = planete_j1
            SystemeSolaire[j+1]= planete_j
# affiche les planetes pour vérifier
affiche(SystemeSolaire)

venus    masse   0.82 distance   108mkm periode   224.7j
terre    masse   1.00 distance   150mkm periode   365.3j
mars     masse   0.11 distance   228mkm periode   687.0j
jupiter  masse 339.42 distance   778mkm periode  4332.6j
saturne  masse 107.45 distance  1427mkm periode 10759.2j

4.4.4. Tracé du système solaire#

4.4.4.1. tracer des planetes#

# adimensionnalisation
dist_terre  = Terre.distance
rayon_terre = Terre.rayon
# couleur des planetes
Couleurs=['b','r','#c0fb2d','#aaa662','#ceb301']
#
plt.figure(figsize=(12,7))
ax = plt.axes(xlim=(-10,10),ylim=(-5,5))
for k,planete in enumerate(SystemeSolaire):
    col= Couleurs[k]
    r1 = planete.distance/dist_terre
    r2 = 0.1*planete.rayon/rayon_terre
    c2 = plt.Circle((r1,0),radius=r2,color=col)
    ax.add_patch(c2)
    plt.text(r1,0.5,planete.nom[0].upper(),fontsize=18)
plt.axis('equal')
plt.xlabel('distance soleil')

Text(0.5, 0, 'distance soleil')

_images/bead28dfef1642ee785dbf11eeaa6b3c4401a7897281cba130d8972009af941c.png

4.4.5. Animation des planetes#

dist_terre  = Terre.distance
rayon_terre = Terre.rayon
Couleurs=['b','r','#c0fb2d','#aaa662','#ceb301']
Cplanetes=[0]*len(SystemeSolaire)
# figure
ax = None
fig = None
#
def init_anim():
    global SystemeSolaire, Cplanetes, fig, ax
    fig=plt.figure(figsize=(6,6))
    ax = plt.axes(xlim=(-10,10),ylim=(-10,10))
    r0 = 0.1
    c0 = plt.Circle((0,0),radius=r0,color='k')
    ax.add_artist(c0)
    for planete in SystemeSolaire:
        r1 = planete.distance/dist_terre
        c1 = plt.Circle((0,0),radius=r1,lw=1,fill=False)
        ax.add_patch(c1)
    plt.axis('equal')
    plt.axis('off')
    return
def init():
    global SystemeSolaire,Cplanetes,ax
    for k,planete in enumerate(SystemeSolaire):
        col= Couleurs[k]
        c2 = plt.Circle((0,0),radius=0,color=col)
        ax.add_patch(c2)
        Cplanetes[k]=c2
    return (Cplanetes[k] for k in range(len(SystemeSolaire)))
#
def animate(i):
    global SystemeSolaire,Cplanetes
    t  = i*20.
    for k,planete in enumerate(SystemeSolaire):
        r1 = planete.distance/dist_terre
        w1 = 2*np.pi/planete.periode
        r2 = 0.12*planete.rayon/rayon_terre
        x  = r1*np.cos(w1*t)
        y  = r1*np.sin(w1*t)
        Cplanetes[k].center = (x,y)
        Cplanetes[k].radius = r2
    return (Cplanetes[k] for k in range(len(SystemeSolaire)))
#

import matplotlib.animation as animation
from matplotlib import rc
init_anim()
anim=animation.FuncAnimation(fig, animate, range(200),interval=100, init_func=init)

_images/869c588f027e9d93468794857751e28144ef174f46fc0e81774ecc55519bf4d5.png

rc('animation', html='jshtml')
anim

4.5. Fin#

print("\t\tSystème utilisé")
import sys
print("Système :\t\t",sys.platform)
import platform
print(platform.platform())
print("Ordinateur:\t\t",platform.machine())
print("Version de Python:\t",sys.version)
import IPython
print("Version de IPython:\t",IPython.__version__)
import numpy
print("Version de numpy:\t",numpy.version.version)
import scipy
print("Version de scipy:\t",scipy.version.version)
import matplotlib
print("Version de matplotlib:\t",matplotlib.__version__)

		Système utilisé
Système :		 linux
Linux-5.15.0-173-generic-x86_64-with-glibc2.35
Ordinateur:		 x86_64
Version de Python:	 3.10.12 (main, Mar  3 2026, 11:56:32) [GCC 11.4.0]
Version de IPython:	 8.37.0
Version de numpy:	 1.26.4
Version de scipy:	 1.15.2
Version de matplotlib:	 3.10.7

5. Gestion de fichiers sous python#

Marc BUFFAT, dpt mécanique, Université Claude Bernard Lyon 1

%matplotlib inline
import numpy as np
import scipy as sp
import matplotlib.pyplot as plt
from matplotlib import rcParams
rcParams['font.family'] = 'serif'
rcParams['font.size'] = 14
from IPython.core.display import HTML
from IPython.display import display
from matplotlib import animation
#css_file = 'style.css'
#HTML(open(css_file, "r").read())

5.1. Fichier informatique#

définition (Wikipédia) un fichier informatique est une collection d’informations numériques (séquences d’octets) réunies sous un même nom, enregistrées sur un support de stockage tel qu’un disque dur, un CD-ROM, une clé USB …

En vue de faciliter leur organisation, les fichiers sont disposés dans des systèmes de fichiers qui permettent de placer les fichiers dans des emplacements appelés répertoires ou dossiers eux-mêmes organisés selon le même principe de manière à former une hiérarchie arborescente

Sur le système de l’ordinateur, un fichier est repéré par son nom avec éventuellement une extension et le répertoire dans lequel il se trouve.

En calcul scientifique, on considère 2 types principaux de fichiers:

des fichiers de programmes qui contiennent un code informatique exécutable par l’ordinateur, sous 2 formes principales:
- fichiers exécutables qui contiennent un code binaire directement exécutable par le processeur (extension .exe sous windows). Ces programmes ne sont pas modifiables et sont propre à l’ordinateur utilisé (programmes commerciaux ou propriétaires)
- fichiers sources qui contiennent le code source du programme dans un langage de programmation (extension .py pour les programmes Python). Ces programmes peuvent être modifiés et être exécuter sur n’importe quel ordinateur.
des fichiers de données qui contiennent des données (data) crées et manipulées par les programmes, sous 2 formes principales:
- des fichiers binaires qui contiennent l’information brute, utilisant le codage de l’ordinateur (le plus efficace pour le stockage)
- des fichiers textes qui contiennent les données écrites sous forme de lignes de texte (format le plus simple à manipuler)

Dans la suite on ne considérera que les fichiers programmes Python (extension .py) et des fichiers de données textes (extension .dat).

Pour éviter des problèmes de compatibilités entre systèmes informatiques, on choisira des noms de fichiers sans caractères accentués, sans espace ni caractéres spéciaux (autre que . ou _ ou - )

5.2. Fichier de programme (ou script) sous Python#

un fichier contenant un programme Python a par convention une extension .py. Un simple éditeur de texte (notepad sous Windows) suffit pour créer et modifier un programme Python

gedit « voir le site gedit » un éditeur simple et efficace
vim et variante « soir le site vim » un éditeur de référence pour les programmeurs
spyder « voir le site spyder » un environnement de développement Python à la matlab pour les applications scientifiques. Sous Spyder, on utilisera plutôt la console Ipython pour exécuter les programmes, que la console par défaut qui utilise des conventions matlab différentes d’un interpréteur Python classique.
jupyter lab extension jupyter pour un environnement de programmation

5.3. Fichiers de données sous Python#

Comme dans la pluspart des langages informatiques, lire ou écrire dans un fichier, on va associer à un fichier (à son nom) une variable informatique de type file, qui posséde des fonctions (ou méthodes) permettant de lire ou écrire des données dans le fichier.

5.3.1. fonction de lecture/écriture#

ouverture d’un fichier mon_fichier.dat en lecture: f=open(« nom_fichier.dat »,”r”)
ou en écriture: f=open(« nom_fichier.dat »,”w”)
lecture / ecriture (caractères) f.read(n) ou f.readline() f.write(chaine)
lecture / écriture de tableaux (numpy) A=loadtxt(file (ou nom de fichier)) savetxt(file (ou nom de fichier), A)
fermeture du fichier f.close()
attention caractères accentués (codage UTF8) encodage des chaines string.encode ouverture fichier (bytes) open(« mon fichier », »wb »)

5.4. Exemple manipulation de courbes de Lissajous#

Une courbe de Lissajous (d’après le physicien français « Jules Antoine Lissajous ») est une courbe paramétrique du plan dont les composantes sont des fonctions périodiques du paramètre (en général le temps en physique, et que l’on peut observer avec un oscilloscope).

On veut manipuler et tracer les courbes de Lissajous suivantes

(5.1)#\[\begin{eqnarray} x(t) &=& \sin{(\frac{2\pi}{p} t)}\\ y(t) &=& \sin{(\frac{2\pi}{q} t + \phi)} \end{eqnarray}\]

Le rapport $n=\frac{p}{q}$ est le paramètre de la courbe et $\phi$ le déphasage. Le temps $T$ de parcours de la courbe est le plus petit commun multiple de $p$ et $q$:

\[T=ppcm(p,q)=\frac{pq}{pgcd(p,q)}\]

5.4.1. Création d’une petite bibliothèque d’application#

lissajous.py : fichier contient des fonctions utiles, et une partie main qui n’est exécutée que si l’on exécute directement le fichier avec un interpréteur python. Cette partie sert en général à tester la bibliothèque.

Cette partie n’est pas exécutée lorsque l’on importe la bibliothèque avec import .

def pgcd(a,b):
    
def Lissajous(p,q,phi,N):

def LissajousFigure(p,q,phi):

# programme principal (main) du module
if __name__ == "__main__":
    # test de la bibliothéque

Attention pour l’utiliser comme bibliothèque, il ne faut écrire aucune instruction python en dehors des fonctions et de la partie main (dans le if)

5.4.2. fichier lissajous.py #

fichier source Python

%%bash --err /dev/null
ls -al lissajous.py
#cat lissajous.py

-rwxrwxr-x 1 buffat buffat 1580 févr. 10  2025 lissajous.py

5.4.3. exécution du programme#

%%bash --err /dev/null
echo "execution avec interpreteur python"
python3 lissajous.py
ls -al lissajous.png

execution avec interpreteur python

-rw-rw-r-- 1 buffat buffat 97884 avril  9 13:57 lissajous.png

5.4.4. utilisation de Jupyter lab#

http://localhost:8888/lab

5.5. module (ou librairie)#

module fichier contenant des fonctions et des définitions en python (bibliothèque)

nom du module nom du fichier sans l’extension .py

utilisation du module mon_module.py

import mon_module
import mon_module as mn
from mom_module import ma_fonction
from mon_module import *

from lissajous import pgcd,LissajousFigure
print(pgcd(21,9))
LissajousFigure(5,3,np.pi/2)

_images/1e9571cca9d5a962fb4ae687b1e7558cbef10862778a06745a0b6a26060e1c6d.png

rem si votre bibliothèque n’est pas dans le répertoire courant, il faudra alors spécifier son chemin d’accès (path) (voir documentation Python)

5.5.1. Exemple: écriture de données#

écriture sur fichiers de $n$ points d’une courbe de Lissajous avec le format suivant

# ligne de commentaire
# n
t0 x0 y0
.......
ti xi yi
.......
tn-1 xn-1 yn-1

5.5.1.1. Programme Python#

import numpy as np
import matplotlib.pyplot as plt
from lissajous import Lissajous
# calcul des points
p=1; q=2; 
N=50
t,x,y=Lissajous(p,q,np.pi/2,N)
# mise sur fichier
F=open("lissajous.dat","w")
F.write("# courbe de lissajous avec p=%d q=%d\n"%(p,q))
F.write("# %d\n"%(N))
np.savetxt(F,np.transpose([t,x,y]))
F.close()

%%bash --err /dev/null
ls -al lissajous.dat
cat lissajous.dat

-rw-rw-r-- 1 buffat buffat 3839 avril  9 13:57 lissajous.dat

# courbe de lissajous avec p=1 q=2
# 50
0.000000000000000000e+00 0.000000000000000000e+00 1.00000000

0000000000e+00
4.081632653061224164e-02 2.536545839095073474e-01 9.917900138232461638e-01
8.16326530

6122448328e-02 4.907175520039378513e-01 9.672948630390294511e-01
1.224489795918367319e-01 6.95682550

6034863826e-01 9.269167573460217469e-01
1.632653061224489666e-01 8.551427630053460849e-01 8.71318704

1233895203e-01
2.040816326530612013e-01 9.586678530366605777e-01 8.014136218679566159e-01
2.44897959

1836734637e-01 9.994862162006878936e-01 7.183493500977276014e-01
2.857142857142856984e-01 9.74927912

1818236193e-01 6.234898018587335944e-01
3.265306122448979331e-01 8.865993063730001067e-01 5.18392568

3105251592e-01
3.673469387755101678e-01 7.402779970753157190e-01 4.047833431223940570e-01
4.08163265

3061224025e-01 5.455349012105490392e-01 2.845275866310327251e-01
4.489795918367346372e-01 3.15108218

0236212671e-01 1.595998950333795963e-01
4.897959183673469274e-01 6.407021998071323055e-02 3.20515775

7165516541e-02
5.306122448979591066e-01 -1.911586287013718743e-01 -9.602302590768126145e-02
5.714285

714285713969e-01 -4.338837391175580094e-01 -2.225209339563142819e-01
6.122448979591835760e-01 -6.482

283953077878635e-01 -3.453650544213070495e-01
6.530612244897958663e-01 -8.201722545969556410e-01 -4.

625382902408350372e-01
6.938775510204081565e-01 -9.384684220497602203e-01 -5.721166601221693293e-01

7.346938775510203357e-01 -9.953791129491981193e-01 -6.723008902613164528e-01
7.755102040816326259e-0

1 -9.871817834144501758e-01 -7.614459583691344235e-01
8.163265306122448051e-01 -9.144126230158128310

e-01 -8.380881048918403797e-01
8.571428571428570953e-01 -7.818314824680299147e-01 -9.009688679024190

350e-01
8.979591836734692745e-01 -5.981105304912169851e-01 -9.490557470106684157e-01
9.3877551020408

15647e-01 -3.752670048793745883e-01 -9.815591569910653291e-01
9.795918367346938549e-01 -1.2787716168

45066350e-01 -9.979453927503363353e-01
1.020408163265306145e+00 1.278771616845061354e-01 -9.97945392

7503363353e-01
1.061224489795918213e+00 3.752670048793733115e-01 -9.815591569910655512e-01
1.1020408

16326530503e+00 5.981105304912150977e-01 -9.490557470106688598e-01
1.142857142857142794e+00 7.818314

824680295816e-01 -9.009688679024192570e-01
1.183673469387755084e+00 9.144126230158122759e-01 -8.3808

81048918410459e-01
1.224489795918367152e+00 9.871817834144499537e-01 -7.614459583691352007e-01
1.265

306122448979442e+00 9.953791129491983414e-01 -6.723008902613174520e-01
1.306122448979591733e+00 9.38

4684220497606644e-01 -5.721166601221699954e-01
1.346938775510204023e+00 8.201722545969559741e-01 -4.

625382902408353702e-01
1.387755102040816313e+00 6.482283953077891958e-01 -3.453650544213082152e-01
1

.428571428571428381e+00 4.338837391175600633e-01 -2.225209339563155031e-01
1.469387755102040671e+00 

1.911586287013736507e-01 -9.602302590768251045e-02
1.510204081632652962e+00 -6.407021998071230073e-0

2 3.205157757165480459e-02
1.551020408163265252e+00 -3.151082180236208230e-01 1.595998950333792354e-

01
1.591836734693877320e+00 -5.455349012105463746e-01 2.845275866310311152e-01
1.632653061224489610e

+00 -7.402779970753142758e-01 4.047833431223928913e-01
1.673469387755101900e+00 -8.86599306372999551

6e-01 5.183925683105244930e-01
1.714285714285714191e+00 -9.749279121818235083e-01 6.2348980185873337

23e-01
1.755102040816326481e+00 -9.994862162006878936e-01 7.183493500977270463e-01
1.795918367346938

549e+00 -9.586678530366613549e-01 8.014136218679558388e-01
1.836734693877550839e+00 -8.5514276300534

70841e-01 8.713187041233888541e-01
1.877551020408163129e+00 -6.956825506034870488e-01 9.269167573460

215248e-01
1.918367346938775420e+00 -4.907175520039380734e-01 9.672948630390293401e-01
1.95918367346

9387710e+00 -2.536545839095085686e-01 9.917900138232460527e-01
2.000000000000000000e+00 -4.898587196

589412829e-16 1.000000000000000000e+00

5.5.2. Lecture des données avec `loadtxt`#

F=open("lissajous.dat","r")
t,x,y =np.loadtxt(F,unpack=True)
F.close()
print("lecture des données taille=",x.shape, y.shape)
# tracer
plt.plot(x,y,lw=2)

lecture des données taille= (50,) (50,)

[<matplotlib.lines.Line2D at 0x7feee906ceb0>]

_images/a28fec1c7d564210519a0e6b431978631c82c4c20800c6d454d13679ac8d321c.png

5.6. Entrée Sortie en python#

5.6.1. Lecture au clavier#

lecture au clavier input d’une valeur sans conversion

val = input("message")

ensuite, on doit faire une conversion

val = int(input("entier ="))
print(val,type(val))
chaine = input("chaine =")
print(chaine,type(chaine))

5.6.2. Écriture avec formatage (mise en forme)#

print(chaine%(val1,val2,..))

chaine est une chaîne de caractères contenant des champs %[n][t]

n = entier optionnel spécifiant la largueur (en colonnes) du champ
t spécifie le type du champ: d=entier (décimal), f=réel (float) , e=réel avec exposant, g=réel (format général), s=chaîne (string)

format avec python 3

string.format(var1,var2,..)

où string est une chaine avec des {} pour spécifier la position des variables (pour le formattage on remplace % par : )

i=1425
x=3242.627
ch="chaine"
print(i,x,ch)

1425 3242.627 chaine

print("i={:} {:08d}  x: {:6.2f} ={:8.1e} titre={}".format(i,i,x,x,ch))

i=1425 00001425  x: 3242.63 = 3.2e+03 titre=chaine

5.6.3. Fonctions utiles de gestion des fichiers#

# liste des fichiers dans un repertoire
import glob
fichiers = glob.glob('*')
print(fichiers)

['Numpy-Matplotlib.ipynb', 'Programmation_questions.ipynb', 'lissajous.png', 'data', 'Variables.ipynb', 'CalculSymbolique.pdf', 'PenduleDouble.ipynb', 'README.md', 'alunissage.py', 'lissajous.dat', 'Algorithme.ipynb', 'lissajous.py', 'images', 'Introduction_questions.ipynb', 'CalculSymbolique.ipynb', 'mafig.pdf', 'anim_pendule.py', 'HistoriqueInfo.md', 'planetes.dat', 'Introduction2.ipynb', 'StructureDonnees.ipynb', 'caslin.png', 'BasePython.ipynb', 'Cours.ipynb', 'BasePython.md', 'BibliothequeScientifique.ipynb', 'BaseProgrammation.ipynb', 'output_35_0.png', 'periode.png', 'README.anim', 'Fichiers.ipynb', 'serie', 'MGC1061M', '__pycache__', 'PythonOpt.ipynb', 'Pendule.ipynb', 'mabib.py', 'monCR.tex', 'treillis.py', 'Introduction.pdf', 'draw_pendule.py', 'cours.tplx', 'lunar_landing.py', 'Pendule.py', 'Introduction.ipynb', 'Cours.pdf', 'StructureDonnees.pdf', 'Pendule.pdf', 'planete.py', 'HistoriqueInfo.ipynb', 'mafig.png', 'casnonlin.png', 'PythonScientifique.ipynb']

5.7. FIN#

1. Modélisation mécanique d’un treillis#

Ce premier cas d’application concerne la modélisation mécanique d’un treillis par éléments finis en statique et dynamique avec une application sur une ferme de toiture industrielle.

Il comprend 4 parties:

la modélisation par éléments finis d’un treillis
un notebook pour le calcul en statique
un notebook pour le calcul en dynamique
un notebook d’application au cas d’une ferme de toiture

1.1. Modélisation d’un treillis par éléments finis#

1.1.1. Introduction#

On se propose d’étudier la modélisation d’un treillis par éléments finis. Un treillis est un assemblage de poutres métalliques (voir figure ci-dessous).

Treillis

On suppose que l’assemblage subit uniquement des efforts de compression (ou extension) sans effet de flexion.

1.1.2. Etude d’une poutre seule#

Pour une poutre de longueur L , de coefficient d’élasticité E et de section S, l’équation d’équilibre statique s’écrit (en notant $u(x)$ le déplacement)

\[\frac{\partial}{\partial x}(ES\frac{\partial u}{\partial x})=0\]

A cette équation, il faut ajouter 2 conditions aux limites:

soit des conditions d’encastrement (condition de Dirichlet):
$u=0$
soit des conditions de force imposée (condition de Neuman):
$ES\frac{\partial u}{\partial x}=F$

Dans la suite on supposera que l’on impose 2 conditions de force en $x=0$ et $x=l$.

1.1.2.1. Formulation éléments finis#

Soit $\delta u$ un déplacement virtuel licite (i.e. vérifiant les C.L. de Dirichlet), la formulation faible s’écrit:

\[ \underbrace{\int_{0}^{l}ES\frac{\partial u}{\partial x}\frac{\partial\delta u}{\partial x}dx}_{\mbox{travail des forces internes}} =\underbrace{F_{2}\delta u(l)-F_{1}\delta u(0)}_{\mbox{travail des forces extérieures}} \]

On cherche une solution approchée $u^{h}$ en utilisant une approximation avec des éléments finis $P^{1}$ sur un maillage d’un seul élément $e_{1}=[0,l]$

\[u^{h}(x)=u_{1}\Phi_{1}(x)+u_{2}\Phi_{2}(x)\]

où $\Phi_{1}(x)$ et $\Phi_{2}(x)$ sont les fonctions base associées aux deux noeuds du maillage $P_{1}(x=0)$ et $P_{2}(x=l)$

La formulation Galerkin s’écrit, en prenant pour fonction test $\delta u^{h}$, les fonctions de base $\Phi_{1}(x)$ et $\Phi_{2}(x)$ :

\[\begin{split}\begin{aligned} u_{1}\int_{0}^{l}ES\frac{\partial\Phi_{1}}{\partial x}\frac{\partial\Phi_{1}}{\partial x}dx+u_{2}\int_{0}^{l}ES\frac{\partial\Phi_{2}}{\partial x}\frac{\partial\Phi_{1}}{\partial x}dx & = & -F_{1}\\ u_{1}\int_{0}^{l}ES\frac{\partial\Phi_{1}}{\partial x}\frac{\partial\Phi_{2}}{\partial x}dx+u_{2}\int_{0}^{l}ES\frac{\partial\Phi_{2}}{\partial x}\frac{\partial\Phi_{2}}{\partial x}dx & = & +F_{2} \end{aligned}\end{split}\]

c’est un système linéaire 2x2

\[\begin{split}\mathbf{K}\left[\begin{array}{c} u_{1}\\ u_{2} \end{array}\right]=\mathbf{B}\end{split}\]

$\mathbf{K}$ est la matrice de rigidité qui vaut pour des éléments $P^{1}$ de longueur $h$

\[\begin{split}\mathbf{K}=\frac{ES}{h}\left[\begin{array}{cc} 1 & -1\\ -1 & 1 \end{array}\right]\,\end{split}\]

1.1.2.2. Cas d’une poutre incliné d’un angle $\alpha$#

poutre

Dans les axes $(\chi,\xi)$ liés à la barre, le déplacement $u$ est donné par les équations précédentes et le déplacement $v$ n’est associé à aucune équation. Le système linéaire sur $u$ et $v$ s’écrit:

\[\begin{split}\left[\begin{array}{cccc} K_{11} & 0 & K_{12} & 0\\ 0 & 0 & 0 & 0\\ K_{21} & 0 & K_{22} & 0\\ 0 & 0 & 0 & 0 \end{array}\right]\left[\begin{array}{c} u_{1}\\ v_{1}\\ u_{2}\\ v_{2} \end{array}\right]=\left[\begin{array}{c} -F_{1}\\ 0\\ F_{2}\\ 0 \end{array}\right]\end{split}\]

Pour passer dans le repère fixe, il faut introduire la matrice de changement de repère:

\[\begin{split}\mathbf{R}=\left[\begin{array}{cc} \cos\alpha & \sin\alpha\\ -\sin\alpha & \cos\alpha \end{array}\right]\end{split}\]

ce qui permet d’écrire le déplacement dans $(\chi,\xi)$ en fonction du déplacement dans $(x,y)$

\[\begin{split}\left[\begin{array}{c} u_{1}\\ v_{1}\\ u_{2}\\ v_{2} \end{array}\right]=\mathbf{T}\left[\begin{array}{c} u_{1x}\\ u_{1y}\\ u_{2x}\\ u_{2y} \end{array}\right]\,\,\mbox{avec\,}\,\,\mathbf{T=}\left[\begin{array}{cccc} R_{11} & R_{12} & 0 & 0\\ R_{21} & R_{22} & 0 & 0\\ 0 & 0 & R_{11} & R_{12}\\ 0 & 0 & R_{21} & R_{22} \end{array}\right]\end{split}\]

et de projeter les forces de $(\chi,\xi)$ dans $(x,y)$

\[\begin{split}\left[\begin{array}{c} F_{x1}\\ F_{y1}\\ F_{x2}\\ F_{y2} \end{array}\right]=\mathbf{T^{t}}\left[\begin{array}{c} F_{1}\\ 0\\ F_{2}\\ 0 \end{array}\right]\,\end{split}\]

ce qui conduit au système linéaire dans $(x,y)$

\[\begin{split}\mathbf{T^{t}}\left[\begin{array}{cccc} K_{11} & 0 & K_{12} & 0\\ 0 & 0 & 0 & 0\\ K_{21} & 0 & K_{22} & 0\\ 0 & 0 & 0 & 0 \end{array}\right]\mathbf{T}\left[\begin{array}{c} u_{1x}\\ u_{1y}\\ u_{2x}\\ u_{2y} \end{array}\right]=\left[\begin{array}{c} F_{1x}\\ F_{1y}\\ F_{2x}\\ F_{2y} \end{array}\right]\end{split}\]

1.1.3. BE étude du treillis#

Un treillis est un assemblage de $Ne$ poutres reliées en $Nn$ noeuds. On veut déterminer les déplacements $(u_{x},u_{y})$ de chaque noeuds, ainsi que les contraintes dans les poutres.

1.1.3.1. étude statique#

Le modèle discret consiste à écrire le système linéaire

\[\mathcal{K}U=\mathcal{B}\]

où $U$ est le vecteur déplacement de $(2\,Nn)$ réels, et $\mathcal{K}$ la matrice de rigidité et $\mathcal{B}$ le second membre (forces appliquées)

1.1.3.1.1. Algorithme:#

Assemblage de la matrice $\mathcal{K}$ (boucle sur les barres)
Calcul du second membre $\mathcal{B}$ (forces appliquées)
Conditions aux limites (Dirichlet, Neuman)
Résolution du système linaire

1.1.3.1.2. Structure de données#

Description du treillis de $Ne$ barres et $Nn$ noeuds:

tableau de coordonnées des noeuds $X$ vecteur $(2Nn)$ réels
graphe du treillis (table de connection) $G$ tableau $(2,Ne)$ entiers
conditions aux limites aux noeuds $CL$ tableau $(Nn)$ entiers
(code = type de CL 0 rien, 1 u=0, 2 v=0, 3 u,v=0, 4 Fe)
Forces extérieures $Fe$ tableau $2Nn$ réels

1.1.3.1.3. Implémentation en python#

utilisation de numpy et de structure (classe)

vecteur réel
X=array([0.2,0.3,0.5])
X=array([0.2,0.3,0.5],dtype=float)
vecteur entier
I=array([1,3,4],dtype=int)
matrice
A=array([[1,2],[3,4]])
concatenation
concatenate((A,B))
concatenate((A,B),axis=1)
hstack((A,B))
vstack((A,B))
vstack((hstack((A,B)),hstack((B,A))))
manipulation
transpose(A)
adressage indirecte
A[ix_([1,0],[2,1])]
resolution système lineaire
linalg.solve(A,B)
valeurs propres
D,V=linalg.eig(A)
inverse
linalg.inv(A)

1.1.3.2. étude dynamique#

Pour une poutre de longueur L , de coefficient d’élasticité E et de section S, l’équation d’équilibre dynamique s’écrit (en notant $u(x,t)$ le déplacement)

\[\rho S\frac{\partial^{2}u}{\partial t^{2}}=\frac{\partial}{\partial x}(SE\frac{\partial u}{\partial x})\]

Matrice élémentaire de masse

\[\begin{split}M=\rho SL\left[\begin{array}{cc} \frac{1}{3} & \frac{1}{6}\\ \frac{1}{6} & \frac{1}{3} \end{array}\right]\end{split}\]

\[\begin{split}\mathcal{M}_{e}=\left[\begin{array}{cccc} M_{11} & 0 & M_{12} & 0\\ 0 & 0 & 0 & 0\\ M_{21} & 0 & M_{22} & 0\\ 0 & 0 & 0 & 0 \end{array}\right]\end{split}\]

Pour le treillis, en vibration libre, on doit résoudre le système d’EDO linéaire:

\[\mathcal{M}\ddot{U}=\mathcal{K}U\]

Ce qui conduit à un problème aux valeurs propres pour la matrice $A$:

\[A\Lambda=\lambda\Lambda\mbox{ avec }A=\mathcal{M}^{-1}\mathcal{K}\]

1.2. Modélisation d’un treillis en statique#

Marc BUFFAT, Département mécanique UCB Lyon 1 _images/treillisfig1.png

%matplotlib inline
# bibliothéque
import numpy as np
import scipy as sp
import matplotlib.pyplot as plt
from matplotlib import rcParams
rcParams['font.family'] = 'serif'
rcParams['font.size'] = 14

1.2.1. Introduction#

On se propose d’étudier la modélisation d’un treillis par éléments finis.

Le treillis étudié est un assemblage de poutres métalliques (voir figure) soumis à des efforts de traction compression _images/treillisfig1.png

### Modélisation du treillis

$N_e=2$ barres et $N_n=3$ noeuds
Assemblage des poutres: en statique système linéaire de dimension $2N_n = 6$ $$ K [u] = B $$
- $K$ matrice de rigidité $2N_n\times 2N_n$
- $[u]$ vecteur déplacements des noeuds $2N_n$

1.2.2. Structure de données#

Description du treillis de $N_e$ barres et $N_n$ noeuds:

tableau de coordonnées des noeuds X : vecteur $2N_n$
graphe du treillis (table de connection) G : tableau $(2,N_e)$
conditions aux limites de Dirichlet aux noeuds (code = type de CL)
numérotation à partir de 0 (convention python)

CL vecteur $N_n$ - 0 libre (pas de CL) - 1 déplacement nul suivant x - 2 déplacement nul suivant y - 3=1+2 déplacement nul suivant x et y

FCL tableau $N_n,2$ - composante de la force extérieure imposée au noeud

propriétés des barres: section, module d’Young, masse volumique

S, E, rho

# données brutes du problème
# ===================
# coordonnées des 3 nds
X = np.array([[0.0,0.0],[0.5,0.5],[1.0,0.0]])
# table de connection (numérotation naturelle à partir de 1)
Tbc = np.array([[1,2],[2,3]])
# parametres des barres en unité SI
E = 200.e9
S = 0.000025
rho = 8000
# Force appliquée en unité SI
F = 10000

1.2.3. Traitement des données#

A partir des données brutes, construire les données utiles pour la résolution du problème

numérotation python Tbc -> G
tableau des CL et des forces FCL

G = None
Nn = None
Ne = None
CL = None
FCL = None
### BEGIN SOLUTION
### END SOLUTION

1.2.4. Structure de données#

Regrouper tous les tableaux précédents dans une structure unique de type classe

 class Treillis():

approche générique
initialisation avec les tableaux génériques :
- table de connexion et coordonnées

# structure de donnees treillis: données table de connexion G et coordonnées de nds
class Treillis():
    def __init__(self,G,X):
        '''
           initialisation d un treillis a partir de la table de connexion G
           et des coordoonees X des noeuds. Attention: numérotation à partir de 0
        '''
        return

# utilisation: création du treillis, tableaux des CL et parametres 
tr = Treillis(G,X)
### BEGIN SOLUTION
### END SOLUTION

1.2.5. Affichage de la structure#

fonction print_treillis(tr,titre)
fonction trace_treillis(tr,titre)

def print_treillis(treillis,titre):
    '''affiche les informations du treillis avec un titre optionel'''

# vérification
print_treillis(tr,"mon treillis")

def trace_treillis(tr,titre):
    '''trace d un treillis tr avec un titre'''
    return

# vérification
trace_treillis(tr,"treillis initial")

1.2.6. Amélioration du tracé#

on souhaite améliorer la fonction précédente, en ajoutant la possibilité de tracer le numéro des noeuds.

utilisation d’un paramètre optionnel numero

def trace_treillis(tr,titre,numero=False):
    '''trace d un treillis tr avec un titre'''
    return

# vérification
trace_treillis(tr,"treillis initial",numero=True)

1.2.7. Algorithme#

structure de donnees: Treillis

initialisation des données du treillis
fonctions pour calculer des données nécessaires:
- longueur des barres,
- matrice de rotation (cosinus directeur)
assemblage de la matrice avec une boucles sur les barres
application des conditions aux limites
résolution du système linéaire
analyse du résultat

démarche

chaque partie correspond à l’écriture d’une fonction
on valide à chaque étape
validation de l’algorithme global sur un cas simple
application au cas étudié

1.2.8. Fonctions utiles#

bibliothéque treillis.py

longueur_treillis(tr) calcul des longueurs des barres
matrice_rotation(tr,l) calcul de matrice de rotation de la barre l

Soit $R$ la matrice de rotation:

\[\begin{split} \mathbf{R}=\left[\begin{array}{cc} \cos\alpha & \sin\alpha\\ -\sin\alpha & \cos\alpha \end{array}\right] \end{split}\]

la matrice de changement de repère $T$ s’écrit: $$ \mathbf{T=}\left[\begin{array}{cccc} R_{11} & R_{12} & 0 & 0\\ R_{21} & R_{22} & 0 & 0\\ 0 & 0 & R_{11} & R_{12}\\ 0 & 0 & R_{21} & R_{22} \end{array}\right] $$

def longueur_treillis(tr):
    """ calcul longueur des barres """
    # longueur des barres

# vérification
L = longueur_treillis(tr)
print("Longueur des barres L=",L)

# matrice de transfert T tq Ae = Tt.K.T
def matrice_rotation(tr,l):
    """ calcul matrice de transfert de l'element l (de 0 a ne-1) """

# vérification
T0 = matrice_rotation(tr,0)
print("T0=",T0)
T0 @ T0.transpose() 

1.2.9. Assemblage de la matrice#

Pour une barre le système $(4x4)$ s’écrit

\[\begin{split} \mathbf{T^{t}}\left[\begin{array}{cccc} K_{11} & 0 & K_{12} & 0\\ 0 & 0 & 0 & 0\\ K_{21} & 0 & K_{22} & 0\\ 0 & 0 & 0 & 0 \end{array}\right]\mathbf{T}\left[\begin{array}{c} u_{1x}\\ u_{1y}\\ u_{2x}\\ u_{2y}\end{array}\right] =\left[\begin{array}{c} F_{1x}\\ F_{1y}\\ F_{2x}\\ F_{2y} \end{array}\right] \end{split}\]

qu’il faut assembler dans un système global $2 N_n x N_n$

\[ A U = F \]

assemblage(tr) fonction d’assemblage pour le calcul de la matrice A (sans CL)

1.2.9.1. Algorithme assemblage#

 A matrice globale dimension (2nn,2nn)
 boucle sur les barres : k de 0 a ne-1    
     n1,n2 numéros des 2 sommets de la barre k
     2n1, 2n1+1, 2n2, 2n2+1  numéros des 4 ddl de la barre
     calcul de la matrice élémentaire (4*4) Ae = Tt K T
     ajout de la contribution Ae dans la matrice globale A:
         lignes 2n1,2n1+1,2n2,2n2+1  colonnes 2n1,2n1+1,2n2,2n2+1

# assemblage de la matrice globale
def assemblage(tr):
    """ assemblage de la matrice de rigidite"""
    # matrice de rigidite elementaire
    Ke=np.array([[1,0,-1,0],[0,0,0,0],[-1,0,1,0],[0,0,0,0]]);
    A=np.zeros((2*tr.nn,2*tr.nn));

# validation
# Assemblage de la matrice
A=assemblage(tr)
# verification A= At
print("Symétrie:",np.allclose(A,A.transpose()))
# somme des lignes
print("Somme des lignes:",np.sum(A,axis=1))
# déterminant
print("déterminant ",np.linalg.det(A))
#

1.2.10. Applications de CL#

sur la matrice A et le second membre B: - climites(tr,AA,BB) - si déplacement nul - met un 1 sur la diagonale, 0 sur la ligne et la colonne et dans B - forcage extérieure dans B

def climites(tr,AA,BB):
    '''imposition des CL dans la matrice A et des forces dans le 2nd membre B'''
    # attention on ne modifie pas les arguments mais on renvoie une nouvelle matrice A et B
    A = AA.copy()
    B = BB.copy()

# verification
# second membre
B=np.zeros((2*tr.nn))
# CLimites
A,B=climites(tr,A,B)
print("déterminant ",np.linalg.det(A))
print("B=\n",B)

1.2.11. Résolution#

utilisation de la bibliothèque linalg pour résoudre le système linéaire

fonction solve
transformation du vecteur solution $2*Nn$ en un tableau de déplacement $(Nn,2)$

# calcul de la solution
U=np.linalg.solve(A,B)
U=U.reshape((tr.nn,2))
print("Deplacement U=\n",U)
print("RDM: ",F/(E*S/L[0]))

1.2.12. Forces volumiques#

pour une poutre L suportant une largueur de toiture H avec masse par unité de surface $\rho$ $$ F = \rho g LH$$
répartition sur les 2 noeuds suivant force suivant y : $F/2$
charge supplémentaire du vent ou à la neige (attention direction)

\[p \approx \frac{1}{2} \rho_{air} Ue^2\]

1.2.13. Fin#

1.3. Vibration libre d’un treillis#

Marc BUFFAT, Département mécanique UCB Lyon 1

%matplotlib inline
import numpy as np
import scipy as sp
import matplotlib.pyplot as plt
from matplotlib import rcParams
rcParams['font.family'] = 'serif'
rcParams['font.size'] = 14
from IPython.core.display import HTML
from IPython.display import display,Image
from matplotlib import animation
#from JSAnimation import IPython_display
#css_file = 'style.css'
#HTML(open(css_file, "r").read())

1.3.1. Vibrations d’un système a 1ddl#

\[ M \ddot{u} + K u = f(t)\]

étude en fréquence $f(t)=F \sin{\omega t} $
solution $u(t) = U \sin{\omega t}$

\[ U = \frac{F}{K -\omega^2 M} \]

raisonnance pour $\omega = \sqrt{K/M} $

dans ce cas la solution s’écrit

\[ u(t)=-\frac{F}{2\omega} t \cos{\omega t} \]

1.3.2. Modes de vibration libre#

Modèle analytique $u(x,t)$

\[\rho S \frac{\partial^2 u}{\partial t^2} - E S \frac{\partial^2 u}{\partial x^2} = 0\]

Modèle discret: vecteur $U(t)$

\[ M \frac{ d^2 U }{dt^2 } + K U = 0\]

modes propres de pulsation

\[ U = \cos{\omega t} \;\Lambda \]

problème aux valeurs propres généralisées $\omega=\sqrt{\lambda}$

\[ \lambda M \Lambda = K \Lambda \]

matrice dimension $2Nn \times 2Nn $ , donc $2Nn$ valeurs propres

attention élimination des modes rigides (associés au CL en déplacement)

1.3.3. Matrice élémentaire de rigidité#

\[\begin{split} Ke = \frac{E S }{L} \left[\begin{array}{cc} 1 & -1 \\ -1 & 1 \\ \end{array}\right] \end{split}\]

1.3.4. Matrice élementaire de Masse#

formulation EF P1 $$ Me = \frac{\rho S L }{6} \left[\begin{array}{cc} 2 & 1 \\ 1 & 2 \\ \end{array}\right] $$
avec condensation de masse $$ Me = \frac{\rho S L }{2} \left[\begin{array}{cc} 1 & 0 \\ 0 & 1 \\ \end{array}\right] $$

1.3.5. Equations Elements finis#

La formulation EF conduit à un système différentielle $$ M \ddot{U} + K U = 0$$

1.3.5.1. Cas d’une barre encastrée en $x=0$#

avec un seul élèment $P^1$ à 2 ddl $U=\{u_1,u_2\}$

\[u_h(x,t) = u_1(t) \Phi_1(x) + u_2(t) \Phi_2(x) \]

La condition aux limites impose $u_1=0$, donc un seul ddl $u_2$

\[u_h(x,t) = u_2(t) \Phi_2(x) \]

Equation différentielle ordinaire

\[\frac{\rho S L }{3} \ddot{u}_2 + \frac{E S}{L} u_2 = 0 \]

Solution sinusoidale $u_2=\cos{\omega_1 t}$

\[ \omega_1 = \sqrt{3} \; \sqrt{\frac{E S}{\rho S L^2}} \mbox{ ou } \omega_1 = \sqrt{2} \; \sqrt{\frac{E S}{\rho S L^2}}\]

mode discret

\[ u_1(x,t) = \frac{x}{L} cos(\sqrt{3} \; \sqrt{\frac{E S}{\rho S L^2}} t )\]

1.3.6. Modele analytique ( 1 barre)#

\[ \begin{align}\begin{aligned}\rho S \frac{\partial^2 u}{\partial t^2} - E S \frac{\partial^2 u}{\partial x^2} = 0\\soit\\$$\frac{\partial^2 u}{\partial t^2} - c^2 \frac{\partial^2 u}{\partial x^2} = 0\end{aligned}\end{align} \]

avec $c = \sqrt{\frac{E S}{\rho S}}$ et les conditions aux limites $u_{(x=0)}=0$ et ${\frac{\partial u}{\partial x}}_{(x=L)}= 0 $

solution variables séparées $ u(x,t)= X(x)*Y(t) $

\[ c^{2}\frac{\ddot{X}}{X} = \frac{\ddot{Y}}{Y} = \omega^2\]

avec les conditions aux limites

\[\omega_k = (2k+1)\frac{\pi}{2L} c \]

kième mode propre

\[u^k (x,t) = sin((2k+1) \frac{\pi x}{2 L }) cos((2k+1)\frac{\pi}{2}\sqrt{\frac{E S}{\rho S L^2}} t)\]

premier mode de vibration

\[ u^1 (x,t) = sin(\frac{\pi x}{2 L }) cos(\frac{\pi}{2}\sqrt{\frac{E S}{\rho S L^2}} t)\]

1.3.7. Simulation numérique: mode discret#

1.3.7.1. treillis 1 barre encastré 1ddl#

solution du problème aux valeurs propres généralisées

\[ A \Lambda = \lambda B \Lambda \]

avec

\[\begin{split} A = M = \left[\begin{array}{cccc} 1 & 0 & 0 & 0\\ 0 & 1 & 0 & 0\\ 0 & 0 & \frac{\rho S L }{3} & 0\\ 0 & 0 & 0 & 1\\ \end{array}\right] \mbox{ et } B = K = \left[\begin{array}{cccc} 1 & 0 & 0 & 0\\ 0 & 1 & 0 & 0\\ 0 & 0 & \frac{E S}{L} & 0\\ 0 & 0 & 0 & 1\\ \end{array}\right] \;\; \end{split}\]

# A.N. resolution 
from scipy.linalg import eig
rho=8000; S=0.05**2; L=1.0; E=2.e08; 
M=np.array([[1,0,0,0],[0,1,0,0],[0,0,rho*S*L/3,0],[0,0,0,1]])
print("M=",M)
K=np.array([[1,0,0,0],[0,1,0,0],[0,0,E*S/L,0],[0,0,0,1]])
print("K=",K)
vp,UP=eig(K,M)
print("vp=",vp) 
print(UP.shape)

# 1 mode de vibration (3ieme)
n1=2
omega =np.sqrt(vp[n1].real)
omega1=np.sqrt(3*E*S/(rho*S*L**2))
omegae=np.pi/2*np.sqrt(E*S/(rho*S*L**2))
print("omega calcul=%g mode 1=%g exacte=%g"%(omega,omega1,omegae))
print("VP=",UP[:,n1])

U0 = 0.1
X=np.linspace(0,L,20)
Y=np.sin(np.pi*X/(2*L))
YP1 = UP[2,n1]*X
plt.plot(X,U0*YP1,'-',label='EF P1')
plt.plot(X,U0*Y,'-',label='exacte')
plt.legend(loc=0)
plt.xlabel("X")
plt.ylabel("U")
plt.title("1er mode de vibration");

1.3.8. vibration de la barre#

I = [0,5,10,15,19]
T = np.linspace(0,4*np.pi/omegae,100)
plt.figure(figsize=(12,6))
plt.subplot(1,2,1)
for i in I:
    U = YP1[i]*np.cos(omega*T)
    plt.plot(T,U,label="x={:.1f}".format(X[i]))
plt.legend()
plt.xlabel("t")
plt.ylabel("U(t)")
plt.title("vibration mode 1 (P1)");
plt.subplot(1,2,2)
for i in I:
    U = Y[i]*np.cos(omegae*T)
    plt.plot(T,U,label="x={:.1f}".format(X[i]))
plt.legend()
plt.xlabel("t")
plt.ylabel("U(t)")
plt.title("vibration mode 1 (exacte)");

1.3.9. Avec des éléments finis $P^2$#

\[\begin{split} \mathbf{K} = \frac{2ES}{L}\left[\begin{array}{ccc} +\frac{7}{6} & +\frac{1}{6} & -\frac{4}{3}\\ +\frac{1}{6} & +\frac{7}{6} & -\frac{4}{3}\\ -\frac{4}{3} & -\frac{4}{3} & +\frac{8}{3} \end{array}\right] \mbox{ et } \mathbf{M} = \frac{\rho S L}{2}\left[\begin{array}{ccc} +\frac{4}{15} & -\frac{1}{15} & +\frac{2}{15}\\ -\frac{1}{15} & +\frac{4}{15} & +\frac{2}{15}\\ +\frac{2}{15} & +\frac{2}{15} & +\frac{16}{15} \end{array}\right] \end{split}\]

M=rho*S*L/30.*np.array([[4.,-1.,2.],[-1.,4.,2.],[2.,2.,16.]])
M[0,:]=0; M[:,0]=0; M[0,0]=1.0
print("M=",M)
K=E*S/L/3.*np.array([[7.,1.,-8.],[1.,7.,-8.],[-8.,-8.,16.]])
K[0,:]=0; K[:,0]=0; K[0,0]=1.0
print("K=",K)
vp,UP=eig(K,M)
print("vp=",vp.real)
print("VP=",UP)

omega =np.sqrt(vp[0].real)
omega1=np.sqrt(3*E*S/(rho*S*L**2))
omegae=np.pi/2*np.sqrt(E*S/(rho*S*L**2))
print(omega,omega1,omegae)
print(UP[:,0])

X  =  np.linspace(0,L,20)
Ye = UP[1,0]*np.sin(np.pi*X/(2*L));
Y1 = UP[1,0]*2*X*(X-L/2) + UP[2,0]*4*X*(L-X)
plt.plot(X,-U0*Y1,'-',label='EF P2',lw=2)
plt.plot(X,-U0*Ye,'-',label='exacte')
plt.legend(loc=0)
plt.xlabel("X")
plt.ylabel("U")
plt.title("1er mode de vibration");

X  =  np.linspace(0,L,20)
Ye = -UP[1,1]*np.sin(3*np.pi*X/(2*L));
Y1 =  UP[1,1]*2*X*(X-L/2) + UP[2,1]*4*X*(L-X)
plt.plot(X,U0*Y1,'-',label='EF P2',lw=2)
plt.plot(X,U0*Ye,'-',label='exacte')
plt.legend(loc=0)
plt.xlabel("X")
plt.ylabel("U")
plt.title("2nd mode de vibration");

1.3.10. Application à un treillis#

import sys
from scipy import linalg
from Treillis.treillis import *
# resolution du treillis en statique
# ==================================
fichier = "Treillis/mon_treillis1.dat"
# creation du treillis
tr=lecture_treillis(fichier)
tr.E = 200*1.e9
tr.S = 0.000025
tr.rho = 8000
# CL
tr.CL[0]=3
tr.CL[3]=3
tr.FCL[5,:]=[0.,-10000]
# trace du treillis
U = np.zeros((tr.nn,2))
# Assemblage de la matrice
AA=assemblage(tr)
# second membre
BB=np.zeros((2*tr.nn))
# CLimites
A,B=climites(tr,AA,BB)
# solution
U=np.linalg.solve(A,B)
U=U.reshape((tr.nn,2))
print("Deplacement U=\n",U)
# trace de la deformee
plt.figure(figsize=(12,6))
plt.subplot(1,2,1)
trace_treillis(tr,"Treillis initial")
plt.subplot(1,2,2)
trace_treillis(tr,"Deformee",U*100)

# en dynamique
AA=assemblage(tr)
# second membre
BB=np.zeros((2*tr.nn))
# CLimites
K,B=climites(tr,AA,BB)
# masse
AA=assemblage_masse(tr)
# second membre
BB=np.zeros((2*tr.nn))
# CLimites
M,B=climites(tr,AA,BB)
#
vp,VP=linalg.eig(K,b=M)
vp=np.real(vp)
ip=np.argsort(vp)
print("IP=",ip)
print("VP=",vp[np.ix_(ip)])
# elimination modes rigides (nbre CL =1,2,3)
nr = np.count_nonzero(tr.CL == 1) + np.count_nonzero(tr.CL == 2) 
nr = nr + 2* np.count_nonzero(tr.CL == 3)
print("Nbre de modes rigides ",nr)
fp=np.sqrt(vp[np.ix_(ip)])/(2*np.pi)
fp=fp[nr:]
print("Frequences propres ",fp)
# tracer 1er mode
mode=nr
print("mode ",mode," f=",fp[mode-nr])
U1=VP[:,ip[mode]]/linalg.norm(VP[:,ip[mode]])
U1=U1.reshape((tr.nn,2))
print(U1)
# 2nd mode
mode=nr+1
print("mode ",mode," f=",fp[mode-nr])
U2=VP[:,ip[mode]]/linalg.norm(VP[:,ip[mode]])
U2=U2.reshape((tr.nn,2))
print(U2)

# tracer du mode
plt.figure(figsize=(12,6))
plt.subplot(1,2,1)
#trace_treillis(tr,"",U*0.0,False)
trace_treillis(tr,"1er mode vib",U1*0.3,False)
# tracer du mode
plt.subplot(1,2,2)
#trace_treillis(tr,"",U*0.0,False)
trace_treillis(tr,"2nd mode vib",U2*0.3,False)

# animation
from matplotlib import animation
from IPython.core.display import HTML
def anim(mode):
    w = 2*np.pi*fp[mode]
    tt=np.linspace(0,2*np.pi/w,20)
    U1=VP[:,ip[mode]]/linalg.norm(VP[:,ip[mode]])
    U1=U1.reshape((tr.nn,2))
    # tracer
    fig = plt.figure(figsize=(8,4))      
    ax  = plt.axes()
    plt.axis('equal')
    plt.axis('off')
    Xd = tr.X + 0.5*U1
    pts, = ax.plot(Xd[:,0],Xd[:,1],'o',markersize=10,color='b')
    bars = [None]*tr.ne
    for i in range(tr.ne):
        n1=tr.G[i,0]
        n2=tr.G[i,1]
        bars[i],= ax.plot([Xd[n1,0],Xd[n2,0]],[Xd[n1,1],Xd[n2,1]],'-g',lw=3)

    def plot_mode(t):
        U11 = U1*np.cos(w*t)
        Xd = tr.X + 0.5*U11
        pts.set_xdata(Xd[:,0])
        pts.set_ydata(Xd[:,1])
        for i in range(tr.ne):
            n1=tr.G[i,0]
            n2=tr.G[i,1]
            bars[i].set_xdata([Xd[n1,0],Xd[n2,0]])
            bars[i].set_ydata([Xd[n1,1],Xd[n2,1]])
        return
    an=animation.FuncAnimation(fig, plot_mode, frames=tt)
    plt.title("vibration mode {} f={:.0f}hz".format(mode,fp[mode]))
    plt.show()
    return an

# selection mode
mode=nr+2
an = anim(mode)

HTML(an.to_html5_video())

1.3.11. FIN#

2. Modélisation de l’effet Magnus#

Ce cas d’application concerne la modélisation de l’effet magnus, et son application à la propulsion.

Il comprend 4 parties:

la modélisation de l’effet magnus
la téorie dans le cadre des écoulements potentiels
la simulation de l’écoulement autour d’un cylindre en rotation
l’application au cas d’une turbovoile pour la propulsion des navires

2.1. Modélisation de l’effet Magnus#

2.1.1. Introduction#

effet magnus

Effet Magnus (wikipedia)

: L’effet Magnus, découvert par Heinrich Gustav Magnus (1802-1870), physicien allemand, est un principe physique qui explique la force tangentielle subie par un objet en rotation se déplaçant dans un fluide. C’est cette force qui explique la modification de trajectoire d’une balle de tennis ou un ballon de football lifté . Cet effet est également utilisé comme moyen de sustentation ou de propulsion.

2.1.1.1. Principe#

Le cylindre en rotation dans l’écoulement de fluide incompressible ci-dessous va accélérer le fluide en dessus du cylindre (car la vitesse de rotation du cylindre s’ajoute à celle du fluide), alors qu’il décélère le fluide en dessous (car la vitesse de rotation du cylindre est opposée à celle du fluide).

En négligeant la dissipation visqueuse, et en utilisant la relation de Bernoulli, cette variation de vitesse induit une dépression en dessus du cylindre et une surpression en dessous, d’où l’apparition une force perpendiculaire à l’écoulement et dirigée vers le haut.

principe effet magnus

2.1.1.2. Application: Propulsion de navire#

Flettner rotor ship (1920) par Anton Flettner (aérodynamicien allemand)

E-ship (2010) transport de pales de turbine (Enercon GmbH.)

2.1.2. Analyse dimensionnelle#

5 paramètres:

: $\rho$ $U_{0}$ $R$ $\omega$ et $F_{p}$ (force de portance)

3 unités:

: $m$ , $s$ , $kg$

2 nombres sans dimension:

\[C_{p}=\frac{F_{p}}{\rho U_{0}^{2} R}\mbox{ et }\alpha=\frac{\omega R}{U_{0}}\]

Forme sans dimension de la loi de portance (effet magnus)

\[C_{p}=f(\alpha)\]

2.1.3. Définition du problème#

On considère l’écoulement d’un fluide incompressible (air à des faibles vitesses d’écoulement), de masse volumique $\rho$ constante, autour d’un obstacle cylindrique de rayon R et d’axe Oz.

L’écoulement est à deux dimensions (vitesses parallèles au plan xOy et indépendantes de z) et stationnaire. Un point M du plan xOy est repéré par ses coordonnées polaires $(r,\theta)$. L’obstacle, dans son voisinage, déforme les lignes de courant ; loin de l’obstacle, le fluide est animé d’une vitesse uniforme $U_{0}$ suivant Ox.

Repère sur le cylindre {width= »60% »}

L’écoulement est supposé irrotationnel, et le fluide est supposé parfait (viscosité nulle).

Dans ce cas l’écoulement est un écoulement potentiel, c.a.d que la vitesse découle d’un champ potentiel $\Phi(x,y)$:

\[\begin{split}\overrightarrow{U}=\overrightarrow{grad}\,\Phi=\left[\begin{array}{c} \frac{\partial\Phi}{\partial x}\\ \frac{\partial\Phi}{\partial y} \end{array}\right]\end{split}\]

qui vérifie une équation de Laplace:

\[\Delta\Phi=\frac{\partial^{2}\Phi}{\partial x^{2}}+\frac{\partial^{2}\Phi}{\partial y^{2}}=0\]

Les conditions aux limites sont telles que loin de l’obstacle la vitesse est égale à $U_{0}$, et que le fluide glisse sur l’obstacle, c.a.d. que la vitesse normale est nulle $\overrightarrow{U}.\overrightarrow{n}=0$.

2.1.3.1. Solution analytique#

La solution analytique de ce problème peut être calculée avec la théorie des fonctions conformes (champ uniforme plus dipôle):

\[\Phi=U_{0}(\frac{R^{2}}{r}+r)\cos\theta\]

Cette solution est tracée sur la figure ci-dessous. Le champ de vitesse est perpendiculaire aux lignes de potentiel ( $\overrightarrow{U}$ est parallèle à $\overrightarrow{grad}\Psi$ ). Sur le cercle $(r=R)$, la vitesse normale est nulle $(\frac{\partial\Phi}{\partial r}=0)$ et la vitesse tangentielle vaut:

\[U_{\theta}=\frac{1}{r}\frac{\partial\Phi}{\partial\theta}=-2U_{0}\sin\theta\]

potentiel psi

Cette fonction potentielle $\Phi$ est associée à une fonction orthogonale $\Psi$ (tracée ci-dessus), telle que

\[\begin{split}\overrightarrow{U}=\left[\begin{array}{c} \,\frac{\partial\Psi}{\partial y}\\ -\frac{\partial\Psi}{\partial x} \end{array}\right]\end{split}\]

Cette fonction $\Psi$ est la fonction de courant de l’écoulement, i.e. les lignes $\Psi=cste$ sont des lignes tangentes au vecteur vitesse ( $\overrightarrow{U}.\overrightarrow{grad}\,\Psi=0$). L’écoulement étant stationnaire, ces lignes correspondent aux trajectoires des particules fluides. On montre facilement que $\Psi$ vérifie aussi une équation de Laplace:

\[\Delta\Psi=0\]

Les conditions aux limites associées sont:

(2.1)#\[\Psi(r=\infty)=yU_{0}\,\,,\,\,\,\Psi(r=R)=0\]

2.1.3.2. Simulation par éléments finis#

Pour résoudre ce problème par élément finis, on choisit un domaine de calcul $\Omega$ fini: i.e. la frontière $r=\infty$ est ramenée à une distance finie $r=10R$ du cylindre (à cette distance la perturbation due au cylindre est négligeable): $\Omega=[R,10R]x[0,2\pi]$. On note $\Gamma_{0}(r=10R)$ la frontière extérieure et $\Gamma_{1}(r=R)$ la frontière du cylindre.

On découpe le domaine $\Omega$ en $N_{e}$ éléments triangulaires (triangulation):

maillage éléments finis

A partir de ce maillage, on construit ensuite une solution approchée $\Psi^{h}(x,y)$ comme une approximation polynomiale par morceaux. Ainsi avec une approximation quadratique (élément finis $P^{2}$), sur chaque élément $e_{k}$, la solution est un polynôme de degré 2 en x et y:

\[\Psi^{h}(x,y)=a_{0}+a_{1}x+a_{2}y+a_{3}xy+a_{4}x^{2}+a_{5}x^{6}\]

et possède donc 6 degrés de liberté, qui sont les 3 valeurs aux sommets de l’élément et les 3 valeurs sur les milieux de cotés

element P2

Compte tenu des conditions aux limites et de la continuité inter-éléments, la solution $\Psi^{h}$ possède $N$ degrés de liberté (ce nombre dépend du nombre d’éléments et du nombre de noeuds sur les frontières) et s’écrit:

\[\Psi^{h}(x,y)=\sum_{j=1}^{N}\Psi_{j}\phi_{j}(x,y)+\sum_{k\in\Gamma_{0}\cup\Gamma_{1}}\Psi_{0}\Phi_{k}(x,y)\]

Le terme en $\Psi_{0}$ permet de vérifier les conditions aux limites:

\[\Psi_{0}|_{\Gamma_{0}}=U_{0}y\,\,,\,\,\Psi_{0}|_{\Gamma_{1}}=0\]

Pour déterminer les valeurs inconnues $\Psi_{j}$, on écrit la formulation faible de la méthode des résidues pondérés, en multipliant l’équation [(2.1)] par une fonction test $w_{i}(x,y)$, puis en intégrant sur le domaine $\Omega$:

\[\int_{\Omega}\Delta\Psi^{h}\,w_{i}\,dxdy=0\]

On intègre par partie, en utilisant la formule de Green, pour obtenir la relation suivante:

\[\int_{\Omega}\overrightarrow{grad}\Psi^{h}.\overrightarrow{grad}w_{i}\,d\Omega-\int_{\Gamma}\overrightarrow{grad}\Psi^{h}.\overrightarrow{n}\,w_{i}\,d\Gamma=0\]

En utilisant la méthode de Galerkin, les fonctions tests $w_{i}$ sont les fonctions de base $\phi_{i}$ associées au degré de liberté de $\Psi^{h}$:

\[w_{i}(x,y)=\frac{\partial\Psi^{h}}{\partial\Psi_{i}}=\phi_{i}(x,y)\]

Ces fonctions tests s’annulent sur la frontière $\Gamma$ de $\Omega$, et donc l’intégrale de bord sur $\Gamma$ est nulle. En remplaçant $\Psi^{h}$ par son expression éléments finis, la formulation faible s’écrit:

\[\sum_{j=1}^{N}\Psi_{j}\int_{\Omega}\overrightarrow{grad}\phi_{j}.\overrightarrow{grad}\Phi_{i}\,d\Omega=-\sum_{k\in\Gamma_{1}}\Psi_{0}\int_{\Omega}\overrightarrow{grad}\phi_{k}.\overrightarrow{grad}\Phi_{i}\,d\Omega\,\,\,\forall i=1,N\]

qui est un système linéaire de $N$ équations pour $N$ inconnues $\{\Psi_{j}\}_{j=1,N}$ qu’il suffit de résoudre pour obtenir la solution approchée $\Psi^{h}$.

2.1.3.3. Résolution avec COMSOL#

COMSOL est un programme qui permet de calculer l’approximation par éléments finis d’un problème physique. Les différentes étapes de la résolution sont:

choix du modèle physique (EDP) et choix du degré d’approximation (éléments $P^{2}$ par défaut).
création de la géométrie
définition des coefficients du modèle dans $\Omega$
définition des conditions aux limites sur la frontière $\Gamma=\partial\Omega$
maillage du domaine de calcul
résolution
analyse et tracé de la solution

2.1.3.4. Validation du résultat numérique#

La solution analytique de ce problème s’écrit:

\[\Psi_{e}(x,y)=U_{0}r\,\sin\theta\,(1-\frac{R^{2}}{r^{2}})\]

On peut alors calculer l’erreur entre cette solution exacte et la solution par élément finis, et on trouve

\[-9.3\,10^{-3}\leq\Psi_{e}-\Psi^{h}\leq9.3\,10^{-3}\,\,\,\mbox{avec}\,\,\,\left|\psi_{e}\right|\leq0.93\]

En analysant cette erreur, on peut s’apercevoir qu’une partie de l’erreur est due à la condition aux limites sur $\Gamma_{1}$, car avec les conditions aux limites imposées on a $\Psi^{h}\neq\Psi_{e}$ sur $\Gamma_{1}$. En imposant exactement sur $\Gamma_{1}$$\Psi^{h}=\Psi_{e}$, on obtiens une erreur encore beaucoup plus faible:

\[-1.24\,10^{-4}\leq\Psi_{e}-\Psi^{h}\leq9.5\,10^{-5}\]

Attention

: sensibilité de la solution numérique aux conditions aux limites à l’infini

2.1.4. Application: étude de l’effet Magnus#

2.1.4.1. Description du modèle#

On reprend l’exemple précédent, en supposant que le cylindre tourne à vitesse angulaire $\omega_{0}\overrightarrow{e}_{z}$. On suppose maintenant que le fluide est visqueux, pour pouvoir être entraîné par le cylindre, mais reste irrotationnel. L’écoulement reste un écoulement à potentiel, vérifiant une équation de Laplace:

\[\Delta\Psi=0\]

Les conditions aux limites sont telles que loin de l’obstacle la vitesse est égale à $U_{0}$. Pour la condition aux limites sur l’obstacle, la condition physique est que la vitesse du fluide doit être égale à la vitesse du cylindre $U_{r}=\omega_{0}R$, Cependant, cette condition aux limites n’est pas compatible avec les hypothèses de fluide irotationnel. On impose donc une condition moins forte: l’égalité de l’intégrale de la vitesse tangentielle sur le cylindre, et la condition d’imperméabilité, soit

\[\oint_{r=R}\overrightarrow{U}.\overrightarrow{dl}=2\pi\omega_{0}R^{2},\,\,\,\,\overrightarrow{U}.\overrightarrow{n}=0\]

La première condition impose la circulation de vitesse $\kappa$ autour du cylindre:

\[\kappa=\oint_{r=R}\overrightarrow{U}.\overrightarrow{dl}=\int_{0}^{2\pi}u_{\theta}(R,\theta)\,d\theta=2\pi\omega_{0}R^{2}\]

qui peut s’interpréter comme la condition qu’en moyenne la vitesse du fluide $u_{\theta}$ sur le cylindre est égale à la vitesse du cylindre $\omega_{0}R$.

La mise en rotation du fluide par le cylindre engendre autour du cylindre une répartition de pression qui n’est plus symétrique par rapport à $Ox$. Cette répartition de pression crée une force de portance sur le cylindre $\Gamma_{1}$ (suivant Oy):

\[F=\int_{\Gamma_{1}}p\,n_{y}ds\]

où $n_{y}$ est la composante suivant y de la normale $\overrightarrow{n}$ à la surface $\Gamma$. La théorie des écoulements potentiels montre que cette force est proportionnelle à la circulation $\kappa$

\[F=-\rho_{0}U_{0}\kappa\]

Cette force de portance, perpendiculaire à la vitesse du fluide loin du cylindre, est appelée effet Magnus et explique les mouvements de lift des balles de tennis ou de golf.

2.1.4.2. Mise en équation#

Pour résoudre ce problème par éléments finis, on calcule une approximation $\Psi^{h}$ de la fonction de courant $\Psi$ solution de:

\[\Delta\Psi=0\,\,\mbox{\,\,\ dans\,}\,\Omega\]

associée aux conditions aux limites:

\[\Psi|_{\Gamma_{0}}=U_{0}y,\,\,\,\Psi|_{\Gamma_{1}}=\Psi_{1}\]

La seconde conditions aux limites implique que le cylindre $\Gamma_{1}$ est une ligne de courant et la valeur $\Psi_{1}\neq0$ imposée est telle que l’on ait une circulation $\kappa$ autour du cylindre.

2.1.4.2.1. solution sans rotation#

\[\Delta\Psi_{1}=0\,\,\mbox{\,\,\ dans\,}\,\Omega\]

associée aux conditions aux limites:

\[\Psi_{1}|_{\Gamma_{0}}=U_{0}y,\,\,\,\Psi_{1}|_{\Gamma_{1}}=0\]

2.1.4.2.2. solution avec rotation pure#

\[\Delta\Psi_{2}=0\,\,\mbox{\,\,\ dans\,}\,\Omega\]

associée aux conditions aux limites:

\[\Psi_{2}|_{\Gamma_{0}}=0,\,\,\,\Psi_{2}|_{\Gamma_{1}}=1\]

2.1.4.2.3. solution générale#

c’est une combinaison linéaire des 2 solutions de base

\[\psi=\psi_{1}+\beta\psi_{2}\]

Pour déterminer $\beta$, il faut une condition supplémentaire. Dans notre cas on va imposer la position du point d’arrêt $P_{a}$ sur le cylindre, ce qui va fixer $\beta$. En effet par définition, la vitesse doit être nulle au point d’arrêt $P_{a}$. Hors $P_{a}$ est sur le cylindre, donc sa vitesse normale est nulle (imperméabilité). La condition impose donc la nullité de la vitesse tangentielle en $P_{a}$:

\[\overrightarrow{U}(P_{a}).\overrightarrow{t}=0\]

En notant $\theta_{a}$ l’angle du point d’arrêt $P_{a}=[R\cos\theta_{a},R\sin\theta_{a}]$, la condition s’écrit en fonction de $\psi$:

\[-\frac{\partial\psi}{\partial y}\sin\theta_{a}-\frac{\partial\psi}{\partial x}\cos\theta_{a}=0\]

ce qui fournit la relation pour calculer $\beta$:

\[\frac{\partial\psi_{1}}{\partial y}\sin\theta_{a}+\frac{\partial\psi_{1}}{\partial x}\cos\theta_{a}=\beta\left(\frac{\partial\psi_{2}}{\partial y}\sin\theta_{a}+\frac{\partial\psi_{2}}{\partial x}\cos\theta_{a}\right)\]

2.1.4.3. Résolution avec COMSOL#

Les différentes étapes de la résolution sont:

choix du modèle physique (EDP) et choix du degré d’approximation (éléments $P^{2}$ par défaut).
création de la géométrie
maillage du domaine de calcul
définition des coefficients du modèle dans $\Omega$
définition des conditions aux limites sur la frontière $\Gamma=\partial\Omega$
résolution
analyse et tracé de la solution

Un exemple de résultat de simulation de l’effet Magnus avec Comsol est accessible sous le lien suivant

Simulation de l’effet Magnus avec COMSOL

2.1.4.4. Analyse du résultat#

Pour analyser le résultat, on calcule tout d’abord le champ de pression en utilisant Bernoulli:

\[p+\frac{1}{2}\rho u^{2}=cste\]

soit en fonction de $\Psi$ (à une constante près)

\[p=-\frac{1}{2}\rho(\left(\frac{\partial\Psi}{\partial x}\right)^{2}+\left(\frac{\partial\Psi}{\partial x}\right)^{2})\]

On en déduit la force de portance $F$

\[F=\int_{\Gamma_{1}}p.n_{y}\,d\Gamma\]

répartition pression avec COMSOL

Avec COMSOL, on obtiens pour $\Psi_{0}=0.4$, $U_{0}=1$, $\rho_{0}=1$ et $R=0.1$ :

\[F=-1.098\]

Le calcul de la circulation $\kappa$ en fonction de $\Psi$ s’écrit:

\[\kappa=\int_{\Gamma_{1}}\left(\frac{\partial\Psi}{\partial x}n_{x}+\frac{\partial\Psi}{\partial y}n_{y}\right)d\Gamma\]

puisque $u_{\theta}=\overrightarrow{U}.\overrightarrow{t}=\frac{\partial\Psi}{\partial x}n_{x}+\frac{\partial\Psi}{\partial y}n_{y}$.

Le calcul COMSOL donne:

\[\kappa=1.090\]

ce qui est en très bon accord avec la théorie (solution exacte), qui pour une circulation positive, prédit une force de portance négative. Compte tenu de cette relation, on peut aussi en déduire la vitesse de rotation $\omega_{0}$ du cylindre pour la valeur de $\Psi_{1}$ imposée:

\[\omega_{0}=17.34\,rd/s\]

2.1.5. Précision d’une simulation numérique#

La précision d’un calcul numérique dépend de plusieurs facteurs, parmi lesquels:

le modèle mathématique
la discrétisation numérique du modèle mathématique
la précision des calculs sur l’ordinateur

Le dernier point est en général négligeable si on prend le soin de faire des calculs en double précision.

L’erreur due au choix du modèle mathématique dépend évidemment du problème physique.

L’erreur de discrétisation dépend du choix de la méthode numérique utilisée. Dans la cas des éléments finis, elle dépend:

du choix de l’interpolation : $P^{1}$, $P^{2}$ …
du choix et de la qualité du maillage $\mathcal{T}^{h}$
de la façon d’imposer les conditions aux limites

Pour des éléments finis, l’erreur de discrétisation est en général majorée par une erreur d’interpolation, qui pour une approximation $P^{q}$ est d’ordre $\theta(h^{1+q})$, où $h$ est la dimension caractéristique de l’élément.

Pour un triangle $e_{k}$, on définit la taille $h_{k}$ comme la longueur du plus grand coté (figure ci-dessous).

triangle

On définit aussi le diamètre $d_{i}$ du cercle inscrit et le diamètre $d_{e}$ du cercle circonscrit, ainsi que le rapport d’aspect du triangle $R=\frac{d_{i}}{d_{e}}$. Pour un triangle équilatéral, on a $d_{e}=h$ et $d_{i}=\frac{\sqrt{2}}{3}h$, donc $R=\frac{\sqrt{2}}{3}$. Par contre pour un triangle de plus en plus aplati, le rapport d’aspect tend vers 0 (car $d_{i}\rightarrow0$) , ce qui indique la dégénérescence du triangle.

Avec ces définitions, on montre que l’erreur d’interpolation $\mathcal{P}^{1}$ sur un triangle vérifie une relation du type:

\[\left\Vert f-f^{h}\right\Vert =\sqrt{\int_{e_{k}}(f-f^{h})^{2}\,dxdy}\,\le\,C\,h_{k}^{2}\,\sqrt{\int_{e_{k}}\left((\frac{\partial^{2}f}{\partial x^{2}})^{2}+(\frac{\partial^{2}f}{\partial y^{2}})^{2}\right)dxdy}\]

Cependant la constante $C$ dépend du rapport d’aspect $R$ et l’estimation d’erreur dégénère lorsque $R\rightarrow0$, i.e. lorsque le triangle devient de plus en plus aplatis.

Sous COMSOL, la qualité du maillage permet d’avoir une idée sur la valeur de ce rapport d’aspect $R$. Pour cela, on calcul sur chaque élément la valeur suivante, fonction de l’aire $A$ et de la longueur des cotés $h_{1},h_{2},h_{3}$:

\[q=\frac{4\sqrt{3}A}{h_{1}^{2}+h_{2}^{2}+h_{3}^{2}}\]

qui varie de 0 (si le triangle est dégénérée) à 1 (pour un triangle équilatéral $h_{1}=h_{2}=h_{3}=h$ et $A=\frac{\sqrt{3}}{4}h^{2}$). Un bon maillage nécessite en générale une valeur $q\ge0,3$.

2.2. Notebook: écoulement potentiel autour d’un cylindre#

Marc Buffat département mécanique, Univerisité Lyon 1 MagnusForce

%matplotlib inline
import numpy as np
import scipy as sp
import matplotlib.pyplot as plt
from matplotlib import rcParams
rcParams['font.family'] = 'serif'
rcParams['font.size'] = 14

2.2.1. Modèle fluide parfait incompressible en 2D#

la solution se calcule avec la fonction de courant $\Psi(x,y)$

\[\begin{split} \overrightarrow{U}=\left[\begin{array}{c} \,\frac{\partial\Psi}{\partial y}\\ -\frac{\partial\Psi}{\partial x} \end{array}\right] \end{split}\]

solution de l’équation de Laplace :

\[ -\Delta \Psi = 0 \]

Les conditions aux limites associées sont :

\[\Psi(r=\infty)=yU_{0}\,\,,\,\,\,\Psi(r=R)=cste \]

En coordonnées polaires, l’équation s’écrit:

\[ \frac{1}{r}\frac{\partial }{\partial r} \left(r \frac{\partial \Psi}{\partial r} \right) + \frac{1}{r^2}\frac{\partial^2 \Psi}{\partial r^2} = 0 \]

avec

\[ U_r = \frac{1}{r} \frac{\partial \Psi}{\partial \theta} \mbox{ et } U_\theta = -\frac{\partial \Psi}{\partial r} \]

et

\[ U_x = U_r \cos\theta - U_\theta \sin \theta \mbox{ et } U_y = U_r \sin\theta + U_\theta \cos \theta\]

La solution analytique peut être obtenue avec la théorie des fonctions conformes comme une combinaison linéaire de solutions élémentaires du type :

solution uniforme
source/puit
vortex
dipôle / doublet

2.2.2. Solutions élémentaires#

2.2.2.1. champ uniforme vitesse constante $U_0$#

\[\Psi_1 = U_0 y\]

d’où le champ de vitesse $$ \overrightarrow{U_1}=\left[\begin{array}{c} U_x=U_0\\ U_y=0 \end{array}\right] $$

2.2.2.2. source d’intensité C à l’origine#

\[\Psi_2 = C \theta\]

d’où le champ de vitesse

\[\begin{split} \overrightarrow{U_2}=\left[\begin{array}{c} U_r=\frac{C}{2\pi r}\\ U_\theta=0 \end{array}\right] \end{split}\]

Le flux de vitesse à travers un cercle de rayon r est constant et vaut :

\[ \int_0^{2\pi} U_r r d\theta = C \]

2.2.2.3. vortex d’intensité $\Gamma$ à l’origine#

\[\Psi_3 = \frac{\Gamma}{2\pi} \ln r\]

d’où le champ de vitesse

\[\begin{split} \overrightarrow{U_3}=\left[\begin{array}{c} U_r=0\\ U_\theta=\frac{\Gamma}{2\pi r} \end{array}\right] \end{split}\]

La circulation de vitesse sur un cercle de rayon r est constante et vaut:

\[ \int_0^{2\pi} U_\theta r d\theta = \Gamma \]

2.2.2.4. doublet à l’origine d’intensité $\kappa$#

\[\Psi_4 = \frac{1}{2\pi}\frac{\kappa}{r} \sin\theta\]

d’où le champ de vitesse

\[\begin{split} \overrightarrow{U_4}=\left[\begin{array}{c} U_r=\frac{1}{2\pi}\frac{\kappa}{r^2} \cos\theta\\ U_\theta= \frac{1}{2\pi}\frac{\kappa}{r^2} \sin\theta \end{array}\right] \end{split}\]

2.2.3. Trace des solutions élémentaires#

calcul des solutions élémentaires avec des fonctions Python
tracer des solutions avec matplotlib
attention aux singularités (à l’origine)

# prise en compte de singularité
def Filtre(U0,U,V):
    '''filtrage du champ vitesse >U0: attention U,V modifié'''
    nx,ny = U.shape
    for i in range(nx):
        for j in range(ny):
            um = np.sqrt(U[i,j]**2 + V[i,j]**2)
            if um>U0:
                U[i,j] = U[i,j]*U0/um
                V[i,j] = V[i,j]*U0/um
    return

# domaine de calcul et maillage (grille) pour le calcul de psi et de la vitesse
L = 2
N = 100
# on elimine le pt singulier x=0,y=0
x = np.linspace(-L,L,N)
y = np.linspace(-L/2,L/2,N)
# pour les lignes de courant psi
X, Y = np.meshgrid(x, y)
# et le champ de vitesse U
XX,YY = np.meshgrid(x[::4],y[::4])

2.2.3.1. champ uniforme#

\[\begin{split} \Psi_1 = U_0 y \mbox{ et } \overrightarrow{U_1}=\left[\begin{array}{c} U_x=U_0\\ U_y=0 \end{array}\right] \end{split}\]

def Psi_uniforme(U0,x,y):
    return U0*y
def Vit_uniforme(U0,x,y):
    return U0,0

U0   = 1.0
Lpsi = np.linspace(-U0*L/2,U0*L/2,21)
Psi1 = Psi_uniforme(U0,X,Y)
U1,V1= Vit_uniforme(U0,XX,YY)      

# tracer
fig, ax = plt.subplots(figsize=(12,6))
ax.contourf(X,Y,Psi1,levels=21)
ax.contour(X,Y,Psi1,levels=Lpsi,colors='r')
ax.quiver(XX,YY,U1,V1)
plt.axis('equal')
plt.title("Champ uniforme");

_images/04b86a07848bab7798275621c85291fbaf44b78ae009b9296c449c6069d3d71b.png

2.2.3.2. Source / puit#

\[\begin{split} \Psi_2 = C \theta \mbox{ et } \overrightarrow{U_2}=\left[\begin{array}{c} U_r=\frac{C}{2\pi r}\\ U_\theta=0 \end{array}\right] \end{split}\]

def Psi_source(C,x,y):
    return C*np.arctan2(x,y)
def Vit_source(C,x,y):
    r2 = x**2+y**2
    u = C*2*np.pi*x/r2
    v = C*2*np.pi*y/r2
    return u,v

C = 0.1
Psi2  = Psi_source(C,X,Y)
U2,V2 = Vit_source(C,XX,YY)
Filtre(U0,U2,V2)

fig, ax = plt.subplots(figsize=(12,6))
ax.contour( X,Y,Psi2,levels=21,colors='r')
ax.contourf(X,Y,Psi2,levels=21)
ax.quiver(XX,YY,U2,V2)
plt.axis('equal')
plt.title("Source");

_images/10114c14a193b2dba0fe4f6d5e355c0411386b388360773cafa901ee52712162.png

2.2.3.3. Vortex#

\[\begin{split} \Psi_3 = \frac{\Gamma}{2\pi} \ln r \mbox{ et } \overrightarrow{U_3}=\left[\begin{array}{c} U_r=0\\ U_\theta=\frac{\Gamma}{2\pi r} \end{array}\right] \end{split}\]

def Psi_vortex(G,x,y):
    r = np.sqrt(x**2+y**2)
    return (G/(2*np.pi))*np.log(r)
def Vit_vortex(G,x,y):
    r2 = x**2 + y**2
    u =  (G/(2*np.pi))*y/r2
    v = -(G/(2*np.pi))*x/r2
    return u,v

G = 2.0
Psi3  = Psi_vortex(G,X,Y)
U3,V3 = Vit_vortex(G,XX,YY)
Filtre(U0,U3,V3)

fig, ax = plt.subplots(figsize=(12,6))
ax.contour( X,Y,Psi3,levels=Lpsi,colors='r')
ax.contourf(X,Y,Psi3,levels=21)
ax.quiver(XX,YY,U3,V3)
plt.axis('equal')
plt.title("Vortex");

_images/7edced43162491ccaf0ec9d55906ea5d77d53a88757a06a59e3ecaa920814f89.png

2.2.3.4. Doublet#

\[\begin{split} \Psi_4 = \frac{1}{2\pi}\frac{\kappa}{r} \sin\theta \mbox{ et } \overrightarrow{U_4}=\left[\begin{array}{c} U_r=\frac{1}{2\pi}\frac{\kappa}{r^2} \cos\theta\\ U_\theta= \frac{1}{2\pi}\frac{\kappa}{r^2} \sin\theta \end{array}\right] \end{split}\]

def Psi_doublet(K,x,y):
    r2  = x**2+y**2
    val = (K/(2*np.pi))*y/r2
    return val
def Vit_doublet(K,x,y):
    r4  = (x**2+y**2)**2
    u =  (K/(2*np.pi))*(x**2-y**2)/r4
    v =  (K/(2*np.pi))*2*x*y/r4
    return u,v

K = 2.0
Psi4  = Psi_doublet(K,X,Y)
U4,V4 = Vit_doublet(K,XX,YY)
Filtre(2*U0,U4,V4)

fig, ax = plt.subplots(figsize=(12,6))
ax.contour( X,Y,Psi4,levels=Lpsi,colors='r')
ax.contourf(X,Y,Psi4,levels=21)
ax.quiver(XX,YY,U4,V4)
plt.axis('equal')
plt.title("Doublet");

_images/e4229696b954ab1b6727660d5074c2f1c0650aec5e428b311f317c688bb45ffc.png

2.2.3.5. combinaison linéaire de 2 solutions#

Ecoulement autour d’un cylindre immobile

Psi = Psi1 - Psi4
U   = U1 - U4
V   = V1 - V4

fig, ax = plt.subplots(figsize=(12,6))
ax.contour( X,Y,Psi,levels=Lpsi,colors='r')
ax.contourf(X,Y,Psi,levels=Lpsi)
ax.quiver(XX,YY,U,V)
plt.axis('equal')
plt.title("Ecoulement autour d'un cylindre immobile");

_images/77703a2601b24a3295f073524d6ef758f419246c06a73da88365d42d9f4c5083.png

2.2.3.6. combinaison linéaire de 3 solutions#

Ecoulement autour d’un cylindre immobile en rotation

Psi = Psi1 - Psi4 + Psi3
U   = U1 - U4 + U3
V   = V1 - V4 + V3

fig, ax = plt.subplots(figsize=(12,6))
ax.contour( X,Y,Psi,levels=Lpsi,colors='r')
ax.contourf(X,Y,Psi,levels=Lpsi)
ax.quiver(XX,YY,U,V)
plt.axis('equal')
plt.title("Ecoulement autour d'un cylindre en rotation");

_images/1923c43441b38562da2da17eeb52584c865b601d11003ed2c097081c270883ed.png

2.2.4. Ecoulement autour d’un cylindre de rayon R fixé#

cas du cylindre immobile
détermination de la bonne combinaison linéaire

R = 0.5
# condition d'arret
u0,v0 = Vit_doublet(1.0,R,0)
alpha = -U0/u0
print("alpha=",alpha)

alpha= -1.5707963267948966

Psi5   = Psi_doublet(alpha,X,Y)
U5, V5 = Vit_doublet(alpha,XX,YY)
Psi = Psi1 + Psi5
U   = U1 + U5
V   = V1 + V5
Filtre(2*U0,U,V)

fig, ax = plt.subplots(figsize=(12,6))
ax.contour( X,Y,Psi,levels=Lpsi,colors='r')
ax.contourf(X,Y,Psi,levels=Lpsi)
ax.quiver(XX,YY,U,V)
plt.axis('equal')
plt.title("Ecoulement autour d'un cylindre immobile de rayon R={}".format(R));

_images/8cd3c33088274b450de7521c8b13b3bb3065ca6107e914000243b354cd31d426.png

2.2.4.1. Analyse de l’écoulement#

solution sur le cylindre (vitesse et $\Psi$)
calcul de la répartition de pression sur le cylindre
bilan des forces
- trainée
- portance

Nc = 101
#Theta = np.linspace(0,2*np.pi,Nc)
Theta = np.linspace(-np.pi,np.pi,Nc)
Xc = R*np.cos(Theta)
Yc = R*np.sin(Theta)
# calcul solution sur le cylindre
Psic1 = Psi_uniforme(U0,Xc,Yc)
Psic2 = Psi_doublet(alpha,Xc,Yc)
Psic  = Psic1 + Psic2
Uc1,Vc1 = Vit_uniforme(U0,Xc,Yc)
Uc2,Vc2 = Vit_doublet(alpha,Xc,Yc)
Uc = Uc1 + Uc2
Vc = Vc1 + Vc2
Ut = (-Uc*Yc + Vc*Xc)/R
Un = ( Uc*Xc + Vc*Yc)/R
# calcul pression
rho = 1.0
Pr  = 0.5*rho*(U0**2 - Uc**2 - Vc**2)

plt.figure(figsize=(10,6))
plt.plot(Theta,Ut,label="$U_t$")
plt.plot(Theta,Un,label="$U_n$")
plt.plot(Theta,Psic,'--',label="$\Psi$")
plt.legend()
plt.xlabel("$\\theta$")
plt.ylabel("Vitesse U")
plt.title("Répartition de Vitesse et $\Psi$ sur le cylindre");

_images/5e6cad2e9d8b25251478b3289634cfd40554f22ae4ade597bcfc94a1590b0f52.png

plt.figure(figsize=(10,6))
plt.plot(Xc[:Nc//2],Pr[:Nc//2],label="$Pr+$")
plt.plot(Xc[Nc//2:],Pr[Nc//2:],'--',label="$Pr-$")
plt.xlabel("$x$")
plt.ylabel("Pression")
plt.legend()
plt.title("Répartition sur le cylindre");

_images/6570c89c45569310c7b121333f098a8a366ce22a7dd59556e62635ed44887eed.png

2.2.5. Cas du cylindre en rotation#

cas du cylindre rayon R
détermination de la bonne combinaison linéaire

# on fixe le point d'arret
theta0 = -np.pi/4
x0 = R*np.cos(theta0)
y0 = R*np.sin(theta0)
u0,v0 = Vit_doublet(alpha,x0,y0)
u0  = u0 + U0
ut0 = (-u0*y0 + v0*x0)/R
# calcul vitesse du vortex
u1,v1 = Vit_vortex(1.0,x0,y0)
ut1 = (-u1*y0 + v1*x0)/R
beta = -ut0/ut1
print("beta=",beta,ut0*2*np.pi*R)

beta= 4.442882938158365 4.442882938158366

Psi6   = Psi_vortex(beta,X,Y)
U6, V6 = Vit_vortex(beta,XX,YY)
Psi = Psi1 + Psi5 + Psi6
U   = U1 + U5 + U6
V   = V1 + V5 + V6
Filtre(2*U0,U,V)
psi0 = Psi_uniforme(U0,x0,y0) + Psi_doublet(alpha,x0,y0) + Psi_vortex(beta,x0,y0)
print("Psi sur le cylindre ",psi0)
Lpsi0=np.linspace(psi0-U0*L/2,psi0+U0*L/2,21)

Psi sur le cylindre  -0.49012907173427345

fig, ax = plt.subplots(figsize=(12,6))
ax.contour( X,Y,Psi,levels=Lpsi0,colors='r')
ax.contourf(X,Y,Psi,levels=Lpsi0)
ax.quiver(XX,YY,U,V)
plt.axis('equal')
plt.title("Ecoulement autour d'un cylindre en rotation de rayon R={}".format(R));

_images/a45e8520ea4c462be2ccdcca0c4f2b4ae59267f7a9d9c4996af7e9945e4599ac.png

2.2.5.1. Analyse de l’écoulement#

calcul de la solution vitesse et $\Psi$ sur le cylindre
calcul de la répartition de pression sur le cylindre
bilan des forces
- trainée
- portance

# solution sur le cylindre
Psic3 = Psi_vortex(alpha,Xc,Yc)
Psic  = Psic1 + Psic2 + Psic3
Uc3,Vc3 = Vit_vortex(beta,Xc,Yc)
Uc = Uc1 + Uc2 + Uc3
Vc = Vc1 + Vc2 + Vc3
Ut = (-Uc*Yc + Vc*Xc)/R
Un = ( Uc*Xc + Vc*Yc)/R
Pr  = 0.5*rho*(U0**2 - Uc**2 - Vc**2)

plt.figure(figsize=(10,6))
plt.plot(Theta,Ut,label="$U_t$")
plt.plot(Theta,Un,label="$U_n$")
plt.plot(Theta,Psic,'--',label="$\Psi$")
plt.legend()
plt.xlabel("$\\theta$")
plt.ylabel("Vitesse U")
plt.title("Répartition de Vitesse et $\Psi$ sur le cylindre en rotation");

_images/a9ad035f1f952a04cb8b07ab9ffe07e70d4d247e0f843ee505fd4790d6a638db.png

plt.figure(figsize=(10,6))
plt.plot(Xc[:Nc//2+1],Pr[:Nc//2+1],label="$Pr+$")
plt.plot(Xc[Nc//2:],Pr[Nc//2:],label="$Pr-$")
plt.legend()
plt.xlabel("$x$")
plt.ylabel("Pression")
plt.title("Répartition sur le cylindre");

_images/0c093ca3583e89a6af32a2dc442be1f47efde9f650b2c9a9bb1db50827626bc5.png

# calcul des forces de pression
Fp =  np.trapz(Pr[:Nc//2+1],Xc[:Nc//2+1])
Fm =  np.trapz(Pr[Nc//2:],  Xc[Nc//2:])
print("portance=",Fp+Fm,"beta=",beta)

portance= 4.439960215340383 beta= 4.442882938158365

2.2.6. FIN#

2.3. Notebook: Application à l’étude d’une TurboVoile#

Marc BUFFAT, département mécanique, UCB Lyon 1

import sympy as sp
import numpy as np
import matplotlib.pyplot as plt

2.3.1. Analyse et paramètres du problème#

$U_0$ vitesse relative du vent d’angle $\beta$
$V_b$ vitesse du bateau
$R$ rayon du cylindre
$H$ hauteur du cylindre
$\omega$ vitesse de rotation
$C_f$ coefficient de traînée du bateau ,$S_m$ surface mouillée
$\rho_0$ masse volumique de l’air et $\rho$ masse volumique de l’eau
$Fp$ force de portance
$Ft$ force de traînée

# parametres du problème
Cf, U0, R, H, omega, rho0, rho, L = sp.symbols('C_f U_0 R H omega rho_0 rho L',positive=True)

2.3.2. Modèle fluide parfait: solution analytique#

cylindre rayon R et de longueur L
solution potentiel

2.3.2.1. Solution potentiel autour cylindre#

solution sans rotation

\[ \psi_1 = U_0 r \sin\theta - U_0\frac{R^2\sin\theta}{r}\]

vortex de circulation $\Gamma$

\[ \psi_2 = \frac{\Gamma}{2\pi} \ln r \]

solution de l’équation de Laplace

\[\Delta \psi = 0\]

Le champ de vitesse est obtenu par

\[ U_r = \frac{1}{r} \frac{\partial \Psi}{\partial \theta} \mbox{ et } U_\theta = -\frac{\partial \Psi}{\partial r} \]

et

\[ U_x = U_r \cos\theta - U_\theta \sin \theta \mbox{ et } U_y = U_r \sin\theta + U_\theta \cos \theta\]

Détermination de la relation entre l’angle du point d’arrêt et la circulation.

# calcul de la circulation en fct de omega
Gamma = 2*sp.pi*R*omega*R
display("Circulation Gamma",Gamma)

'Circulation Gamma'

\[\displaystyle 2 \pi R^{2} \omega\]

r,theta = sp.symbols('r theta')
psi1 = U0*r*sp.sin(theta) - U0*R**2*sp.sin(theta)/r
psi2 = Gamma/(2*sp.pi)*sp.log(r)
display("psi = ",psi1+psi2)

'psi = '

\[\displaystyle - \frac{R^{2} U_{0} \sin{\left(\theta \right)}}{r} + R^{2} \omega \log{\left(r \right)} + U_{0} r \sin{\left(\theta \right)}\]

Ut1 = -psi1.diff(r).subs(r,R)
Ut2 = -psi2.diff(r).subs(r,R)
Ut = Ut1 + Ut2
display("cdt arret Ut=0 soit ",Ut,"=0")

'cdt arret Ut=0 soit '

\[\displaystyle - R \omega - 2 U_{0} \sin{\left(\theta \right)}\]

'=0'

# condition au pt d'arret -> theta en fct de omega
display("Ut=",Ut)
Theta=sp.solve(Ut,theta)[1]
display(sp.Eq(theta,Theta))

'Ut='

\[\displaystyle - R \omega - 2 U_{0} \sin{\left(\theta \right)}\]

\[\displaystyle \theta = - \operatorname{asin}{\left(\frac{R \omega}{2 U_{0}} \right)}\]

Condition au pt d’arrêt sur le cylindre: $\theta$ < 20 degré (analogie profil d’aile)
d’où la vitesse de rotation max du cylindre avec $\sin\theta \approx \theta$

# w_max
theta_m = sp.symbols('theta_m')
omega_m = theta_m*2*U0/R
display("omega max:",omega_m)

'omega max:'

\[\displaystyle \frac{2 U_{0} \theta_{m}}{R}\]

2.3.2.2. Force de portance#

la force de portance pour un cylindre de rayon R et de longueur L est avec ce modèle proportionnelle à la circulation de vitesse $\Gamma$

L’analyse dimensionnelle donne

Fp = rho0*U0*Gamma*L
display("Force de portance ",Fp)

'Force de portance '

\[\displaystyle 2 \pi L R^{2} U_{0} \omega \rho_{0}\]

2.3.3. Analyse du calcul comsol#

simulation pour différentes valeurs de $\theta_1$
calcul répartition de pression autour du cylindre
analyse des résultats

polaire

calcul portance

portance

calcul circulation

d’où la loi $F_p = F(\Gamma)$

calcul de la vitesse de rotation $ \omega$

\[\Gamma = 2 \pi R \omega R\]

2.3.4. Bilan des forces#

turbo_voile

Calcul dans le référentiel du bateau

$U_r$ = vitesse du vent / bateau (correspond à $U_0$)
$V_b$ = vitesse absolue du bateau
$V $ = vitesse absolue du vent d’angle $\alpha$ / $V_b$
$\beta$ = angle d’incidente de U0 dans le repère lié au bateau

Pour les forces :

la force de portance = fonction de Ur et $\beta$
la force de trainée = fonction de Vb

Analyse dimensionnelle

force de portance $$ F_p = F_p(\rho,\omega,U_0,R,H)$$
force de traînée (surface mouillée du bateau) $$ F_t = F_t(\rho,C_f,V_b,S_m)$$

# Force de portance
Fp = rho0*U0*Gamma*H
display("Portance: Fp=",Fp)

'Portance: Fp='

\[\displaystyle 2 \pi H R^{2} U_{0} \omega \rho_{0}\]

# alpha:angle du vent, V vitesse vent, Vb vitesse bateaux 
Cf, alpha = sp.symbols('C_f alpha')
V, Vb = sp.symbols('V V_b')

# angle beta de la vitesse relative
beta = sp.atan2(V*sp.sin(alpha),V*sp.cos(alpha)+Vb)
display("Angle du vent: beta =",beta)

'Angle du vent: beta ='

\[\displaystyle \operatorname{atan}_{2}{\left(V \sin{\left(\alpha \right)},V \cos{\left(\alpha \right)} + V_{b} \right)}\]

# module vitesse relative
Ur = sp.sqrt((V*sp.sin(alpha))**2+(V*sp.cos(alpha)+Vb)**2)
display("Vitesse relative du vent: Ur =",Ur)

'Vitesse relative du vent: Ur ='

\[\displaystyle \sqrt{V^{2} \sin^{2}{\left(\alpha \right)} + \left(V \cos{\left(\alpha \right)} + V_{b}\right)^{2}}\]

display("Angle de la force de portance:",sp.cos(sp.pi/2-beta))

'Angle de la force de portance:'

\[\displaystyle \frac{V \sin{\left(\alpha \right)}}{\sqrt{V^{2} \sin^{2}{\left(\alpha \right)} + \left(V \cos{\left(\alpha \right)} + V_{b}\right)^{2}}}\]

# d'où la force motrice
Fm = sp.cos(sp.pi/2-beta)*Fp.subs({U0:Ur})
display("Force motrice Fm=",Fm)

'Force motrice Fm='

\[\displaystyle 2 \pi H R^{2} V \omega \rho_{0} \sin{\left(\alpha \right)}\]

# force de trainee fonction de la surface mouillée S
S = sp.symbols('S')
Ft = rho*Cf*Vb**2*S/2
display("Force de trainée: Ft=",Ft)

'Force de trainée: Ft='

\[\displaystyle \frac{C_{f} S V_{b}^{2} \rho}{2}\]

2.3.5. Equilibre des forces pour un nbre de cylindres Nc#

On suppose que le bateau n’a pas d’autre moyens de propulsion que la turbovoile

Bilan des forces => portance + trainée = 0
Calcul de la surface mouillée S du bateau (pour la trainée)
Loi de la vitesse V en fonction de omega

# cas avec Nc cylindres
Nc = sp.symbols('N_c')
eq1 = sp.Eq(Nc*Fm,Ft)
display("Equilibre du bateau ",eq1)

'Equilibre du bateau '

\[\displaystyle 2 \pi H N_{c} R^{2} V \omega \rho_{0} \sin{\left(\alpha \right)} = \frac{C_{f} S V_{b}^{2} \rho}{2}\]

# taille du bateau L longeur , La largueur, h tirant d'eau
L, La, h = sp.symbols('L L_a h')
Sm = L*sp.sqrt(La**2+4*h**2)
# surface mouillée
display("Surface mouillée ",Sm)
eq3 = eq1.subs({S:Sm})
# equilibre => Vb en fonction de omega
display("Equilibre du bateau:",eq3)

'Surface mouillée '

\[\displaystyle L \sqrt{L_{a}^{2} + 4 h^{2}}\]

'Equilibre du bateau:'

\[\displaystyle 2 \pi H N_{c} R^{2} V \omega \rho_{0} \sin{\left(\alpha \right)} = \frac{C_{f} L V_{b}^{2} \rho \sqrt{L_{a}^{2} + 4 h^{2}}}{2}\]

2.3.6. Application numérique#

conteneur Enercon L=130 m x La=22.5 m avec h~5 m tirant d’eau
vent V=30 nds (x 1.85 km/h) de travers ($\alpha = \pi/2$)
Nc=4 cylindres
hauteur cylindre H = 2*h
diametre La/5
ordre de grandeur trainee CF=0.00133
$\theta_m$=20 degrés

# parametres
vals = { theta_m:20*np.pi/180, Nc:4, V:30*1.85*1000/3600, H:2*h, R:La/10, Cf:0.00133, alpha:sp.pi/2 , L:130, La:22.5, h:5, 
         rho0:1, rho:1000, sp.pi:np.pi }
# vitesse de rotation en tr/min
N = sp.symbols("N")
# rotation max
Wmax = omega_m.subs(vals).subs({U0:Ur,Vb:0}).subs(vals)
print('omega_max=',Wmax,' N_max=',Wmax*60/(2*np.pi)," tr/min")
# vitesse bateau fct de omega
eq4=eq3.subs(vals).subs(vals).subs(omega,N*2*sp.pi/60)
display("Equation d'équilibre ",eq4)
# vitesse bateau en nds
VB=sp.solve(eq4,Vb)[1]*3600/1000/1.85
display("Vitesse du bateau en nds",sp.Eq(Vb,VB))

omega_max= 4.78349498694742  N_max= 45.6790123456790  tr/min

"Equation d'équilibre "

\[\displaystyle 653.843971028376 \pi N = 2128.58439241318 V_{b}^{2}\]

'Vitesse du bateau en nds'

\[\displaystyle V_{b} = 1.91160137298645 \sqrt{N}\]

plt.rcParams['figure.figsize'] = 10, 10
sp.plot(VB,(N,0,Wmax*60/(2*np.pi)),title="vitesse du bateau pour un vent de travers de V=30 [nds]",
        xlabel="N en [tr/min]",ylabel="$V_b$ en [nds]",lw=4);

_images/db1f38b3e49a04861a8b6916eb1e798174bdb399ff648f539cd94c3aac084155.png

2.3.7. Analyse paramétrique#

On calcule la taille de la turbovoile en fonction de la taille du bateau L variant de 80 a 200m

N = 5
LL = np.linspace(80,200,N)
for i in range(N):
    # parametres
    vals = { theta_m:20*np.pi/180, Nc:1, V:60*1.85*1000/3600, H:2*h, R:La/10, Cf:0.00133, alpha:sp.pi/2 , 
            La:22.5*LL[i]/130, h:5*LL[i]/130, L:LL[i], rho0:1, rho:1000, sp.pi:np.pi }
    # rotation max
    Wmax = omega_m.subs(vals).subs({U0:Ur,Vb:0}).subs(vals)
    # vitesse bateau fct de omega
    eq4=eq3.subs(vals).subs(vals).subs({omega:Wmax})
    # vitesse bateau en nds
    VB=sp.solve(eq4,Vb)[1]*3600/1000/1.85
    #
    val = []
    for p in [L,La,h,R]:
        val.append(p.subs(vals).subs(vals))
    # 
    print("L={:.1f}\t turbovoile  La={:.1f} h={:.1f} R={:.1f} wmax={:.1f} \t vitesse Vb={:.1f} [nds]".
          format(val[0],val[1],val[2],val[3],Wmax,VB))

L=80.0	 turbovoile  La=13.8 h=3.1 R=1.4 wmax=15.5 	 vitesse Vb=12.9 [nds]
L=110.0	 turbovoile  La=19.0 h=4.2 R=1.9 wmax=11.3 	 vitesse Vb=12.9 [nds]
L=140.0	 turbovoile  La=24.2 h=5.4 R=2.4 wmax=8.9 	 vitesse Vb=12.9 [nds]

L=170.0	 turbovoile  La=29.4 h=6.5 R=2.9 wmax=7.3 	 vitesse Vb=12.9 [nds]
L=200.0	 turbovoile  La=34.6 h=7.7 R=3.5 wmax=6.2 	 vitesse Vb=12.9 [nds]

2.3.8. The END#

3. Modélisation de l’écoulement autour d’un profil d’aile#

Ce cas d’application concerne la modélisation numérique de l’écoulement autour d’un profil d’aile. L’objectif est de déterminer les efforts aérodynamiques sur le profil.

On étudie ensuite son application au design d’une pale d’éolienne à axe horizontal.

Il comprend les parties suivantes:

La théorie de l’aile
L’écoulement autour d’un profil d’aile NACA
Le calcul par simulation de l’écoulement autour d’un profil Naca
L’analyse des résultats de la simulation
L’application au design d’une pale d’éolienne en utilisant XFoil

3.1. Théorie de l’aile#

Marc BUFFAT département mécanique, Lyon 1

%matplotlib inline
import numpy as np
import sympy as sp
import matplotlib.pyplot as plt

from metakernel import register_ipython_magics
register_ipython_magics()

3.1.1. Théorie de la portance#

lien vers le sondage
- Questions portance

portance

les différentes théories de la portance (sur internet):

Appuie de l’aile sur l’air (principe de réaction de Newton)
Bernoulli + égalité des temps de parcours
Analyse globale: bilan de quantité de mouvement, déviation de l’écoulement
Théorie de la Circulation de vitesse

3.1.2. Caractéristique aérodynamique d’un profil#

Analyse dimensionnelle:

Etude d’un profil 2D

étude en fonction de l’angle $\alpha$
fonction du nombre de Reynolds
$C_L$ coefficiant de portance: projection suivant $\vec{N}=[-\sin\alpha, \cos\alpha]$ $\perp$ à $\vec{U}_0$

\[ C_L = \int_{C} Pr\, \vec{n}.\vec{N}\,ds \]

\[ C_L = \int_{C} Pr\, ds \mbox{ avec } ds = dx \cos{\alpha} + dy \sin{\alpha} \]

$C_m$ coefficient de moment (moment pression) par rapport à un point de référence $P=[x_{ref}=0.25, y_{ref}=0]$ , qui n’est donc pas nécessairement sur la ligne moyenne à 1/4 de corde

\[ C_m = \int_{C} Pr\, \vec{n} \land \vec{PM} \,ds \]

\[C_m = \int_{C} -Pr \left((x-x_{ref}) dx + (y-y_{ref}) dy \right)\]

on calcule alors le moment en un point Q quelconque par la relation de transport

\[ \mathcal{M}_Q = \mathcal{M}_P + \vec{QP} \land \vec{F}_p \mbox{ avec } \vec{F}_p = C_L [-\sin\alpha, \cos\alpha]\]

d’où la position du centre de poussée exacte / à $P$ $x,y$ (erreur)

\[ \frac{C_m}{C_L} = (x- x_{ref})\cos\alpha + (y- y_{ref})\sin\alpha\]

3.1.3. Analyse par bilan de quantité de mouvement#

Modèle simplifié des efforts exercés par un profil d’aile de longueur (de corde) $L_c$ en incidence d’angle $\alpha$ sur un écoulement stationnaire d’un fluide parfait incompressible de masse volumique $\rho$ et de vitesse $U_1$ (schéma ci-dessous)

portance

On constate expérimentalement que l’effet de l’aile sur l’écoulement, aux faibles incidences, est de dévier les lignes de courant d’un angle $\alpha$ correspondant à l’angle d’incidence du profil, c.a.d l’angle d’inclinaison du bord de fuite.

ligne courant

d’après « How do wings work? » H. Babinsky, Physics Education 38(6) 2003

A 1-minute video released by the University of Cambridge sets the record straight on a much misunderstood concept - how wings lift. I start by giving the wrong explanation and asking who has heard it and every time 95% of the audience puts their hand up. Only a handful will know that it is wrong.Professor Holger BabinskyIt’s one of the most tenacious myths in physics and it frustrates aerodynamicists the world over. Now, …

from IPython.display import YouTubeVideo

YouTubeVideo('e0l31p6RIaY', width=600, height=400)

3.1.3.1. bilan de masse et de quantité de mouvement#

portance

\[ U_1 S_1 = U_2 S_2 \]
\[ p_1 = p_2 = p_0 \]
\[ S_2 = L_c W \]

on en déduit par bilan de quantité de mouvement la force de portance et le coefficient de portance $C_p$

\[ Fp = \rho U_1^2 \sin\alpha W L_c \]

soit à faible incidence une proportionnalité du coefficient de portance $C_p$ avec l’angle d’incidence $\alpha$

\[ C_p = \frac{F_p}{\rho U_1^2 W L_c} \approx \alpha\]

3.1.4. Potentiel complexe#

Pour un écoulement incompressible 2D de fluide parfait irrotationnel le potentiel $\phi(x,y)$ et la fonction de courant $\psi(x,y)$ vérifie l’équation de Laplace:

\[ \Delta \psi = \Delta \phi = 0 \]

dont le champ de vitesse $\vec{U} = [U_x,U_y]$ en cartésien est donné par

\[ U_x = \frac{\partial \phi}{\partial x} = \frac{\partial \psi}{\partial y} \]

\[ U_y = \frac{\partial \phi}{\partial y} = -\frac{\partial \psi}{\partial x} \]

On définit alors un potentiel complexe $\Phi(z)$:

\[\Psi(z) = \phi(x,y) + i \psi(x,y)\]

qui est une fonction analytique en z, dont la dérivée est indépendante de dz.

Dans ce cas $\Phi(z)$ est différentiable en z et vérifie:

(3.1)#\[\begin{eqnarray} \frac{d \Phi}{d z} &=& \frac{\partial \phi}{\partial x} + i \frac{\partial \psi}{\partial x}\\ &=& -i \frac{\partial \phi}{\partial y} + \frac{\partial \psi}{\partial y}\\ &=& U_x - i U_y \\ \end{eqnarray}\]

On le vérifie en prenant respectivement $\delta z = \delta x $ et $\delta z = i\delta y$ et en utilisant les relations de Cauchy-Riemann entre $\phi$ et $\psi$

On en déduit le module de la vitesse:

\[\| U \| = \sqrt{U_x^2+U_y^2} = \left|\frac{d \Psi}{d z}\right|\]

et son angle $\theta$

\[ \theta = -arg\left(\frac{d \Psi}{d z} \right)\]

3.1.4.1. Coordonnées polaires#

En coordonnées polaires, l’équation de Laplace s’écrit:

\[ \frac{1}{r}\frac{\partial }{\partial r} \left(r \frac{\partial \Psi}{\partial r} \right) + \frac{1}{r^2}\frac{\partial^2 \Psi}{\partial r^2} = 0 \]

et le champ de vitesse en polaire

\[ U_r = \frac{\partial \Phi}{\partial r} \mbox{ et } U_\theta = \frac{1}{r} \frac{\partial \Phi}{\partial \theta} \]

\[ U_r = \frac{1}{r} \frac{\partial \Psi}{\partial \theta} \mbox{ et } U_\theta = -\frac{\partial \Psi}{\partial r} \]

soit en coordonnées cartésiennes $$ U_x = U_r \cos\theta - U_\theta \sin \theta \mbox{ et } U_y = U_r \sin\theta + U_\theta \cos \theta$$

3.1.4.2. Potentiel complexe en polaire#

puisque $z = r e^{i\theta}$ on a donc :

(3.2)#\[\begin{eqnarray} \frac{d \Phi}{d z} &=& \frac{\partial \phi}{\partial r} e^{-i\theta} + i \frac{\partial \psi}{\partial r} e^{-i\theta}\\ &=& -\frac{i}{r}\frac{\partial \phi}{\partial \theta}e^{-i\theta} + \frac{1}{r}\frac{\partial \psi}{\partial \theta}e^{-i\theta}\\ &=& U_x -i U_y = (U_r - i U_\theta) e^{-i\theta} \\ \end{eqnarray}\]

avec $$ U_r = U_x \cos \theta + U_y \sin \theta \mbox{ et } U_\theta = -U_x \sin \theta + U_y \cos \theta$$

On en déduit aussi

(3.3)#\[\begin{eqnarray} \frac{\partial \Phi}{\partial r} &=& \frac{\partial \phi}{\partial r} + i \frac{\partial \psi}{\partial r}\\ &=& U_r - i U_\theta \end{eqnarray}\]

3.1.5. Potentiel complexe de l’écoulement autour d’un cylindre#

On a vu que l’écoulement potentiel autour d’un cylindre de rayon $R$ en rotation est la somme d’un écoulement uniforme, d’un doublet et d’un vortex de circulation $\Gamma$

Pour une vitesse $U_0$ d’incidence $\alpha$ , le potentiel complexe s’écrit

\[\Psi(z) = U_0 z e^{-i\alpha} + U_0 e^{i\alpha} \frac{R^2}{z} - i \frac{\Gamma}{2\pi}\log(\frac{z e^{-i\alpha}}{R})\]

et la vitesse complexe $W = U- i V$

\[ W = U_0 e^{-i\alpha} - U_0 e^{i\alpha} \frac{R^2}{z^2} - i \frac{\Gamma}{2\pi z}\]

# parametres
R, U0 = sp.symbols('R U_0',real=True, positive=True)
Gamma,alpha = sp.symbols('Gamma alpha',real=True)
# variables
theta,x,y = sp.symbols('theta x y',real=True)
r         = sp.symbols('r',real=True,Positive=True)
# fonctions pour manipuler les complexes
from sympy import re,im,I

Pour une vitesse $U_0$ d’incidence $\alpha$ , le potentiel complexe s’écrit

\[\Psi(z) = U_0 z e^{-i\alpha} + U_0 e^{i\alpha} \frac{R^2}{z} - i \frac{\Gamma}{2\pi}\log(\frac{z e^{-i\alpha}}{R})\]

# potentiel complexe Phi
z = sp.symbols('z')
Phi = U0*z*sp.exp(-I*alpha) + U0*sp.exp(I*alpha)*R**2/z \
    - (I*Gamma/(2*sp.pi))*sp.log(z*sp.exp(-I*alpha)/R)
display("Phi(z)=",Phi)

'Phi(z)='

\[\displaystyle - \frac{i \Gamma \log{\left(\frac{z e^{- i \alpha}}{R} \right)}}{2 \pi} + \frac{R^{2} U_{0} e^{i \alpha}}{z} + U_{0} z e^{- i \alpha}\]

la vitesse complexe $W = U- i V$ est alors donnée par $$ W = \frac{d \Phi}{d z}$$

# calcul de la vitesse complexe W
W = 0
### BEGIN SOLUTION
W = sp.diff(Phi,z)
### END SOLUTION
display("W(z)=",W)

'W(z)='

\[\displaystyle - \frac{i \Gamma}{2 \pi z} - \frac{R^{2} U_{0} e^{i \alpha}}{z^{2}} + U_{0} e^{- i \alpha}\]

# calcul fonction de courant psi (en fonction de (x,y)
psi = 0
### BEGIN SOLUTION
psi = im(Phi.subs(z,x+sp.I*y))
### END SOLUTION
display('psi(x,y)=',psi)

'psi(x,y)='

\[\displaystyle - \frac{\Gamma \log{\left(\sqrt{\frac{x^{2}}{R^{2}} + \frac{y^{2}}{R^{2}}} \right)}}{2 \pi} + R^{2} U_{0} \operatorname{im}{\left(\frac{e^{i \alpha}}{x + i y}\right)} - U_{0} x \sin{\left(\alpha \right)} + U_{0} y \cos{\left(\alpha \right)}\]

# en déduire les composantes de vitesse dans Ux et Uy
Ux = 0
Uy = 0
### BEGIN SOLUTION
Ux =  re(W.subs(z,x+sp.I*y))
Uy = -im(W.subs(z,x+sp.I*y))
### END SOLUTION
display("Ux=",Ux)
display("Uy=",Uy)

'Ux='

\[\displaystyle - \frac{\Gamma y}{2 \pi \left(x^{2} + y^{2}\right)} - R^{2} U_{0} \operatorname{re}{\left(\frac{e^{i \alpha}}{\left(x + i y\right)^{2}}\right)} + U_{0} \cos{\left(\alpha \right)}\]

'Uy='

\[\displaystyle \frac{\Gamma x}{2 \pi \left(x^{2} + y^{2}\right)} + R^{2} U_{0} \operatorname{im}{\left(\frac{e^{i \alpha}}{\left(x + i y\right)^{2}}\right)} + U_{0} \sin{\left(\alpha \right)}\]

3.1.5.1. Visualisation des champs#

utilisation d’une technique de masque pour éliminer l’écoulement à l’intérieur du cercle

# parametres numeriques: choix arbitraire de Gamma et alpha ! 
R0   = 0.5
Uinf = 1.0
vals = [(Gamma,-2*np.pi*R*U0),(alpha,np.deg2rad(10)),(R,R0),(U0,Uinf)]
#vals = [(Gamma,0),(alpha,np.deg2rad(10)),(R,R0),(U0,Uinf)]
display("parametres:",vals)

'parametres:'

[(Gamma, -6.28318530717959*R*U_0),
 (alpha, 0.17453292519943295),
 (R, 0.5),
 (U_0, 1.0)]

# domaine de calcul et maillage (grille) pour le calcul de psi et de la vitesse
L = 3
N = 201
pas = 8 # pas pour les vitesses
# points de calcul
xg = np.linspace(-L,L,N)
yg = np.linspace(-L/2,L/2,N)
# pour les lignes de courant psi
X, Y = np.meshgrid(xg, yg)
# utilisation d'un masque
Z = X + 1j*Y
Z = np.ma.masked_where(np.absolute(Z)<0.95*R0,Z)
X = Z.real
Y = Z.imag
# et le champ de vitesse U
XX = X[::pas,::pas]
YY = Y[::pas,::pas]

# calcul des valeurs numeriques
psi1 = sp.lambdify([x,y],psi.subs(vals),'numpy')
u1   = sp.lambdify([x,y],Ux.subs(vals),'numpy')
v1   = sp.lambdify([x,y],Uy.subs(vals),'numpy')
with np.errstate(divide='ignore'):
    Psi1 = psi1(X,Y)
    Psi0 = psi1(0.5,0)
    Levs = np.linspace(Psi0-L/5,Psi0+L/5,21)
    U1   = u1(XX,YY)
    V1   = v1(XX,YY)

<lambdifygenerated-1>:2: RuntimeWarning: invalid value encountered in sqrt
  return -0.17364817766693*x + 0.984807753012208*y + 1.5707963267949*log(2.0*sqrt(x**2 + y**2))/pi + 0.25*imag(exp(0.174532925199433*1j)/(x + 1j*y))

# tracer
fig, ax = plt.subplots(figsize=(12,6))
#ax.contourf(X,Y,Psi1,levels=21)
ax.contour(X,Y, Psi1,levels=Levs,colors='r')
ax.quiver(XX,YY,U1,V1)
cercle = plt.Circle((0.,0.),R0,color='yellow')
ax.add_artist(cercle)
plt.axis('equal')
plt.axis('off')
plt.title("Ecoulement autour du cylindre");

_images/8a2d99f674ee59706fd32e294a705e7275828f119a95b35e27c3eb2de62ff76e.png

3.1.6. Transformation conforme#

On définit une transformation conforme $Z=F(z)$ qui préserve les angles dans le plan complexe.

3.1.6.1. Transformation de Joukovski#

\[ Z = z + \frac{C^2}{z} \]

avec $ z = x + i y $ (plan d’origine) et $Z = X + i Y$ (plan transformé)

C = sp.symbols('C',real=True,positive=True)
F = lambda z: z + C**2/z
display("Z=F(z)",F(z))

'Z=F(z)'

\[\displaystyle \frac{C^{2}}{z} + z\]

3.1.6.2. transformation conforme#

transformation de l’écoulement autour d’un cercle de centre $x0,y0$ de rayon $R=R_0$

$c$ paramétre de la transformation $0\le c\le 1$
$\beta$ angle du point centre / horizontal

# selection cas 1,2,3,4 (choix 4 par defaut)
cas  = 4
beta = 0
if cas==0:
    # cercle
    x0 = 0
    y0 = 0
    R0 = float(R.subs(vals))
    c  = float(R.subs(vals)*0.0)
elif cas==1:
    # parametres pour une ellipse
    x0 = 0
    y0 = 0
    R0 = float(R.subs(vals))
    c  = float(R.subs(vals)*0.9)
elif cas==2:
    # parametre pour une plaque
    x0 = 0
    y0 = 0
    R0 = float(R.subs(vals))
    c  = float(R.subs(vals)*1.0)
elif cas==3:
    # parametres pour un profil symetrique
    R0 = float(R.subs(vals))
    c  = float(R.subs(vals)/1.2) 
    x0 = c - R0
    y0 = 0
else:
    # parametres pour un profil cambre
    beta = np.deg2rad(8)
    R0 = float(R.subs(vals))
    c  = float(R.subs(vals)/1.2)
    x0 = c - R0*np.cos(beta)
    y0 = R0*np.sin(beta)

# calcul de la transformation en coordonnées cartésienne
xc = re(F(x+sp.I*y))
yc = im(F(x+sp.I*y))
display("xc,yc=",xc,yc)

'xc,yc='

\[\displaystyle \frac{C^{2} x}{x^{2} + y^{2}} + x\]

\[\displaystyle - \frac{C^{2} y}{x^{2} + y^{2}} + y\]

# transformation du cercle
Xc = sp.lambdify([x,y],xc.subs(C,c),'numpy')
Yc = sp.lambdify([x,y],yc.subs(C,c),'numpy')
Theta = np.linspace(0,2*np.pi,410)
XC = Xc(x0+R0*np.cos(Theta),y0+R0*np.sin(Theta))
YC = Yc(x0+R0*np.cos(Theta),y0+R0*np.sin(Theta))

plt.figure(figsize=(8,6))
plt.plot(XC,YC)
plt.title("Transformation de Joukovski")
plt.axis('equal');

_images/729f9e3975e457cf633efacbc6b3f19cd4065943fec2c8fef918d9a671e94b9b.png

3.1.6.3. calcul de l’écoulement par transformation de Joukovski#

calcul pour un profil
transformation numérique

def TransJ(Z,lam):
    '''transformation Joukovski (complexe)'''
    return Z + lam**2/Z
def Cercle(C0,R):
    '''pts du cercle de centre C0 (complexe) et de rayon R'''
    Theta = np.linspace(0,2*np.pi,200)
    return C0 + R*np.exp(1j*Theta)

def solutionJ(PHI,lam,C0,R,xg,yg):
    '''calcul de la solution par transformation de Joukovski PHI'''
    X,Y  = np.meshgrid(xg,yg)
    Z    = X+1j*Y
    Z    = np.ma.masked_where(np.absolute(Z-C0)<=R,Z)
    Zc   = Z - C0
    #Phiz = PHI(Z)
    # BUG numpy: calcul 
    # np.log(Z) cannot be calculated correctly due to a numpy bug np.log(MaskedArray);
    Phiz = Zc.copy()
    with np.errstate(divide='ignore'):        
        for m in range(Zc.shape[0]):
            for n in range(Zc.shape[1]):
                Phiz[m,n] = PHI(Zc[m,n])
    # Joukovski transformation 
    J  = TransJ(Z, lam)  
    cercle = Cercle(C0, R)
    airfoil= TransJ(cercle, lam) 
    return J, Phiz.imag, airfoil

#valeur des  parametres
display("parametres ",vals)
alpha0 = np.rad2deg(float(alpha.subs(vals)))
print("alpha={} deg. beta={} deg. c={} C0={}".format(alpha0,np.rad2deg(beta),c,x0+1j*y0))

'parametres '

[(Gamma, -6.28318530717959*R*U_0),
 (alpha, 0.17453292519943295),
 (R, 0.5),
 (U_0, 1.0)]

alpha=10.0 deg. beta=8.0 deg. c=0.4166666666666667 C0=(-0.0784673677041185+0.06958655048003272j)

display(Phi)
# calcul du potentiel complexe 
PhiJ = Phi.subs(vals)
display("Phi(z)=",PhiJ)
# conversion fonction python
PHI  = sp.lambdify(z,PhiJ)

\[\displaystyle - \frac{i \Gamma \log{\left(\frac{z e^{- i \alpha}}{R} \right)}}{2 \pi} + \frac{R^{2} U_{0} e^{i \alpha}}{z} + U_{0} z e^{- i \alpha}\]

'Phi(z)='

\[\displaystyle 1.0 z e^{- 0.174532925199433 i} + \frac{1.5707963267949 i \log{\left(2.0 z e^{- 0.174532925199433 i} \right)}}{\pi} + \frac{0.25 e^{0.174532925199433 i}}{z}\]

# calcul solution ! pble mauvaise valeur de la circulation 
J, Psi, airfoil = solutionJ(PHI,c,x0+1j*y0,R0,xg,yg)
# tracer
fig=plt.figure(figsize=(12,8))
ax=fig.add_subplot(111)
# this means that the flow is evaluated at Juc(z) since c_flow(Z)=C_flow(csi(Z))
cp=ax.contour(J.real, J.imag, Psi,levels=Levs, colors='blue', linewidths=1,
                linestyles='solid')
ax.plot(airfoil.real, airfoil.imag)
ax.fill(airfoil.real, airfoil.imag,color='yellow')
plt.axis('off')
plt.title("Problème: eclt au bord de fuite car circulation non correcte")
ax.set_aspect('equal');
#plt.savefig('eclt_cylindre.png')

_images/8d26746706d09ceaca5393f82cd865f34dcf285b443c622c75c8be20a794ebbd.png

3.1.6.4. circulation de Kutta-Joukovsky#

condition de Kutta au bord de fuite !

\[ \Gamma = -4 \pi U_0 R \sin{(\alpha + \beta)}\]

# valeur de Joukovski
GammaJ = -4*sp.pi*U0*R*sp.sin(alpha+beta)
display("GammaJ=",GammaJ)
# potentiel Joukovsky
PhiJ = Phi.subs(Gamma,GammaJ).subs(vals)
display("PhiJ=",PhiJ)
# conversion fonction python
PHI  = sp.lambdify(z,PhiJ)

'GammaJ='

\[\displaystyle - 4 \pi R U_{0} \sin{\left(\alpha + 0.139626340159546 \right)}\]

'PhiJ='

\[\displaystyle 1.0 z e^{- 0.174532925199433 i} + 0.309016994374947 i \log{\left(2.0 z e^{- 0.174532925199433 i} \right)} + \frac{0.25 e^{0.174532925199433 i}}{z}\]

# calcul solution
J, Psi, airfoil = solutionJ(PHI,c,x0+1j*y0,R0,xg,yg)
# tracer
fig=plt.figure(figsize=(12,8))
ax=fig.add_subplot(111)
# this means that the flow is evaluated at Juc(z) since c_flow(Z)=C_flow(csi(Z))
cp=ax.contour(J.real, J.imag, Psi,levels=Levs, colors='blue', linewidths=1,
                linestyles='solid')
ax.plot(airfoil.real, airfoil.imag)
ax.fill(airfoil.real, airfoil.imag,color='yellow')
plt.axis('off')
ax.set_aspect('equal');

_images/a9c1cb9ff1ed5bedaaec6d52a465cbea1ffe54b373a1e241971f0d5bceb9cd52.png

3.1.7. Calcul de la Portance#

portance

3.1.7.1. Condition de Kutta - Joukovsky#

création circulation vitesse $\Gamma$ telle que la vitesse reste parallèle au bord de fuite (pas de contournement)

condition de kutta-joukovski

circulation

\[\Gamma = 4 \pi U_0 R_0 \sin(\alpha + \beta)\]

portance

\[ F_p = \rho \Gamma U_0 \]

en faible incidence

\[ F_p \propto (\alpha + \beta)\]

3.1.7.2. Expérience de Prandtl (~1918)#

création de la circulation

from IPython.display import YouTubeVideo

YouTubeVideo('VcggiVSf5F8', width=600, height=400)

3.1.8. Annexe 1: Joukovski-airfoil#

Ce notebook est inspiré d’un notebook de calcul de l’écoulement autour d’un profil de Joukovski

3.1.9. FIN#

3.2. Écoulement autour d’un profil NACA#

Écoulement autour d’un profil

Bilan des forces sur un avion

3.2.1. Théorie de l’aile#

Les différentes (fausses) théories

théorie du temps de transit (Bernoulli)
théorie effet Venturi

Oui mais l’écoulement autour d’une plaque plane génère une portance !!

3.2.1.1. théorie de kutta-joukowsky#

Génération d’une circulation de vitesse $\Gamma=\oint\,\overrightarrow{u}.\overrightarrow{dl}$ autour du profil pour que l’écoulement autour du bord de fuite ne le contourne pas (condition de Kutta)!!!.

_images/solpot0.png — Fig. 3.1 solution sans condition de Kutta Joukowsky#

_images/solpot1.png — Fig. 3.2 solution avec condition de Kutta Joukowsky#

Cette circulation $\Gamma$ s’écrit pour un profil symétrique ($a$ rayon du cercle générateur):

\[\Gamma=4\pi U_{0}a\sin\alpha\approx4\pi U_{0}a\alpha\]

_images/TheoTourbillonInitiateur.png — Fig. 3.3 génération du tourbillon de bout d’aile#

_images/TheoTourbillonNouveau.png — Fig. 3.4 tourbillon de bout d’aile#

3.2.1.2. Force de portance#

Théorie de la portance

\[F_{p}=\rho U_{0}\Gamma\]

Analyse dimensionnelle

\[\frac{F_{p}}{\frac{1}{2}\rho U_{0}^{2}L}=f(R_{e},\alpha)\]

3.2.2. Définition du profil#

On veut étudier l’écoulement autour d’un profil NACA pour calculer la portance du profil en fonction de l’angle d’incidence $\alpha$.

3.2.2.1. Définition d’un profil NACA#

Un profil NACA est un profil d’aile standardisé définie à partir de 3 paramètres:

l’épaisseur relative $t$
la cambrure maximale $m=y_{a}$ au point $p=x_{a}$ de la ligne moyenne

3.2.2.1.1. profil symétrique#

Dans ce cas $y_{a}=x_{a}=0$ et le seul paramètre est l’épaisseur relative $t$.

L’extrados est donné par la courbe $y_{e}=F(x)$ suivante et l’intrados par la courbe symétrique $y_{i}=-F(x)$ pour $0\leq x\leq L_{c}$

\[F(x)=5.0\,t\left(0.2969\,\sqrt{x}-0.1260\,x-0.3516\,x^{2}+0.2843\,x^{3}-0.1036\,x^{4}\right)\,\]

La longueur de la corde $L_{c}$ suivant $x$ est donnée par la racine de $F(L_{c})=0$, i.e. $L_{c}=1$ (calcul sans dimension);

Un profil NACA symétrique à 4 chiffres est donné par la nomenclature NACA00yz où les 2 chiffres $yz$ représente l’épaisseur relative du profil $t=yz/100$. Les profils NACA suivants sont des profils symétriques: NACA0009 et NACA0012.

3.2.2.1.2. profil cambré#

Dans ce cas on donne la ligne moyenne du profil en fonction de $x_{a}$ et $y_{a}$. Cette ligne moyenne a pour équation en fonction de $x$

\[\begin{split}y_{m}=\left\{ \begin{array}{l} y_{a}\frac{(2x_{a}-x)\,x}{x_{a}^{2}}\,\mbox{ si }0\leq x\leq x_{a}\\ y_{a}\frac{(1-x)(1+x-2x_{a})}{(1-x_{a})^{2}}\mbox{ si }x\leq x\leq L_{c} \end{array}\right.\end{split}\]

L’équation de l’extrados $\{y_{e},x_{e}\}$ est obtenu à partir de cette ligne moyenne, en introduisant l’angle $\theta$

\[\theta(x)=\frac{dy_{m}}{dx}\]

et la fonction $F(x)$précédente:

\[\begin{split}\begin{align} x_{e}(x) & = x-F(x)\,\sin(\theta(x))\\ y_{e}(x) & = y_{m}+F(x)\,\cos(\theta(x))\,\,\,\,\mbox{ pour }0\leq x\leq L_{c}\end{align}\end{split}\]

L’équation de l’intrados $\{x_{i},y_{i}\}$ correspond à l’équation symétrique en $y$

\[\begin{split}\begin{align} x_{i}(x) & = x+F(x)\,\sin(\theta(x))\\ y_{i}(x) & = y_{m}-F(x)\,\cos(\theta(x))\,\,\,\,\mbox{ pour }0\leq x\leq L_{c}\end{align}\end{split}\]

La longueur de la corde $L_{c}$ vérifie $y_{e}(L_{c})=y_{i}(L_{c})$

Un profil NACA cambré à 4 chiffres est donnée par la nomenclature NACAwxyz où les 2 chiffres $yz$ représentent l’épaisseur relative du profil ($t=yz/100$), le chiffre $w$ la courbure relative en % ($y_{a}=w/100$) et le chiffre $x$ la position de la flèche de la ligne moyenne en $1/10$ de corde ($x_{a}=x/10$). Les profils NACA suivants sont des profils cambrés : NACA6512, NACA6412, NACA2412, NACA4412, NACA4415.

Un profil NACA cambré à 5 chiffres est donnée par la nomenclature NACAwxxyz où les 2 chiffres $yz$ représentent l’épaisseur relative du profil ($t=yz/100$), le chiffre $w$ la courbure relative en % ($y_{a}=w/100$) et les 2 chiffres chiffres $xx$ la position de la flèche de la ligne moyenne en $2/100$ de corde $(x_{a}=2xx/100$). Les profils NACA suivants sont des profils cambrés à 5 chiffres: NACA23112 et NACA16020.

3.2.2.2. Modélisation de l’écoulement autour d’un profil NACA symétrique#

L’écoulement autour d’un profil symétrique en incidence est représenté sur la figure ci-dessous:

_images/schema.png — Fig. 3.5 écoulement autour d’un profil symétrique#

L’écoulement avec incidence positive présente un point d’arrêt sous le bord d’attaque, et un point d’arrêt sur le bord de fuite. La répartition de pression est donnée sur la figure ci-dessous ou on a tracé le coefficient de pression $C_{p}$ :

\[C_{p}=\frac{p-p_{0}}{\frac{1}{2}\rho U_{0}^{2}}\]

en fonction de la corde. On voit qu’à l’extrados il est négatif (zone de dépression) et positif à l’intrados (zone de surpression) sauf tout près du bord d’attaque où il y a une forte dépression. En l’absence de frottement visqueux, la force de pression qui s’exerce sur le profil est une force de portance.

_images/cp.png — Fig. 3.6 coefficient de portance $C_p$#

Pour modéliser cet écoulement, on utilise un modèle potentiel. Le champ de vitesse $\overrightarrow{U}$ découle d’un potentiel $\Phi$ , et aussi d’une fonction de courant $\Psi$. Les conditions aux limites sont:

à l’infini (frontière $\Gamma_{0}$), la vitesse est égale à $U_{0}$ avec un angle d’incidence $\alpha$: $\overrightarrow{U}=U_{0}[\cos\alpha,\sin\alpha]$
sur le profil (frontière $\Gamma_{1}$), la vitesse est tangente au profil (glissement): $\overrightarrow{U}.\overrightarrow{n}=0$

Mais sur le profil, on doit imposer une condition supplémentaire, qui est la circulation de vitesse $\kappa$ autour du profil $\Gamma_{1}$:

\[\kappa=\oint_{\Gamma_{1}}\overrightarrow{U}.\overrightarrow{dl}\]

C’est cette circulation qui crée la portance. La valeur de cette circulation $\kappa$ doit être telle que l’écoulement au bord de fuite soit régulier (i.e. que le point d’arrêt soit sur le bord de fuite). C’est la condition de Kutta Joukovski. Sur les figures ci-dessous, on a représenté à gauche une solution potentielle avec $\kappa=0$ (non physique) et à droite la solution avec la valeur de $\kappa$ satisfaisant la condition de Kutta Joukovski.

3.2.2.3. Calcul du point d’application#

_images/portance.jpg — Fig. 3.9 Portance sur un profil#

Le point d’application $P\left(x_{P},y_{P}\right)$ (aussi appelé centre de poussée) de la résultante des forces aérodynamiques $\overrightarrow{F}=\overrightarrow{F_{p}}+\overrightarrow{F_{t}}$ doit vérifier la relation suivante sur le moment de la résultante par rapport à un point $O$ quelconque:

\[\overrightarrow{OP}\wedge\overrightarrow{F}=\int_{S}\overrightarrow{OM}\wedge\overrightarrow{dF}\]

avec $$\overrightarrow{F}={\displaystyle \int_{S}-p\overrightarrow{n}dS} \mbox{ et } \overrightarrow{dF}=-p\overrightarrow{n}dS$$

$M\left(x_{M},y_{M}\right)$ est un point du profil et $\overrightarrow{n}$ le vecteur normal sortant. En appliquant cette relation, on obtient l’équation de la droite qui porte la résultante, ce qui n’est pas suffisant pour définir la position du centre de poussée. On calcul donc sa position sur la ligne moyenne.

3.3. TP Écoulement autour d’un profil NACA avec COMSOL#

On considère l’écoulement potentiel autour d’un profil NACA en incidence (angle $\alpha$).

profil NACA en incidence

Pour effectuer la simulation sous COMSOL, il faut générer un maillage autour du profil. Pour cela on doit définir le profil NACA à l’aide d’un fichier de points au format DXF (CAO).

3.3.1. Définition d’un profil NACA#

Un profil NACA est un profil d’aile standardisé défini à partir de 3 paramètres:

l’épaisseur relative $t$
la cambrure relative $y_{a}$ au point $x_{a}$de la ligne moyenne

3.3.1.1. profil symétrique#

Dans ce cas $y_{a}=x_{a}=0$ et le seul paramètre est l’épaisseur relative $t$. L’extrados est donné par la courbe sans dimension $y_{e}==F(x)$ suivante et l’intrados par la courbe symétrique $y_{i}=-F(x)$. Les longueurs sont adimensionalisées par la longueur de corde $L_c$ du profil:

\[F(x)=5.0\,t\left(0.2969\,\sqrt{x}-0.1260\,x-0.3516\,x^{2}+0.2843\,x^{3}-0.1015\,x^{4}\right)\mbox{ pour }0\leq x\leq 1\]

La fonction $F(x)$ ne vérifie pas exactement $F(1)=0$ et donc une épaisseur nulle au bord de fuite. Pour assurer une épaisseur nulle au bord de fuite pour la simulation numérique, on modifie le dernier coefficient par $-0.1036$, ce qui induit un petit changement de la surface portante.

\[F(x)=5.0\,t\left(0.2969\,\sqrt{x}-0.1260\,x-0.3516\,x^{2}+0.2843\,x^{3}-0.1036\,x^{4}\right)\mbox{ pour }0\leq x\leq 1\]

Un profil NACA symétrique à 4 chiffres est donnée par la nomenclature NACA00yz où les 2 chiffres $yz$ représente l’épaisseur relative du profil $t=yz/100$.

3.3.1.2. profil cambré Naca à 4 chiffres#

Dans ce cas on donne la ligne moyenne du profil en fonction de $x_{a}$ et $y_{a}$. Cette ligne moyenne a pour équation en fonction de $x$

\[\begin{split}y_{m}=\left\{ \begin{array}{l} \frac{y_{a}(2x_{a}-x)\,x}{x_{a}^{2}}\mbox{ si }0\leq x\leq x_{a}\\ \frac{y_{a}(1-x)(1+x-2x_{a})}{(1-x_{a})^{2}}\mbox{ si }x_{a}\leq x\leq 1 \end{array}\right.\end{split}\]

L’équation de l’extrados $\{y_{e},x_{e}\}$ est obtenu à partir de cette ligne moyenne, en introduisant l’angle $\theta$

\[\theta(x)=\frac{dy_{m}}{dx}\]

et la fonction $F(x)$précédente:

\[\begin{split}\begin{aligned} x_{e}(x) & = & x-F(x)\,\sin(\theta(x))\\ y_{e}(x) & = & y_{m}+F(x)\,\cos(\theta(x))\mbox{ pour }0\leq x\leq 1\end{aligned}\end{split}\]

L’équation de l’intrados $\{x_{i},y_{i}\}$ correspond à l’équation symétrique en $y$

\[\begin{split}\begin{aligned} x_{i}(x) & = & x+F(x)\,\sin(\theta(x))\\ y_{i}(x) & = & y_{m}-F(x)\,\cos(\theta(x))\mbox{ pour } 0\leq x\leq 1\end{aligned}\end{split}\]

Un profil NACA cambré à 4 chiffres est donnée par la nomenclature NACAwxyz où les 2 chiffres $yz$ représentent l’épaisseur relative du profil ($t=yz/100$), le chiffre $w$ la courbure relative en % ($y_{a}=w/100$) et le chiffre $x$ la position de la flèche de la ligne moyenne en $1/10$ de corde ($x_{a}=x/10$).

3.3.1.3. profil cambré Naca 5 chiffres#

Un profil NACA cambré à 5 chiffres est donnée par la nomenclature NACAabcyz où les 2 derniers chiffres $yz$ représentent l’épaisseur relative du profil ($t=yz/100$), et abc vaut 210,220,230,240,250 pour des cambrure simples et 211,221,231,241,251 pour des cambrures doubles. Les équations sont les mêmes sauf pour l’équation de la ligne moyenne $y_{m}$ que l’on la trouvera sur le site Wikipédia ici :

https://fr.wikipedia.org/wiki/Profil_NACA

3.3.2. Format CAO DXF#

DXF, sigle de Drawing eXchange Format, est un format créé par la société Autodesk servant à échanger des fichiers DAO ou CAO entre systèmes CAO n’utilisant pas le même format de fichier natif. Il a été conçu à l’origine pour représenter les modèles 3D créés avec AutoCAD.

Pour décrire une courbe, on la décrit par une liste de segments en donnant les coordonnées des 2 points extrémités et décrit dans le sens trigonométrique.

Dans l’exemple en annexe`, on décrit un triangle en définissant ces 3 cotés après les entêtes LINE et les codes 10,20,30 pour les coordonnes x,y,z du 1er point et 11,21,31 pour les coordonnées du 2nd point. Toutes les autres entêtes doivent être conservées.

    # programme python pour l'écriture d'une courbe au format DXF dans un fichier
    # la courbe est définie par les pts X,Y,Z 
    # attention si la courbe est fermée, le dernier pt = premier pt 
    def writedxf(nom,X,Y,Z):     
        fid=open(nom,'w')     
        fid.write('999\nCourbe au format DXF\n 0\nSECTION\n 2\nENTITIES\n 0\n')
        for i in range(len(X)-1):       
          fid.write('LINE\n 8\n 0\n')       
          # ecriture du segment Pi-Pi+1         
          fid.write('10\n %.4f\n 20\n %.4f\n 30\n %.4f\n'%(X[i],Y[i],Z[i]))
          fid.write('11\n %.4f\n 21\n %.4f\n 31\n %.4f\n'%(X[i+1],Y[i+1],Z[i+1]))
          fid.write(' 0\n')
        fid.write('ENDSEC\n 0\nEOF\n')     
        fid.close()     
        return 

3.3.3. Simulations avec COMSOL#

L’objectif de l’étude est de simuler avec COMSOL l’écoulement autour du profil pour déterminer la force de portance et son point d’application pour différentes incidences $-10^{\circ}\leq\alpha\leq20^{\circ}$.

En séance on traitera le cas d’un profil symétrique. Puis vous aurez à étudier le profil NACA qui vous a été attribué.

3.3.3.1. Préparation#

Pour effectuer la simulation avec COMSOL, il faut créer un fichier permettant de définir la géométrie du profil et mailler le domaine de calcul, et un second fichier pour extraire la répartition de pression sur le profil. Pour valider la simulation, on a aussi besoin de données de comparaison.

En utilisant le profil naca nacawxyz qui vous a été attribué, vous devez, en utilisant un programme python, créer 2 fichiers et donner la position du point J pour l’application de la condition de Kutta-Joukovsky.

Ecrire un fichier nacawxyz.dxf décrivant le profil avec un format CAO DXF (décrit en 2.2). On utilisera la fonction writedxf permettant d’écrire un fichier dxf qui vous est fournie et documentée en section 2.3. Nous considérerons ici des valeurs de Z nulles (étude 2D). Afin de décrire une courbe fermée, le dernier point du fichier devra être strictement équivalent au premier, et le premier point doit correspondre au bord de fuite.
Ecrire un fichier nacawxyz.txt contenant sur 2 colonnes les coordonnées des points sur le profil en partant du bord de fuite, avec un parcours dans le sens trigonométrique. Afin de décrire une courbe fermée, le dernier point du fichier devra être strictement équivalent au premier.
définir les coordonnée $x_{j},x_{j}$ du point J où la condition de kutta-Jukowski doit être appliquée. Ce point doit se situer sur la tangente au bord de fuite. Il n’est numériquement pas possible d’appliquer cette condition exactement au bord de fuite. Le point J sera donc situé à une distance d du bord de fuite, tel que $d<corde/20$. On donnera aussi les coordonnées de la tangente au bord de fuite.
Le site http://airfoiltools.com/ propose un référencement des caractéristiques aérodynamique de profil cambré. En utilisant leur moteur de recherche http://airfoiltools.com/search/index, recherchez le profil étudier. Grâce à l’onglet Airfoil detail, vous pouvez obtenir les caractéristiques (coeff traînée et portance en fonction de l’angle) pour différentes valeur du nombre de Reynolds. En choisissant une valeur de Re pertinente , extraire les valeur de $C_{p}$ en fonction de l’angle $\alpha$, que vous reporterez dans un fichier nacawxyz_portance.dat.

3.3.4. Bibliothèque d’analyse#

On va construire une bibliothèque d’analyse de la solution écrite en Python. On définit une structure de données Naca, correspondant à un profil naca d’épaisseur e, de point de cambrure (xa,ya), pour stocker N points X,Y (vecteurs numpy) sur le Naca (en partant du bord de fuite), les valeurs de la pression Pr (matrice numpy (N,nv) ) et de la vitesse (U, V) en ces N points fonction des nv valeurs alpha (vecteur numpy) de l’angle $\alpha$ de la vitesse $U_{0}$ par rapport à l’horizontal.

    # donnees sur le NACA
    class Naca():
        def __init__(self,E,X0,Y0):
            # caracteristique du profil (epaisseur, cambrure)
            self.e  = E         
            self.xa = X0
            self.ya = Y0
            # points sur le profil
            self.X=None
            self.Y=None
            # donnees fct de l'angle d'incidence alpha
            self.alpha=None
            self.Pr=None
            self.U =None
            self.V =None
            return 

Écrire une fonction lecture(nom_fichier,E,Xa,Ya) qui lit les données de pression dans le fichier de nom nom_fichier pour un profil d’épaisseur E et de point de cambrure Xa,Ya, et qui renvoie une structure naca remplie. Pour l’extraction des valeurs de alpha stockées sur une ligne du fichier, on pourra utiliser la fonction read_alpha(line) fournie.
Écrire une fonction portance(naca) qui renvoie un vecteur numpy contenant la force de portance (par unité de longueur transverse) calculée avec les données de la structure naca pour toutes les valeurs de $\alpha$. On pourra utiliser la fonction trapz de la bibliothèque numpy pour calculer les intégrales. On pourra calculer la portance en projettant la résultante de pression suivant x et y.
Écrire une fonction lecture_vitesse(nom_fichier,naca) qui lit les données de vitesse dans le fichier de nom nom_fichier et qui initialise la vitesse dans la structure naca. La fonction ne renvoie rien directement, mais modifie la structure naca (vitesse U et V).
Écrire une fonction circulation(naca) qui renvoie un vecteur numpy contenant la circulation de vitesse calculée avec les données de la structure naca pour toutes les valeurs de $\alpha$. On pourra utiliser la fonction trapz de la bibliothèque numpy pour calculer les intégrales. On pourra calculer la circulation en décomposant $\overrightarrow{U}.\overrightarrow{dl}$ suivant x et y.

On fournit un canevas python bib_naca.py.

3.3.5. Analyse des données#

En utilisant la bibliothèque précédente, écrire un programme python pour analyser les deux fichiers de résultats (pression et vitesse). Pour le profil naca étudié, on calculera l’angle $\theta_{0}$ de la ligne moyenne au bord d’attaque et l’angle $\theta_{1}$ de la ligne moyenne au bord de fuite. On calculera la portance et la circulation de vitesse en fonction de l’angle d’incidence $\beta = \alpha-\theta_{0}$ par rapport à la ligne moyenne. Pour déterminer la précision des simulations, on calcule la vitesse au bord de fuite, et on définit l’erreur comme le rapport $err=|un/ut|$ entre la vitesse normale $un$ et la vitesse tangentielle $ut$ au bord de fuite. On justifiera ce choix.

Une caractéristique importante d’un profil est la position du centre de pousée. C’est le point d’application de la résultante aérodynamique; il est situé sur la corde du profil et se déplace en fonction de l’incidence (il se situe en général à 1/4 de la corde du bord d’attaque). Ce point d’application $P\left(x_{P},y_{P}\right)$ (aussi appelé centre de poussée) de la résultante des forces aérodynamiques $\overrightarrow{F}=\overrightarrow{F_{p}}+\overrightarrow{F_{t}}$ doit vérifier la relation suivante sur le moment de la résultante par rapport à un point $O$ quelconque:

\[\begin{split}\begin{aligned} \mathcal{M}_{O} & =\mathcal{M}_{P}+\overrightarrow{OP}\wedge\overrightarrow{F}\,\,\Rightarrow\overrightarrow{OP}\wedge\overrightarrow{F}=\int_{S}\overrightarrow{OM}\wedge\overrightarrow{dF}\mbox{ car }\mathcal{M}_{P}=0\\\end{aligned}\end{split}\]

avec $\overrightarrow{F}={\displaystyle \int_{S}-p\overrightarrow{n}dS}$, $\overrightarrow{dF}=-p\overrightarrow{n}dS$, $M\left(x_{M},y_{M}\right)$ un point du profil et $\overrightarrow{n}$ le vecteur normal sortant. En appliquant cette relation, on obtient l’équation de la droite qui porte la résultante, ce qui n’est pas suffisant pour définir la position du centre de poussée. Mais en calculant le centre de poussée sur la ligne moyenne on lève cette indétermination.

On demande de sortir les courbes au format pdf pour le compte rendu:

la courbe des polaires de pression sur le profil calculées avec comsol pour différentes valeurs de $\beta=\alpha-\theta_{0}$
la portance et la circulation en fonction de $\beta=\alpha-\theta_{0}$ calculées avec comsol
l’estimation de l’erreur en fonction de $\beta=\alpha-\theta_{0}$
Calcul de la position du centre de poussée en fonction de $\beta$. On discutera de la stabilité du profil en comparant cette position de $P$ par rapport au centre de gravité du profil.

3.3.6. Compte rendu#

Écrire un compte rendu de 3 pages maximum, pour déterminer la précision de la simulation comsol, donner l’évolution de la portance en fonction de l’angle d’incidence, ainsi que la loi de portance en fonction de la circulation.

3.3.7. Annexe#

3.3.7.1. exemple de fichier DXF pour un triangle ([0,0,0],[1,0,0],[0.5,1,0])#


    triangle

    SECTION  

    ENTITIES  

    LINE



0

0  

0

0  

0  

0  

    LINE



0
     20
0  

0

5  

0  

0  

    LINE



5

0  

0

0  

0  

0  

    ENDSEC
    0
    EOF

3.4. Ecoulement autour d’un profil NACA#

profil naca

%matplotlib inline
import numpy as np
import matplotlib.pyplot as plt

3.4.1. Définition d’un profil NACA#

Un profil NACA est un profil d’aile standardisé défini à partir de 3 paramètres:

l’épaisseur relative $t$
la cambrure relative $y_{a}$ au point $x_{a}$ de la ligne moyenne

3.4.1.1. profil Naca symétrique#

Dans ce cas $y_{a}=x_{a}=0$ et le seul paramètre est l’épaisseur relative $t$. L’extrados est donné par la courbe $y_{e}==F(x)$ suivante et l’intrados par la courbe symétrique $y_{i}=-F(x)$

\[F(x)=5.0\,t\left(0.2969\,\sqrt{x}-0.1260\,x-0.3516\,x^{2}+0.2843\,x^{3}-0.1036\,x^{4}\right)\,\,\,\mbox{ pour }\,\,0\leq x\leq L_{c}\]

La longueur de la corde $L_{c}$ suivant x est donnée par la racine de $F(L_{c})=0$, i.e. $L_{c}=1$ (sans dimension)

Un profil NACA symétrique à 4 chiffres est donnée par la nomenclature NACA00yz où les 2 chiffres yz représente l’épaisseur relative du profil $t=yz/100$.

3.4.1.2. profil cambré Naca à 4 chiffres#

Dans ce cas on donne la ligne moyenne du profil en fonction de $x_{a}$ et $y_{a}$. Cette ligne moyenne a pour équation en fonction de x

\[\begin{split} y_{m}=\left\{ \begin{array}{l} \frac{y_{a}(2x_{a}-x)\,x}{x_{a}^{2}}\,\,\,\mbox{ si }\,\,\,\,0\leq x\leq x_{a}\\ \frac{y_{a}(1-x)(1+x-2x_{a})}{(1-x_{a})^{2}}\,\,\,\mbox{ si }\,\,\,\,x_{a}\leq x\leq L_{c} \end{array}\right. \end{split}\]

L’équation de l’extrados ${y_{e},x_{e}}$ est obtenu à partir de cette ligne moyenne, en introduisant l’angle $\theta$

\[\theta(x)=\frac{dy_{m}}{dx}\]

et la fonction F(x)précédente:

\[x_{e}(x) = x-F(x)\,\sin(\theta(x))\]

\[y_{e}(x) = y_{m}+F(x)\,\cos(\theta(x))\,\,\,\,\mbox{ pour }0\leq x\leq L_{c}\]

L’équation de l’intrados ${x_{i},y_{i}}$ correspond à l’équation symétrique en y

\[x_{i}(x) = x+F(x)\,\sin(\theta(x))\]

\[y_{i}(x) = y_{m}-F(x)\,\cos(\theta(x))\,\,\,\,\mbox{ pour }0\leq x\leq L_{c}\]

La longueur de la corde $L_{c}$ vérifie $y_{e}(L_{c})=y_{i}(L_{c})$

Un profil NACA cambré à 4 chiffres est donnée par la nomenclature NACAwxyz où les 2 chiffres yz représentent l’épaisseur relative du profil ($t=yz/100$), le chiffre w la courbure relative en % ($y_{a}=w/100$) et le chiffre x la position de la flèche de la ligne moyenne en 1/10 de corde ($x_{a}=x/10$).

3.4.1.3. profil cambré Naca à 5 chiffres#

Un profil NACA cambré à 5 chiffres est donnée par la nomenclature NACAabcyz où les 2 derniers chiffres yz représentent l’épaisseur relative du profil (t=yz/100), et abc vaut 210,220,230,240,250 pour des cambrure simples et 211,221,231,241,251 pour des cambrures doubles. Les équations sont les mêmes sauf pour l’équation de la ligne moyenne y_{m} que l’on la trouvera sur le site Wikipédia ici :

WIKIPEDIA

3.4.2. liste des profils étudiés#

fichier='listenaca.txt'
f = open(fichier,'r')
profils = f.readlines()
f.close()
# remove /n
profils = [p.strip() for p in profils]
print(profils)

3.4.3. Tracé de profils NACA#

def traceProfil(Ep,Xa,Ya):
    def F(x):
        return 5*Ep*(0.2969*np.sqrt(x)-0.1260*x-0.3537*x**2+0.2843*x**3-0.1015*x**4)
    def Ym(x):
        if (x<Xa): 
            return Ya*(2*Xa-x)*x/Xa**2
        else :
            return Ya*(1-x)*(1+x-2*Xa)/(1-Xa)**2
    def Theta(x):
        if (x<Xa): 
            return Ya*(2*Xa-x)/Xa**2-Ya*x/Xa**2
        else :
            return Ya*(1-x)/(1-Xa)**2-Ya*(1+x-2*Xa)/(1-Xa)**2
    #
    X = np.linspace(0,1,100)
    FX= np.zeros(X.size); THETA= np.zeros(X.size); YM= np.zeros(X.size);
    for i,x in enumerate(X):
        FX[i]=F(x);
        THETA[i]=Theta(x);
        YM[i]=Ym(x);
    # nbre de points
    n=len(X);
    XX=np.zeros((2*n-1));
    YY=np.zeros((2*n-1));
    ZZ=np.zeros((2*n-1));
    XX[:n]=X  - FX*np.sin(THETA);
    YY[:n]=YM + FX*np.cos(THETA);
    XX[n::]= X[n-2::-1]  + FX[n-2::-1]*np.sin(THETA[n-2::-1])
    YY[n::]= YM[n-2::-1] - FX[n-2::-1]*np.cos(THETA[n-2::-1])
    # tracer profil et ligne moyenne
    plt.plot(XX,YY,'-')
    plt.plot(X,YM,'--')
    return

3.4.3.1. Profil NACA wxyz#

Un profil NACA cambré à 4 chiffres est donnée par la nomenclature NACAwxyz où les 2 chiffres yz représentent l’épaisseur relative du profil ($t=yz/100$), le chiffre w la courbure relative en % $y_{a}=w/100$ et le chiffre x la position de la flèche de la ligne moyenne en 1/10 de corde $x_{a}=x/10$.

$NACAy_ax_ae_p$

# profil Naca4412 (xyee)
plt.figure(figsize=(12,6))
traceProfil(12./100,4./10,4./100.)
plt.title("Naca 4412")
plt.axis('equal');

# profil Naca2412 (xyee)
plt.figure(figsize=(12,6))
traceProfil(12./100,4./10,2./100.)
plt.title("Naca 2412")
plt.axis('equal');

3.4.4. Simulation avec XFOIL naca0012#

from xfoil import XFoil
xf = XFoil()

from xfoil.test import naca0012
xf.airfoil = naca0012

# naca 0012 xfoil
X = xf.airfoil.coords[:,0]
Y = xf.airfoil.coords[:,1]
# data COMSOL
X1,Y1 = np.loadtxt('Naca0012.txt',unpack=True)

plt.plot(X,Y,label="Comsol")
plt.plot(X1,Y1,'--',label="XFoil")
plt.title("Naca 0012")
plt.legend()
plt.axis('equal');

3.4.5. Lecture des résultats de simulation COMSOL#

from sol_naca import Naca,lecture,lecture_vitesse, portance
naca12 = lecture("PressionNaca0012.dat",12,0,0)
Fp = portance(naca12)
for num in range(0,naca12.alpha.size,2):
    plt.plot(naca12.X,naca12.Pr[:,num],label="$\\alpha$={:.1f}".format(np.degrees(naca12.alpha[num])))
plt.legend()
plt.xlabel("$x/L_c$")
plt.ylabel("$(p-p_0)/E_0$")
plt.title("Répartition de pression (calcul COMSOL)");

# solution pour un angle donnee
num=3
alpha = naca12.alpha[num]
Alpha = np.degrees(alpha)
# Xfoil
Cl, Cd, Cm, Cp = xf.a(Alpha)
# polaire: pression adim
Xp,Pr = xf.get_cp_distribution()
Pr = Pr*0.5
print("Cp Portance à alpha={:.1f}° COMSOL={:.3f} Xfoil={:.3f}".format(Alpha,Fp[num]/0.5,Cl))

plt.plot(naca12.X,naca12.Pr[:,num],label="COMSOL $\\alpha$={:.1f}".format(Alpha))
plt.plot(Xp,Pr,label="xfoil")
plt.title("Repartition pression")
plt.legend();

Cl,Cd,Cd,Cm,Cp,alpha

xref=0.25; yref=0.0
CL =  (np.trapz(naca12.Pr[:,num],naca12.X)*np.cos(alpha)   + np.trapz(naca12.Pr[:,num],naca12.Y)*np.sin(alpha))/0.5
CM = -(np.trapz(naca12.Pr[:,num]*(naca12.X-xref),naca12.X) + np.trapz(naca12.Pr[:,num]*(naca12.Y-yref),naca12.Y))/0.5
print(CM," ecart centre de poussée / 1/4 corde =",CM/CL)

Par définition le coefficiant de moment = moment des forces de pression / à P = xref=0.25,yref=0.0 (i.e. 1/4 cordre)

\[ Cm = \int p\vec{n} \wedge \vec{MP} ds = -\int p \left((x-x_{ref}) dx + (y-y_{ref}) dy\right)\]

Donc le centre de poussée = 1/4 cordre / bord d’attaque est une bonne approximation (erreur 1/100)

3.4.6. FIN#

1. Références#

1.1. Bibliographie#

[BUF18]

Marc BUFFAT. Inpros: introduction à la programmation scientifique. 2018. URL: https://perso.univ-lyon1.fr/marc.buffat/2022/BOOK_INPROS/index.html.

[Dow16]

Allen B. Downey. How to think like a computer scientist. 2016. URL: https://www.greenteapress.com/wp/think-python-2e/.

[Fou01]

Python Software Foundation. Python official tutorials. 2001. URL: https://docs.python.org/.

[Fou11]

Python Software Foundation. Scientific python. 2011. URL: https://scipy.org/.

[Fou21]

Python Software Foundation. Symbolic python. 2021. URL: https://sympy.org/.

[fou23]

Python Software foundation. Jupyter notebook. 2023. URL: https://jupyter.org/.

[Mar18]

Nicolas Martin. La méthode scientifique: le meilleur des mondes informatiques ( gérard berry ). 2018. URL: https://www.radiofrance.fr/franceculture/podcasts/la-methode-scientifique/gerard-berry-le-meilleur-des-mondes-informatiques-4408389.

(C) Marc BUFFAT, 2024

Cours Ateliers Numériques en Mécanique

Contenu

Cours Ateliers Numériques en Mécanique#

1. Introduction#

1.1. Objectif du cours#

1.1.1. Outils informatiques:#

1.1.2. Applications#

2. Histoire et enjeux de l’informatique scientifique#

2.1. Historique de l’informatique#

2.1.1. Premiers Algorithmes#

2.1.2. Origine du nom Algorithme#

2.1.3. Prémices de l’informatique#

2.1.4. Les pionniers: machine de Turing#

2.1.5. Premiers ordinateurs électronique#

2.1.6. Naissance de la micro-informatique ~ 1970#

2.1.7. L’ère internet (à partir 1990)#

2.1.8. L’ère actuelle: mobilité , données#

2.1.9. Conclusion#

2.2. Utilité d’un ordinateur en science#

2.2.1. Modélisation numérique#

2.2.2. Méthode#

2.3. Utilisation de l” informatique en science#

2.3.1. modélisation en mécanique de système complexe#

2.3.2. expérience numérique fondamentale#

2.3.3. CFD = Computational Fluid Dynamics#

2.3.4. CFD = Color Fluid Dynamics ?#

2.3.5. Attention au mirage de la simulation numérique#

2.4. FIN de la leçon#

3. Python pour les scientifiques#

3.1. Introduction au calcul scientifique avec Python#

3.1.1. Objectifs#

3.1.1.1. Apprendre:#

3.1.1.2. À la fin du cours, vous pourrez:#

3.1.2. Fonctionnement d’un ordinateur#

3.1.2.1. composants d’un ordinateur#

3.1.2.2. programmes / processus#

3.1.2.3. Système de fichiers (Unix)#

3.1.3. Execution d’un programme#

3.1.3.1. Calcul de la somme d’une série:#

3.1.3.2. Algorithme de base#

3.1.3.3. Execution sur un ordinateur#

3.1.3.4. Langage machine (binaire)#

3.1.4. Langages de programmation#

3.1.4.1. Langage assembleur#

3.1.4.2. Langage compilé: compilateur#

3.1.4.3. Langage interprété:#

3.1.4.3.1. Matlab (Matrix Laboratory)#

3.1.4.3.2. langage python#

3.1.5. Historique des langages de programmation#

3.1.5.1. Langages les plus utilisées par les scientifiques#

3.1.6. Exigences du calcul scientifique#

3.1.6.1. Démarche du calcul scientifique#

3.1.7. Documentation#

3.1.7.1. Installation (logiciel libre)#

3.1.8. Fin de la leçon#

3.2. Base de programmation en Python#

3.2.1. Questions préliminaires#

3.2.2. Historique#

3.2.3. Pourquoi Python#

3.2.3.1. Avantages du langage Python#

3.2.3.2. Défauts du langage Python#

3.2.3.3. Approche générale de développement#

3.2.4. Philosophie de Python (Zen of python)#

3.2.5. Concepts de base#

3.2.5.1. Algorithme#

3.2.5.2. Programme informatique#

3.2.6. Les chaînes de caractères en Python#

3.2.6.1. exercice#

3.2.6.2. Indexation#

3.2.6.3. sous-chaînes#

3.2.6.4. manipulation des chaînes de caractères#

3.2.6.4.1. commentaires#

3.2.6.4.2. application#

3.2.6.4.3. count()#

3.2.6.4.4. find() & index()#

3.2.6.4.5. index()#

3.2.6.4.6. exercices:#

3.2.6.4.7. split()#

3.2.7. Les listes en Python#

3.2.7.1. exercice:#

3.2.6.4. manipulation des chaînes de caractères#

3.2.6.4.3. `count()`#

3.2.6.4.4. `find()` & `index()`#

3.2.6.4.5. `index()`#

3.2.6.4.7. `split()`#

3.2.7.2. `append()`: ajout d’éléments#

3.2.7.3. `in` : test si dans la liste#

3.2.8.1. Itération avec `for`#

3.2.8.4. Test avec `if`#