梅 杰，張博文，張春云，彭海峰，崔 苗

(大連理工大學(xué)，工業(yè)裝備結(jié)構(gòu)分析國家重點實驗室，遼寧省空天飛行器前沿技術(shù)重點實驗室，遼寧 大連 116024)

冪強化彈塑性材料在工程領(lǐng)域中的應(yīng)用十分廣泛，這種材料的受力及變形行為分析具有重要的理論價值和實際工程意義。冪強化彈塑性材料的特點是材料的變形規(guī)律會根據(jù)應(yīng)力狀態(tài)而改變，正是這種性質(zhì)導(dǎo)致材料在達到屈服點后呈現(xiàn)非線性的本構(gòu)關(guān)系，使得相關(guān)的力學(xué)問題難以解析求解。在工程問題中，材料的本構(gòu)參數(shù)和結(jié)構(gòu)的邊界條件往往不容易確定，在這種情況下，反問題為確定這些參數(shù)提供了一種新思路。

反問題研究主要分為兩個部分：正問題求解和反問題求解。冪強化材料的平面應(yīng)變正問題的求解方法主要包括解析方法和數(shù)值方法，解析解[1-3]比較難獲得，而數(shù)值解[4]主要是借助商業(yè)軟件或者編寫源程序。編寫源程序需要深厚的理論基礎(chǔ)并且花費大量時間和精力，不太適合工程技術(shù)人員，相比之下，商用軟件的精度與效率都比較高，工程應(yīng)用范圍更廣。

關(guān)于反問題的研究，國內(nèi)外學(xué)者已經(jīng)做了不少的研究工作，提出很多方法。何永勇等[5]基于遺傳算法，實現(xiàn)了轉(zhuǎn)子裂紋的識別。郭紅玲等[6]基于蟻群算法求解彈性本構(gòu)參數(shù)區(qū)間。姜紹飛等[7]對傳統(tǒng)多粒子群協(xié)同優(yōu)化(MPSCO)算法進行分布式并行化改進，開發(fā)了基于云計算的 PMPSCO 算法，提出了基于 PMPSCO 算法的框架結(jié)構(gòu)物理參數(shù)辨識方法。Guo等[8]提出了一種基于貓群優(yōu)化(CSO)算法的單、雙二極管模型參數(shù)估計優(yōu)化技術(shù)，結(jié)果表明該方法性能良好，CSO算法可以有效地解決太陽能電池模型參數(shù)辨識的優(yōu)化問題。楊海天等[9]采用Levenberg-Marquardt(LM)算法對彈性支撐與本構(gòu)參數(shù)進行單一/組合反演，以及基于梯度法求解平面雙模問題確定本構(gòu)參數(shù)[10]。郭力等[11]針對非線性參數(shù)反演過程中靈敏度計算困難的問題，提出了靈敏度分析和彈塑性參數(shù)反演的一種新方法。Astroza等[12]提出了一種新的框架，將基于高保真度力學(xué)的非線性(滯后)有限元模型和非線性隨機濾波技術(shù)相結(jié)合，用于框架結(jié)構(gòu)中未知材料參數(shù)的估計。韓陽等[13]基于LM算法對向家壩巖石的三軸壓縮蠕變試驗曲線進行擬合及參數(shù)識別。Cui等[14]基于最小二乘法，進行瞬態(tài)熱傳導(dǎo)反問題中的參數(shù)辨識，他們還基于改進的LM算法[15]求解瞬態(tài)非線性熱傳導(dǎo)的多重邊界參數(shù)。薛齊文等[16]引入Bregman距離函數(shù)及其加權(quán)函數(shù)，應(yīng)用正則化技術(shù)，建立了一種非定常熱力耦合反問題的數(shù)值求解模式。韓雯雯等[17]基于共軛梯度法，利用圓管外壁面可測的有限溫度信息，同時反演出第三類邊界條件。大體上看，這些方法可以分為兩種：隨機法和梯度法。隨機法可以實現(xiàn)全局搜索，但是計算時間長、收斂慢、效率低，通常在實際工程中很少采用；而梯度法效率高、計算時間短，在研究中被廣泛使用。

