企业官网建设流程全解析

做齿轮动力学的人，大概率都被“时变啮合刚度”这个词折腾过。它是齿轮系统振动噪声分析里最核心的激励源，算不准，后面的动态响应基本就是空中楼阁。工程上算时变啮合刚度，方法无非那么几类：有限元精度高但建模麻烦、计算慢；而石川公式法作为经典解析方法，用悬臂梁模型把轮齿变形拆解成几部分，配合MATLAB可以在几秒钟内算出一条刚度曲线，非常适合方案设计阶段的参数扫描和规律研究。这篇文章就从一个完整项目出发，把石川公式的原理、MATLAB实现步骤、算例结果和实际踩过的坑全部梳理一遍，给正在做齿轮设计或者科研入门的同学一个可以直接照抄的参考。

1. 项目背景与整体设计思路

1.1 时变啮合刚度在齿轮系统分析中的位置

齿轮正常啮合时，由于重合度通常大于1，啮合区内会出现单齿承载和双齿承载交替的状态。双齿啮合时载荷由两对齿分担，整体刚度高；单齿啮合时只有一对齿受力，刚度明显下降。这个刚度变化随啮合相位周期性出现，就是“时变啮合刚度”。它本质上是一个周期性的参数激励，直接决定了齿轮系统的动载荷、振动响应和辐射噪声。

在做齿轮箱故障诊断或者减振降噪设计时，很多人上来就直接建动力学模型、算频谱，但忽略了激励源本身的准确性。如果时变啮合刚度的幅值和相位不对，后面算出来的边频带、共振峰值全都会偏掉。所以一套能够快速、可重复地计算时变啮合刚度的工具，是齿轮动力学研究的基础设施。

时变啮合刚度计算的价值还不止于动力学。齿面修形量设计、齿根弯曲强度校核、胶合与磨损分析，很多环节都要用到啮合刚度这个基础量。能够快速算出不同参数组合下的刚度曲线，工程师就能在设计阶段把振动风险提前规避掉，而不是等样机出来再改。

1.2 为什么选择石川公式法而不是有限元

我刚接触这个课题时，第一反应也是“直接上有限元不就行了”。齿轮啮合刚度用有限元确实可以算得很准，二维平面应力模型加接触单元，网格细一点，结果跟试验对得上。但它的问题也很明显：每对啮合位置都要重新建立接触，一对齿轮在整个啮合周期里要算几十个工况，每个工况做一次非线性接触求解，时间成本极高。做参数研究的时候，模数、齿数、变位系数稍微一改就要重新建模，效率完全跟不上。

石川公式法的思路完全不一样。它把单个轮齿看作受载的变截面悬臂梁，通过材料力学公式直接计算齿面载荷作用点处的变形，再换算成刚度。整个过程是纯解析的，没有网格、没有接触迭代，计算速度以毫秒计。虽然精度和有限元相比有一定差距，但在工程初步设计和规律性研究层面完全够用。更重要的是，解析公式把每个几何参数都显式地暴露出来，参数对刚度的影响规律一目了然，这是有限元黑盒计算很难做到的。

石川公式也不是没有缺点。它对齿根过渡曲线、齿圈弹性做了大量简化，算出的绝对刚度值可能跟实测差10%甚至更多。但它的相对变化趋势、单双齿交替的形态、参数敏感度规律，跟试验和有限元都是一致的。作为工程快速评估工具，这个性价比很划算。

1.3 MATLAB实现的技术路线

整个项目的实现路线我分成四步走：

第一步，定义齿轮的基本几何参数，包括模数、齿数、压力角、变位系数、齿宽、弹性模量和泊松比，计算基圆、齿顶圆、齿根圆、分度圆齿厚这些基础量。

第二步，基于石川公式的悬臂梁模型，写一个单齿柔度计算函数。输入某个接触点的位置和载荷角，输出该齿轮轮齿在法向载荷下的柔度，也就是单位力产生的变形量。这里要处理弯曲、剪切、基础变形三部分，接触变形单独再算。

第三步，计算啮合几何关系，包括实际啮合线长度、基节、重合度。有了这些量，就能把一个啮合周期内每个时刻参与啮合的齿对数量和啮合位置确定下来。

