Advances of Phase Field Modeling of Martensitic Phase Transformation

界面摩擦是材料的固有属性之一^[1]，普遍存在于各种材料和结构中。摩擦界面的宏观运动遵循经典摩擦定律。然而，若探究摩擦的物理机制，需在微细层面研究界面接触的结构特征和界面运动的微过程。研究界面静/动转化的传统方法一般通过静摩擦系数（f_s）和动摩擦系数（f_d）的变化来描述，然而这种描述过于粗糙，无法刻画界面运动状态转化的细节和机制。事实上，界面由静到动变化的瞬间，界面存在显著的波动现象和波的精细结构变化^[2-4]。弄清这些波动规律，探索波动的起源和物理机制，将有助于深入理解摩擦的微观本质，拓展对界面摩擦动力学行为的认知，从而为工程中构件的摩擦设计提供理论支撑，为自然界各种摩擦运动（如地震等）的预判提供新的波动方法。

目前，对界面摩擦行为已有较多研究，大多集中在界面的摩擦特征，研究主要包括两个方向：(1) 探索各种物理过程中不同尺度下界面摩擦系数、摩擦运动的规律^[5-11]；(2) 关注特定尺度界面的摩擦运动对更大尺度物理过程的影响，如界面粗糙度对撞击或侵彻的影响^[12-17]、纤维/基体界面强度对整体复合材料性能的影响^[18]、颗粒接触摩擦对颗粒群剪切流动变形的影响^[19-20]以及摩擦运动对地震的影响^[21-23]等。上述研究工作涉及的时间尺度相对较长，贯穿整个相关物理过程。因此，不论是研究界面摩擦运动规律还是其影响，都没有过多关注界面摩擦运动早期由于摩擦运动导致的界面波动现象。事实上，不存在理想的准静态荷载，扰动的外荷载可诱发界面波动的发生和传播，这些波动信号携带了界面的实际微接触几何、物理特征和摩擦微运动信息。实验上实时检测该波动信号，理论上建立波结构和界面接触微结构的关系，从而分析其传播及演化规律，将有助于全过程理解界面的动态摩擦力学行为。

实验方面，直接观察摩擦界面上的波动现象是极具挑战性的，目前只有少数实验做到了这一点。Pyrak-Nolte 等^[24]首先在断裂面上观察到一类新的弹性界面波，这类界面波包含两类，即快波（fast waves）和慢波（slow waves），其产生机制同界面的断裂刚度有关。Xia等^[25]在有机玻璃摩擦界面滑动过程中，通过光学方法直接观测到界面上的亚瑞利波和超剪切波。Rubinstein等^[2-3]也开展了有机玻璃准静态剪切实验，通过光学方法精细测量了界面微滑移过程中实际接触面的几何形态变化，结果表明，在界面由静到动的微转化瞬间，接触界面产生并传播了3种应力扰动，即快瑞利波、超剪切波和慢波。Ferrer等^[26]实时测量了界面摩擦运动瞬时来自于界面的声辐射信号，认为这些信号幅值和数量的变化可以用来界定界面由静到动转化的特征时间节点。上述实验揭示了一个重要现象：在摩擦界面由静到动转化的瞬时，界面的运动和变形行为是个微动态过程，并伴随显著的波动特征。

理论方面，Braun等^[27]建立了弹簧-滑块模型，采用接触刚度综合描述界面的接触特征，结果表明，在界面运动尚未转为整体摩擦滑移之前，由于界面的波动作用，界面初始未受扰动区已经转化为高度非均匀应力状态。Svetlizky等^[28]提出了粗糙界面有限断裂的理论模型，认为界面精细波结构的产生与界面上大量实际微接触群被瞬时剪断事件有关，基于模型模拟计算的波结构与实验中的慢波^[2-3]吻合较好，但与其他两种波有一定差异。Bartolomeo等^[29]建立了界面波结构特征和界面运动转化的联系，认为弄清该波结构特征和传播规律将有助于理解界面摩擦运动的机制。Kammer等^[30]建立了断裂能与界面摩擦强度的关系模型，其三维有限元模拟计算表明，界面存在亚瑞利波和超声波，这两类波的产生和转化取决于两种尺度的相互作用，即界面黏性区域宽度和非均匀区特征尺度。事实上，各种尺度的界面均是粗糙不平的，这种粗糙性在更小尺度上表现为界面的实际自然接触具有更精细的微尺度起伏性，这些起伏的结构特征和演化规律必然深刻影响界面的摩擦强度和动力学行为^{[1, 5]}。界面摩擦运动诱发的界面波动本质上来源于界面实际接触中的几何结构变化，也必然蕴含了这些接触的特征信息。然而，困难在于难以对界面上波的精细结构进行定量描述和解释，从而建立界面动力学行为和波结构的联系。

本工作主要研究界面微接触断裂事件引起的波动效应，为此建立了简单界面接触模型，以模拟界面微滑动过程中微接触的断裂行为以及界面上相应的波结构特征。该模型将展示摩擦滑动过程中界面微凸起断裂的演变过程及其相应的波结构，以揭示波的精细结构特征和产生机理。

1.   平面摩擦滑动模型

采用ABAQUS软件建立了平面摩擦滑动模型，如图1所示。上基体（Part-1）、下基体（Part-2）的宽度均为10 mm，高度分别为40和10 mm。为表征界面的粗糙性，宏观上赋予界面摩擦系数，细观上在界面中心设计了边长为0.1 mm的等边三角形微接触凸起。考虑到摩擦界面的地震背景，摩擦界面的样本材料选择弹脆性玄武岩体，材料的破坏采用D-P准则描述，相关弹性参数和断裂参数见表1。应力脉冲荷载作用下，摩擦界面的动力学行为是一个微过程的波传播问题，为展现波传播的细节和精细结构，单元网格的划分非常小，尤其是在微凸起近区域，最小网格尺寸5 μm，时间步长10⁻⁹ s。

图  1  摩擦界面模型

Figure  1.  Friction interface model

下载: 全尺寸图片 幻灯片

t=0时刻在上部滑块左侧施加瞬态荷载（σ^*）10 MPa，下滑块底面为固定约束，其他边界设置为无反射类型。如图1所示，为更好地揭示波动过程和应力波的精细结构，选取两类单元：(1) 在界面上以微凸起为中心，上下、左右对称选取8组单元，形成4个区域，通过这些单元的应力历程展现界面上的波动效应和精细波结构；(2) 在微凸起的正上方，旨在揭示基体内的波结构。

表  1  计算材料参数

Table  1.  Material parameters for calculation

Density/ (kg·m−3)	Elastic modulus/GPa	Shear modulus/GPa	P wave velocity/(m·s–1)	S wave velocity/(m·s–1)	Friction coefficient	Internal friction angle/(°)
2300	62.8	24.1	5225	3237	0.1	44
Expansion angle/(°)	Hardening coefficient	Fracture strain	Tensile strength/MPa	Cohesion strength/MPa	Shear stress ratio	Absolute plastic strain
0	6.98	0.0075	3.5	8	0.33	0

下载: 导出CSV 

| 显示表格

2.   计算结果及分析

2.1   界面上的波结构特征

图2为上界面单元的3种应力波结构特征。从上到下依次为σ₁₁、σ₂₂和σ₁₂，从左至右依次为Zone-1、Zone-2、Zone-3和Zone-4。σ₁₁和σ₂₂分别为图1中x和y方向的正应力，σ₁₂为单元的剪应力。图2中σ₁₁的分布演化表明，在3 μs微过程中，σ₁₁的变化呈现出显著的3波结构特征。最先到达的波速度约为5200 m/s，幅值为10 MPa，显然该波来自于加载端面的平面纵波；后继两个波的平均速度约为5200和3000 m/s，可能是纵波和Rayleigh界面波。由于Zone-2和Zone-3、Zone-1和Zone-4均关于微凸起对称分布，从波的到达时间来看，这两个扰动来自于微凸起的断裂，即微凸起断裂瞬时，沿接触界面将传播纵波和界面波，速度约为0.93C_s。从σ₂₂和σ₁₂波结构来看，来自于加载端面的纵波同样会引起界面单元竖直方向应力分量σ₂₂和剪应力σ₁₂的变化，幅值强度分别为1.22和0.12 MPa，但二者跳跃变化的机制不同，σ₂₂的变化源于泊松效应，而σ₁₂的微扰动与界面的微摩擦相对滑动有关。另外，也可以看到，微接触凸起断裂诱发的纵波和界面波不会引起Part-1部分界面单元σ₂₂应力分量变化，但Rayleigh波会导致该类界面单元应力σ₁₂的扰动。值得注意的是，在σ₂₂的波形结构中，左端面纵波未到之前已存在一个微扰动（图2和图3中σ₂₂应力波动结构图中灰色矩形框），并且这些微扰动信号的幅值一致，时间同步（到达时约为0.3 μs）。对比相应的σ₁₁和σ₁₂的波结构发现，并不存在类似波形，说明该扰动并不是沿着界面传播的，而是起源于界面，进一步的物理机制需要结合微凸起的断裂过程来深入分析。

图  2  Part-1上界面单元的波动特征

Figure  2.  Wave structure of interface elements at Part-1

下载: 全尺寸图片 幻灯片

图  3  Part-2下界面单元的波动特征

Figure  3.  Wave structure of interface at Part-2

下载: 全尺寸图片 幻灯片

图3为下界面单元的3种应力波结构特征。在图3所示σ₁₁的波动信号中，也呈现类似的三波结构：其中后到达的两个波同界面上方对应单元一致，均来自于界面微凸起断裂产生的纵波和界面波；与上部对称单元相比，最先到达的波的速度和到达时间是一致的，说明两个波同源，但其幅值差别较大，这是因为上部单元位于Part-1上，承受主动荷载，而下部单元在Part-2上，其σ₁₁的变化主要在于界面的摩擦作用。对比图2和图3的波形结构，可以看出，除个别点外，界面上下对称单元的剪应力σ₁₂和垂直应力σ₂₂的扰动变化几乎是完全一致的，这源于扰动起源的一致性和所取单元位置的对称性。另外，在σ₂₂的变化过程中，相同时刻主动加载纵波未到达之前也出现了相似的微扰信号。

综合图2和图3展示的界面波动信号可以看出，在施加载荷的瞬间，界面存在显著的波动现象和精细的波结构，这种精细结构特征源自于加载波和微凸起的相互作用造成的微接触断裂。

2.2   应力波与界面微缺陷相互作用过程

图4以应力云图的形式给出了加载弹性波的传播图像及与界面微凸起的相互作用过程。由于图4主要展示了波传播过程及波的精细结构特征，因此并没有绘出反映应力大小的标度柱状图，这样使得图4中波结构更简明清晰。t=0时，加载弹性波包含3个应力扰动，即σ₁₁、σ₂₂和σ₁₂以平面纵波的形式向基体内传播。图4中，t=0.905 μs时，加载纵波传播至界面微凸起位置，并与之相互作用，由于应力集中引起微凸起断裂，从而以断裂点为中心，形成球面P波和S波向上、下基体内传播，同时沿着界面形成Rayleigh波，图4(b)中σ₂₂的变化清晰展示了该过程。对于界面上的单元，只能感应到P波和Rayleigh波，因此在图2和图3中，在主动加载波后，σ₁₁为双波结构（P波和Rayleigh波），而σ₂₂和σ₁₂只有Rayleigh波，即界面上的单元无法反映由于微凸起断裂而形成的S波。图4(b)中，自t=0.311 μs开始，可以清晰看到两个P波自界面分别向上、下基体中传播，如图2和图3中σ₂₂信号中的灰色区域所示。剪应力σ₁₂的扰动波阵面的演化如图4(c)所示。整体来看，剪应力σ₁₂的扰动阵面是以界面为中心呈近似对称平面圆锥形，随着主平面纵波σ₁₁阵面一起传播，在传播过程中，圆锥面逐渐向上、下基体扩大。

