刘超，石艺娜，秦承森，梁仙红

(北京应用物理与计算数学研究所，北京100088)

0 引言

冲击相变是当前冲击波与爆轰物理研究领域的热点问题之一，由于其对于高温高压极端条件下材料物性研究的科学意义及在高新技术领域的应用背景，冲击相变研究受到广泛关注。同时，冲击相变还是一个非线性、非平衡的复杂物理过程，具有时间相关效应;并且冲击相变研究既是一个涉及物理、力学、材料科学及冶金学的交叉学科问题，又是一个典型的多个尺度问题。因此，冲击相变也是冲击波物理研究的难点问题之一。自1956 年Bancroft 等采用炸药加载首次观察到铁的冲击相变［1］以来，众多研究者在冲击相变领域积累了丰富的文献资料，同时也发现了许多有趣的实验现象。

时至今日，冲击相变研究仍集中于实验和唯象理论模型方面，对于其细观机理研究相对较少，其中一个很重要的原因是缺乏细观尺度数值模拟手段。尽管多尺度算法研究发展十分迅猛，但多集中于宏观尺度计算方法与分子动力学等微观尺度算法耦合的搭桥类算法，细观尺度的计算方法很少。目前，研究多晶材料冲击响应的细观尺度数值模拟方法，原则上可分为两类:一类是以离散元为代表的粒子类方法，另一类是以连续介质力学为基础的传统数值模拟方法。与传统数值模拟方法相比，离散元具有模型构建方法易行、晶粒间各相异性特性表征便捷、算法实现简单等优点。

离散元法是20 世纪70 年代初由美国学者cundall 首先提出的［2］，最初主要应用于岩石力学、颗粒态群体及土壤力学问题分析中。离散元方法允许单元间的相对运动，而且不象连续介质模型那样的依赖于高度简化、规定性的本构关系，并具有算法简单、易于实现的优点。20 世纪90 年代初，Sawamoto等［3］首先将离散元方法成功地用于混凝土动态冲击破坏等非线性大变形问题的数值模拟研究。刘凯欣等在这一领域做了大量的研究工作［4-5］。1999年以来，Yano 等［6］、Case 等7］利用离散元方法研究了铜、铁等多晶金属的冲击响应。2000 年以后，王文强［8］、于继东等［9］将离散元方法应用到非均质材料炸药在冲击作用下细观损伤的研究。这些工作展示了离散元方法模拟细观非均质材料动力学问题的能力。由于可以方便地表征晶粒间取向的分布特征，离散元方法在模拟细观非均质材料的冲击响应方面具有独特的优势。

本文利用离散元方法结合基于无扩散相变的两相模型，对于α 铁的冲击相变过程进行了数值模拟研究，给出了铁的相边界及冲击Hugoniot 关系。并在此基础上对于多晶铁的冲击相变过程进行了模拟，研究了晶粒大小及传播距离对于冲击波前沿不规则度的影响，并对冲击相变过程中相变特征量随局部冲击压力的变化规律进行了统计。

1 计算方法与相变模型

离散元方法通过求解多体运动的牛顿力学方程组，跟踪全部单元的运动轨迹，来揭示系统与外界的相互作用和自身响应、演化规律。与传统的数值模拟方法相比，离散元具有处理冲击载荷作用下单元间常见的大变形、断裂等问题方便，算法实现简单等优点。

1.1 单元间相互作用力模型

一般单元间可能包括以下相互作用力［10］:1)中心势力;2)中心阻尼;3)弹塑性剪切力;4)切向阻尼;5)干摩擦力。如图1 所示，图中v 为线速度，ω为角速度。

图1 单元间相互作用力的示意图Fig.1 Interaction model of an element-pair

利用图1 描述的单元间相互作用力模型，可将单元i 与j 间的作用力合力表示为

根据所研究的问题选取合适的单元间作用力模型是离散元方法研究的核心内容，本文采用适合于描述材料冲击响应的单元间作用力模型，模型中单元间相互作用力主要包括中心势力与中心阻尼，两相单元间相互作用力模型参数取值见文献［11］。

1.2 温度与相变模型

