基于工程實例的非線性問題數值軟件選取分析

趙亞鵬， 孔 亮

(1.青島理工大學土木工程學院， 青島 266033； 2.青島理工大學理學院， 青島 266033)

不同的工程材料往往具有不同的應力應變關系，或者說對應不同的強度準則[1-3]，一般可將工程材料分為金屬材料和非金屬材料，金屬材料的力學行為較為簡單，如廣泛采用的Tresca準則，認為金屬材料表現為靜水壓力的非相關性[4]，而在非金屬材料中，典型代表則為巖土材料，主要表現為壓硬性[5]，或者稱為靜水壓力相關性。巖土材料的應力應變關系較為復雜，現代土力學理論一般采用增量的彈塑性理論來解決巖土材料進入塑性階段后的力學問題[6-8]，而這種非彈性的彈塑性力學關系即為典型的非線性關系。工程材料的非線性力學問題如果按照普通的解析方法求解較為困難且不實際，在計算機發展早期，由于硬件的局限性，數值解法的應用并不突出，隨著計算機的革命，計算機性能越來越高，數值解法的優越性逐漸凸顯并被廣泛應用，并隨之出現了一系列的基于不同算法的數值模擬軟件[9]。根據不同研究角度，廣大學者開展了一系列的數值模擬研究工作，并取得了大量研究成果，但是針對非線性問題的數值模擬軟件選取及其建模過程的分析相對較少，而模擬軟件的選取及建模又對分析結果有著重要影響，因此有必要對模擬軟件選取的相關內容作進一步分析。首先對數值模擬進行概括性分類總結及介紹，并對其發展趨勢進行簡要分析，在此基礎上結合具體工程(采礦領域)實例就非線性問題的數值軟件選取相關內容展開討論，可為相似地質條件下的非線性問題數值軟件選取及建模提供參考。

1 數值模擬分類

依據不同的劃分原則，可將數值模擬劃分為不同類別，例如:根據求解微分方程方法的不同，可以劃分為有限元法、有限差分法、邊界元法、有限體積法等[10]；根據研究對象的連續與否，則可劃分為連續介質分析方法與非連續介質方法[11]。

1.1 有限單元法

有限單元法(finite element method ，FEM)以變分原理和加權余量法為基礎，其基本原理是將連續的研究對象離散成有限個單元，各單元之間由有限個節點相聯系[12]，即將一系列的微分方程離散化，再編制相應的計算機程序，借助計算機進行求解。而根據更為具體的分析方法，有限元法又可分為里茲法、伽遼金法等[13]。隨著計算方法的不斷發展，一種有別于傳統有限單元法，而基于廣義有限單元法的多尺度單元法被提出[14]，由于該方法的基函數具有多尺度特性，因而被廣泛應用于解決多孔介質滲流問題[15]。目前基于有限單元法的數值軟件眾多，且多用于解決固體力學及結構力學問題，常用的模擬軟件有ANSYS、ABAQUS、COMSOL等。有限單元法從本質上來說屬于連續介質分析方法。

1.2 離散單元法

離散單元法(discrete element method ，DEM)在求解方式上屬于顯示求解方法，其與有限單元法類似，都是將連續的研究對象劃分成有限個單元，但是各單元之間的節點可以分離，這與有限法元截然不同。單元之間的作用力則由力與位移的關系求出，個別單元的運動規律可由牛頓第二定律得出[16]。由于離散單元法所劃分單元之間的離散特性，該方法可以對二維及三維空間內的連續對象進行離散顆粒之間的力學行為研究，而當顆粒的尺寸足夠小時，則可反映研究對象的微觀屬性，因此離散元法又被認為是宏微觀研究相結合的有效橋梁[17]。離散單元法在非連續介質領域如散體力學、斷裂力學、爆炸力學方面得到了廣泛應用[18-19]，此外也被廣泛應用于結構工程、巖土工程等領域[20]。離散單元法數值軟件的典型代表為顆粒流程序(particle flow code ，PFC)，此外還包括UDEC(universal distinct element code)、3DEC(3 dimension distinct element code)等。離散單元法屬于非連續介質分析方法。

1.3 有限體積法

