半圓形初始裂紋在交變拉伸載荷下的裂尖應(yīng)力場數(shù)值模擬

周曉松，張焱冰，張岳林

1中國人民解放軍軍事科學(xué)院 國防科技創(chuàng)新研究院，北京100071

2海軍工程大學(xué) 艦船與海洋學(xué)院，湖北武漢430033

3中國人民解放軍91404部隊，河北秦皇島 066001

0 引 言

裂紋普遍存在于航空航天及海洋結(jié)構(gòu)物中，例如，因材料缺陷和加工引起的初始裂紋以及疲勞載荷引起的萌生裂紋。因此，開展裂紋擴展規(guī)律及裂尖應(yīng)力場的研究，對于準(zhǔn)確評估板結(jié)構(gòu)的剩余承載能力具有重要意義，尤其是對于承受交變波浪載荷的船體結(jié)構(gòu)。為此，國內(nèi)外學(xué)者開展了此類結(jié)構(gòu)的極限強度研究。例如，Cui等［1］針對一系列不同長度、位置、方向角的裂紋，對鋼板極限強度影響的非線性進行了有限元分析；Bayatfar等［2］對在單軸壓縮載荷作用下的非加筋板和加筋板單元的極限強度特性進行了研究，發(fā)現(xiàn)裂紋長度及其位置對極限強度有著重要影響；閆小順等［3］提出了一種坐標(biāo)變換建模方法，用來解決球柱結(jié)合殼焊趾表面裂紋的有限元建模問題；吳建國等［4］基于能量差率基本原理，以非穿透裂紋深度的相對值作為控制裂紋虛比例擴展的無量綱幾何參量，來構(gòu)造求解裂紋張開位移幅值的伯努利微分方程，從而導(dǎo)出以裂紋絕對尺寸和相對尺寸為參數(shù)的裂紋張開位移幅值的表達(dá)式，得到了有限大體非穿透裂紋三維應(yīng)力強度因子的閉合解。在上述研究中，大多數(shù)是將含裂紋的非加筋板或加筋板簡化為二維薄殼結(jié)構(gòu)，即認(rèn)為裂紋厚度沿板截面是均勻分布的。然而，在實際工程應(yīng)用中，裂紋在板截面上的形狀往往是不規(guī)則的，故應(yīng)用實體單元對裂紋擴展和裂尖應(yīng)力場進行分析有望提高計算的精度。

本文將以半圓形三維初始裂紋為例，基于虛擬裂紋閉合技術(shù)（Virtual Crack Closure Technique，VCCT），運用Marc有限元程序，對裂紋在交變拉伸載荷作用下的裂尖應(yīng)力場進行分析，從而為考慮初始裂紋的板結(jié)構(gòu)疲勞強度計算提供參考。

1 Marc中的虛擬裂紋閉合技術(shù)

虛擬裂紋閉合技術(shù)最先由Rybicki和Kanninen［5］提出，主要用于計算二維裂紋的問題。該方法通過對有限元分析結(jié)果進行后處理，得到所需裂紋擴展的能量釋放率。Shivakumar等［6］后來將該方法推廣到了三維裂紋的計算中，Narayana等［7］則進一步發(fā)展了針對面狀裂紋的一般形式的計算公式。在上述方法中，由于所需的斷裂參數(shù)是對有限元分析結(jié)果進行后處理得到的，因此比較繁雜，也削弱了有限元分析和參數(shù)計算的有機聯(lián)系。為此，Xie和 Biggers［8］提出了啞節(jié)點斷裂單元（Fracture element with dummy nodes）方法。圖1所示為一個九節(jié)點塊體單元的虛擬裂紋閉合法示意圖［4］。圖中：只有位于裂紋尖端的節(jié)點1和2對單元的剛度矩陣產(chǎn)生影響；節(jié)點3和4用于計算裂紋張開位移；節(jié)點5～9用于確定局部坐標(biāo)系并計算虛擬裂紋的擴展面積。為計算裂紋尖端處的節(jié)點力，在節(jié)點1和2之間放置剛度很大的強彈簧。

圖1 九節(jié)點塊體單元的虛擬裂紋閉合法示意圖［4］Fig.1 Sketch of 9-node solid elements using VCCT［4］

由此，可得到3種斷裂形式的應(yīng)變能釋放率，分別表示如下：

式中：GⅠ，GⅡ和GⅢ分別為3種斷裂形式的應(yīng)變能釋放率；為由全局坐標(biāo)系轉(zhuǎn)換為局部坐標(biāo)后節(jié)點 1上的節(jié)點力；Δu3，4，Δv3，4和Δw3，4為局部坐標(biāo)下節(jié)點3和4之間的相對位移；A為裂紋深度。

對于復(fù)合型斷裂，當(dāng)斷裂判據(jù)滿足式（2）時，

