bstart-trb 0.2.0

Add your description here

圆柱滚子轴承切片计算详解

本文档详细介绍 BSTART 系统中圆柱滚子轴承（CRB - Cylindrical Roller Bearing）的切片计算原理和实现逻辑。

---

目录

• 1. 切片计算概述
• 2. 载荷分布计算
• 3. 滚子切片应力计算
• 4. 输入输出格式
• 5. 使用示例
• 6. 与寿命计算的关系
• 7. 边缘效应与凸度修形
• 8. 滚子倾斜角 tilt 详解
• 9. f2py 封装（Python 直接调用）
• 10. 完整计算流程
• 参考文献

---

1. 切片计算概述

1.1 什么是切片计算

圆柱滚子轴承的接触区域是线接触，滚子与滚道的接触应力沿滚子长度方向分布不均匀。为了精确计算应力分布，将滚子沿轴向分成多个切片（条带），分别计算每个切片的接触应力。

滚子切片示意图：
┌─────────────────────────────────────────────┐
│  1  │  2  │  3  │ ... │ n-1 │  n  │         │ ← 滚子
└─────────────────────────────────────────────┘
  ↑      ↑      ↑           ↑      ↑
  σ₁     σ₂     σ₃         σₙ₋₁   σₙ          ← 各切片应力

1.2 切片计算的必要性

传统方法	切片方法
假设应力均匀分布	计算真实应力分布
忽略边缘效应	捕捉边缘应力集中
保守估计寿命	更精确的寿命预测

1.3 相关文件

文件	功能
CRB_main.py	圆柱滚子轴承计算主程序
CRB_C.py	C 语言 DLL 接口，计算载荷分布
TAPER_FORTRAN.py	Fortran 程序接口，计算切片应力
CRB_load.dll	C 语言编译的载荷分布计算库
TAPER.exe	Fortran 编译的应力计算程序

---

2. 载荷分布计算

2.1 计算原理

在径向载荷作用下，各滚子承受的载荷不同。最大载荷出现在载荷方向（通常为 0° 位置），随角度增加载荷逐渐减小。

载荷分布公式： $$Q_j = Q_{max} \cdot \left(1 - \frac{1}{2\varepsilon}\left(1 - \cos\psi_j\right)\right)^n$$

其中：

• $Q_j$：第 j 个滚子的载荷（N）
• $Q_{max}$：最大滚子载荷（N）
• $\varepsilon$：载荷分布参数
• $\psi_j$：第 j 个滚子的方位角（rad）
• $n$：载荷-变形指数（线接触 $n = 10/9$）

2.2 代码实现

载荷分布通过 C 语言 DLL 计算（CRB_C.py）：

def c_wfk(z: int, pd: float, l: float, fr: float, r1: float, r2: float):
    """
    计算圆柱滚子轴承的载荷分布
    
    参数:
        z: 滚子个数
        pd: 径向游隙 (mm)
        l: 滚子长度 (mm)
        fr: 径向载荷 (N)
        r1: 滚子半径 (mm)
        r2: 滚道半径 (mm)
    
    返回:
        dict: {角度: [接触变形, 接触载荷]}
    """
    clib = ctypes.CDLL('./CRB_load.dll')
    clib.CRB.restype = TestStruct
    
    value = clib.CRB(
        ctypes.c_double(z), 
        ctypes.c_double(pd),
        ctypes.c_double(l), 
        ctypes.c_double(fr),
        ctypes.c_double(r1), 
        ctypes.c_double(r2)
    )
    
    # 解析结果
    c_dict = {}
    for i in range(z):
        angle = value.array[i*3] * 180 / math.pi  # 方位角（度）
        deformation = value.array[i*3 + 1]         # 接触变形
        load = value.array[i*3 + 2]                # 接触载荷
        c_dict[angle] = [deformation, load]
    
    return c_dict

2.3 输出示例

