基于QGA-SVM的堆石料離散元細觀參數標定模型

楊 杰,馬春輝,程 琳,冉 蠡,黃志鴻

(1.西安理工大學水利水電學院,陜西 西安 710048;2.西安理工大學省部共建西北旱區生態水利國家重點實驗室,陜西 西安 710048)

隨著世界范圍內水資源開發力度的不斷加大,堆石壩已成為眾多高壩大庫的推薦壩型,目前正朝向300 m級高壩發展。由于受試驗原理、試驗環境、試驗費用等因素的制約,堆石料室內三軸試驗難以全面、準確、真實地反應筑壩材料力學特性,是造成堆石壩運行期沉降超出預期的原因之一[1]。近年來,離散元憑借其可從細觀尺度探究宏觀現象、物理力學關系明確等優勢,在土工試驗數值模擬方面得到廣泛應用,但離散元仍存在細觀參數標定困難等問題,需開展進一步研究。如何快速、準確地標定堆石料離散元細觀參數,對研究堆石料材料力學性能、確保水利水電工程安全具有重要意義。

顆粒接觸模型及細觀參數標定是離散元數值模擬土工試驗的重要環節,對模擬結果有著直接、明顯的影響,也是制約離散元方法發展的瓶頸之一。由于實際巖土體細觀結構的復雜性和當前硏究的局限性,尚沒有形成一套完善的力學理論用于描述土體細觀與宏觀參數之間的定量關系,細觀參數標定多數采用反復嘗試和人為調試的方式,存在較大的盲目性,且效率低、精度差。當前關于離散元細觀與宏觀參數間互相影響的研究多集中于定性分析,在定量分析方面的研究相對較少。Yang等[2]研究了巖石平行黏結模型中,顆粒數量、細觀參數與宏觀參數以及單軸抗壓強度的關系;Yoon[3]采用中心合成設計對接觸黏結模型的細觀參數進行標定,并研究了各細觀參數與宏觀力學特性的相關性;徐小敏等[4]建立了線性接觸模型的顆粒法向剛度、顆粒剛度比等細觀彈性常數與顆粒材料宏觀彈性常數的經驗公式;趙彥國等[5]系統研究了平行黏結模型中細觀參數對宏觀特性的影響,并提出了細觀與宏觀特性的理論公式;周喻[6]采用BP神經網絡建立了巖土體細觀與宏觀力學參數間的關系;周博[7]擬合了黏性材料內摩擦角、黏聚力的多元非線性公式,定量地描述細觀參數和宏觀剪切強度參數的聯合關系,并給出了切向與法向黏結強度比的建議值;Cheng等[8]采用連續準蒙特卡洛法對土顆粒細觀參數進行了標定。在堆石料研究方面,李守巨[9]以堆石料室內三軸試驗的實測應力-應變曲線為目標,采用響應面法反演離散元三軸試驗模型中的接觸剛度、切向剛度、摩擦系數,取得了良好的效果。以上研究為離散元細觀參數標定提供了新思路、新方法,但標定模型多對細觀參數數目進行簡化處理,且模型缺乏處理小樣本、多參數、復雜非線性關聯的能力;另外上述研究多以巖石的摩爾-庫倫強度參數作為宏觀參數開展研究,相關成果難以推廣至雙曲線模型等更為復雜、參數更多的堆石料本構模型中。

近年來,機器學習算法迅速發展,其在處理小樣本、非線性、多輸出方面具有明顯優勢,為上述問題的求解提供了新手段。本文針對堆石料離散元三軸試驗中存在的細觀參數標定困難問題,以室內三軸試驗實測應力-應變曲線為標定目標,采用機器學習理論建立細觀與宏觀參數間關系,以期迅速、準確的完成堆石料細觀參數標定。

1 堆石料離散元接觸模型

