損傷分析的精確增量迭代算法

趙冰，宋文浩，陳健，彭暉，彭旭龍

(長沙理工大學(xué) 土木工程學(xué)院，湖南 長沙 410004)

損傷非線性分析已應(yīng)用于眾多領(lǐng)域，如：材料疲勞壽命的預(yù)測、構(gòu)件裂紋的萌生與擴展及巖土強度的分析等。對于損傷力學(xué)，其結(jié)構(gòu)應(yīng)力依賴于變形歷史，全量理論已不再適用。必須采用增量理論，跟蹤結(jié)構(gòu)的平衡路徑。因此，增量平衡方程在非線性分析中，得到了廣泛的應(yīng)用。唐雪松[1?2]等人采用損傷力學(xué)-有限元法，研究了加載次序?qū)ζ趬勖挠绊懽饔?，并對標準疲勞試驗?gòu)件進行了疲勞壽命預(yù)測。張杰毅[3]等人采用損傷力學(xué)-有限元法，預(yù)估了球軸承的接觸疲勞壽命。李聰成[4]等人采用損傷力學(xué)有限元法，分析了蠕變疲勞交互作用下，影響裂紋萌生壽命的因素。張我華[5]等人將損傷力學(xué)引入巖石類介質(zhì)的非線性本構(gòu)模型，分別建立了脆性動力和彈黏塑性動力的損傷模型，并編制相應(yīng)的有限元程序，分析巖石的動力學(xué)行為。損傷非線性分析已運用到醫(yī)學(xué)領(lǐng)域。Taylor[6]等人采用損傷力學(xué)-有限元法，模擬了人體骨骼在人造環(huán)境內(nèi)的疲勞行為，應(yīng)用于人體腦骨骼樣本的疲勞失效次數(shù)、彈性模量的退化及骨骼內(nèi)永久應(yīng)變的累積狀況?？紤]這些傳統(tǒng)損傷有限元分析中的增量迭代算法，會產(chǎn)生迭代無法消除的漂移誤差。作者擬提出一種全新的增量迭代方法：正割剛度-附加荷載法(Secant stiffness– additional load method,簡稱為 SSALM)，以校正傳統(tǒng)增量迭代法的全部漂移誤差，尋求其真實的結(jié)構(gòu)響應(yīng)，以提高算法的穩(wěn)定性。

1 傳統(tǒng)損傷算法中漂移誤差的成因

損傷非線性有限元格式可簡單分為兩類：

1) 第一類是利用力的平衡條件(內(nèi)力與外力平衡)或虛功原理得到的全量形式[7]：

式中：KS(u)為結(jié)構(gòu)在某種狀態(tài)(λ,u)下的正割剛度矩陣；u為結(jié)構(gòu)位移矢量；B為應(yīng)變矩陣；σ為全量應(yīng)力張量；Ωe為單元體積；G為組裝矩陣。

外力由荷載因子λ與基準荷載矢量f的乘積形式表示。當結(jié)構(gòu)的外力與內(nèi)力處于非平衡狀態(tài)時，則體系內(nèi)存在殘余力矢量r為：

全量應(yīng)力張量σ可由損傷本構(gòu)表達式為[8]：

式中：ε為全量應(yīng)變張量；CD為損傷本構(gòu)張量。

將式(3)代入式(1)，全量平衡方程為：

2) 第二類為增量平衡方程，其有限元基本式[9]為：

式中：Δσ為增量應(yīng)力張量；Δu為增量位移矢量；Δλ為增量荷載因子；KT(u)為切線剛度矩陣。

在實際應(yīng)用過程中，常采用增量的傷本構(gòu)[10]：

對式(6)進行線性化處理[11]，可得：

將式(6)和(9)稱為傳統(tǒng)增量迭代法，而式(9)是傳統(tǒng)增量平衡方程。該方法已應(yīng)用于許多領(lǐng)域，如：疲勞壽命預(yù)測[1?3]和高階梯度損傷理論[12]等。對比式(5)，(9)可以看出，式(9)中的傳統(tǒng)方法是采用KS(u)代替了式(5)的KT(u)。然而，切線剛度矩陣[13]由式(2)導(dǎo)出：