# 14 个滚子，径向载荷 4450N
c_dict = {
    0.0: [0.0023, 1928.5],      # 0° 滚子：最大载荷
    25.7: [0.0019, 1456.2],     # 25.7° 滚子
    51.4: [0.0012, 892.1],      # 51.4° 滚子
    77.1: [0.0005, 324.8],      # 77.1° 滚子
    102.9: [0.0, 0.0],          # 102.9° 滚子：无接触
    # ... 其他滚子
}

---

3. 滚子切片应力计算

3.1 切片参数

参数	说明	默认值	限制
n	切片数（条带数）	30	最大 50（f90）/ 30（原版 FOR）
m1	凸度类型：0-直线, 1-圆弧, 2-对数	2 (对数)	-
m2	滚道类型：0-平面, 1-内圈, 2-外圈	1 或 2	-

⚠️ 切片数限制说明：

• 原版 TAPER_A.FOR 使用固定数组 dimension p(30),a(30),d(30,30)，最大 30 片
• 新版 taper_slice.f90 使用 MAX_SLICES = 50，最大 50 片
• 如需更多切片，需修改源码中的数组大小限制并重新编译

3.2 TAPER 类参数说明

class TAPER:
    def __init__(self, d1, d2, l, alfa, q, tilt, m1, n, m2, dr1, fai, curveQ):
        """
        Fortran 程序接口
        
        参数:
            d1: 滚子小端直径 (mm)
            d2: 滚子大端直径 (mm)（圆柱滚子 d1=d2）
            l: 滚子有效长度 (mm)
            alfa: 滚子半圆锥角 (°)（圆柱滚子 alfa=0）
            q: 滚子载荷 (N)
            tilt: 滚子倾斜角度 (°)
            m1: 凸度类型 - 0:直线, 1:圆弧, 2:对数
            n: 切片数（条带数） ← 关键参数
            m2: 滚道类型 - 0:平面, 1:内圈, 2:外圈
            dr1: 滚道直径 (mm)
            fai: 滚道半圆锥角 (°)
            curveQ: 设计载荷 (N)（用于计算对数凸度）
        """

3.3 应力计算原理

3.3.1 简化的 Hertz 公式（参考）

传统的 Hertz 椭圆接触应力公式为：

$$\sigma_i = \frac{2 \cdot q_i}{\pi \cdot a_i \cdot b_i}$$

其中：

• $\sigma_i$：第 i 个切片的最大接触应力（MPa）
• $q_i$：第 i 个切片的线载荷（N/mm）
• $a_i, b_i$：接触椭圆半轴（mm）

但此公式假设各切片独立，不考虑切片间的弹性耦合。

3.3.2 Fortran 实际实现：影响系数法

TAPER.exe（源码 fortran/TAPER_A.FOR）使用更精确的弹性半空间影响系数法：

核心算法步骤：

1. 建立影响系数矩阵 [d]

   • infmat 子程序（第324-347行）
   • 计算各切片间的弹性耦合关系
   • 使用高斯积分计算影响系数（calf, f 子程序）

2. 求解线性方程组 [d]{p}={s}

   • gauss 子程序（第387-440行）
   • 使用高斯消元法
   • 其中 {s} 为间隙数组，{p} 为接触应力数组

3. 迭代收敛

   • contac 子程序中的双重迭代
   • 外层：载荷平衡迭代（调整变形 dc）
   • 内层：接触区域收敛迭代

关键代码片段：

c     计算影响系数矩阵并求解应力 p()
      call infmat(n,h,a,d)           ! 建立影响系数矩阵
      call gauss(t0,n,ierror,d,p,s)  ! 求解方程组，得到应力

c     根据应力计算新的接触半宽度
      do 24 i=1,n
24    a(i)=1.767e-5*r(i)*p(i)        ! Hertz 关系

其中 t0 = 4*(1-v²)/(π*E) 为材料常数（v=0.3 泊松比，E=206GPa 弹性模量）。

与简化公式的区别：

特性	简化 Hertz 公式	影响系数法
切片耦合	不考虑	通过影响系数矩阵耦合
接触区判断	固定	迭代确定接触/非接触区
精度	近似	精确（弹性半空间理论）
计算复杂度	O(n)	O(n³)

