輪軌法向接觸斑計算 ，同時輸出 接觸斑面積、長軸 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 / 有限元驗證
- 當輪緣貼靠或磨耗出現共形/多點接觸時，Hertz 假設失效，需采用
– CONTACT（Kalker 三維程序）
– ABAQUS / ANSYS 顯式有限元（彈塑性、大變形） - 在 MATLAB 中可調用 Python-Abaqus 批處理，再把接觸斑面積讀回做對比。
