碳纤维增强复合材料机械连接的可靠性高，承载能力强，是主承力结构的重要连接手段^[1-2]。孔壁作为材料断口和连接表面，是连接件上的应力集中部位和薄弱环节。孔壁表面形貌在微观层面是粗糙的，具有接触的不均匀性，虽然结构尺寸较小，但在很大程度上影响着机械系统的性能及机器的服役寿命等^[3-5]。此外，由于复合材料各向异性和纤维束导向结构，加工形成的孔壁表面几何形貌具有不确定性，表面几何形貌特征复杂且伴随着大量断崖式突降区域和不规则凸起^[6-7]。由于复合材料连接孔壁微观结构尺寸属于微米级且具有不确定性，在连接强度研究中难以使用现有的模型，需要根据实测微观数据进行精准表征重构。

目前的重构方法大多是基于逆向建模技术，根据实体模型，扫描数据，重新生成数字模型用于后续的工程问题研究。Li等^[8]利用三维数据扫描仪测量物理样品的点云数据，借助逆向工程软件Geomagic Studio实现零件的逆向建模和数字模型的误差检测分析。Wang等^[9]利用实验测量得到的数据点云重构具有随机形貌的柔性关节界面，建立2个粗糙随机表面之间的接触模型。但对于复合材料孔壁形貌的重构，上述方法并不适用。现有的逆向建模重构技术是基于插值函数对采集点云进行拟合重构得到连续曲面^[10]。相较于传统插值方法，分形插值根据整体与局部相似的原理，将插值数据点的变化特征映射到相邻点之间的局部区域。复合材料孔壁复杂微观表面上相邻2个信息点之间的连线并不是线性或者光滑的，而是存在局部的变化特征^[11]。传统插值重构方法难以对其进行准确的拟合表征，从而影响后续物理模型的准确性。然而这种2个信息点之间更精细一级的特征，用分形插值方法得到的结果更加符合实际。

分形理论作为研究自然界复杂形状和结构的工具，在复杂结构研究中起着至关重要的作用^[12-13]。随着分形理论研究的深入，针对分形插值应用于粗糙表面的研究，出现了一些实用性的方法。Yang等^[14]分别使用了中点位移法(MD)和Weierstrass-Mandelbrot分形函数法(WM)对一维和二维分形表面进行模拟。该研究突出了2种分形模型的形态特征、统计特性和适用性，为今后分形模型的发展与应用提供了科学依据。翟晓雅等^[15]提出了一种基于分形模型的3D打印路径规划算法，保证边缘细节的精度条件下最大限度地提高打印效率。赵林等^[16]以实测激光切割表面轮廓数据为基础，应用三维表面分形插值算法，对激光切割表面形貌进行分形插值模拟，分析模拟精度。此外，分形模型在海底地形^[17]、地震研究^[18]、风速预测^[19]等领域均有广泛应用。由于周期性的钻孔路径，复合材料连接孔壁具有周期性的形貌特性，较适合使用自仿射的分形插值方法来研究。然而在现有的分形理论中，给定纵向压缩比的分形插值曲面在任意小的局部都与整体的结构相似，只能得到具有周期变化形态的自相似结构^[20]。而复合材料粗糙表面的粗糙度变量是一种区域化的变量，既有随机性，也有相关性。因此需要对现有的自仿射分形插值方法进行改进。

针对复合材料孔壁复杂表面形貌难以快速表征重构的问题，本文以碳纤维增强复合材料连接表面为研究对象，基于自仿射分形理论，提出了一种分形插值重构复杂表面的方法。利用连接表面少量测量单元，对复合材料孔壁表面微观形貌进行了重构；同时利用大量测量单元建立了准确的原始实测表面模型，通过对比实测表面模型与分形插值重构表面模型特征，讨论了采样间隔对插值重构结果的影响，验证重构表面的几何有效性；最后开展了有限元仿真验证重构表面的力学有效性。

1 分形曲面与分形插值 1.1 分形曲面与分形插值

分形是指具有自相似性的几何图形，它们在工业设计领域和自然界中是普遍存在的，例如车轮轮廓、闪电外形等几何形状都具有分形元素。分形曲面在计算机图形学、地震学、岩土学等许多实际领域有广泛的应用。常见的分形曲面主要有：Koch分形曲面、布朗分形曲面、随机网格分形曲面等^[20]。

分形插值曲面是根据分形曲面的自相似性原理和迭代函数，将已知数据插值成具有自相似结构的曲面。与分形插值曲线类似，分形插值曲面任何一个局部都与整体自相似。

1.2 分形插值数学模型

