基于機(jī)器學(xué)習(xí)和遺傳算法的非局部晶體塑性模型參數(shù)識(shí)別1)

周 瑞 熊宇凱 儲(chǔ)節(jié)磊 闞前華 康國(guó)政 張 旭 ,2)

* (西南交通大學(xué)力學(xué)與航空航天學(xué)院,應(yīng)用力學(xué)與結(jié)構(gòu)安全四川省重點(diǎn)實(shí)驗(yàn)室,成都 610031)

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

引言

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

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

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

本文基于非局部晶體塑性本構(gòu)模型,借助DAMASK平臺(tái)[23]生成鎳基高溫合金的含孔洞模型單軸拉伸模擬數(shù)據(jù),進(jìn)而建立從參數(shù)到應(yīng)力響應(yīng)的機(jī)器學(xué)習(xí)模型.在此基礎(chǔ)上,將機(jī)器學(xué)習(xí)模型與遺傳算法進(jìn)行耦合,通過(guò)已有實(shí)驗(yàn)數(shù)據(jù),確定鎳基高溫合金的非局部晶體塑性模型參數(shù).最后通過(guò)SHAP 工具,分析參數(shù)變化的影響.

1 模型與方法

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

圖1 參數(shù)確定優(yōu)化過(guò)程流程圖Fig.1 Flowchart of the optimization procedure for constitutive parameters

1.1 非局部晶體塑性本構(gòu)模型

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

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

式中,Lp是位錯(cuò)運(yùn)動(dòng)引起的塑性速度梯度,可以表示為

位錯(cuò)運(yùn)動(dòng)時(shí)會(huì)遇到障礙,其平均速度vα可表示為

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

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

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

其中,G是剪切模量,為滑移系之間位錯(cuò)的相互作用強(qiáng)度系數(shù).

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

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

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

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

其中,V是體元體積,an為對(duì)應(yīng)的面法向矢量,An為編號(hào)為n的面元面積.

1.2 遺傳算法模型

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

(1) 基因是評(píng)價(jià)公式中的可調(diào)參數(shù),即所選取的非局部晶體塑性模型參數(shù).

(2) 染色體表示個(gè)體的所有基因集合,即多個(gè)本構(gòu)參數(shù)的集合.染色體代表單獨(dú)個(gè)體,即一組本構(gòu)參數(shù).

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

(4) 適應(yīng)度是對(duì)應(yīng)于評(píng)價(jià)公式得到的數(shù)值,適應(yīng)度越大表明染色體對(duì)應(yīng)的評(píng)價(jià)結(jié)果越優(yōu).

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

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

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

本文將應(yīng)力-應(yīng)變曲線中的屈服應(yīng)力和最終應(yīng)力作為評(píng)價(jià)指標(biāo),通過(guò)確定二者的值即可確定應(yīng)力-應(yīng)變曲線.本文最大模擬應(yīng)變?yōu)?.1,最終應(yīng)力表示的是應(yīng)變?yōu)?.1 時(shí)的應(yīng)力.二者的值可以通過(guò)機(jī)器學(xué)習(xí)模型得到.本文依據(jù)常用的L2范數(shù)[27],建立了對(duì)應(yīng)的評(píng)價(jià)公式評(píng)估適應(yīng)度

其中,Dobj為評(píng)價(jià)值,和 σN分別對(duì)應(yīng)機(jī)器學(xué)習(xí)計(jì)算和實(shí)驗(yàn)得到的應(yīng)變?yōu)?.1 處的應(yīng)力值,和 σ0分別對(duì)應(yīng)機(jī)器學(xué)習(xí)計(jì)算得到的屈服應(yīng)力和實(shí)驗(yàn)得到的屈服應(yīng)力.評(píng)價(jià)值描述了屈服應(yīng)力和最終應(yīng)力與實(shí)驗(yàn)值的最大差值,當(dāng)評(píng)價(jià)值最小時(shí)說(shuō)明計(jì)算結(jié)果與實(shí)驗(yàn)值最為接近.個(gè)體的適應(yīng)度為Fobj=Dm-Dobj,其中Dm表示每代中Dobj的最大值,此時(shí)適應(yīng)度越大,對(duì)應(yīng)的評(píng)價(jià)值越小.

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

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

1.3 可解釋機(jī)器學(xué)習(xí)模型SHAP

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