綜上所述，考慮到ABAQUS有限元商用軟件具有強大的建模和計算能力，并且非常適于求解非線性問題，本文應(yīng)用ABAQUS軟件對冪強化材料的平面應(yīng)變問題進行數(shù)值求解，然后采用梯度法對反問題進行求解。在梯度法中，靈敏度矩陣的準(zhǔn)確確定非常重要，文獻[18]驗證了復(fù)變量求導(dǎo)法相對差分法的優(yōu)越性。然而在本文以前，能夠準(zhǔn)確確定靈敏度矩陣的復(fù)變量求導(dǎo)法很難與ABAQUS有限元軟件相結(jié)合，這是因為普通用戶無法獲得ABAQUS有限元軟件中的全部源程序代碼。為此，本文將ABAQUS中的用戶單元子程序(UEL)與復(fù)變量求導(dǎo)法相結(jié)合[19]，首次用于力學(xué)反問題的求解。既保證正問題求解的高效率與高精度，同時采用復(fù)變量求導(dǎo)法準(zhǔn)確確定靈敏度矩陣。最后，采用最小二乘法對反問題進行求解。據(jù)筆者目前所知，并未有其他文獻涉及本文研究內(nèi)容。

1 正問題的有限元推導(dǎo)

1.1 基本控制方程

本文研究的是冪強化彈塑性材料的平面應(yīng)變問題的反問題分析。在正問題中，模型的變形存在兩個階段：線彈性變形階段與塑性變形階段。

有限單元法中，單元需要滿足靜力平衡方程[20]：

式中：B為應(yīng)變位移矩陣；D為單元的應(yīng)力應(yīng)變矩陣；Km為單元剛度陣；U為單元節(jié)點的位移向量；F為單元的節(jié)點載荷列向量。

位移邊界條件：

式中：us、vs為單元節(jié)點的已知位移；為已知邊界的位移分量。

應(yīng)力邊界條件：

1.2 彈性變形階段

線彈性階段的本構(gòu)矩陣表達式為：

1.3 塑性變形階段

當(dāng)物體內(nèi)某一點的Von Mises應(yīng)力超過了材料的初始屈服應(yīng)力σ0，材料開始進入屈服狀態(tài)，物體的變形開始出現(xiàn)塑性變形。此時，材料的屈服應(yīng)力可以表示為：

式中：E1為材料的強化系數(shù)；p為等效塑性應(yīng)變；n為材料的強化指數(shù)，一般是介于0～1的正數(shù)。

由于材料本構(gòu)關(guān)系的非線性，在材料發(fā)生塑性變形時，應(yīng)力-應(yīng)變關(guān)系很難通過解析解求得。因此，本文利用Newtown迭代法求解單元的等效塑性應(yīng)變增量Δp，在屈服階段，單元要滿足屈服方程[21]：

將屈服方程進行泰勒展開，可得：

則可以得到：

式中，H=npn-1為當(dāng)前屈服過程中應(yīng)變強化系數(shù)。

在Newtown迭代完成后，需要對各種變量進行更新，更新的方法如式(10)～式(18)所示：

式中：為試探應(yīng)力σtr的應(yīng)力偏量；Δεp為塑性應(yīng)變增量；Δεe為彈性應(yīng)變增量；εp為塑性應(yīng)變；εe為彈性應(yīng)變；I為單位矩陣。

計算塑性變形下的Jacobian矩陣：

Jacobian矩陣的表達式為：

求得Jacobian矩陣后，將J中與z方向無關(guān)的元素提取出來，組成新的矩陣D。

將D代入式(1)中，即可求解當(dāng)前載荷下單元的靜力平衡方程。

2 反問題理論基礎(chǔ)

2.1 反問題中的目標(biāo)函數(shù)

在反問題中，材料的楊氏模量、泊松比、位移邊界條件與載荷邊界條件的大小是未知的，其他條件與正問題一樣給定，所需要的額外信息是物理量測量值，這些測量值可以通過實驗儀器測量所得。在本文中，這些測量信息是測點的位移值，通過數(shù)值模擬得到。

二維平面應(yīng)變的反問題的優(yōu)化目標(biāo)函數(shù)可以表示為：

式中：M為測量點的數(shù)量；N為待反演參數(shù)的個數(shù)；z=(z1,z2,z3,…,zN)T為待反演的參數(shù)向量；和ui分別為測點位移的測量值和計算值，i=1,2,…,M。

2.2 最小二乘法

全局余量矢量為R，其表達如式(23)所示。

式中，R為殘余矢量，代表u*和u的差距。

經(jīng)過推導(dǎo)，殘余矢量與參數(shù)增量的關(guān)系可以表示為如下矩陣形式。具體推導(dǎo)見文獻[14]。

式中，[?u/?z]為靈敏度矩陣，由位移對各反演參數(shù)的一階偏導(dǎo)數(shù)組成，詳見文獻[14]。

由于式(24)為病態(tài)方程組，應(yīng)用最小二乘法可以得到穩(wěn)定的解，表達式如下所示：

