岩石是一种天然的损伤材料,其内部含有大量的微裂纹、微缺陷,这些天然缺陷在外力的作用下逐渐劣化并最终引起材料的破坏。岩石在冲击荷载作用下性质逐渐劣化并最终引起材料破坏的过程称为冲击损伤演化。冲击演化规律的研究在岩石损伤力学的研究中占有十分重要的地位。

数值模拟早已成为目前科学研究中的一个重要组成部分,且由于其具有明显的可重复性和低成本,在科学研究领域中发挥着越来越重要的作用。同样,在冲击损伤演化规律的研究中,数值模拟也占有很重要的地位。目前也有众多学者利用数值方法对载荷作用下的损伤演化规律进行模拟[1 ,2] ,但所采用的数值分析方法基本都是有限单元法;笔者为了探索其他数值分析方法在损伤演化计算中的应用,试采用一种较新的数值分析方法———数值流形方法来对冲击损伤演化规律进行模拟,并开发了相应的计算子程序。

数值流形方法是石根华博士继提出不连续变形分析方法(DDA)之后,又提出的一种新的数值分析方法[3 ,4] 。由于具有明显的优越性,所以一开始他就吸引了很多人的注意。到目前为止,该方法已在岩土工程分析中得到了较为广泛的应用[5] 。由于提出的时间比较短,应用范围相对来说还比较狭窄,不象有限元和离散元等在解决实际问题中应用的那么普遍。但是由于其在计算的过程中是采用两套网格———物理网格和数学网格,因此在计算材料变形过程中裂纹的产生、扩展等方面具有明显的优越性,因而也在这类问题中得到了广泛的应用。笔者在早期开发的仅考虑线弹性本构关系的数值流形方法程序的基础上加入以应力波衰减为基础的冲击损伤模型,并利用扩展后的程序实现了对岩石长杆在冲击载荷作用下的损伤演化过程的模拟。

《1   岩石冲击损伤应力波衰减模型》

1   岩石冲击损伤应力波衰减模型

《1.1 损伤变量的定义》

1.1 损伤变量的定义

岩石中损伤的存在和发展会引起应力波波速的衰减,所以 Robin 和 Ahrens 采用了如下的损伤参量表达式[8]

为受载损伤后岩石的纵波波速, 为受载前岩石的纵波波速。由于(2 是受载损伤前后岩石弹性模量之比,所以 D 表述的是岩石的相关模量的变化。

1.2 损伤演化方程及本构方程[6]

在冲击载荷的作用下岩石的破坏主要与微裂纹的拉伸激活有关。岩石冲击实验和超声波实验研究揭示了应力波在损伤岩石中传播时的衰减特性,即随着损伤的不断增加,应力波在其内传播的衰减愈大,表现为衰减系数的不断增加,损伤耗散能增加,因此可以认为衰减系数 α 和损伤耗散能 YD 之间的关系式为:

式(2)反映了损伤演化的基本规律。其中 Kα 是拟合常数,α0 可以认为是岩石介质的初始衰减系数,随着应力波的传播和损伤的演化,介质的衰减系数增加,其损伤程度越来越大。损伤参数与衰减系数间的关系为线性关系:

将式(2)和式(3)写成导数形式:

式中 C 是材料常量(C = 1/B)。式(4)和式(5)为体积拉伸下的损伤演化方程。由于岩石在压缩状态下,损伤演化一般很小,所以一般不考虑岩石在压缩状态下的损伤累计。

损伤材料的本构关系为:

将上式写成偏量部分和体积部分的率形式为:

式中, KG 为损伤岩石的体积模量和剪切模量, D 为损伤参数, 为应变偏量张量, 为单位张量。式(4),式(5)和式(7)构成了岩石应力波衰减损伤模型的封闭方程。该模型具有输入参数少、实验参数易得的优点,所以可以考虑将该模型加入到数值流形方法的计算程序中与弹性本构关系一起构成复合本构关系,来对冲击问题进行模拟研究。这里采用的弹性本构关系也应该是考虑损伤后的弹性本构关系,即:

《2   冲击损伤演化的数值流形方法模拟》

2   冲击损伤演化的数值流形方法模拟

《2.1 数值流形方法计算程序功能拓展》

2.1 数值流形方法计算程序功能拓展

由于数值流形方法原程序,仅适合于计算静载作用下的弹性变形问题,就所研究的问题而言,需要改进的地方主要集中在两个方面:

1)原程序中作用于材料的载荷是恒定的静载,而这与冲击载荷特性有着明显的不同,因此载荷的作用形式需要扩展。笔者将冲击载荷简化为三角波载荷的形式,即有明显的起止作用时间和应力峰值。