有限體積法(finite volume method ，FVM)首先將研究區域劃分為有限個體積，然后再進行方程式的構建，區別于有限單元法和離散單元法，有限體積法基于積分形式的守恒方程而不是微分方程。該方法從物理觀點來構造離散方程[21]，每一個離散方程都是有限大小體積上某種物理量守恒的表示式，物理意義明確。有限體積法被廣泛應用于流體力學研究領域[22-23]以及用以解決傳熱傳質學等問題[24]。有限體積法的常用數值軟件有FLUENT、FLOTHERM、STAR-CD等。

1.4 邊界單元法

邊界單元法(boundary element method ，BEM)是在有限元法之后衍生并發展起來的一種數值方法，其主要特點在于僅在連續介質的邊界上劃分單元，用滿足控制方程的函數去逐漸逼近邊界條件[25]，這與有限元在連續介質內劃分單元有著很大區別。正是由于邊界元法的特殊性，使得對于某些特殊的研究對象，可以達到化整為零、維度降低的效果，在保持較高精度的同時大幅簡化計算過程[26]。邊界單元法在斷裂力學、運動邊界尤其是解決無限與半無限問題方面具有十分廣泛的應用[27-28]。與其他方法相比，目前專門的邊界元軟件相對較少，有LINFLOW等。該方法也屬于連續介質分析方法。

1.5 有限差分法

有限差分法(finite difference method，FDM)首先對連續介質的研究對象進行網格劃分，然后基于所劃分的網格點，以差商代替微商[29]，本質是將研究問題離散為差分格式進而求解，其又可具體分為顯示差分與隱式差分。有限差分法具有廣泛的適用性，易于編程且容易利用計算機實現，但是對研究區域的連續性要求較高。在傳統有限差分法基礎上，近期又發展出一種無網格形式的廣義有限差分法[30-31]，進一步擴展了有限差分的適用性。有限差分法的典型模擬軟件為FLAC3D，其被廣泛應用于巖土工程、地下工程及采礦工程[32-33]，尤其適用于解決大變形問題。

2 數值模擬發展趨勢

從20世紀50年代數值模擬技術誕生起，其發展歷程就與計算機技術緊密聯系在一起，計算機技術與數學等多學科的交叉融合則共同促進了數值模擬領域的快速發展，相應地，數值模擬又對各研究領域起到了反哺促進的作用。經過數十年的發展，數值模擬方法已經成為各個研究領域較為成熟的常規分析方法，通過對大量數值模擬研究成果的總結，可以認為當前數值模擬研究有著如下的三大發展趨勢。

2.1 多相及多物理場耦合

數值模擬在誕生初期受限于計算機硬件以及算法的不足，僅能對低維、低相以及單一物理場進行模擬，其研究成果也僅能作為特定前提條件下的結論參考。隨著計算機發展以及各種算法的提出，同時基于工程實際的需求，實現三維狀態下的多相及多物理場耦合影響研究成為廣大學者們的一致訴求。目前對于多相及多物理場模擬研究已經取得了一定成果，如應力-滲流耦合研究[34]、溫度-應力耦合研究[35]、溫度-應力-滲流-化學耦合研究等[36]。

2.2 計算方法耦合

根據前述可知，數值模擬方法可以分為有限元、離散元、邊界元及有限差分等不同類別，而各種方法都有其自身優缺點，即不同方法具有不同的適用性。經過長期的研究工作，發現某些研究對象往往具有多重屬性，或者同一研究對象在不同區域具有特定特性，又或者研究對象本身是由兩種不同屬性結構組合而成(如對含斷層地層進行分析時，斷層及完整巖層則具有不同屬性)，此時如果將計算方法相耦合則可以在集不同方法優點的同時彌補單一方法的不足。較為常見的耦合計算方法如有限元與離散元耦合[37]、有限差分與離散元耦合[38]、邊界元與離散元耦合等[39]。

2.3 二次開發數值模擬

