"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"############ Si vous faites ce TP pour la première fois, il faut télécharger le fichier Trajectoire.dat ###########\n",
" #################### A ne faire que la première fois ######################\n",
"# Décommentez la ligne suivante et executer la cellule\n",
"#! wget hpp.education/Lessons/CFD_part2/Files/Trajectoire.dat"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"import numpy as np\n",
"import matplotlib.pyplot as plt\n",
"import scipy\n",
"import time\n",
"#Option pour afficher les figures dans le notebook et eviter le plt.show():\n",
"%matplotlib inline "
]
},
{
"cell_type": "markdown",
"metadata": {
"collapsed": true
},
"source": [
"On se propose dans ce TP de calculer numériquement la trajectoire d'un projectile sphérique dans l'air lancé à une vitesse initiale $U_0$ dans la direction $\\theta_0$ et soumis à son poid $\\mathbf{P}=m\\mathbf{g}$ et aux forces aérodynamique de trainée et de portance $\\mathbf{F_D}$ et $\\mathbf{F_L}$.\n",
"\n",
"\n",
"# I - Formalisation du problème\n",
"\n",
"## Mise en équation:\n",
"\n",
"Le principe fondamentale de la dynamique appliqué à ce projectil donne:\n",
"\n",
"$$\n",
"m\\dfrac{d\\mathbf{U}}{dt}=\\mathbf{P}+\\mathbf{F_D}+\\mathbf{F_L}\n",
"$$\n",
"\n",
"Avec $\\mathbf{U}=\\left(\\begin{array}{c} u_x \\\\ u_y \\end{array}\\right)$ et $U=\\lVert{\\mathbf{U}}\\lVert$.\n",
"\n",
"En projetant cette équation sur les deux axes du repère $(Oxy)$ on obtient:\n",
"\n",
"$$\n",
"m\\dfrac{du_x}{dt}=-F_D\\cos(\\theta)-F_L\\sin(\\theta)\n",
"$$\n",
"\n",
"$$\n",
"m\\dfrac{du_y}{dt}=-F_D\\sin(\\theta)+F_L\\cos(\\theta)-mg\n",
"$$\n",
"\n",
"\n",
"Avec:\n",
"$$F_D=\\dfrac{1}{2}\\rho U^2 S C_d$$\n",
"$$F_L=\\dfrac{1}{2}\\rho U^2 S C_l$$\n",
"$$\\theta=arctan\\left(\\dfrac{u_y}{u_x}\\right)$$\n",
"\n",
"\n",
"Le système ainsi présenté est donc parfaitement décrit à l'aide de son vecteur position / vitesse défini par:\n",
"\n",
"$$\n",
"\\mathbf{W}=\\left(\\begin{array}{c} x \\\\ y \\\\ u_x \\\\ u_y \\end{array}\\right)\n",
"$$\n",
"\n",
"Le problème du calcul de trajectoire (mettre à jour $\\mathbf{W}$ à chaque instant) se ramène donc à un problème du type:\n",
"$$\n",
"\\dfrac{d\\mathbf{W}}{dt}=K(\\mathbf{W},t)\n",
"$$\n",
"\n",
"avec \n",
"\n",
"$$\n",
"K(\\mathbf{U},t)=\\left(\\begin{array}{c} u_x \\\\ u_y \\\\ \\dfrac{-F_D\\cos(\\theta)-F_L\\sin(\\theta)}{m} \\\\ \\dfrac{-F_D\\sin(\\theta)+F_L\\cos(\\theta)-mg}{m}\\end{array}\\right)\n",
"$$"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"On défini les paramètres suivants qui serviront de référence pendant le TP:"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"m=1. # Masse de 1kg\n",
"D=0.1 # Diametre de la sphere de 10cm\n",
"g=9.81 # Acceleration de la pesanteur\n",
"S=np.pi*D**2/4. # Surface de référence\n",
"rho=1.24 # Masse volumique standard de l'air\n",
"Cd=0.4 # Coefficient de trainée\n",
"Cl=0.01 # Coefficient de portance"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Initialisation\n",
"\n",
"Avant de commencer le calcul, il faut initialiser le problème c'est à dire définir l'état initial du vecteur position/vitesse $\\mathbf{W}$. On prendra comme position initiale le point $(0,1)$."
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"# Etat initial:\n",
"U0=60. # Norme de vitesse initiale\n",
"theta0=np.pi/3. # Orientation du vecteur vitesse au temps 0\n",
"W=np.array([0.,0.,0.,0.]) \n",
"W[0]=?? # position x initiale\n",
"W[1]=?? # position y initiale\n",
"W[2]=?? # vitesse ux initiale\n",
"W[3]=?? # vitesse uy initiale"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Définir la fonction K qui va permettre la mise à jour de $\\mathbf{W}$ comme définie plus haut:"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"def K(W,t):\n",
" theta=...\n",
" Fd=...\n",
" Fl=...\n",
" K1=...\n",
" K2=...\n",
" K3=...\n",
" K4=...\n",
" return np.array([K1,K2,K3,K4])"
]
},
{
"cell_type": "markdown",
"metadata": {
"collapsed": true
},
"source": [
"
II - Résolution numérique d'ordre 1: Schéma d'Euler
\n",
"\n"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Le schéma d'euler explicite s'écrit de la façon suivante (Cf chap4):\n",
"\n",
"$$\n",
"\\mathbf{W}^{n+1}=\\mathbf{W}^{n}+\\Delta t K(\\mathbf{W}^{n},t)\n",
"$$\n",
"\n",
"Implémenter le schéma d'Euler pour le calcul de la vitesse à chaque instant, puis en déduire le calcul de la position à chaque instant. \n",
"\n",
"Comparer vos résultats avec les valeurs expérimentales fourni dans le fichier Trajectoire.dat."
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"# Chargement des résultats théoriques:\n",
"th=np.loadtxt('Trajectoire.dat',dtype='f')\n",
"# !!! Attention, il faut télécharger ce fichier avec la toute première cellule"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"dt=0.2\n",
"t=0\n",
"while W[1]>0:\n",
" t+=dt\n",
" W=....\n",
" plt.plot(...,...,'.r')\n",
"xmax=... # Valeur de x lors de l'impact au sol\n",
"plt.title(str(t))\n",
"plt.xlabel('xmax: '+str(xmax)+'m')\n",
"plt.plot(th[:,0],th[:,1],'k')"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Discussion:\n"
]
},
{
"cell_type": "markdown",
"metadata": {
"collapsed": true
},
"source": [
"
III - Résolution numérique d'ordre élevé: Schéma de Runge-Kutta
\n",
"\n",
"## RK2\n",
"\n",
"Le schéma de Runge_Kutta permet d'augmenter l'ordre d'intégration et consiste à calculer des étapes intermédiaires entre 2 itérations. (On parle de sous-iterations). \n",
"\n",
"\n",
"Implémenter le schéma de Runge-Kutta d'ordre 2 vu en cours pour la mise à jour de $\\mathbf{W}$. On veillera à bien réinitialiser les grandeurs initiales:"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"# Etat initial:\n",
"U0=60. # Norme de vitesse initiale\n",
"theta0=np.pi/3. # Orientation du vecteur vitesse au temps 0\n",
"W=np.array([0.,0.,0.,0.]) \n",
"W[0]=?? # position x initiale\n",
"W[1]=?? # position y initiale\n",
"W[2]=?? # vitesse ux initiale\n",
"W[3]=?? # vitesse uy initiale"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"dt=0.2\n",
"t=0\n",
"while y>0:\n",
" t+=dt\n",
" k1=...\n",
" k2=...\n",
" W=...\n",
" plt.plot(...,'.r')\n",
"xmax=...\n",
"plt.title(str(t))\n",
"plt.xlabel('xmax: '+str(xmax)+'m')\n",
"plt.plot(th[:,0],th[:,1],'k')"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Discussions:\n",
"\n",
"\n",
"## RK4\n",
"\n",
"Implémenter l'ordre 4 du schéma de Runge-Kutta et commenter les écarts obeservés avec le schéma d'Euler et la version d'ordre 2."
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"# Etat initial:\n",
"U0=60. # Norme de vitesse initiale\n",
"theta0=np.pi/3. # Orientation du vecteur vitesse au temps 0\n",
"W=np.array([0.,0.,0.,0.]) \n",
"W[0]=?? # position x initiale\n",
"W[1]=?? # position y initiale\n",
"W[2]=?? # vitesse ux initiale\n",
"W[3]=?? # vitesse uy initiale"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"dt=0.2\n",
"t=1.\n",
"while y>0:\n",
" t+=dt\n",
" k1=...\n",
" k2=...\n",
" k3=...\n",
" k4=...\n",
" W=...\n",
" plt.plot(x,y,'.r')\n",
"xmax=...\n",
"plt.title(str(t))\n",
"plt.xlabel('xmax: '+str(xmax)+'m')\n",
"plt.plot(th[:,0],th[:,1],'k')"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Discussions..."
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"# Conclusion\n",
"\n",
"Présentez une synthèse de vos résultats et des étapes de calcul de votre TP."
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"from IPython.core.display import HTML\n",
"style=open('notebooks.css', \"r\").read()\n",
"HTML(style)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"\n",
"\n",
"