对数凸度修形可以优化应力分布，减少边缘应力集中。

3.3.3 重要澄清：载荷不是均分到各切片

常见误解：如果总载荷 q=3000N，切成 n=100 片，是否每片分 30N？

答案：不是！

Fortran 程序的输入是总载荷 q，而不是每片的载荷。实际计算流程如下：

输入: 总载荷 q=3000N, 切片数 n=100

步骤1: 估算整体变形
       dc = 3.84e-5 × q^0.9 / rl^0.8

步骤2: 计算各切片间隙
       s(i) = dc - z(i)    (z 取决于凸度修形和倾斜)

步骤3: 建立 n×n 影响系数矩阵 [d]
       d(i,j) 表示切片 j 对切片 i 的弹性影响

步骤4: 求解方程组 [d]{p} = {s}
       得到各切片应力 p(i)

步骤5: 计算平衡载荷
       q1 = π × Σ(p(i) × a(i)) × h

步骤6: 如果 |q1 - q| / q > 0.5%
       调整 dc，返回步骤2

输出: 各切片应力 p(i) (MPa) —— 非均匀分布！

为什么各切片应力不均匀？

1. 影响系数矩阵的耦合效应：每个切片的变形受其他切片载荷的影响
2. 边缘效应：边缘切片只有单侧邻居"帮忙承载"，应力更高
3. 凸度修形：故意让边缘间隙变小来补偿边缘效应
4. 倾斜：tilt≠0 时一侧接触紧密、另一侧可能脱离

3.4 代码实现（旧版 - 文件传参）

# CRB_main.py (旧版实现，通过 INPUT_TAPER.dat 文件传参)

class CRB:
    def __init__(self, z, Pd, l, l_total, fr, r1, r2, r3, curveQ):
        """
        圆柱滚子轴承计算类
        
        参数:
            z: 滚子数量
            Pd: 径向游隙 (mm)
            l: 滚子有效长度 (mm)
            l_total: 滚子总长度 (mm)
            fr: 径向载荷 (N)
            r1: 滚子半径 (mm)
            r2: 内滚道半径 (mm)
            r3: 外滚道半径 (mm)
            curveQ: 设计载荷 (N)
        """
        self.n_slice = 30  # 切片数
    
    def CRB_Main(self):
        """主计算流程"""
        # 1. 计算载荷分布
        c_dict, a_list = c_wfk(self.z, self.Pd, self.l, self.fr, self.r1, self.r2)
        
        # 2. 内滚道切片应力计算
        roller_inner = TAPER(
            d1=2*self.r1, d2=2*self.r1,  # 滚子直径
            l=self.l,                     # 滚子有效长度
            alfa=0,                       # 半圆锥角（圆柱=0）
            q=0,                          # 载荷（后续赋值）
            tilt=0,                       # 倾斜角
            m1=2,                         # 对数凸度
            n=30,                         # 切片数 ← 关键参数
            m2=1,                         # 内圈
            dr1=2*self.r2,               # 内滚道直径
            fai=0,                        # 滚道半圆锥角
            curveQ=self.curveQ           # 设计载荷
        )
        
        # 对每个滚子计算切片应力
        line_dict_inner = {}
        for angle, values in c_dict.items():
            load = values[1]  # 该滚子的载荷
            roller_inner.q = load
            roller_inner.write_dat()
            roller_inner.run_exe()
           
            # 读取 30 个切片的应力值
            with open('TAPER.txt', 'r') as f:
                stresses = []
                for line in f.readlines()[1::2]:
                    stresses += line.strip().split()
                line_dict_inner[angle] = stresses
        
        # 3. 外滚道切片应力计算（类似）
        roller_outer = TAPER(..., m2=2, dr1=2*(-self.r3), ...)
        # ...
        
        return c_dict, line_dict_inner, line_dict_outer

---

4. 输入输出格式

4.1 输入文件（INPUT_TAPER.dat）