影响单元温度的力学过程可以分为可逆过程与不可逆过程，可逆的力学过程如中心势力的作用过程，不可逆的力学过程为耗散力的作用过程［11］:

采取与Forbes［12］类似的推导方法，可以得到温度的可逆部分，即等熵过程中温度可表示为

式中:T0为初始温度;p 为静水压力;B0与等温体模量KT相关，B0= KT/β，β 为状态方程参数;γ0为Gruneisen 参数。

不可逆部分仅考虑了由热传导和粘性力带来的能量耗散过程，不可逆过程带来的温升是上述两类过程的累计效应:

连接与接触单元间考虑了基于傅里叶定律的热传导过程:

式中:ΔQ 为在Δt 时间内由单元j 向单元i 传导的热;κ 为热传导系数;Ti和Tj分别为单元i 与单元j的温度;dij为单元i 与单元j 间的距离;Aij为单元i与单元j 间的接触面积。

由热传导导致的单元温度变化，可表示为

由粘性力带来的温升:

式中:m 为单元质量;cv为等容比热;Cn为粘性系数为单元i 与单元j 间的相对速度在其中心连线方向的投影。

本文采用基于无扩散相变的两相模型，该模型基于如下假设:1)每个单元中的相变均匀产生，并由局部能量条件控制;2)单元中的压力和温度满足局部平衡条件。相变模型包括3 个重要的组成部分:1)平衡的相边界;2)局部阈值条件;3)相变动力学方程。

首先，平衡的相边界，即相边界处两相的Gibbs自由能相等Gα(p，T)= Gε(p，T)). 此处i 相的Gibbs 自由能的表达式为

式中:Hi和Si分别为比焓与熵，下标i 分别表示α与ε 相。

在压力不太高的情况下，冲击压缩产生的熵增不大，等熵过程与冲击压缩过程差别不大。因此，对于本文所研究的压力范围内，可采用Murnaghan 状态方程，并采取与Forbes［12］类似的推导方法，得到以下比焓与熵的表达式:

式中:下标0 为初始状态;V0i为零压下的比体积;具体温度及相变模型参数取值见表1.

表1 温度及相变模型参数Tab.1 Values of temperature and phase transition model parameters

其次，局部阈值条件，即当自由能差额ΔG=Gα-Gε超过某一阈值时相变立即被触发。计算中相变和逆相变的激活能分别为ΔGf和ΔGb.

再次，相变动力学方程，用于描述相变份额λE的变化率，本文分别采用1 阶与2 阶动力学方程，具体形式如下:

式中:λE为单元中ε 相的质量份额;1 阶相变动力学方程中λE∈［0，1］，初始时刻取λE=0;2 阶相变动力学方程中λE∈［0.001，0.999］，初始取λE=0.001;τE为单元相变松弛时间，由单元直径与马氏体相变速率之比进行估算，文中相变速率近似取常数1 km/s.

α-ε 相变伴随着约5%的体积变化，数值模拟中通过依据单元中ε 相的质量份额改变单元半径来模拟这一效应:

式中:r0为单元初始半径;V0α、V0ε分别为零压状况下a 和ε 相的比体积。

本文在计算中仅考虑了静水压力而未考虑偏应力对于相变过程的影响。

2 α 铁的数值模拟结果与分析

计算模型飞片与靶板均为6 400 μm×27 μm 的α 铁，单元直径为9 μm，计算单元总数约5 000 个。如图2 所示，计算模型的上、下边界采用周期性边界条件，左、右边界采用自由边界条件。初始温度300 K，飞片不同的速度从左端撞击靶板。

图3 给出了准静态加载条件下铁的相图。从图3 可以看出，本文计算结果与文献［16 -19］实验结果符合较好，证明本文采用的无扩散两相模型能够正确反应α 铁的相变特性。

图2 计算模型示意图Fig.2 Initial model

图3 准静态加载条件下铁的相图Fig.3 Phase diagram of iron under quasi-static loading