根据分形曲面的理论，学者已经研究了通过迭代函数系统(IFS)或递归迭代函数系统(RIFS)来构造分形曲面^[21]。迭代函数系统以仿射压缩变换为框架，根据几何对象的整体与局部具有自相似的结构，将总体形状以一定的概率按不同的仿射压缩变换迭代下去，直至得到满意的分形图形。在迭代函数系统的基础上，分形插值曲面有多种构造方法。根据复合材料表面几何形貌自仿射特征，本文选用矩形区域分形插值曲面的自仿射分形插值方法。

令I=[a, b], J=[c, d], 设区域D=I×J={(x, y)|a≤x≤b, c≤y≤d}, 以Δx, Δy为步长, 将D剖分为网格

(1)

给定一组网格点上的数据(x_i, y_j, z_{i, j})(i=0, 1, …, N; j=0, 1, …, M), 使用常见的二元分形插值函数f: D→R, 并且满足f(x_i, y_j)=z_{i, j}(i=0, 1, …, N; j=0, 1, …, M)。

在区域(x_n-1, y_m-1, z_{n-1, m-1})→(x_n, y_m, z_{n, m})中, 令x方向的压缩变换为φ_n(x)=a_nx+b_n, y方向的压缩变换为ψ_m(y)=c_my+d_m, 并且需要满足条件

(2)

由此可得

(3)

将结果带入压缩变换中, 得

(4)

令z方向的压缩变换为

(5)

式中, e_{n, m}, f_{n, m}, g_{n, m}, s_{n, m}, k_{n, m}均为待求解的参数且满足

(6)

令s_{n, m}(n∈{1, 2, …, N}, m∈{1, 2, …, M})为自由参数, 称为纵向比例因子, 表示大范围内的数据变化特征映射到小范围区间的缩小倍数, 故需要满足|s_{n, m}| < 1。联立(5)~(6)式得到方程组, 求解得

(7)

将(7)式代入(5)式中得到z方向的压缩变换表达式。其中参数e_{n, m}, f_{n, m}, g_{n, m}, k_{n, m}根据(7)式求得, s_{n, m}为纵向压缩比, 控制分形插值曲面形状, 根据插值形态人为给定。

2 微观表面分形插值重构 2.1 微观表面分形插值 2.1.1 整体框架

本文算法总体框架如图 1所示, 主体分为三部分。主程序中读取实验测量输出的原始数据, 并对原始数据进行填补, 再转换坐标矩阵, 将三维数据文件转换为对应x, y坐标的绝对z值矩阵; 再通过对分形插值参数进行计算, 得到三维压缩比; 最后对z值矩阵进行插值, 得到插值矩阵, 完成分形插值过程。

2.1.2 几何特征提取

实验使用的测量仪器是由奥地利Alicona公司生产的全自动刀具测量仪, 如图 2a)所示, 能够测量零部件的微观表面, 其水平光学分辨率为0.44 μm, 垂直光学分辨率为10 nm, 最小测量高度0.01 μm。实验所用螺栓为航空级钛合金螺栓, 材料为Ti6Al4V, 规格M6。T700级复合材料层合板连接孔名义直径5.98 mm, 铺层顺序[-45/90/45/0]_4s, 纤维体积分数为60%。实验前将层合板沿着孔的中轴线进行横向剖切, 便于后期的测量观察。

首先对螺栓杆杆壁进行测量, 测量区域如图 2b)所示。由于金属的各向同性, 钛合金螺栓只需选取杆壁的一个区域进行测量, 命名为Ti-region。对复合材料层合板的孔壁进行测量, 由于层合板的叠层结构, 需要在孔的入口处至中间层选取4个具有代表性区域进行多次测量, 分别命名为CFRP-region-1, CFRP-region-2, CFRP-region-3, CFRP-region-mid。具体测量区域如图 2c)所示。测量区域为284 μm×216 μm的类矩形区域, 每个测量区域选取采样间隔为4个数据点, 共7 815个测量点。

2.1.3 丢失数据填补

主程序中首先读取原始测量文件, 将测量数据变换为矩阵形式。由于在进行原始测量数据采集时, 所使用的光学显微镜容易受到样本反光率、环境光和景深等因素的影响, 导致表面形貌的个别数据丢失。由于丢失数据点的位置随机, 需要对整体数据进行遍历, 寻找丢失数据点位置并对其进行填补。当丢失点位于中心位置时, 如图 3a)所示, 利用周围8个点的z值数据取平均进行缺失点的填补, 即缺失点z值为

(8)

式中, z_{(i, j)}表示第i行j列的数据点数值。同理, 当丢失点位于边角位置时, 如图 3b)所示, 对丢失的z值进行填补。这种方法符合粗糙表面变化规律, 较好地解决光学测量丢失少量数据点的问题。

将矩阵中三维测量数据转换为对应x, y坐标的z值。由于光学显微镜测量得到的数据是基于全局坐标系, 无法直接进行使用, 需要得到相对坐标后对其进行全局坐标系与局部坐标系转换。

