基于機器學習和遺傳算法的非局部晶體塑性模型參數識別1)

周 瑞 熊宇凱 儲節磊 闞前華 康國政 張 旭 ,2)

* (西南交通大學力學與航空航天學院,應用力學與結構安全四川省重點實驗室,成都 610031)

? (西南交通大學計算機與人工智能學院,成都 610031)

引言

晶體塑性有限元方法基于材料在離散滑移系的塑性滑移來研究單晶和多晶材料的塑性變形行為[1].由于在微米尺度下,材料非均勻變形相關的應變梯度效應不可忽略,需要建立可描述幾何必需位錯(geometrically necessary dislocations,GND)對力學行為的影響的非局部晶體塑性本構模型[2-5].非局部晶體塑性模型考慮的位錯機制較多,描述位錯行為的本構模型參數較多.合適的材料參數是晶體塑性模型擁有良好預測能力的先決條件,如何準確確定非局部晶體塑性本構模型的參數極具挑戰.

傳統的本構模型參數的確定方法采用的是“試錯法”,即通過調整材料參數,直至模擬結果與實驗結果或參考數據相匹配[6-7].在參數較少且參數影響較為明顯的情況下,通過試錯法是易于得到晶體塑性參數的.但是對于具有大量參數,且參數間可能會存在非線性作用的復雜本構模型,通過“試錯法”調控參數極為困難.基于梯度的優化方法常常用于求解逆問題,已有學者借助優化算法確定本構模型參數[8-9].Mahnken 等[10]基于梯度的優化方法識別黏塑性材料本構模型的參數,但整個優化過程的收斂性和收斂速度極大地依賴于初始點的選取,仍無法避免不斷的試錯.相比之下,遺傳算法(genetic algorithm,GA)采用了進化程序和交叉重組算法,并不需要有關目標函數的梯度信息.因此,有的學者考慮通過遺傳算法確定本構模型參數[11-13].Savage 等[14]基于遺傳算法,結合了微觀結構數據與應力數據,成功識別了本構模型參數.Frydrych 等[15]基于壓痕實驗數據,通過將遺傳算法的結果帶入有限元方法中進行模擬,識別了晶體塑性模型的參數.然而遺傳算法每次迭代都需要執行一次有限元模擬,提取輸出結果,并通過與實驗數據比較計算出適應度.晶體塑性有限元的模擬計算成本過高,會導致整個優化流程效率低下.因此,遺傳算法是一種更為穩健的優化算法[16],可成功用于參數的識別,但是這種收斂的穩定性是以犧牲效率為代價的[17].

這里本文引入了機器學習方法來避免實際模擬計算成本過高的問題.機器學習模型的計算速度,遠快于有限元的模擬過程,并且將機器學習模型用于替代模型的技術已經比較成熟,許多學者將其用于替代本構模型[18-19],用于建立參數到應力響應的映射.此外,已有學者將機器學習算法與本構模型計算過程相結合,通過機器學習算法加速有限元的計算.Heider 等[20]在回退-映射算法中加入了機器學習模型,加速了各向異性與非彈性材料響應的本構模型的計算過程.Wen 等[19]使用梯度提升樹建立模型,當輸入等效應力、等效應變、溫度和等效應變率時,可以直接預測出塑性應變率和應變硬化模量.他們將這個模型加入有限元計算過程中,實現了計算的加速.其次,已有學者將機器學習方法用于識別晶體塑性模型的參數[21].例如,Veasna 等[22]結合Pareto最優點和機器學習方法,建立了多目標的優化模型,該模型能夠有效確定晶體塑性本構參數.然而,結合機器學習和遺傳算法進行非局部晶體塑性本構模型參數確定的研究還未見報道.

本文基于非局部晶體塑性本構模型,借助DAMASK平臺[23]生成鎳基高溫合金的含孔洞模型單軸拉伸模擬數據,進而建立從參數到應力響應的機器學習模型.在此基礎上,將機器學習模型與遺傳算法進行耦合,通過已有實驗數據,確定鎳基高溫合金的非局部晶體塑性模型參數.最后通過SHAP 工具,分析參數變化的影響.

1 模型與方法

對于簡單的本構模型,可以通過機器學習從應力-應變曲線得到參數.但復雜的本構模型不滿足一一對應關系,一條應力-應變曲線可能對應多組參數,所以難以直接確定本構參數.因此,本文采用了一種新穎的耦合機器學習模型的遺傳算法進行參數優化,能夠有效地規避上述問題.本文采用的算法流程分為3 個階段:第1 階段調控非局部晶體塑性模型的參數,通過晶體塑性有限元模擬生成大量數據;第2 階段是根據已生成的數據,建立3 種不同的機器學習模型;第3 階段是將機器學習模型與遺傳算法相結合,通過不斷的優化迭代獲取最優參數.整個優化過程的流程圖如圖1 所示.

