， ， Khemais Saanouni，

(1.大連理工大學 工業(yè)裝備結(jié)構(gòu)分析國家重點實驗室，大連 116024；2.ICD/LASMIS，University of Technology of Troyes，Troyes 10000)

1 引 言

一般塑性材料在受力變形達到拉伸極限強度后，會出現(xiàn)明顯的應變軟化直至材料最終斷裂。由材料微缺陷、損傷或溫度等引起的應變軟化現(xiàn)象導致了材料失穩(wěn)，進而失去承載能力，同時使得采用局部連續(xù)介質(zhì)損傷力學模型的準靜態(tài)問題控制方程橢圓性喪失，成為數(shù)學上病態(tài)的問題，且數(shù)值響應病態(tài)依賴于有限元網(wǎng)格，而收斂到不正確的沒有物理意義的有限元解上[1]。

針對上述問題，已展開大量研究，但有效的且根本性的補救措施是通過在本構(gòu)方程中引入特定尺寸內(nèi)物質(zhì)點間的相互作用，以保證控制方程的橢圓性，即正則化機制。當前引入正則化機制的途徑主要有四條。一是將經(jīng)典連續(xù)體模型擴展為率相關材料模型[2-4]，該粘性效應隱式包含了一種內(nèi)長度參數(shù)的作用，可以減弱數(shù)值結(jié)果對網(wǎng)格尺寸的依賴性。二是采用梯度理論模型[5,6]，最初是由Mindlin[7]提出的應變的二階梯度線彈性模型。Triantafyllidis等[8]提出在超彈性材料的應變能密度函數(shù)中引入二階變形梯度張量的依賴項，以保證控制方程的橢圓性。De Borst等[9]提出了廣義梯度塑性模型，并在屈服函數(shù)中引入了等效塑性應變的高階空間梯度。三是引入額外自由度以表征材料的非局部行為，Cosserat介質(zhì)[10]引入了旋轉(zhuǎn)自由度并相應產(chǎn)生了微曲率，進而引入了與微曲率能量共軛的對偶應力和具有特征長度意義的內(nèi)部長度參數(shù)。De Borst[11]推導了Mises屈服準則下的彈塑性Cosserat連續(xù)體模型。李錫夔等[12,13]提出了壓力相關彈塑性Cosserat連續(xù)體模型，考慮了非關聯(lián)Drucker-Prager屈服準則，并對其本構(gòu)參數(shù)的取值進行了分析。Micromorphic介質(zhì)最初由Eringen等[14-16]提出，視每個質(zhì)點為一個微小物體，并允許微小物體變形和做剛體轉(zhuǎn)動。Forest等[17,18]發(fā)展了包含新的微態(tài)自由度的基于廣義虛功原理和擴展的亥姆霍茲自由能的微態(tài)方法。同時，Eringen的微態(tài)連續(xù)體模型可視為微應變模型與Cosserat模型的結(jié)合體，分別對應于微態(tài)方法中微態(tài)變量取為局部應變和旋轉(zhuǎn)張量時的模型。四是采用完全非局部的連續(xù)體理論[19]，其基本思想不再基于局部行為假設，而是假定整個域內(nèi)的物質(zhì)點均處于相互作用的狀態(tài)。Ba?ant等[1,20,21]以此建立了針對控制應變硬化所有變量的積分形式的正則化機制，并應用于混凝土的應變軟化分析中。

2 本構(gòu)模型

從虛功原理出發(fā)，通過引入與損傷內(nèi)變量相關聯(lián)的微態(tài)韌性損傷，以反映損傷演化過程中的非局部效應；在熱力學定律和克勞修斯不等式的幫助下，通過亥姆霍茲自由能求解得到各狀態(tài)變量的表達式；應用廣義法向法則根據(jù)彈塑性耗散勢能獲得耗散內(nèi)變量的演化方程。

2.1 平衡方程

(1)

相應地，外力虛功和慣性力虛功可表示為

(2)

應用廣義的虛功原理Pint+Pext=Pa,可以得到經(jīng)典的位移平衡方程及邊界條件：

(3)

額外的與新引入的微態(tài)損傷相關的廣義平衡方程及邊界條件：

(4)

2.2 狀態(tài)方程

微分形式的熱力學第一定律可擴展為

(5)

