巖層運動并行計算系統StrataKing及應用

王學濱,余保健,馬立強,李小帥,張欽杰,杜 軒

(1.遼寧工程技術大學 計算力學研究所,遼寧 阜新 123000;2.遼寧工程技術大學 力學與工程學院,遼寧 阜新 123000;3.中國礦業大學 礦業工程學院,江蘇 徐州 221116;4.新疆工程學院 礦業工程與地質學院,新疆 烏魯木齊 830023)

科學計算是繼理論研究和實驗研究之后出現的第三種科學研究手段。隨著社會的進步,科學計算的地位變得愈發重要,而其他研究手段的重要性則有所降低?？茖W計算可用于結構優化、工藝篩選和科學現象探索等諸多方面,發揮著至關重要甚至不可替代的作用。

科學計算可用于土木、交通、能源、材料、國防、機械、兵器、航空和航天等諸多工程領域,由此帶來了有關產品的研制效率和可靠性大幅提高,有關工程的建設更加經濟、安全和高效。在歐美,普遍認為“基于模擬的工程應成為工程與科學領域國家優先發展項目”“從競爭中勝出,就是從計算中勝出”?？茖W計算需要一定的計算機硬件和軟件或程序的支撐。目前,科學計算的作用在一些領域還未充分發揮出來。根源之一是科技人員平時接觸到的科學計算僅是在微機上利用國外的串行商業軟件進行的。這些計算的規模小(通常僅幾十萬個單元),效率低,模型較粗糙。隨著各界計算機軟件或程序開發水平和使用水平的提升,上述情況會逐漸得到改善。

力學分析軟件或程序主要基于連續方法、非連續方法和連續-非連續方法。后者兼具前兩者的優點(例如,連續方法的計算精度高和非連續方法的適用性廣)正在快速發展。在國際上,有限元與離散元耦合方法已取得了諸多的進展[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)在靜水壓力條件下,當單元數目較多時,巷道圍巖的裂紋以更加彎曲、柔美而非突兀、生硬狀向圍巖深部發展,這與塑性力學中圓環形的各處均破壞的塑性區有所不同。在采動條件下,工作面煤壁超前支承壓力和水平應力均具有峰值,后者與巖層之間的黏結作用有關,這與基于連續介質假定的厚壁筒徑向應力特征存在不同。