MOUVEMENT BALISTIQUE AVEC FROTTEMENT QUADRATIQUE CCP MP 2017

Recherche de la portée maximale

On cherche à déterminer la partée maximale d'un tir avec frottements quadratiques. Pour cela on commence par calculer la portée du tir étudié dans le fichier Tartaglia_CCP_MP_2017 , puis on cherchera à déterminer l'angle qui maximise cette portée.

In [1]:
%display latex
In [1]:
g=9.81; vo=380; alpha=n((18*pi)/180);    # angle initial 18°
m_1=n(11350*4*pi*(2e-3)^3/3)             #masse des plombs n°1, 5 et 10 en kg
m_5=n(11350*4*pi*(1.5e-3)^3/3)
m_10=n(11350*4*pi*(0.875e-3)^3/3)
beta_1=3.40046e-6   #  beta force en -beta v^2  en unite SI
beta_5=1.912759e-6
beta_10=6.508693e-7

Le détail de la méthode de résolution est dans le fichier Tartaglia_CCP_MP_2017.

In [2]:
var('t')
(x,z,vx,vz)=var('x,z,vx,vz')
sol1 = desolve_system_rk4([vx,vz,-(beta_1/m_1)*sqrt(vx^2+vz^2)*vx,-g-(beta_1/m_1)*sqrt(vx^2+vz^2)*vz],
                          [x,z,vx,vz],
                          ics=[0,0,0,vo*cos(alpha),vo*sin(alpha)],
                          end_points=n(2*vo*sin(alpha)/g),step=n(2*vo*sin(alpha)/g)/100,ivar=t) 

# première estimation de la durée à partir du temps de retombée du mvt sans frottements

La solution renvoie une liste dont chaque élément est de la forme: $[t_i,x(t_i),z(t_i),vx(t_i),vz(t_i)]$.

On a tracé la courbe dans le fichier Tartaglia_CCP_MP_2017 la portée est entre 330 et 335 m. On précise le tracé de la courbe.

In [3]:
graph1=line([(sol1[i][1],sol1[i][2]) for i in range(len(sol1))],color='blue')
show(graph1,xmin=332,xmax=333,ymin=-5,ymax=5,  axes_labels=['$x$','$z$'])   
#même si la variable s'appelle z sur la verticale l'option de tracé s'écrit avec ymin et ymax 

Calcul direct de la portee du tir (valeur de $x$ quand le projectile retombe dans le plan $z=0$.)

In [4]:
i=0   #initialisation 
while sol1[i][2]>=0:
    i+=1
portee_inf=sol1[i-1][1],sol1[i-1][2]  #couple (x,z) pour la dernière valeur de z>0
portee_sup=sol1[i][1],sol1[i][2]  #couple (x,z) pour la première valeur de z<0
print(portee_inf,portee_sup)
((332.3458053838148, 0.4015398108540875), (335.5340487878719, -5.68547863544882))

On obtient un encadrement de $x_p$.

Pour l'évaluer, on approxime la courbe à une droite joignant ces deux points. Elle coupe la droite $z=0$ pour $$x=x_{i-1}-z_{i-1}\frac{x_i-x_{i-1}}{z_i-z_{i-1}}$$ qu'on approximera à la portée $x_p$.

In [5]:
xp=portee_inf[0]-portee_inf[1]*(portee_sup[0]-portee_inf[0])/(portee_sup[1]-portee_inf[1])
print(xp)
332.5561229086778

La valeur trouvée coïncide avec la valeur lue graphiquement.


But: déterminer xp en fonction de $\alpha$ afin de chercher la valeur de l'angle de tir qui maximise la portée.

In [6]:
xp=[] #initialisation de la liste xp
graph=Graphics()  #initialise un graphique vide
liste_theta=[]    #initialisation de la liste des valeurs de theta
pas_theta=n(pi/45)  # pas de 4°
for j in range(1,15):
    theta=j*pas_theta
    sol = desolve_system_rk4([vx,vz,-(beta_1/m_1)*sqrt(vx^2+vz^2)*vx,-g-(beta_1/m_1)*sqrt(vx^2+vz^2)*vz],
                          [x,z,vx,vz],
                          ics=[0,0,0,vo*cos(theta),vo*sin(theta)],
                          end_points=20,step=(20/100),ivar=t) 
    graph+=line([(sol[i][1],sol[i][2]) for i in range(len(sol))], color=hue(float(j)/15))
    i=0   #initialisation 
    while sol[i][2]>=0:
        i+=1
    portee_inf=sol[i-1][1],sol[i-1][2]  #couple (x,z) pour la dernière valeur de z>0
    portee_sup=sol[i][1],sol[i][2]  #couple (x,z) pour la première valeur de z<0
    portee=portee_inf[0]-portee_inf[1]*(portee_sup[0]-portee_inf[0])/(portee_sup[1]-portee_inf[1])
    xp.append(portee)
    liste_theta.append(n((theta*180)/pi))  #liste des theta en degre
