从质点运动到连续介质：用Python模拟刚体旋转与变形（附完整代码）

从质点运动到连续介质用Python模拟刚体旋转与变形附完整代码刚体旋转与变形模拟是计算力学和工程仿真中的基础课题。想象你手中握着一块橡皮当你旋转它时整体形状保持不变而当你挤压它时内部各点的相对位置会发生变化——这两种运动模式分别对应刚体运动与变形运动。本文将用Python带你实现这两种运动的数值模拟通过代码直观理解连续介质力学的核心概念。1. 环境准备与基础概念在开始编码前我们需要明确几个关键概念。刚体运动包含平移和旋转变形运动则涉及质点间相对位置的变化。拉格朗日描述追踪特定质点的运动轨迹而欧拉描述则关注固定空间点上流动的质点。安装必要的Python库pip install numpy sympy matplotlib ipywidgets基础数学工具配置import numpy as np import sympy as sp from sympy import Matrix, symbols, cos, sin, diff import matplotlib.pyplot as plt from matplotlib.animation import FuncAnimation from IPython.display import HTML定义坐标系变换的核心参数# 定义符号变量 X1, X2, X3 symbols(X1 X2 X3) # 物质坐标 x1, x2, x3 symbols(x1 x2 x3) # 空间坐标 t symbols(t) # 时间变量 theta symbols(theta) # 旋转角度2. 刚体旋转的数值实现2.1 二维刚体旋转理论考虑边长为b的正方形平板绕原点逆时针旋转30°的情况。刚体运动方程可表示为def rigid_motion_2d(X, Y, angle_deg): angle_rad np.deg2rad(angle_deg) rotation_matrix np.array([ [np.cos(angle_rad), -np.sin(angle_rad)], [np.sin(angle_rad), np.cos(angle_rad)] ]) return rotation_matrix np.array([X, Y])验证旋转矩阵的正交性Q Matrix([ [cos(theta), -sin(theta), 0], [sin(theta), cos(theta), 0], [0, 0, 1] ]) QT Q.T # 验证Q^T Q I sp.simplify(QT * Q) # 应输出单位矩阵2.2 可视化实现创建正方形平板并应用旋转def plot_rigid_rotation(b2, angle30): # 创建初始正方形 corners np.array([[0,0], [b,0], [b,b], [0,b], [0,0]]) # 应用旋转 rotated np.array([rigid_motion_2d(p[0]-b/2, p[1]-b/2, angle) for p in corners]) rotated np.array([b/2, b/2]) # 平移回中心 # 绘图 fig, ax plt.subplots(figsize(8,8)) ax.plot(corners[:,0], corners[:,1], b-, label初始构形) ax.plot(rotated[:,0], rotated[:,1], r--, label旋转后构形) ax.set_aspect(equal) ax.legend() plt.show()执行旋转可视化plot_rigid_rotation(angle30) # 30度旋转3. 变形运动的数值模拟3.1 简单剪切变形考虑以下变形运动def shear_deformation(X, Y, gamma0.5): x X gamma * Y y Y return x, y对应的变形梯度张量# 定义变形梯度F dx/dX X_vec Matrix([X1, X2, X3]) x_vec Matrix([X1 0.5*X2, X2, X3]) # 剪切变形 F x_vec.jacobian(X_vec) print(变形梯度张量F:) sp.pprint(F)3.2 可视化对比def plot_deformation(b2, gamma0.5): # 创建初始网格 x np.linspace(0, b, 10) y np.linspace(0, b, 10) X, Y np.meshgrid(x, y) # 应用变形 x_def, y_def shear_deformation(X, Y, gamma) # 绘图 fig, ax plt.subplots(1, 2, figsize(12,6)) ax[0].scatter(X, Y, cb, s5) ax[0].set_title(初始构形) ax[1].scatter(x_def, y_def, cr, s5) ax[1].set_title(变形后构形) for axi in ax: axi.set_aspect(equal) plt.show()4. 动态可视化与交互探索4.1 旋转动画实现def animate_rotation(): b 2 corners np.array([[0,0], [b,0], [b,b], [0,b], [0,0]]) center np.array([b/2, b/2]) fig, ax plt.subplots(figsize(8,8)) line_initial, ax.plot([], [], b-, label初始) line_rotated, ax.plot([], [], r--, label旋转) ax.set_xlim(-b, 2*b) ax.set_ylim(-b, 2*b) ax.legend() def init(): line_initial.set_data(corners[:,0], corners[:,1]) return line_initial, def update(frame): angle frame % 360 rotated np.array([rigid_motion_2d(p[0]-center[0], p[1]-center[1], angle) for p in corners]) rotated center line_rotated.set_data(rotated[:,0], rotated[:,1]) return line_rotated, ani FuncAnimation(fig, update, framesnp.arange(0, 360, 2), init_funcinit, blitTrue) plt.close() return HTML(ani.to_jshtml())4.2 变形动画实现def animate_shear(): b 2 x np.linspace(0, b, 10) y np.linspace(0, b, 10) X, Y np.meshgrid(x, y) fig, ax plt.subplots(figsize(8,8)) scat ax.scatter([], [], cr, s5) ax.set_xlim(-1, 3) ax.set_ylim(-1, 3) def init(): scat.set_offsets(np.c_[X.ravel(), Y.ravel()]) return scat, def update(frame): gamma frame / 100 x_def, y_def shear_deformation(X, Y, gamma) scat.set_offsets(np.c_[x_def.ravel(), y_def.ravel()]) return scat, ani FuncAnimation(fig, update, framesnp.arange(0, 100, 2), init_funcinit, blitTrue) plt.close() return HTML(ani.to_jshtml())5. 拉格朗日与欧拉描述的代码对比5.1 拉格朗日实现追踪特定质点的运动轨迹def lagrangian_trajectory(X0, Y0, t_max5, dt0.1): 模拟质点在简单剪切流场中的轨迹 times np.arange(0, t_max, dt) x_traj, y_traj [], [] for t in times: x X0 0.3 * Y0 * t # 随时间增加的剪切 y Y0 x_traj.append(x) y_traj.append(y) return times, x_traj, y_traj5.2 欧拉实现固定空间点观察流动def eulerian_observation(x_obs, y_obs, t_max5, dt0.1): 观察固定点(x_obs,y_obs)上不同质点的属性变化 times np.arange(0, t_max, dt) particles [] for t in times: # 求解逆映射: X x - γYt, Y y Y y_obs X x_obs - 0.3 * Y * t particles.append((X, Y)) return times, particles5.3 对比可视化def compare_descriptions(): # 拉格朗日描述 t, x_lag, y_lag lagrangian_trajectory(1, 1) # 欧拉描述 t_eul, particles eulerian_observation(1.5, 1) x_eul [p[0] for p in particles] y_eul [p[1] for p in particles] plt.figure(figsize(10,5)) plt.subplot(121) plt.plot(t, x_lag, b-, labelx(t)) plt.plot(t, y_lag, r-, labely(t)) plt.title(拉格朗日描述质点轨迹) plt.legend() plt.subplot(122) plt.plot(t_eul, x_eul, g--, labelX(t)) plt.plot(t_eul, y_eul, m--, labelY(t)) plt.title(欧拉描述通过固定点的质点) plt.legend() plt.show()6. 工程应用实例柔性机构分析结合上述技术我们分析一个简单柔性机构的变形def compliant_mechanism_analysis(): # 创建柔性铰链模型 L 10 # 长度 w 1 # 宽度 h 0.2 # 厚度 E 2e6 # 弹性模量 F 10 # 作用力 # 简化的变形计算 (实际工程中需用有限元分析) def deformation(x, y): delta F * L**3 / (3 * E * w * h**3) return x, y delta * (3*(x/L)**2 - 2*(x/L)**3) # 可视化 x np.linspace(0, L, 20) y_bottom np.zeros_like(x) y_top np.ones_like(x) * w # 未变形形状 plt.plot(x, y_bottom, k-, linewidth2) plt.plot(x, y_top, k-, linewidth2) # 变形形状 x_def_b, y_def_b deformation(x, y_bottom) x_def_t, y_def_t deformation(x, y_top) plt.plot(x_def_b, y_def_b, r--) plt.plot(x_def_t, y_def_t, r--) plt.title(柔性机构变形分析) plt.axis(equal) plt.show()在实现这些模拟时有几个常见陷阱需要注意旋转矩阵的正交性验证、变形梯度行列式保持正值保证物理合理性、以及拉格朗日与欧拉描述的正确转换。通过逐步调试和可视化验证可以确保模拟结果的物理正确性。