多晶铁冲击相变的离散元方法研究

多晶铁冲击相变的离散元方法研究

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

（北京应用物理与计算数学研究所，北京，8009信箱，100088）

摘要 采用离散元方法结合基于无扩散相变的两相模型、热力学相容的Helmoholtz自由能函数与有限速率相变动力学方程，模拟了α铁的冲击相变过程，计算得到了铁的相边界及冲击Hugoniot关系。对于多晶铁的冲击相变过程进行了模拟，研究了冲击波前沿不规则程度随冲击波传播距离及晶粒大小的变化规律，对于多晶冲击相变过程进行了统计，得到了相变特征量随局部冲击压力的变化规律。

关键词 冲击相变；离散元方法；多晶材料

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

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

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

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

1

计算方法与相变模型

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

1

1.1 单元间相互作用力模型

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

图1. 单元间相互作用力的示意图

Fig.1 Interaction model of an element-pair

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

F=P+fvt+fvn+fd+fsh (1.1)

其中，P为中心势力，fvt为切向摩擦力，fvn为中心阻尼，fd为干摩擦力，fsh为弹塑性剪切力。单元间相互作用力的形式、单元间作用状态判定方法及等效应力与应变的计算等参见文献[10]。

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

1.2 温度与相变模型

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

ijijijijijijijijijijijT=Trev+∑ΔTirrt>0 (1.2)

利用（1.10）式，可以得到温度的可逆部分，即等熵过程中温度可表示为：

Trev1⎧⎡⎤⎫β⎛⎞P⎪⎪=T0exp⎨γ0⎢⎜+1⎟−1⎥⎬⎢⎥⎪⎢⎝B0⎠⎪⎥⎦⎭ (1.3)

⎩⎣其中，T0为初始温度；P为静水压力；B0与等温体模量KT相关B0=KT态方程参数；γ0为Gruneisen参数。

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

ΔTirr=ΔTcond+ΔTvis (1.4)

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

Ti−TjijAΔtΔQ=−κijd (1.5)

ΔQ为在Δt时间内由单元j向单元i传导的热，κ为热传导系数，Ti和Tj分别为单元iijij与单元j的温度，d为单元i与单元j间的距离，A为单元i与单元j间的接触面积。

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

2

ΔTcond=由粘性力带来的温升：

ΔQMCv (1.6)

ijij1(Cnvn)⋅vnΔTvis=ΔtMC2v

(1.7)

ij其中，M为单元质量，Cv为等容比热，Cn为粘性系数，vn为单元i与单元j间的相对速度在其中心连线方向的投影。

本文采用基于无扩散相变的两相模型，该模型基于如下假设：

（1） 每个单元中的相变均匀产生，并由局部能量条件控制；

（2） 单元中的压力和温度满足局部平衡条件。

相变模型包括三个重要的组成部分：

（1）平衡的相边界；（2）局部阈值条件；（3）相变动力学方程。

首先，平衡的相边界，即相边界处两相的Gibbs自由能相等（Gα(P,T)=Gε(P,T)）。此处i相的Gibbs自由能的表达式为：

Gi(P,T)=Hi(P,T)−TSi(P,T)其中，H和S分别为比焓与熵。

(1.8)

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

1⎧⎫β⎞i⎪⎪⎛PHi(P,T)=H0i+Cvi(T−T0)−γ0iCviT0⎨⎜+1⎟−1⎬⎠⎪⎝B0i⎪⎩⎭11−⎧⎫⎞βi⎪βiV0iB0i⎪⎛P++1⎟−1⎬⎨⎜βi−1⎪⎝B0i⎠⎪⎩⎭ (1.9)

1⎧⎫βi⎛T⎞⎛⎞P⎪⎪Si(P,T)=S0i+Cviln⎜⎟−γ0iCvi⎨⎜+1⎟−1⎬⎝T0⎠⎠⎪⎝B0i⎪⎩⎭ (1.10)

此处，下标“i”的取值为1或2，分别表示α与ε相；下标“0”为初始状态；V0i为零压下的比体积；具体温度及相变模型参数取值见表1。

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

第三，相变动力学方程，用于描述相变份额λE的变化率，本文分别采用一阶与二阶动力学方程，具体形式如下：

dλE1−λEλE(1−λE)=or (1.11)

dtτEτE式中，λE为单元中ε相的质量份额，一阶相变动力学方程中λE∈[0.0,1.0]，初始时刻取 3

λE=0.0；二阶相变动力学方程中λE∈[0.001,0.999]，初始取λE=0.001；τE为单元相变松弛时间，由单元直径与马氏体相变速率之比进行估算，文中相变速率近似取常数1kms。

表1、 温度及相变模型参数

Table.1 Values of temperature and phase transition model parameters

ρ(gcm3)

γ0

Cv(J/g/K)H0(J/g)S0(J/g/K)KT(GPa)

β

α相 7.87[14] 1.68[11] 0.4414[11] 79.3[12] 0.489[15] 166[11] 4.25[11]ε相 8.28[14] 2.24[11] 0.4460[11] 169a 0.539[11] 139.6[13] 4.4[13]a、根据P=11.15GPa，T=300K条件下，两相自由能相等（即Gα=Gε），确定的参数。

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

⎛V⎞r=r0⎜1−λE+λE0α⎟V0ε⎠⎝ （1.12）

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

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

−122

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

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

图2 计算模型示意图

Fig.2 Initial model sketch

图3给出了准静态加载条件下铁的相图，图中菱形为Gilen等人的实验结果[16]，圆圈为Bundy 等人的实验结果[17]，方块为 kenney等人的实验结果[18]，三角型为Clogherty等人的实验结果[19]，实线为本文的计算结果，计算结果与实验结果符合较好，证明本文采用的无扩散两相模型能够正确反应α铁的相变特性。

