考慮混凝土塑性耗散的CDM-XFEM 裂縫計算方法1)

金 浩 *, 余 朔

* (東南大學交通學院,南京 211189)

? (同濟大學道路與交通工程教育部重點實驗室,上海 201804)

引言

混凝土結構是建筑、交通、水利等領域的重要基礎材料,其可塑性高、耐久性好是受廣泛應用的主要原因.但混凝土抗拉強度遠低于抗壓強度[1],結構在服役期間受外界載荷的影響容易產生裂縫[2],裂縫的出現不僅會降低結構剛度,還為外部侵蝕介質的侵入提供了快捷通道,從而加速內部鋼筋銹蝕[3-4]、降低結構承載力.因此,采用準確的計算方法預測混凝土的裂縫發展是治理裂縫的基本前提,也是保障結構安全的重要手段.

基于非線性斷裂力學,混凝土開裂大致分為以下三個過程[5-7]:(1)當材料拉伸應力接近抗拉強度時,在斷裂過程區會形成密集分布的微裂紋,載荷-變形的關系不再具有線性斜率,但材料的宏觀響應仍然保持穩定;(2)由于載荷增大引起微裂紋的合并和交叉,以及材料基體中骨料的脫黏,材料的宏觀響應變得不穩定,載荷-變形的關系位于軟化段中,局部范圍出現較大變形,在應變場中形成了弱不連續性;(3)隨著載荷繼續增大,斷裂過程區黏聚力逐漸減小至零,最終形成一個無黏聚應力裂紋的強不連續區域.

為模擬混凝土開裂的上述過程,出現了較多的數值計算方法,目前主流方法有三大類:連續損傷力學方法(continuou damage method,CDM)、擴展有限單元法(mechanics-extended finite element method,XFEM)及兩者相結合的方法(CDM-XFEM).其中CDM 方法是假定當等效應力達到抗力準則時,開始發生初始斷裂,即在損傷單元中,應力-應變關系會被有效應力-應變關系所取代,因此,CDM 方法只更新裂紋擴展過程中材料的本構關系,其單元網格保持不變[8-9].如Vilppo 等[10]從熱力學的角度提出了各向異性損傷的本構模型,分析了混凝土單軸受拉及受壓情況下的破壞規律.Poliotti 等[11]提出了考慮剪脹參數變化的塑性損傷模型,分析了不同強度的混凝土受剪切作用下的受力狀態.Pereira 等[12]提出了一種與有效率相關的非局部損傷模型,模擬了單切口混凝土受拉出現分叉裂縫的情況.而XFEM 方法是擴充帶有不連續性質的形函數來代表計算區域內的間斷,在計算過程中,不連續場的描述完全獨立于網格邊界[13-14].如Agathos 等[15]在XFEM 的基礎上,提出采用向量水平集的算法更新裂縫面的演化,研究了三點受彎作用下含單切口混凝土的裂縫形態.Schatzer 等[16]將顯-隱式結合的水平集函數融入到XFEM 方法中,分析了含邊緣裂口混凝土的裂縫張開度.為提高裂縫的計算效率及穩定性,Haghani等[17]開發了基于XFEM 和α時間積分結合的混凝土裂縫擴展程序,研究了振動作用下混凝土大壩的裂縫擴展規律.

相比而言,CDM 和XFEM 這兩種方法各有利弊.CDM 方法能夠很好地描述裂縫擴展的第一個階段,但該方法不能描述離散的開裂面,同時存在網格誘導偏差及虛假應力傳遞的弊端[18-20].XFEM 方法雖然能夠很好的描述宏觀裂紋的擴展,但不能較好地描述第一階段中密集分布的微裂紋,計算出的裂紋分布與實際差異較大[21].因此,部分學者提出將CDM 及XFEM 相結合的方法,該方法基于能量等效原理,建立CDM 與XFEM 黏聚裂縫的能量轉化關系,當損傷達到臨界值時實現轉換[22].如Comi 等[23]采用CDM-XFEM 方法研究了含雙切口的試件受豎向拉伸作用下的開裂情況,Jin 等[24]采用該方法研究了含單切口混凝土梁受彎作用下的裂縫形態,Pandey等[25]采用該方法研究了高頻疲勞載荷作用下混凝土裂縫的擴展規律.

