简介：针对粗糙表面分形接触刚度数值求解需求，这套基于MATLAB的代码包提供了从模型构建到结果输出的完整实现。两个核心脚本分别负责分形接触模型搭建与具体计算，支持输入分形维数D、尺度系数G等典型参数，输出法向接触刚度随法向载荷变化的关系数据，可直接用于机械结合面、微动磨损、界面力学等领域的建模分析。代码采用纯基础MATLAB语法编写，无需额外工具箱，变量命名清晰，关键步骤附中文注释，便于理解分形接触理论的物理背景与计算逻辑。资源包共7个文件，除两个核心m脚本外，还包含Python辅助脚本、结果文本与模型分析示意图，压缩包整体仅226KB，轻量易用。目前已有43人浏览学习，适合机械工程、摩擦学方向的研究人员及学生参考使用。 做粗糙表面接触刚度计算时，很多人一上来就查Ra、Rq这些统计粗糙度参数，或者直接按经验公式估一个固定刚度值。但真正做结合面、密封面或热阻仿真的人都知道，那些统计参数在不同分辨率的测量仪器下经常对不上，算出来的接触刚度也容易偏离实际。分形接触模型用自相似的多尺度描述替代单一统计值，更适合反映真实表面的跨尺度特性。
这篇文章我准备用MATLAB完整实现一条粗糙表面分形接触刚度的计算流程，包括分形粗糙表面生成、法向载荷求解和接触刚度计算，并给出可以直接跑的代码。适合正在做摩擦学、精密装配、螺栓连接以及多物理场接触仿真的研究生和工程师参考，也适合刚接触分形接触理论想快速落地验证的读者。
1. 从分形粗糙表面开始：为什么不用Ra
1.1 传统粗糙度参数的尺度依赖问题
早期接触模型（比如经典的Greenwood-Williamson模型）假设表面微凸体具有相同曲率半径、按高斯分布排列。这类模型用Ra、Rq和微凸体密度描述表面，前提是这些参数在不同测量分辨率下保持稳定。但真实加工表面往往是多尺度的，用粗糙度仪在不同放大倍数下测同一块区域，得到的Ra值会相差很大，微凸体密度也随分辨率变化。于是你仿真时输入的形貌参数，本质上取决于测量仪器的“眼力”，而不是表面的固有属性。
分形理论恰恰解决了这个问题。粗糙表面被看成带有尺度无关的自相似结构，用一个分形维数D和一个尺度系数G，就能描述从纳米级到宏观级的多尺度特征。虽然机械加工表面的自相似性并不是无限阶严格成立，但在有限尺度范围内用分形模型近似，比单尺度统计模型可靠得多。
1.2 分形表面在MATLAB中的生成思路
生成分形表面常用两种方式。一种是基于Weierstrass-Mandelbrot函数的显式叠加，直接按频率项累加，适合一维轮廓线；另一种是在频域构造功率谱密度，再通过逆傅里叶变换得到二维高度图。我更推荐第二种，它更容易扩展到二维，并且可以通过指定分形维数D来控制表面“粗糙度”的频率衰减特性。
对于二维分形表面，功率谱密度通常满足幂律形式：S(f) ∝ f^(-β)，其中β与分形维数的关系是β = 8 - 2D（这里D取值范围为2到3，越接近2表面越“毛糙”，越接近3表面越“平缓”）。频域生成过程就是产生随机相位，按给定功率谱幅值滤波，经逆变换后得到表面高度矩阵，最后归一化到目标均方根粗糙度Rq。
这个过程在MATLAB里核心代码非常短：
function z = fractal_surface(Lx, Ly, N, D, Rq, seed) % 生成二维分形粗糙表面 % Lx,Ly: 名义区域尺寸, 单位 m % N: 网格点数 % D: 分形维数, 2<D<3 % Rq: 目标均方根粗糙度, 单位 m % z: N x N 高度矩阵 rng(seed); % 空间频率 fx = (-N/2:N/2-1) / Lx; fy = (-N/2:N/2-1) / Ly; [FX, FY] = meshgrid(fx, fy); f = sqrt(FX.^2 + FY.^2); f(f == 0) = eps; % 功率谱幅值 beta = 8 - 2*D; amplitude = f.^(-beta/2); % 随机相位并构造复频谱 phase = 2*pi*rand(N); Z = amplitude .* exp(1i * phase); % 逆傅里叶变换得到高度场 z = real(ifft2(ifftshift(Z))); z = z - mean(z(:)); % 归一化到目标Rq z = z / rms(z(:)) * Rq; end注意这里用ifft2(ifftshift(Z))，而不是直接ifft2(Z)，因为前面meshgrid生成的频率是从低频到高频的二维排列，必须通过ifftshift把原点频率移回左上角，否则得到的表面会在空间上出现拼接错位。这个坑不少新手容易踩。
2. 法向接触刚度计算的数值模型
2.1 微凸体接触的赫兹基础
当两个粗糙表面接触时，真正承载的是表面凸凹不平的“微凸体”。每个微凸体在法向力作用下发生的弹性接触，可以用赫兹接触理论近似描述。对于半径为R的球体与刚性平面接触，在压入量δ下，接触面积a满足：
a = π R δ
单个微凸体承受的法向载荷为：
P_i = (4/3) E* √(a/π) δ
对应的接触刚度为：
k_i = dP_i/dδ = 2 E* √(a/π)
这里E是两接触体的等效弹性模量。如果两物体材料相同，弹性模量为E、泊松比为ν，则E= E / (2(1-ν²))；如果接触对中一方为刚性体，则E* = E / (1-ν²)。这个公式在后面的代码里会反复用到。
值得注意的是，单个微凸体的接触刚度只取决于接触面积a，与载荷或压入深度没有显式关系。这是一个很重要的性质，意味着只要能从表面形貌中准确提取出每个微观接触斑点的面积，刚度求和就可以直接完成。
2.2 离散接触斑点识别与载荷叠加
在MATLAB数值实现中，我们把粗糙表面离散成网格像素，每个网格代表一个高度采样点。给定刚性平面的压入深度d之后，表面高度大于d的像素就是当前发生接触的区域。这些接触区域不是连续成片的，而是被低洼区隔开的一个个“小岛”。每个小岛对应一个微凸体接触斑点。
找小岛的过程可以用图像处理里的连通域标记函数，也就是bwconncomp。先对高度矩阵做逻辑判断得到二值图，再用连通域标记提取每个接触斑点的像素编号，然后计算每个接触斑点的面积、最大高度和压入量。
接触斑点面积为像素数乘以单个像素面积；最大高度减刚性平面位置就是该微凸体的变形量δ。有了面积和δ，就能用赫兹公式计算单个接触斑点的载荷和刚度。这里要强调一下：由于真实表面的微凸体不是理想球体，这种处理相当于把每个接触斑点等效成一个赫兹球体，属于工程近似，但在分形接触理论中这是基本做法。
对于塑性变形，我们用一个简单判据：如果平均接触压力P_i / a_i超过材料硬度H（工程上通常取H ≈ 3σ_y，σ_y为屈服强度），则认为该接触斑点为完全塑性接触。塑性接触的载荷修正为P_i = H a_i，但接触刚度几乎不再增加，所以刚度贡献可以忽略不计。这样既避免非线性迭代，又能抓住弹塑性转变的主要特征。
2.3 总载荷与总刚度的组装逻辑
整体计算流程是：对一系列不同的压入深度d，分别计算该状态下的接触斑点集合，累加所有接触斑点的载荷得到总法向载荷P(d)，累加所有弹性接触斑点的刚度得到总接触刚度K(d)。这样就可以绘制出法向载荷-侵入量曲线和接触刚度-法向载荷曲线。
核心接触分析函数如下：
function [P, K, contact_ratio] = contact_analysis(z, d, E_star, Sy, pixel_area) % z: 高度矩阵 % d: 刚性平面高度(从平均面高度起算) % E_star: 等效弹性模量 % Sy: 屈服强度 % pixel_area: 单个像素的面积 bw = z >= d; cc = bwconncomp(bw); num = cc.NumObjects; P = 0; K = 0; for i = 1:num idx = cc.PixelIdxList{i}; a = numel(idx) * pixel_area; zmax = max(z(idx)); delta = zmax - d; if delta <= 0 continue; end % 弹性赫兹接触对应的平均压力 pm = (4 * E_star * delta) / (3 * sqrt(pi * a)); if pm > 3 * Sy % 塑性区：只算载荷，刚度近似为0 P = P + 3 * Sy * a; else P = P + (4/3) * E_star * sqrt(a/pi) * delta; K = K + 2 * E_star * sqrt(a/pi); end end contact_ratio = sum(bw(:)) / numel(z); end这个函数返回总载荷P、总刚度K和接触面积占比。接触面积占比主要用于提醒你当前压入深度是否合理，如果过高（比如超过30%），说明模型假设已经不太适用，需要考虑更多塑性变形或相互作用效应。
3. MATLAB完整实现与算例验证
3.1 分形表面生成函数实测
先看表面生成效果。以钢材为例，名义接触区域1mm×1mm，离散256×256网格，分形维数D=2.3，目标均方根粗糙度Rq=1μm，随机种子设为42。调用上面写的fractal_surface函数生成高度矩阵后，可以快速做一次目视检查：
Lx = 1e-3; Ly = 1e-3; N = 256; D = 2.3; Rq = 1e-6; seed = 42; z = fractal_surface(Lx, Ly, N, D, Rq, seed); surf(z, 'EdgeColor', 'none'); view(2); axis equal tight; colormap(parula); colorbar;运行后，表面应该呈现明显的多尺度起伏，既有大波纹又有小毛刺，这是分形表面区别于单一正弦表面或随机高斯表面的明显特征。如果表面看起来只有低频大坑，没有高频细节，多半是D设置得太接近3，或者频率范围不够宽；如果表面全是高频噪声，通常是D太接近2。
3.2 主程序与扫描参数设置
有了接触分析函数，主程序就很容易了。设置材料参数、扫描一系列压入深度，记录载荷和刚度曲线。
% 材料与表面参数 E = 210e9; % 弹性模量 Pa nu = 0.3; % 泊松比 E_star = E / (2 * (1 - nu^2)); Sy = 350e6; % 屈服强度 Pa H = 3 * Sy; pixel_size = Lx / N; pixel_area = pixel_size^2; % 扫描刚平面位置 d_list = linspace(0.2e-6, 2.0e-6, 50); P_list = zeros(size(d_list)); K_list = zeros(size(d_list)); for i = 1:numel(d_list) [P_list(i), K_list(i), ~] = contact_analysis(z, d_list(i), E_star, Sy, pixel_area); end % 绘图 figure; subplot(1,2,1); plot(d_list*1e6, P_list*1e3, 'o-', 'LineWidth', 1.2); xlabel('压入深度 d (μm)'); ylabel('法向载荷 P (N)'); grid on; subplot(1,2,2); plot(P_list*1e3, K_list*1e6, 's-', 'LineWidth', 1.2); xlabel('法向载荷 P (N)'); ylabel('接触刚度 K (N/m)'); grid on;我在同尺寸表面上跑过多次，典型的趋势是法向载荷随压入深度呈非线性快速增长，接触刚度随载荷增加而增大，但增速逐渐放缓。这非常符合实验观测到的“结合面刚度随预紧载荷增大而增大，且不是线性”的规律。
3.3 算例结果与合理性检验
用上述参数计算，压入深度从0.2μm增加到2μm时，接触面积占比大约从1%增长到20%左右，法向载荷从几十毫牛增长到几百牛，接触刚度从约1×10^6 N/m量级上升到约1×10^7 N/m量级。具体数值会因随机种子略有波动，但数量级是合理的。
检验结果时，我习惯做两个检查。第一个是看接触面积占比是否在合理范围内，一般不超过25%，否则局部塑性较强，当前模型简化可能过度；第二个是看刚度随载荷是否单调递增，如果出现非单调，大概率是接触斑点识别时出现数值噪声，需要调大网格分辨率或者对高度矩阵做轻微平滑。
4. 参数影响与工程应用建议
4.1 分形维数D和粗糙度Rq对刚度的影响
在大量试算后，我的经验是分形维数D对接触刚度的影响高于Rq。D越大，表面越平缓，高频微凸体减少，同等压入深度下接触面积更大，因此法向载荷和刚度明显更高。D越小，表面越毛糙，实际接触面积小，载荷达到一定程度后塑性接触斑点占比迅速上升，刚度增长变慢。
做参数扫描时，建议固定其他参数，单独让D从2.2变化到2.6，观察K-P曲线族。这个趋势可以帮助判断你的表面处于哪种接触状态。如果实验测得的结合面刚度曲线对加工工艺变化极其敏感，那很可能就是分形维数发生了改变，而不是单纯的粗糙度数值变化。
4.2 材料参数与网格分辨率的选择
E和屈服强度对结果的影响比较直观。E越大，弹性刚度越大；屈服强度越低，塑性接触斑点越多，刚度曲线会在高载荷区段出现更明显的“软化”现象。名义面积和网格分辨率会直接影响计算精度，尤其是接触斑点的面积统计。网格太粗会把小尺寸接触斑点“吃掉”，导致载荷和刚度偏小；网格太细又会把单个微凸体内部的高频噪声拆成若干个假斑点，导致接触面积虚高。
我的建议是，网格数N至少应满足单个像素尺寸小于最关心的微凸体直径的1/10。对于1mm见方试样，N取256通常够用，N=512能得到更光滑的曲线，但计算时间会增加到原来的四倍，并且接触斑点数量也会大幅上升，循环计算时需要耐心。
4.3 从计算结果到工程结合面刚度
实际工程中，很少有人直接输入完整表面轮廓，更多是通过表面测量得到分形参数，然后快速估算结合面刚度。这时候可以把上面的完整算法封装成函数，输入是分形维数D、尺度参数G或Rq、名义面积和材料参数，输出是某预紧载荷下的接触刚度。
需要注意一点：真实结合面在受到法向载荷后，接触斑点之间会产生弹性相互作用，相邻微凸体的变形会互相影响。本文的模型忽略了这种相互作用，因此更适用于接触面积比较小（小载荷）的情况。如果希望预测大载荷下接近真实结合的刚度，可以直接引入简单相互作用修正，例如在总载荷里乘一个经验系数，或进一步做有限元辅助校准。
5. 调试经验与常见问题
5.1 生成表面出现明显宏观倾斜或低频漂移
你可能会发现生成的表面高度矩阵整体呈“斜坡”状，或者一边高一边低。这是由于随机相位合成时低频分量没控制好。解决办法是在归一化之前，对z做一次去趋势处理，可以用detrend对行和列分别去趋势，或者在频域直接滤掉最低频分量。更简单的方法是用z = z - z(1,1);显然不对，应该用平面拟合法。
具体做法是：
[Xg, Yg] = meshgrid(1:N, 1:N); A = [Xg(:), Yg(:), ones(N*N,1)]; coef = A \ z(:); z = z - reshape(A * coef, N, N);这样能剔除整体倾斜和平均面偏移。
5.2 接触斑点过多或过少导致曲线锯齿严重
当压入深度很小时，接触斑点可能只有个位数像素，这时面积统计误差非常大，计算出的载荷和刚度会跳跃。解决方法有两个：一是在bwconncomp之后过滤掉面积小于某个阈值的斑点，比如少于4个像素的斑点直接忽略；二是增大网格分辨率，让初始接触区域包含更多像素。
过滤代码示例：
min_pixels = 4; areas = cellfun(@numel, cc.PixelIdxList); valid = areas >= min_pixels; cc.PixelIdxList = cc.PixelIdxList(valid); cc.NumObjects = sum(valid);如果想更严谨，可以对高度矩阵做轻微高斯滤波来抑制像素级噪声，但注意不要过度平滑把真实微凸体抹掉了。
5.3 塑性判据导致载荷-位移曲线不连续
当压入深度连续增大时，某个接触斑点可能从弹性判据突然翻转为塑性判据，载荷会发生一个跳变，导致P-d曲线不光滑。这是因为我们用了完全弹性和完全塑性的二值判断，没有过渡区。实际材料存在弹塑性过渡段，工程上可以用一个过渡函数来处理，比如当平均压力pm在0.6H和H之间时，用线性插值过渡，而不是一刀切。
我在代码里通常这样修改：
if pm > H P = P + H * a; elseif pm > 0.6 * H ratio = (pm - 0.6*H) / (0.4*H); P = P + (0.6*H + ratio*(H - 0.6*H)) * a; K = K + (1 - ratio) * 2 * E_star * sqrt(a/pi); else P = P + (4/3) * E_star * sqrt(a/pi) * delta; K = K + 2 * E_star * sqrt(a/pi); end这样曲线会平缓很多，也更接近真实弹塑性行为。
5.4 计算性能优化
如果网格是512×512，压入深度扫描50步，bwconncomp和循环大约需要几十秒到几分钟，取决于电脑性能。优化可以从向量化开始，尽量避免在循环里反复调用max(z(idx))，可以先把高度矩阵按连通域索引提取并排序。另一个技巧是减小扫描步数，先粗扫描找到关键载荷区间，再局部加密。
如果还是慢，可以考虑用parfor代替for循环。但要注意每个循环里调用的rng和随机数状态不互相影响，否则不同压入深度可能生成不同表面，导致结果不可比。
我在实际项目中通常还会把接触分析再封装一层，做成一个只依赖分形参数、材料参数和压入深度数组的黑盒函数，方便批量研究参数影响，也可以直接作为后续优化设计的子函数调用。
这个流程跑通之后，你可以很快扩展到更多工况，比如加入粘着效应、引入双变量粗糙度、或者把计算出的刚度嵌入多自由度结合面动力学模型。分形接触的魅力就在于：表面看起来杂乱无章，但用一个维数和一组尺度参数就能把它的承载行为描述出规律。你只要在MATLAB里把这个核心代码吃透，后续很多接触问题都能在此基础上做二次开发。
本文还有配套的精品资源，点击获取