圖1 參數確定優化過程流程圖Fig.1 Flowchart of the optimization procedure for constitutive parameters

1.1 非局部晶體塑性本構模型

本文所采用的本構模型為基于位錯機制的非局部晶體塑性模型,該模型考慮了金屬材料變形過程中多種位錯的力學行為.在有限變形框架下,變形梯度F可以通過乘法分解為彈性變形梯度Fe和塑性變形梯度Fp

式中,Fe描述晶格彈性變形,而Fp則描述由位錯滑移引起的塑性變形.塑性變形梯度的演化率可以表示為

式中,Lp是位錯運動引起的塑性速度梯度,可以表示為

位錯運動時會遇到障礙,其平均速度vα可表示為

其中,cat為原子濃度,vT為黏滯速度,其可以表示為,η為黏滯系數,tS和tP分別表示位錯克服Peierls 障礙和固溶原子障礙需要的時間.位錯在障礙前的等待時間與位錯嘗試跨越障礙的頻率和其成功跨越障礙的概率相關[25],可以表達為

其中,f為位錯嘗試跨越障礙的頻率,冪指數項為位錯成功跨越障礙的概率,kB為玻爾茲曼常數,τS和τP分別為固溶原子障礙和Peierls 障礙的強度,dobst為固溶原子直徑,wk為雙扭折寬度,p和q為描述障礙能壘的系數,τeff為滑移面上的有效分切應力,其可以表示為 α 滑移系上驅動力 τα與位錯運動需要克服的林位錯交互作用力 τcr之差

位錯之間的相互作用力 τcr可以表示為

其中,G是剪切模量,為滑移系之間位錯的相互作用強度系數.

位錯的增殖、湮滅以及位錯極性(單個位錯與位錯偶)的轉變對材料的力學行為有重要的影響.假定不同類型位錯的增殖率相同,則位錯增殖率可以描述為

位錯偶的湮滅類型有兩種:其一為攀移造成的刃型位錯偶湮滅,刃型位錯偶的演化率可以表示為

對于非局部模型,關鍵在于考慮了位錯流動,即位錯信息在材料點之間的交換.滑移系 α 上的位錯流量可以定義為可動位錯密度和其運動速度的乘積,即fα=ραvα.由材料點的位錯信息交換,可以給出位錯密度的演化[25]

位錯的流動通過穿透因子χ控制.使用有限體積離散的方法,位錯流入鄰近體積單元可以表示為

其中,V是體元體積,an為對應的面法向矢量,An為編號為n的面元面積.

1.2 遺傳算法模型

遺傳算法作為一種直接搜索方法,其優化結果通常會得到全局的最優值而不是局部最優值[26].遺傳算法中的部分概念定義如下:

(1) 基因是評價公式中的可調參數,即所選取的非局部晶體塑性模型參數.

(2) 染色體表示個體的所有基因集合,即多個本構參數的集合.染色體代表單獨個體,即一組本構參數.

(3) 種群表示染色體的集合,即多組非局部晶體塑性模型參數的集合.

(4) 適應度是對應于評價公式得到的數值,適應度越大表明染色體對應的評價結果越優.

(5) 選擇階段是根據適應度得到染色體被保留的概率,然后淘汰適應度低的個體.

(6) 交叉遺傳表示種群中兩個染色體間有概率發生基因互換.

(7) 變異表示染色體的基因可能發生隨機變化.

本文將應力-應變曲線中的屈服應力和最終應力作為評價指標,通過確定二者的值即可確定應力-應變曲線.本文最大模擬應變為0.1,最終應力表示的是應變為0.1 時的應力.二者的值可以通過機器學習模型得到.本文依據常用的L2范數[27],建立了對應的評價公式評估適應度

其中,Dobj為評價值,和 σN分別對應機器學習計算和實驗得到的應變為0.1 處的應力值,和 σ0分別對應機器學習計算得到的屈服應力和實驗得到的屈服應力.評價值描述了屈服應力和最終應力與實驗值的最大差值,當評價值最小時說明計算結果與實驗值最為接近.個體的適應度為Fobj=Dm-Dobj,其中Dm表示每代中Dobj的最大值,此時適應度越大,對應的評價值越小.