采用CDM-XFEM 準確計算裂縫的關鍵在于實現CDM 與XFEM 之間能量轉化的平衡性.首先要采用合適的牽引-分離函數模型描述裂縫張開的變化模式,如直線型[26]、雙折線型[27]、指數型函數[28]等,其次是選取合適的能量轉化點,將宏觀裂縫出現時損傷模型中的剩余能量全部轉移到黏結裂縫模型中.為此,Roth 等[29]假定能量轉化時,兩者能量等效的分割線為垂線,且轉化點的應力相等,以此建立了能量平衡方程并求解了牽引-分離函數.Bobinski等[30]假定開裂時僅出現不可恢復變形,不考慮材料剛度發生退化,并自定義能量轉化時的臨界位移,從而構建了兩者能量的轉化方程.Wang 等[31]基于彈性損傷理論,假定材料發生開裂時彈性能全部回彈,從而建立了對應的能量平衡方程.

然而,相關試驗研究證明[32-33],在實際情況下混凝土出現宏觀裂紋時既存在剛度退化現象,也存在不可逆變形的特征.因此,本文考慮CDM 與XFEM之間斷裂能轉化時的塑性耗散,重新構建能量轉化方程,采用廣義逆的最小二乘法求解臨界位移及牽引-分離函數的系數,最終形成考慮混凝土塑性耗散的CDM-XFEM 計算方法.通過與雙切口混凝土受剪拉開裂的試驗進行對比驗證,結果表明采用考慮塑性耗散CDM-XFEM 方法能夠有效地預測混凝土裂縫的發展規律.

1 考慮塑性耗散的CDM-XFEM 計算方法

1.1 CDM-XFEM 方法的計算過程

CDM-XFEM 裂縫計算方法的過程為:當損傷因子d小于臨界值時,基于連續損傷力學描述微裂縫的發展,當損傷因子d大于臨界值dcrit時,將損傷模型所需耗散的剩余能量轉移到黏結裂縫模型中,采用擴展有限元單元法描述宏觀裂縫的張開規律,如圖1 所示,其計算過程主要包含4 個部分.

圖1 CDM-XFEM 方法計算過程示意圖Fig.1 Calculation process CDM-XFEM method

1.2 混凝土受拉損傷及塑性應變的計算式

受拉彈性損傷因子D的演化準則的計算式[35]為

式中,R為損傷演化參數,其計算為[29]

其中lrve為單元特征長度,對于C3D8R 單元,一般取值0.02 m[36].r0為混凝土的初始剛度,計算式為

對于單次加載-卸載后出現塑性損傷的材料,Hatzigeorgioiu 等[37]根據試驗結果將彈性損傷因子D修正為Dk,則變化后的彈性模量Ek表示為

式中,E為初始彈性模量,則總應變ε可由彈性應變εe和塑性應變εp表示為

對于混凝土材料,Hatzigeorgioiu 等[37]通過數值計算給出了最優值k=2,本文采用該數值,代入式(6)可得出塑性應變的簡化計算公式

1.3 黏結裂縫的牽引-分離函數形式

目前CDM-XFEM 計算模型中均不考慮黏結裂縫單元的切向牽引力ts,僅考慮黏結裂縫的法向牽引力tn,在本文研究中同樣只考慮裂縫間的法向牽引力,其牽引-分離函數體現的是裂縫張開量與裂縫間牽引力的相互關系,常用的牽引-分離函數主要包括:線型、雙折線型,多項式型及指數型,其計算為

由式(8)可知,第1 種函數求取的參數最少,但相關研究表明,采用線性函數表達牽引-分離法則計算得出的裂縫與實際裂縫相差較大[38];第2 種和第3 種函數需求取的未知變量數分別為3 和4,需引入較多的邊界條件才能求取相應系數[39];第4 種采用了指數型函數來表示,既體現了裂縫張開與牽引力之間的非線性關系,求解的未知變量也較少,因此,在本文的研究中也選擇指數型函數來表達黏結裂縫的牽引-分離關系.

1.4 能量轉化的基本假定條件

本文假定能量轉化時 (D=Dcr) CDM 模型中單元的最大拉應力等于XFEM 中裂縫張開至ucr時的牽引力,即 σ(εcr)=t(ucr),此時產生的塑性應變為εds,損傷消耗的能量 (當＞r0時開始出現損傷) 等于黏結裂縫張開至ucr的能量,即圖2 中藍色的面積等于圖3 中藍色的面積;同時,損傷模型的剩余能量均全部轉移到XFEM 粘結裂縫模型中,即圖2 中綠色和黃色的累計面積等于圖3 中黃色的面積,其中‖u‖代表裂縫的總張開量.

圖2 CDM 的應力-應變曲線Fig.2 Stress-strain curve of CDM

圖3 XFEM 粘結裂縫的牽引-分離曲線Fig.3 Traction-separation curve of XFEM cohesive crack

1.5 能量轉化方程

根據圖3 所示,基于能量等效原理,可得出關系式(9)和式(10)

式(9)中,能量Ψ1可由r0到εcr的積分面積減去能量Ψ2的面積,即