2)原程序中材料的应力应变关系仅仅采用了弹性本构关系;笔者将其修改为在不同的载荷作用阶段和作用情况下采用不同的本构关系,即采用弹性本构关系和应力波衰减模型相结合的复合本构关系。当块体单元处于拉伸状态并且应变变化率大于给定值(取应变率临界值为1 000 ,即认为当应变率大于1 000时,按动态问题处理,否则按静态问题处理[8] )时,采用应力波衰减模型;而应变变化率小于给定值时和当块体单元处于压缩状态时采用弹性本构关系。

《2.2 程序计算流程》

2.2 程序计算流程

将上述冲击损伤本构模型加入到改进的数值流形方法程序中,可实现对冲击损伤模型的数值模拟,详细的算法流程如下:

(a)若一单元损伤变量 D = 1 ,则不再对该单元进行计算而直接进入下一单元的计算,否则进入(b);

(b)计算体积应变 ,如果 < 0 或 = 0 则进入(c),否则进入(d);

(c)按照弹性模型,由应变计算应力,并且保持 D 值不变(不考虑体积压缩时的损伤积累);

(d)计算体积应变变化率痹ε,并计算新的 D 值,如果 > 阀值,则进入(e),否则进入(f);

(e)按照动态(应力波衰减模型)由应变计算应力值,进入(g);

(f)按照静态(弹性模型)由应变计算应力值,进入(g);

(g)如果 D > 1 或 D = 1 ,则令 D = 1 ,并置当前所计算块体的压力和偏应力为零,进入下一步的计算。

在每步的计算中需要对所有单元都进行计算,并保存当前步的应力、应变计算结果以供下步计算使用。

《2.3 应力波冲击损伤演化过程模拟》

2.3 应力波冲击损伤演化过程模拟

材料在冲击载荷作用下的破坏,是由于应力波在材料中传播时,由于波的反射和叠加,使在材料中传播的应力波强度超过了材料的强度极限而导致的,所以研究应力波在材料中的传播过程具有重要的意义。同时在岩石爆破等工程实际中,应力波在破岩过程中也起着非常重要的作用,然而目前的理论都认为岩石的破坏过程是爆炸应力波与爆生气体共同作用的结果,因而为了能够研究应力波在岩石爆破过程中对破岩的贡献,目前也大都采用霍布金森杆等可以忽略爆生气体的作用而仅考虑应力波作用下材料的破坏过程来研究,因此在本文中,对应力波在一岩石长杆中的传播过程进行模拟研究。

2.3.1   计算模型与计算参数计算模型为如图 1所示,模型长为 6 m ,宽为 0.3 m (长细比为 20),整个杆是一个自由杆,不受任何约束,在杆的左端面受一均匀分布的冲击载荷 p 的作用,其中 bcde 分别为位于杆中轴线上的左端面、中部、右端面及 1/4 处和 3/4 处的 5 个测点。 p 为作用在左端面的三角波冲击载荷,其参数为:载荷峰值时间 t1 = 80 μs ,载荷结束时间 t2 = 300 μs ,载荷峰值 Fmax = 2.9 × 107 N 。由于在冲击损伤模型中,还要用到材料本身的材料参数和应力波模型中的衰减参数,而这些参数都需要通过相应的实验才能测定,所以笔者根据文献[7]中相关试验采用表 1 所示的计算参数。