第四步，在MATLAB中对一个啮合周期做离散扫描，对每个啮合相位计算参与啮合的所有齿对的柔度并求和取倒数，得到该时刻的总啮合刚度，最后画成曲线。

这套流程从公式到代码是逐行对应的，每一步都有明确的物理含义，调试起来也很方便。下面我把每一步的原理和关键细节展开讲。

2. 石川公式法的原理与公式体系

2.1 悬臂梁模型怎么建

石川公式法的核心思想，是把轮齿简化成固定在齿圈上的悬臂梁。玩过跳水板的人都懂，人在板头一踩，板子会弯，弯多少取决于板子伸出多长、截面多厚、材质多硬。轮齿受力弯曲是一个道理，载荷作用点越高、齿根越薄、材料越软，变形越大。

但轮齿不是等截面的规则长方体，从齿根到齿顶厚度是逐步变小的。石川公式的处理办法很朴素：把轮齿沿高度方向切成很多层，每一层近似成一个矩形截面，宽度的变化按线性处理，这样整个齿就变成一个梯形截面的悬臂梁组合。这个梯形假设是石川公式的“灵魂”，它让弯曲变形的积分变得非常直观。

实际操作中，我并没有完全照搬石川论文里那几个复杂的闭合表达式，而是按变截面悬臂梁的通用思路，在MATLAB里把轮齿沿齿高离散成几十层，逐层计算弯曲变形和剪切变形，最后累加。这样做的好处是代码逻辑清晰，也能方便地处理不同齿廓形状，跟很多文献中“离散化石川公式”的处理方式是一致的。如果你手头有齿轮手册里石川公式的标准系数，可以在我的代码基础上直接替换公式，不影响整体框架。

2.2 轮齿变形的四个组成

一对齿在啮合点受到法向载荷作用时，总变形量不是单一的弯曲造成的，拆开来看至少有四部分：

弯曲变形。这是最主要的部分，就是悬臂梁在载荷作用下的挠曲。计算时把轮齿离散成n层，第i层的高度位置已知，该层截面对应的齿厚由梯形假设线性插值得到，该层中心到载荷点的力臂也容易求。然后按材料力学中悬臂梁受横向力作用的挠曲公式积分累加，得到总的弯曲柔度。

剪切变形。齿轮轮齿的截面相对较“矮胖”，剪切变形不能忽略。这部分由切应力引起，跟载荷在齿面法向上的分力和截面积有关。计算时同样按层次累加，注意剪切修正系数，通常取1.2左右。

齿根基础变形。齿根并不是固定在绝对刚性的基础上，齿圈本身在载荷下也会发生局部弹性变形，导致齿根产生附加位移。这部分在石川公式中常用经验公式处理。我在代码中采用的是一组基于工程经验的系数表达式，用齿根处截面尺寸和载荷作用高度归一化处理。不同文献给出的基础变形系数略有差异，实际使用时可以结合有限元结果或者试验数据标定。

接触变形。齿面接触区域在载荷下会产生赫兹局部压陷，两齿面相互趋近，这也算进啮合总变形里。计算时把接触点附近的两个齿面近似成半径不同的两平行圆柱体，用赫兹线接触公式求解接触半宽和趋近量。

把四部分柔度加在一起，就是一对齿在该啮合位置的总柔度，取倒数就是这对齿的啮合刚度。

2.3 啮合位置几何与多齿啮合合成

一对齿从进入啮合到退出啮合，接触点在齿面上是不断移动的。在啮合线这个参考系里看，接触点从啮合线的一端匀速移动到另一端。这个移动距离叫作实际啮合线长度，它与基节的比值就是重合度。

重合度大于1意味着什么？在某段时间里，前一齿对还没退出，后一齿对已经进入，两对齿同时分担载荷。这时候整体啮合刚度不是某对齿单独算出来的刚度，而是并联关系，总柔度等于各齿对柔度之和，总刚度等于总柔度的倒数。在单齿啮合区，只有一对齿工作，刚度自然就低。这个“高-低-高-低”的周期性波动，正是齿轮振动的主要激励。