此時，節(jié)點1和節(jié)點2處的強彈簧完全失效，裂紋擴展。其中，GⅠC，GⅡC和GⅢC分別為3種斷裂形式應(yīng)變能釋放率的門檻值。

2 裂紋擴展的數(shù)值模擬

2.1 模型概述

圖2所示為本文構(gòu)建的數(shù)值模型。具體尺寸及相關(guān)參數(shù)如下：板尺寸為200 mm×100 mm×20 mm，單元數(shù)劃分為26 660個，節(jié)點為4 743個，裂紋半徑為5 mm。材料為DQSK36，材料密度為7.85×10-9t/mm3，楊氏模量為2.07×105MPa，泊松比為0.28。圖3所示為DQSK36材料的塑性力學(xué)特性曲線。本文模型采取左側(cè)固支，右側(cè)施加隨時間t變化的載荷P，并由下式計算：

圖2 模型網(wǎng)格劃分Fig.2 Meshing of the model

圖3 DQSK36材料塑性力學(xué)特性Fig.3 Plasticity mechanical characteristics of material DQSK36

2.2 初始裂紋定義

本文將模型視為接觸體，并基于VCCT方法指定在疲勞加載模式下裂紋擴展的相關(guān)參數(shù)。其中，疲勞載荷周期取為1 s，裂紋擴展增量取為4 mm，裂紋擴展方向采用最大Hoopstress方法。

在構(gòu)建初始裂紋尖端和裂紋擴展尖端周圍的網(wǎng)格時，采用全局網(wǎng)格重劃分技術(shù)自動完成，并使用了Patran四面體網(wǎng)格生成器。通過距離控制單元尺寸的變化，在距離裂紋尖端20 mm以內(nèi)的區(qū)域，對網(wǎng)格進行細(xì)化，接近裂紋尖端部位采用單元長度為3 mm的網(wǎng)格重新劃分，逐漸遠(yuǎn)離裂紋尖端時采用過渡網(wǎng)格，其他部位采用均勻的10 mm單元尺寸劃分。

采用Crack initators定義初始裂紋尖端自動建模的參數(shù)，通過小面構(gòu)成的曲面在接觸體中插入模板裂紋。然后，分別定義預(yù)載荷和往復(fù)加載分析工況，其中預(yù)載荷工況采用1個增量步，在0.5 s內(nèi)完成；往復(fù)加載工況采用10個增量步，在5 s內(nèi)完成，激活全局網(wǎng)格重劃分和初始裂紋建模選項，提交計算。

3 裂紋水平面的應(yīng)力場

打開Marc有限元程序的結(jié)果后處理文件，讀取計算結(jié)果，得到如圖4所示的裂紋水平面的應(yīng)力云圖。由圖可以看出，由于模型的對稱性，裂尖應(yīng)力場關(guān)于裂紋中心點對稱分布。

圖4 裂紋水平面的應(yīng)力云圖Fig.4 Stress contours of crack in horizontal plane

圖5所示為以下面的裂尖為原點，在極坐標(biāo)系下，在不同時刻（t=1.5，2 s）時的von Mises等效應(yīng)力隨r的變化曲線。圖中，r為計算點離開裂紋尖端的距離，a為裂紋長度，θ為裂紋與載荷方向的夾角。

當(dāng)t=1.5 s時，在θ=0或 π位置，計算點的von Mises等效應(yīng)力隨著離開裂尖距離的增大而呈線性增大，當(dāng)r=0.08a時達(dá)到最大值，為155 MPa，此時的載荷為100 MPa。用該值無量綱化計算的點應(yīng)力為1.55P，此后計算點的von Mises等效應(yīng)力隨離開裂尖距離的增大而呈線性減小。

當(dāng)t=2 s時，載荷由100 MPa降低到50 MPa，降幅達(dá) 50%。在θ=0，π，π/6，5π/6，-π/6，-5π/6，-π/3，-2π/3，-π/2 位置，隨著載荷的減小，計算點的von Mises等效應(yīng)力比t=1.5 s時刻的有所減小，但幅度小于載荷減小的幅度（50%）。這說明計算點發(fā)生了塑性變形，計算點的von Mises等效應(yīng)力隨r的變化規(guī)律與t=1.5 s時刻的相似，即基本上呈線性或雙線性關(guān)系。而在θ=π/3，2π/3，π/2位置，在載荷變小的情況下，當(dāng)0≤r＜0.35a時，von Mises等效應(yīng)力隨r的增大而迅速減小，當(dāng)r=0.35a時達(dá)到最小。這是因為在整個加載過程中r=0.35a附近區(qū)域的von Mises等效應(yīng)力水平均未超過材料的屈服應(yīng)力，其未發(fā)生塑性應(yīng)變。

圖5 不同時刻von Mises等效應(yīng)力隨r的變化Fig.5 Variation of von Mises equivalent stress with r at different times

4 裂紋橫斷面的應(yīng)力場