《图 1》

图 1 计算模型

Fig.1 Calculation model

《表 1》

表 1 模型计算参数

Table 1 Calculation parameters of the model

2.3.2   损伤演化模拟结果与分析在损伤模型的研究中,材料损伤演化规律的研究占有着重要的地位,他反映了材料在外力作用下其性质的逐渐劣化过程。由于数值流形方法也是通过划分单元来进行计算的,所以对于不同的单元其损伤是不同的,也就是说在外力的作用下,材料内部的损伤将呈现一个渐变的趋势。为了更清楚地理解材料的损伤情况,笔者采用了对材料内部的不同位置处的损伤情况和整个材料的损伤演化趋势两种情况分析。

1)不同位置点处的损伤演化分析在距离载荷点的不同位置处,材料的损伤程度是不同的,一般而言,在距离载荷点比较近的位置处,由于作用力的强度很大,很容易对材料造成严重的损伤,而随着距作用点距离的增加,应力波的强度不断减少,因此也就渐渐地导致了材料内部相应点处损伤程度的缓解。所以对材料内部各不同位置点处损伤程度的分析,对把握材料内部不同位置点处的受力情况和材料性质的劣化情况也有着很重要的意义。下面就对上例中不同位置点处的损伤情况进行分析。由于不同点处的损伤程度和该点的应力波强度是密切相关的,所以很有必要对各点在不同时刻的应力强度进行分析,但是基于侧重点和篇幅的限制,在此仅给出应力波第一次到达各点时的应力强度值,按照模型中的顺序,dbec 5 点处的应力强度峰值分别为: 94.0,67.0,58.0,52.0 和 3.8 。(MPa)

在本算例中,材料的损伤是用应力波波速来定义的,而由弹性理论可知,以应力波波速定义的损伤变量和以弹性模量变化来定义的损伤变量是一致的,所以可把材料的损伤定义为:

其中: Ei 为第 i 个单元在初始状态时的弹性模量; 为第 i 个单元在某一损伤状态下的弹性模量; i 为单元号。根据该式的定义可求出材料内部不同点处的损伤演化曲线如图 2 所示。

《图 2》

图 2 ~ e 点处的损伤演化曲线图

Fig.2 Damage evolution curves of bcde points

从这 5 个点处的损伤演化曲线图可以看出,在不同位置点处对材料的损伤是不同的。首先从损伤的严重程度来说,依次是: 点> d 点> b 点> e 点> c 点,与载荷作用点距离由近至远的规律,也就是说距载荷作用点越近,损伤就越严重。随着距载荷点距离的增加,应力波在传播过程中就会发生衰减,造成应力峰值下降,强度减弱,因而对材料的损伤程度也就相应减小。其次,从损伤开始增加的时间来说,依次是: 点 < d 点 < b 点 < e 点 < c 点,这对于 点来说,当载荷作用到该点上以后,其损伤几乎是立刻就开始增加,并迅速达到最大值,而对 d 点来说,损伤开始增加的时间就稍稍向后推移了一点,这是因为应力波的传播是需要一个过程的,只有当应力波从点 传到 d 点的时候, d 点损伤才开始增加,这对于其他几个点来说也是同样的道理,所以损伤开始增加的时间是依次延长的。第三,从损伤增加的速率来说,依次是: 点 > d 点 > b 点 > e 点 > c 点。这一点由图 2 可以明显看出,这些点的损伤由初始损伤增加到曲线上的一个极值时,所需的时间基本上都是相等的,但是损伤的峰值却相差很大,对于 点来说,它是由 0.2 增加到了 0.412 ;对于 d 点来说,它是由 0.2 增加到了 0.306 ;对于 b 点来说,它是由 0.2 增加到了 0.270 ;对于 e 点来说,它是由 0.2 增加到了 0.255 ;对于 c 点来说,它是由 0.2 增加到了 0.200 3 。所以明显可以看到损伤的增加速率相差是很大的。而从应力波传播的角度来说,应力波的强度峰值在开始点处也是很强的,然而随着应力波传播距离的增加,其峰值降低,峰值的上升坡度也有所缓和。第四,从各点的损伤随时间的变化曲线还可以看出,材料的损伤由初始值增加到最终值并不是一次完全完成的,而是经历了几个相应的周期,这说明当第一次应力波传播过后,随后反射回来的应力波又继续对材料作用,而且其强度也较大并对材料造成了相应的损伤。最后,从材料的损伤演化曲线中还可以看出,在应力波反复对材料作用并引起材料损伤的过程中,每两次损伤增加的中间阶段都有一段平直的线段,而且平直线段是由短逐渐变长的,这也充分说明了应力传播的速度是在衰减的。因为点的位置是固定的,即应力波传播的距离是不变的,但从由前到后的几个不同周期中,应力波从该点开始传播到又传播回该点的一个周期时间是逐渐增加的,这充分说明应力波的传播速度是下降的,这也从另一方面说明了应力波衰减损伤模型的正确性。