考慮等溫過程，結(jié)合式(5)可得到Clausius-Duhem不等式為

(6)

為方便起見，在缺少實驗支撐的前提下，本文假設微態(tài)損傷變量引起的耗散可以忽略不計，則通過將亥姆霍茲自由能函數(shù)代入不等式(6)中，可得狀態(tài)方程為

(7)

與經(jīng)典彈塑性損傷模型一致的機械耗散項為

(8)

式中R和Xi j分別為與各向同性硬化應變和隨動硬化應變相共軛的應力，Y為與損傷共軛的應力，也稱為損傷能量釋放率。

考慮如下形式的亥姆霍茲自由能：

(9)

通過狀態(tài)方程(7)，可得各廣義應力變量的表達式為

σ=(1-d)Λεe，R=(1-d)Qr，X=(1-d)Cα

(10)

損傷能量釋放率可進一步分解為兩部分之和：

(11)

式中Yloc>0是經(jīng)典彈塑性損傷模型中的損傷能量釋放率，Yn l來自于損傷和微態(tài)損傷的耦合作用，同時Y恒為正。

2.3 演化方程

采用式(12)定義的Von Mises屈服函數(shù)和耗散勢能。

(12)

采用廣義法向非關聯(lián)流動法則得到耗散內(nèi)變量的演化方程為

(13)

(14)

式中Hp為彈塑性硬化模量,

(15)

進而可求得應變率張量

(16)

其中，四階彈塑性切線模量可表示為

Y(nΛ)?[σ/(1-d)]}

(17)

3 有限元離散

3.1 應變空間的平衡方程

利用式(10)的狀態(tài)方程，將式(4)的微態(tài)損傷平衡方程變換到微態(tài)損傷應變空間，并忽略其廣義體力與接觸力，可以得到簡潔形式的偏微分方程和邊界條件：

(18)

式(18)與Engelen等在隱式的梯度損傷模型中提出的亥姆霍茲方程是一致的。

3.2 弱形式的初始邊界值問題

采用典型的加權(quán)余量法對式(3，18)的偏微分方程和Neumann邊界條件進行弱化，得到其對應的積分形式的初始邊界值問題：

(19)

3.3 有限元空間離散

采用基于Hu-Washizu變分原理的假設應變四邊形單元以減少沙漏和閉鎖現(xiàn)象的影響，則單元離散的弱形式初始邊界值問題為

(20)

式中 上標e代表參考單元。與位移相關的阻尼質(zhì)量矩陣和內(nèi)外力向量分別為

(21)

與微態(tài)損傷對應的阻尼質(zhì)量矩陣和內(nèi)外力向量分別為

(22)

由于動力顯式求解過程是條件穩(wěn)定的，須確定穩(wěn)定的最大時間步長Δt≤min(Le/Cd),其中，Le為單元的特征長度，Cd為等效的材料聲速，可通過Lame常數(shù)近似確定。

3.4 本構(gòu)積分過程

為了計算內(nèi)外力向量以求解式(20)的代數(shù)方程組，需要根據(jù)已知tn時刻的各狀態(tài)變量得到tn +1時刻的應力張量、硬化變量、局部損傷和微態(tài)損傷因子。本文采用完全隱式的彈性預測與塑性修正的Simo圖形返回算法更新該非關聯(lián)彈塑性模型的狀態(tài)變量。

彈性預測

(23)

判斷屈服函數(shù)：

(24)

塑性修正

在tn +1時刻的塑性應變與應力張量：

(25)

損傷及損傷能量釋放率：

(26)

各項同性硬化應變及應力：

(27)

隨動硬化應變及應力：

(28)

Δλ>0時，滿足tn +1時刻的屈服函數(shù)：

(29)

4 算例與討論

4.1 DP1000材料的拉伸實驗

將對偶雙相鋼DP1000材料在Zwick 250拉伸機中以0.1 mm/s拉伸速度進行準靜態(tài)單向拉伸實驗。樣件厚度為1.5 mm，其幾何形狀尺寸如圖1所示。擬合的材料參數(shù)列入表1。