判断某个时刻有几对齿参与啮合，方法很巧妙。把当前刚进入啮合的齿对在啮合线上的位置记作x，它从0逐渐增加到基节长度pb，完成一个啮合周期。与此同时，它的“前任”齿对位置是x + pb，如果这个位置仍然小于实际啮合线长度，说明前任还没退出，还有两对齿。等x + pb超过实际啮合线长度，就只剩当前这一对齿了。沿一个周期扫描x，就能得到完整的时变啮合刚度曲线。

这里还有一个几何细节需要处理。接触点在啮合线上移动时，它在两个齿轮齿廓上的位置是反方向变化的：在齿轮1上从齿根滑向齿顶，在齿轮2上从齿顶滑向齿根。因此，对齿轮1和齿轮2，接触点处的接触半径、压力角都要分别由啮合位置x和(L_act - x)推导，不能混用。

3. MATLAB程序实现细节

3.1 齿轮几何参数模块

我习惯把所有齿轮几何参数封装成结构体，这样函数传参非常方便。下面这段代码定义了齿轮的基本参数并计算出各特征圆半径和齿厚函数：

%% 齿轮参数定义 g1.m = 2; % 模数 mm g1.z = 20; % 齿数 g1.alpha = 20*pi/180; % 分度圆压力角 rad g1.ha = 1; % 齿顶高系数 g1.c = 0.25; % 顶隙系数 g1.x = 0; % 变位系数 g1.b = 20; % 齿宽 mm g1.E = 2.06e5; % 弹性模量 MPa g1.nu = 0.3; % 泊松比 g2 = g1; g2.z = 40; % 从动轮齿数 %% 几何预处理 function [g1, g2, geo] = gear_geometry(g1, g2) % 分度圆、基圆、齿顶圆、齿根圆半径 g1.r = g1.m * g1.z / 2; g1.rb = g1.r * cos(g1.alpha); g1.ra = g1.r + g1.m * (g1.ha + g1.x); g1.rf = g1.r - g1.m * (g1.ha + g1.c - g1.x); g2.r = g2.m * g2.z / 2; g2.rb = g2.r * cos(g2.alpha); g2.ra = g2.r + g2.m * (g2.ha + g2.x); g2.rf = g2.r - g2.m * (g2.ha + g2.c - g2.x); % 分度圆弧齿厚 g1.st = g1.m * pi / 2 + 2 * g1.x * g1.m * tan(g1.alpha); g2.st = g2.m * pi / 2 + 2 * g2.x * g2.m * tan(g2.alpha); % 中心距（标准安装） geo.a = g1.r + g2.r; % 啮合角（标准安装时等于分度圆压力角） geo.alpha_w = g1.alpha; % 实际啮合线长度 geo.L_act = 0.5 * (sqrt(g1.ra^2 - g1.rb^2) + sqrt(g2.ra^2 - g2.rb^2)) ... - geo.a * sin(geo.alpha_w); % 基节 geo.pb = pi * g1.m * cos(g1.alpha); % 重合度 geo.eps = geo.L_act / geo.pb; end

代码里最需要注意的是单位统一。我全程用毫米和牛顿，弹性模量单位是MPa，也就是N/mm²，这样最后算出来的刚度单位是N/mm。如果混用了米和帕斯卡，结果会差好几个数量级。

任意半径处的齿厚计算要用到渐开线函数invα = tanα - α。这个函数在MATLAB里没有内置，需要自己写。齿厚公式是：

function s = gear_st(g, r) inv = @(a) tan(a) - a; a_r = acos(min(g.rb ./ r, 1)); % 该半径处的压力角 s = 2 * r .* (g.st / (2 * g.r) + inv(g.alpha) - inv(a_r)); end

这里有一个要注意的细节：当r小于基圆半径时，acos函数会出问题，因为g.rb/r大于1。我在代码里用了min(g.rb./r, 1)做截断，避免复数结果。但这也意味着基圆以下的齿厚计算带有近似性质，后面我会专门讲这个坑。

3.2 单齿柔度计算函数

单齿柔度是整个程序的核心，输入该齿轮的接触点半径和压力角，输出该齿轮单位法向载荷产生的变形量，也就是柔度。弯曲和剪切部分用分层累加实现，齿根基础变形用经验公式：

