如何调整我的篮球-FT（罚球）模拟以找到关节上的必要扭矩

Question

我正在模拟篮球投掷，它以角运动开始并进入抛射运动。 Physics set-up of simulation

我的目标是深入了解投掷篮球时施加在关节（如肘部和肩部）上的扭矩量。在我的模拟中，扭矩和释放角度是输入，轨迹是输出。我想对输入进行“tweek”处理，以获得可以扫网（得分）的球轨迹。 现在我将它的范围缩小到肘关节，所以基本上是一个弹射器，如下面的代码所示，而且篮子还没有在那里。

但在此扩展之前，我想针对多个“释放角度”和不同的“肘部扭矩”运行模拟。 如您所见，我尝试在第 11 行和第 17 行（已注释掉）中创建一个包含扭矩和释放角度的附加列表，并且我想添加另一个循环，以便整个模拟使用不同的 Angles[re_angl_deg] 列表和 Torques[torq_elb ] 列表作为输入。 不幸的是，这并没有真正奏效。 现在我的问题是：还有其他方法可以让模拟运行多次，每次都有不同的角度和扭矩（就像我制作的列表一样）。

希望有人能给我一些建议！问候.TF

import numpy as np
import matplotlib.pyplot as plt
import math

# Define settings.  
endTime = 60 # The time (seconds) that we simulate.  
dt = 0.001 # The time step (seconds) that we use in the discretization.  
w0 = 0 # The initial velocity [m/s].  
tet0 = 15 # The initial position [m].  
torq_elb = 70     #[Nm]
#torq_elb = np.array([50,55,60,65,70,75,80,85,90,95,100])       #[Nm]
I = 0.16        #[kgm^2]  inertia
x_pos = 0       #initial x position
y_pos = 0       # initial y Position
r = 0.3         # lower arm length [m]
re_angl_deg = 50  # release angle degrees
#re_angl_deg = np.array([30,32,34,36,38,40,42,44,45,47,49,50,52,54])
re_angl_rad = math.radians(re_angl_deg)
g = 9.81 # The gravitational acceleration [m/s^2].    
m_ball = 0.6 # The mass of the ball [kg]. 
rho_air = 1.2041 #[kg/m3]
c_drag =  0.3       # drag coefficient
area_b = (4*math.pi*r**2)/2 # surface area of the ball

 
# Set up variables.  
time = np.arange(0, endTime + dt, dt) # A list with all times we want to plot at.  
w = np.zeros(len(time)) # A list for the angular velocity. [rad/s]
w_deg = np.zeros(len(time))   # list for angle of angular velocity 
tet = np.zeros(len(time)) # A list for the distance.  
x_pos = np.zeros(len(time)) # A list for the xpos
y_pos = np.zeros(len(time)) # A list for the ypos.
x_pos_released = np.zeros(len(time)) # A list for the xpos when released
y_pos_released = np.zeros(len(time)) # A list for the ypos when released
v_x = np.zeros(len(time)) # A list for the xpos
v_y = np.zeros(len(time)) # A list for the ypos.
v_t = np.zeros(len(time)) # List for tangent speed
v_ang = np.zeros(len(time))  # list for angle of the tangent speed
air_res = np.zeros(len(time))    # list for total air resistance
air_res_x = np.zeros(len(time)) # list for X vector of air resistance
air_res_y = np.zeros(len(time))     # list for Y vector of air resistance

w[0] = w0  #start w
tet[0] = math.radians(tet0)  
x_pos[0] = 0   # xpos start
y_pos[0] = 1.8  # approximate hight of a person that is throwing
released=False
  
# Run simulation.
for i in range(1, len(time)):  
    if (released==False):
        w[i] = w[i-1] + torq_elb / I *dt  
        tet[i] = tet[i-1] + w[i] * dt  
        x_pos[i] = x_pos[i-1] + r * math.cos(w[i] - w[i-1]) * dt
        y_pos[i] = y_pos[i-1] + r * math.sin(w[i] - w[i-1]) * dt
    if (tet[i]  > re_angl_rad) and (released==False):
        v_release = w[i] * r
        tet_release = tet[i]
        v_x[i-1] = v_release * math.cos(tet_release)
        v_y[i-1] = v_release * math.sin(tet_release)
        x_pos_released[i-1]=x_pos[i]
        y_pos_released[i-1]=y_pos[i]
        released=True
        w[i] = 0       
    if (released==True):
        v_x[i]=v_x[i-1]
        v_y[i]=v_y[i-1]
        v_t[i] = math.sqrt(v_x[i-1]**2+v_y[i-1]**2)
        v_ang[i] = math.atan(v_y[i-1]/v_x[i-1])
        air_res[i] = 0.5*rho_air*v_t[i]**2*c_drag*area_b ##Force airresitance    
        air_res_x[i] = (air_res[i] * math.cos(v_ang[i]))/m_ball
        air_res_y[i] = (air_res[i] * math.sin(v_ang[i]))/m_ball
        v_x[i]= v_x[i-1] - air_res_x[i] * dt
        v_y[i]= v_y[i-1] - (g + air_res_y[i]) * dt 
        x_pos_released[i] = x_pos_released[i-1] + ((v_x[i-1] + v_x[i])/2) * dt
        y_pos_released[i] = y_pos_released[i-1] + ((v_y[i-1] + v_y[i])/2) * dt
          