在離散元三軸試驗數值模擬過程中,模型的顆粒數目、粒徑級配、組裝方式、加載速度、接觸模型、細觀參數等眾多因素均會對試件的應力-應變曲線造成影響,其中接觸模型及細觀參數是影響最大、最難調控的因素。接觸模型是離散元描述顆粒變形、受力、運動等狀態的基礎,離散元接觸模型可分為剛度模型、滑動模型、黏結模型三類[6-7]。其中,剛度模型建立了顆粒間接觸力和相對位移間的關系;滑動模型建立法向力和切向力之間的關系,判別兩個接觸體是否發生相對運動;黏結模型是在剛度模型的基礎上,明確了法向拉力和切向力的上限值。黏結模型又可分為接觸黏結模型和平行黏結模型,當顆粒黏結破壞后將遵循滑動模型。接觸黏結模型通過黏結點將球與球連接,能夠設定法向與切向黏結力,其原理如圖1所示,圖中Fn為顆粒的法向接觸力,Un為法向位移,kn為法向接觸剛度,Fs為顆粒的切向接觸力,Fsc為顆粒切向黏結破裂時的接觸力,Fsmax為顆粒發生滑動時的接觸力,Us為切向位移,ks為切向接觸剛度。平行黏結模型假定球體接觸后接觸處被其他黏結性材料填充,黏結性材料的有效剛度與球體接觸點的剛度并聯連接,可設定黏結力和黏結力矩。

圖1 接觸黏結模型原理[7]

根據接觸黏結模型基本原理,若顆粒間運動趨勢為法向擠壓,其壓力與位移關系服從線性剛度模型;若顆粒間運動趨勢為法向脫離,其法向拉力與法向位移成正比,且法向拉力達到設定法向黏結力時黏結破裂,顆粒不再有法向拉力;若顆粒間運動趨勢為切向擠壓,其切向力與切向位移成正比,且切向力達到設定切向黏結力時黏結破裂,顆粒將遵循滑動模型。因此,黏結模型是對線性剛度模型的完善與發展。

近年來,國內外關于離散元堆石料三軸試驗研究的接觸模型選擇與細觀參數取值統計情況如表1所示,涉及到的細觀參數包括:顆粒的法向接觸剛度kn,切向接觸剛度ks,摩擦系數μ;顆粒的法向黏結力bn,切向黏結力bs;側墻的法向接觸剛度knw1,上下加壓板的法向接觸剛度knw2,側墻的切向接觸剛度ksw1、上下加壓板的切向接觸剛度ksw2;孔隙率n等。對于采用接觸黏結模型的堆石料離散元三軸試驗,通常是顆粒間仍采用線性剛度模型,對于粒徑較大的塊石由采用接觸黏結模型的顆粒簇代替生成,以模擬塊石的復雜形狀、破碎過程。綜合分析可知:①堆石料離散元三軸試驗研究所采用的試樣,已逐步由簡單的二維模型向更為復雜的三維模型發展;②早期的堆石料細觀研究多采用簡單的線性剛度模型,隨著研究的發展,接觸黏結模型被廣泛用于構建堆石料中的顆粒簇,以深入研究堆石料破碎發展過程等問題;③由于受到眾多因素的影響,細觀模型參數變化幅度較大甚至存在數量級的差別,難以總結可推廣的堆石料細觀參數標定準則;④除個別案例外,接觸模型若采用接觸黏結模型,其接觸剛度參數數值較采用線性剛度模型有大幅度減小;⑤雖然接觸黏結模型更符合堆石料力學特性,但相比于線性剛度模型,其涉及參數數量更多、取值范圍更廣、細觀與宏觀參數間作用機理更為復雜,因此應對接觸黏結模型的細觀參數確定開展進一步的研究。

表1 堆石料離散元三軸試驗接觸模型與細觀參數取值統計

注：二維試樣尺寸為寬度×長度;三維圓柱試樣尺寸為直徑×高度；三維立方體試樣尺寸為寬度×長度×高度。

不同于巖石、土等材料,堆石料由形狀不規則、多棱角、排列緊密的顆粒組成,具有咬合力大、抗剪強度高等特點。若僅用圓球模擬堆石料,勢必造成顆粒排列形式單一,咬合力較弱。結合上述分析結論,為了更好地模擬堆石料顆粒物理力學特性,提高離散元三軸試驗模擬精度,本研究將建立cluster堆石料顆粒簇,顆粒簇將采用接觸黏結模型模擬堆石料力學特性。