10.0      10.0      9.6       0.0       1340.859878 
0.0       2         30        2         -150.0    0.0       14100.0   

行	参数	说明
1	d1, d2, l, alfa, q	滚子小端直径, 大端直径, 有效长度, 半圆锥角, 载荷
2	tilt, m1, n, m2, dr1, fai, curveQ	倾斜角, 凸度类型, 切片数, 滚道类型, 滚道直径, 滚道半圆锥角, 设计载荷

为什么有两个半圆锥角（alfa 和 fai）？

1. 滚子半圆锥角 alfa

• 描述滚子自身的圆锥形状
• 圆柱滚子：alfa = 0
• 圆锥滚子：alfa > 0（典型 10°-15°）

2. 滚道半圆锥角 fai

• 描述滚道的圆锥形状
• 圆柱滚子轴承：fai = 0（滚道是圆柱面）
• 圆锥滚子轴承：fai > 0（滚道是圆锥面）

为什么需要分开？

理论上，在圆锥滚子轴承中，滚子和滚道的半圆锥角应该相同（alfa = fai），这样才能实现完美的线接触。

但在实际设计和计算中，分开它们可以：

1. 处理制造误差 - 滚子和滚道的实际半圆锥角可能略有不同
2. 更通用的计算模型 - 同一套代码可以处理圆柱滚子轴承和圆锥滚子轴承
3. 处理特殊接触 - 还可以处理滚子与平面的接触（fai = 0, m2 = 0）

对于圆柱滚子轴承（CRB）：

• alfa = 0（滚子是圆柱形）
• fai = 0（滚道是圆柱面）

4.2 输出文件（TAPER.txt）

  1269.  1092.  1061.  1041.  1028.  1018.  1011.  1005.  1001.   998.
   995.   993.   991.   990.   990.   990.   990.   991.   993.   995.
   998.  1001.  1005.  1011.  1018.  1028.  1041.  1061.  1092.  1269.

30 个数值，对应 30 个切片的接触应力（MPa）。

可以看到：

• 两端应力高（边缘效应）：1269 MPa
• 中间应力低且均匀：约 990 MPa
• 对数凸度修形有效减缓了边缘应力集中

4.3 应力分布可视化

计算完成后生成应力分布图：

• inner.png：内滚道接触应力分布
• outer.png：外滚道接触应力分布

def make_report(self):
    """生成报告和图表"""
    # ... 计算 ...
    
    # 绘制内滚道应力分布
    plt.figure(figsize=(6, 4), dpi=200)
    for angle in df_inner.index:
        plt.plot(
            df_inner.loc[angle, :].keys(),  # X: 切片位置
            df_inner.loc[angle, :].values,   # Y: 应力值
            label=f'{angle}°'
        )
    plt.xlabel("沿滚子的距离(mm)")
    plt.ylabel("接触应力(MPa)")
    plt.title("内滚道接触应力")
    plt.savefig('inner.png')

---

5. 使用示例

5.1 完整计算流程

from CRB_main import CRB

# 创建圆柱滚子轴承对象
crb = CRB(
    z=14,           # 14 个滚子
    Pd=0.041,       # 径向游隙 0.041mm
    l=9.6,          # 滚子有效长度 9.6mm
    l_total=10,     # 滚子总长度 10mm
    fr=4450,        # 径向载荷 4450N
    r1=5,           # 滚子半径 5mm
    r2=55,          # 内滚道半径 55mm
    r3=75,          # 外滚道半径 75mm
    curveQ=14100    # 设计载荷 14100N
)

# 执行计算并生成报告
crb.make_report()

5.2 参数调整

如需修改切片数，编辑 CRB_main.py 中的 n 参数：

# 默认 30 切片
roller_inner = TAPER(..., n=30, ...)

# 更精细的计算（如 50 切片）
roller_inner = TAPER(..., n=50, ...)

---

6. 与寿命计算的关系

切片计算得到的应力分布可用于更精确的寿命预测：

1. 传统方法：使用平均应力计算 L10 寿命
2. 切片方法：考虑应力分布，基于最大应力或累积损伤计算寿命

