劉鈞玉，張?zhí)煊?，蘇艷，寧寶寬

沈陽工業(yè)大學 建筑與土木工程學院，遼寧 沈陽 110870

在Williams[1]給出的裂紋尖端奇異應力場解析表達式中，奇異參數(shù)包括張開型（I型裂紋）應力強度因子（stress intensity factor）KI和剪切型（II型裂紋）應力強度因子KII，除此之外還包括被稱為裂尖參數(shù)的T應力（T-stress）和高階奇異參數(shù)。隨著時代進步，計算機數(shù)值計算手段日趨發(fā)展，針對T應力和高階項等奇異參數(shù)的相關研究內容逐漸增多[2-6]。裂紋尖端奇異場的參數(shù)受很多因素影響，包括結構的幾何構型、材料屬性、荷載類別以及裂紋傾角等，理論分析方法可以求解簡單情況下的參數(shù)，但對于不同因素條件下的實際問題，一般使用數(shù)值方法來對參數(shù)進行計算。目前，計算裂紋尖端奇異應力場的數(shù)值方法主要有權函數(shù)法、有限差分法、有限元法、邊界元法以及比例邊界有限元法等。

有限元法需要對裂紋尖端部位進行離散，邊界元法需要對裂紋面進行離散，且二者都是基于無法近似奇異點周圍解析解的分片光滑函數(shù)[7]。比例邊界有限元法具有有限元法和邊界元法的優(yōu)勢。在計算裂紋尖端奇異應力場時，在裂紋面通過相似中心的情況下，可以不用離散裂紋面，這為前處理提供了便捷，使計算量減少。除此之外，比例邊界有限元法相比于有限元法不需要基本解，且計算的數(shù)值結果在徑向是精確的，同樣也是環(huán)向收斂于有限元意義的解，裂紋尖端奇異場的這些參數(shù)由此可根據(jù)定義直接進行提取。因此比例邊界有限元法在當下已在很多方面得到應用，例如對各向異性材料[7]、溫度應力作用下的奇異應力場[8]的計算，以及對動應力強度因子[9]及不同材料高階奇異性的求解[5]。文獻[10]計算了劈拉試件裂紋尖端的奇異應力場，文獻[11]對比例邊界有限元法近年來的發(fā)展以及T應力和裂紋尖端奇異參數(shù)的計算進行了總結。

本文基于比例邊界有限元法對半無限大彈性板裂紋尖端的奇異應力場進行計算，提取了裂紋尖端處的T應力以及高階奇異項等參數(shù)。與數(shù)值求解結果進行對比驗證了其精度與有效性。且對半無限大彈性板進行了斷裂分析，通過改變裂紋長度和加載方式等影響因素，得出了裂紋尖端T應力以及高階奇異參數(shù)的變化規(guī)律。在對比后可以看出T應力以及高階奇異參數(shù)對于判斷裂紋的斷裂特性具有一定價值。

1 比例邊界有限元法基本原理

彈性力學平面問題的平衡方程：

式中：σ和u分別為應力和位移，ρ為質量密度，L為微分算子。比例邊界坐標變換如圖1所示。

圖1 比例邊界坐標變換

式(2)為比例邊界坐標變換公式：

式中：N(η)為插值形函數(shù)，為相似中心橫縱坐標，ξ、η為比例邊界有限元坐標。將邊界ξ=1上結點的坐標(x?,y?)用(x,y)表示。對式(1)按式(2)將坐標轉換為比例邊界坐標ξ、η，可得：

該方程組為二階線性常微分方程組，沿 ξ方向的位移場{u(ξ)}可以解析求解。

系數(shù)矩陣E0、E1、E2和M0可以用數(shù)值積分方法直接計算，但是僅需在邊界上進行離散計算[12-14]。

式中：M0項在動力荷載問題中存在，Nu(η)為位移插值函數(shù)。

1.1 基本方程求解

靜力情況下，M0項為零。引入新變量：

式中Q(ξ)代表域內一點的等效節(jié)點力。

式(3)轉化為一階常微分方程：

式中Z為哈密頓系數(shù)矩陣，具體形式為

首先求解Z的特征值問題

式中Λ為對角矩陣。然后將Λ與Φ進行分塊排列

可得：

式中λi、-λi均為Λ的一組特征值（λi的實部為負）。由此式(4)的解可設定為

式中c1、c2為積分常數(shù)。由邊界處位移u(ξ=1)求解。

對于有限域來說，有限域剛度為

對于無限域來說，無限域剛度為

1.2 裂紋尖端奇異應力場

如圖2所述，有裂紋的子塊內保持比例邊界有限元的特點。

圖2 帶裂紋子結構塊（超單元）

式(1)中位移場的解為

式中：n為方陣Φ11的維數(shù)，Φ為特征向量矩陣，?i是矩陣Φ11的第i列，ci是積分常數(shù)矢量c1的第i個元素。插值函數(shù)Nu(η)可以計算單元位移場。

應力σ(ξ)為

式中：D為彈性矩陣，B1(η) 、B2(η)計算過程見文獻[7]和文獻[9] 。將位移式(5)代入到式(6)得到子結構內部點的應力：

