半無限大彈性板裂紋尖端奇異應(yīng)力場(chǎng)計(jì)算

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

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

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

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

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

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

彈性力學(xué)平面問題的平衡方程：

式中：σ和u分別為應(yīng)力和位移，ρ為質(zhì)量密度，L為微分算子。比例邊界坐標(biāo)變換如圖1所示。

圖1 比例邊界坐標(biāo)變換

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

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

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

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

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

1.1 基本方程求解

靜力情況下，M0項(xiàng)為零。引入新變量：

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

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

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

首先求解Z的特征值問題

式中Λ為對(duì)角矩陣。然后將Λ與Φ進(jìn)行分塊排列

可得：

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

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

對(duì)于有限域來說，有限域剛度為

對(duì)于無限域來說，無限域剛度為

1.2 裂紋尖端奇異應(yīng)力場(chǎng)

如圖2所述，有裂紋的子塊內(nèi)保持比例邊界有限元的特點(diǎn)。

圖2 帶裂紋子結(jié)構(gòu)塊（超單元）

式(1)中位移場(chǎng)的解為

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

應(yīng)力σ(ξ)為

式中：D為彈性矩陣，B1(η) 、B2(η)計(jì)算過程見文獻(xiàn)[7]和文獻(xiàn)[9] 。將位移式(5)代入到式(6)得到子結(jié)構(gòu)內(nèi)部點(diǎn)的應(yīng)力：

式(5)和式(6)為子塊內(nèi)的位移和應(yīng)力計(jì)算公式。為了計(jì)算裂紋尖端奇異應(yīng)力場(chǎng)，將子塊內(nèi)的位移和應(yīng)力的表達(dá)式表示為極坐標(biāo)的形式比較方便。此時(shí)，徑向坐標(biāo)為

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

將式(8) 和式(9)組合可形成Williams 以極坐標(biāo)表達(dá)類似的應(yīng)力場(chǎng)。Williams[1]裂紋尖端奇異應(yīng)力場(chǎng)對(duì)應(yīng)關(guān)系參見文獻(xiàn)[5] 。

2 數(shù)值算例

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

圖3為半無限大彈性板受均布拉伸荷載的計(jì)算模型。其中有限域第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塊的相似中心坐標(biāo)是(0,0)，無限域（第2塊）的相似中心坐標(biāo)是(-1,0)，第3塊的相似中心坐標(biāo)是(-1,2)。本文結(jié)果與解析解提取的高階奇異參數(shù)a2、a3見表1。

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

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

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

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

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

圖5為半無限大彈性板受集中剪切荷載的計(jì)算模型。其中寬度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，中心坐標(biāo)取在裂紋尖端（0,0）。

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

本文結(jié)果與文獻(xiàn)[15]給出的解析解對(duì)比見表2。

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

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

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

圖6為半無限大彈性板受均布剪切荷載的計(jì)算模型。計(jì)算中寬度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。中心取在裂紋尖端坐標(biāo)（0，0）。文獻(xiàn)[15] 給出了解析解，本文結(jié)果與解析解對(duì)比見表3。

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

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

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

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

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

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

表4 半無限大彈性板裂紋邊受集中拉伸荷載計(jì)算結(jié)果

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

3 結(jié)論

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