function c = tooth_compliance(g, r_c, alpha_c) E = g.E; nu = g.nu; b = g.b; h_x = r_c - g.rf; % 接触点距齿根的高度 if h_x < 1e-6 c = 0; return; end n_layer = 50; % 沿齿高离散层数 dh = h_x / n_layer; S_f = gear_st(g, g.rf); % 齿根处齿厚 S_x = gear_st(g, r_c); % 接触点处齿厚 G = E / (2 * (1 + nu)); % 剪切模量 sum_b = 0; % 弯曲累加项 sum_s = 0; % 剪切累加项 for i = 1:n_layer y = (i - 0.5) * dh; % 第i层距齿根高度 S_i = S_f + (S_x - S_f) * y / h_x; % 梯形截面假设 L_i = h_x - y; % 该层中心到载荷点的力臂 I_i = b * S_i^3 / 12; % 矩形截面惯性矩 sum_b = sum_b + L_i^2 / I_i * dh; sum_s = sum_s + 1 / (b * S_i) * dh; end % 弯曲柔度 delta_b_over_F = 12 * cos(alpha_c)^2 / (b * E) * sum_b; % 剪切柔度 delta_s_over_F = 1.2 * sin(alpha_c)^2 / (G * b) * sum_s; % 齿根基础变形（经验公式） delta_f_over_F = cos(alpha_c)^2 / (b * E) * 3.5 * (h_x / S_f); c = delta_b_over_F + delta_s_over_F + delta_f_over_F; % mm/N end

离散层数取50层，这个值是我试出来的平衡点。层数太少，结果对网格密度敏感，曲线不光滑；层数太多，计算时间虽然也不长，但没必要。如果你想让曲线更平滑，可以取100层。

齿根基础变形这个经验式，不同文献差异真的很大。我用的3.5这个系数是一个参考文献里的推荐值，代码里你要根据自己验证的结果去标定。如果你的计算对象是齿圈很薄的齿轮，这个系数可能需要加大。

3.3 时变啮合刚度主程序

主程序的逻辑很简单：在一个基节周期内均匀取若干个相位点，对每个相位点判断有几对齿在啮合，分别算柔度，合成总刚度。接触变形单独用一个函数算：

function ch = contact_compliance(g1, g2, r1c, r2c) % 渐开线齿面法向曲率半径 = 基圆切线长度 R1 = sqrt(r1c^2 - g1.rb^2); R2 = sqrt(r2c^2 - g2.rb^2); Rstar = (1/R1 + 1/R2)^(-1); % 等效曲率半径 F = 100; % 假定的法向载荷，用于接触半宽 w = F / g1.b; % 线载荷 N/mm Estar = 1 / ((1 - g1.nu^2)/g1.E + (1 - g2.nu^2)/g2.E); aH = sqrt(4 * w * Rstar / (pi * Estar)); % 赫兹接触半宽 % 两平行圆柱接触趋近量换算的柔度 ch = 2 / (pi * Estar * g1.b) * (log(4 * Rstar / aH) - 1); % mm/N end

赫兹接触变形对总刚度的影响相对较小，但它的量级要大致对。注意这里假定了100N的参考载荷，因为赫兹接触有轻微的非线性，刚度跟载荷的对数相关，影响不大，但要在代码里说明清楚。

单齿对刚度的组装函数：

function K = single_pair_stiffness(g1, g2, geo, x) % x: 该齿对在啮合线上距A点的位置 (0~L_act) N1A = sqrt(g2.ra^2 - g2.rb^2); % 齿轮1基圆切点N1到啮入点A的距离 N2B = sqrt(g1.ra^2 - g1.rb^2); % 齿轮2基圆切点N2到啮出点B的距离 r1c = sqrt(g1.rb^2 + (N1A + x)^2); % 齿轮1接触点半径 r2c = sqrt(g2.rb^2 + (N2B - x)^2); % 齿轮2接触点半径 alpha1 = acos(g1.rb / r1c); % 齿轮1接触点压力角 alpha2 = acos(g2.rb / r2c); c1 = tooth_compliance(g1, r1c, alpha1); c2 = tooth_compliance(g2, r2c, alpha2); ch = contact_compliance(g1, g2, r1c, r2c); K = 1 / (c1 + c2 + ch); % N/mm end