圖6為在不同時刻裂紋橫斷面的應(yīng)力云圖。由圖可以看出，在裂紋橫斷面上，裂紋的擴展方向總是沿著裂紋形狀的法向，即以裂紋中心為原點，在極坐標(biāo)下，a為任意值時計算點的von Mises等效應(yīng)力隨r的變化均一致。由此，本文繪制了任意a位置時的裂紋橫斷面應(yīng)力隨r變化的曲線，分別如圖7和圖8所示。

圖6 裂紋橫斷面不同時刻的應(yīng)力云圖Fig.6 Stress contours of crack cross-section at different times

圖7 不同時刻裂紋橫斷面應(yīng)力隨r的變化曲線Fig.7 Variation of crack stress in cross section with r at different times

圖8 不同時刻裂紋橫斷面應(yīng)力隨r的變化曲線Fig.8 Variation of stress at crack cross section with r at different times

當(dāng)t=1.5 s時，計算點的von Mises等效應(yīng)力隨r的增大而增大，在r=0.5a時達(dá)到最大，此位置恰好在裂紋邊緣，其應(yīng)力值為1.45P。此后，von Mises等效應(yīng)力隨r的增大而減小，當(dāng)r=0.8a時，von Mises等效應(yīng)力-r/a曲線逐漸呈平穩(wěn)的直線，應(yīng)力值趨向于P，即距原點的距離大于0.8a的位置不再產(chǎn)生應(yīng)力集中。

當(dāng)t=2 s時，計算點的von Mises等效應(yīng)力隨r的增大先增后減，在r=0.9a時達(dá)到峰值2P。此后，von Mises等效應(yīng)力隨r的增大而減小，在r=1.4a時，應(yīng)力值趨向于P。

對于t=2.5s時的情況，當(dāng)r＜0.9a時，計算點的von Mises等效應(yīng)力隨r的變化仍然呈現(xiàn)出高度非線性，在r=0.9a時，von Mises等效應(yīng)力達(dá)到峰值2P，此后，von Mises等效應(yīng)力隨r的增大而減??；在r=1.5a時，von Mises等效應(yīng)力值趨向于P。

綜上所述并結(jié)合圖8，得出了如表1所示的不同時刻von Mises等效應(yīng)力出現(xiàn)的位置（瞬時裂紋半徑rt）、應(yīng)力集中開始消失的位置rd以及裂紋邊緣應(yīng)力峰值σmax。其中，比值rd/rt表示裂紋位置的嚴(yán)重程度，值越大，應(yīng)力集中越嚴(yán)重。

表1 裂紋橫斷面應(yīng)力場參數(shù)隨時間變化Table 1 Parameters of the stress field in cross section of crack varying with different times

由表1可知：在任意時刻下，應(yīng)力集中總是從1.38～1.67倍裂紋半徑的位置開始消失；當(dāng)載荷為2/3倍材料屈服應(yīng)力時，裂紋邊緣應(yīng)力峰值總是等于1.45P；當(dāng)載荷為1/3倍材料屈服應(yīng)力時，裂紋邊緣應(yīng)力峰值總是等于2P。

需要指出的是，由于在t=4 s后裂紋擴展到板厚近末端，故整個模型中不再出現(xiàn)應(yīng)力集中消失的位置，尤其是在t=5 s時，結(jié)合圖6中t=5 s時的裂紋橫斷面不同時刻的應(yīng)力云圖，裂紋擴展到板厚邊緣貫穿了整個板厚。

5 結(jié) 論

本文基于三維虛擬裂紋閉合技術(shù)，以半圓形初始裂紋為例，研究了其在交變載荷作用下的裂尖應(yīng)力場。根據(jù)本文的研究，得到以下結(jié)論：

1）對于裂紋水平面的應(yīng)力場，在發(fā)生塑性變形的區(qū)域，其von Mises等效應(yīng)力隨計算點距裂尖距離的增大而減小，兩者之間基本上呈線性或雙線性關(guān)系；在未發(fā)生塑性變形的區(qū)域，基本上是裂尖距離的三次函數(shù)。

2）對于裂紋橫斷面，裂紋總是沿裂紋邊緣形狀的法向方向擴展，應(yīng)力最大值總是出現(xiàn)在裂紋邊緣。在裂紋內(nèi)部，應(yīng)力值隨著距裂紋中心的距離而呈現(xiàn)高度的非線性分布；在裂紋外部，應(yīng)力值隨著距裂紋中心距離的增大而減小。當(dāng)距裂紋中心的距離大于裂紋半徑1.67倍時，應(yīng)力集中消失。

3）同樣，對于裂紋橫斷面，當(dāng)載荷為材料屈服應(yīng)力的2/3倍時，裂紋邊緣最大應(yīng)力值總是為載荷的1.45倍；當(dāng)載荷為材料屈服應(yīng)力的1/3倍時，裂紋邊緣最大應(yīng)力值總是為載荷的2倍。