简介：本资源为2021年全国大学生物理实验竞赛一等奖获奖项目——威尔伯福斯摆（Wilberforce Pendulum）的完整开源实现，面向物理类本科生、实验课程教师及对振动与耦合动力学感兴趣的科研初学者。项目聚焦共振耦合现象，系统呈现振动能量在纵向与横向模式间周期性转换的物理机制，助力理解振动理论、能量守恒、傅里叶分析等核心概念。压缩包共94个文件，含39个JavaScript前端交互模块（用于数据可视化与参数调节）、5个Python脚本（承担数值模拟与数据分析）、16个XML配置与界面定义文件、7个CSS/JS样式资源及配套文档（README.md、LICENSE、ss.md等），整体仅912KB，轻量易部署。目前已有293人学习下载，内容结构清晰：backend含Flask服务与模型逻辑，src为前端工程，cal.py与app.py体现关键算法与控制流程，配套说明文档详述实验原理与复现实操路径，可直接用于课程演示、创新实验复现或毕业设计参考。
1. 威尔伯福斯摆不是“两个摆的简单组合”，而是耦合振动系统的教科书级实证载体
2021年全国大学生物理实验竞赛一等奖作品.zip 解压后核心是wilberforce_pendulum_publish-master，它并非一个演示动画或仿真课件，而是一套完整闭环的物理实验系统：从硬件信号采集（含光电门/加速度计原始数据流）、实时数据处理（Python 后端服务）、动态可视化（React 前端界面），到理论建模验证（含参数辨识与相图重构）。这个项目真正解决的是高校物理实验中长期存在的“现象可观、规律难量、模型难验”三重断层——学生能看见摆动，但无法精确捕捉能量在纵向拉伸模态与横向扭转模态间的周期性转移；教师有理论公式，却缺乏可复现、可调节、可对比的实测数据支撑。它适合两类人：一是需要将《力学》《振动与波》课程实验升级为“可编程物理实验”的高校教师；二是正在准备创新实验设计、需同时体现硬件搭建、数据建模与可视化表达能力的本科生团队。项目中cal.py的相位差计算逻辑、models.py的双自由度微分方程离散化实现、以及app.py中 WebSocket 实时数据管道的设计，共同构成了一个可拆解、可替换、可扩展的物理实验数字底座。
2. 威尔伯福斯摆的物理建模必须从双自由度耦合方程出发，而非单摆近似
2.1 为什么不能直接套用单摆公式？——耦合项的物理意义与数学表征
威尔伯福斯摆的核心特征在于其纵向（轴向）振动与横向（扭转）振动存在强动力学耦合。若忽略耦合项，仅将两部分视为独立单摆，会导致振幅衰减预测偏差超40%，相位关系完全失真。真实系统需用如下双自由度微分方程组描述：
$$ \begin{cases} m\ddot{z} + k_z z + \kappa (\theta - \alpha z) = 0 \ I\ddot{\theta} + k_\theta \theta + \kappa (\theta - \alpha z) = 0 \end{cases} $$
其中 $z$ 为纵向位移，$\theta$ 为扭转角，$m$ 和 $I$ 分别为等效质量与转动惯量，$k_z$、$k_\theta$ 为各自刚度，$\kappa$ 是耦合刚度系数，$\alpha$ 表征几何耦合比例（由连杆长度与摆臂偏心距决定）。项目models.py中的关键实现正是对这一方程组的无量纲化与四阶龙格-库塔离散化：
# models.py 片段：双自由度系统状态空间建模 def wilberforce_ode(t, y, params): """ y = [z, dz/dt, theta, dtheta/dt] params = [m, I, kz, ktheta, kappa, alpha, c_z, c_theta] # 含阻尼项 """ z, vz, theta, vtheta = y m, I, kz, ktheta, kappa, alpha, cz, ctheta = params # 耦合项明确体现在两个方程的右侧：kappa*(theta - alpha*z) dzdt = vz dvzdt = (-kz * z - kappa * (theta - alpha * z) - cz * vz) / m dthetadt = vtheta dvthetadt = (-ktheta * theta - kappa * (theta - alpha * z) - ctheta * vtheta) / I return [dzdt, dvzdt, dthetadt, dvthetadt]提示：
kappa参数在实际标定时不可设为零。项目cal.py中通过扫频实验获取共振峰分裂间距 $\Delta f = f_2 - f_1$，再代入公式 $\kappa \approx \frac{1}{2} m (2\pi \Delta f)^2$ 进行初值估计，这是避免数值求解发散的关键前置步骤。
2.2 实验数据如何反推模型参数？——基于最小二乘的在线辨识流程
理论模型需与实测数据对齐。项目采用分阶段参数辨识策略：先固定几何参数（$\alpha$ 由结构测量获得），再通过三组独立实验数据联合优化剩余6个参数。cal.py中的identify_parameters()函数执行此过程：
# cal.py 片段：多目标参数辨识主函数 def identify_parameters(raw_data_list): """ raw_data_list: [ (t1, z1, theta1), (t2, z2, theta2), ... ] 返回最优参数向量及各目标函数残差 """ def cost_func(params): total_error = 0 for t, z_obs, theta_obs in raw_data_list: # 使用当前params积分ODE，得到模拟z_sim, theta_sim sol = solve_ivp( wilberforce_ode, [t[0], t[-1]], [z_obs[0], 0, theta_obs[0], 0], # 初值取首帧观测值 args=(params,), t_eval=t, method='RK45', rtol=1e-6 ) z_sim, theta_sim = sol.y[0], sol.y[2] # 加权误差：z通道权重0.6，theta通道权重0.4（因光电门精度更高） error_z = np.mean((z_obs - z_sim)**2) error_theta = np.mean((theta_obs - theta_sim)**2) total_error += 0.6 * error_z + 0.4 * error_theta return total_error # 初始猜测：基于文献值与粗略测量 x0 = [0.15, 0.002, 85.0, 12.0, 3.2, 0.8, 0.05, 0.03] # m, I, kz, ktheta, kappa, alpha, cz, ctheta result = minimize(cost_func, x0, method='L-BFGS-B', bounds=[(0.1,0.2), (0.001,0.005), (70,100), (8,15), (1,5), (0.5,1.2), (0.01,0.1), (0.01,0.1)]) return result.x, result.fun # 执行辨识（示例调用） optimal_params, min_error = identify_parameters([ (t_exp1, z_exp1, theta_exp1), (t_exp2, z_exp2, theta_exp2), (t_exp3, z_exp3, theta_exp3) ])该代码逻辑说明：solve_ivp在每次迭代中重新积分微分方程，将模拟输出与三组实测数据比对；bounds参数强制物理合理性（如质量不能为负、阻尼系数需在合理量级）；最终返回的optimal_params可直接写入config.py供后续仿真使用。失败时常见原因包括：初始猜测偏离过远导致minimize收敛至局部极小，此时应检查raw_data_list中时间序列是否对齐、初值是否取自同一时刻。
2.3 模型验证必须通过相图与能量轨迹——而不仅是时域波形拟合
仅比对时域波形（z-t, θ-t）易掩盖相位误差。项目采用相图（Phase Portrait）与模态能量轨迹（Energy Trajectory）双重验证。cal.py中plot_phase_energy()函数生成关键诊断图：
| 验证维度 | 计算方法 | 物理意义 | 正常表现 |
|---|---|---|---|
| z-vz 相图 | plt.plot(z, vz) | 纵向振动能量守恒性 | 封闭椭圆，长轴方向反映 $k_z/m$ |
| θ-vθ 相图 | plt.plot(theta, vtheta) | 扭转振动能量守恒性 | 封闭椭圆，长轴方向反映 $k_\theta/I$ |
| 能量交换轨迹 | E_z = 0.5*m*vz² + 0.5*kz*z²,E_θ = 0.5*I*vθ² + 0.5*kθ*θ²,plt.plot(E_z, E_θ) | 耦合强度与能量守恒 | 近似直线，斜率 ≈ -1，总能量 $E_z+E_θ$ 波动 < 5% |
# cal.py 片段：能量轨迹绘制（关键诊断） def plot_energy_trajectory(z, vz, theta, vtheta, params): m, I, kz, ktheta, *_ = params E_z = 0.5 * m * vz**2 + 0.5 * kz * z**2 E_theta = 0.5 * I * vtheta**2 + 0.5 * ktheta * theta**2 E_total = E_z + E_theta plt.figure(figsize=(12,4)) plt.subplot(131) plt.plot(z, vz, 'b-', alpha=0.7) plt.xlabel('z (m)'); plt.ylabel('vz (m/s)') plt.title('Longitudinal Phase Portrait') plt.subplot(132) plt.plot(theta, vtheta, 'r-', alpha=0.7) plt.xlabel('θ (rad)'); plt.ylabel('vθ (rad/s)') plt.title('Torsional Phase Portrait') plt.subplot(133) plt.plot(E_z, E_theta, 'g-', alpha=0.7) plt.xlabel('E_z (J)'); plt.ylabel('E_θ (J)') plt.title('Energy Exchange Trajectory') plt.grid(True) # 添加总能量波动标注 plt.figtext(0.02, 0.02, f'Total energy std: {np.std(E_total):.3e} J', fontsize=10) plt.tight_layout() plt.show() # 调用示例（使用辨识后的参数） plot_energy_trajectory(z_exp1, vz_exp1, theta_exp1, vtheta_exp1, optimal_params)参数辨识成功的标志是第三幅图中能量轨迹呈高线性度（R² > 0.995）且总能量标准差低于 $10^{-4}$ J。若出现明显弯曲，说明耦合项 $\kappa$ 或阻尼系数估计不足；若总能量持续下降过快，则需增大cz,ctheta。
3. 实时数据管道设计：从传感器原始信号到前端动态相图的低延迟链路
3.1 后端服务如何承载高频振动数据？——WebSocket 与异步任务的协同架构
威尔伯福斯摆典型振动频率在 1–3 Hz，但为精确捕捉相位关系，采样率需达 100 Hz 以上。app.py采用 Flask-SocketIO 构建全双工通信，避免 HTTP 轮询的延迟累积。其核心是分离「数据接收」与「数据广播」两个异步任务：
# app.py 片段：异步数据管道主干 from flask_socketio import SocketIO, emit, join_room import asyncio from threading import Lock socketio = SocketIO(app, async_mode='eventlet', cors_allowed_origins="*") data_buffer = [] # 存储最近1000帧原始数据 buffer_lock = Lock() sampling_rate = 100 # Hz @socketio.on('connect') def handle_connect(): join_room('pendulum_data') # 任务1：模拟传感器数据流（实际项目中替换为串口/USB读取） async def sensor_reader(): """每10ms生成一帧数据：z, theta, timestamp""" t = 0.0 while True: # 模拟真实传感器噪声（高斯白噪声+量化误差） z_raw = 0.02 * np.sin(2*np.pi*1.8*t + 0.1) + np.random.normal(0, 0.0005) theta_raw = 0.015 * np.sin(2*np.pi*1.9*t - 0.3) + np.random.normal(0, 0.0003) frame = { 't': round(t, 3), 'z': round(z_raw, 6), 'theta': round(theta_raw, 6), 'ts': int(time.time() * 1000) # 毫秒级时间戳 } with buffer_lock: data_buffer.append(frame) if len(data_buffer) > 1000: data_buffer.pop(0) await asyncio.sleep(0.01) # 10ms间隔 → 100Hz t += 0.01 # 任务2：向所有客户端广播最新数据（每50ms推送一次，降低前端压力） async def broadcaster(): while True: await asyncio.sleep(0.05) with buffer_lock: if data_buffer: latest = data_buffer[-1] socketio.emit('pendulum_update', latest, room='pendulum_data') # 启动异步任务（Flask-SocketIO eventlet 模式下） @socketio.on('start_stream') def start_stream(): socketio.start_background_task(sensor_reader) socketio.start_background_task(broadcaster)注意：
sensor_reader中await asyncio.sleep(0.01)是硬性时间约束，确保采样间隔稳定。实际部署时需替换为serial.Serial.readline()或pyusb读取，但必须保证循环内无阻塞操作（如文件I/O、数据库查询），否则会拖慢整个事件循环。
3.2 前端如何实现毫秒级相图渲染？——WebGL 加速的 Canvas 动画
src/pages/PendulumView.tsx使用react-konva（基于 HTML5 Canvas）实现轻量级实时绘图，避免 DOM 重排开销。关键优化点在于：只重绘新增点，不刷新整图；使用requestAnimationFrame同步屏幕刷新率：
// src/pages/PendulumView.tsx 片段：高效相图渲染 import { Stage, Layer, Line, Circle } from 'react-konva'; const PendulumView = () => { const [phasePoints, setPhasePoints] = useState<{x: number, y: number}[]>([]); // WebSocket 监听器：仅追加新点，不重建数组 useEffect(() => { const socket = io(); socket.on('pendulum_update', (data: {z: number, theta: number}) => { // 坐标归一化：z∈[-0.03,0.03]→x∈[100,500], theta∈[-0.02,0.02]→y∈[100,400] const x = 100 + ((data.z + 0.03) / 0.06) * 400; const y = 100 + ((0.02 - data.theta) / 0.04) * 300; // Y轴翻转 setPhasePoints(prev => [ ...prev.slice(-199), // 保留最近200点 { x, y } ]); }); return () => socket.close(); }, []); return ( <Stage width={600} height={500}> <Layer> {/* 动态相图曲线 */} <Line points={phasePoints.flatMap(p => [p.x, p.y])} stroke="#3b82f6" strokeWidth={2} tension={0.5} // 样条插值平滑 /> {/* 当前点高亮 */} {phasePoints.length > 0 && ( <Circle x={phasePoints[phasePoints.length-1].x} y={phasePoints[phasePoints.length-1].y} radius={4} fill="#ef4444" /> )} </Layer> </Stage> ); }; export default PendulumView;该实现每帧仅更新points属性，Konva 库内部自动进行 Canvas 路径重绘，实测在 Chrome 中可稳定维持 60 FPS。若需更高性能（如添加频谱图），可切换至regl（WebGL）方案，但本项目复杂度下 Canvas 已足够。
3.3 数据校准与异常过滤：前端预处理保障可视化可靠性
原始传感器数据含直流偏移与脉冲噪声。src/utils/dataProcessor.ts在前端完成轻量级校准，避免后端重复计算：
// src/utils/dataProcessor.ts 片段：前端实时校准 export class DataProcessor { private zOffset: number = 0; private thetaOffset: number = 0; private zHistory: number[] = []; private thetaHistory: number[] = []; // 启动时采集1秒静止数据计算零点偏移 calibrateZero(zStream: number[], thetaStream: number[]): void { this.zOffset = zStream.reduce((a, b) => a + b, 0) / zStream.length; this.thetaOffset = thetaStream.reduce((a, b) => a + b, 0) / thetaStream.length; } // 滑动窗口中位数滤波（抗脉冲噪声） filterOutlier(value: number, history: number[], windowSize: number = 5): number { history.push(value); if (history.length > windowSize) history.shift(); const sorted = [...history].sort((a, b) => a - b); return sorted[Math.floor(sorted.length / 2)]; } process(raw: {z: number, theta: number}): {z: number, theta: number} { let zClean = raw.z - this.zOffset; let thetaClean = raw.theta - this.thetaOffset; // 应用中位数滤波（z通道更易受冲击干扰） zClean = this.filterOutlier(zClean, this.zHistory); thetaClean = this.filterOutlier(thetaClean, this.thetaHistory); // 限幅：超出±5倍标准差则截断（防止传感器饱和） const zStd = Math.sqrt(this.zHistory.reduce((sum, x) => sum + (x - this.zOffset)**2, 0) / this.zHistory.length); const thetaStd = Math.sqrt(this.thetaHistory.reduce((sum, x) => sum + (x - this.thetaOffset)**2, 0) / this.thetaHistory.length); zClean = Math.max(-5*zStd, Math.min(5*zStd, zClean)); thetaClean = Math.max(-5*thetaStd, Math.min(5*thetaStd, thetaClean)); return { z: parseFloat(zClean.toFixed(6)), theta: parseFloat(thetaClean.toFixed(6)) }; } } // 使用示例 const processor = new DataProcessor(); processor.calibrateZero(initialZ, initialTheta); // 静止期调用 socket.on('pendulum_update', (raw) => { const cleaned = processor.process(raw); // 更新相图... });此校准逻辑使相图中心稳定在原点，消除因温度漂移导致的缓慢偏移；中位数滤波有效抑制开关机瞬态、桌面震动等脉冲干扰，避免相图出现异常飞点。
4. 从竞赛作品到教学工具：参数扫描与对比实验的快速构建方法
4.1 一键生成多组对比实验——修改 config.py 即可切换物理场景
项目将所有可调参数集中于config.py，教师无需改代码即可开展探究式教学。例如，研究耦合强度对能量交换周期的影响，只需修改KAPPA值并重启服务：
# config.py 关键参数节选（教学场景常用调整项） # ———————————————————————————————————————————————— # 【基础物理参数】 MASS = 0.15 # kg, 摆锤质量 INERTIA = 0.002 # kg·m², 摆锤转动惯量 K_Z = 85.0 # N/m, 纵向刚度 K_THETA = 12.0 # N·m/rad, 扭转刚度 KAPPA = 3.2 # N·m/rad, 耦合刚度 ← 修改此处！ ALPHA = 0.8 # 无量纲, 几何耦合系数 # 【阻尼参数】 C_Z = 0.05 # N·s/m, 纵向阻尼 C_THETA = 0.03 # N·m·s/rad, 扭转阻尼 # 【仿真与显示】 SIMULATION_DT = 0.005 # s, 数值积分步长 PLOT_WINDOW_SEC = 10.0 # s, 实时绘图时间窗 PHASE_PLOT_SCALE = 200 # px/unit, 相图缩放因子修改KAPPA = 1.0后重启app.py，前端相图将立即显示能量交换周期显著延长（理论周期 $T_{exchange} \propto 1/\kappa$）；设KAPPA = 0.0则两模态完全解耦，相图退化为两个独立椭圆。这种即时反馈极大提升课堂演示效果。
4.2 快速导出科研级数据包——JSON 格式兼容 MATLAB/Python 分析
所有实时采集与仿真数据均按标准 JSON Schema 输出，便于导入主流分析工具。views.py中/api/export接口提供三种导出模式：
# views.py 片段：标准化数据导出 from flask import jsonify, send_file import json import io @app.route('/api/export', methods=['POST']) def export_data(): data_type = request.json.get('type') # 'raw', 'simulated', 'both' duration = request.json.get('duration', 30) # 秒 # 从内存缓冲区或数据库提取指定时长数据 if data_type == 'raw': export_data = get_raw_buffer_last(duration) elif data_type == 'simulated': export_data = run_simulation_for(duration) else: export_data = { 'raw': get_raw_buffer_last(duration), 'simulated': run_simulation_for(duration) } # 生成符合物理实验数据规范的JSON output = { "metadata": { "experiment_id": "WP_2021_AWARD", "timestamp": datetime.now().isoformat(), "sampling_rate_hz": 100, "units": {"z": "m", "theta": "rad", "t": "s"}, "parameters_used": current_config_dict() # 读取config.py当前值 }, "data": export_data } # 内存中生成JSON文件并返回 json_str = json.dumps(output, indent=2) json_bytes = io.BytesIO(json_str.encode('utf-8')) json_bytes.seek(0) return send_file( json_bytes, mimetype='application/json', as_attachment=True, download_name=f'wilberforce_data_{data_type}_{duration}s.json' )教师在课堂演示后，可立即导出wilberforce_data_raw_30s.json，学生用 Pythonpandas.read_json()或 MATLABjsondecode()直接加载，无缝接入傅里叶分析、李雅普诺夫指数计算等高阶实验。
4.3 教学提示卡：三个必做对比实验与预期现象
为降低教学实施门槛，项目附带ss.md（教学提示文档），明确列出三个核心对比实验及其现象判据：
| 实验编号 | 操作步骤 | 预期现象 | 教学要点 |
|---|---|---|---|
| EXP-01 | 将KAPPA从 3.2 降至 0.8，保持其他参数不变 | 能量交换周期延长约 2.2 倍；相图中能量轨迹斜率绝对值减小 | 耦合强度 $\kappa$ 直接决定模态间能量转移速率，验证公式 $T_{exchange} = \frac{2\pi}{\sqrt{(k_z/m - k_\theta/I)^2 + 4\kappa^2}}$ |
| EXP-02 | 将C_Z从 0.05 增至 0.2，C_THETA不变 | z-vz 相图椭圆迅速坍缩为点，θ-vθ 相图仍保持较完整椭圆；总能量衰减中 z 分量占比超 70% | 阻尼不对称性导致能量单向耗散，引申讨论非保守系统的哈密顿量破缺 |
| EXP-03 | 将K_Z与K_THETA设为相等（如均为 50.0），KAPPA=3.2 | 出现简并模态：z 与 θ 振动频率相同，相图呈现完美圆形；能量轨迹变为严格直线 | 简并条件下系统具有额外对称性，是理解量子简并、光子晶体带隙的基础类比 |
每个实验均配有curl命令示例，教师可直接在终端执行触发参数变更，无需打开编辑器：
# EXP-01 示例：降低耦合强度 curl -X POST http://localhost:5000/api/update_config \ -H "Content-Type: application/json" \ -d '{"KAPPA": 0.8}'此设计将复杂的物理概念转化为可触摸、可测量、可重复的操作指令，真正实现“做中学”。
本文还有配套的精品资源，点击获取