式(5)和式(6)為子塊內的位移和應力計算公式。為了計算裂紋尖端奇異應力場，將子塊內的位移和應力的表達式表示為極坐標的形式比較方便。此時，徑向坐標為

將式(8)代入到式(7)得到

將式(8) 和式(9)組合可形成Williams 以極坐標表達類似的應力場。Williams[1]裂紋尖端奇異應力場對應關系參見文獻[5] 。

2 數(shù)值算例

2.1 半無限大板裂紋邊受均布拉伸荷載

圖3為半無限大彈性板受均布拉伸荷載的計算模型。其中有限域第1塊（相似中心O1）的寬度取W=1.5，高度H=4，裂紋長度a=1，均布力加載長度b分別取0.35、0.5、0.75和0.8。彈性模量E=1，泊松比ν=1.5。圖4給出了比例邊界有限元法離散網(wǎng)格的示意圖，豎邊共離散單元數(shù)N=16，第1塊的相似中心坐標是(0,0)，無限域（第2塊）的相似中心坐標是(-1,0)，第3塊的相似中心坐標是(-1,2)。本文結果與解析解提取的高階奇異參數(shù)a2、a3見表1。

圖3 半無限大彈性板受均布拉伸荷載

圖4 半無限大彈性板網(wǎng)格離散示意

表1 半無限大彈性板裂紋邊受均布拉伸荷載（W=1.5，N=16）

由表1中數(shù)據(jù)可以看到比例邊界有限元法的計算結果與解析解的最大誤差為3%，這說明比例邊界有限元法能夠較為精確地計算應力強度因子。

2.2 半無限大彈性板裂紋邊受集中剪切荷載

圖5為半無限大彈性板受集中剪切荷載的計算模型。其中寬度W=2 ，高度H=4 ，裂紋長度a分別取0.7、0.8、0.9、1.0、1.1、1.2和1.3。彈性模量E=1，泊松比ν=0.25。共離散單元數(shù)N=19 。比例邊界有限元法離散網(wǎng)格的示意圖見圖4，中心坐標取在裂紋尖端（0,0）。

圖5 半無限大彈性板受集中剪切荷載

本文結果與文獻[15]給出的解析解對比見表2。

表2 半無限大彈性板裂紋邊受集中剪切荷載

由表2中數(shù)據(jù)可以看到本文計算結果與解析解的最大誤差在2%以內，表明比例邊界有限元法能夠較為精確地計算應力強度因子。

2.3 半無限大彈性板裂紋邊受均布剪切荷載

圖6為半無限大彈性板受均布剪切荷載的計算模型。計算中寬度W=1.5，高度H=4，裂紋長度取a=1，均布力加載長度b分別取 0.65、 0.70、0.75、0.80、0.85、0.90和0.95 。彈性模量E=1，泊松比ν=0.25。共離散單元數(shù)取N=8。中心取在裂紋尖端坐標（0，0）。文獻[15] 給出了解析解，本文結果與解析解對比見表3。

圖6 半無限大彈性板受均布剪切荷載

表3 半無限大彈性板裂紋邊受均布剪切荷載

由表3中數(shù)據(jù)可以看到本文計算結果與解析解的最大誤差在4%以內，表明比例邊界有限元法可以較為精確地計算應力強度因子。

2.4 半無限大彈性板裂紋邊受集中拉伸荷載

圖7為半無限大彈性板受集中拉伸荷載的計算模型。計算中寬度W=2，高度H=4，裂紋長度a=1，彈性模量E=1，泊松比ν=0.25。模型受集中力P=1。中心取在裂紋尖端坐標為（0，0）。文獻[15]給出了解析解，本文結果與解析解和數(shù)值結果的對比見表4。

圖7 半無限大彈性板裂紋邊受集中拉伸荷載

表4 半無限大彈性板裂紋邊受集中拉伸荷載計算結果

計算中離散單元數(shù)取N=11。由表4中數(shù)據(jù)可以看到本文結果a1、a2與解析解最大誤差在2%以內，表明了比例邊界有限元法能夠較為精確地計算應力強度因子；而高階奇異參數(shù)a3的計算結果與文獻中結果的最大誤差在4%以內，表明本文方法能夠較為精確地提取高階奇異參數(shù)。

3 結論

本文基于比例邊界有限元法計算了半無限大彈性板裂紋尖端奇異應力場，給出了靜力無限域剛度的計算方程，分析了比例邊界元法位移場、應力場表達式與斷裂力學裂紋尖端奇異場參數(shù)的對應關系。對半無限大彈性板在集中荷載和均布荷載作用下的裂紋尖端奇異場進行了研究，提取了應力強度因子、T 應力以及高階項等奇異參數(shù)，將提取結果與解析解和文獻中結果進行了對比，結果表明比例邊界有限元法能夠較精確地計算半無限大彈性板裂紋尖端應力場的奇異參數(shù)。對于拉伸荷載作用下，隨著裂紋長度的增加，應力強度因子逐漸增大；對于集中剪切荷載作用下，隨著裂紋長度的增加，應力強度因子逐漸降低；均布剪切荷載作用下應力強度因子逐漸增加。