趙亞鵬, 孔 亮
(1.青島理工大學(xué)土木工程學(xué)院, 青島 266033; 2.青島理工大學(xué)理學(xué)院, 青島 266033)
不同的工程材料往往具有不同的應(yīng)力應(yīng)變關(guān)系,或者說對應(yīng)不同的強度準則[1-3],一般可將工程材料分為金屬材料和非金屬材料,金屬材料的力學(xué)行為較為簡單,如廣泛采用的Tresca準則,認為金屬材料表現(xiàn)為靜水壓力的非相關(guān)性[4],而在非金屬材料中,典型代表則為巖土材料,主要表現(xiàn)為壓硬性[5],或者稱為靜水壓力相關(guān)性。巖土材料的應(yīng)力應(yīng)變關(guān)系較為復(fù)雜,現(xiàn)代土力學(xué)理論一般采用增量的彈塑性理論來解決巖土材料進入塑性階段后的力學(xué)問題[6-8],而這種非彈性的彈塑性力學(xué)關(guān)系即為典型的非線性關(guān)系。工程材料的非線性力學(xué)問題如果按照普通的解析方法求解較為困難且不實際,在計算機發(fā)展早期,由于硬件的局限性,數(shù)值解法的應(yīng)用并不突出,隨著計算機的革命,計算機性能越來越高,數(shù)值解法的優(yōu)越性逐漸凸顯并被廣泛應(yīng)用,并隨之出現(xiàn)了一系列的基于不同算法的數(shù)值模擬軟件[9]。根據(jù)不同研究角度,廣大學(xué)者開展了一系列的數(shù)值模擬研究工作,并取得了大量研究成果,但是針對非線性問題的數(shù)值模擬軟件選取及其建模過程的分析相對較少,而模擬軟件的選取及建模又對分析結(jié)果有著重要影響,因此有必要對模擬軟件選取的相關(guān)內(nèi)容作進一步分析。首先對數(shù)值模擬進行概括性分類總結(jié)及介紹,并對其發(fā)展趨勢進行簡要分析,在此基礎(chǔ)上結(jié)合具體工程(采礦領(lǐng)域)實例就非線性問題的數(shù)值軟件選取相關(guān)內(nèi)容展開討論,可為相似地質(zhì)條件下的非線性問題數(shù)值軟件選取及建模提供參考。
依據(jù)不同的劃分原則,可將數(shù)值模擬劃分為不同類別,例如:根據(jù)求解微分方程方法的不同,可以劃分為有限元法、有限差分法、邊界元法、有限體積法等[10];根據(jù)研究對象的連續(xù)與否,則可劃分為連續(xù)介質(zhì)分析方法與非連續(xù)介質(zhì)方法[11]。
有限單元法(finite element method ,F(xiàn)EM)以變分原理和加權(quán)余量法為基礎(chǔ),其基本原理是將連續(xù)的研究對象離散成有限個單元,各單元之間由有限個節(jié)點相聯(lián)系[12],即將一系列的微分方程離散化,再編制相應(yīng)的計算機程序,借助計算機進行求解。而根據(jù)更為具體的分析方法,有限元法又可分為里茲法、伽遼金法等[13]。隨著計算方法的不斷發(fā)展,一種有別于傳統(tǒng)有限單元法,而基于廣義有限單元法的多尺度單元法被提出[14],由于該方法的基函數(shù)具有多尺度特性,因而被廣泛應(yīng)用于解決多孔介質(zhì)滲流問題[15]。目前基于有限單元法的數(shù)值軟件眾多,且多用于解決固體力學(xué)及結(jié)構(gòu)力學(xué)問題,常用的模擬軟件有ANSYS、ABAQUS、COMSOL等。有限單元法從本質(zhì)上來說屬于連續(xù)介質(zhì)分析方法。
離散單元法(discrete element method ,DEM)在求解方式上屬于顯示求解方法,其與有限單元法類似,都是將連續(xù)的研究對象劃分成有限個單元,但是各單元之間的節(jié)點可以分離,這與有限法元截然不同。單元之間的作用力則由力與位移的關(guān)系求出,個別單元的運動規(guī)律可由牛頓第二定律得出[16]。由于離散單元法所劃分單元之間的離散特性,該方法可以對二維及三維空間內(nèi)的連續(xù)對象進行離散顆粒之間的力學(xué)行為研究,而當顆粒的尺寸足夠小時,則可反映研究對象的微觀屬性,因此離散元法又被認為是宏微觀研究相結(jié)合的有效橋梁[17]。離散單元法在非連續(xù)介質(zhì)領(lǐng)域如散體力學(xué)、斷裂力學(xué)、爆炸力學(xué)方面得到了廣泛應(yīng)用[18-19],此外也被廣泛應(yīng)用于結(jié)構(gòu)工程、巖土工程等領(lǐng)域[20]。離散單元法數(shù)值軟件的典型代表為顆粒流程序(particle flow code ,PFC),此外還包括UDEC(universal distinct element code)、3DEC(3 dimension distinct element code)等。離散單元法屬于非連續(xù)介質(zhì)分析方法。
有限體積法(finite volume method ,F(xiàn)VM)首先將研究區(qū)域劃分為有限個體積,然后再進行方程式的構(gòu)建,區(qū)別于有限單元法和離散單元法,有限體積法基于積分形式的守恒方程而不是微分方程。該方法從物理觀點來構(gòu)造離散方程[21],每一個離散方程都是有限大小體積上某種物理量守恒的表示式,物理意義明確。有限體積法被廣泛應(yīng)用于流體力學(xué)研究領(lǐng)域[22-23]以及用以解決傳熱傳質(zhì)學(xué)等問題[24]。有限體積法的常用數(shù)值軟件有FLUENT、FLOTHERM、STAR-CD等。
邊界單元法(boundary element method ,BEM)是在有限元法之后衍生并發(fā)展起來的一種數(shù)值方法,其主要特點在于僅在連續(xù)介質(zhì)的邊界上劃分單元,用滿足控制方程的函數(shù)去逐漸逼近邊界條件[25],這與有限元在連續(xù)介質(zhì)內(nèi)劃分單元有著很大區(qū)別。正是由于邊界元法的特殊性,使得對于某些特殊的研究對象,可以達到化整為零、維度降低的效果,在保持較高精度的同時大幅簡化計算過程[26]。邊界單元法在斷裂力學(xué)、運動邊界尤其是解決無限與半無限問題方面具有十分廣泛的應(yīng)用[27-28]。與其他方法相比,目前專門的邊界元軟件相對較少,有LINFLOW等。該方法也屬于連續(xù)介質(zhì)分析方法。
有限差分法(finite difference method,F(xiàn)DM)首先對連續(xù)介質(zhì)的研究對象進行網(wǎng)格劃分,然后基于所劃分的網(wǎng)格點,以差商代替微商[29],本質(zhì)是將研究問題離散為差分格式進而求解,其又可具體分為顯示差分與隱式差分。有限差分法具有廣泛的適用性,易于編程且容易利用計算機實現(xiàn),但是對研究區(qū)域的連續(xù)性要求較高。在傳統(tǒng)有限差分法基礎(chǔ)上,近期又發(fā)展出一種無網(wǎng)格形式的廣義有限差分法[30-31],進一步擴展了有限差分的適用性。有限差分法的典型模擬軟件為FLAC3D,其被廣泛應(yīng)用于巖土工程、地下工程及采礦工程[32-33],尤其適用于解決大變形問題。
從20世紀50年代數(shù)值模擬技術(shù)誕生起,其發(fā)展歷程就與計算機技術(shù)緊密聯(lián)系在一起,計算機技術(shù)與數(shù)學(xué)等多學(xué)科的交叉融合則共同促進了數(shù)值模擬領(lǐng)域的快速發(fā)展,相應(yīng)地,數(shù)值模擬又對各研究領(lǐng)域起到了反哺促進的作用。經(jīng)過數(shù)十年的發(fā)展,數(shù)值模擬方法已經(jīng)成為各個研究領(lǐng)域較為成熟的常規(guī)分析方法,通過對大量數(shù)值模擬研究成果的總結(jié),可以認為當前數(shù)值模擬研究有著如下的三大發(fā)展趨勢。
數(shù)值模擬在誕生初期受限于計算機硬件以及算法的不足,僅能對低維、低相以及單一物理場進行模擬,其研究成果也僅能作為特定前提條件下的結(jié)論參考。隨著計算機發(fā)展以及各種算法的提出,同時基于工程實際的需求,實現(xiàn)三維狀態(tài)下的多相及多物理場耦合影響研究成為廣大學(xué)者們的一致訴求。目前對于多相及多物理場模擬研究已經(jīng)取得了一定成果,如應(yīng)力-滲流耦合研究[34]、溫度-應(yīng)力耦合研究[35]、溫度-應(yīng)力-滲流-化學(xué)耦合研究等[36]。
根據(jù)前述可知,數(shù)值模擬方法可以分為有限元、離散元、邊界元及有限差分等不同類別,而各種方法都有其自身優(yōu)缺點,即不同方法具有不同的適用性。經(jīng)過長期的研究工作,發(fā)現(xiàn)某些研究對象往往具有多重屬性,或者同一研究對象在不同區(qū)域具有特定特性,又或者研究對象本身是由兩種不同屬性結(jié)構(gòu)組合而成(如對含斷層地層進行分析時,斷層及完整巖層則具有不同屬性),此時如果將計算方法相耦合則可以在集不同方法優(yōu)點的同時彌補單一方法的不足。較為常見的耦合計算方法如有限元與離散元耦合[37]、有限差分與離散元耦合[38]、邊界元與離散元耦合等[39]。
數(shù)值模擬作為一種研究方法,其最終目的仍然是指導(dǎo)工程實踐,這與實驗研究、理論分析等研究方法相一致。而不同研究方法之間以及研究結(jié)論與工程實踐之間往往會有一定出入,正如前述分析,無論是多物理場耦合還是計算方法耦合,都是為使研究結(jié)果更貼近于工程實際。然而大量研究結(jié)果表明,對于某些特殊環(huán)境中的研究對象,采取常規(guī)方式并不能得出理想結(jié)果,此時對數(shù)值模擬所依賴的本構(gòu)模型等進行改造顯然是一種更深層面的解決辦法。例如考慮非共軸影響的二次開發(fā)數(shù)值模擬[40]、基于FLAC3D黏土邊界面模型的隧道施工分析[41]、基于ABAQUS平臺的纖維彈塑性本構(gòu)模型開發(fā)等[42]。
研究區(qū)域包含3上及3下兩層可采煤層,煤厚分別為5 m和4 m,煤層間距4 m,煤層傾角1°~7°,平均4°,屬半亮-半暗型氣煤,色黑、性脆,似玻璃光澤,硬度2.2~2.5。區(qū)域內(nèi)對工作面影響較大斷層主要有田崗斷層和二龍崗斷層,田崗斷層走向169°~186°,傾向79°~96°,傾角∠70°~75°,落差H為300~410 m,屬臺階狀斷層;二龍崗斷層走向180°~188°,傾向270°~278°,傾角∠65°,落差40~83 m,兩斷層均屬不導(dǎo)水正斷層。距離3煤底板26 m為三灰含水層,平均厚度7 m,水壓最大2.9 MPa。工作面布置分別為23上614及23下614工作面。工作面及斷層位置關(guān)系見圖1。
圖1 工作面平面圖Fig.1 Plan of working face
基于地質(zhì)概況,并結(jié)合工程實際,主要研究內(nèi)容如下。
(1)研究在近距離煤層下行式開采過程中的覆巖運動規(guī)律。
(2)針對三灰含水層,研究底板承壓水對回采過程的影響規(guī)律。
(3)對工作面開采過程中,斷層的活化規(guī)律進行研究,尤其是田崗臺階狀斷層的活化規(guī)律。
(4)研究確定田崗及二龍崗斷層的合理煤柱尺寸。
數(shù)值模擬軟件的選取需要考慮諸多因素,基于以上工程概況,主要考慮如下角度。
3.3.1 實際工程地質(zhì)概況
數(shù)值軟件的選取最根本的出發(fā)點,應(yīng)是以實際地質(zhì)概況為基礎(chǔ)。本工程需要考慮的因素主要有:
(1)田崗斷層與二龍崗斷層。即在數(shù)值模擬時,需要考慮斷層的存在,這就需要模擬軟件能夠解決斷層與周圍巖層的接觸問題[43]。
(2)三灰含水層。由于承壓水的存在,在進行普通應(yīng)力場分析時,還需要考慮滲流場的問題,即數(shù)值模擬軟件必須能夠進行流固耦合的多場分析。
3.3.2 研究內(nèi)容及目的
在正式進行一項研究工作之前,首先應(yīng)確定研究目的及研究內(nèi)容,此后的研究工作也應(yīng)以此為依據(jù)而具體展開。不同的研究內(nèi)容及目的,一般所采用的研究手段也不盡相同,因此就試驗與模擬的研究手段來說,無論是試驗設(shè)備還是模擬軟件都應(yīng)根據(jù)具體的研究內(nèi)容及目的進行具體化的選擇。針對本工程,試驗與數(shù)值模擬部分在考慮前述研究內(nèi)容與目的基礎(chǔ)上將具體結(jié)合其他考慮因素進行綜合確定。
3.3.3 與其他研究方法相對應(yīng)
這一因素主要是基于研究方法“整體性”的考慮。本工程除了采用數(shù)值模擬的研究手段外,還采用了“相似材料模擬試驗”的研究方法。圖2所示為相似材料模擬試驗設(shè)計圖及鋪設(shè)完成圖。
圖2 相似材料模擬試驗圖Fig.2 Diagram of similar material simulation experiment
基于研究手段之間的“呼應(yīng)”關(guān)系,數(shù)值模擬部分應(yīng)該將“相似材料模擬試驗”內(nèi)容作為軟件選取及建模時的考慮因素。就此試驗部分而言,數(shù)值模擬部分的模型規(guī)模應(yīng)與試驗對應(yīng),即至少要包括“臺階田崗斷層”、二龍崗斷層、三灰含水層以及上覆和下伏巖層,此外,其所選軟件的模擬結(jié)果應(yīng)盡可能與試驗結(jié)果相對應(yīng)。
3.3.4 研究對象及行業(yè)領(lǐng)域
主要是指整體性的研究領(lǐng)域。就本工程而言,其首先屬于采礦工程領(lǐng)域、其次可將其擴展為巖土工程及地下工程領(lǐng)域。基于此,首先在領(lǐng)域范圍上,對于數(shù)值模擬軟件的選取就已經(jīng)有了一定的傾向性(更為具體的軟件選取考慮角度將在后續(xù)內(nèi)容中提到),如巖土工程領(lǐng)域常用的FLAC3D、ANSYS、PFC、ABAQUS等[44-46]。
3.3.5 其他因素
在此將以上四類之外的考慮因素以及更為具體化的考慮因素統(tǒng)稱為其他因素,主要包括:
(1)本構(gòu)模型。不同的研究對象其應(yīng)力應(yīng)變關(guān)系不同,因而選取的本構(gòu)模型也不同,最為代表性的就是金屬材料和巖土材料。本工程涉及的煤層開挖問題,屬于巖土的本構(gòu)關(guān)系,一般巖土的本構(gòu)關(guān)系,莫爾庫倫準則相對比較準確[47-48]。因此考慮含有莫爾庫倫模型的模擬軟件。
(2)開挖問題。本工程主要研究煤層開挖問題,因此,模擬軟件必須能夠?qū)r體開挖進行模擬,模擬開挖問題的軟件也較多,如FLAC3D、ANSYS、PLAXIS等[49-51]。
(3)大變形問題。大型滑坡問題以及本工程所涉及的煤層開挖問題一般都屬于大變形問題,如本工程涉及兩層煤層8 m巖層的開挖,其在走向長度上工作面的長度可達數(shù)百米,因此可以想象本工程的數(shù)值模擬必然面臨大變形的問題。與有限元主要解決小變形問題不同,在大變形及大面積屈服方面,F(xiàn)LAC3D具有明顯的優(yōu)勢[52]。
(4)模型規(guī)模。模型規(guī)模包含模型尺寸、網(wǎng)格多少、節(jié)點數(shù)量等內(nèi)容。本工程針對的對象較為復(fù)雜,包括多條斷層,此外水平巖層多達20余層,而針對斷層活化,斷層帶需要進行網(wǎng)格的細化,此外還要考慮豎向斷層與水平巖層的對接問題,可想而知,整個模型的節(jié)點及網(wǎng)格數(shù)量龐大,其最終導(dǎo)致計算過程也必然較長。FLAC3D采用的拉格朗日算法在計算中不形成剛度矩陣,不需迭代滿足彈塑性本構(gòu)關(guān)系,而只需通過應(yīng)變來計算應(yīng)力。因此對于大體積、劃分單元及結(jié)點數(shù)較多的非線性巖體的開挖計算來說是較適宜的[53]。
(5)建模問題。前面已經(jīng)分析,本次模擬所需建立的模型規(guī)模較大,前期建模過程也將相對復(fù)雜,在建模方面,ANSYS優(yōu)勢明顯,而FLAC3D在建模方面較為“呆板”,但是針對FLAC3D建模方面的短板,可借助“RHINO”軟件強大的建模能力,其可對包括FLAC3D、ANSYS在內(nèi)的多款模擬軟件進行模型文件輸出[54-55]。
(6)覆巖垮落問題。本工程中的一個研究角度就是對近距離煤層開采時的覆巖垮落規(guī)律進行研究,此內(nèi)容在試驗部分已經(jīng)進行了研究,因此考慮數(shù)值模擬中也可對此內(nèi)容進行部分體現(xiàn)。由于FLAC3D采用顯式算法獲得模型全部運動方程(包括內(nèi)變量)的時間步長解,從而可以追蹤材料的漸近破壞和垮落[56-57],這與本次研究內(nèi)容相對應(yīng)。
基于以上分析,本工程采用FLAC3D作為所選取的模擬軟件,建模的思路(過程)如下。
3.4.1 地質(zhì)巖層匯總
首先對地質(zhì)巖層類別進行匯總與歸納,主要以實際地質(zhì)鉆孔所獲得的鉆孔柱狀圖為依據(jù),獲得模擬范圍內(nèi)的各巖層厚度及巖層分類。
3.4.2 確定模型尺寸
模型高度及寬度尺寸的確定主要以“研究角度范圍”(如本工程研究覆巖垮落,則需要考慮對縱向高度的覆巖垮落范圍進行預(yù)留)、“研究對象范圍”(如本工程涉及的兩斷層、含水層)、“邊界范圍(要考慮模型的邊界效應(yīng))”及“其他因素”(如與試驗研究相對應(yīng))為側(cè)重點。
3.4.3 模型建立
可采用“RHINO”軟件對模型進行建立,主要針對復(fù)雜模型,如本工程的模型或模型表面有起伏等復(fù)雜模型,建立完成后輸出為FLAC3D文件;也可采用ANSYS軟件進行建模,其與FLAC3D可進行對接;或者對于相對簡單的模型采用CAD進行模型的“粗建”,獲得所需各點坐標后,再根據(jù)坐標在FLAC3D中進行建模,此方法雖然較為煩瑣,但是網(wǎng)格劃分效果較好,對于本工程中斷層帶的網(wǎng)格細化較為有利;對于更為簡單的模型,當然應(yīng)直接在FLAC3D中進行模型建立。
本工程綜合考慮,最終采用了首先在CAD中進行“粗建”,獲得坐標后再進行FLAC3D軟件建模的方法,各點坐標在EXCEL表格中采用公式計算,并將FLAC3D的建模命令也轉(zhuǎn)移至EXCEL文件中進行,最終將所有命令導(dǎo)出為“TXT”文件,供FLAC3D直接讀取。其整體的建模過程也并不十分復(fù)雜,可供參考。所建模型尺寸為:長×寬×高=500 m×200 m×220 m,模型節(jié)點、網(wǎng)格數(shù)量分別為317 460和299 750,最終建立完成的模型如圖3所示。
圖3 模型網(wǎng)格劃分及斷層示意圖Fig.3 Model meshing and fault division
模擬軟件的模擬過程通常主要研究應(yīng)力場,其次還涉及滲流場、溫度場等多場的耦合[58-59],而具體來說則主要有應(yīng)力、位移、塑性區(qū)、水壓等各種模型數(shù)據(jù)。除了整體的模型云圖外,常規(guī)監(jiān)測方式還包括測點以及測線的布設(shè),對分析過程中的以上數(shù)據(jù)進行監(jiān)測,而這些測點與測線的設(shè)置應(yīng)當遵循一定的原則。
(1)以研究內(nèi)容、對象及研究目的為設(shè)置依據(jù)。
(2)選定的測點及測線位置應(yīng)當具有一定的區(qū)域代表性或宏觀代表性。即所獲得的數(shù)據(jù)規(guī)律應(yīng)具有一定區(qū)域或類似工程概況的普適性。
(3)避免重復(fù)設(shè)置。應(yīng)當避免對類似測點的重復(fù)性設(shè)置,一方面,重復(fù)設(shè)置獲取的數(shù)據(jù)規(guī)律具有重復(fù)性、類似性;另一方面,會增加計算及后處理的工作量。在具有代表性的基礎(chǔ)上,重復(fù)性設(shè)置顯然也是沒有必要的。
(4)差異性設(shè)置。即當研究角度、側(cè)重點不同時,即使是相同或類似的研究對象,也應(yīng)避免“一概而論”性設(shè)定。
本次模擬在整體云圖分析的基礎(chǔ)上,也布設(shè)了大量的應(yīng)力、位移、水壓測點和測線。通過對本次模型云圖及測點的分析,其呈現(xiàn)極強的規(guī)律性,并與文獻[60-61]結(jié)論一致,這也說明本次模擬過程較為成功。本次模型布設(shè)的相關(guān)測點及測線主要有:兩級臺階斷層帶應(yīng)力、水壓測線、煤層底板不同深度(10、20、50 m)應(yīng)力、位移測線、煤層間位移測線等。
應(yīng)當進一步說明的是:針對斷層帶應(yīng)力測點,在一級臺階斷層及二級臺階斷層的斷層界面處沿煤層頂、底板每隔10 m布設(shè)一應(yīng)力測點,其中頂板測點距離煤層(3上煤)的高度設(shè)定為正值,底板測點距離煤層的深度設(shè)為負值;各坐標方向如圖3所示。提取相關(guān)數(shù)據(jù)后整理的數(shù)據(jù)曲線如圖4~圖6所示。
圖4 3上煤煤柱尺寸為80 m時,兩級斷層界面應(yīng)力變化Fig.4 The stress change law of two-stage fault interface when the pillar size of 3upper coal is 80 m
圖5 3上煤煤柱尺寸為80 m時,底板應(yīng)力位移圖Fig.5 The change law of stress and displacement of floor when the pillar size of 3upper coal is 80 m
圖6 3上、3下煤開采對比Fig.6 Comparison of 3upper coal and 3lower coal
根據(jù)以上繪制的應(yīng)力、位移、水壓變化曲線,可知,數(shù)據(jù)呈現(xiàn)極強的變化規(guī)律,且符合煤層開挖過程中的應(yīng)力、位移變化規(guī)律,這說明無論是軟件選取、模型建立還是開挖過程模擬,都是比較合理的,這也進一步印證了前文所述內(nèi)容的相對合理性。
數(shù)值模擬與試驗研究相結(jié)合的方式已經(jīng)成為各個研究領(lǐng)域常規(guī)的研究構(gòu)架,不同研究方法之間的結(jié)果印證也使得宏觀結(jié)果更具說服力。正是在這一思路下結(jié)合具體工程實例就試驗研究基礎(chǔ)上模擬研究的一般考慮因素或構(gòu)思過程展開討論,文中所述內(nèi)容具有一定的宏觀代表性,同時又具有特定案例專屬性,可為相似工程案例數(shù)值模擬軟件選取及建模提供參考。主要結(jié)論如下:
(1)數(shù)值模擬的發(fā)展趨勢為:多相及多物理場耦合、計算方法耦合以及二次開發(fā)數(shù)值模擬。
(2)數(shù)值模擬軟件選取時考慮因素主要有:實際工程地質(zhì)概況、研究內(nèi)容及目的、與其他研究手段相對應(yīng)、研究對象及行業(yè)領(lǐng)域、其他因素(具體因素)。
(3)建模思路或過程為(結(jié)合本工程):地質(zhì)巖層匯總、模型尺寸確定、模型建立(不同方法建模)。
(4)測點與測線的設(shè)置應(yīng)當遵循的原則為:以研究內(nèi)容、對象及目的為依據(jù)、具有一定的區(qū)域代表性或宏觀代表性、避免重復(fù)設(shè)置、差異性設(shè)置。
(5)數(shù)據(jù)曲線呈現(xiàn)出極強的規(guī)律性,并與相關(guān)文獻結(jié)論一致,印證了文中軟件選取、模型建立、開挖過程以及測點布設(shè)的正確及合理性。