2 QGA-SVM標定模型建立

2.1 標定模型目標函數

在離散元數值模擬的基礎上,建立基于量子遺傳算法(quantum genetic algorithm,QGA)和支持向量機(support vector machine, SVM)的堆石料離散元細觀參數標定模型,其目標函數為

S(x1,x2,…,xM)=

(1)

式中:x1,x2,…,xM為M個待確定的離散元細觀參數;q為應力-應變曲線中提取的定點個數;Fi為QGA-SVM模型的第i個應變值對應的應力值;Ti為室內三軸試驗第i個應變值對應的應力實測值。

2.2 QGA-SVM標定模型

QGA是由Narayanan等[21]基于量子基本理論與遺傳算法思維建立的一種算法,是目前量子衍生算法中發展最為成熟的算法之一,具有優秀的全局搜索能力。作為機器學習算法發展最熱門的算法之一,SVM[22]建立在統計學習理論的Vapnik-chervonenkis(VC)維理論和結構風險最小原理的基礎上。SVM可根據有限的樣本信息,尋求模型復雜性和學習能力間的最佳折中,具有較強的理論基礎,其極值解為全局最優解而非局部最小值,對未知樣本有較好的泛化能力。核函數的類型及參數對機器學習算法的性能有至關重要的影響,為進一步提高SVM模型計算能力,本文采用混合核函數[23]:

K(x,xi)=gexp(-‖x-xi‖2/δ2)+

(1-g)(ηxxi+r)d

(2)

式中:K(·)為核函數;x為訓練樣本的輸入;g為組合核函數的待尋優參數;δ為高斯核參數(帶寬參數);η、r、d為多項式核參數。

基于QGA-SVM的堆石料離散元細觀參數標定模型主要包括訓練機器學習模型和搜索細觀參數兩部分,其計算流程如圖2所示。

本文采用拉丁超立方抽樣(latin hypercube sampling,LHS)構建離散元細觀參數組合,使盡可能少次數的離散元計算能夠表達范圍更廣的細觀參數組合。與隨機取樣法、正交設計法等其他抽樣算法相比,LHS具有適用范圍廣、抽樣估值穩定、樣本具有更好的代表性和均勻性等優點[24]。

圖2 基于QGA-SVM的堆石料細觀參數標定模型流程圖

憑借強大的表達、并行計算能力,QGA計算性能對參數依賴性小,其取值較為固定,同時SVM核函數的變化范圍也較為固定。QGA-SVM細觀參數標定模型具有自動適應不同圍壓、不同試樣的優點,盡可能避免傳統智能算法需人為修改模型參數、造成分析結果差別較大和模型推廣能力差的問題。QGA參數取值和SVM混合核函數參數變化范圍如表2所示。

表2 QGA-SVM參數

3 QGA-SVM標定模型驗證

為驗證模型可行性,采用文獻[10,11,25]中的某水電站筑壩石料粒徑級配及其室內三軸試驗結果進行模型驗證,其粒徑級配如圖3所示。根據粒徑級配,模擬試驗采用擠壓排斥法生成離散元三軸試樣,如圖4所示。為模擬堆石料應力應變特性,采用接觸黏結模型模擬堆石料,其中堆石料顆粒間采用線性剛度模型,對于粒徑較大塊石由隨機顆粒組成的cluster顆粒簇進行代替,cluster顆粒簇內部通過bond鍵連接。

圖3 某水電站筑壩石料粒徑級配曲線

圖4 離散元三軸試驗試樣

離散元三軸試件采用圓柱形,尺寸為?300 mm×650 mm,初始孔隙比為0.35。試件共5 026個顆粒,其中cluster顆粒簇包含508個顆粒。由伺服控制程序控制試驗的等壓固結、加載、卸載,圍壓為800 kPa。設定上下加載壓盤的運動速度為0.05 m/s,以模擬試樣的靜力加載。

