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

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

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

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

1 傳統損傷算法中漂移誤差的成因

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

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

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

外力由荷載因子λ與基準荷載矢量f的乘積形式表示。當結構的外力與內力處于非平衡狀態時，則體系內存在殘余力矢量r為：

全量應力張量σ可由損傷本構表達式為[8]：

式中：ε為全量應變張量；CD為損傷本構張量。

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

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

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

在實際應用過程中，常采用增量的傷本構[10]：

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

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

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

2 SSALM 法

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

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

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

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

式(14)可簡寫為：

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

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

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

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

式中：j為迭代次數；i為增量步。

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

式中：l為規定弧長。

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

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

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

3 試驗算例

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

3.1 損傷桿模型的單軸拉伸

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

式中：ε為單軸拉伸應變；κ=3.50×10?6為損傷參數。

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

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

3.2 損傷梁模型的純彎曲

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

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

式中：為損傷力學中所定義的有效應力，κ=5.0×105Pa 為損傷參數。

該純彎梁的解析解在文獻[17]中已經給出。

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

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

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

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

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

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

4 傳統法和SSALM 法對比分析

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

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

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

傳統增量迭代法與SSALM 的相互關系，如圖6 所示。虛線和空心點表示傳統增量迭代法的結果和路徑；點劃線和實心點表示SSALM 的結果和路徑。

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

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

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

5 結論

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

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

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