式(25)利用高斯消元法求解，可以得到更新的z的值。

式中，w為松馳因子[14]。

2.3 復(fù)變量求導(dǎo)法

1967 年，Lyness和Moler[22]提出了復(fù)變量求導(dǎo)法，該方法一般是將實數(shù)函數(shù)f(x)中的自變量x替換成復(fù)數(shù)變量x+ih，其中h取值很小，一般取10-20，通過泰勒級數(shù)展開后，f(x+ih)可以表示為：

由于h的取值很小，所以f(x)的一階導(dǎo)數(shù)可以表示如下：

式中，Im代表虛部。

在復(fù)變量求導(dǎo)法中，由于直接計算，所以一階偏導(dǎo)數(shù)的計算與空間步長大小無關(guān)，不涉及截斷誤差。因此，復(fù)變量求導(dǎo)法在一階偏導(dǎo)數(shù)的計算精度方面，總體上要優(yōu)于傳統(tǒng)的差分法[23]。本文的各靈敏度系數(shù)均采用復(fù)變量求導(dǎo)法獲得。

2.4 復(fù)變量求導(dǎo)法與ABAQUS有限元方法的結(jié)合

由于ABAQUS有限元軟件對復(fù)雜不規(guī)則或大型結(jié)構(gòu)的非線性力學(xué)問題的適應(yīng)性強，本文采用ABAQUS有限元軟件對正問題進行求解?，F(xiàn)有的ABAQUS軟件只能實現(xiàn)實數(shù)值的求解，而沒有對復(fù)數(shù)進行運算的求解器。為了采用復(fù)變量求導(dǎo)法獲得反演計算所需的靈敏度系數(shù)，本文采用ABAQUS的UEL子程序，這涉及實數(shù)變量與復(fù)數(shù)變量的形式轉(zhuǎn)換。

對于每個單元，需要將單元的每個節(jié)點變成一個實數(shù)節(jié)點和一個虛數(shù)節(jié)點相結(jié)合的復(fù)數(shù)節(jié)點，則二維復(fù)數(shù)單元的節(jié)點位移矢量具有如下形式：

式中：uRe=(ux,uy)T代表傳統(tǒng)單元節(jié)點的位移值；uIm代表虛數(shù)節(jié)點的位移。

通常可以將uIm表示為參數(shù)θ的擾動h對節(jié)點位移uRe的影響：

式中，u′為位移uRe對于參數(shù)θ的偏導(dǎo)數(shù)。

由于將傳統(tǒng)實單元轉(zhuǎn)換成了復(fù)數(shù)單元，那么相應(yīng)地，靜力平衡方程(1)的復(fù)數(shù)形式為：

式中，B*和D*分別是B和D的復(fù)數(shù)形式。

為了利用ABAQUS求解器求解這些方程，對于形如n×n的復(fù)數(shù)矩陣就必須化成2n×2n的實數(shù)矩陣，例如復(fù)數(shù)單元剛度陣的實數(shù)化如下所示：

靜力平衡方程的實數(shù)形式為：

將式(33)傳遞給ABAQUS求解，就能求出單元上虛值節(jié)點的位移，從而求出該點位移對被擾動參數(shù)的靈敏度系數(shù)，組成反演過程中的靈敏度矩陣。

2.5 反問題的收斂準(zhǔn)則

當(dāng)反演計算過程中目標(biāo)函數(shù)的值小于一個給定的足夠小的正數(shù)，或者前后兩次的目標(biāo)函數(shù)幾乎不變化時，迭代結(jié)束，可以表示為：

式中，ξ為給定的足夠小的正數(shù)，本文取ξ=10-5。

反問題的計算過程可以總結(jié)為以下幾步：

1) 假設(shè)待辨識的各參數(shù)的初始值為z0，假設(shè)初始迭代步k=0。

2) 將當(dāng)前的參數(shù)值zk傳遞進被調(diào)用的正問題求解程序，獲得該迭代步下的計算值uk。

3) 檢查目標(biāo)函數(shù)是否滿足收斂準(zhǔn)則，如果滿足，則迭代過程結(jié)束，計算停止；如果不滿足，繼續(xù)以下步驟。

4) 采用復(fù)變量求導(dǎo)法計算各參數(shù)的靈敏度系數(shù)[?u/?z]，測量值與計算值的差。通過最小二乘法求解得到Δz。

5) 計算得到zk+1。

6) 當(dāng)前迭代步為k=k+1，返回第(2)步，繼續(xù)計算。

3 數(shù)值算例1

3.1 復(fù)數(shù)有限元方法用于求解正問題的精度驗證