这里要注意r1c和r2c的计算。N1A是从动轮齿顶圆切点长度，它决定了齿轮1在啮入点A时的接触半径；N2B是主动轮齿顶圆切点长度，决定了齿轮2在啮出点B时的接触半径。x从0走到L_act，齿轮1的接触半径先小后大，对应齿根到齿顶；齿轮2反过来，先大后小，对应齿顶到齿根。逻辑上是自洽的。

最后的扫描循环：

%% 时变啮合刚度扫描 n_samp = 200; phi = linspace(0, 1, n_samp); K_total = zeros(n_samp, 1); K_cur = zeros(n_samp, 1); K_pre = zeros(n_samp, 1); for i = 1:n_samp x_cur = phi(i) * geo.pb; % 当前齿对在啮合线上的位置 Kc = single_pair_stiffness(g1, g2, geo, x_cur); K_cur(i) = Kc; K_total(i) = Kc; x_pre = x_cur + geo.pb; % 前一齿对的位置 if x_pre <= geo.L_act Kp = single_pair_stiffness(g1, g2, geo, x_pre); K_pre(i) = Kp; K_total(i) = K_total(i) + Kp; end end

扫描点数取200个每周期，画出来的曲线已经很光滑，计算时间也就几毫秒级别。如果需要做参数扫描，这个速度可以轻松跑上千组工况。

4. 算例结果与参数影响规律

4.1 典型刚度曲线长什么样

用上面代码跑一组标准齿轮副，模数2mm，主动轮20齿，从动轮40齿，压力角20度，齿宽20mm，得到的时变啮合刚度曲线具有非常鲜明的特征：在一个啮合周期内，双齿啮合区刚度明显高出一截，单齿啮合区刚度掉下去，形成“山峰-山谷-山峰-山谷”的周期性形态。

从机理上看，双齿区的“山峰”并不是因为每个齿变硬了，而是两个弹簧并联，总柔度是两个齿对柔度之和的倒数，必然比单个齿对的刚度高。这个高低交替的幅值差，正是齿轮啮合冲击的主要来源。在设计时，可以调节重合度来改变单双齿区的宽度比例，从而控制振动激励的强度。

另一个容易忽略的现象是，双齿啮合区内两个齿对分别承担的载荷并不相等。先进入啮合的齿对在齿根位置接触，柔度较大，分配到的载荷比例较小；后进入的齿对接触位置更靠近齿顶，柔度较小，承担的载荷反而更大。这对齿根应力分析和修形设计都有直接参考价值，因为实际齿面载荷分布并不是均等的。

新进入啮合的齿对在x=0处，理论上轮齿齿根处的齿厚很小，柔度很大，但x=0对应的是从动轮齿顶与主动轮齿根接触的位置，齿面曲率半径的几何关系让这个位置的刚度并不是全周期最低点。真正的最低点出现在单齿啮合区的中间位置，因为此时既没有第二对齿分摊载荷，当前齿对的接触点又离齿根有一段距离，弯曲力臂较大。

4.2 参数变化对刚度曲线的影响规律

用这个程序做参数扫描非常方便，我这里总结几个我实测过的规律：

模数的变化对刚度的影响最直接。模数增大，齿厚近似线性增大，弯曲柔度大幅下降，啮合刚度显著上升。但同时轮齿的尺寸变大，基节增大，啮合周期的“宽度”也变了。模数的影响基本是一次方的强正相关。

齿数增加时，在模数不变的情况下，齿轮直径变大，齿根圆的曲率变平缓，单齿柔度会下降，但重合度也会增加，双齿啮合区的比例上升。整体上刚度曲线变得更高更平缓，波动幅度减小。齿数增多对降低振动激励是有利的，这就是为什么高速齿轮箱倾向选多齿数方案的原因之一。

变位系数的影响比较微妙。正变位会让齿根厚度增加、齿顶变尖，单齿刚度在齿根附近有所提升，但齿顶附近的齿厚变小，接触点靠近齿顶时刚度反而下降。同时变位改变了实际啮合线长度和重合度。这个规律必须在具体参数下扫描才能看得清楚，不建议凭经验拍脑袋。

齿宽对刚度是线性正比关系，这个很好理解，相当于悬臂梁的宽度增加。但要注意，齿宽增大后接触线的长度也增加，接触柔度同样会下降，所以刚度跟齿宽接近线性关系而不是严格线性。在对比不同方案时，如果只关心啮合刚度的形态，可以用单位齿宽刚度来归一化比较。