# Display results.  
plt.plot(time, tet, color='pink',label='Angular Displacement') 
plt.plot(time, w, color='yellow', label='Angular Velocity')  
plt.plot(time, y_pos_released,color='blue',label='Y_position')
plt.plot(time, x_pos_released,color='purple',label='X_position')
plt.plot(time, v_t, color='black',label='Tangent Velocity')
plt.plot(time, air_res, color='cyan',label='Force Air_res')  
plt.xlabel('Time [s]')  
plt.ylabel('YYYY') 
plt.legend()
plt.xlim(0,1)
plt.ylim(0,10) 
plt.show()  

plt.plot(x_pos_released,y_pos_released,':',color='green' ,label='Trajectory ball')
plt.xlabel('Distance[m]')  
plt.ylabel('Height[m]') 
plt.xlim(0,4)
plt.ylim(0,4)
plt.legend() 
plt.show()

Answer 1

我可能会将每个模拟的值存储在字典中，这样您就可以将其写入文件，然后只需调用特定的释放角度和扭矩来显示，而不是尝试绘制每一个。但是你只需要在这里计算出嵌套循环的逻辑即可。

然后您可以遍历结果以获取图表或将其保留为输入，以便您可以调用特定组合来输出。您需要pip install choice 来实现我的解决方案的输入部分。

代码：

import numpy as np
import matplotlib.pyplot as plt
import math

# Define settings.  
endTime = 60 # The time (seconds) that we simulate.  
dt = 0.001 # The time step (seconds) that we use in the discretization.  
w0 = 0 # The initial velocity [m/s].  
tet0 = 15 # The initial position [m].  
#torq_elb = 70     #[Nm]
torq_elb_list = np.array([50,55,60,65,70,75,80,85,90,95,100])       #[Nm]
I = 0.16        #[kgm^2]  inertia
x_pos = 0       #initial x position
y_pos = 0       # initial y Position
r = 0.3         # lower arm length [m]
#re_angl_deg = 50  # release angle degrees
re_angl_deg_list = np.array([30,32,34,36,38,40,42,44,45,47,49,50,52,54])
g = 9.81 # The gravitational acceleration [m/s^2].    
m_ball = 0.6 # The mass of the ball [kg]. 
rho_air = 1.2041 #[kg/m3]
c_drag =  0.3       # drag coefficient
area_b = (4*math.pi*r**2)/2 # surface area of the ball


