Metadata-Version: 2.4
Name: bstart-trb
Version: 0.2.0
Summary: Add your description here
License-File: LICENSE
Requires-Python: >=3.11
Requires-Dist: numpy>=2.4.2
Requires-Dist: openpyxl>=3.1.5
Requires-Dist: pandas>=3.0.0
Requires-Dist: plotly>=6.5.2
Description-Content-Type: text/markdown

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

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

---

## 目录

- [1. 切片计算概述](#1-切片计算概述)
- [2. 载荷分布计算](#2-载荷分布计算)
- [3. 滚子切片应力计算](#3-滚子切片应力计算)
- [4. 输入输出格式](#4-输入输出格式)
- [5. 使用示例](#5-使用示例)
- [6. 与寿命计算的关系](#6-与寿命计算的关系)
- [7. 边缘效应与凸度修形](#7-边缘效应与凸度修形)
- [8. 滚子倾斜角 tilt 详解](#8-滚子倾斜角-tilt-详解)
- [9. f2py 封装（Python 直接调用）](#9-f2py-封装python-直接调用)
- [10. 完整计算流程](#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`）：

```python
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 输出示例

```python
# 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 类参数说明

```python
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）
   - 内层：接触区域收敛迭代

**关键代码片段：**

```fortran
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 代码实现（旧版 - 文件传参）

```python
# 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`：外滚道接触应力分布

```python
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 完整计算流程

```python
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` 参数：

```python
# 默认 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` 子程序中，各切片的间隙计算为：

```fortran
! 对于直线滚子 (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 编译方法

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

### 9.4 Python 调用示例

```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 代码示例

```python
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）
- 非均匀应力分布
- 不同凸度修形（直线、圆弧、对数）
- 接触区域的动态判断