本文的算例由正問題和反問題組成，其中正問題是利用ABAQUS用戶單元子程序(UEL)求解冪強化彈塑性材料的非線性平面應(yīng)變問題。本文選取的算例模型為帶橢圓孔的六邊形，如圖1所示。其中模型的尺寸如圖1所示，橢圓孔位于模型中心，模型一共由2346個單元組成(圖2)，單元類型為8節(jié)點的完全積分平面應(yīng)變單元。

圖1 算例1的模型尺寸和測點位置 /mmFig.1 Model dimensions and measurement points in example 1

a、b邊為固支邊，c邊的初始位移條件為初始的載荷邊界條件為橢圓孔墻壁d受到了與外法線方向相反的均布載荷P|d=20MPa，材料的彈性模量為52.1×10MPa，泊松比為0.3，初始屈服應(yīng)力為100 MPa。

圖2 算例1模型的網(wǎng)格劃分形式Fig.2 The grid division form of model in example 1

材料的屈服應(yīng)力表達式為：

為了驗證分別使用ABAQUS的自帶CPE8單元與使用本文建立的自定義復(fù)變量UEL計算彈塑性問題的結(jié)果一致性，從模型上隨機提取10個節(jié)點，節(jié)點位置如圖1所示，節(jié)點的位移值如表1所示。

表1 兩種單元下節(jié)點位移的對比以及誤差Table 1 Comparison of node displacements and errors between the two elements

此外，本文還比較了應(yīng)力結(jié)果。隨機選取5個單元，單元中心積分點的y方向正應(yīng)力的結(jié)果如表2所示。

從表1和表2可以看出，無論是將節(jié)點位移還是將積分點應(yīng)力分量進行比較，利用ABAQUS傳統(tǒng)單元計算的結(jié)果與利用復(fù)變量UEL計算的結(jié)果非常接近，兩者計算結(jié)果相對誤差均小于0.05%。可見，將復(fù)變量引入到UEL中以后，對計算結(jié)果的實數(shù)部分沒有影響。

表2 單元中心積分點應(yīng)力的比較Table 2 Comparison of stresses at the central integral point of the element

采用應(yīng)力或者位移作為測量物理量進行反演，都是有效的，本文將位移作為測量物理量進行后文的反演。

3.2 本構(gòu)參數(shù)與邊界條件的反演

當(dāng)結(jié)構(gòu)處于塑性變形階段，材料的本構(gòu)參數(shù)很難確定，因此，反演方法提供了一種新手段。在本文的反問題中，彈性模量、泊松比以及邊界條件的值未知，需要進行反演，其他的計算條件均已知，與3.1節(jié)正問題模型相同。將正問題中一些節(jié)點x方向上的位移作為測量值，測量點選取的個數(shù)為100個，選取的方法為，以單元編號為參考，初始單元編號為40，此后每隔20個單元再選為下一個單元，依此類推，最后一個被選中的單元編號為2020。將被選中單元的第五個節(jié)點的x方向位移作為測點測量值。需要反演的各參數(shù)的初始值為：彈性模量E=2×105MPa，泊松比ν=0.2，邊界位移為[0.2,0.2]Tmm，均布載荷P為2 MPa。

圖3顯示了各參數(shù)的收斂過程，經(jīng)過13次迭代后，各參數(shù)分別收斂于真實值，收斂時目標(biāo)函數(shù)為1.35×10-6，可見本文提出的方法對于力學(xué)參數(shù)的反演具有較高的精度。同時，從圖4可以觀察到，優(yōu)化目標(biāo)函數(shù)的值總體上下降得比較快，這說明本文使用的算法不僅具有高精度，而且具有較高的效率。

圖3 各參數(shù)的收斂過程Fig.3 The convergence process of each parameter

圖4 目標(biāo)函數(shù)的變化過程Fig.4 Change process of objective function

3.3 反演初值的影響

為了研究參數(shù)反演初值的選取對反演過程的影響，選取4組參數(shù)初值，除初值的選取不同以外，其他條件均與3.2節(jié)相同。研究在不同初值情況下，迭代的收斂情況以及迭代次數(shù)，如表3所示。

表3 不同初值對收斂過程的影響Table 3 The influence of different initial values on the process of convergence

如表3顯示，收斂情況與初值有關(guān)。如果初值選取不當(dāng)，如第4組，可能會導(dǎo)致反演過程不收斂，而當(dāng)初值選取合適，如第1組～第3組，收斂迭代的總次數(shù)與初值沒有明顯的關(guān)系。所以在反演過程中，已知辨識參數(shù)的大致范圍對成功反演有幫助。

4 數(shù)值算例2

4.1 反演效率和精度分析

為了驗證本文反演算法對不同模型和不同參數(shù)的適應(yīng)性，考慮另一個模型，即一個帶孔的圓環(huán)模型，如圖5所示，圓環(huán)的內(nèi)半徑為40 mm，外半徑為100 mm，在圓環(huán)的周向方向上，存在著4個均勻分布的半徑為10 mm的圓孔。圓環(huán)內(nèi)壁受到P1=100 MPa的壓強，圓孔受到P2=30 MPa的壓強。由于模型的對稱性，在有限元的數(shù)值分析中將模型的1/4取出作為分析對象。

圖5 算例2模型的尺寸 /mmFig.5 Model size of example 2

有限元模型的網(wǎng)格劃分如圖6所示，模型一共由404個單元組成，單元類型為8節(jié)點的完全積分平面應(yīng)變單元。模型兩條直邊為對稱邊界條件。材料的楊氏模量55×10MPa，泊松比為0.35，彈塑性參數(shù)與算例1保持一致。測點為隨機選取的20個單元的第五個節(jié)點，測量物理量為該測點的x方向位移。

圖6 算例2模型的網(wǎng)格劃分形式Fig.6 The grid division form of model in example 2

在反問題中，待反演的參數(shù)為彈性模量E，泊松比ν，壓強P1和P2。初值的的選取E=2.5MPa ×106MPa，ν=0.25，P1=200 MPa，P2=50 MPa，松弛因子w取1。圖7為在無誤差測量值下的參數(shù)反演過程。

4.2 測量誤差的影響

在實驗過程和實際的工程中，位移的大小是通過儀器測量得到的。不可避免地，測量中會存在誤差，為了研究測量誤差對于最后反演結(jié)果的影響，定義相對誤差Erel為：

圖7 各參數(shù)的收斂過程Fig.7 The convergence process of each parameter

式中：z為待反演的參數(shù)；zestimated和zexact分別為在有/無測量誤差時參數(shù)反演的結(jié)果。

為了定量地分析測量誤差的大小對于反演參數(shù)的影響，假設(shè)測點的誤差分別為1%、5%、10%，所以測量的位移值可以表示為：

式中：umeasured代表位移的測量值；ζ為測量誤差；γ為-1～1的隨機數(shù)。

計算條件與4.1節(jié)相同，分別在20個測點上施加隨機誤差，測量誤差分別為1%、5%、10%。

表4 三種測量誤差對反演結(jié)果的影響Table 4 The influence of three kinds of measurement errors on inversion results

從表4可以看出，測量誤差對反演結(jié)果有影響，測量誤差越大，反演的參數(shù)與真實值的差距也越大。在一定的測量誤差內(nèi)，應(yīng)用本文的方法仍然可以反演出較為準(zhǔn)確的結(jié)果，同時可以發(fā)現(xiàn)Erel的值都小于相應(yīng)給定的測量誤差。

5 結(jié)論

本文首次將ABAQUS的用戶單元子程序(UEL)與復(fù)變量求導(dǎo)法結(jié)合起來，用于冪強化彈塑性材料的平面應(yīng)變反問題。該方法可以準(zhǔn)確求出反演所需要的計算值與靈敏度矩陣，并且驗證了復(fù)變量UEL求解力學(xué)問題的精度。在反問題中，以計算值與測量值之間的殘差作為優(yōu)化目標(biāo)函數(shù)，在上一步得到了靈敏度矩陣后，利用最小二乘法求解優(yōu)化參數(shù)的增量，進行多次迭代直到收斂。可以得出如下結(jié)論：

(1) ABAQUS中的UEL與復(fù)變量求導(dǎo)法相結(jié)合是確定求解反問題的梯度法中靈敏度矩陣的一種有效方法，所建立的復(fù)單元既不影響實數(shù)部分的計算，又能得到精確的靈敏度系數(shù)。

(2)本文建立的反問題求解算法是有效的，可以同時反演多個不相關(guān)的參數(shù)，包括材料本構(gòu)參數(shù)和邊界條件。

(3)反演參數(shù)初值影響反演結(jié)果：在合適的范圍內(nèi)，不同的初值對于總迭代次數(shù)和辨識結(jié)果沒有明顯影響；超過該范圍，可能會導(dǎo)致反演發(fā)散。

(4)當(dāng)存在一定的測量誤差時，本文建立的反演方法仍然具有較高的精度。

本文為解決復(fù)雜結(jié)構(gòu)非線性材料的力學(xué)反問題提供了一種新方法。