數值模擬作為一種研究方法，其最終目的仍然是指導工程實踐，這與實驗研究、理論分析等研究方法相一致。而不同研究方法之間以及研究結論與工程實踐之間往往會有一定出入，正如前述分析，無論是多物理場耦合還是計算方法耦合，都是為使研究結果更貼近于工程實際。然而大量研究結果表明，對于某些特殊環境中的研究對象，采取常規方式并不能得出理想結果，此時對數值模擬所依賴的本構模型等進行改造顯然是一種更深層面的解決辦法。例如考慮非共軸影響的二次開發數值模擬[40]、基于FLAC3D黏土邊界面模型的隧道施工分析[41]、基于ABAQUS平臺的纖維彈塑性本構模型開發等[42]。

3 實例分析

3.1 工程概況

研究區域包含3上及3下兩層可采煤層，煤厚分別為5 m和4 m，煤層間距4 m，煤層傾角1°～7°，平均4°，屬半亮-半暗型氣煤，色黑、性脆，似玻璃光澤，硬度2.2～2.5。區域內對工作面影響較大斷層主要有田崗斷層和二龍崗斷層，田崗斷層走向169°～186°，傾向79°～96°，傾角∠70°～75°，落差H為300～410 m，屬臺階狀斷層；二龍崗斷層走向180°～188°，傾向270°～278°，傾角∠65°，落差40～83 m，兩斷層均屬不導水正斷層。距離3煤底板26 m為三灰含水層，平均厚度7 m，水壓最大2.9 MPa。工作面布置分別為23上614及23下614工作面。工作面及斷層位置關系見圖1。

圖1 工作面平面圖Fig.1 Plan of working face

3.2 研究內容

基于地質概況，并結合工程實際，主要研究內容如下。

(1)研究在近距離煤層下行式開采過程中的覆巖運動規律。

(2)針對三灰含水層，研究底板承壓水對回采過程的影響規律。

(3)對工作面開采過程中，斷層的活化規律進行研究，尤其是田崗臺階狀斷層的活化規律。

(4)研究確定田崗及二龍崗斷層的合理煤柱尺寸。

3.3 軟件選取

數值模擬軟件的選取需要考慮諸多因素，基于以上工程概況，主要考慮如下角度。

3.3.1 實際工程地質概況

數值軟件的選取最根本的出發點，應是以實際地質概況為基礎。本工程需要考慮的因素主要有：

(1)田崗斷層與二龍崗斷層。即在數值模擬時，需要考慮斷層的存在，這就需要模擬軟件能夠解決斷層與周圍巖層的接觸問題[43]。

(2)三灰含水層。由于承壓水的存在，在進行普通應力場分析時，還需要考慮滲流場的問題，即數值模擬軟件必須能夠進行流固耦合的多場分析。

3.3.2 研究內容及目的

在正式進行一項研究工作之前，首先應確定研究目的及研究內容，此后的研究工作也應以此為依據而具體展開。不同的研究內容及目的，一般所采用的研究手段也不盡相同，因此就試驗與模擬的研究手段來說，無論是試驗設備還是模擬軟件都應根據具體的研究內容及目的進行具體化的選擇。針對本工程，試驗與數值模擬部分在考慮前述研究內容與目的基礎上將具體結合其他考慮因素進行綜合確定。

3.3.3 與其他研究方法相對應

這一因素主要是基于研究方法“整體性”的考慮。本工程除了采用數值模擬的研究手段外，還采用了“相似材料模擬試驗”的研究方法。圖2所示為相似材料模擬試驗設計圖及鋪設完成圖。

圖2 相似材料模擬試驗圖Fig.2 Diagram of similar material simulation experiment

基于研究手段之間的“呼應”關系，數值模擬部分應該將“相似材料模擬試驗”內容作為軟件選取及建模時的考慮因素。就此試驗部分而言，數值模擬部分的模型規模應與試驗對應，即至少要包括“臺階田崗斷層”、二龍崗斷層、三灰含水層以及上覆和下伏巖層，此外，其所選軟件的模擬結果應盡可能與試驗結果相對應。

3.3.4 研究對象及行業領域

主要是指整體性的研究領域。就本工程而言，其首先屬于采礦工程領域、其次可將其擴展為巖土工程及地下工程領域。基于此，首先在領域范圍上，對于數值模擬軟件的選取就已經有了一定的傾向性(更為具體的軟件選取考慮角度將在后續內容中提到)，如巖土工程領域常用的FLAC3D、ANSYS、PFC、ABAQUS等[44-46]。