数值模拟结果显示室温(300K)下平衡态相边界通过11.15GPa压力点。Barker等[20]通过实验发现铁的冲击相变压力阈值在(12.9GPa, 13.7GPa)范围内，逆相变压力阈值在(9.4GPa,

10.2GPa)范围内。因此，文中取相变的激活能为13.4J/g，对应室温条件下相变的临界压力13.37GPa；逆相变的激活能取为-8.24J/g，对应室温条件下逆相变的临界压力为9.8GPa。

4

图3 铁的相图

Fig.3 Phase diagram of iron

图4(a)给出了Fe的P-V/V0Hugonoit曲线；图4(b) 给出了Fe的冲击波波速D与波后粒子速度up的关系图。图4(a)中实心方块为Bancroft等人的实验数据[1]，空心圆圈为Barker等人的实验数据[20]，实线为本文计算结果。图4(b)中实心方块为Barker等人的实验数据[20]，实心三角型为Brown等人的实验数据[21]，空心圆圈为本文计算结果。计算结果与实验结果符合较好，证明计算所用单元间作用力模型及相变模型能够正确反应铁的冲击压缩及相变特性。

5

(a) 铁的P-V/V0 Hugoniot关系 (b) 铁的D-up关系

(a) P-V/V0

Hugonoit of iron (b) D-upHugonoit of iron

图4. 铁的冲击Hugoniot关系

Fig.4 Hugoniot data of Iron

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

3

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

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

图5 两种不同晶粒尺度多晶铁模型局部放大图

Fig.5 An enlarged view of initial model of polycrystalline iron

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

表2、两种晶粒尺度多晶铁模型参数

Table.2 parameter of polycrystalline iron model

模型

模型尺寸μm2

晶粒总数1

2

晶粒平均尺寸μm

136.0×136.0

67.5×67.5

2单元直径μm

9

9

63240×540

3240×540

94

384

图6 两种不同晶粒尺度多晶铁中的冲击波前沿位置

Fig.6 The shock wave front in polycrystalline iron

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

图7 冲击波前沿位置中线的定义

Fig.7 The definition of shock wave front midline

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

图8 冲击波前沿位置的标准偏差随冲击波传播距离的变化

Fig.8 Variation of standard deviation in the shock wave front position

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

7

图9 冲击波前沿附近的质量分额与压力

Fig.9 Tansformaed mass fraction and pressure fields

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

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

图10 局部平均压力与相变质量份额

Fig.10 Local average pressure and transformed mass fraction

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

13GPa）与实验结果符合更好。

4

结论

数本文利用离散元方法结合基于无扩散相变的两相模型，模拟了α铁的冲击相变过程，值模拟结果表明：（1）计算得到铁的相边界、冲击Hugoniot关系实验结果符合较好，证实采用离散元方法结合两相模型模拟α铁冲击相变过程是可行的。（2）对于多晶铁的冲击响应过程进行了模拟，初步的数值模拟结果显示冲击波的波面不规则程度随传播距离的增加而增加，大晶粒模型内冲击波前沿的不规则程度更高。（3）对于多晶铁冲击相变过程进行了统计，结果表明冲击相变的局部压力-相变质量份额曲线均逐渐逼近准静态压缩实验曲线上的点；与一阶动力学方程的相比，二阶动力学方程得到的局部平均压力-相变质量分额曲线与实验结果符合更好。

参 考 文 献

1 Bancroft D, Peterson E L, Minshall S 1956 J. Appl. Phys. 27:291.

2 P.A. Cundall, in Proc. Symp. Int. Soc. Rock Mech., (1971), II, Art. 8.

3 Y. Sawamoto, H. Tsubota, Y. Kasai, et al. Nuclear Engineering and Design. 1998, 179:157–177.

4 刘凯欣,高凌天，郑文刚，工程力学（增刊）2000：470-474.

5 刘凯欣,高凌天，固体力学学报2004, 25(2): 181-185.

6 K. Yano, Y. Horie. Phys Rev B, 1999, 59: 13672-03680.

7 S Case, Y. Horie. J. Mech. Phys. Solids2007,55: 589–614.

8 王文强，博士学位论文，中国科学技术大学，2000.

9 于继东,王文强，爆炸与冲击2008，28(6):488–493.

10 Tang Z P, Horie Y, Psakhie S G. High Pressure Shock Compression of Solids IV, 1997, 143–176.

11 Yano K, Horie Y. Mesomechanics of the α–ε Transition in Iron. 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. NSWC/WOL TR 77-137.

13 Jephcoat A P, Mao H K, Bell P M. The static compression of iron to 78GPa with rare gas solids as

pressure-transmitting media. J. Geophys. Res. 1986, 91: 4677 – 4684.

14 Anderson O L. Equations of State of Solids for Geophysics and Ceramic Science. Oxford University, New

York. , 1995, P195.

15 Lide D R, Kehiaian H V. CRC Handbook of Thermophysical and Thermodynamical Data. CRC, Ann Arbor.,

1994, P136.

16 Giles P M, Longenbach M H, Marder A R. High-Pressure α–ε Martensitic Transformation in Iron. J.

Appl. Phys., 1971, 42(11): 4290–4295.

17 Bundy F P. Pressure-temperature phase diagram of iron to 200 kbar, 900 _C. J. Appl. Phys., 1965, 36(2): 616–620.

18 Kennedy G C, Newton R C. Solids Under Pressure. Editted by Paul W, Warschauer D M, P163.

19 Clougherty E V, Kaufman L. High Pressure Measurement. Edited by Giadini A A, Lloyd (Butterworths,

Wshington, D. C.), P152.

20 Barker L M, Hollenbach R E, Shock wave study of the α