print(xp)    
show(graph,aspect_ratio=1,ymin=-100)
[248.43962496387306, 293.06578179600615, 316.3682184999193, 329.79067656981596, 336.94682054292383, 339.44775010098766, 338.11086112538084, 333.4713172561592, 325.8733694645093, 315.5684156597507, 302.70906659268803, 287.4410792783512, 269.8614904921159, 250.04659530692982]
In [19]:
list_plot([[liste_theta[i],xp[i]] for i in range(len(xp))],axes_labels=[r'$\theta ^{\circ}$',r'$x_p\  {\rm (m)}$'],
         gridlines='minor',xmin=0,ymax=350)   
Out[19]:

On reprend les mêmes calculs pour les plombs n°5.

In [12]:
xp5=[] #initialisation de la liste xp
graph5=Graphics()  #initialise un graphique vide
liste_theta=[]    #initialisation de la liste des valeurs de theta
pas_theta=n(pi/45)  # pas de 4°
for j in range(1,15):
    theta=j*pas_theta
    sol5 = desolve_system_rk4([vx,vz,-(beta_5/m_5)*sqrt(vx^2+vz^2)*vx,-g-(beta_5/m_5)*sqrt(vx^2+vz^2)*vz],
                          [x,z,vx,vz],
                          ics=[0,0,0,vo*cos(theta),vo*sin(theta)],
                          end_points=16,step=(16/100),ivar=t) 
    graph5+=line([(sol5[i][1],sol5[i][2]) for i in range(len(sol5))], color=hue(float(j)/15))
    i=0   #initialisation 
    while sol5[i][2]>=0:
        i+=1
    portee5_inf=sol5[i-1][1],sol5[i-1][2]  #couple (x,z) pour la dernière valeur de z>0
    portee5_sup=sol5[i][1],sol5[i][2]  #couple (x,z) pour la première valeur de z<0
    portee5=portee5_inf[0]-portee5_inf[1]*(portee5_sup[0]-portee5_inf[0])/(portee5_sup[1]-portee5_inf[1])
    xp5.append(portee5)
    liste_theta.append(n((theta*180)/pi))  #liste des theta en degre
print(xp5)    
show(graph5,ymin=-100)
[200.58318770105993, 233.73966487386394, 250.8824735398078, 260.5942568457922, 265.56343038195286, 266.9663922709362, 265.46637989061753, 261.43592604717406, 255.14947423952142, 246.78688942135668, 236.48375812161194, 224.33294785107955, 210.42091828049604, 194.80268852750757]
In [17]:
list_plot([[liste_theta[i],xp5[i]] for i in range(len(xp5))],axes_labels=[r'$\theta ^{\circ}$',r'$x_p\  {\rm (m)}$'],
         gridlines='minor',xmin=0,ymax=280)   
Out[17]:

et pour les plombs n°10

In [14]:
xp10=[] #initialisation de la liste xp
graph10=Graphics()  #initialise un graphique vide
liste_theta=[]    #initialisation de la liste des valeurs de theta
pas_theta=n(pi/45)  # pas de 4°
for j in range(1,15):
    theta=j*pas_theta
    sol10 = desolve_system_rk4([vx,vz,-(beta_10/m_10)*sqrt(vx^2+vz^2)*vx,-g-(beta_10/m_10)*sqrt(vx^2+vz^2)*vz],
                          [x,z,vx,vz],
                          ics=[0,0,0,vo*cos(theta),vo*sin(theta)],
                          end_points=14,step=(14/100),ivar=t) 
    graph10+=line([(sol10[i][1],sol10[i][2]) for i in range(len(sol5))], color=hue(float(j)/15))
    i=0   #initialisation 
    while sol10[i][2]>=0:
        i+=1
    portee10_inf=sol10[i-1][1],sol10[i-1][2]  #couple (x,z) pour la dernière valeur de z>0
    portee10_sup=sol10[i][1],sol10[i][2]  #couple (x,z) pour la première valeur de z<0
    portee10=portee10_inf[0]-portee10_inf[1]*(portee10_sup[0]-portee10_inf[0])/(portee10_sup[1]-portee10_inf[1])
    xp10.append(portee10)
    liste_theta.append(n((theta*180)/pi))  #liste des theta en degre
print(xp10)    
show(graph10,ymin=-100)
[127.73130088283992, 146.88262454039545, 156.6258342527939, 162.0224350458136, 164.60946208772134, 165.07997242789756, 163.8171022918599, 161.0498214068402, 156.93501150579922, 151.57903672343974, 145.06503455811273, 137.4508588293882, 128.78581228615957, 119.10571581545973]
In [18]:
list_plot([[liste_theta[i],xp10[i]] for i in range(len(xp10))],axes_labels=[r'$\theta ^{\circ}$',r'$x_p\  {\rm (m)}$'],
         gridlines='minor',xmin=0,ymax=170)   
Out[18]:

Calcul des $v_{\infty}$

In [26]:
v1_inf=n(sqrt(m_1*g/beta_1));v1_inf
Out[26]:
33.1247719729984
In [27]:
v5_inf=n(sqrt(m_5*g/beta_5));v5_inf
Out[27]:
28.6868921484768
In [29]:
v10_inf=n(sqrt(m_10*g/beta_10));v10_inf
Out[29]:
21.9099771664711
In [30]:
m_1
Out[30]:
0.000380342150594604
In [32]:
m_5
Out[32]:
0.000160456844782099
In [33]:
m_10
Out[33]:
0.0000318499408334837
In [ ]: