论文标题：基于数字体相关和有限元模型修正的金属基复合材料损伤机制识别
摘要
本文介绍了一种对金属基复合材料在压缩载荷下损伤机制的三维实验与数值研究。所研究的复合材料是重型刹车片中常用材料的简化版本，仅包含九种组分。该分析结合了实验室X射线计算机断层扫描、数字体相关技术，以及基于图像的有限元方法和被称为有限元模型修正的反向识别方法。原位CT与DVC的结合使我们能够理解控制该非均质金属基体内应变局部化和损伤累积的机制。损伤主要表现为位于大尺寸石墨G2颗粒之间的薄“胞壁”区域材料断裂。此外，本文开发的FEMU方法能够识别控制宏观压缩行为的两个关键组分（即“基体”——其中合并了九种组分中的八种，以及G2颗粒）的弹性模量，并预测基体内的损伤区域。研究强调，在FEMU中引入损伤机制对于可靠地预测复合材料的压缩响应至关重要。

1. 引言
   重型车辆对制动性能的持续高要求，推动了新型金属刹车片衬片材料设计方法的发展。这些材料需要同时满足高摩擦系数、高导热性、耐用性、耐腐蚀性、低噪音和成本效益等多重苛刻要求。因此，对材料配方的精确控制至关重要。传统的“试错法”研发模式效率低下且成本高昂。利用逼真的数值模拟可以显著提高开发过程中的预测能力，实现高通量的材料筛选和设计。然而，有效的数值模拟需要深入理解摩擦材料在实际工况下的力学、热学、摩擦学和动力学行为。

刹车片衬片通常包含多达20种不同组分，可分为磨料、摩擦改性剂、填料和粘结剂。以往研究虽然探索了配方和服役引起的组分性能变化，但将这些变化与材料行为的预测性建模联系起来仍存在空白。例如，材料在宏观和微观尺度上的力学响应直接影响其摩擦性能。此外，实践中有些组分的性能很难通过实验测定。因此，开发一种整体的实验/数值方法，以识别复合材料中有效组分的性能，并模拟损伤局部化的机制至关重要。

原位X射线计算机断层扫描已被广泛用于三维分析材料在微观结构层面的变形方式。除了直接观察图像，CT还能用于创建基于真实微观结构的有限元模型，并应用DVC技术来相对直接地评估位移场和应变局部化。具体来说，基于DVC提供的全场应变数据，可以运用有限元模型修正等反演方法，迭代更新FEM中未知的材料参数，使计算应变场与测量应变场相匹配。DVC与反演方法的结合已被用于计算球墨铸铁、云杉木材、弹性体泡沫和核石墨等多种材料的宏观本构特性。然而，目前仍缺乏旨在确定细观和微观尺度组分性能，以及开发简化损伤描述方法的研究。巨大的计算成本是其广泛应用的主要障碍。

为此，本文引入了一种结合原位CT、DVC和基于图像的FEM的混合方法。由于宏观应变局部化主要由G2颗粒控制，本文将3D图像基FEM简化为双材料模型。随后，开发了一种基于FEMU的反向识别方法，并引入损伤准则来真实地模拟复合材料内的损伤累积。

Fig. 1. (a) 典型火车刹车片图像，刹车片衬里材料被分成5个独立的闸瓦。 (b) 金属刹车片衬里微观结构的SEM-BSE图像，显示了低倍率下的主要组分。高倍率插图详细展示了金属基体，其中较亮的元素是铜。

Fig. 2. (a) 宏观压缩测试装置图，使用从刹车片衬里中提取的立方体样品。 (b) 右侧应变计记录的应力-应变曲线。蓝线对应应力-应变曲线加载段的线性回归拟合。 © 弹性模量随应力水平下降的量化曲线图。

2. 材料与方法
   2.1. 刹车片衬里材料
   研究对象是一种半金属刹车片衬里复合材料。为便于阐明损伤机制，其成分被简化为仅包含九种对刹车性能至关重要的主要组分：铁粉和铜粉等金属粉末、填料、陶瓷颗粒和两种不同类型的石墨颗粒。填料理化性质属于专利，无法公开。因此，该材料被视为一种铁基材料，各组分的尺寸分布和体积百分比如表1所示。

(此处插入 Table 1)

Table 1
主要组分的平均尺寸和体积分数汇总表。

制造过程包括粉末混合、压制和烧结。G1是次生合成石墨，G2是原生合成石墨。由于两者在显微镜下灰度相近，难以区分，因此分类主要依据尺寸和形状：G1呈片状，G2球形度更高且尺寸更大。由于它们在力学性能上扮演的角色截然不同，本研究的建模方法将明确分割出G2颗粒，而将G1颗粒与陶瓷、填料和金属基体一并归入“基体”相。

2.2. 宏观压缩试验
使用INSTRON-5500万能试验机进行宏观压缩试验。加载方向平行于Z轴，即压制方向。试验采用立方体试样，应变速率为0.01 mm/s。实验装置包括力传感器、两个应变片和一个用于拍摄表面图像并进行数字图像相关分析的相机。我们先前的研究结果表明，宏观应变局部化受G2颗粒周围及其内部的应变控制。

压缩试验由一系列在递增载荷下进行的加卸载循环组成。进行这些间歇性加卸载循环的目的在于：(1) 追踪弹性模量随施加应力增加而逐渐退化的过程，(2) 捕捉非线性和迟滞行为，这为损伤起始、累积和应变不可逆性等微观结构机制提供了见解，(3) 识别损伤定律的参数，因为刚度随应力的变化是校准所选损伤模型参数所必需的。应力水平的选择基于以下考虑：(i) 远低于预期损伤阈值开始，(ii) 保持在材料的典型服役范围内，(iii) 略微超过屈服点以诱发可观察的损伤并确保力学响应的可测量变化，(iv) 避免整体开裂或宏观失效。

不同载荷水平下的应力-应变行为遵循相似模式：初始加载为线弹性阶段，随后是准线性阶段，该阶段主要发生拉伸变形，累积应变。卸载时，材料表现出非线性行为，重新加载时又恢复线性。卸载-再加载分支的响应具有显著的迟滞性。这种迟滞以及卸载时的非线性，主要由G2颗粒控制。

应力-应变曲线的另一个重要特征是弹性模量随载荷增加而降低。20 MPa时的弹性模量值几乎是3 MPa时的一半，且从5到20 MPa，这种下降趋势近似线性。弹性模量的逐渐降低归因于基体中的损伤累积，尤其是在颗粒间薄区域。这种效应即使在中等载荷下也会出现，这已被DVC和SEM分析所证实。

2.3. X射线计算机断层扫描与原位压缩试验
实验室X射线显微CT在里尔的LaMcube实验室ISIS4D平台上进行。使用了两种尺寸的样品以获得不同的体素尺寸。直径15 mm、高24 mm的样品用于观察能代表宏观响应的压缩行为。原位加载装置置于CT旋转台上。记录了0 MPa和15 MPa两个载荷下的CT数据，用于后续的DVC分析。体素尺寸为(22.7 μm)³，该分辨率只能分辨G2颗粒，而G1颗粒几乎无法区分。因此，在此配置下，复合材料可简化为双材料模型。

直径5 mm、高15 mm的样品用于分辨大多数组分并观察损伤机制。该样品先在0 MPa下扫描，然后通过试验机进行60 MPa的非原位加载，之后再次扫描。高分辨率X射线管产生的体素尺寸为(4 μm)³。这两个数据集并非用于直接的力学响应比较，而是互补使用：一个用于定量建模，另一个用于定性观察损伤机制。

Fig. 3. XCT实验示意图，(a) Ø5 mm样品，(b) Ø15 mm样品。 © Ø5 mm样品的重建断层图，体素尺寸为(4 µm)³。 (d) Ø15 mm样品的重建断层图，体素尺寸为(22.7 µm)³。

使用Fiji软件处理Ø15 mm样品的CT数据。采用局部厚度法计算基体材料的厚度，该方法为每个体素计算能恰好包含该点的、内切于物体的最大球体的直径。

2.4. 数字体相关
使用YaDICs软件计算Ø15 mm样品的位移和应变场。选择了具有良好空间相干性的光流有限元方法。微观结构提供的天然散斑图案允许使用最小为16×16×16体素的单元尺寸。总体残差水平较低，表明配准成功。使用了掩膜来减少边界条件的影响。YaDICs软件能实现亚体素精度，位移不确定度通常低于2 µm。G2颗粒与基体间的对比度足够清晰，因此CT噪声的影响被认为是次要的。

应变基于DVC获得的位移场U进行计算：

e_{ij} = \frac{1}{2}(\frac{\partial u_i}{\partial x_j} + \frac{\partial u_j}{\partial x_i})
e
ij

2
1
(
∂x
j

∂u
i

∂x
i

∂u
j

)

3.1. 有限元模型修正
将Ø15 mm样品在0 MPa下的分割数据导入Avizo软件，生成微观结构的表面网格。原始表面网格被迭代简化，并重新划分网格以提高三角形质量。使用GMSH软件生成约250万个节点的双材料体网格。进行了网格无关性验证，确保识别出的参数和应变局部化模式不受离散化显著影响。

FEMU使用Abaqus商业软件进行。边界条件源自DVC获得的位移场，作为狄利克雷边界条件施加于上下表面。这确保了模型反映实验中观察到的实际运动学，避免了对应力或接触载荷的假设。

Fig. 4. (a) Ø15 mm样品CT数据的3D渲染图。 (b) 微观结构的表面网格。 © 表面网格细节，显示三角形质量。 (d) Abaqus生成的3D渲染图，显示施加在顶面和底面的边界条件。

3.1.1. 双材料的反向识别
反问题的目标是识别未知的本构参数θ。解决此类问题常用三类方法：基于平衡方程的方法、最小化本构方程误差的方法，以及最小化数值量与实验位移场之间差异的有限元模型修正方法。本文将采用FEMU。

Fig. 5. 本研究中开发的反向识别过程示意图。

FEMU的核心是最小化代价函数，该函数评估每个网格节点上实验（DVC）与数值（FEM）应变场之间的差异：

\text{Cost function} = e(\theta) = \sqrt{\sum_i \sum_j (\Delta e_{ij})^2}
Cost function=e(θ)=
i
∑

j
∑
(Δe
ij
)
2

使用Levenberg-Marquardt算法迭代更新材料参数，直到达到收敛准则。

3.1.2. 考虑基体损伤定律的反向识别
为模拟损伤，引入一个标量损伤变量D。基于连续介质热力学，结合弹性和损伤的自由能表达式为：

\Psi = \frac{1}{2} (1 - D) (2\mu \text{tr}(\varepsilon : \varepsilon) + \lambda (\text{tr}(\varepsilon))^2)
Ψ=
2
1
(1−D)(2μtr(ε:ε)+λ(tr(ε))
2
)

本构关系为：

\sigma = \frac{\partial \Psi}{\partial \varepsilon} = (1 - D) (2\mu \varepsilon + \lambda \text{tr}(\varepsilon)I)
σ=
∂ε
∂Ψ
=(1−D)(2με+λtr(ε)I)

损伤驱动力Y为弹性应变能释放率。损伤演化由下式控制，其中X是Y在加载历史中的最大值。

为最小化待识别参数，选择了简单的损伤模型，其函数F(Y)定义为：

F(Y) = \frac{1}{1 + b_1 (Y - Y_0) + b_2 (Y - Y_0)^2}
F(Y)=
1+b
1
(Y−Y
0
)+b
2
(Y−Y
0
)
2

1

参数Y0, b1, b2通过宏观压缩试验结果（图2c）使用最小二乘法确定。表2总结了模拟中使用的所有参数。

(此处插入 Table 2)

Table 2
FEM模拟中使用的识别和假设的材料参数汇总表。

泊松比取0.15是合理的，因为：(i) 反向识别的主要目标是有效弹性模量，在单轴加载下，应变场对泊松比的敏感性远低于对杨氏模量的敏感性；(ii) 实验数据主要提供轴向位移，限制了同时稳定识别泊松比的能力。虽然未进行完整的不确定性量化，但对损伤参数±10%变化的敏感性测试显示，其对局部化模式和预测模量的影响是适度的。

尽管所选损伤演化律相对简单，但其选择是出于以下考虑：(1) 最小化待识别参数，确保数值稳定性和避免过拟合；(2) 宏观拟合目标旨在捕捉损伤累积的整体效应；(3) 有效性已得到证实。模型未包含显式塑性，因为主要的非线性与损伤累积相关。

3.2. 损伤机制的确定
与DIC结果一致，DVC结果显示应变场是非均匀的，在某些区域出现应变带。通常认为，基体承受的应变水平最低，而G2颗粒控制着应变的局部化。有趣的是，并非所有G2颗粒都会导致高应变局部化，基体本身也可能局部化相当程度的应变。

因此，G2颗粒的存在是应变带出现的必要条件。基体局部厚度的计算显示，其变化范围从最薄区域的20-50 μm到400 μm。这些最薄区域被称为“胞壁”区域。沿Z轴每层切片的平均厚度计算表明，应变带倾向于出现在平均厚度局部极小值的区域。因此，G2颗粒之间薄胞壁区域的出现被认为是应变带发展的另一个必要条件。

Fig. 6. Ø15 mm样品在15 MPa名义压缩载荷下的XCT重建切片，叠加了DVC结果。细节显示应变作为G2颗粒函数的局部化情况，© G2颗粒未导致高应变局部化，(d) G2诱导了应变带。

Fig. 7. (a) Ø15 mm样品在15 MPa下对应XY平面的XCT重建切片。 (b) 相应的局部厚度计算结果。 © ZY平面的XCT切片，叠加了相应的DVC计算结果。 (d) 沿Z轴每个XY切片计算的平均厚度和平均Ezz值。

文献报道，裂纹数量增多通常导致材料弹性模量更严重地下降。在本研究中，损伤倾向于从最薄的胞壁处开始逐渐累积。这种损伤主要是断裂主导的。在自由表面，裂纹可能表现为胞壁向外翘曲、壁的一部分被破碎并向外移位，或壁只是向外弯曲而没有明显开裂。基体在某些区域可能被严重破碎。此外，位于裂纹路径上的陶瓷和石墨颗粒可能发生脱粘。

Fig. 8. (a) Ø5 mm样品的XCT切片，显示其复杂微观结构。 (b-e) 加载至60 MPa后，不同胞壁的损伤机制。 (f) 加载后样品表面的SEM图像。

基体内部的损伤机制更难观察，因为石墨的存在限制了胞壁的屈曲和弯曲。无论如何，由于应变带也发生在材料内部，预计内部壁也会因相邻材料的位移而产生裂纹。

3.3. 损伤机制的建模
经过几轮FEMU迭代后，识别出的G2弹性模量为0.2 GPa，基体为10 GPa。由于实验压缩载荷并非与轴线完全平行，导致样品一侧承受更高的压缩位移。FEMU计算高估了材料的抗压强度，因为它无法正确复制基体内的应变局部化。

Fig. 9. (a) 使用未考虑损伤的FEMU计算的Ezz云图。 (b) 引入损伤定律后的相同计算结果。在©和(d)两个计算的细节中，渲染了对应于-0.2% Ezz水平的等高线。 (e) 对应于©和(d)所示区域的DVC结果细节。

通过在反向识别中引入损伤定律，重新计算出的G2颗粒压缩模量为0.27 GPa，而基体的杨氏模量约为9.5 GPa。这使得Ezz应变场现在能够预测基体内的应变局部化，计算结果与DVC结果吻合良好。从本质上讲，重新计算得到的9.5 GPa杨氏模量更准确地代表了基体易于损伤的非均质压缩行为。

Fig. 10. 损伤参数图，叠加了DVC应变等高线。

当评估由考虑损伤的反向识别计算所预测的损伤区域时，可以发现这些区域倾向于位于相邻的G2颗粒之间，并在空间上排列成垂直带。这些计算出的损伤区域很可能与最薄胞壁的出现相关。该方法的准确性得到了以下因素的支持：计算与实验应变局部化模式的良好吻合、优化过程中代价函数的收敛，以及使用基于真实CT图像生成的网格。这些结果证明了在反向识别框架中加入损伤定律的价值。从更广泛的角度来看，这种方法不仅能识别有效的材料参数，还能对易损区域进行空间预测。

4. 结论
   本研究结合原位X射线计算机断层扫描、数字体相关和有限元模型修正，来表征和模拟金属基复合材料在压缩载荷下的损伤累积。通过所提出的方法，可以合理地解释和模拟导致宏观应变带形成的损伤机制。主要结论如下：

原位CT压缩试验表明，所研究的金属基复合材料的宏观压缩行为受最大颗粒控制，这些颗粒诱导了宏观应变局部化。因此，复合材料可简化为双材料模型。

应变带的形成与石墨颗粒紧密相邻的情况相关。这导致基体局部厚度显著减薄的“胞壁”形成。

胞壁中的损伤累积，主要表现为材料断裂，被确定为控制宏观弹性模量随载荷增加而逐渐退化的关键因素。

未考虑损伤的FEMU方法确定了双材料的压缩弹性模量，但由于未计入非均质基体的应变局部化效应，该方法高估了材料的承载能力。

在反向识别过程中引入损伤定律，能够更真实地模拟双材料复合材料的压缩响应。通过这种方法，我们能够预测基体内的损伤易发区域。该方法有望为复合材料的微观结构设计提供指导。

CRediT作者贡献声明
Itziar Serrano-Munoz: 撰写-审阅与编辑，撰写-初稿，可视化，方法论，调研，正式分析，概念构思。Vincent Magnier: 撰写-初稿，可视化，方法论，概念构思。Philippe Dufrenoy: 监督，项目管理，经费获取。

利益冲突声明
作者声明，他们没有任何已知的可能影响本文所报告工作的竞争性经济利益或个人关系。

致谢
本研究工作得到了法国国家研究署、国际运输安全与联运校园、上法兰西大区、欧盟共同体、高等教育与研究部以及国家科学研究中心的支持。作者衷心感谢所有这些机构的支持。感谢FLERTEX SINTER公司提供研究材料。

数据可用性
数据将根据要求提供。

参考文献
(此处省略，以保持译文简洁，原文参考文献列表已完整提供)