基于人工神經網絡的顆粒材料本構關系及邊值問題研究

張廣江, 楊德澤, 楚錫華

(武漢大學 工程力學系, 武漢 430072)

0 引 言

顆粒材料通常是指大量相互接觸的離散顆粒與周圍流體介質一起組成的復雜體系,廣泛存在于日常生產生活的許多領域．以沙石、泥土等典型的顆粒材料為例,其被大量地運用于巖土工程、道路工程、地質工程等工程實踐中,還與山體滑坡、土地沙漠化、泥石流等許多地質災害密切相關．為了保證工程的安全性、預測并防止相關地質災害的發生,需要提高對顆粒材料的認識,理解并預測其在不同外荷載作用下的力學響應．

目前,研究顆粒材料物理力學性質的主要數值方法有基于離散途徑的離散單元法和基于連續體途徑的有限單元法．離散單元法[1-2]基于Newton運動定律和力-位移運動定律,通過搜索并計算顆粒之間的相互作用,確定所有顆粒的運動狀態．它能夠有效地反映顆粒的排列方式、形狀、粒度等細觀結構信息對顆粒材料力學性質的影響,在計算過程中更加真實地體現顆粒材料的本質,但是由于顆粒數量眾多,需要消耗大量的計算資源和時間．有限單元法則是從宏觀上將顆粒材料視作一種連續體,采用連續介質理論對其進行分析研究,與離散單元法相比計算效率較高,但分析的結果十分依賴于顆粒材料本構模型的準確性．由于顆粒數量眾多、形狀各異[3]且內部的細觀結構復雜,顆粒材料呈現出非連續性、非均勻性、各向異性[4-5],使其具有了一些十分獨特且復雜的物理力學性質,如“糧倉效應”[6-7]、“分離效應”[8-9]、“壓力凹陷現象”[10]等,導致建立一個可以精確預測顆粒材料力學行為的本構模型十分困難．現有的顆粒材料本構模型往往是基于唯象假設建立起來的,許多模型十分復雜,本構參數的確定存在一定的困難,或者是僅基于材料的某一現象建立,適用范圍存在一定的局限性[11]．

為了解決兩者在模擬顆粒材料時的困難,Munjiza和Owen等[12-13]將離散單元法和有限單元法結合,創建了有限元-離散元耦合方法(FDEM),通過有限單元法模擬顆粒材料的變形,內部顆粒的運動則用離散單元法描述,充分發揮了兩者的優勢．

近年來,隨著計算機科學技術的快速發展,算力迅猛提升,新的算法不斷涌現,機器學習再次引起了人們的關注．機器學習通過“學習算法”,從現有數據中產生模型,并能通過產生的模型對新的情況進行判斷和預測,被認為是一個可能解決材料本構問題的途徑[14]．人工神經網絡作為機器學習中極具代表性的一類算法,具備強大的學習能力,推動了各相關領域的變革．為此,許多研究人員致力于用其解決顆粒材料數值模擬所面臨的一些難題:Benvenuti等[15]利用人工神經網絡將宏觀實驗結果與微觀數值參數聯系起來,提出了一種識別離散元模擬參數的方法; Zhang 等[16]使用長短期記憶神經網絡提出了一種模擬顆粒材料循環特性的方法;Wang等[17]和Qu等[18]分別利用卷積神經網絡和循環神經網絡成功預測了顆粒材料的本構規律．本文嘗試應用人工神經網絡算法,將基于離散顆粒模型的離散單元法與基于連續介質模型的有限單元法有機結合,訓練一個可以在主空間上準確描述顆粒材料本構關系的模型并運用于求解實際邊值問題,形成一套新的、完整的求解方案．