5. 实操中的常见问题与排查技巧

5.1 几何计算的坑

先说单位。我在代码里统一用mm和N，但很多人从文献里抄公式时会遇到把弹性模量写成2.06e11的情况，那是帕斯卡单位，如果直接跟mm混算，结果会差1000倍。建议在最开始加一行注释，明确所有输入输出的单位，免得过了两个月回来自己都看懵。

其次是渐开线齿厚公式在基圆以下失效的问题。标准齿轮z=20时，基圆半径大于齿根圆半径，齿根圆这一段其实已经不是渐开线了。代码里用min(g.rb./r, 1)截断虽然避免了报错，但它默认了齿根圆处仍然按渐开线延伸计算齿厚，这会让齿根圆齿厚偏大。好在石川公式本身就是近似方法，这个误差在工程上可以接受。如果你需要更高的精度，建议把齿根过渡曲线用圆弧近似，再重新推导根部齿厚。

再来说说啮合点位置判断。很多初学者会在主循环里写if x_pre <= geo.L_act，却忽略了当x_cur本身就接近基节时，x_pre可能已经超过L_act很多，这时只有一个齿对。这个判断必须放在循环体内部，对每个相位点单独判断，不能在整个周期外统一判断。

还有一个隐蔽的坑是啮合角。如果齿轮是变位安装或者中心距有调整，啮合角不再等于分度圆压力角，实际啮合线长度的公式要改用啮合角。代码里我默认了标准安装，你在实际项目中要加上啮合角的反算逻辑，用渐开线函数方程fzero求解。

5.2 结果验证与误差控制

程序写完了，最重要的一步是验证。我一般分三个层次检查：

第一层，量级判断。齿轮啮合刚度的典型量级是每毫米宽度在10⁵到10⁶ N/mm之间。如果算出来只有几千或者上千万，基本就是单位或者公式系数出了问题。20mm齿宽、2模数齿轮副，总刚度在2×10⁵ N/mm量级是比较合理的。

第二层，用有限元对比。如果你手头有ANSYS或者Abaqus，选取几个特征啮合位置做二维接触分析，跟石川公式结果对比。对比时重点关注变化趋势，而不是绝对数值。石川公式通常算出的刚度比有限元偏高或者偏低，取决于基础变形系数的标定，这个偏差范围在20%以内都算正常。

第三层，跟文献数据对比。齿轮传动领域关于时变啮合刚度的论文非常多，很多都提供了具体参数下的刚度曲线。找参数一致或者接近的文献，把曲线叠在一起看形状是否吻合，特别是单双齿交替的位置有没有错位。如果错位，大概率是啮合线长度或者基节计算有误。

5.3 代码性能与工程化建议

虽然解析法已经很快，但做大批量参数优化时，循环里面的分层累加还是会成为瓶颈。我有一个优化建议：把50层离散和循环累加改成矩阵向量化计算，用数组运算替代for循环，速度能快5到10倍。做法是把齿高分层向量、齿厚向量、力臂向量都用列向量表示，然后用sum函数一次性累加。

另一个工程化建议是把整个计算封装成一个函数，输入参数和输出刚度曲线封装成结构体，方便后续和动力学模型联动。比如做齿轮箱振动响应计算时，你只需要把这个函数放在循环里，每步更新啮合相位，输出当前刚度值，就能嵌入到Runge-Kutta积分器里。

最后想说一下代码的可读性。这种分析程序写完以后，过半年你大概率要回头改参数。我给每个函数开头都写了输入输出注释，关键公式旁边标注了物理意义和单位。尤其是石川公式这种经验性很强的内容，建议把公式来源文献和版本号也写在注释里，方便追溯。

我在实际项目里用这套代码做过的最大规模一次扫描，是连续跑了5000组不同的齿轮参数组合，研究变位系数和齿数对啮合刚度波动幅值的影响，整个过程不到一分钟。这种批量分析能力，是有限元方案很难想象的。也是从那次以后，我更加确信石川公式法在工程前期设计中的价值，关键是你得把它的边界条件搞清楚，知道它哪里准、哪里近似、哪里可以用、哪里不能用。只要心里有数，这套方法就是齿轮振动分析工具箱里最顺手的一把刀。