單純地采用增量法，必然引起解的漂移[13]。如果沒有模型誤差，這種解的漂移通常是通過在增量步內(nèi)設(shè)置迭代來校正的[13]。因為傳統(tǒng)增量迭代法忽略掉了。所以，導(dǎo)致增量平衡方程所表征的平衡路徑失真，而引起漂移誤差，并且不能通過迭代技術(shù)消除。采用KT(u)表達的增量迭代算法，在峰值點附近，常常存在剛度矩陣病態(tài)化、奇異化及收斂性差[14]等問題，其穩(wěn)定性在很大程度上取決于KT(u)的條件數(shù)。

2 SSALM 法

對式(2)進行全微分展開，可得：

將式(10)代入式(11)，可得：

對于損傷力學(xué)，正割剛度矩陣KS(u)的表達式為：

將式(13)代入式(12)，可得：

式(14)可簡寫為：

式中：i為當前平衡狀態(tài)；i-1 為上一個平衡狀態(tài)。

其中，λi>λi?1。利用式(17)減去式(18)，則可得增量平衡方程為：

式(15)，(21)具有相同的形式，均由KS(u)和Df表達。因此，正割剛度-附加荷載法的增量平衡方程，已從2 個不同的角度完成推導(dǎo)。

為了提高SSALM 法的求解精度，可以在增量步內(nèi)加入平衡校正的迭代步[13]。通過將式(15)或式(21)與迭代技術(shù)相結(jié)合，得到新的增量迭代方程：

式中：j為迭代次數(shù)；i為增量步。

由于引入了荷載因子λ，式(22)多出的一個未知量需要補充一個約束方程。該弧長形式的約束方程[15]為：

式中：l為規(guī)定弧長。

Crisfield[16]等人提出過其他類型的約束方程。本研究采用最經(jīng)典的弧長迭代法。

SSALM 法由KS(u)來表達的，避免了在峰值點附近KT(u)病態(tài)化的現(xiàn)象。

將SSALM 法編入FORTRAN 有限元程序中，每個迭代步的收斂條件，采用位移判斷準則。

3 試驗算例

依托損傷桿的單軸拉伸和損傷梁的純彎曲2 個試驗算例，驗證SSALM 法的精確性。

3.1 損傷桿模型的單軸拉伸

采用8 節(jié)點平面應(yīng)變單元，建立損傷桿單軸拉伸的有限元模型，共12 609 個節(jié)點和4 096 個單元。桿左端完全固定，桿右端施加平行于桿軸的均布荷載q，如圖1 所示。模型材料參數(shù)為：彈性模量E0=21.0 GPa，泊松比υ=0。采用與應(yīng)變相關(guān)的損傷演化律：

式中：ε為單軸拉伸應(yīng)變；κ=3.50×10?6為損傷參數(shù)。

圖1 單軸拉伸損傷桿的有限元網(wǎng)格和邊界條件Fig. 1 FEM mesh and boundary conditions of the damage bar

損傷桿的單軸拉伸有解析解。分別使用傳統(tǒng)方法、SSALM 法和解析法，求出損傷桿的荷載-位移曲線，如圖2 所示。從圖2 中可以看出，傳統(tǒng)方法出現(xiàn)了明顯的漂移誤差，即使改變增量步的步長，漂移誤差并未改善。而采用SSALM 法，所模擬的結(jié)果，與解析解非常吻合。表明：傳統(tǒng)增量迭代法發(fā)生了無法迭代消除的漂移誤差，而SSALM 法可以跟蹤到結(jié)構(gòu)的真實響應(yīng)。

3.2 損傷梁模型的純彎曲

圖2 單軸拉伸損傷桿模型的載荷-位移曲線Fig. 2 Load-displacement curves of the uniaxial tension of damage bar

懸臂梁的純彎曲有限元模型如圖3 所示，采用與損傷桿單軸拉伸的有限元模型相同的網(wǎng)格。梁的左端完全固定，梁的右端施加了大小為M的彎矩。模型材料參數(shù)為：彈性模量E0=21.0 GPa，泊松比υ=0。采用了與應(yīng)力相關(guān)的損傷演化律。損傷變量定義(損傷演化方程式)為：