根據應變等效原則[29],即無損傷試件產生的應變ε等價于含損傷試件產生的應變,則等效應力可表示為

將式(2)、式(4) 和式(12) 代入到式(11) 中,則r0到εcr的積分面積可表示為

式(10)中,能量Ψ2可由為直角三角形面積表示,即

其中,εcr-εds可表示為

將式(15)代入式(14),則能量Ψ2的計算式表示為

為簡化表示,將損傷因子Dcr表示為關于εcr的函數A(εcr),則式(16)為

將式(13)和式(17)代入到式(11),可得到能量Ψ1的計算公式

式中,Γ(x)為伽馬函數,裂縫的牽引分離函數形式采用指數型函數,即式(8)的tn4,則能量Ψ4可表示為

將式(18)和式(19)代入到式(9),則第1 個能量等效方程可表示為

同理,能量Ψ3和能量Ψ5的方程可表示為

將式(17)、式(21) 和式(22) 代入到式(10),則第2 個能量等效方程可表示為

根據1.4 節提出的應力等效假定條件,除了滿足能量等效,轉換時的應力條件也應該相等,即

聯立方程式(20)、式(23)及式(24),最終求出參數α,β和ucr.

2 CDM-XFEM 的數值求解算法

2.1 數值求解流程

本文采用ABAQUS 軟件的用戶自定義擴展模塊,編寫相應的用戶單元子程序 UEL(user defined element).計算流程總體分3 個部分,即

(1) 計算八節點六面體單元的應力及損傷因子,以0.7 作為臨界損傷標準[29],判定模型是否達到能量轉化條件.

(2) 求解能量轉化的臨界位移,計算牽引-分離函數的系數,確定裂縫張開量.

(3) 基于主應力確定裂縫的發展方向,即裂縫的發展方向垂直于最大主應力方向,更新裂縫尖端及裂縫面的水平集函數,具體如圖4 所示.

圖4 考慮混凝土塑性耗散的CDM-XFEM 方法計算流程Fig.4 Calculation flow of CDM-XFEM method considered concrete plastic dissipation

2.2 能量轉化系數求解方法

對于牽引-分離函數參數α,β和臨界位移ucr,采用廣義逆的最小二乘法進行求解[40],步驟為:先將能量等效方程及應力等效方程轉為3 個待求的非線性方程,即式(25),然后對臨界應變εcr設置初始值,代入到式(26)和式(27)進行迭代求算,最終求出最小范數解,這里設置的最小范數的約束條件[29]為式(28)

2.3 裂縫尖端及裂縫面水平集更新

裂縫水平集更新的關鍵是將速度場從裂紋前緣擴展到整個水平集子域[41],裂縫尖端和裂縫面水平集的速度分量如圖5 所示,各個變量的關系式為

圖5 裂縫水平集函數速度分量示意圖Fig.5 Velocity component of crack level set function

裂縫水平集更新步驟主要如下:

(1) 對裂尖水平集進行初值化,使其滿足以下方程

(2) 定義裂縫的初始擴展速度,即式(32),本文θ取值45°,R取值0.005 m[42],再將1 維的速度場擴展到3 維的水平集子域,進行正交延展,即式(33)所示

根據Peng[43]的方法,采用Hamilton–Jacobi 方程對上述方程進行穩態求解,使其滿足方程式(34),計算表達式為

(3) 修正速度分量中裂紋面水平集的速度,修正目的是保證裂縫面內(ψ ＜0)的為0,使已生成裂縫不受水平集更新影響,保證裂紋前緣平滑[44],其中H(ψ)為階躍函數,修正方程為

(4) 采用修正后的速度,更新裂紋面水平集函數,再對裂紋面水平集函數進行重初值化,即

(5) 更新裂尖水平集函數,使其滿足以下條件

(6) 對水平集函數進行正交化,使其滿足?φ·?ψ=0,并進行再初值化處理,然而循環至第(2)步,即

3 應用算例

3.1 雙切口混凝土受剪拉作用開裂試驗

基于Nooru-Mohamed[45]的試驗,對比分析不同方法得出的混凝土受剪拉作用下的開裂情況.如圖6(a)所示,混凝土試件的長度和高度均為200 mm,厚度為50 mm,切口位于混凝土兩側中心,切口深度均為25 mm,寬度為5 mm,混凝土的強度為C30,混凝土底部和右下側部固定約束,左上側和頂部與鋼框架連接,加載鋼框架的長度為220 mm,寬度為120 mm,如圖6(b)所示.

圖6 雙切口混凝土受剪拉開裂試驗[45]Fig.6 Shear tensile cracking test of concrete contained double notched[45]