2.1.4 改进的自仿射分形插值重构

在全局统一纵向压缩比的分形插值曲面中, 在任意小的局部都与整体的结构相似, 只能得到具有周期变化形态的自相似结构。复合材料微观表面的粗糙度变量是一种区域化的变量, 即既有随机性, 也有相关性。因此需要使用改进的自仿射分形插值方法, 即让每个局部子区域总生成自己的纵向压缩比。

每个子区间上的压缩比可以通过该方法得到: 首先根据观测值数据, 运用最小二乘法, 得到一次拟合面方程, 该拟合面反映区域观测值的平均平面; 再利用此方程, 求出每个插值点位置上的预测值; 用各点上的观测值减去相应的预测值, 得到偏差值, 偏差值反映观测值偏离平均平面的程度。根据变差函数原理, 数据具有相关性, 某个观测点的变异特征能够反映附近较小领域内的变化特征。因此以观测点偏差值作为纵向压缩比的依据是合理的。

每个子区间上的一次拟合面方程表达式为

(9)

式中，b₀, b₁, b₂为待求系数。根据最小二乘法, 使得观测值与预测值的离差平方和Q最小, 即

(10)

分别求Q对b₀, b₁, b₂的偏导数, 并令其为0。求解整理可得

(11)

求解(11)式, 可得到系数b₀, b₁, b₂。同时将n个观测点数据代入一次拟合面方程式(8)中可得到n个预测值, i=1, 2, …, n。用观测值z减去预测值得到偏差值

(12)

在矩形区域的分形插值中, 规则网格点数据的表示形式为{x_n, y_m, z_{n, m}}(n=0, 1, …, N, m=0, 1, …, M)。类似地, 将预测值、偏差值分别表示为和e_{n, m}。为了保证|s_{n, m}| < 1, 令全局区间偏差值最大值为

(13)

取每个子区间上压缩比s_{n, m}为

(14)

按照改进的自仿射分形插值方法, 对网格点数据进行分形插值重构。对于给定一组网格点的数据{x_i, y_j, z_{i, j}}(i=0, 1, …, N, j=0, 1, …, M), 网格剖分为(1)式所示。x方向和y方向的插值公式由(4)式求出; z方向的插值公式由(5)式求出, 每个子区间的纵向压缩比s_{n, m}由观测值结合(9)~(14)式求出。最后对每组网格数据进行插值重构, 得到插值重构数据。根据x, y, z的值, 重构出表面轮廓。

2.2 微观表面分形插值重构结果

分别对螺栓杆表面区域与复合材料孔壁表面4个特征区域进行分形插值, 初始采样间隔均为4个数据点, 得到多组分形插值数据。其中CFRP-region-mid的表面轮廓重构结果如图 4所示, 图 4a)为实验采集的原始表面轮廓, 图 4b)为分形插值重构表面。分形插值表面与实测表面的几何结构相似且重要几何特征基本相同, 如图 4中的红圈所示。分形插值重构表面更平滑, 相较于原始实测表面更加的简化。

通过上述分析发现, 分形插值后的粗糙表面几何形貌和原始实测粗糙表面几何形貌基本相似, 但在某些细节特征上略有不同。图 5左侧为原始实测表面, 中间为采样表面, 右侧为插值重构表面。分析发现局部细微差异主要为两部分, 首先如图 5a)所示, 当采样间隔设置过大时, 部分包含在采样间隔内的细节特征数据点会被筛选出去, 出现细节特征丢失。其次如图 5b)所示, 当采样间隔过大时, 采样点与采样点间的距离过大, 会导致插值点所依据的采样点超出邻域范围, 不能很好地表征插值点的原始形貌, 导致细节特征畸变。因此, 选取适当的采样间隔对于分形插值重构的结果至关重要, 将在3.1.2节进行详细说明。

3 分形插值重构表面有效性 3.1 几何有效性验证 3.1.1 三维表面轮廓表征与误差分析

三维表面轮廓表征参数主要分为几何参数和性能参数^[22]。由于本研究主要是对复合材料连接孔壁进行表征, 几何参数选用常用的三维轮廓算术平均偏差S_a; 性能参数选用表面接触面积比S_dr。其中, 三维轮廓算术平均偏差S_a表示被测表面上所有像素点到基准面上对应点的直线距离的算数平均值。

(15)

式中: l_{(i, j)}表示像素点到基准面的直线距离; M和N分别为基准面的边长。表面接触面积比S_dr描述采样面积内表面在水平面上的投影面积和表面积的相差程度。

(16)

式中: A表示总展开面积; Δx, Δy分别为x和y方向相邻像素点的间隔距离。表征结果如图 6所示。

对比原始实测表面与插值重构表面的三维表面轮廓参数表征结果发现, 复合材料孔壁4个代表性区域的三维轮廓算术平均偏差S_a的误差均在2%以内, 接触面积比S_dr的误差均在10%以内。通过对重构模型与实测模型表面几何参数分析, 认为分形插值模型具有较好的几何精确性。复合材料表面的粗糙度起伏较大, 并且在分形插值过程中出现细节特征丢失与细节特征畸变现象, 从而导致粗糙表面总展开面积A产生较大波动, 带来接触面积比S_dr的误差。通过上述分析, 结果符合精度要求。

3.1.2 采样间隔对分形插值的影响

由上文发现采样间隔对分形插值结果影响较大。针对此问题进行研究, 又分别进行了相同区域, 采样间隔为2~8的原始采样表面测量实验, 使用上文相同的方法, 得到采样间隔分别为2~8的原始采样表面形貌与分形插值结果并分析。在某些微观细节特征上, 原始实测表面与插值重构表面有细微差别, 并且随着采样间隔的增大, 相似性降低, 会出现原始形貌所不具有的表面特征。使用三维轮廓算术平均偏差S_a进行表征, 如图 7所示。

随着采样间隔的增大, 原始实测表面与插值重构表面三维轮廓算术平均偏差之间的误差增大。当采样间隔超过5个数据点时, 误差量超过预估范围, 分形插值结果与原始表面产生较大差异, 精度无法保证。当采样间隔小于等于5个数据点时, 误差值均小于5%, 符合精度要求。采样间隔越小, 原始实测表面与分形插值表面的差异越小, 难以显示出等效表征的效果。因此, 取采样间隔为3~5个数据点, 既可以保证分形插值精度, 也可以对原始实测表面进行有效的等效表征。

综上, 通过对比重构结果发现, 选取合适的采样间隔时, 分形插值后三维模型的粗糙表面几何形貌和原始三维模型的粗糙表面几何形貌基本相似, 分形插值重构结果具有几何有效性。

3.2 力学有效性验证 3.2.1 有限元模型

利用有限元软件进行力学挤压仿真模拟, 根据挤压力学性能所关注的物理量, 对比原始表面模型与重构表面模型的应力应变以及相互作用面积, 验证分形插值重构结果的力学有效性, 从而验证分形插值重构方法的正确性。

将逆向建模得到的三维实体模型导入有限元软件中, 分别对复合材料板与钛合金板进行材料属性设置、接触设置和边界条件设置, 给予挤压载荷，模拟复合材料装配连接挤压过程。最终得到的有限元模型如图 8所示。

3.2.2 结果分析

将每组原始表面模型与重构表面模型进行有限元仿真, 图 9为CFRP-region-mid区域的仿真结果。根据有限元模拟结果图可以发现原始表面轮廓挤压模型和插值重构表面轮廓挤压模型的接触区域大致相同, 部分细节特征有所不同。对每组模型最大应力值与最大变形量进行分析, 如表 1~2所示。结果发现, 4组模型的最大应力值和最大变形量差距均在10%以内。

利用二值化程序分别对4组有限元变形云图的变形区进行处理计算, 得出4组变形区域与整个面积的百分比, 如表 3所示。通过分析发现, CFRP-region-1的原始表面轮廓与插值表面轮廓的变形区域占比相差较大, 达到9.5%, 4组的差距均在10%以内, 整体趋势与图 6b)中接触面积比误差结果一致。分析重构模型与实测模型力学性能, 认为分形插值模型具有较好的力学预测精确性。分析原因发现, 由于CFRP-region-1区域的接触表面粗糙度较大, 加剧细节特征丢失与细节特征畸变现象。此外, 在进行曲面逆向重构拟合时, 急剧起伏的曲面拟合误差较大, 拟合的效果不太理想, 进而导致原始轮廓模型与插值轮廓模型的误差。

通过对比重构前后力学挤压模拟结果发现，分形插值重构预测模型的力学性能误差均在10%以内，分形插值重构结果具有力学有效性。

4 结论

本文主要提出了一种利用分形插值重构研究复合材料孔壁微观表面等效表征的方法，并使用重构误差分析和有限元仿真验证其有效性。分析得出以下结论：

1) 分形插值过程中，由于细节特征丢失与细节特征畸变，原始表面与插值表面有细微差别，并随采样间隔增大，相似性降低。当采样间隔为3~5时，既可以保证分形插值误差，又可以对原始实测表面进行有效的等效表征。

2) 通过重构的曲面模型，对比原始表面与插值表面几何形貌相似性，三维轮廓算术平均偏差的误差均在2%以内，接触面积比的误差均在10%以内，证明了分形插值重构结果具有几何有效性。

3) 通过有限元仿真，在挤压载荷下，实测原始表面与插值重构表面的最大应力、最大变形量以及接触变形区域面积预测误差均在10%以内，证明了分形插值重构结果具有力学有效性。