本文采用的遺傳算法程序流程主要包含兩個階段,第1 個階段是確定非局部晶體塑性模型參數并且隨機生成種群染色體,即隨機化初始參數.第2 階段主要是種群演化過程,首先通過機器學習模型計算種群中每個染色體對應的評價值,然后根據評價值得到適應度.此時種群中的最優值若是滿足精度要求即可停止演化,若不滿足則進入選擇階段.選擇階段將根據適應度淘汰部分個體,并將剩余個體依次進行交叉遺傳和變異,最終得到新一代種群.然后計算新種群的適應度,反復迭代直至滿足要求,如圖2所示.

圖2 非局部晶體塑性本構模型參數確定的遺傳算法流程示意圖Fig.2 Flowchart of genetic algorithm for determing the parameters of nonlocal crystal plasticity model

1.3 可解釋機器學習模型SHAP

本文采用的可解釋的機器學習模型來自Lundberg等[28]提出的SHAP 模型.SHAP 源于博弈論,是一種解釋機器學習模型預測的方法.SHAP 是一種加性特征屬性方法,即其模型的輸出是輸入變量的線性疊加[29].每個特征的貢獻通過Shapley 值[30]表示,因此可解釋模型g(x′) 可以表示為

其中,x′是從原始數據集x簡化得到的輸入特征,M是數據集中的特征數量,φ0是一個常數,φi是每個特征i對應的權重.

2 參數優化過程

2.1 非局部晶體塑性模型參數選擇

本文以鎳基高溫單晶合金作為研究材料.非局部晶體塑性模型中參數眾多,大部分參數可以通過實驗曲線或者材料本征的物理特性來確定,如彈性模型、泊松比以及位錯Burgers 矢量大小等.因此,僅僅與塑性變形相關的本構參數,特別是與位錯密度演化有關的模型參數被選為待擬合參數.待擬合參數會在一定范圍內變化,它們的大小變化和組合在很大程度上影響了位錯演化和塑性屈服與流動行為.參數的范圍參考了文獻中給出的鎳基高溫合金參數取值[31-33].其次在選擇上下限時需滿足相應理論要求,如p和q的取值需分別滿足 0

表1 待擬合的晶體塑性模型參數Table 1 Unfit crystal plasticity parameter

表2 鎳基高溫合金的非局部晶體塑性模型參數Table 2 Nonlocal crystal plasticity parameter for Ni-based superalloys

2.2 數據生成

本文的優化過程以鎳基高溫合金實驗應力-應變曲線的屈服應力和最終應力作為目標結果,其數據源于參考文獻[37].為體現非局部晶體塑性本構的特點,選取了含孔洞模型進行單軸拉伸模擬.在訓練機器學習模型前,需獲得足夠的數據.因此,本文依據實驗中的試樣尺寸,建立了大小為0.5 mm×3 mm×3 mm 的晶體塑性模型,其包含有直徑為0.5 mm 的穿透孔.模型以晶粒取向[001]作為拉伸方向,單元數量為2880,單元類型為C3D8R,如圖3 所示.模型約束了x-y平面x和y方向的位移,并在x-y面施加沿z軸的應變控制載荷,應變率為1.0×10-4s-1,拉伸位移為0.3,因此最終應變為10%.根據建立的晶體塑性模型,在已確定的參數范圍內,隨機調整參數,并提交計算.最后進行后處理,每組數據取200 個應力-應變點組成應力-應變曲線,并提取相應的屈服應力和最終應力.將以上過程重復多次,完成數據生成,共生成1200 組數據.

圖3 鎳基合金有限元模型示意圖(單位:mm)Fig.3 Finite element model of Ni-based alloy (unit:mm)

本文建立了3 種機器學習模型,分別是支持向量回歸模型(support vector regression,SVR)、長短期記憶神經網絡(long short-term memory,LSTM)和前饋神經網絡(feed forward neural network,FFNN).在SVR 模型和FFNN 模型中,將選定的6 個非局部晶體塑性模型參數作為輸入參數,屈服應力和最終應力作為輸出參數.而在LSTM模型中,輸入參數為6 個非局部晶體塑性模型參數和應變,每個參數序列長度均設置為200,輸出結果為應力.在SVR 模型中采用多項式核函數,其多項式次數為3.在神經網絡中,使用均方誤差(mean squared error,MSE)作為損失函數,其可以表示為

其中,n為樣本數量,yi為樣本標簽,為機器學習的預測值.同時使用決定系數R2來描述機器學習的回歸效果,其表示為

圖4 最終應力的(a)訓練集和(b)測試集的數據分布;屈服應力的(c)訓練集和(d)測試集的數據分布Fig.4 The data distribution of final stress in (a) the train dataset and(b) the test dataset;The data distribution of yield stress in (c) the train dataset and (d) the test dataset