式中：為損傷力學(xué)中所定義的有效應(yīng)力，κ=5.0×105Pa 為損傷參數(shù)。

該純彎梁的解析解在文獻[17]中已經(jīng)給出。

圖3 純彎曲損傷梁的有限元網(wǎng)格和邊界條件Fig. 3 FEM mesh and boundary conditions of the damage cantilever beam

考慮到圣維南原理(遠離右端的橫截面，受邊界效應(yīng)的影響不大)和純彎曲的變形特點，遠離右端的任意截面，具有相同的損傷分布。本算例取中間橫截面，分析其損傷分布和變形。

根據(jù)SSALM 法和解析法，得到梁的中間橫截面的損傷分布，如圖4 所示。從圖4 中可以看出，受壓區(qū)無損傷發(fā)生，受拉區(qū)的損傷發(fā)展與梁的深度d成線性關(guān)系。梁的中性層隨著外加彎矩M的增大，逐漸下移。最大損傷Dmax出現(xiàn)在梁頂。這些規(guī)律與式(25)一致，且與實際結(jié)果相吻合。采用SSALM 法得到的結(jié)果與解析解非常吻合。

圖4 梁中間橫截面損傷分布Fig. 4 The distribution of damage in the cross section of the beam

圖5 彎矩M 與最大損傷Dmax 的曲線關(guān)系Fig. 5 The bending moment M-maximum damage Dmax curves

梁頂最大損傷Dmax與外部施加彎矩M之間的關(guān)系，如圖5 所示。從圖5 中可以看出，傳統(tǒng)增量迭代法追蹤的平衡路徑與解析解偏離的非常嚴重。而SSALM 法精確地捕捉到了結(jié)構(gòu)的真實平衡路徑，不會出現(xiàn)迭代無法消除的漂移誤差。

4 傳統(tǒng)法和SSALM 法對比分析

通過對比分析傳統(tǒng)增量迭代法和SSALM 法，分析引起傳統(tǒng)增量迭代法誤差的原因和SSALM 法如何修正該誤差。

式(6)的增量應(yīng)力-應(yīng)變關(guān)系，由式(3)得到表達式(26)或(27)為：

比較式(7)與式(27)可知，ε·ΔCD項在式(7)中被忽略。傳統(tǒng)的損傷有限元法，由于忽略掉了 ε·ΔCD項，其精度必然會下降。而在SSALM 中，被忽略掉的ε·ΔCD項，被轉(zhuǎn)化成附加荷載f D，引入到增量平衡方程之中，克服了傳統(tǒng)方法所產(chǎn)生的漂移誤差。

傳統(tǒng)增量迭代法與SSALM 的相互關(guān)系，如圖6 所示。虛線和空心點表示傳統(tǒng)增量迭代法的結(jié)果和路徑；點劃線和實心點表示SSALM 的結(jié)果和路徑。

在增量步內(nèi)利用迭代技術(shù)，可消除部分漂移誤差，將圖6 中空心方點校正到空心圓點。

圖6 傳統(tǒng)法和SSALM 法的誤差對比分析Fig. 6 The error contrastive analysis in two different incremental-iterative methods

由式(29)可以知，傳統(tǒng)增量迭代法中，由于忽略 ε·ΔCD，產(chǎn)生了漂移誤差(模型誤差)，該漂移誤差可以通過引入fD有效消除，使得損傷有限元分析更貼近真實解。表明：SSALM 法將KT對應(yīng)的增量位移，分解為KS對應(yīng)的增量位移和fD對應(yīng)的增量位移2 個部分。

5 結(jié)論

1) 在進行損傷非線性分析時，傳統(tǒng)的增量迭代算法，會產(chǎn)生迭代無法消除的漂移誤差。

2) SSALM 法可以精確地捕捉到結(jié)構(gòu)的真實平衡路徑，不會出現(xiàn)迭代無法消除的漂移誤差。

3) SSALM 法將KT對應(yīng)的增量位移分解為KS對應(yīng)的增量位移和fD對應(yīng)的增量位移2 個部分。