2022IYPT/CUPT第十一题 旋风球（橡皮筋上的球）

题目回顾

• 哈哈，我们拿了一个不太雅观的国三
• 本人是最佳反方。。
• 团队配合不太行的确不太行，好多题没做。关键是学校也没支持。经典的一句概括，别的学校保研大战，我们大学娱乐开怀。别人比赛热情满满，我们比赛炸鸡热腾腾的。比人记住的都是怎么科研，我们记住的是怎么卷北京烤鸭。
• 总结：烤鸭好吃
• 我最后做的是第11题，第12题，第16题。当的队长
• 12题做的的确不咋的。
• 16题引用的是德国佬的一篇文章。
• 11题嘛，力学分析，写个和计算甲有关的文章吧！！
• Connect two metal balls with an elastic band, then twist the elastic band and put the balls on a table. The balls will begin to spin in one direction, then in the other. Explain this phenomenon and investigate how the behaviour of such a “pendulum” depends on the relevant parameters.

力学理论分析

对单个小球做动力学分析

• 采用极坐标系，以橡皮筋中点为坐标原点

 

 

符号	物理意义
A	橡皮筋横截面积
E	橡皮筋弹性模量
L	橡皮筋原长的一半
R	皮筋某时刻长度的一半
𝐺	金属球速度
T	橡皮筋的扭转角
	金属球的滚动摩擦系数

数值计算

• 根据测得的参量，设置步长为0.01， 利用四阶龙格库塔法数值解微分方程
• 对于方程需要的六个初始条件

 代码段

import numpy as np
import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import Axes3D
from math import *
 
def runge_kutta4(df,a,b,h,y0):
    num = len(y0)
    t = np.arange(a,b+h,h)
    w = np.zeros((t.size,num))
    w[0,:] = y0
 
    for i in range(t.size - 1):
        s0 = df(t[i], w[i,:])
        s1 = df(t[i] + h/2.,w[i,:] + h * s0 / 2.)
        s2 = df(t[i] + h/2.,w[i,:] + h * s1 / 2.)
        s3 = df(t[i+1], w[i,:] + h * s2)
        w[i+1,:] = w[i,:] + h * (s0 + 2*(s1+s2) + s3) / 6.

    return t,w
 
def df(t, variables):
    x,y,z,x1,y1,z1 = variables
    
    A = np.zeros((3,3))
    b = np.zeros(3)
    
    A[0,0] = 1
    A[0,1] = 0
    A[0,2] = 0
    
    A[1,0] = 0
    A[1,1] = 1
    A[1,2] = 0
 
    A[2,0] = 0
    A[2,1] = 0
    A[2,2] = 1
 

    r1 = 2 * 10**-3
    r = 14.5 * 10**-3
    m = 88.25 * 10**-3
    G = 53.4
    E = 268153
    k = 0.5
    u = 0.35
    A1 = r1**2
    L = 10 * 10**-2
    I = 0.4 * m * r**2
    g = 9.8
    
    b[0] = -u*g*sin(atan(x1/x/y1))-E*A1/m/L*2*(x-L/2)+(x+r)*x1**2
    b[1] = (-u*g*cos(atan(x1/x/y1))-2*x1*y1)/x 
    b[2] = -G*3.14* r1**3 * r * z/32/I/x + r*u*m*g*cos(atan(x1/x/y1))/I - k*m*g/I
    
    x2,y2,z2 = np.linalg.solve(A, b)
 
    return np.array([x1,y1,z1,x2,y2,z2])
 
a,b = 0.0,15.0
h = 0.01
 
x = 1/2 * 10**-2
y = 0
z = 100*2*3.14
x1 = 0.1
y1 = 0.000001
z1 = 0.000001
 
y0 = np.array([x,y,z,x1,y1,z1])
t,w = runge_kutta4(df, a, b, h, y0)
t = t
x = w[:,0]
y = w[:,1]
z = w[:,2]
x1 = w[:,3]
y1 = w[:,4]
z1 = w[:,5]
 
fig,axes = plt.subplots(3,2)
axes[0][0].plot(t,x)
axes[0][1].plot(t,y)
axes[1][0].plot(t,z)
axes[1][1].plot(t,x1)
axes[2][0].plot(t,y1)
axes[2][1].plot(t,z1)
plt.show()

•  显然这个结果不太让人满意

从“单摆”出发，拟合这一物理过程

 思路

 

 

直径1.5cm	扭转次数	钟摆运动 持续时间	直径1.8cm	扭转次数	钟摆运动持续时间
光滑的玻璃面	30	13.5s	光滑的玻璃面	30	14.3s
	60	13.8s		60	14.6s
	80	14.1s		80	14.8s
较光滑的木地板	30	12.1s	较光滑的木地板	30	12.7s
	60	12.3s		60	12.9s
	80	12.6s		80	13.2s
粗糙的水泥地面	30	6.8s	粗糙的水泥地面	30	7.8s
	60	6.9s		60	8.0s
	80	7.0s		80	8.2s