本文基于離散單元法獲取顆粒材料的主應力、主應變以及對應的應力-應變矩陣等數據,并利用人工神經網絡算法建立起可以在主空間上準確描述顆粒材料本構關系的人工神經網絡模型;在假設材料的主應力空間與主應變空間重合的前提下,結合轉軸公式,通過UMAT子程序將模型導入ABAQUS中,并用于求解顆粒材料的邊值問題．數值算例首先考察了平板雙軸壓縮問題,結果表明了該方案的可行性和可靠性,并通過進一步考慮邊坡問題,結果顯示訓練后的人工神經網絡模型能夠在一定程度上預測復雜工況的發展．

1 人工神經網絡算法

人工神經網絡是由具有適應性的簡單單元組成的廣泛并行互連的網絡,能夠模擬生物神經系統對真實世界物體所作出的交互反應[19]．它通常由一個輸入層、多個隱藏層和一個輸出層組成,且每一層都包含多個神經元,其中輸入層通過接收輸入值傳遞給隱藏層,而輸出層則產生預測的結果．

從數學上看,人工神經網絡本質上是一個復合函數,通過層與層之間的線性變換以及神經元上的非線性激活函數,使其具備了強大的擬合以及預測能力,被廣泛地運用于分類以及回歸等任務．對于顆粒材料而言,通過唯象假設建立精確的顆粒材料本構模型十分復雜且困難,在已知有關顆粒材料的應力、應變等數據時, 人工神經網絡能夠直接從數據中學習, 建立應力-應變的內在聯系, 十分有效、 便捷地形成顆粒材料的本構模型．

為建立起一個可以準確描述輸入數據與輸出數據內在關系的人工神經網絡模型,需要對人工神經網絡不斷地進行訓練,即通過向前傳播和向后傳播不斷地調整網絡中神經元之間的連接權以及閾值,直到損失函數在容許誤差范圍內．向前傳播是指神經網絡在訓練過程中,數據從輸入層傳播到輸出層的過程．向后傳播則是指將神經網絡訓練后的預測值與訓練樣本真實值進行比較,將兩者之間的誤差利用神經網絡反向傳播調整連接權以及閾值,使得預測值和真實值之間的誤差不斷下降．

圖1 人工神經網絡訓練過程Fig. 1 The training process of the artificial neural network

2 數據準備與人工神經網絡模型的訓練

2.1 基于離散元法的數值雙軸試驗

離散元法是以Newton運動定律為基礎、接觸力計算為核心,用來求解非連續介質材料問題的一種數值模擬方法,能夠有效地對顆粒材料的物理力學行為進行模擬．與基于唯象模型,通過有限元法對顆粒材料的宏觀力學行為進行模擬相比,離散元法能夠直接把顆粒的排列方式、形狀、粒度、孔隙度等細觀結構參數以及顆粒的接觸本構考慮在內,更加接近顆粒材料的本質屬性．為了獲取材料在一系列不同外荷載作用下的響應,對于標準試樣不論是真實試驗還是數值模擬試驗,相較于獲取其主方向下的應力和應變,獲取其一般方向下的應力和應變通常更加困難．因此,本文基于PFC2D,采用離散元法模擬雙軸試驗,產生顆粒材料的主應力、主應變以及對應的應力-應變矩陣等相關數據,用于訓練人工神經網絡模型．

離散元數值雙軸試驗的標準試樣如圖2所示,其中不同的顏色代表不同的顆粒半徑,樣本生成模型參數見表1．在一個區域生成約7 200個、半徑為2～4 mm隨機分布的顆粒后,通過伺服方法,移動試樣的四個邊界來施加雙向等壓壓力,直到試樣穩定達到100 kPa的圍壓,伺服后通過改變上下、左右墻面的移動速率,來控制雙向產生不同的應變比．根據顆粒材料的特性,本文主要考慮以壓為主的情況,所有工況的最大軸向應變為7.5%,一共產生了202組數據作為訓練集或驗證集、40組數據作為測試集,訓練集、驗證集、測試集的應變加載路徑如圖3所示．

表1 樣本模型參數

圖2 離散元數值雙軸試樣 圖3 加載路徑示意圖 Fig. 2 The discrete element numerical biaxial sample Fig. 3 Schematic diagram of the loading path

