周曉松，張焱冰，張岳林

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

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

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

0 引 言

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

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

1 Marc中的虛擬裂紋閉合技術

虛擬裂紋閉合技術最先由Rybicki和Kanninen［5］提出，主要用于計算二維裂紋的問題。該方法通過對有限元分析結果進行后處理，得到所需裂紋擴展的能量釋放率。Shivakumar等［6］后來將該方法推廣到了三維裂紋的計算中，Narayana等［7］則進一步發(fā)展了針對面狀裂紋的一般形式的計算公式。在上述方法中，由于所需的斷裂參數(shù)是對有限元分析結果進行后處理得到的，因此比較繁雜，也削弱了有限元分析和參數(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用于確定局部坐標系并計算虛擬裂紋的擴展面積。為計算裂紋尖端處的節(jié)點力，在節(jié)點1和2之間放置剛度很大的強彈簧。

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

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

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

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

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

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

2.1 模型概述

圖2所示為本文構建的數(shù)值模型。具體尺寸及相關參數(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材料的塑性力學特性曲線。本文模型采取左側固支，右側施加隨時間t變化的載荷P，并由下式計算：

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

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

2.2 初始裂紋定義

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

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

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

3 裂紋水平面的應力場

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

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

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

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

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

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

4 裂紋橫斷面的應力場

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

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

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

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

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

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

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

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

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

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

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

5 結 論

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

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

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

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