試驗時先在混凝土側部施加水平力Ph分別到5 kN,10 kN,27.5 kN,然后在頂端施加拉伸力Ps,試驗時在切口位置安裝了位移計測試了位移變化情況,即CTOD,型號為WPS 拉繩式位移計,測試精度為0.01 mm.

3.2 裂縫分布對比

圖7 為試驗及不同計算方法得出的裂縫分布對比圖,其中各種計算方法采用的最小網格尺寸為4 mm,可看出采用CDM 方法計算得出的裂縫形態在Ph=5 kN 時出現了分叉裂縫,與試驗差別較大,且裂縫后期的走向還受到網格依耐的影響;采用XFEM 方法計算得出的裂縫僅為一條,且均位于左端切口附近;采用原有的CDM-XFEM 方法計算得出的裂縫條數與試驗相同,但不同載荷階段的彎曲程度與試驗差距較大,如Ph=27.5 kN 時試驗從右切口延伸出的裂縫的拐點位于50 mm 的位置,但模擬得出裂縫的拐點卻位于150 mm 的位置.采用考慮塑性耗散的CDM-XFEM 方法計算得出的裂縫條數及走向與試驗均較為接近,雖然裂縫的彎曲程度與試驗有部分差異,但整體發展趨勢一致.以切口部位為中軸線,裂縫距中軸線的最大垂直距離為H,當Ph加載到27.5 kN 時,試驗、CDM 方法、XFEM 方法、既有的CDM-XFEM 方法、考慮塑性耗散的CDM-XFEM 方法算出的H值分別為46.5 mm,76.5 mm,69.3 mm,58.2 mm,48.2 mm,由此可知,采用考慮塑性耗散CDM-XFEM 計算出的混凝土裂縫形態與試驗差異最小.

圖7 試驗與各種方法算出的裂縫分布對比Fig.7 Comparison of crack distribution calculated by various methods and experiment

圖8 為考慮混凝土塑性耗散的CDM-XFEM 方法算出的不同載荷下的三維裂縫面云圖,對比圖7(a)可知,數值計算出的裂縫形態在混凝土前部和后部的裂縫走向相同,而試驗的前部和后部的裂縫走向略有差異,這是由于混凝土的開裂形態還受到內部骨料分布、骨料形態及ITZ 界面力學性能等因素的影響[47-49].

圖8 考慮塑性耗散的CDM-XFEM 計算出的三維裂縫分布云圖Fig.8 3D crack distribution calculated by CDM-XFEM considered plastic dissipation

3.3 拉力與切口位移關系曲線對比

圖9 為側向載荷Ph為5 kN 和10 kN 時,試驗與數值計算得出的混凝土切口處張開位移CTOD 與豎向拉力Ps的關系曲線,由圖可知,試驗曲線出現拐點的位移分別為15 kN 和12 kN,采用CDM 方法計算出的初始剛度大于試驗數據,而采用XFEM 方法計算出的初始剛度遠大于試驗數據,在后續的軟化階段,CDM 或XFEM 方法計算出的曲線均與試驗差異較大,原有的CDM-XFEM 方法計算出的曲線雖然與試驗走向基本一致,但在軟化階段與試驗的吻合程度不如考慮塑性耗散的CDM-XFEM 方法計算的曲線.在整個加載過程中,試驗與CDM 方法、XFEM 方法、原有的CDM-XFEM 方法、考慮塑性耗散的CDM-XFEM 方法計算得出的最大位移差值占比分別為26%,20%,14%,8%,說明采用考慮塑性耗散的CDM-XFEM 方法分析混凝土的開裂規律是合理可行的.

圖9 試驗與各種方法算出的拉力-張開位移曲線對比Fig.9 Comparison of the tension-open displacement curves obtained by experiment and various calculation methods

4 結論

本文通過考慮能量轉化時的塑性耗散,選取裂縫的牽引-分離函數模式,重新構建了能量轉化方程,最終提出了考慮混凝土塑性耗散的CDM-XFEM裂縫計算方法.通過與雙切口混凝土受剪拉開裂試驗進行對比分析,結果表明:

(1) 采用考慮塑性耗散的CDM-XFEM 計算方法能有效模擬混凝土的開裂過程,與其他計算方法相比,在裂縫形態方面與試驗差異最小;

(2) CDM 和XFEM 方法計算出的拉力-張開位移曲線均與試驗差異較大,原有的CDM-XFEM 方法計算出的曲線雖然與試驗走向基本一致,但在軟化階段與試驗的吻合程度不如考慮塑性耗散的CDM-XFEM 方法.說明采用考慮塑性耗散的CDMXFEM 方法分析混凝土的受力開裂是合理可行的.