為了使訓練后的人工神經網絡模型能夠替代顆粒材料的本構關系用于求解邊值問題,離散元數值試驗輸出的數據主要是主應力、主應變以及對應的應力-應變矩陣．其中,主應力為顆粒樣本施加在墻面上的平均反作用力,主應變為墻體的相對位移,應力-應變矩陣[20]為

Cijkl=(kn-ks)aijkl+ksbijkl,

(1)

其中,kn,ks分別為顆粒之間的法向接觸剛度以及切向接觸剛度;aijkl,bijkl分別為表征元的組構張量．

平面情況下,主應力與主應變的增量關系用矩陣形式表示如下:

(2)

2.2 人工神經網絡模型訓練過程

為了訓練出完整可靠的人工神經網絡模型,便于通過ABAQUA子程序UMAT求解實際邊值問題,采用了兩條神經網絡[20]分別訓練應變與應力、應力-應變矩陣之間的關系．經過多次試驗,人工神經網絡結構形式如表2所示．兩條神經網絡均基于全連接的人工神經網絡,采用標準BP算法[21],其中激活函數為sigmoid函數,損失函數為均方誤差函數．

表2 神經網絡結構參數

為了更快、更精準地訓練出人工神經網絡模型,所有樣本的輸入數據、輸出數據在進行訓練前根據下面公式進行歸一化處理,使得所有樣本的輸入數據、輸出數據均位于[-1,1]之間．

(3)

兩條神經網絡的訓練誤差如圖4所示,為了避免過擬合,在指定訓練次數的前提下,取驗證集誤差最小時的參數為神經網絡模型的參數．從圖5(a)、(b)可以看出,訓練后的人工神經網絡模型均能在99%以上解釋應變與應力、應力-應變矩陣之間的關系,表明了模型在一定應變范圍內具有良好的預測能力．值得注意的是,在后續實際求解邊值問題的過程中,訓練后的人工神經網絡1模型無論是對求解的精度還是求解過程的收斂性都至關重要,而人工神經網絡2模型并不影響求解的精度,只對求解過程的收斂性產生重大影響．

(a) 神經網絡1 (b) 神經網絡2 (a) Neural network 1 (b) Neural network 2

(a) 神經網絡1 (b) 神經網絡2 (a) Neural network 1 (b) Neural network 2

3 基于人工神經網絡模型的UMAT子程序開發

為了將訓練后的人工神經網絡模型運用于求解邊值問題,本文通過UMAT子程序將訓練后的人工神經網絡模型導入ABAQUS,并用于求解邊值問題．以平面問題為例,UMAT子程序流程圖如圖6所示．

圖6 UMAT程序流程Fig. 6 The UMAT program flow chart

為了使得訓練后的人工神經網絡模型在UMAT中參與計算,假定主應力空間與主應變空間重合,主軸坐標系與一般直角坐標系之間需進行轉換．

平面問題轉軸公式為:Pij=PmnLimLjn,令σ1應力主方向與x軸之間的夾角為θ．由于主方向上剪應力為0,即

γ12=-2sinθcosθ×εx+2sinθcosθ×εy+(cos2θ-sin2θ)γxy=0,

(4)

得出

(5)

已知θ后,分別通過坐標轉換1將應變轉換為主應變,坐標轉換2將主應力轉換為應力:

(6)

在對主方向下的D矩陣進行維度轉換后,通過坐標轉換3將D矩陣從主應力與主應變的關系轉換為一般直角坐標系下應力與應變的關系,即

(7)

(8)

其中

4 數 值 算 例

4.1 雙軸壓縮試驗

為了通過有限元法求解顆粒材料的邊值問題,傳統方法往往需要將離散元模型中的細觀材料參數轉換成宏觀材料參數,建立以經典彈性模型和Mohr-Coulomb屈服準則為基礎形成的彈塑性模型,其函數表達式分別為

σij=Cijklεkl,

(9)

(10)

(11)