3 結果與討論

3.1 遺傳算法計算結果

基于以上數據,按照上述模型設置分別對3 個機器學習模型進行訓練.并依據MSE對模型的訓練結果進行評估,經過機器學習模型參數調整后,完成模型的訓練,其模型最終性能可以通過MSE和決定系數R2表示,如表3 所示.與其他模型相比,LSTM的MSE最小,表明其模型性能最好.這表示了LSTM能夠較好地預測應力-應變響應過程,且具有良好的泛化能力.這是因為LSTM 能夠有效地學習數據的時序特征,該特點與輸入數據更為契合.因此,在3 種模型的復雜程度相差不大的情況下,其MSE最小.而LSTM 對應的R2也最小,但是考慮到R2更適合用于評價線性回歸結果,而LSTM 計算的標簽是非線性的應力-應變曲線,因此R2難以反映出LSTM 的性能.SVR 和FFNN 使用的標簽符合線性回歸,其R2也反映了二者能夠對屈服應力和最終應力進行準確的預測.

表3 不同機器學習模型的性能表現Table 3 Performance of different machine learning model

基于訓練完成的機器學習模型,通過前文敘述的遺傳算法模型進行本構模型參數確定.將遺傳算法的迭代次數設置為500,種群大小為1000,交叉率和變異率分別為65%和10%.以評價公式的結果作為誤差,用以求解非局部晶體塑性模型參數.在算法運算過程中,每代種群中最優個體的誤差變化過程如圖5 所示,最終的誤差值見表4.評價值越小,表明其結果與實驗結果越接近.從圖中可得,使用了不同機器學習模型的遺傳算法模型均可收斂,并且最終的誤差均較小.遺傳算法模型在迭代100 次左右時已接近收斂,而使用FFNN 的遺傳算法模型收斂速度最快,且其最終精度最高.

表4 遺傳算法的最終誤差Table 4 The final error of genetic algorithm

圖5 遺傳算法中每代誤差變化Fig.5 The variation of the error in every generation of genetic algorithm

遺傳算法得到的晶體塑性模型參數見表5.不同機器學習模型的遺傳算法優化結果有一定差異,主要有兩個原因:一是不同模型的精度影響了收斂區域的大小;二是同一條應力-應變曲線可能對應多組晶體塑性模型參數,導致了多解現象,且每個解的收斂區域也不相同.其次,部分參數在使用不同機器學習模型的情況下取值范圍比較接近,例如參數k2在設定最大值為60 的情況下,各個模型下取值都超過50;參數p和參數q的取值都接近1.為了驗證這3 個參數的分布情況,本文采用FFNN 作為機器學習模型,用遺傳算法進行了50 次求解,其參數取值分布如圖6 所示.這3 個參數的數值分布比較集中,其中參數p和參數q的平均值都接近1.參數分布的穩定性說明它們對屈服應力和最終應力有較大的影響,因此模型求解得到的參數接近.

表5 根據不同機器學習模型得到的晶體塑性模型參數Table 5 Crystal plasticity parameter obtained by different machine learning model

圖6 遺傳算法得到的參數最優解分布圖Fig.6 The distribution of the optimal parameters obtained by the genetic algorithm model

將遺傳算法模型得到的參數代入已建立的有限元模型中進行計算,并將計算結果與實驗結果進行對比.如圖7 所示,3 種模型得到的應力-應變曲線整體上均與實驗曲線接近,并且二者的曲線趨勢一致,只有硬化階段存在較小差異.由于評價公式主要關注屈服點和應變最大處的應力,而由遺傳算法計算得到的最終誤差值均不超過0.02,如表4 所示,因此在這兩點處模擬結果與實驗曲線重合程度高.該模型得到的參數能夠滿足要求,參數對應的應力-應變曲線符合實驗趨勢.

圖7 基于優化程序得到參數的晶體塑性有限元模擬結果與實驗結果的對比Fig.7 The crystal plasticity results of finite element simulation with the parameters obtained by the optimization program are compared with the experimental results

注意到模擬結果與實驗結果在塑性階段存在一定偏差.考慮到大多數金屬的應變硬化行為,可以用指數曲線來描述:.其中 σT和 εT分別表示真應力和真應變,K表示強度系數,n表示應變硬化指數.硬化指數n可以描述金屬的應變硬化程度,因此可以將其作為評價指標之一,引入到式(14) 中.LSTM 可以輸出應力-曲線從而計算曲線的硬化指數,故選擇使用LSTM 作為機器學習模型.在遺傳算法中引入硬化指數n作為評價指標,并使用LSTM作為機器學習模型,此新模型記作LSTM-1.