$$L_{10} = \left(\frac{Q_c}{Q_{eq}}\right)^{10/3}$$

其中 $Q_{eq}$ 可以基于切片应力分布计算等效载荷。

---

7. 边缘效应与凸度修形

7.1 边缘效应 (Edge Loading)

对于理想的圆柱滚子（无凸度修形），如果滚子和滚道都是完美的直线贝母面，理论上应力应该是均匀分布的。但实际上会出现边缘应力集中（边缘应力趋于无穷大）：

无凸度修形的应力分布：

应力
  ↑        ∞           ∞
  |       /|           |\
  |      / |           | \
  |     /  |___________|  \
  |    /                   \
  +-----------------------------→ 滚子长度
       边缘   中间均匀   边缘

这是因为在滚子端部，接触区域突然终止，导致应力急剧上升。

7.2 凸度修形的作用

为了减少边缘应力集中，采用凸度修形（Crowning）：

凸度类型	m1	特点
直线	0	无修形，边缘应力高
圆弧	1	简单修形，需要指定圆弧半径
对数	2	最优修形，理论上可实现均匀应力分布

对数凸度的设计原理：

$$z = t_1 \cdot \ln\left(\frac{1}{1-(2y/l)^2}\right)$$

其中：

• $z$：滚子端部的倒角量（凸度修形量），单位 mm
• $y$：沿滚子轴向的位置（以滚子中心为原点），单位 mm
• $l$：滚子有效长度（即 rl），单位 mm
• $t_1$：对数凸度幅度参数，单位 mm

t₁ 与 curveQ 的计算公式

核心公式（源自 TAPER_A.FOR 第 58 行）：

$$t_1 = \frac{1}{2} \cdot t_0 \cdot \frac{curveQ}{rl}$$

其中 $t_0$ 是材料弹性常数：

$$t_0 = \frac{4(1-\nu^2)}{\pi \cdot E}$$

参数	含义	典型值	单位
$t_1$	对数凸度幅度参数	计算结果	mm
$t_0$	材料弹性常数	5.626×10⁻⁶	mm²/N
$\nu$	泊松比	0.3 (钢)	无
$E$	弹性模量	2.06×10⁵	MPa
$curveQ$	设计载荷	用户输入	N
$rl$	滚子有效长度	用户输入	mm

数值计算示例：

对于钢材（ν=0.3, E=206000 MPa）：

$$t_0 = \frac{4 \times (1-0.3^2)}{3.14159 \times 206000} = \frac{4 \times 0.91}{647118} \approx 5.626 \times 10^{-6} \text{ mm}^2/\text{N}$$

简化形式：

$$t_1 = 2.813 \times 10^{-6} \cdot \frac{curveQ}{rl} \quad \text{(mm)}$$

验证示例：

假设 curveQ = 14100 N, rl = 9.6 mm：

$$t_1 = 2.813 \times 10^{-6} \times \frac{14100}{9.6} = 4.13 \times 10^{-3} \text{ mm} = 4.13 \text{ μm}$$

对数凸度曲线的特性

位置	y 值	z(y) 值	说明
中心	0	0	无修形
1/4 处	±rl/4	0.134×t₁	微小修形
端部附近	±0.4995rl	6.2×t₁	快速增大
极端端部	±0.49995rl	8.5×t₁	趋于无穷

当实际载荷 q = curveQ 时，应力分布均匀。

7.3 设计载荷 curveQ 的作用

7.3.1 curveQ 的物理含义

curveQ 是「用于设计对数凸度修形的基准载荷」。它决定了滚子端部倒角（凸度修形量）的大小。

设计目标：对数凸度的目的是在特定载荷下，让接触应力沿滚子长度均匀分布。

无凸度修形（理想情况）：     无凸度修形（实际情况）：
应力                          应力    ∞           ∞
  ↑  ________________              ↑   /|           |\
  |  |              |             |  / |           | \
  |  |              |             | /  |___________|  \
  +------------------→            +--------------------→
     滚子长度                       滚子长度
   (理论均匀)                  (边缘应力集中)