為實現基于QGA-SVM的堆石料離散元細觀參數標定,采用離散元軟件PFC(particle flow code)進行離散元三軸試驗模擬,使用FISH語言實現細觀參數的讀入、離散元三軸試驗和應力-應變曲線的輸出;隨后采用Matlab對應力-應變曲線進行插值,以獲得0～15%中每0.5%間隔的軸向應變對應應力值,共計31個應力值;最后采用Matlab完成QGA-SVM算法的訓練與計算,實現堆石料離散元細觀參數標定。

綜合分析相關研究,并在一定的試算基礎上確定細觀參數取值范圍,如表3所示。其中側墻法向剛度約為球體顆粒剛度的1/10,墻體的切向剛度設為0。通過LHS在細觀參數取值范圍內構建40組參數組合,利用離散元計算相應的應力-應變曲線,并提取其中31個應力值;將40組材料參數組合作為SVM的輸入數據,相應應力值作為SVM的輸出數據,采用QGA搜索確定SVM核參數,使其性能達到最佳狀態;最后,以室內三軸試驗值與計算應力值誤差最小為目標,發揮SVM快速計算應力-應變曲線的能力,采用QGA全局搜索細觀參數。

表3 堆石料離散元三軸試驗細觀參數取值范圍及標定結果

在SVM訓練中,QGA計算得SVM性能最佳時的混合核參數組合為:δ=1.377 4、η=0.943 3、r=0.163 6、d=2.227 8、g=0.444 5。SVM與離散元計算值的平均絕對誤差(mean absolute error,MAE)為0.015 3,表明SVM對訓練數據的擬合精度非常高,訓練已完全能夠代替離散元實現應力應變計算。隨后采用QGA搜索細觀參數,最終細觀參數標定如表3所示。

通過離散元正算細觀參數標定結果,得室內三軸試驗、SVM與離散元計算值對比如圖5所示。由圖5可知:①SVM計算值與室內三軸試驗的MAE為0.25,兩者在應力-應變曲線的前部與中部吻合度較高,在中前部及后部略有誤差,模型整體表現良好,較好地反映了堆石料變形情況;②SVM與離散元正算值的MAE為0.28,兩者在應力-應變曲線的前部、后部存在一定誤差,表明SVM在前部、后部的訓練需加強;③離散元正算值與室內三軸試驗的MAE為0.19,兩者在應力-應變曲線的前部存在一定誤差,在中部、后部擬合效果較好。計算結果中,離散元正算值相較于SVM計算值更為接近室內三軸試驗值,表明SVM在擬合應力-應變曲線中仍存在一定誤差,致使正算值存在優于SVM計算值的可能。上述分析表明:基于QGA-SVM的堆石料離散元細觀參數標定模型是可行的。

圖5 各類試驗的應變應力對比

為驗證QGA優化速度,本文分別采用QGA與遺傳算法(GA)標定細觀參數,其收斂速度如圖6所示。相比于GA,QGA迭代收斂速度快,計算精度高。因此,本文所建立的QGA-SVM標定模型在計算精度、速度方面優勢明顯。

圖6 優化模型收斂速度對比

4 結 語

通過總結國內外關于堆石料細觀模型選擇與參數選取的研究進展,認為接觸黏結模型更適合作為堆石料的離散元細觀模型。針對堆石料離散元三軸試驗數值模擬中存在的細觀參數標定影響因素多、耗時長、成本高等問題,建立基于QGA-SVM的堆石料細觀參數標定模型。模型采用QGA算法優化SVM核函數完成機器學習的訓練,使其擬合精度、預測精度達到代替離散元計算的要求。隨后,模型以SVM計算應力-應變曲線與室內三軸試驗實測值差值最小為目標,采用QGA搜索離散元細觀參數,完成堆石料離散元三軸試驗細觀參數標定。通過實例證明,QGA-SVM可快速、精確地標定離散元材料細觀模型參數,可將模型進一步擴展至考慮體應變或本構模型參數等方面,具有良好的學術研究、工程應用推廣價值。