其中,x′是從原始數(shù)據(jù)集x簡(jiǎn)化得到的輸入特征,M是數(shù)據(jù)集中的特征數(shù)量,φ0是一個(gè)常數(shù),φi是每個(gè)特征i對(duì)應(yīng)的權(quán)重.

2 參數(shù)優(yōu)化過(guò)程

2.1 非局部晶體塑性模型參數(shù)選擇

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

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

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

2.2 數(shù)據(jù)生成

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

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

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

其中,n為樣本數(shù)量,yi為樣本標(biāo)簽,為機(jī)器學(xué)習(xí)的預(yù)測(cè)值.同時(shí)使用決定系數(shù)R2來(lái)描述機(jī)器學(xué)習(xí)的回歸效果,其表示為

圖4 最終應(yīng)力的(a)訓(xùn)練集和(b)測(cè)試集的數(shù)據(jù)分布;屈服應(yīng)力的(c)訓(xùn)練集和(d)測(cè)試集的數(shù)據(jù)分布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 結(jié)果與討論

3.1 遺傳算法計(jì)算結(jié)果

基于以上數(shù)據(jù),按照上述模型設(shè)置分別對(duì)3 個(gè)機(jī)器學(xué)習(xí)模型進(jìn)行訓(xùn)練.并依據(jù)MSE對(duì)模型的訓(xùn)練結(jié)果進(jìn)行評(píng)估,經(jīng)過(guò)機(jī)器學(xué)習(xí)模型參數(shù)調(diào)整后,完成模型的訓(xùn)練,其模型最終性能可以通過(guò)MSE和決定系數(shù)R2表示,如表3 所示.與其他模型相比,LSTM的MSE最小,表明其模型性能最好.這表示了LSTM能夠較好地預(yù)測(cè)應(yīng)力-應(yīng)變響應(yīng)過(guò)程,且具有良好的泛化能力.這是因?yàn)長(zhǎng)STM 能夠有效地學(xué)習(xí)數(shù)據(jù)的時(shí)序特征,該特點(diǎn)與輸入數(shù)據(jù)更為契合.因此,在3 種模型的復(fù)雜程度相差不大的情況下,其MSE最小.而LSTM 對(duì)應(yīng)的R2也最小,但是考慮到R2更適合用于評(píng)價(jià)線性回歸結(jié)果,而LSTM 計(jì)算的標(biāo)簽是非線性的應(yīng)力-應(yīng)變曲線,因此R2難以反映出LSTM 的性能.SVR 和FFNN 使用的標(biāo)簽符合線性回歸,其R2也反映了二者能夠?qū)η?yīng)力和最終應(yīng)力進(jìn)行準(zhǔn)確的預(yù)測(cè).

表3 不同機(jī)器學(xué)習(xí)模型的性能表現(xiàn)Table 3 Performance of different machine learning model

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

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

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

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

表5 根據(jù)不同機(jī)器學(xué)習(xí)模型得到的晶體塑性模型參數(shù)Table 5 Crystal plasticity parameter obtained by different machine learning model

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

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

圖7 基于優(yōu)化程序得到參數(shù)的晶體塑性有限元模擬結(jié)果與實(shí)驗(yàn)結(jié)果的對(duì)比Fig.7 The crystal plasticity results of finite element simulation with the parameters obtained by the optimization program are compared with the experimental results

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

將優(yōu)化求解得到的參數(shù)代入有限元模型中進(jìn)行計(jì)算,結(jié)果如圖8 所示.本文使用了應(yīng)力-應(yīng)變曲線下面積差值的絕對(duì)值作為指標(biāo)來(lái)評(píng)價(jià)曲線與實(shí)驗(yàn)數(shù)據(jù)的擬合情況,各個(gè)模型結(jié)果在表6 中列出.機(jī)器學(xué)習(xí)模型的性能越好,其應(yīng)力-應(yīng)變曲線擬合越好,LSTM,FFNN 和SVR 3 個(gè)模型得到的曲線均能夠很好地與實(shí)驗(yàn)數(shù)據(jù)吻合.在引入硬化指數(shù)作為評(píng)價(jià)指標(biāo)后,擬合效果得到了明顯的提升.這表明,本文引入的硬化指數(shù)能夠有效地反映出實(shí)驗(yàn)數(shù)據(jù)的特征,從而提高了優(yōu)化模型的準(zhǔn)確性.