对数凸度修形的作用就是补偿这个边缘应力集中现象。

7.3.2 curveQ 与实际载荷 q 的关系

关系	应力分布特征	说明
q = curveQ	均匀分布	最优状态 - 凸度完美补偿边缘效应
q < curveQ	中间高、边缘低	凸度过大，边缘可能脱离接触
q > curveQ	边缘高、中间低	凸度不足，边缘应力集中未完全消除

7.3.3 工程应用建议

在实际设计中：

1. 如果已知工作载荷范围：curveQ 通常选择为预期最常用载荷
2. 保守设计：curveQ 可选择较大值，确保在高载荷时边缘不会应力过高
3. 轴承样本：例如示例中 curveQ=14100N，意味着该滚子的对数凸度是按 14100N 载荷设计的

7.3.4 一句话总结

curveQ 决定了「在哪个载荷下应力分布最均匀」。

7.4 滚道直径 dr1 的符号约定

在 Hertz 接触理论中，采用符号约定区分凸面和凹面：

• 正曲率半径 → 凸面（如滚子）
• 负曲率半径 → 凹面（如外圈滚道）

因此，INPUT_TAPER.dat 中 dr1=-150.0 表示：

• 这是一个凹面滚道（外圈）
• 直径的绝对值是 150mm
• 曲率半径 aa = 0.5 × (-150) = -75mm

等效曲率半径计算（外圈接触 m2=2）：

$$r = \frac{a \cdot aa}{aa - a} = \frac{5 \times (-75)}{(-75) - 5} = \frac{-375}{-80} = 4.6875 \text{ mm}$$

---

8. 滚子倾斜角 tilt 详解

8.1 tilt 的物理含义

滚子倾斜角 tilt 表示滚子轴线与滚道轴线之间的夹角（单位：度）。

滚子倾斜示意图：

                 正常状态 (tilt=0)          倾斜状态 (tilt>0)
               ┌─────────────────┐       ┌─────────────────┐
               │     roller      │       │     roller   ╱  │
               │  ═══════════    │       │  ═══════════╱   │
               │                 │       │            ╱    │
               └─────────────────┘       └───────────╱─────┘
                                                   tilt角

8.2 tilt 对应力分布的影响

8.2.1 无倾斜 (tilt=0) 且无凸度修形

即使没有倾斜，由于边缘效应，应力分布仍然不均匀：

应力分布（tilt=0, m1=0）：

应力
  ↑        ∞           ∞      ← 边缘应力理论上趋于无穷大
  |       /|           |\
  |      / |           | \
  |     /  |___________|  \
  |    /                   \
  +-----------------------------→ 滚子长度
       边缘   中间均匀   边缘

原因：影响系数矩阵的边界特性导致边缘切片应力集中。

8.2.2 有倾斜 (tilt>0)

倾斜会导致不对称的应力分布：

应力分布（tilt>0, m1=0）：

应力
  ↑   
  |    ∞
  |   /|
  |  / |
  | /  |
  |/   |_______________        ← 大端可能不接触（应力=0）
  +----------------------------→ 滚子长度
   小端           大端
   （载荷集中）    （载荷很小/无）

8.2.3 数学分析

在 gaps 子程序中，各切片的间隙计算为：

! 对于直线滚子 (m1=0):
z = (y + rl/2.) * tilt
s(i) = dc - z

位置	y	z	s(i) = dc - z
小端	-rl/2	0	dc（正常接触）
中间	0	rl/2 × tilt	dc - rl/2 × tilt
大端	+rl/2	rl × tilt	dc - rl × tilt（可能不接触）

当 rl × tilt > dc 时，大端将脱离接触。

8.3 tilt 的来源

滚子倾斜的主要原因：

来源	说明
轴挠曲	轴在载荷下弯曲，导致内圈相对外圈倾斜
安装误差	内外圈安装不同心或不平行
壳体变形	载荷导致的壳体变形
轴承游隙	大游隙可能导致运行中产生倾斜

