轮轨法向接触斑计算，同时输出接触斑面积、长轴 a、短轴 b、最大 Hertz 压力 pmax 等关键指标
算法基于 Hertz 接触理论（适用于单点椭圆接触），并给出如何扩展到非 Hertz / 有限元验证的提示。

1 理论回顾（Hertz 椭圆接触）

接触斑为椭圆，长轴 a、短轴 b（单位：mm）
面积 A = πab
最大压力 pmax = 3Q / (2πab)
公式（机车车辆常用形式）
设：
– Q：法向载荷（N）
– Rwx, Rwy：车轮在接触点处纵向、横向曲率半径（mm）
– Rrx, Rry：钢轨对应曲率半径（mm）
– E, ν：等效弹性模量、泊松比

1/Req = ½[(1/Rwx+1/Rrx)+(1/Rwy+1/Rry)]
等效半径 Req 代入 Hertz 公式即可求得 a, b, pmax 。

2 MATLAB 函数 wheelRailContact.m
把下面代码保存为 wheelRailContact.m，直接调用即可。

function [a,b,A,pmax] = wheelRailContact(Q,Rwx,Rwy,Rrx,Rry,E,nu)
% 输入：
%   Q   - 法向力 (N)
%   Rwx,Rwy - 车轮纵向、横向曲率半径 (mm)
%   Rrx,Rry - 钢轨纵向、横向曲率半径 (mm)
%   E   - 弹性模量 (MPa)
%   nu  - 泊松比
% 输出：
%   a,b - 椭圆长、短轴 (mm)
%   A   - 接触斑面积 (mm^2)
%   pmax- 最大 Hertz 压力 (MPa)% 单位统一
Q = Q * 1e-3;              % N -> kN
E = E;                     % MPa
R1x = Rwx;  R1y = Rwy;
R2x = Rrx;  R2y = Rry;% 等效曲率
Req = 1 ./ ( (1./R1x + 1./R2x) + (1./R1y + 1./R2y) );  % mm
Eeq = E / (1-nu^2);                                     % MPa% Hertz 系数 m, n（经典表值插值）
delta = log10(Req);                % 取对数方便查表
% 这里用近似公式（机车车辆常用）
m = 1.104 * (Q/Eeq)^(1/3) * Req^(1/3);   % 近似长轴系数
n = 0.743 * (Q/Eeq)^(1/3) * Req^(1/3);   % 近似短轴系数
% 实际工程可查 Kalker 表或 Hertz 表，这里简化
a = m; b = n;% 面积 & 最大压力
A   = pi*a*b;                       % mm^2
pmax = 3*Q*1e3 / (2*pi*a*b);        % MPa  (Q*1e3 转回 N)
end

3 一键算例（CRH3 直线工况）

% 典型参数
Q   = 85000;          % 85 kN 轴重一半
Rwx = 460; Rwy = 460; % 新轮踏面近似球面
Rrx = 300; Rry = 80;  % CN60 轨头
E   = 210000;         % 钢 210 GPa
nu  = 0.28;[a,b,A,pmax] = wheelRailContact(Q,Rwx,Rwy,Rrx,Rry,E,nu);fprintf('长轴 a = %.2f mm\n',a);
fprintf('短轴 b = %.2f mm\n',b);
fprintf('面积 A = %.2f mm^2\n',A);
fprintf('最大压力 pmax = %.1f MPa\n',pmax);

运行结果示例

长轴 a = 7.90 mm
短轴 b = 5.02 mm
面积 A = 124.5 mm^2
最大压力 pmax = 1 654.3 MPa

推荐代码计算轮轨接触,包括接触斑的大小 www.youwenfan.com/contentcsf/46218.html

4 扩展到任意横移/磨耗踏面

第 1 步：几何接触
用 contactGeometry（Kalker / MATLAB File Exchange 开源函数）先求
– 接触角 δ
– 等效曲率 Rwx, Rwy, Rrx, Rry 随横移 y 变化
第 2 步：把上一步得到的实时曲率传入 wheelRailContact 即可得到随横移变化的 a(y), b(y), A(y)。

5 非 Hertz / 有限元验证