{ "cells": [ { "cell_type": "markdown", "id": "3190d34c", "metadata": {}, "source": [ "## La bille dans l'étonnoir" ] }, { "cell_type": "markdown", "id": "051190d2", "metadata": {}, "source": [ "Une bille, assimilée à un point matériel, se déplace sur la surface d'un cône renversé sous l'influence de la gravité. Tout frottement est négligeable. La vitesse et la position initiale à $t=0$ sont données; on souhaite déterminer la trajectoire. \n", "\n", "La masse de la bille est $m$ et l'angle d'ouverture du cône est $2\\alpha$. \n", "\n", "\n", "\n", "Un problème avec des contraintes comme celui-ci se traite le plus facilement dans le formalisme lagrangien ou hamiltonien de la mécanique classique. Mais pour ceux et celles qui ne connaissent pas bien ces approches, on utilisera ici le formalisme newtonien. \n", "\n", "Tout d'abord, choisissons les coordonnées: en vue de la symétrie du problème, on utilisera des coordonnées cylindriques $(r, \\phi, z)$, avec l'axe des $z$ = l'axe de symétrie verticale, $r=\\sqrt{x^2+y^2}$ et $\\tan\\phi=y/x$. Dans le plan des $(x,y)$, vu d'en haut, on a alors les coordonnées et vecteurs de base suivants:\n", "\n", "\n", "\n", "Les équations de mouvement sont\n", "$$\\vec F=m\\ddot{\\vec x}$$\n", "avec $\\vec F$ la somme des forces et\n", "$$\n", "\\vec x=x\\;\\vec e_x+y\\;\\vec e_y+z\\;\\vec e_z=r\\;\\vec e_r+z\\;\\vec e_z\n", "$$\n", "la position de la bille. On va écrire l'accélération $\\ddot{\\vec x}$ en coordonnées cylindriques, sachant que les vecteurs de base, pour une base mobile, sont des fonctions des coordonnées qui évoluent avec le temps. La vitesse et l'accélération peuvent s'écrire\n", "$$\\dot{\\vec x}=\\dot r\\;\\vec e_r+r\\;\\dot{\\vec e_r}+\\dot z\\;\\vec e_z$$\n", "$$\\ddot{\\vec x}=\\ddot r\\;\\vec e_r+2\\dot r\\;\\dot{\\vec e_r}+r\\;\\ddot{\\vec e_r}+\\ddot z\\;\\vec e_z$$\n", "Avec $\\vec e_r=\\cos(\\phi)\\;\\vec e_x+\\sin(\\phi)\\;\\vec e_y$ et $\\vec e_\\phi=-\\sin(\\phi)\\;\\vec e_x+\\cos(\\phi)\\;\\vec e_y$, on a\n", "$$\\dot{\\vec e_r}=-\\dot\\phi\\sin(\\phi)\\;\\vec e_x+\\dot\\phi\\cos(\\phi)\\;\\vec e_y=\\dot\\phi\\;\\vec e_\\phi$$\n", "$$\\ddot{\\vec e_r}=-\\dot\\phi^2\\cos(\\phi)\\;\\vec e_x-\\ddot\\phi\\sin(\\phi)\\;\\vec e_x-\\dot\\phi^2\\sin(\\phi)\\;\\vec e_y+\\ddot\\phi\\cos(\\phi)\\;\\vec e_y=-\\dot\\phi^2\\;\\vec e_r+\\ddot\\phi\\;\\vec e_\\phi$$\n", "Tout ensemble, on peut donc exprimer l'accélération en coordonnées cylindriques comme\n", "$$\n", "\\ddot{\\vec x}=(\\ddot r-r\\dot\\phi^2)\\;\\vec e_r+(r\\ddot\\phi+2\\dot r\\dot\\phi)\\;\\vec e_\\phi+\\ddot z\\;\\vec e_z\\,.\n", "$$\n", "On définit le moment cinétique en direction des $z$,\n", "$$\n", "\\ell=mr^2\\dot\\phi\\,,\n", "$$\n", "ce qui permet d'écrire\n", "$$\n", "\\ddot{\\vec x}=(\\ddot r-r\\dot\\phi^2)\\;\\vec e_r+\\frac{1}{mr}\\dot\\ell\\;\\vec e_\\phi+\\ddot z\\;\\vec e_z\\,.\\qquad\\qquad(\\star)\n", "$$\n", "\n", "Ayant réécrit l'accélération en coordonnées cylindriques, il reste à établir un bilan des forces pour obtenir $\\vec F$ dans les équations de mouvement. Dans un plan qui contient l'axe des $z$, on a la force normale $\\vec F_n$ et la force gravitationnelle $\\vec F_g$ selon le schéma\n", "\n", "\n", "\n", "La force gravitationnelle est $\\vec F_g=-mg\\;\\vec e_z$. La force normale est orthogonale à la surface du cône et prend donc la forme $\\vec F_n = F_n\\;(\\sin(\\alpha)\\;\\vec e_z-\\cos(\\alpha)\\;\\vec e_r)$; on va laisser sa magnitude $F_n$ indéterminée pour l'instant. La somme des forces est alors\n", "$$\n", "\\vec F=\\vec F_n+\\vec F_g=-F_n\\cos(\\alpha)\\;\\vec e_r+(F_n\\sin(\\alpha)-mg)\\;\\vec e_z\\,.\\qquad\\qquad(\\star\\star)\n", "$$\n", "Avec le principe fondamental de la dynamique $\\vec F=m\\ddot{\\vec x}$ et éqs.$(\\star)$ et $(\\star\\star)$, on obtient finalement trois équations de mouvement pour les trois composantes:\n", "$$\\ddot r=r\\dot\\phi^2-\\frac{F_n}{m}\\cos(\\alpha)\\,,$$\n", "$$\\dot\\ell = 0\\,,$$\n", "$$\\ddot z=\\frac{F_n}{m}\\sin(\\alpha)-g\\,.$$\n", "* La deuxième équation $\\dot\\ell=0$ indique que le moment cinétique est conservé, ce qui est une conséquence de l'invariance par rotations autour de l'axe des $z$: $\\ell = mr^2\\dot\\phi = $ const.\n", "* Les coordonnées $r$ et $z$ ne sont pas indépendantes, car la bille doit rester sur la surface du cône qui est donnée par $r=z\\tan(\\alpha)$. Si on injecte $\\ddot r=\\ddot z\\tan(\\alpha)$ dans l'équation de mouvement, cela détermine la magnitude de la force normale $F_n$ car\n", "$$\n", "-\\frac{F_n}{m}\\cos(\\alpha)+r\\dot\\phi^2=\\ddot r=\\ddot z\\tan(\\alpha)=\\left(\\frac{F_n}{m}\\sin(\\alpha)-g\\right)\\tan(\\alpha)\n", "$$\n", "d'où\n", "$$\\frac{F_n}{m}\\left(\\sin(\\alpha)\\tan(\\alpha)+\\cos(\\alpha)\\right)=r\\dot\\phi^2+g\\tan(\\alpha)$$\n", "ou bien\n", "$$ F_n=m\\left(r\\dot\\phi^2\\cos(\\alpha)+g\\sin(\\alpha)\\right)\\,. $$\n", "* Si on injecte cette expression encore dans l'équation de mouvement pour $r$, on obtient enfin\n", "$$\\ddot r=r\\dot\\phi^2-\\frac{F_n}{m}\\cos(\\alpha)=r\\dot\\phi^2\\sin^2(\\alpha)-g\\sin(\\alpha)\\cos(\\alpha)$$\n", "ou bien, avec $\\dot\\phi=\\frac{\\ell}{mr^2}$,\n", "$$\n", "\\boxed{\\ddot r=\\frac{\\ell^2}{m^2 r^3}\\sin^2(\\alpha)-g\\sin(\\alpha)\\cos(\\alpha)\\,.}\n", "$$\n", "\n", "C'est cette équation qu'il faudra résoudre numériquement. Pour ce faire, il convient d'écrire l'équation de second ordre pour $r$ comme deux équations de premier ordre en introduisant la vitesse radiale $v_r=\\dot r$. Ensemble avec l'équation pour $\\phi$, on a alors les trois équations différentielles\n", "$$\\boxed{\\dot r=v_r\\,,\\qquad \\dot v_r=\\frac{\\ell^2}{m^2r^3}\\sin^2(\\alpha)-g\\sin(\\alpha)\\cos(\\alpha)\\,,\\qquad \\dot\\phi=\\frac{\\ell}{mr^2}\\,.}\n", "$$\n", "\n", "Après tout ce travail préliminaire, on peut enfin procéder à la solution numérique. On prend une bille de $m=1$ g (même si, en fait, la solution finale sera indépendante de la masse), un angle $\\alpha=45^\\circ$, un point de départ à $r=10$ cm et $\\phi=0$, une vitesse radiale initiale de $0.2$ m/s et une vitesse angulaire initiale de $\\omega_0=\\dot\\phi(0)=4/$s. On calcule la solution pendant $5$ s avec un incrément de $0.1$ ms." ] }, { "cell_type": "code", "execution_count": null, "id": "aa03beae", "metadata": {}, "outputs": [], "source": [ "import numpy as np\n", "\n", "# Constantes physiques:\n", "g = 9.81 # accélération gravitationnelle en m/s^2\n", "alpha = np.pi/4 # 1/2 * angle d'ouverture\n", "m = 1.E-3 # masse de la bille en kg\n", "# Constantes numériques: \n", "h = 1.E-4 # pas d'incrément en s\n", "N = 50000 # nombre de pas à calculer\n", "# Conditions initiales: \n", "r0 = .1 # rayon\n", "v0 = .2 # vitesse radiale\n", "phi0 = 0.0 # angle\n", "omega0 = 4. # vitesse angulaire\n", "# Constantes dépendantes:\n", "L = m * r0**2 * omega0 # moment cinétique\n", "sa, ca = np.sin(alpha), np.cos(alpha)" ] }, { "cell_type": "markdown", "id": "71addc98", "metadata": {}, "source": [ "Afin de résoudre les équations de mouvement, on utilise d'abord la méthode d'Euler:" ] }, { "cell_type": "code", "execution_count": null, "id": "e061e855", "metadata": {}, "outputs": [], "source": [ "def euler(f, u0, h, N):\n", " # Solution de l'EDO u'(t) = f(u(t)) par la méthode d'Euler explicite\n", " # f = une fonction vectorielle avec 3 composantes\n", " # u0 = un tableau 1-dimensionnel contenant 3 conditions initiales\n", " # h = incrément, N = nombre de pas\n", " # Renvoie un tableau de (N+1)*3 valeurs de u entre t=0 et t=Nh \n", " u = np.empty([N+1, 3]) # variables dynamiques\n", " u[0] = u0 # u[0] = première ligne du tableau u = conditions initiales\n", " for n in range(N):\n", " u[n+1] = u[n] + h * f(u[n])\n", " return u" ] }, { "cell_type": "markdown", "id": "a9e9bdbf", "metadata": {}, "source": [ "Il faut encore spécifier la fonction f des membres de droite du système d'EDO." ] }, { "cell_type": "code", "execution_count": null, "id": "417ca698", "metadata": {}, "outputs": [], "source": [ "# La fonction vectorielle des membres de droite des e.d.m.\n", "# u = un vecteur à 3 composantes (r, v et phi au temps t)\n", "def f(u):\n", " r, v, phi = u\n", " fr = v\n", " fv = L**2 * sa**2 / (m**2 * r**3) - g * sa * ca\n", " fphi = L / (m * r**2)\n", " return np.array([fr, fv, fphi])" ] }, { "cell_type": "markdown", "id": "f008bd50", "metadata": {}, "source": [ "Calcul de la solution:" ] }, { "cell_type": "code", "execution_count": null, "id": "ba8fb199", "metadata": {}, "outputs": [], "source": [ "u0 = np.array([r0, v0, phi0])\n", "solution_euler = euler(f, u0, h, N)" ] }, { "cell_type": "markdown", "id": "5ad029fa", "metadata": {}, "source": [ "On trace la solution pour $r(t)$:" ] }, { "cell_type": "code", "execution_count": null, "id": "b45c0940", "metadata": {}, "outputs": [], "source": [ "import matplotlib.pyplot as plt\n", "\n", "t_euler = np.linspace(0, (N+1)*h, N+1)\n", "\n", "plt.plot(t_euler, solution_euler[:, 0])\n", "plt.xlabel(\"t [s]\")\n", "plt.ylabel(\"r(t) [m]\")\n", "plt.title(\"Coordonnée radiale en fonction du temps (Euler)\")\n", "plt.show()" ] }, { "cell_type": "markdown", "id": "74f7e182", "metadata": {}, "source": [ "La croissance de l'amplitude visible indique que la précision de la méthode n'est pas suffisante. On a le choix soit de calculer plus de points intermédiaires, soit d'utiliser une méthode plus précise, comme par exemple la méthode RK4:" ] }, { "cell_type": "code", "execution_count": null, "id": "b3a969bf", "metadata": {}, "outputs": [], "source": [ "def rk4(f, u0, h, N):\n", " # Solution de l'EDO u'(t) = f(u(t)) par la méthode de Runge-Kutta du 4ème ordre\n", " # Toutes les variables sont les mêmes que pour la fonction euler() ci-dessus\n", " u = np.empty([N+1, 3]) # variables dynamiques\n", " u[0] = u0 # conditions initiales\n", " for n in range(N):\n", " k1 = f(u[n])\n", " k2 = f(u[n] + h/2 * k1)\n", " k3 = f(u[n] + h/2 * k2)\n", " k4 = f(u[n] + h * k3)\n", " u[n+1] = u[n] + h * (k1/6 + k2/3 + k3/3 + k4/6)\n", " return u\n", "\n", "# La méthode RK4 permet de prendre 100-fois moins de points intermédiaires et toujours avoir une précision supérieure\n", "h = 1.E-2\n", "N = 500\n", "t_rk4 = np.linspace(0, (N+1)*h, N+1)\n", "\n", "solution_rk4 = rk4(f, u0, h, N)\n", "\n", "plt.plot(t_rk4, solution_rk4[:, 0])\n", "plt.xlabel(\"t [s]\")\n", "plt.ylabel(\"r(t) [m]\")\n", "plt.title(\"Coordonnée radiale en fonction du temps (RK4)\")\n", "plt.show()" ] }, { "cell_type": "markdown", "id": "89c9dcaa", "metadata": {}, "source": [ "On trace le rapport des deux (un point sur 100 pour la méthode d'Euler):" ] }, { "cell_type": "code", "execution_count": null, "id": "311b0527", "metadata": {}, "outputs": [], "source": [ "plt.plot(t_rk4, 1 - solution_euler[::100, 0] / solution_rk4[:, 0])\n", "plt.xlabel(\"t [s]\")\n", "plt.ylabel(\"Δr / r\")\n", "plt.title(\"Écart relatif\")\n", "plt.show()" ] }, { "cell_type": "markdown", "id": "343ba3a0", "metadata": {}, "source": [ "L'angle $\\phi(t)$ peut aussi être tracé:" ] }, { "cell_type": "code", "execution_count": null, "id": "e64f2096", "metadata": {}, "outputs": [], "source": [ "plt.plot(t_rk4, solution_rk4[:, 2], label=\"ϕ\")\n", "phimod2pi = np.ma.array(solution_rk4[:, 2] % (2*np.pi))\n", "masked = np.ma.masked_where(phimod2pi - np.roll(phimod2pi, -1) > 0, phimod2pi) # astuce pour tracer discontinuités\n", "plt.plot(t_rk4, masked, label=\"ϕ mod 2π\")\n", "plt.xlabel(\"t [s]\")\n", "plt.ylabel(\"ϕ(t)\")\n", "plt.legend()\n", "plt.title(\"Angle polaire\")\n", "plt.show()" ] }, { "cell_type": "markdown", "id": "bc8515d8", "metadata": {}, "source": [ "La représentation 3D sert pour illustrer la précession de la trajectoire:" ] }, { "cell_type": "code", "execution_count": null, "id": "55dc2700", "metadata": {}, "outputs": [], "source": [ "import mpl_toolkits.mplot3d\n", "ax = plt.figure().add_subplot(projection='3d')\n", "t = np.linspace(0, h * N, N + 1)\n", "r = solution_rk4[:, 0]\n", "phi = solution_rk4[:, 2]\n", "z = ca / sa * r\n", "x = r * np.cos(phi)\n", "y = r * np.sin(phi)\n", "X = np.arange(-.12, .12, 0.01)\n", "Y = np.arange(-.12, .12, 0.01)\n", "X, Y = np.meshgrid(X, Y)\n", "R = np.sqrt(X**2 + Y**2)\n", "Z = ca / sa * R\n", "surf = ax.plot_surface(X, Y, Z, linewidth=0, antialiased=False, alpha=0.2)\n", "ax.plot(x, y, z)\n", "plt.title(\"Trajectoire en 3D\")\n", "plt.show()" ] }, { "cell_type": "markdown", "id": "f854041c", "metadata": {}, "source": [ "On termine avec une animation:" ] }, { "cell_type": "code", "execution_count": null, "id": "7d59ed20", "metadata": {}, "outputs": [], "source": [ "import matplotlib.animation as animation, matplotlib\n", "\n", "fig = plt.figure()\n", "ax_anim = fig.add_subplot(projection='3d')\n", "\n", "ax_anim.plot_surface(X, Y, Z, linewidth=0, antialiased=False, alpha=0.2)\n", "ax_anim.view_init(elev=50, azim=30, roll=15)\n", "\n", "def init():\n", " return ax_anim.plot([], [], [])\n", "\n", "def update(frame) :\n", " return ax_anim.plot([x[frame]], [y[frame]], [z[frame]], 'ro')\n", "\n", "ani = animation.FuncAnimation(fig, update, 500, blit=True, interval=10)\n", "\n", "from IPython.display import HTML\n", "matplotlib.rcParams['animation.embed_limit'] = 100.\n", "HTML(ani.to_jshtml())" ] }, { "cell_type": "code", "execution_count": null, "id": "54d35395", "metadata": {}, "outputs": [], "source": [] } ], "metadata": { "kernelspec": { "display_name": "Python 3 (ipykernel)", "language": "python", "name": "python3" }, "language_info": { "codemirror_mode": { "name": "ipython", "version": 3 }, "file_extension": ".py", "mimetype": "text/x-python", "name": "python", "nbconvert_exporter": "python", "pygments_lexer": "ipython3", "version": "3.12.3" } }, "nbformat": 4, "nbformat_minor": 5 }