数值模拟结果显示室温300 K 下平衡态相边界通过11.15 GPa 压力点。Barker 等［20］通过实验发现铁的冲击相变压力阈值在(12.9 GPa，13.7 GPa)范围内，逆相变压力阈值在(9.4 GPa，10.2 GPa)范围内。因此，文中取相变的激活能为13.40 J/g，对应室温条件下相变的临界压力13.37 GPa;逆相变的激活能取为-8.24 J/g，对应室温条件下逆相变的临界压力为9.80 GPa.

图4(a)给出了Fe 的p-V/V0Hugonoit 曲线;图4(b)给出了Fe 的冲击波波速D 与波后粒子速度up的关系图。由图4 可以看出，本文计算结果与文献［1，20 -21］实验结果符合较好，证明计算所用单元间作用力模型及相变模型能够正确反应铁的冲击压缩及相变特性。

此外，图4(a)的模拟结果显示，冲击相变压力阈值约12.95 GPa，比室温下的结果13.37 GPa 低了约0.42 GPa，这说明冲击导致靶板内温度升高，从而致使相变的压力阈值降低。图4(b)结果表明:当冲击压力小于相变压力阈值时，样品内波形为单波结构;当冲击压力在相变压力阈值至约36 GPa 之间，波形为双波结构(冲击波与相变波);当冲击压力继续增大时，相变波与冲击波汇合成稳定的冲击波。当冲击压力低于相变压力阈值时，冲击波速度随波后粒子速度近似呈线性增长;当冲击压力高于相变压力阈值时冲击波近似以定常速度传播，而相变波波速随波后粒子速度增加较快。

图4 铁的冲击Hugoniot 关系Fig.4 Hugoniot data of Iron

3 多晶铁的数值模拟结果与分析

冲击压缩过程中多晶金属，在晶粒尺度上是一个各向异性、非平衡过程，此时的冲击波结构不能为传统的连续介质力学描述［22］。Meyers 和Carvalho应用晶粒取向的概然论模型研究了多晶镍的冲击响应，结果表明传入多晶中的冲击波前沿变得不规则，其波面不规则度与晶粒尺寸同量级［23］。本文构建了两种不同晶粒尺度的多晶铁模型(如图5)，采用离散元方法研究了晶粒尺度对于冲击波波面不规则程度的影响。

图5 两种不同晶粒尺度多晶铁模型局部放大示意图Fig.5 Enlarged view of initial model of polycrystalline iron

多晶铁模型中的晶粒取向呈随机分布，以此表征晶粒间取向的分布差异。各晶粒的取向为(0° ～60°)之间的随机数，晶粒内部单元为规则的密排六边形，晶界处的间隙由小尺度单元填满，图5 中不同灰度代表不同晶粒取向，两种模型的具体参数参见表2 所示。为避免边界效应的影响，计算模型的上、下边界为周期型边界条件，左侧为固壁边界，右侧为自由边界，初始时刻模型以一定的速度撞击固壁。

表2 两种晶粒尺度多晶铁模型参数Tab.2 Parameters of polycrystalline iron models

图6 为两种不同晶粒大小的多晶铁以700 m/s的速度撞击固璧后，某时刻模型中的冲击波前沿位置图。此处，定义波后粒子速度等于0.5 倍撞击速度的位置为冲击波前沿位置，图中深色部分为冲击压缩区，浅色部分为未受冲击区域。数值模拟结果表明，冲击波前沿分布并不均匀，并且大晶粒模型内冲击波前沿的不规则程度更高。

图7 为冲击波前沿位置中线定义示意图，图中的黑色粗实线代表冲击波前沿位置，阴影部分为波后的冲击压缩区。如图7 所示，沿y 向在冲击波轮廓线的波峰(y 向最大值)与波谷(y 向最小值)之间选取一条直线(平行于x 轴)，使得直线上方的阴影面积之和等于直线下方的空白面积之和，并定义该直线位置为冲击波前沿位置中线。

图6 两种不同晶粒尺度多晶铁中的冲击波前沿位置Fig.6 The shock wave front in polycrystalline iron

图7 冲击波前沿位置中线的定义Fig.7 The definition of shock wave front midline

统计不同时刻两种不同晶粒大小的多晶铁模型中，冲击波前沿位置距离中线位置的标准偏差。图8 给出了冲击波前沿位置的标准偏差随冲击波传播距离的变化关系。从统计结果看，冲击波的波面不规则程度随传播距离的增加而增加;晶粒大小对于冲击波前沿的不规则程度有一定影响，大晶粒模型内冲击波前沿的不规则程度更高。

图8 冲击波前沿位置的标准偏差随冲击波传播距离的变化Fig.8 Variation of standard deviation in the shock wave front position

图9为模型1 以600 m/s 速度撞击固壁后，冲击波前沿附近的质量份额与压力图，图中Ve 代表相变质量份额，p 为压力。从图可以看到:由于晶粒间取向分布差异及晶粒边界效应的影响，压力及质量份额场的分布很不均匀，应力远高于或低于波后平均应力的区域集中于晶界附近;从冲击波前沿位置可见，由于应力集中效应的影响，相变首先发生在晶界处，然后穿透到晶粒内部;在冲击波前沿处的一些晶粒中，可观察到指状传播的相边界。

图9 冲击波前沿附近的质量份额与压力Fig.9 Transformed mass fraction and pressure fields

室温下静压实验的测量结果表明，即使在准静态条件下，铁的α→ε 相变也并非沿平衡面进行［24］。Boettger 等的理论计算结果表明，冲击相变过程中得到的p-λ 曲线与静压实验的测量结果非常相近［25］。本文采用宽度27 μm 的采样窗，沿着波的传播方向统计相变质量份额与局部平均压力。

图10 为采用以上统计方法得到的局部平均压力与相变质量份额曲线。图中的离散点为静压实验的测量结果［24］，两条曲线分别为采用1 阶与2 阶动力学方程得到的统计结果。从图10 中可见，冲击相变的局部压力-相变质量份额曲线均逐渐逼近准静态压缩实验曲线上的点。本文的统计结果与Boettger 等［25］的结果可以互相佐证。采用两种动力学方程得到结果的重要不同点在于相变临界压力，1 阶模型的相变临界压力约为8 GPa，2 阶模型的相变临界压力为10 GPa;与1 阶动力学方程的统计结果相比，2 阶动力学方程得到的局部平均压力-相变质量分额曲线在低压段(小于13 GPa)，与实验结果符合更好。

4 结论

本文利用离散元方法，结合基于无扩散相变的两相模型，模拟了α 铁的冲击相变过程，数值模拟结果表明:

1)计算得到铁的相边界、冲击Hugoniot 关系与实验结果符合较好，证实采用离散元方法结合两相模型模拟α 铁冲击相变过程是可行的。

图10 局部平均压力与相变质量份额Fig.10 Local average pressure and transformed mass fraction

2)对多晶铁的冲击响应过程进行了模拟，初步的数值模拟结果显示冲击波的波面不规则程度随传播距离的增加而增加，大晶粒模型内冲击波前沿的不规则程度更高。

3)对多晶铁冲击相变过程进行了统计，结果表明冲击相变的局部压力-相变质量份额曲线均逐渐逼近准静态压缩实验曲线上的点;与1 阶动力学方程的相比，2 阶动力学方程得到的局部平均压力-相变质量分额曲线在低压段与实验结果符合更好。

References)

［1］Bancroft D，Peterson E L，Minshall S. Polymorphism of iron at high pressure［J］. J Appl Phys，1956，27(3):291 -298.

［2］Cundall P A. A computer model for simulating progressive large scale movement in block rock system［C］∥Symposium ISRM.Nancy，France:International Society of Rock Mechnics，1971:129 -136.

［3］Sawamoto Y，Tsubota H，Kasai Y，et al. Analytical studies on local damage to reinforced concrete structures under impact loading by discrete element method ［J］. Nucl Eng Des，1998，179:157 -177.

［4］刘凯欣，高凌天，郑文刚. 混凝土动态破坏过程的数值模拟［J］. 工程力学，2000(增刊):470 -474.LIU Kai-xin，GAO Ling-tian，ZHENG Wen-gang. Numerical simulation for the concrete failure process under shock loading［J］.Engineering Mechanics，2000(Suppl)470 -474. (in Chinese)