其中p為等效應力,q為Mises應力,c為黏聚力,φ為內摩擦角,θ為Lode角．

宏觀參數標定方法如下:

1) 在圍壓達到100 kPa后,保持左右伺服應力不變,在上下兩側緩慢地施加位移荷載,在材料變形為線性變化時,通過下面公式分別計算彈性模量和Poisson比:

(12)

2) 通過改變圍壓的大小,分別為100 kPa,200 kPa,300 kPa,之后保持左右伺服應力不變,在上下兩側緩慢地施加位移荷載,得到材料的應力-應變曲線(圖7(a)),再根據Mohr-Coulomb理論繪制Mohr-Coulomb包絡線(圖7(b)),并計算相關的材料參數．

(a) 應力-應變曲線 (b) Mohr-Coulomb包絡線(a) Stress-strain curves (b) The Mohr-Coulomb envelope

材料的宏觀力學參數最終標定為:彈性模量E=6.5×107Pa,Poisson比μ=0.3,內摩擦角32°,剪脹角8°,黏聚力250 Pa．

為了驗證離散元法產生的數據訓練后的人工神經網絡模型是否可以成功替代顆粒材料的本構模型,通過有限元法用于求解實際邊值問題,分別運用經典彈塑性模型與人工神經網絡模型模擬一個5 m×5 m的平板雙軸壓縮試驗,對比兩者的數值計算結果．

離散元模擬過程中通過伺服機制產生了大小為100 kPa的初始應力,使得訓練后的人工神經網絡模型隱含了該初始應力條件．為了使運用經典彈塑性模型求解與人工神經網絡模型求解時的實際加載路徑一致,在運用經典彈塑性模型進行求解時,先施加初始應力(圖8(a)),再通過位移加載(圖8(b))的方式進行雙軸壓縮,而運用人工神經網絡模型求解時,由于訓練后的人工神經網絡模型隱含了該初始條件,則直接施加相同的位移荷載(圖8(b))．

圖8 平板壓縮試驗Fig. 8 The plate compression test

人工神經網絡模型與經典彈塑性模型的計算結果如圖9所示,分別為位移云圖、應變云圖和應力云圖．在施加相同位移荷載的情況下,應變與應力的計算結果基本一致,其中應變最大值與應力最大值的相對誤差分別為9.78%,10.44%．從應變云圖可以看出,兩者的變形模式相一致,均會出現明顯的應變局部化現象,應變主要集中在一個“X”形的剪切帶中．為了進一步分析人工神經網絡模型與經典彈塑性模型在雙軸壓縮試驗中的計算結果,我們統計了所有單元共計1 600個Gauss積分點處的Mises應力值并計算相對誤差,如圖10所示,其中柱狀圖表示位于指定誤差范圍內的Gauss積分點數目,曲線圖表示在指定誤差下Gauss積分點數目所占的比例．兩個模型計算雙軸壓縮試驗的平均相對誤差為11.8%,驗證了利用人工神經網絡算法建立起可以在主空間上有效描述顆粒材料本構關系的人工神經網絡模型,并運用于求解邊值問題的可行性．

(a) 經典彈塑性模型(a) The classical elastoplastic model

圖10 平板應力誤差曲線Fig. 10 The stress error curve of plate

由于經典彈塑性模型和人工神經網絡模型的構建方式不同,經典彈塑性模型是基于唯象假設建立的具備一定適用范圍的本構模型,在運用于具體材料時,需進行參數標定．而人工神經網絡模型則是直接從具體材料的應力、應變等數據中學習,建立了一個適用于該特定材料的本構模型,導致兩者在描述材料的本構關系時存在著一定的區別．