results = {}
count = 1
tot = len(torq_elb_list) * len(re_angl_deg_list)
for torq_elb in torq_elb_list:
    for re_angl_deg in re_angl_deg_list:
        results[(torq_elb, re_angl_deg)] = {}
 
        re_angl_rad = math.radians(re_angl_deg ) 
 
        # Set up variables.  
        time = np.arange(0, endTime + dt, dt) # A list with all times we want to plot at.  
        w = np.zeros(len(time)) # A list for the angular velocity. [rad/s]
        w_deg = np.zeros(len(time))   # list for angle of angular velocity 
        tet = np.zeros(len(time)) # A list for the distance.  
        x_pos = np.zeros(len(time)) # A list for the xpos
        y_pos = np.zeros(len(time)) # A list for the ypos.
        x_pos_released = np.zeros(len(time)) # A list for the xpos when released
        y_pos_released = np.zeros(len(time)) # A list for the ypos when released
        v_x = np.zeros(len(time)) # A list for the xpos
        v_y = np.zeros(len(time)) # A list for the ypos.
        v_t = np.zeros(len(time)) # List for tangent speed
        v_ang = np.zeros(len(time))  # list for angle of the tangent speed
        air_res = np.zeros(len(time))    # list for total air resistance
        air_res_x = np.zeros(len(time)) # list for X vector of air resistance
        air_res_y = np.zeros(len(time))     # list for Y vector of air resistance
        
        w[0] = w0  #start w
        tet[0] = math.radians(tet0)  
        x_pos[0] = 0   # xpos start
        y_pos[0] = 1.8  # approximate hight of a person that is throwing
        released=False

  
        # Run simulation.
        for i in range(1, len(time)):  
            if (released==False):
                w[i] = w[i-1] + torq_elb / I *dt  
                tet[i] = tet[i-1] + w[i] * dt  
                x_pos[i] = x_pos[i-1] + r * math.cos(w[i] - w[i-1]) * dt
                y_pos[i] = y_pos[i-1] + r * math.sin(w[i] - w[i-1]) * dt
            if (tet[i]  > re_angl_rad) and (released==False):
                v_release = w[i] * r
                tet_release = tet[i]
                v_x[i-1] = v_release * math.cos(tet_release)
                v_y[i-1] = v_release * math.sin(tet_release)
                x_pos_released[i-1]=x_pos[i]
                y_pos_released[i-1]=y_pos[i]
                released=True
                w[i] = 0       
            if (released==True):
                v_x[i]=v_x[i-1]
                v_y[i]=v_y[i-1]
                v_t[i] = math.sqrt(v_x[i-1]**2+v_y[i-1]**2)
                v_ang[i] = math.atan(v_y[i-1]/v_x[i-1])
                air_res[i] = 0.5*rho_air*v_t[i]**2*c_drag*area_b ##Force airresitance    
                air_res_x[i] = (air_res[i] * math.cos(v_ang[i]))/m_ball
                air_res_y[i] = (air_res[i] * math.sin(v_ang[i]))/m_ball
                v_x[i]= v_x[i-1] - air_res_x[i] * dt
                v_y[i]= v_y[i-1] - (g + air_res_y[i]) * dt 
                x_pos_released[i] = x_pos_released[i-1] + ((v_x[i-1] + v_x[i])/2) * dt
                y_pos_released[i] = y_pos_released[i-1] + ((v_y[i-1] + v_y[i])/2) * dt
                
        # Store simulation results
        results[(torq_elb, re_angl_deg)]['time'] = time
        results[(torq_elb, re_angl_deg)]['tet'] = tet
        results[(torq_elb, re_angl_deg)]['w'] = w
        results[(torq_elb, re_angl_deg)]['y_pos_released'] = y_pos_released
        results[(torq_elb, re_angl_deg)]['x_pos_released'] = x_pos_released
        results[(torq_elb, re_angl_deg)]['v_t'] = v_t
        results[(torq_elb, re_angl_deg)]['air_res'] = air_res
        
        print(f'Finished simulation: torque: {torq_elb}    release: {re_angl_deg}\t{count} of {tot}')
        count+=1
          
# Display results for a particular (torq, re_angl_deg).

#pip install choice
import choice

print('Choose torque.')
torque_choice = int(choice.Menu([str(x) for x in torq_elb_list.tolist()]).ask())

print('Choose release degree.')
release_choice = int(choice.Menu([str(x) for x in re_angl_deg_list.tolist()]).ask())


time = results[(torque_choice, release_choice)]['time']
tet = results[(torque_choice, release_choice)]['tet']
w = results[(torque_choice, release_choice)]['w']
y_pos_released = results[(torque_choice, release_choice)]['y_pos_released']
x_pos_released = results[(torque_choice, release_choice)]['x_pos_released']
v_t = results[(torque_choice, release_choice)]['v_t']
air_res = results[(torque_choice, release_choice)]['air_res']

  
plt.plot(time, tet, color='pink',label='Angular Displacement') 
plt.plot(time, w, color='yellow', label='Angular Velocity')  
plt.plot(time, y_pos_released,color='blue',label='Y_position')
plt.plot(time, x_pos_released,color='purple',label='X_position')
plt.plot(time, v_t, color='black',label='Tangent Velocity')
plt.plot(time, air_res, color='cyan',label='Force Air_res')  
plt.xlabel('Time [s]')  
plt.ylabel('YYYY') 
plt.title('Torque: %s - Release: %s°' %(torque_choice,release_choice))
plt.legend()
plt.xlim(0,1)
plt.ylim(0,10) 
plt.show()  

plt.plot(x_pos_released,y_pos_released,':',color='green' ,label='Trajectory ball')
plt.xlabel('Distance[m]')  
plt.ylabel('Height[m]') 
plt.title('Torque: %s - Release: %s°' %(torque_choice,release_choice))
plt.xlim(0,4)
plt.ylim(0,4)
plt.legend() 
plt.show()

输入：

Choose torque.
Make a choice:
 0: 50
 1: 55
 2: 60
 3: 65
 4: 70
 5: 75
 6: 80
 7: 85
 8: 90
 9: 95

Enter number or name; return for next page

? 0

Choose release degree.
Make a choice:
 0: 30
 1: 32
 2: 34
 3: 36
 4: 38
 5: 40
 6: 42
 7: 44
 8: 45
 9: 47

Enter number or name; return for next page

? 4

输出：