［5］刘凯欣，高凌天. 离散元法在求解三维冲击动力学问题中的应用［J］. 固体力学学报，2004，25(2):181 -185.LIU Kai-xin，GAO Ling-tian. The application of discrete element method in solving three-dimensional impact dynamics problems［J］. Acta Mechanica Solida Sinica，2004，25(2):181 -185.(in Chinese)

［6］Yano K，Horie Y. Discrete-element modeling of shock compression of polycrystalline copper ［J］. Phys Rev B，1999，59:13672 -13680.

［7］Case S，Horie Y. Mesomechanics of the α-δ transition in iron［J］. J Mech Phys Solids，2007，55:589 -614.

［8］王文强. 离散元方法及其在材料和结构力学响应分析中的应用［D］. 合肥. 中国科学技术大学，2000.WANG Wen-qiang. Discrete element method and its use in analysis of response of materials and structures［D］. Hefei:University of Science and Technology of China，2000. (in Chinese)

［9］于继东，王文强，刘仓理，等. 炸药冲击响应的二维细观离散元模拟［J］. 爆炸与冲击，2008，28(6):488 -493.YU Ji-dong，WANG Wen-qiang，LIU Cang-li，et al. Two-dimensional mesoscale discrete element simulation of shock response of explosives［J］. Explosion and Shock Waves，2008，28(6):488 -493. (in Chinese)

［10］Tang Z P，Horie Y，Psakhie S G. Discrete meso-element modeling of shock processes in powders［C］∥High Pressure Shock Compression of Solids Ⅳ，Response of Highly Porous Solid to Shock Loading. N Y:Springer，1997:143 -176.

［11］Yano K，Horie Y. Mesomechanics of the α-ε transition in iron［J］. International Journal of Plasticity，2002，18:1427 -1446.

［12］Forbes J W. Experimental investigation of the kinetics of the shock-induced alpha to epsilon phase transformation in armco iron，WSU-SDL 76-01［R］. Pullman，Washington:Washington State University，1976:112 -120.

［13］Anderson O L. Equations of state of solids for geophysics and ceramic science［M］. New York:Oxford University，1995:195.

［14］Lide D R，Kehiaian H V. CRC handbook of thermophysical and thermodynamical data［M］. US:CRC，1994:136.

［15］Jephcoat A P，Mao H K，Bell P M. The static compression of iron to 78 GPa with rare gas solids as pressure-transmitting media［J］. J Geophys Res，1986，91:4677 -4684.

［16］Giles P M，Longenbach M H，Marder A R. High-pressure α-ε martensitic transformation in iron ［J］. J Appl Phys，1971，42(11):4290 -4295.

［17］Bundy F P. Pressure-temperature phase diagram of iron to 200 kbar，900 ℃［J］. J Appl Phys，1965，36(2):616 -620.

［18］Kennedy G C，Newton R C. Solids under pressure［M］. New York，McGrawHill，1963:163.

［19］Clougherty E V，Kaufman L. High pressure measurement［M］.Washington，Giadini A A，Lloyd E C，1963:152.

［20］Barker L M，Hollenbach R E，Shock wave study of the αTMε phase transition in iron［J］. J Appl Phys，1974，45:4872 -4887.

［21］Brown J M，Fritz J N，Hixson R S. Hugoniot data for iron［J］.J Appl Phys，2000，88:5496 -5498.

［22］Wallace D C. Irreversible thermodynamics of flow in solids［J］.Phys Rev B，1980，22(4):1477 -1486.

［23］Meyers M A ，Carvalho M S. Shock-front irregularities in polycrystalline metals［J］. Mater Sci Eng，1976，24:131 -135.

［24］Taylor R D，Pasternak M P，Jeanloz R. Hysteresis in the high pressure transformation of bcc-to hcp-iron ［J］. J Appl Phys，1991，69(8):6126 -6128.

［25］Boettger J C，Wallace D C. Metastability and dynamics of the shock-induced phase transition in iron［J］. Phys Rev B，1997，55:2840 -2849.