采用局部彈塑性損傷模型計算三個不同網(wǎng)格尺寸(0.4 mm，0.8 mm和1.6 mm)下的結(jié)構(gòu)位移-荷載曲線,如圖2所示?？梢钥闯?，網(wǎng)格尺寸為1.6 mm的系統(tǒng)響應曲線與實驗的位移-荷載曲線吻合，拉伸試件斷裂時的位移約為9.0 mm。同時可以觀察到明顯的網(wǎng)格依賴性，隨著網(wǎng)格尺寸的減小(0.8 mm和0.4 mm)，局部彈塑性損傷模型計算的系統(tǒng)響應曲線下降時刻越早，斷裂位移相應越小(分別對應為8.5 mm和8.2 mm附近)。隨著網(wǎng)格尺寸以1/2的比率遞減，相鄰網(wǎng)格尺寸所對應的斷裂位移間的偏差也越來越小(約為0.5 mm和0.3 mm)。

圖1 DP1000單向拉伸實驗樣件

Fig.1 Geometry of the DP1000 uniaxial tension specimen

圖2 拉伸實驗局部損傷模型下的位移-荷載曲線

Fig.2 Displacement-force curves of the local damage model of the uniaxial tensile test

圖4分別給出了網(wǎng)格尺寸為0.4 mm的結(jié)構(gòu)拉伸至剛出現(xiàn)破壞單元時，局部和微態(tài)損傷模型計算的局部損傷分布圖?？梢钥闯?，微態(tài)損傷模型的損傷集中帶比局部損傷模型預測的損傷集中帶更寬，損傷分布更加平滑。圖5的曲線分布趨勢與之類似。自峰值處(約18.0 mm)向兩側(cè)延伸方向，與局部損傷模型計算的有明顯集中現(xiàn)象的局部損傷分布相比，在微態(tài)損傷模型中的局部損傷變量，由于微態(tài)損傷的正則化作用，其分布更加平滑且均有不同程度的抬升(10 mm～25 mm)。

圖3 拉伸實驗微態(tài)損傷模型下的位移-載荷曲線

Fig.3 Displacement-force curves of the micromorphic damage model of the uniaxial tensile test

圖4 0.4 mm單元的損傷分布

Fig.4 Distributions of damage for mesh size is 0.4 mm

4.2 十字形零件的沖壓實驗

圖6為十字形零件的沖壓成形模具的幾何信息，材料為1.5 mm厚的DP1000對偶雙相鋼，壓邊力為400 kN，摩擦系數(shù)約為0.15?？紤]微裂紋在壓縮狀態(tài)下的閉合效應，式(11)的局部損傷能量釋放率修改為式(30)，并取h的值為0.25。

(30)

圖5 0.4 mm單元的損傷變量沿軸向的分布

Fig.5 Distribution of damage along thexaxial for mesh size is 0.4 mm

圖6 十字形零件模具的幾何信息[22]

Fig.6 Geometric information of the cross-section drawing part[22]

圖7 十字形零件成形局部損傷模型的位移-載荷曲線

Fig.7 Displacement-force curves of the local damage model of the cross-section drawing part

圖8 十字形零件成形微態(tài)損傷模型的位移-載荷曲線

Fig.8 Displacement-force curves of the micromorphic damage model of the cross-section drawing part

5 結(jié) 論

本文根據(jù)Forest的微態(tài)方法，提出引入微態(tài)損傷因子及其一階梯度，以考慮應變集中的非局部效應，建立了新的適用于有限變形框架的熱力學一致性的微態(tài)損傷彈塑性本構(gòu)模型。通過在擴展的亥姆霍茲自由能中定義局部損傷與微態(tài)損傷的耦合方式，進而得到了包含材料內(nèi)部特征尺度的類似于隱式梯度理論中補充的亥姆霍茲方程，同時獲得了新的包含局部損傷和微態(tài)損傷的損傷能量釋放率。由微態(tài)損傷變量對損傷能量釋放率進行正則化作用，進而影響了局部損傷的演化過程，最終達到了一致性的系統(tǒng)位移-載荷響應。對DP1000材料的單向拉伸和十字形零件的沖壓形實驗進行數(shù)值分析，結(jié)果表明，該微態(tài)損傷模型可以有效地避免材料應變軟化行為導致的數(shù)值響應依賴于網(wǎng)格的問題，并弱化損傷變量的局部集中現(xiàn)象。并且該模型可進一步擴展包含混合硬化和塑性應變等非局部效應，以完善對材料應變軟化行為的描述。