岩层运动并行计算系统StrataKing及应用

王学滨,余保健,马立强,李小帅,张钦杰,杜 轩

(1.辽宁工程技术大学 计算力学研究所,辽宁 阜新 123000;2.辽宁工程技术大学 力学与工程学院,辽宁 阜新 123000;3.中国矿业大学 矿业工程学院,江苏 徐州 221116;4.新疆工程学院 矿业工程与地质学院,新疆 乌鲁木齐 830023)

科学计算是继理论研究和实验研究之后出现的第三种科学研究手段。随着社会的进步,科学计算的地位变得愈发重要,而其他研究手段的重要性则有所降低。科学计算可用于结构优化、工艺筛选和科学现象探索等诸多方面,发挥着至关重要甚至不可替代的作用。

科学计算可用于土木、交通、能源、材料、国防、机械、兵器、航空和航天等诸多工程领域,由此带来了有关产品的研制效率和可靠性大幅提高,有关工程的建设更加经济、安全和高效。在欧美,普遍认为“基于模拟的工程应成为工程与科学领域国家优先发展项目”“从竞争中胜出,就是从计算中胜出”。科学计算需要一定的计算机硬件和软件或程序的支撑。目前,科学计算的作用在一些领域还未充分发挥出来。根源之一是科技人员平时接触到的科学计算仅是在微机上利用国外的串行商业软件进行的。这些计算的规模小(通常仅几十万个单元),效率低,模型较粗糙。随着各界计算机软件或程序开发水平和使用水平的提升,上述情况会逐渐得到改善。

力学分析软件或程序主要基于连续方法、非连续方法和连续-非连续方法。后者兼具前两者的优点(例如,连续方法的计算精度高和非连续方法的适用性广)正在快速发展。在国际上,有限元与离散元耦合方法已取得了诸多的进展[1-4]。本文作者团队历时10余年,自主开发了与有限元与离散元耦合方法相平行的一种新的连续-非连续方法[5-7],即拉格朗日元与离散元耦合方法。以拉格朗日元取代有限元法的目的是该方法适于求解大变形问题,不受负刚度问题困扰。近年来,通过引入GPU并行计算技术,上述方法的计算效率和计算规模得以极大提升,为进一步应用奠定了坚实的基础。该方法的程序已被命名为StrataKing(王之岩层,岩层之王)[8-9]。StrataKing的功能较为强大,主要面向矿业领域,特色较为鲜明。

本文介绍了StrataKing的计算流程、基本原理和功能。通过模拟三点弯梁的开裂过程,在一定程度上检验了StrataKing的正确性。通过模拟静水压力条件下圆形巷道围岩的开裂过程和采动条件下采场岩层的运动过程,呈现了区别于基于传统理论的新颖结果,进一步展示了StrataKing的大规模计算能力、丰富功能和矿业鲜明特色。

1 岩层运动并行计算系统StrataKing

1.1 计算流程

岩层运动并行计算系统StrataKing是在Visual-Studio平台上基于C++语言开发的。在该系统中,采用面向对象的程序设计理念,以更加直观地反映事物属性及事物之间的关系;采用统一计算设备架构(CUDA)对该系统进行并行化处理;采用跨平台的C++图形用户界面应用程序框架Qt对计算结果进行可视化。在该系统中,CPU和GPU各司其职,以发挥CPU擅长控制而GPU擅长繁重计算的各自优势。CPU是主处理器,称为主机端;GPU是协处理器,称为设备端。

图1为该系统的流程图。该流程主要包含前处理部分、计算部分和后处理部分。在上述3个部分执行前,首先应声明头文件和类。

在前处理中,首先,对主机端模型类Model进行实例化,声明主机端模型对象M(以下简称主机端模型,用于存储主机端模型数据);然后,建立主机端模型,其包括单元、节点、内棱(单元的边被称为棱。对于任1条棱,若被2个单元所拥有,则该棱被称为内棱)、外棱(对于任1条棱,若仅被1个单元所拥有,则该棱被称为外棱)、接触点、载荷、约束和锚杆等;最后,声明设备端对象Md,并将主机端模型数据传输至设备端。通过CUDA的内存操作函数cudaMemcpy(),进行主机端与设备端之间的数据传输。

计算部分包含运行在主机端的主体循环函数和运行在设备端的核函数。主体循环函数用于控制计算进程并调用核函数;核函数用于指定各线程所需执行的任务。主体循环的具体流程如下:

(1)当时步数目N=0时,计算开始。当N

(2)每次循环完毕后,需检测是否达到数据或云图显示和保存的间隔。若达到,则通过cudaMemcpy()将数据由设备端拷贝回至主机端,以便于对数据或云图进行输出和可视化。

(3)上述流程被依次执行一次后,N加1,直到N>Nt(总时步数目)时,循环结束。在后处理中,进行计算结果的显示和保存。

1.2 功能