將優化求解得到的參數代入有限元模型中進行計算,結果如圖8 所示.本文使用了應力-應變曲線下面積差值的絕對值作為指標來評價曲線與實驗數據的擬合情況,各個模型結果在表6 中列出.機器學習模型的性能越好,其應力-應變曲線擬合越好,LSTM,FFNN 和SVR 3 個模型得到的曲線均能夠很好地與實驗數據吻合.在引入硬化指數作為評價指標后,擬合效果得到了明顯的提升.這表明,本文引入的硬化指數能夠有效地反映出實驗數據的特征,從而提高了優化模型的準確性.

表6 不同模型得到的應力-應變曲線與實驗數據的擬合情況Table 6 Fitting results of stress-strain curves obtained by different models and experimental data

圖8 引入硬化指數n 后,有限元模擬結果與實驗結果的對比Fig.8 Comparison of the finite element simulation results and experimental results after introducing the hardening index n

3.2 基于SHAP 的機器學習模型分析

本文借助SHAP 工具,基于Shapley 值的原理對機器學習模型進行解釋.通過機器學習模型計算得到了晶體塑性模型參數在屈服應力和最終應力下的Shapley 值,并根據大小進行了排序,如圖9 所示.Shapley 值的大小表示了參數對結果的邊際效應影響,發現p和q對屈服應力和最終應力均有較大的影響.值得注意的是,k2對屈服應力影響較小,而對最終應力影響較大.

圖9 對結果影響最大的特征Fig.9 The most important features affect

為了更加明確地表示各個特征取值變化對結果的影響和作用,我們通過SHAP 工具,基于數據集的數據分布,繪制了特征取值對屈服應力和最終應力的影響,如圖10 所示.圖中一個點代表一個樣本,顏色越紅說明特征本身數值越大,顏色越藍說明特征本身數值越小.這里根據數據集來說明參數數值變化時對屈服應力和最終應力造成的影響.在屈服應力圖中,p越大,q和f越小,屈服應力越大;而在最終應力圖中,p越大,q和k2越小,最終應力越大.

圖10 特征的取值產生的影響Fig.10 Effects of features’ values

我們將相應的參數在范圍內變化,并帶入有限元模型中進行驗證,結果如圖11 所示.由圖可知,當p值增大時,屈服應力和最終應力均增加;當q增大時,二者均減小;當k2增大時,屈服應力僅發生微小的變化,最終應力出現下降的趨勢.有限元計算結果與預測結果一致,證明了Shapley 值能夠有效地說明特征變化產生的影響.

圖11 基于有限元模擬,參數(a) p,(b) q,(c) k2 增大時最終應力和屈服應力的變化Fig.11 Based on finite element simulation,the changes of final stress and yield stress when parameters (a) p,(b) q,and (c) k2 increase

關于該現象可以從非局部晶體塑性本構模型中得到解釋.p和q主要是通過式(6)影響位錯速度vα,進而影響 α 滑移系上的塑性滑移率,因此當p增大或q減小時,對應的 α 滑移系上的塑性滑移率減小,最終應力值也會增大.由于選取的屈服點是根據應力-應變曲線進行選取的,在該點狀態下實際上已經出現了部分滑移系的開動,因此p和q對屈服應力的影響趨勢與最終應力一致.參數k2用于調控位錯運動的平均自由程的大小,所以當k2變大時,位錯運動的平均自由程變大,意味著位錯運動遇到的障礙強度變弱,因此最終應力會減小.因此,根據SHAP 分析得到參數變化產生的影響結果與理論一致,這表明其能夠有效地展示參數變化導致的影響.

4 結論

本文將鎳基高溫合金材料的晶體塑性有限元數據用于訓練機器學習模型,并在遺傳算法中耦合了機器學習模型,確定了鎳基高溫合金的晶體塑性模型參數,得到主要結論如下.

(1) 本文建立耦合機器學習的遺傳算法本構模型參數求解模型,該方法可以根據應力-應變曲線確定對應的非局部晶體塑性模型參數,并在精度上滿足要求.

(2) 通過計算Shapley 值,可以得到各個參數對屈服應力和最終應力的邊際影響相對大小.根據Shapley 值的變化,可以從數據分布中得出特征變化對結果產生的影響:p增大或q減小時,屈服應力和最終應力均增加;而k2增大時,屈服應力僅出現微小變化,最終應力出現下降趨勢.

(3) 基于SHAP 框架得到本構參數變化影響的分析結果可以通過晶體塑性本構理論和有限元模擬進行驗證,并且其變化趨勢與有限元模擬和理論分析一致.