下载: 全尺寸图片 幻灯片

图  4  应力波的传播及演化

Figure  4.  Propagation pattern of stress wave in space-time

下载: 全尺寸图片 幻灯片

2.3   微接触的断裂过程及波动效应

如图4所示，在应力波演化过程中的一个重要事件是界面上微接触发生断裂，但图4并没有清晰地描绘该过程中的断裂过程以及相应波的精细结构。图5通过σ₂₂应力云图给出了微凸起断裂的微过程。图5中，当t=1.272 μs时，由于应力集中，三角形微凸起的左下角首先开始起裂，接着裂纹向下部基体和沿凸起根部传播（t=1.562 μs和t=1.604 μs）。与此同时，微凸起的右下角位置也开始出现裂纹（t=1.660 μs），两个角裂纹沿着界面运动（t=1.724 μs），最终微凸起从根部剪断（t=1.740 μs）。在该过程中，主要以微凸起的左、右下角为中心形成新的应力扰动。

图  5  微凸起的断裂过程云图

Figure  5.  Stress nephogram of micro bulge fracture process

下载: 全尺寸图片 幻灯片

微凸起断裂形成新的次声波结构如图6所示。在t=1.780 μs时的应力云图中，可以清晰地观察到纵波、横波和界面波三波结构。通过图2可知，该界面波主要引起界面σ₁₁和σ₁₂的变化，其传播速度约为0.93C_s。显然，只有上、下基体内的单元能感应到纵波和横波扰动。为证实这一点，提取微接触的正上方单元的应力扰动变化信息，如图7所示。图7中，当断裂引起的球面纵波传播至捡取单元时，主要引起σ₂₂和σ₁₂的较大变化，σ₁₁有微小扰动，而横波主要体现在σ₁₂的变化上。因此，从σ₁₂的扰动信号中可以清晰地观察到凸起断裂形成的双波结构，其波速分别为5225和3239 m/s。

图  6  微凸起断裂的波动效应

Figure  6.  Wave effect induced by micro bulge fracture

下载: 全尺寸图片 幻灯片

图  7  微凸起正上方单元的应力扰动

Figure  7.  Stress disturbance of elements directly above the micro bulge

下载: 全尺寸图片 幻灯片

3.   讨　论

有限元模拟结果表明，在施加载荷的瞬间，摩擦界面存在显著的波动效应。界面上波的精细结构主要包含3部分：主动加载脉冲、界面微凸起断裂引起的球面扰动和起源于界面的微扰动。基于传统的固体波动理论和特征线分析方法^[31]，可以在x-t平面上绘出前两种波的传播及结构特征，如图8所示。由于模拟对象处于近似一维应变状态，因此初始时刻的加载弹性波包含3个应力扰动，σ₁₁、σ₂₂和σ₁₂独立向基体内传播。图9给出了t=0.315 μs时3个应力波动的云图。宏观上，3个应力扰动的前沿阵面传播速度是一致的，但由于摩擦界面效应，可以清晰地看到，界面近区域波前沿阵面的形状发生了显著变化。具体而言，σ₁₁在波阵面的右下尖端有应力集中带并延伸至下基体，σ₂₂在界面近区域形成关于界面对称的三角锥，剪应力σ₁₂的波阵面整体上是平面纵波，但在界面近区域呈现出较大的对称圆锥，这种近似圆锥形的剪应力阵面同Rubino等^[23]的实验结果是一致的，但其形成机制需要进一步深入研究。

图  8  基于特征线理论的波结构

Figure  8.  Wave structure based on characteristic line theory

下载: 全尺寸图片 幻灯片

图  9  t=0.315 μs时3个应力扰动的波阵面形状

Figure  9.  Wave fronts of three stress disturbances （t=0.315 μs）

下载: 全尺寸图片 幻灯片

