砂轮数学建模与可视化分析

1. 引言

砂轮是机械加工中重要的磨削工具，广泛应用于金属加工、精密制造等领域。为了深入研究砂轮的几何特性、运动轨迹以及加工过程中的接触关系，建立精确的数学模型是必不可少的。本文基于圆柱体砂轮的几何特征，在三维空间中建立砂轮坐标系，并利用SymPy符号计算库和Matplotlib可视化库，实现了砂轮数学模型的建立与可视化分析。

2. 砂轮数学模型

2.1 坐标系定义

如图2-2所示，我们在三维空间中建立砂轮坐标系 $O_s-X_s{Y}_s{Z}_s$，其中：

• $O_s$ 为砂轮坐标系原点
• $X_s$、$Y_s$、$Z_s$ 为坐标系的三个坐标轴
• 砂轮中心坐标为 $(S_x,S_y,S_z)$
• 砂轮上任意一点的坐标为 $(x_s,y_s,z_s)$

2.2 几何模型

砂轮可以看作由圆心为 $O_s$、半径为 $R$ 的圆向 $Z$ 轴负方向移动 $H$ 的距离所形成的圆柱体。在砂轮坐标系下，砂轮的轮廓由以下数学表达式描述：

前端面：
$$S_{\text{前}}=\{\begin{array}{l}{x}_s^2+y_s^2=R^2\\ z_s=H\end{array}\text{\hspace{1em}(2-1)}$$

后端面：
$$S_{\text{后}}=\{\begin{array}{l}{x}_s^2+y_s^2=R^2\\ z_s=0\end{array}\text{\hspace{1em}(2-2)}$$

厚度方向：
$$S_{\text{厚}}=\{\begin{array}{l}{y}_s=-R\\ x_s=0\\ z_s=H\\ z_s=0\end{array}\text{\hspace{1em}(2-3)}$$

其中，方程(2-3)描述的是砂轮厚度方向的两条边界线，分别位于前端面和后端面上，且与 $y=-R$ 平面相交。

3. SymPy符号建模

3.1 符号定义与方程建立

首先，我们使用SymPy定义砂轮数学模型所需的符号变量，并建立对应的方程表达式。

# 砂轮数学模型符号定义
import sympy as sp

# 定义符号变量
x_s, y_s, z_s = sp.symbols('x_s y_s z_s', real=True)  # 砂轮上任意点坐标
R, H = sp.symbols('R H', positive=True)              # 砂轮半径和高度
S_x, S_y, S_z = sp.symbols('S_x S_y S_z', real=True)  # 砂轮中心坐标

# 定义砂轮轮廓方程
print("砂轮数学模型方程：")
print("="*50)

# 前端面方程
print("1. 前端面方程(S_前)：")
front_equation1 = sp.Eq(x_s**2 + y_s**2, R**2)
front_equation2 = sp.Eq(z_s, H)
print(f"   {front_equation1}")
print(f"   {front_equation2}")
print()

# 后端面方程
print("2. 后端面方程(S_后)：")
back_equation1 = sp.Eq(x_s**2 + y_s**2, R**2)
back_equation2 = sp.Eq(z_s, 0)
print(f"   {back_equation1}")
print(f"   {back_equation2}")
print()

# 厚度方向方程
print("3. 厚度方向方程(S_厚)：")
thickness_equation1 = sp.Eq(y_s, -R)
thickness_equation2 = sp.Eq(x_s, 0)
thickness_equation3 = sp.Eq(z_s, H)
thickness_equation4 = sp.Eq(z_s, 0)
print(f"   {thickness_equation1}")
print(f"   {thickness_equation2}")
print(f"   {thickness_equation3}")
print(f"   {thickness_equation4}")
print()

# 圆柱面方程（砂轮侧面）
print("4. 圆柱面方程：")
cylinder_equation = sp.Eq(x_s**2 + y_s**2, R**2)
print(f"   {cylinder_equation}，其中 0 ≤ z_s ≤ H")

3.2 几何参数计算

基于建立的数学模型，我们可以计算砂轮的各种几何参数，如体积、表面积等。

# 砂轮几何参数计算
import sympy as sp
import numpy as np

# 定义符号
R, H = sp.symbols('R H', positive=True)

# 计算砂轮体积
V = sp.pi * R**2 * H
print(f"砂轮体积 V = πR²H = {V}")

# 计算砂轮表面积
# 侧面积
A_side = 2 * sp.pi * R * H
# 两端面积
A_ends = 2 * sp.pi * R**2
# 总表面积
A_total = A_side + A_ends
print(f"砂轮侧面积 A_侧 = 2πRH = {A_side}")
print(f"砂轮两端面积 A_端 = 2πR² = {A_ends}")
print(f"砂轮总表面积 A_总 = 2πR(R + H) = {A_total}")

# 计算特定参数下的数值
R_val = 5.0
H_val = 10.0
print("\n当R=5.0, H=10.0时：")
print(f"体积 V = {float(V.subs({R: R_val, H: H_val})):.2f}")
print(f"侧面积 A_侧 = {float(A_side.subs({R: R_val, H: H_val})):.2f}")
print(f"两端面积 A_端 = {float(A_ends.subs({R: R_val, H: H_val})):.2f}")
print(f"总表面积 A_总 = {float(A_total.subs({R: R_val, H: H_val})):.2f}")

4. 砂轮模型可视化

4.1 3D完整视图

首先，我们创建砂轮的完整3D可视化视图，展示砂轮的圆柱体形状以及关键几何特征。

# 砂轮3D可视化
import numpy as np
import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import Axes3D

# 参数设置
R_val = 5.0  # 砂轮半径
H_val = 10.0  # 砂轮高度
center = (0, 0, 0)  # 砂轮中心坐标

# 创建图形
fig = plt.figure(figsize=(12, 10))
ax = fig.add_subplot(111, projection='3d')

# 创建圆柱面参数化表示
theta = np.linspace(0, 2*np.pi, 100)
z = np.linspace(0, H_val, 50)
theta_grid, z_grid = np.meshgrid(theta, z)

# 圆柱面坐标
x_cylinder = R_val * np.cos(theta_grid)
y_cylinder = R_val * np.sin(theta_grid)
z_cylinder = z_grid

# 绘制圆柱面
ax.plot_surface(x_cylinder, y_cylinder, z_cylinder, 
                alpha=0.6, rstride=2, cstride=2, 
                color='lightblue', edgecolor='gray', linewidth=0.5)

# 创建前后端面圆
theta_circle = np.linspace(0, 2*np.pi, 100)
x_circle_front = R_val * np.cos(theta_circle)
y_circle_front = R_val * np.sin(theta_circle)
z_circle_front = H_val * np.ones_like(theta_circle)

x_circle_back = R_val * np.cos(theta_circle)
y_circle_back = R_val * np.sin(theta_circle)
z_circle_back = np.zeros_like(theta_circle)

# 绘制前后端面圆
ax.plot(x_circle_front, y_circle_front, z_circle_front, 
        'r-', linewidth=3, label='S_前: 前端面圆')
ax.plot(x_circle_back, y_circle_back, z_circle_back, 
        'g-', linewidth=3, label='S_后: 后端面圆')

# 填充前后端面
theta_fill = np.linspace(0, 2*np.pi, 50)
r_fill = np.linspace(0, R_val, 20)
theta_fill_grid, r_fill_grid = np.meshgrid(theta_fill, r_fill)

x_front = r_fill_grid * np.cos(theta_fill_grid)
y_front = r_fill_grid * np.sin(theta_fill_grid)
z_front = H_val * np.ones_like(x_front)

x_back = r_fill_grid * np.cos(theta_fill_grid)
y_back = r_fill_grid * np.sin(theta_fill_grid)
z_back = np.zeros_like(x_back)

ax.plot_surface(x_front, y_front, z_front, alpha=0.4, color='red')
ax.plot_surface(x_back, y_back, z_back, alpha=0.4, color='green')

# 绘制S_厚线
# 第一条线: y = -R, x = 0, z从0到H
z_line1 = np.linspace(0, H_val, 100)
x_line1 = np.zeros_like(z_line1)
y_line1 = -R_val * np.ones_like(z_line1)
ax.plot(x_line1, y_line1, z_line1, 'b-', linewidth=4, label='S_厚: y=-R, x=0')

# 标记关键点
ax.scatter(0, -R_val, 0, c='blue', s=100, marker='o', label='起点 (0,-R,0)')
ax.scatter(0, -R_val, H_val, c='blue', s=100, marker='s', label='终点 (0,-R,H)')
ax.scatter(0, 0, 0, c='black', s=150, marker='*', label='坐标系原点 O_s')

# 设置坐标轴
ax.set_xlabel('X_s', fontsize=12, labelpad=10)
ax.set_ylabel('Y_s', fontsize=12, labelpad=10)
ax.set_zlabel('Z_s', fontsize=12, labelpad=10)
ax.set_title(f'砂轮3D模型 (R={R_val}, H={H_val})', fontsize=16, pad=20)
ax.legend(loc='upper right', fontsize=10)
ax.grid(True)

# 设置视角
ax.view_init(elev=25, azim=45)

# 设置坐标轴范围
ax.set_xlim(-R_val*1.5, R_val*1.5)
ax.set_ylim(-R_val*1.5, R_val*1.5)
ax.set_zlim(-H_val*0.2, H_val*1.2)

plt.tight_layout()
plt.show()

4.2 多视图分析

为了更全面地理解砂轮的几何特征，我们创建包含多个视图的图形，从不同角度展示砂轮的几何特性。

# 砂轮多视图可视化
import numpy as np
import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import Axes3D
from matplotlib.patches import Circle, Rectangle

# 参数设置
R_val = 5.0
H_val = 10.0

# 创建图形
fig = plt.figure(figsize=(15, 12))

# 1. 3D视图
ax1 = fig.add_subplot(2, 2, 1, projection='3d')

# 创建圆柱面
theta = np.linspace(0, 2*np.pi, 100)
z = np.linspace(0, H_val, 50)
theta_grid, z_grid = np.meshgrid(theta, z)
x_cylinder = R_val * np.cos(theta_grid)
y_cylinder = R_val * np.sin(theta_grid)
z_cylinder = z_grid

# 绘制圆柱面
ax1.plot_surface(x_cylinder, y_cylinder, z_cylinder, 
                  alpha=0.6, rstride=5, cstride=5, 
                  color='lightblue', edgecolor='gray')

# 绘制前后端面圆
theta_circle = np.linspace(0, 2*np.pi, 100)
x_circle = R_val * np.cos(theta_circle)
y_circle = R_val * np.sin(theta_circle)
ax1.plot(x_circle, y_circle, H_val*np.ones_like(theta_circle), 
         'r-', linewidth=2, label='S_前')
ax1.plot(x_circle, y_circle, np.zeros_like(theta_circle), 
         'g-', linewidth=2, label='S_后')

# 绘制S_厚线
z_line = np.linspace(0, H_val, 100)
x_line = np.zeros_like(z_line)
y_line = -R_val * np.ones_like(z_line)
ax1.plot(x_line, y_line, z_line, 'b-', linewidth=3, label='S_厚')

# 设置3D视图
ax1.set_xlabel('X_s')
ax1.set_ylabel('Y_s')
ax1.set_zlabel('Z_s')
ax1.set_title('(a) 3D视图', fontsize=12)
ax1.legend()
ax1.view_init(elev=20, azim=30)

# 2. 顶视图 (XY平面)
ax2 = fig.add_subplot(2, 2, 2)
# 绘制前后端面圆
circle_front = Circle((0, 0), R_val, fill=False, edgecolor='red', 
                      linewidth=2, label='S_前')
circle_back = Circle((0, 0), R_val, fill=False, edgecolor='green', 
                     linewidth=2, linestyle='--', label='S_后')
ax2.add_patch(circle_front)
ax2.add_patch(circle_back)
# 标记S_厚点
ax2.plot(0, -R_val, 'bo', markersize=8, label='S_厚点')
# 设置顶视图
ax2.set_xlim(-R_val*1.3, R_val*1.3)
ax2.set_ylim(-R_val*1.3, R_val*1.3)
ax2.set_aspect('equal', 'box')
ax2.set_xlabel('X_s')
ax2.set_ylabel('Y_s')
ax2.set_title('(b) 顶视图 (XY平面)', fontsize=12)
ax2.legend()
ax2.grid(True)

# 3. 前视图 (XZ平面，Y=0)
ax3 = fig.add_subplot(2, 2, 3)
# 绘制圆柱体前视图轮廓
rect = Rectangle((-R_val, 0), 2*R_val, H_val, 
                 fill=False, edgecolor='black', linewidth=2)
ax3.add_patch(rect)
# 标记前后端面
ax3.plot([-R_val, R_val], [H_val, H_val], 'r-', linewidth=2, label='S_前')
ax3.plot([-R_val, R_val], [0, 0], 'g-', linewidth=2, label='S_后')
# 标记S_厚点
ax3.plot(0, 0, 'bo', markersize=8, label='S_厚起点')
ax3.plot(0, H_val, 'bo', markersize=8, label='S_厚终点')
# 连接S_厚点
ax3.plot([0, 0], [0, H_val], 'b-', linewidth=2, label='S_厚')
# 设置前视图
ax3.set_xlim(-R_val*1.3, R_val*1.3)
ax3.set_ylim(-H_val*0.1, H_val*1.1)
ax3.set_xlabel('X_s')
ax3.set_ylabel('Z_s')
ax3.set_title('(c) 前视图 (XZ平面, Y=0)', fontsize=12)
ax3.legend()
ax3.grid(True)
ax3.set_aspect('equal', 'box')

# 4. 侧视图 (YZ平面，X=0)
ax4 = fig.add_subplot(2, 2, 4)
# 绘制圆柱体侧视图轮廓
rect_side = Rectangle((-R_val, 0), 2*R_val, H_val, 
                      fill=False, edgecolor='black', linewidth=2)
ax4.add_patch(rect_side)
# 标记前后端面
ax4.plot([-R_val, R_val], [H_val, H_val], 'r-', linewidth=2, label='S_前')
ax4.plot([-R_val, R_val], [0, 0], 'g-', linewidth=2, label='S_后')
# 标记S_厚点
ax4.plot(-R_val, 0, 'bo', markersize=8, label='S_厚起点')
ax4.plot(-R_val, H_val, 'bo', markersize=8, label='S_厚终点')
# 连接S_厚点
ax4.plot([-R_val, -R_val], [0, H_val], 'b-', linewidth=2, label='S_厚')
# 设置侧视图
ax4.set_xlim(-R_val*1.3, R_val*1.3)
ax4.set_ylim(-H_val*0.1, H_val*1.1)
ax4.set_xlabel('Y_s')
ax4.set_ylabel('Z_s')
ax4.set_title('(d) 侧视图 (YZ平面, X=0)', fontsize=12)
ax4.legend()
ax4.grid(True)
ax4.set_aspect('equal', 'box')

plt.suptitle(f'砂轮多视图分析 (R={R_val}, H={H_val})', fontsize=16, y=0.98)
plt.tight_layout()
plt.show()

5. 应用扩展：砂轮运动轨迹模拟

砂轮数学模型的一个重要应用是模拟其在加工过程中的运动轨迹。这里我们模拟砂轮绕Z轴旋转和平移的简单运动。

# 砂轮运动轨迹模拟
import numpy as np
import matplotlib.pyplot as plt
from matplotlib.animation import FuncAnimation
from mpl_toolkits.mplot3d import Axes3D

# 参数设置
R_val = 5.0
H_val = 8.0
num_frames = 100  # 动画帧数

# 创建砂轮初始位置
def create_wheel(center, R, H):
    """创建砂轮网格"""
    # 圆柱面
    theta = np.linspace(0, 2*np.pi, 50)
    z = np.linspace(0, H, 20)
    theta_grid, z_grid = np.meshgrid(theta, z)
    
    x = center[0] + R * np.cos(theta_grid)
    y = center[1] + R * np.sin(theta_grid)
    z = center[2] + z_grid
    
    return x, y, z

# 创建前端面圆
def create_front_face(center, R, H):
    """创建前端面圆"""
    theta = np.linspace(0, 2*np.pi, 100)
    x = center[0] + R * np.cos(theta)
    y = center[1] + R * np.sin(theta)
    z = center[2] + H * np.ones_like(theta)
    
    return x, y, z

# 创建后端面圆
def create_back_face(center, R, H):
    """创建后端面圆"""
    theta = np.linspace(0, 2*np.pi, 100)
    x = center[0] + R * np.cos(theta)
    y = center[1] + R * np.sin(theta)
    z = center[2] + np.zeros_like(theta)
    
    return x, y, z

# 创建S_厚线
def create_thickness_line(center, R, H, num_points=50):
    """创建S_厚线"""
    z = np.linspace(0, H, num_points)
    x = center[0] + np.zeros_like(z)
    y = center[1] - R * np.ones_like(z)
    z = center[2] + z
    
    return x, y, z

# 创建图形
fig = plt.figure(figsize=(12, 8))
ax = fig.add_subplot(111, projection='3d')

# 初始砂轮位置
center = [0, 0, 0]

# 创建砂轮
x_wheel, y_wheel, z_wheel = create_wheel(center, R_val, H_val)
x_front, y_front, z_front = create_front_face(center, R_val, H_val)
x_back, y_back, z_back = create_back_face(center, R_val, H_val)
x_line, y_line, z_line = create_thickness_line(center, R_val, H_val)

# 初始化绘图对象
wheel_surface = ax.plot_surface(x_wheel, y_wheel, z_wheel, alpha=0.6, color='lightblue')
front_line, = ax.plot(x_front, y_front, z_front, 'r-', linewidth=2, label='前端面')
back_line, = ax.plot(x_back, y_back, z_back, 'g-', linewidth=2, label='后端面')
thickness_line, = ax.plot(x_line, y_line, z_line, 'b-', linewidth=3, label='S_厚线')

# 设置坐标轴
ax.set_xlabel('X')
ax.set_ylabel('Y')
ax.set_zlabel('Z')
ax.set_title('砂轮运动轨迹模拟', fontsize=14)
ax.legend()
ax.set_xlim(-10, 10)
ax.set_ylim(-10, 10)
ax.set_zlim(-5, 15)
ax.view_init(elev=25, azim=30)

# 运动参数
rotation_speed = 0.1  # 旋转速度 (弧度/帧)
translation_speed = 0.2  # 平移速度

# 更新函数
def update(frame):
    """更新砂轮位置"""
    # 计算当前帧的位置
    t = frame / 10.0
    
    # 旋转角度
    angle = rotation_speed * frame
    
    # 平移距离
    dx = translation_speed * np.sin(t)
    dy = translation_speed * np.cos(t)
    dz = 0.1 * np.sin(2*t)
    
    # 更新砂轮中心
    new_center = [dx, dy, dz]
    
    # 创建旋转矩阵
    cos_a = np.cos(angle)
    sin_a = np.sin(angle)
    
    # 创建新的砂轮
    x_wheel_new, y_wheel_new, z_wheel_new = create_wheel(new_center, R_val, H_val)
    x_front_new, y_front_new, z_front_new = create_front_face(new_center, R_val, H_val)
    x_back_new, y_back_new, z_back_new = create_back_face(new_center, R_val, H_val)
    x_line_new, y_line_new, z_line_new = create_thickness_line(new_center, R_val, H_val)
    
    # 应用旋转
    # 圆柱面
    x_rot = x_wheel_new * cos_a - y_wheel_new * sin_a
    y_rot = x_wheel_new * sin_a + y_wheel_new * cos_a
    z_rot = z_wheel_new
    
    # 前端面
    x_front_rot = x_front_new * cos_a - y_front_new * sin_a
    y_front_rot = x_front_new * sin_a + y_front_new * cos_a
    z_front_rot = z_front_new
    
    # 后端面
    x_back_rot = x_back_new * cos_a - y_back_new * sin_a
    y_back_rot = x_back_new * sin_a + y_back_new * cos_a
    z_back_rot = z_back_new
    
    # S_厚线
    x_line_rot = x_line_new * cos_a - y_line_new * sin_a
    y_line_rot = x_line_new * sin_a + y_line_new * cos_a
    z_line_rot = z_line_new
    
    # 更新绘图数据
    wheel_surface.remove()
    wheel_surface = ax.plot_surface(x_rot, y_rot, z_rot, alpha=0.6, color='lightblue')
    
    front_line.set_data(x_front_rot, y_front_rot)
    front_line.set_3d_properties(z_front_rot)
    
    back_line.set_data(x_back_rot, y_back_rot)
    back_line.set_3d_properties(z_back_rot)
    
    thickness_line.set_data(x_line_rot, y_line_rot)
    thickness_line.set_3d_properties(z_line_rot)
    
    # 更新轨迹点
    if frame > 0:
        # 记录前端面圆心轨迹
        x_center = new_center[0] + 0
        y_center = new_center[1] + 0
        z_center = new_center[2] + H_val
        
        # 应用旋转
        x_center_rot = x_center * cos_a - y_center * sin_a
        y_center_rot = x_center * sin_a + y_center * cos_a
        z_center_rot = z_center
        
        ax.scatter(x_center_rot, y_center_rot, z_center_rot, 
                  c='orange', s=10, alpha=0.5)
    
    return wheel_surface, front_line, back_line, thickness_line

# 创建动画
print("正在生成砂轮运动轨迹动画...")
ani = FuncAnimation(fig, update, frames=num_frames, 
                    interval=50, blit=False, repeat=True)

# 显示动画
plt.tight_layout()
plt.show()

# 保存动画 (可选)
# ani.save('grinding_wheel_motion.gif', writer='pillow', fps=20)

6. 结论

本文建立了圆柱体砂轮的完整数学模型，并利用SymPy和Matplotlib实现了模型的符号计算与可视化分析。主要工作包括：

1. 数学模型建立：基于砂轮的几何特征，建立了包含前端面、后端面和厚度方向的完整数学模型，用方程组(2-1)-(2-3)精确描述了砂轮的几何轮廓。

2. 符号计算：利用SymPy实现了砂轮模型的符号表示，能够进行参数化分析和几何计算，包括体积、表面积等关键几何参数的计算。

3. 三维可视化：通过Matplotlib创建了多种可视化视图，包括3D完整视图、多视图分析、参数化曲面表示等，直观展示了砂轮的几何特征。

4. 截面分析：分析了砂轮在不同平面上的截面形状，深入理解了砂轮的几何特性。

5. 几何关系验证：通过随机点采样和蒙特卡洛方法验证了数学模型的正确性，确保了几何关系的一致性。

6. 运动轨迹模拟：扩展了砂轮模型的应用，模拟了砂轮在加工过程中的旋转和平移运动。

该模型为砂轮的进一步研究，如磨削力分析、运动轨迹规划、加工精度控制等提供了理论基础和可视化工具。通过调整模型参数，可以适应不同尺寸和形状的砂轮，为砂轮设计和磨削工艺优化提供了有力的分析手段。