📚 课程目录第4章参数化齿轮建模与齿数变化批量仿真

4. 参数化齿轮建模与齿数变化批量仿真

任务定义

本案例创建一个参数化直齿圆柱齿轮模型，以齿数 z、模数 m、压力角 α、齿宽 b 为输入参数，自动生成渐开线齿廓。随后批量改变齿数（z = 18, 25, 32）与模数（m = 2, 3），生成 6 组不同齿轮副，对每组齿轮副执行啮合仿真，提取齿面接触应力峰值，对比分析齿数变化对接触应力的影响。

参数化齿轮几何模型

齿轮渐开线齿廓由基圆展开生成。对于齿数 z、模数 m、压力角 α 的标准直齿轮：分度圆直径 d = m·z，基圆直径 d_b = d·cos α，齿顶圆直径 d_a = d + 2m，齿根圆直径 d_f = d − 2.5m。渐开线上任一点的压力角 θ 满足：

r(θ) = r<sub>b</sub> / cos(θ)，θ ∈ [0, θ<sub>max</sub>]

其中 θ_max = arccos(r_b/r_a)。齿廓由基圆两侧对称渐开线构成，齿根处用圆弧过渡。以下 Python 代码生成单个齿轮的二维截面顶点坐标，并输出为 DXF 格式供仿真前处理使用。

import numpy as np
import math

def gear_profile(z, m, alpha_deg=20.0):
    """生成直齿轮渐开线齿廓顶点坐标（笛卡尔系，单位mm）"""
    alpha = math.radians(alpha_deg)
    r = m * z / 2.0                # 分度圆半径
    rb = r * math.cos(alpha)       # 基圆半径
    ra = r + m                     # 齿顶圆半径
    rf = r - 1.25 * m              # 齿根圆半径
    # 齿厚半角（分度圆上）
    inv_a = math.tan(alpha) - alpha
    half_thick = math.pi / (2*z) + inv_a
    # 渐开线展角范围
    theta_max = math.acos(rb / ra)
    n_pts = 30                     # 每侧渐开线采样点数
    thetas = np.linspace(0, theta_max, n_pts)
    # 右侧渐开线（从基圆到齿顶）
    right_side = []
    for th in thetas:
        rr = rb / math.cos(th)
        # 极角：齿中心线偏置 + 渐开线展开角
        phi = half_thick + th - math.atan2(math.sin(th), math.cos(th))  # 简化：实际为 th - inv
        # 标准渐开线极角公式
        phi = half_thick + th - (math.tan(th) - th) if th > 0 else half_thick
        x = rr * math.cos(phi)
        y = rr * math.sin(phi)
        right_side.append((x, y))
    # 左侧渐开线（镜像）
    left_side = [(-x, y) for (x, y) in reversed(right_side)]
    # 齿根圆弧（连接左右渐开线底部）
    root_pts = 8
    root_arc = []
    for i in range(root_pts + 1):
        ang = math.pi + (math.pi / root_pts) * i  # 从π到2π
        root_arc.append((rf * math.cos(ang), rf * math.sin(ang)))
    # 组装单个齿（从齿根左侧→渐开线→齿顶→右侧→齿根）
    tooth = left_side + [(ra*math.cos(0), ra*math.sin(0))] + right_side + root_arc
    # 旋转复制z个齿
    all_pts = []
    for k in range(z):
        rot = 2 * math.pi * k / z
        c, s = math.cos(rot), math.sin(rot)
        for (x, y) in tooth:
            all_pts.append((x*c - y*s, x*s + y*c))
    return np.array(all_pts)

# 示例：生成 z=20, m=2 的齿轮轮廓并保存
pts = gear_profile(20, 2.0)
with open('gear_z20_m2.dxf', 'w') as f:
    f.write('0\nSECTION\n2\nHEADER\n0\nENDSEC\n0\nSECTION\n2\nENTITIES\n')
    for (x, y) in pts:
        f.write(f'0\nLINE\n8\n0\n10\n{x:.6f}\n20\n{y:.6f}\n11\n{x:.6f}\n21\n{y:.6f}\n')
    f.write('0\nENDSEC\n0\nEOF\n')
print(f'生成齿轮 z=20, m=2，顶点数: {len(pts)}')

该代码输出每个齿的完整闭合轮廓点列。实际仿真前处理时，将相邻顶点用直线连接形成闭合截面，拉伸 b = 20 mm 获得三维实体。齿轮副装配时，主动轮与从动轮中心距 a = m(z₁+z₂)/2，从动轮旋转 180°/z₂ 以保证正确啮合相位。

批量仿真与接触应力分析

以下脚本对 6 组参数组合执行啮合仿真。每组仿真中，主动轮以恒定角速度 ω = 10 rad/s 旋转，从动轮施加 50 N·m 阻力矩。材料设为结构钢（E=210 GPa, ν=0.3），接触算法采用增广拉格朗日法。仿真持续 0.2 s（覆盖约 3 个齿啮合周期），提取齿面最大接触应力。

import subprocess, csv, os

# 参数组合：齿数对 (z1, z2)，模数 m
combos = [
    (18, 18, 2.0), (25, 25, 2.0), (32, 32, 2.0),
    (18, 18, 3.0), (25, 25, 3.0), (32, 32, 3.0)
]

results = []
for z1, z2, m in combos:
    # 生成两个齿轮的几何文件
    pts1 = gear_profile(z1, m)
    pts2 = gear_profile(z2, m)
    # 此处调用有限元求解器（如 CalculiX/ANSYS 批处理）
    # 实际命令因求解器而异，以下为示意接口
    inp_file = f'gear_{z1}_{m}.inp'
    with open(inp_file, 'w') as f:
        f.write(f'*NODE\n')  # 实际需写入节点/单元/接触定义
        # ... 完整有限元模型由前处理脚本生成（略）
    # 运行求解（示意）
    # subprocess.run(['ccx', inp_file], check=True)
    # 读取结果文件中的最大接触应力
    # s_max = parse_result(f'gear_{z1}_{m}.dat')
    # results.append((z1, m, s_max))

# 由于未实际运行求解器，此处不输出伪造的应力数值
print('批量仿真任务已提交，共', len(combos), '组参数组合')

验证方法与判定标准

几何验证：运行 gear_profile() 生成 z=20, m=2 的轮廓，检查齿顶圆直径实测为 44 mm（= (20+2)×2），齿根圆直径 35 mm（= (20−2.5)×2），相邻齿距角 18°。将 DXF 导入 CAD 软件测量分度圆齿厚，理论值为 πm/2 = 3.1416 mm。

仿真验证：每组参数仿真完成后，提取接触应力云图中齿面节点的最大等效应力。判定标准：

接触应力峰值出现在啮合线附近，且随齿数 z 增大而减小（齿数多则重合度大，载荷分布更均匀）；
模数 m 从 2 增至 3 时，接触应力显著下降（齿廓曲率半径增大，赫兹接触应力降低）；
所有模型求解收敛，无穿透或负体积错误。

将 6 组 σ_max 值绘制成以 z 为横轴、m 为分组的折线图，观察两条曲线均呈下降趋势，且 m=3 曲线整体低于 m=2 曲线，即验证参数化批量仿真流程正确。