經典彈塑性模型存在顯式的數學表達式,應力-應變曲線是光滑的,而人工神經網絡模型通過數據訓練獲得,往往不存在顯式的數學表達式,其準確性十分依賴于訓練數據的數量與質量,而由于離散元數值模擬過程中存在波動導致應力-應變曲線并不是十分光滑,使得訓練后的神經網絡也不是十分光滑[22]．此外,經典彈塑性模型中彈性模量并不具備應力相關性,即彈性模量始終為一定值,并不會隨著材料應力狀態而發生改變,而理論上人工神經網絡模型能直接從數據中學習,不僅能夠反映顆粒細觀結構對其宏觀力學性能的影響,也能夠在一定程度上捕捉材料彈性模量的變化．正是由于經典彈塑性模型和人工神經網絡模型本質上存在著這些差異,導致了無論是應力還是應變,兩者的計算結果都不可避免地出現一定程度的數值大小差異．

4.2 邊坡試驗

為了研究訓練后的人工神經網絡模型在面對更加復雜的情況時,是否具備一定的能力去模擬并預測破壞的發生,進一步基于新建立的求解方案模擬邊坡問題．邊坡問題是巖土工程中一個十分重要且復雜的問題,不同的邊坡在不同的加載條件下會出現不同的破壞形式,其中滑動破壞是最常見的一種破壞類型,為了防止破壞的發生往往需要對邊坡的穩定性進行分析．本文模擬邊坡問題的數值模型如圖11所示,其中邊坡頂部為一剛性板,施加偏心位移荷載大小為0.06 m,底部邊界為固定約束,右側邊界約束x向位移．圖12(a)、(b)分別為經典彈塑性模型和人工神經網絡模型的計算結果．

圖11 邊坡問題示意圖Fig. 11 Schematic diagram of the slope

通過對比兩者的計算結果,觀察到人工神經網絡模型與經典彈塑性模型模擬的總體趨勢相一致．從應變云圖中可以看出,它們均呈現出明顯的滑動破壞趨勢,且在該位移加載的情況下破壞均未擴展至整個滑坡．通過進一步地提取邊坡試驗中所有Gauss積分點處的Mises應力并計算相對誤差,如圖13(a)所示,結果表明,兩者計算的相對誤差主要集中在50%附近．與平板雙軸壓縮試驗相比,邊坡試驗由于具備更加復雜的加載情況,其計算誤差也相應地明顯提高．可見該神經網絡模型在面對復雜加載情況時還缺乏足夠的能力去精確地捕獲材料的應力-應變響應,只能進行定性的分析．

圖13 邊坡應力誤差圖Fig. 13 Stress errors curve of slope

與經典彈塑性模型相比,人工神經網絡模型的計算誤差主要來源于訓練誤差以及在面對復雜加載情況時,由于主應力空間與主應變空間不重合所造成的誤差．雖然訓練后的人工神經網絡模型在指定的應力-應變范圍內能夠很好地解釋應變-應力、應力-應變矩陣之間的關系,但由于模型泛化能力的限制,圖10中仍存在極少數數據點處的相對誤差遠大于平均誤差的情況．而針對主應力空間與主應變空間可能出現不重合所造成的影響,通過提取邊坡中應變和Mises應力最大的單元,觀察兩個模型計算的正應力相對誤差隨主應力空間與主應變空間夾角變化的趨勢,如圖13(b)所示,兩者總體上呈正相關,正應力的相對誤差隨著主應力空間與主應變空間夾角的增大而增大,表明了主應力空間與主應變空間的不重合是復雜加載情況下誤差的主要來源之一．

5 結 論

本文通過運用人工神經網絡算法,將基于離散顆粒模型的離散單元法與基于連續介質模型的有限單元法有機結合,構建了一個可以在主空間上預測顆粒材料本構關系的人工神經網絡模型,并運用于求解邊值問題,形成了一套完整的求解方案．通過平板雙軸壓縮試驗,驗證了該求解方案的可行性,表明了該人工神經網絡模型在面對簡單加載情況時,能夠有效地捕捉應力-應變等數據之間的關系,但是在面對復雜加載情況時,由于主應力空間與主應變空間不重合使得求解結果出現較大的數值差異,導致該方案只能在一定程度上進行定性分析．