表6 不同模型得到的應(yīng)力-應(yīng)變曲線與實(shí)驗(yàn)數(shù)據(jù)的擬合情況Table 6 Fitting results of stress-strain curves obtained by different models and experimental data

圖8 引入硬化指數(shù)n 后,有限元模擬結(jié)果與實(shí)驗(yàn)結(jié)果的對(duì)比Fig.8 Comparison of the finite element simulation results and experimental results after introducing the hardening index n

3.2 基于SHAP 的機(jī)器學(xué)習(xí)模型分析

本文借助SHAP 工具,基于Shapley 值的原理對(duì)機(jī)器學(xué)習(xí)模型進(jìn)行解釋.通過(guò)機(jī)器學(xué)習(xí)模型計(jì)算得到了晶體塑性模型參數(shù)在屈服應(yīng)力和最終應(yīng)力下的Shapley 值,并根據(jù)大小進(jìn)行了排序,如圖9 所示.Shapley 值的大小表示了參數(shù)對(duì)結(jié)果的邊際效應(yīng)影響,發(fā)現(xiàn)p和q對(duì)屈服應(yīng)力和最終應(yīng)力均有較大的影響.值得注意的是,k2對(duì)屈服應(yīng)力影響較小,而對(duì)最終應(yīng)力影響較大.

圖9 對(duì)結(jié)果影響最大的特征Fig.9 The most important features affect

為了更加明確地表示各個(gè)特征取值變化對(duì)結(jié)果的影響和作用,我們通過(guò)SHAP 工具,基于數(shù)據(jù)集的數(shù)據(jù)分布,繪制了特征取值對(duì)屈服應(yīng)力和最終應(yīng)力的影響,如圖10 所示.圖中一個(gè)點(diǎn)代表一個(gè)樣本,顏色越紅說(shuō)明特征本身數(shù)值越大,顏色越藍(lán)說(shuō)明特征本身數(shù)值越小.這里根據(jù)數(shù)據(jù)集來(lái)說(shuō)明參數(shù)數(shù)值變化時(shí)對(duì)屈服應(yīng)力和最終應(yīng)力造成的影響.在屈服應(yīng)力圖中,p越大,q和f越小,屈服應(yīng)力越大;而在最終應(yīng)力圖中,p越大,q和k2越小,最終應(yīng)力越大.

圖10 特征的取值產(chǎn)生的影響Fig.10 Effects of features’ values

我們將相應(yīng)的參數(shù)在范圍內(nèi)變化,并帶入有限元模型中進(jìn)行驗(yàn)證,結(jié)果如圖11 所示.由圖可知,當(dāng)p值增大時(shí),屈服應(yīng)力和最終應(yīng)力均增加;當(dāng)q增大時(shí),二者均減小;當(dāng)k2增大時(shí),屈服應(yīng)力僅發(fā)生微小的變化,最終應(yīng)力出現(xiàn)下降的趨勢(shì).有限元計(jì)算結(jié)果與預(yù)測(cè)結(jié)果一致,證明了Shapley 值能夠有效地說(shuō)明特征變化產(chǎn)生的影響.

圖11 基于有限元模擬,參數(shù)(a) p,(b) q,(c) k2 增大時(shí)最終應(yīng)力和屈服應(yīng)力的變化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

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

4 結(jié)論

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

(1) 本文建立耦合機(jī)器學(xué)習(xí)的遺傳算法本構(gòu)模型參數(shù)求解模型,該方法可以根據(jù)應(yīng)力-應(yīng)變曲線確定對(duì)應(yīng)的非局部晶體塑性模型參數(shù),并在精度上滿足要求.

(2) 通過(guò)計(jì)算Shapley 值,可以得到各個(gè)參數(shù)對(duì)屈服應(yīng)力和最終應(yīng)力的邊際影響相對(duì)大小.根據(jù)Shapley 值的變化,可以從數(shù)據(jù)分布中得出特征變化對(duì)結(jié)果產(chǎn)生的影響:p增大或q減小時(shí),屈服應(yīng)力和最終應(yīng)力均增加;而k2增大時(shí),屈服應(yīng)力僅出現(xiàn)微小變化,最終應(yīng)力出現(xiàn)下降趨勢(shì).

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