图2和图3的σ₂₂波结构中还存在一个精细的微小扰动。图4(b)所示的σ₂₂应力云图表明，该扰动在t=0.300 μs开始在界面上出现，此时左端面的加载纵波传播了约1.56 mm，这对于介质的内部单元（尤其是近界面区域）来说，意味着波尚未传至时已产生微应力扰动。该扰动并没有沿界面传播，而是起源于界面且向上、下基体材料内传播。为进一步验证这一点，在微凸起正上方取一系列单元，其σ₂₂扰动信号如图10所示，对灰色框所示的信号放大，从而清晰地观察到一个速度约为5220 m/s的微扰动信号自界面处向上传播，显然该波是起源于界面的平面纵波。

图  10  σ₂₂波的精细结构

Figure  10.  Fine structure of σ₂₂ wave

下载: 全尺寸图片 幻灯片

基于上述分析形成了初步认识，该波起源于界面，是平面纵波，但产生的物理机制仍不清楚。调整界面摩擦系数（μ=0，0.3，0.5）对该波的形成没有影响（图11(a)），推测该波可能是微凸起引起的，但改变微凸起的大小（三角形微接触边长a=0.1，0.5，1.0 mm）和形状甚至将微凸起去掉后，发现该波仍然存在（图11(b)）。基于对该波到达时的分析，推测该波可能与上部基体的重力作用于界面有关，为此进一步设计了两个模型：一个是将界面去掉，上下为一个整体，此时该波不存在；另一个是将原模型顺时针旋转90°，即重力不再直接作用于界面，此时该波亦不存在。上述两种模型说明该波的形成可能与作用在界面上的重力调整有关。如图12(a)所示，根据惠更斯原理，在外部荷载作用下，界面上的每个点都是一个新的波源，产生球面扰动，这些球面波系的前沿波阵面形成一个包络面，该包络面的传播速度与球面P波一致，形式上表现为平面纵波的特征。图12(b)所示的数值模拟结果清楚地展现了界面上球面波系与材料内向上下对称方向传播的新纵波扰动的关系。但这里仍然有些基本问题有待商榷，如：界面上单元最初的扰动从哪里来，是否与重力作用有关系，该波出现的时间是否固定，同哪些因素有关。

图  11  界面性能对新纵波扰动的影响

Figure  11.  Effect of frictional interface properties

下载: 全尺寸图片 幻灯片

图  12  界面新纵波扰动的形成机制

Figure  12.  Mechanism of interfacial longitudinal wave

下载: 全尺寸图片 幻灯片

该波的存在给了我们深刻的启示。目前，地震预测系统建立在地震P波和S波之间的速度差上，当地震发生时，基于本研究中的模型结果，在球面P波抵达地表探测系统之前，应有重力扰动引起的纵波存在，若能捡取、分离、明确该信号扰动，将有助于将地震预报的时间提前。为了验证其可行性，进行了宏观大尺度地震模拟，如图13所示。显然，先于地震P波（图13(b)中蓝色虚线部分）的纵波扰动（图13(b)中红色虚线部分）是存在的，这方面的实验和进一步的数值模拟研究将在后续工作中逐步开展。

图  13  地震波结构模拟

Figure  13.  Simulation of seismic wave profile

下载: 全尺寸图片 幻灯片

4.   结　论

建立了简单的界面摩擦模型，通过数值模拟计算分析外荷载作用瞬时界面上的波动效应，得到如下主要结论：

(1) 在加载瞬时，界面运动的微过程中界面上存在清晰的波的精细结构；

(2) 界面上波的精细结构与界面微接触的断裂有关，断裂将产生纵波、横波和界面波；

(3) 在界面微凸起断裂之前，由于重力扰动的作用，自界面产生微扰动，该微扰动以平面纵波的形式向基体内传播。

客观上，界面的粗糙度和起伏度要复杂得多，但本研究通过一个简单的三角形微凸起模型揭示了一个重要规律，即界面的摩擦动力学行为与界面的粗糙起伏引起的断裂密切相关，这种断裂事件将以波动的形式在基体和界面上传播。当考虑界面的实际粗糙度时，相应波结构产生的变化方面的研究将在后续工作中继续开展。另外，应力扰动σ₂₂和σ₁₂的传播过程中波阵面形状的变化机制也需要进一步深入分析，给出合理的物理解释。