本构模型包括两类:单元内部的和单元之间的。单元内部的本构模型为广义胡克定律。所谓的单元内部的本构模型其实是子单元的。一个四边形单元被剖分成两套三角形子单元,即两种覆盖。四边形单元的应力、应变通过对子单元的应力、应变取平均获得。单元之间的本构模型包括两种:一种是虚拟裂纹模型,另一种是岩层之间界面的黏结模型。前者适于模拟两个相邻单元沿边界开裂的应变软化(甚至脆性),后者可以反映两岩层之间的黏结特性。单元类型通常为正方形单元。

通过删除单元的方式模拟巷道和工作面的开挖。开挖方式包括3种:瞬时开挖、以某一单元为中心向外逐圈删除单元的方式和通过将施加在开挖边界上的力逐渐卸去的方式。对于第2种,卸荷速度通过删除两圈单元的时间步隔来确定。对于第3种,将开挖边界以内的单元删除,与此同时,在开挖边界上施加力等于被删除单元对开挖边界的弹性力,其随后随着时步数目的增加而下降。当应力降至零时,开挖完成。

计算模式包括准静力和动力两种。对于前者,与时间有关的量不是真实的,而是虚拟的,例如,速度、加速度和质量,适于模拟时间较长的过程,例如,煤层开采。此时,更关注最终的结果。对于后者,与时间相关的量是真实的,适于模拟某一时间较短的过程。例如,巷道冲击。

目前,不允许单元内部开裂,仅允许单元沿边界开裂。这样,若干原本离散的单元或从母体中脱离出来的若干单元堆积在一起后,可以承受过高的应力,这与实际情况不符。为此,需要对这些单元的应力进行处理。应力跌落方式包括两种:球应力不变和围压不变。

2 计算模型及结果分析

2.1 计算模型

共计算了4个算例,计算条件均为平面应变、大变形,计算模型的示意图见图2,部分参数见表1。

图2 计算模型 Fig.2 Computational models

算例1用于模拟三点弯梁的开裂过程。模型被划分为100×30个正方形单元。在梁上边界跨中节点上施加向下的速度v,不考虑重力作用。计算模式为准静力。应力跌落方式为球应力不变。法向黏聚力σn和法向张开度w之间的关系为指数形式:

(1)

式中:σt为抗拉强度,Pa;G为Ⅰ型断裂能,G=100 N/m。

算例2—3用于模拟静水压力条件下圆形巷道的开裂过程,单元数分别为25万和400万。对于算例2—3,巷道开挖之前的模型分别被划分为500×500和2 000×2 000个正方形单元。不考虑重力作用。计算模式为准静力。应力跌落方式为球应力不变。计算过程包括下列3步:首先,在压应力作用下,将模型计算至平衡,算例2—3分别消耗6 248和25 000个时步;其次,在模型中部以将施加在开挖边界上的力逐渐卸去的方式开挖巷道(半径为0.8 m),算例2—3分别消耗20 833和83 333个时步;最后,在巷道开挖完毕后,继续计算。

算例4用于模拟采动条件下采场的岩层运动过程。工作面开采之前的模型被划分为800×233个正方形单元。计算模式为准静力。应力跌落方式为球应力不变。计算过程包括下列3步:首先,在压应力作用下,将模型计算至平衡,消耗10 000个时步;其次,从左向右开采煤层。采高为10 m,每间隔4 000个时步瞬时开采1 m,共开采200 m;最后,在煤层开采后,继续计算。

2.2 结果分析

2.2.1 三点弯梁

图3给出了算例1中三点弯梁的载荷-位移曲线(图3(a))和以最大主应力σ3表征的拉裂过程(图3(b)—(e)),其中,黑色、灰色线段分别代表剪裂纹、拉裂纹区段。图3(a)中的4点a—d分别对应图3(b)—(e)。位移是梁上边界跨中节点的。由此可以发现,梁跨中横截面下部的σ3集中达到一定程度后,裂纹垂直向上发展,这与常识相一致。当载荷达到峰值时,梁已发生少许开裂。在峰值稍前,载荷-位移曲线呈上凸特征。当载荷达到峰值之后,载荷-位移曲线先呈上凸后呈下凸特征。上述数值结果能与文献[10]的数值结果基本吻合,由此可在一定程度上检验了StrataKing程序的正确性。

图3 三点弯梁的计算结果Fig.3 Results of the three-point bending beam

2.2.2 静水压力条件下圆形巷道围岩

图4(a)—(d)给出了算例2中以最大主应力σ3表征的巷道围岩的开裂过程。由此可以发现:

图4 粗网格条件下圆形巷道围岩的计算结果Fig.4 Results of the circular tunnel surrounding rock for the coarse mesh

当N=21 000时(图4(a),卸荷进行中),巷道周边的σ3较高,有些位置甚至出现了σ3>0。σ3高值区的轮廓大致成正方形。

当N=30 000时(图4(b),卸荷已完成),巷道的左上、左下、右上和右下4个位置存在多条裂纹,有些裂纹有形成V形坑的趋势。

当N=50 000～200 000(图4(c)—(d)),巷道周边的裂纹不断向深部扩展,最终停止发展。应当指出,巷道周边的裂纹分布并不均匀,相比之下,巷道左帮的裂纹较少。