3.3.5 其他因素

在此將以上四類之外的考慮因素以及更為具體化的考慮因素統稱為其他因素，主要包括：

(1)本構模型。不同的研究對象其應力應變關系不同，因而選取的本構模型也不同，最為代表性的就是金屬材料和巖土材料。本工程涉及的煤層開挖問題，屬于巖土的本構關系，一般巖土的本構關系，莫爾庫倫準則相對比較準確[47-48]。因此考慮含有莫爾庫倫模型的模擬軟件。

(2)開挖問題。本工程主要研究煤層開挖問題，因此，模擬軟件必須能夠對巖體開挖進行模擬，模擬開挖問題的軟件也較多，如FLAC3D、ANSYS、PLAXIS等[49-51]。

(3)大變形問題。大型滑坡問題以及本工程所涉及的煤層開挖問題一般都屬于大變形問題，如本工程涉及兩層煤層8 m巖層的開挖，其在走向長度上工作面的長度可達數百米，因此可以想象本工程的數值模擬必然面臨大變形的問題。與有限元主要解決小變形問題不同，在大變形及大面積屈服方面，FLAC3D具有明顯的優勢[52]。

(4)模型規模。模型規模包含模型尺寸、網格多少、節點數量等內容。本工程針對的對象較為復雜，包括多條斷層，此外水平巖層多達20余層，而針對斷層活化，斷層帶需要進行網格的細化，此外還要考慮豎向斷層與水平巖層的對接問題，可想而知，整個模型的節點及網格數量龐大，其最終導致計算過程也必然較長。FLAC3D采用的拉格朗日算法在計算中不形成剛度矩陣，不需迭代滿足彈塑性本構關系，而只需通過應變來計算應力。因此對于大體積、劃分單元及結點數較多的非線性巖體的開挖計算來說是較適宜的[53]。

(5)建模問題。前面已經分析，本次模擬所需建立的模型規模較大，前期建模過程也將相對復雜，在建模方面，ANSYS優勢明顯，而FLAC3D在建模方面較為“呆板”，但是針對FLAC3D建模方面的短板，可借助“RHINO”軟件強大的建模能力，其可對包括FLAC3D、ANSYS在內的多款模擬軟件進行模型文件輸出[54-55]。

(6)覆巖垮落問題。本工程中的一個研究角度就是對近距離煤層開采時的覆巖垮落規律進行研究，此內容在試驗部分已經進行了研究，因此考慮數值模擬中也可對此內容進行部分體現。由于FLAC3D采用顯式算法獲得模型全部運動方程(包括內變量)的時間步長解，從而可以追蹤材料的漸近破壞和垮落[56-57]，這與本次研究內容相對應。

3.4 建模思路(過程)

基于以上分析，本工程采用FLAC3D作為所選取的模擬軟件，建模的思路(過程)如下。

3.4.1 地質巖層匯總

首先對地質巖層類別進行匯總與歸納，主要以實際地質鉆孔所獲得的鉆孔柱狀圖為依據，獲得模擬范圍內的各巖層厚度及巖層分類。

3.4.2 確定模型尺寸

模型高度及寬度尺寸的確定主要以“研究角度范圍”(如本工程研究覆巖垮落，則需要考慮對縱向高度的覆巖垮落范圍進行預留)、“研究對象范圍”(如本工程涉及的兩斷層、含水層)、“邊界范圍(要考慮模型的邊界效應)”及“其他因素”(如與試驗研究相對應)為側重點。

3.4.3 模型建立

可采用“RHINO”軟件對模型進行建立，主要針對復雜模型，如本工程的模型或模型表面有起伏等復雜模型，建立完成后輸出為FLAC3D文件；也可采用ANSYS軟件進行建模，其與FLAC3D可進行對接；或者對于相對簡單的模型采用CAD進行模型的“粗建”，獲得所需各點坐標后，再根據坐標在FLAC3D中進行建模，此方法雖然較為煩瑣，但是網格劃分效果較好，對于本工程中斷層帶的網格細化較為有利；對于更為簡單的模型，當然應直接在FLAC3D中進行模型建立。