2)材料的总体损伤演化分析通过上面的讨论可以知道,材料内部每一点处的损伤随时间的发展显现出相同的规律。但是利用这种方法也是仅仅只能给出有限个点的损伤演化曲线,很难反映出整个材料在外力作用下的损伤情况。因此,下面就从宏观上对材料的总体损伤进行计算,从整上来把握材料的损伤情况。材料的总体损伤应该定义为:在某一时刻所有单元在初始状态时的弹性模量之和减去所有单元在该时刻的经历损伤后的弹性模量之和再除以所有单元在初始状态时的弹性模量之和,得到的这个无量纲值就是材料在该时刻的损伤值。用公式表示为:

 为所有单元在初始状态时的弹性模量之和; 为所有单元在某一损伤状态下的弹性模量之和; n 为单元总数。据式(10)的定义,就可以求得材料在任一时刻的损伤值,把这些不同时刻的损伤值描绘出来,就可以得到在上述冲击载荷作用下材料的损伤演化曲线,其结果如图 3 所示。从该损伤演化曲线中可以看出在最初的 0.003 s 内,损伤由初始损伤 0.2 增加到 0.260 ,损伤增加幅度很大,而在以后的 0.007 s 内,损伤的增加就很小,仅从 0.260 增加到 0.274 。这主要是由于当应力波在材料中传播时,在初始的阶段,应力波强度很大很容易造成材料的破坏而产生新的损伤,而在后期由于应力波衰减很快,其强度被大大削弱,所以很难再使材料的性质发生劣化而产生新的损伤,所以后期材料损伤值增加很小,这一点与理论分析是很吻合的。同时从损伤的演化曲线中可以看出,损伤曲线是凸形的,这也说明损伤在开始阶段增加很快,而随着时间的推移,损伤增加的越来越小,有可能存在着一个损伤的极限值,即随着损伤的不断演化,应力波峰值降低,到某一时刻,降低了的应力波峰值将有可能不会对材料造成新的损伤,这时损伤将趋于一个稳定的值,所以从这点上来说,从该模型中得出的损伤演化曲线是很符合实际的物理过程。

《图 3》

图 3 整体损伤演化曲线图

Fig.3 Damage evolution curve of total mode

《3   结语》

3   结语

利用扩展后的流形程序很好地模拟岩石在冲击载荷作用下的内部损伤演化过程。损伤演化曲线很好地再现了材料在外力作用下的损伤演化过程,并与理论预测结果十分吻合。这充分说明了把数值流形方法引入岩石冲击损伤计算中是可行的,这有望对目前的岩石冲击损伤的数值模拟研究开创出一条新的途径。