图5(a)—(d)给出了算例3中以最大主应力σ3表征的巷道围岩的开裂过程。由此可以发现:

图5 细网格条件下圆形巷道围岩的计算结果Fig.5 Results of the circular tunnel surrounding rock for the fine mesh

当N=66 000时(图5(a),卸荷进行中),巷道左上、左下、右上和右下4个位置存在4个小V形坑,这与这些位置存在一定的剪应力有关。巷道表面与小V形坑相交处存在若干较短的裂纹。

当N=92 000(图5(b),卸荷进行中),巷道周边存在众多较长的裂纹。巷道拱顶、拱底和两帮的若干较长的裂纹形成了4个大V形坑。其余位置的裂纹不如上述4处发育。

当N=120 000～1 100 000时(图5(c)—(d),卸荷已完成),巷道周边的裂纹不断向围岩深部发展,最终仍在发展。应当指出,巷道周边的裂纹分布较算例2中的更加均匀,裂纹以更加弯曲、柔顺而非突兀、生硬的方式向围岩深部发展。显然,算例3的结果更加合理。

图4(e)和图5(e)分别给出了算例2—3中Ns和Nt随着N的演化规律。其中,黑色、灰色线段分别代表剪裂纹、拉裂纹区段。由此可以发现,随着N的增加,二者的增速减缓,直至增速为零,这代表巷道围岩中的裂纹停止发展,此时,巷道围岩处于平衡状态。相比之下,剪裂纹的发展在先,Ns远大于Nt。

2.2.3 采场

图6(a)—(d)给出了算例4中以最大主应力σ3表征的岩层的运动过程。由此可以发现:

图6 采场岩层运动的计算结果Fig.6 Results of the strata motion of a stope

当N=180 000(图6(a),推进距离L=42 m),采空区前、后煤壁存在半圆形的σ3高值区。除上述区域之外,随着远离煤壁,煤层的σ3的值先增大后减小至定值。在采空区上方,顶煤发生开裂,并向下运动;泥质砂岩层1受拉严重;细砂岩层存在一个“V”形的σ3受拉区。

当N=450 000(图6(b),L=110 m),采空区两侧的细砂岩层1存在σ3低值区,这代表这部分岩石受挤压作用强烈。在采空区上方,若干较薄岩层已冒落,细砂岩层1中部存在向上发展的狭长拉裂纹。

当N=810 000(图6(c),L=200 m),采空区上方粗砂岩层以下的岩层已冒落。采空区两侧存在大量拉裂纹和少量剪裂纹。最明显的离层发生在粗砂岩层和砂岩层之间,其上方的诸多岩层发生了明显的弯曲下沉现象。最高位离层位于细砂岩层2和泥质砂岩层2之间,距煤层上表面的距离约为60 m。

当N=1 200 000(图6(d),L=200 m),采空区上方的岩层进一步冒落,导致大部分采空区闭合;若干离层闭合;采空区上方粗砂岩层及其上方的若干岩层发生了比过去更加明显的弯曲下沉现象。与此同时,在采空区两侧,粗砂岩层至泥质砂岩层3之间的拉裂纹和剪裂纹进一步发展。

上述结果不仅呈现了采空区上方的上底短下底长的梯形的“三带”(冒落带、裂隙带和弯曲变形带),还呈现了更大范围的上底(呈马鞍形)长下底短的梯形裂纹分布形态,这与文献[12-13]中的导水裂隙带定性相符。

图6(e)—(f)分别给出了工作面煤壁超前支承压力和水平应力的时空分布规律。其中,支承压力为采高范围内煤层中部单元的垂直应力的绝对值,水平应力也取了绝对值。由此可以发现,随着L的增加,支承压力的峰值有增加的趋势,这与文献[14]的数值结果有类似之处;水平应力也具有峰值,其随着L的增加呈增加的趋势。水平应力具有峰值的现象与岩石力学中厚壁筒的径向应力解答(随着远离筒壁,径向应力的值增大,直至达到远场应力)不符。这是因为在本文模型中,两相邻岩层之间存在界面黏聚力,这不同于连续介质模型。水平应力的峰值的提升现象亦可用于解释支承压力的峰值随着L的增加而增加的现象。

3 结论

1)在岩层运动并行计算系统StrataKing中,CPU和GPU各司其职,以发挥CPU擅长控制而GPU擅长繁重计算的各自优势。应力-应变模块等计算模块均在GPU上执行。单元内部的本构模型为广义胡克定律。单元之间的本构模型包括两种:一种是虚拟裂纹模型,另一种是岩层之间界面的黏结模型。通过模拟三点弯梁的开裂过程,StrataKing的正确性得到了一定程度上的检验。

2)在静水压力条件下,当单元数目较多时,巷道围岩的裂纹以更加弯曲、柔美而非突兀、生硬状向围岩深部发展,这与塑性力学中圆环形的各处均破坏的塑性区有所不同。在采动条件下,工作面煤壁超前支承压力和水平应力均具有峰值,后者与岩层之间的黏结作用有关,这与基于连续介质假定的厚壁筒径向应力特征存在不同。