本工程綜合考慮，最終采用了首先在CAD中進行“粗建”，獲得坐標后再進行FLAC3D軟件建模的方法，各點坐標在EXCEL表格中采用公式計算，并將FLAC3D的建模命令也轉移至EXCEL文件中進行，最終將所有命令導出為“TXT”文件，供FLAC3D直接讀取。其整體的建模過程也并不十分復雜，可供參考。所建模型尺寸為：長×寬×高=500 m×200 m×220 m，模型節點、網格數量分別為317 460和299 750，最終建立完成的模型如圖3所示。

圖3 模型網格劃分及斷層示意圖Fig.3 Model meshing and fault division

3.5 測點設計

模擬軟件的模擬過程通常主要研究應力場，其次還涉及滲流場、溫度場等多場的耦合[58-59]，而具體來說則主要有應力、位移、塑性區、水壓等各種模型數據。除了整體的模型云圖外，常規監測方式還包括測點以及測線的布設，對分析過程中的以上數據進行監測，而這些測點與測線的設置應當遵循一定的原則。

(1)以研究內容、對象及研究目的為設置依據。

(2)選定的測點及測線位置應當具有一定的區域代表性或宏觀代表性。即所獲得的數據規律應具有一定區域或類似工程概況的普適性。

(3)避免重復設置。應當避免對類似測點的重復性設置，一方面，重復設置獲取的數據規律具有重復性、類似性；另一方面，會增加計算及后處理的工作量。在具有代表性的基礎上，重復性設置顯然也是沒有必要的。

(4)差異性設置。即當研究角度、側重點不同時，即使是相同或類似的研究對象，也應避免“一概而論”性設定。

本次模擬在整體云圖分析的基礎上，也布設了大量的應力、位移、水壓測點和測線。通過對本次模型云圖及測點的分析，其呈現極強的規律性，并與文獻[60-61]結論一致，這也說明本次模擬過程較為成功。本次模型布設的相關測點及測線主要有：兩級臺階斷層帶應力、水壓測線、煤層底板不同深度(10、20、50 m)應力、位移測線、煤層間位移測線等。

應當進一步說明的是：針對斷層帶應力測點，在一級臺階斷層及二級臺階斷層的斷層界面處沿煤層頂、底板每隔10 m布設一應力測點，其中頂板測點距離煤層(3上煤)的高度設定為正值，底板測點距離煤層的深度設為負值；各坐標方向如圖3所示。提取相關數據后整理的數據曲線如圖4～圖6所示。

圖4 3上煤煤柱尺寸為80 m時，兩級斷層界面應力變化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時，底板應力位移圖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

根據以上繪制的應力、位移、水壓變化曲線，可知，數據呈現極強的變化規律，且符合煤層開挖過程中的應力、位移變化規律，這說明無論是軟件選取、模型建立還是開挖過程模擬，都是比較合理的，這也進一步印證了前文所述內容的相對合理性。

4 結論

數值模擬與試驗研究相結合的方式已經成為各個研究領域常規的研究構架，不同研究方法之間的結果印證也使得宏觀結果更具說服力。正是在這一思路下結合具體工程實例就試驗研究基礎上模擬研究的一般考慮因素或構思過程展開討論，文中所述內容具有一定的宏觀代表性，同時又具有特定案例專屬性，可為相似工程案例數值模擬軟件選取及建模提供參考。主要結論如下：

(1)數值模擬的發展趨勢為：多相及多物理場耦合、計算方法耦合以及二次開發數值模擬。

(2)數值模擬軟件選取時考慮因素主要有：實際工程地質概況、研究內容及目的、與其他研究手段相對應、研究對象及行業領域、其他因素(具體因素)。

(3)建模思路或過程為(結合本工程)：地質巖層匯總、模型尺寸確定、模型建立(不同方法建模)。

(4)測點與測線的設置應當遵循的原則為：以研究內容、對象及目的為依據、具有一定的區域代表性或宏觀代表性、避免重復設置、差異性設置。

(5)數據曲線呈現出極強的規律性，并與相關文獻結論一致，印證了文中軟件選取、模型建立、開挖過程以及測點布設的正確及合理性。