8.4 如何获取 tilt

8.4.1 方法一：简化轴挠曲计算（推荐）

基于材料力学的简化模型，计算轴在载荷下的挠曲角：

    F1         F2          F3
    ↓          ↓           ↓
────○──────────○───────────○────
   ▲                       ▲
  轴承1                    轴承2
   (θ₁)                    (θ₂)

挠曲角 θ 即为该处的 tilt 值。

输入参数：

• 轴的几何：各段直径、长度
• 轴承位置
• 载荷位置和大小
• 材料：弹性模量 E

8.4.2 方法二：经验估计

根据 ISO/TS 16281 和 Harris 的建议：

工况	tilt 典型范围
高精度安装	0.5' ~ 2' (0.008° ~ 0.033°)
普通安装	2' ~ 5' (0.033° ~ 0.083°)
载荷导致倾斜	取决于轴刚度

注：1' = 1/60°

8.4.3 方法三：有限元分析

使用 ANSYS/Abaqus 等软件进行完整的轴-轴承系统分析，可获得最精确的 tilt 值。

8.5 tilt 在当前代码中的使用

模块	是否使用 tilt
C 模块 (gl_crb_ext)	❌ 假设 tilt=0
Fortran 模块 (taper_slice)	✅ 作为输入参数
CRB_main.py (新版)	❌ 未调用切片计算
TAPER_FORTRAN.py (旧版)	✅ 需手动输入

8.6 INPUT_TAPER.dat 中的 tilt 位置

10.0      10.0      9.6       0.0       1340.859878 
0.0       2         30        2         -150.0    0.0       14100.0   
↑
tilt=0.0 (第二行第一个参数，单位：度)

---

9. f2py 封装（Python 直接调用）

9.1 概述

为了避免通过文件传递参数，使用 f2py 将 Fortran 代码封装为 Python 可直接调用的模块。

9.2 文件位置

文件	说明
fortran_wrapper/taper_slice.f90	Fortran 90 源代码
fortran_wrapper/_taper_slice.cpython-312-darwin.so	编译后的 Python 扩展

9.3 编译方法

cd fortran_wrapper
uv run python -m numpy.f2py -c taper_slice.f90 -m _taper_slice

9.4 Python 调用示例

import sys
sys.path.insert(0, 'fortran_wrapper')
import _taper_slice

# 计算切片应力
result = _taper_slice.taper_calculate(
    d1=10.0,        # 滚子小端直径 (mm)
    d2=10.0,        # 滚子大端直径 (mm)
    rl=9.6,         # 滚子有效长度 (mm)
    alfa=0.0,       # 滚子半圆锥角 (deg)
    q=1340.86,      # 滚子载荷 (N) - 来自载荷分布计算
    tilt=0.0,       # 滚子倾斜角 (deg)
    m1=2,           # 凸度类型: 0=直线, 1=圆弧, 2=对数
    n_slice=30,     # 切片数
    m2=2,           # 滚道类型: 0=平面, 1=内圈, 2=外圈
    dr1=-150.0,     # 滚道直径 (mm), 负值表示凹面
    fai=0.0,        # 滚道半圆锥角 (deg)
    curveq=14100.0  # 设计载荷 (N)
)

# 解析结果
p, a, q_slice, dc, q1, qm, ierr = result

print(f"错误码: {ierr}")
print(f"滚子变形: {dc:.6f} mm")
print(f"各切片接触力 (N): {[f'{x:.2f}' for x in q_slice]}")
print(f"平衡载荷: {q1:.2f} N")
print(f"各切片应力 (MPa): {[f'{x:.0f}' for x in p]}")

9.5 返回值说明

返回值	类型	说明
p	array(n)	各切片接触中心应力 (MPa)
a	array(n)	各切片接触半宽度 (mm)
q_slice	array(n)	各切片接触力 (N)
dc	float	滚子变形 (mm)
q1	float	平衡载荷 (N)
qm	float	平衡力矩 (N·mm)
ierr	int	错误码: 0=成功, 1=方程求解失败

---

10. 完整计算流程

10.1 两阶段计算

第一阶段: C 模块 (gl_crb_ext) - 载荷分布计算
━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━

输入:
  - z=14 (滚子数)
  - Fr=4450N (径向载荷)
  - Pd, l, r1, r2 (几何参数)

         ↓ Newton 迭代

输出:
  {0°: 1340.86N, 25.7°: 1128.3N, 51.4°: 689.2N, ...}
  └─ 每个滚子的接触载荷


第二阶段: Fortran 模块 (taper_slice) - 切片应力计算
━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━

对每个承载滚子:

  输入:
    - q = 1340.86N (该滚子的载荷)
    - 滚子几何参数
    - n=30 (切片数)

         ↓ 影响系数法迭代

  输出:
    30 个切片的接触应力 (MPa)

10.2 代码示例

import sys
sys.path.insert(0, 'fortran_wrapper')
import gl_crb_ext
import _taper_slice

# 第一阶段: 载荷分布计算
load_dist = gl_crb_ext.calculate_load_distribution(
    z=14, pd=0.041, l=9.6, fr=4450, r1=5, r2=55
)

# 第二阶段: 对每个承载滚子计算切片应力
for angle, (deformation, load) in load_dist.items():
    if load > 0:  # 只计算承载滚子
        result = _taper_slice.taper_calculate(
            d1=10.0, d2=10.0, rl=9.6, alfa=0.0,
            q=load,  # 使用该滚子的实际载荷
            tilt=0.0, m1=2, n_slice=30, m2=2,
            dr1=-150.0, fai=0.0, curveq=14100.0
        )
        p, a, q_slice, dc, q1, qm, ierr = result
        print(f"{angle}° 滚子: 载荷={load:.1f}N, 最大应力={max(p):.0f}MPa, 切片力总和={sum(q_slice):.1f}N")

---

参考文献

标准与规范

1. ISO 281:2007 - Rolling bearings — Dynamic load ratings and rating life
2. ISO/TS 16281:2008 - Rolling bearings — Methods for calculating the modified reference rating life

核心理论著作

3. Harris, T.A., Kotzalas, M.N. - Rolling Bearing Analysis, 5th Edition, CRC Press, 2006

   • 第6章：Contact Stress and Deformation（接触应力与变形）
   • 第8章：Roller Bearing Loads（滚子轴承载荷）
   • 影响系数法的核心参考，详细介绍了弹性接触分析方法

4. Johnson, K.L. - Contact Mechanics, Cambridge University Press, 1985

   • 第4章：Normal Contact of Elastic Solids - Line Loading（弹性固体法向接触 - 线载荷）
   • 第9章：Rolling Contact（滚动接触）
   • 弹性半空间理论的经典著作，TAPER_A.FOR 影响系数矩阵算法的理论基础

5. Lundberg, G. & Palmgren, A. - Dynamic Capacity of Rolling Bearings, 1947

   • 滚动接触疲劳寿命理论的奠基之作
   • 载荷-寿命指数 p=10/3（线接触滚子轴承）的理论来源

中文参考

6. 《滚动轴承分析计算与应用》- 邓四二，机械工业出版社

TAPER_A.FOR 代码来源

TAPER_A.FOR 源代码（Version 8.0, Jun. 4, 2007）由 F.K. Wu（吴飞科） 编写。

该程序采用的影响系数法（Influence Coefficient Method）基于 Johnson 弹性接触力学理论，通过以下步骤求解：

1. 将滚子离散为 n 个条带（切片）
2. 基于弹性半空间假设，建立条带间的弹性耦合关系（影响系数矩阵）
3. 使用高斯消元法求解线性方程组
4. 迭代收敛得到各条带的接触应力分布

这种方法比简化的 Hertz 公式更精确，能够处理：

• 边缘效应（edge loading）
• 非均匀应力分布
• 不同凸度修形（直线、圆弧、对数）
• 接触区域的动态判断