何雨軒 卿光輝

（中國(guó)民航大學(xué)航空工程學(xué)院，天津 300300）

0 引言

裂紋損傷是引發(fā)災(zāi)難事故的主要原因之一，因此裂紋尖端附近的應(yīng)力、應(yīng)變和裂紋的擴(kuò)展規(guī)律一直是國(guó)內(nèi)外學(xué)者的研究熱點(diǎn)問(wèn)題。應(yīng)力強(qiáng)度因子是判斷已有裂紋是否擴(kuò)展的物理量，高精度的應(yīng)力強(qiáng)度因子計(jì)算結(jié)果是裂紋擴(kuò)展規(guī)律分析的重要保證。

在彈性范圍內(nèi)斷裂力學(xué)進(jìn)展中最重要的成就之一是1968年Rice提出的計(jì)算應(yīng)力強(qiáng)度因子的J積分方法［1-2］。文獻(xiàn)［3］采用位移外推法、J 積分法和虛擬閉合技術(shù)計(jì)算裂尖應(yīng)力強(qiáng)度因子，比較了三者的優(yōu)劣。主要結(jié)論是：在有限元網(wǎng)格相同的條件下，J 積分和虛擬閉合技術(shù)的強(qiáng)度因子結(jié)果精度高。另一方面，J 積分法最大的優(yōu)點(diǎn)在于解決了位移外推法和虛擬閉合技術(shù)高度依賴裂紋尖端區(qū)域網(wǎng)格細(xì)分的問(wèn)題，且具有積分路徑的無(wú)關(guān)性，是當(dāng)前求解應(yīng)力強(qiáng)度因子值最常用的計(jì)算方法之一。

基于位移有限元法模型研究復(fù)雜結(jié)構(gòu)中裂紋或結(jié)構(gòu)中復(fù)雜裂紋裂尖附近的應(yīng)力、應(yīng)變和裂紋的擴(kuò)展規(guī)律最為普遍［3-5］。一般情況下，基于有限元模型的J 積分算法的精度取決于位移結(jié)果的精度。事實(shí)上，基于最小勢(shì)能原理的位移元法用有限的結(jié)點(diǎn)代替了真實(shí)結(jié)構(gòu)中無(wú)窮多的點(diǎn)，減少了結(jié)構(gòu)的自由度數(shù)，因而這種模型較真實(shí)模型偏硬，給出的位移數(shù)值解是真解的下界，所以通常需要非常細(xì)的有限元網(wǎng)格模型才能得到比較理想的位移結(jié)果。

HERMANN［6］最早提出混合有限元法分析板殼問(wèn)題?；旌显Ｐ涂梢酝瑫r(shí)引入位移和應(yīng)力邊界條件［7-8］，較位移元模型增加了已知的應(yīng)力邊界條件約束，使得模型的剛度更符合實(shí)際情況，進(jìn)而可以提高位移變量和應(yīng)力變量的精度。最近，Qing 和Tian［9］結(jié)合最小勢(shì)能原理和H-R 變分原理建立了分析靜力學(xué)問(wèn)題的非協(xié)調(diào)的混合元。文獻(xiàn)［10］擴(kuò)展了這種混合單元的應(yīng)用范疇。非協(xié)調(diào)混合元模型的結(jié)果精度明顯優(yōu)于非協(xié)調(diào)位移元。

在商用軟件中，對(duì)平面靜態(tài)裂紋的應(yīng)力強(qiáng)度因子進(jìn)行數(shù)值計(jì)算時(shí)可提供的單元類型有Abaqus 中的非協(xié)調(diào)線性元CPS4I，協(xié)調(diào)非線性元CPS8［11］等，Ansys中的平面單元PLANE42，PLANE82 等。這些單元計(jì)算的應(yīng)力強(qiáng)度因子的方式也有多種，歸納成兩大類：基于節(jié)點(diǎn)位移求解（外推法、J 積分）和基于節(jié)點(diǎn)應(yīng)力求解（虛擬閉合），但由于其計(jì)算精度不高，不能很好地表現(xiàn)裂紋尖端的奇異性或需要額外重新設(shè)置成奇異單元［12］。本文針對(duì)應(yīng)力強(qiáng)度因子的精度問(wèn)題，構(gòu)建了在商用軟件中針對(duì)二維問(wèn)題還未曾應(yīng)用過(guò)的非協(xié)調(diào)混合元用于計(jì)算應(yīng)力強(qiáng)度因子。首先根據(jù)廣義H-R 變分原理建立含參數(shù)的非協(xié)調(diào)廣義混合有限元法的數(shù)學(xué)模型，然后簡(jiǎn)要地給出J積分算法的基本理論和過(guò)程，最后通過(guò)算例驗(yàn)證基于非協(xié)調(diào)的廣義混合有限元模型的J積分法的精確性和可靠性。

1 非協(xié)調(diào)四邊形廣義混合元法

假設(shè)線性彈性體的位移邊界條件事先滿足u-=0，則不考慮體積力的H-R 變分原理可表示為［13］

式中：σ表示應(yīng)力向量；C表示材料的剛度系數(shù)矩陣；?表示微分算子；u表示位移向量；V表示連續(xù)體的體積是作用在邊界表面上的載荷。

含參數(shù)α的廣義混合變分原理：

式中，參數(shù)0≤α≤1。

由式（1）和式（2）可得：

非協(xié)調(diào)的四邊形單元（4 節(jié)點(diǎn)）的應(yīng)力場(chǎng)和應(yīng)變場(chǎng)可近似表示為［13-14］：

式中，N和Nr分別為單元的協(xié)調(diào)部分和非協(xié)調(diào)部分（即單元內(nèi)部節(jié)點(diǎn)）的形函數(shù)矩陣。

將式（4）代入式（3）中，可得：

對(duì)式（5）中pe,qe,re分別進(jìn)行變分，消去re后可得非協(xié)調(diào)的廣義混合元的列式：

根據(jù)式（6）組裝各個(gè)單元后可得到模型的控制方程。通常情況下，控制方程中的參數(shù)α取0.75（參見文獻(xiàn)［15］）。從理論上講，非協(xié)調(diào)廣義混合元模型與混合有限元模型一樣，可以同時(shí)引入位移和應(yīng)力邊界條件。更重要的是模型中參數(shù)α的不同取值可以起到調(diào)節(jié)模型剛度的作用，因而可以提高位移結(jié)果精確度［8］。

2 J積分理論與應(yīng)力強(qiáng)度因子計(jì)算

2.1 J積分理論

如圖1 所示，考慮圍繞裂紋尖端的逆時(shí)針回路Γ，關(guān)于J積分的表達(dá)式如下：

圖1 圍繞裂紋尖端的逆時(shí)針?lè)e分回路Fig.1 Anti clockwise integral circuit around crack tip

式中：ui為位移矢量；ds為積分路徑Γ上的微小增量；w為應(yīng)變能密度因子，其定義為

式中，σij和εij分別為應(yīng)力張量和應(yīng)變張量。

Ti為垂直于回路的應(yīng)力矢量：

式中，nj為Γ的單位法矢量。

基于有限元模型的J 積分?jǐn)?shù)值算法參見文獻(xiàn)［3，16］。

2.2 應(yīng)力強(qiáng)度因子的計(jì)算

式中：J為J積分值；E為材料的楊氏模量；μ為泊松比。

3 數(shù)值算例

3.1 算例1:平板中心水平裂紋問(wèn)題[3]

如圖2所示，給定平面內(nèi)寬度為2b、長(zhǎng)度為2h的均質(zhì)平板，板厚1.0 mm，有一長(zhǎng)度為2a的中心裂縫，受到均勻?qū)ΨQ拉伸作用，平面應(yīng)力狀態(tài)下，材料的彈性模量E=200×103MPa，泊松比為μ=0.3，各向同性、均勻、線彈性，該裂紋板承受均勻應(yīng)力σ=30 MPa。

圖2 平板中心水平裂紋Fig.2 Plate with a central crack

幾何參數(shù)取a=20 mm，b=100 mm，h=200 mm。這是一個(gè)典型的Ι 型直裂紋問(wèn)題，其應(yīng)力強(qiáng)度因子KΙ的計(jì)算公式為［15］

因此，

由式（10）可得到解析解：

四分之一對(duì)稱的有限元模型，如圖3（a）所示。有限元網(wǎng)格劃分?jǐn)?shù)為100×200。

如圖3（b）所示，對(duì)稱邊界左側(cè)水平方向位移約束為零，垂直方向位移自由；底部非裂紋區(qū)域垂直方向位移約束為零，水平方向位移自由。應(yīng)力邊界條件如圖3（c），除了已知均布應(yīng)力σ，對(duì)模型右側(cè)給出水平方向應(yīng)力約束為零，垂直方向待求；底部裂紋區(qū)域豎直方向約束為零，水平方向待求。

圖3 初始邊界條件Fig.3 Initial boundary conditions

積分路徑與文獻(xiàn)［3］的相同，如圖4 所示。本文的計(jì)算結(jié)果為243.59，與解析解的結(jié)果243.6 MPa ?相比，誤差僅為-0.27%，該誤差遠(yuǎn)小于文獻(xiàn)［3］中的-1.3%。顯然，基于非協(xié)調(diào)的廣義混合元模型的精度明顯高于非協(xié)調(diào)位移元模型的精度。

圖4 1/4對(duì)稱模型中等效積分區(qū)域Fig.4 The equivalent integral region in 1/4 symmetric model

3.2 算例2:三點(diǎn)彎曲問(wèn)題[3]

三點(diǎn)彎曲模型是測(cè)試材料斷裂韌度的經(jīng)典模型。其幾何構(gòu)型和加載方式如圖5 所示，其中，L、W、a、S、B和P分別表示模型的長(zhǎng)度、寬度、裂長(zhǎng)、跨度、厚度和載荷。

圖5 三點(diǎn)彎曲模型Fig.5 Three point bending model

取L=220 mm，W=50 mm，S=200 mm，P=30 MPa，B=1 mm，彈性模量200×103MPa，泊松比為0.3。

應(yīng)力強(qiáng)度因子的表達(dá)為［17］

式中，形狀因子函數(shù)為

有限元模型的網(wǎng)格劃分?jǐn)?shù)為220×50，非協(xié)調(diào)位移元和非協(xié)調(diào)廣義混合元模型的結(jié)果與解析解的誤差列于表1。

表1 非協(xié)調(diào)位移元和非協(xié)調(diào)廣義混合元結(jié)果與解析解的誤差（SIF）Table 1 Errors between the results and analytical solutions of nonconforming displacement element and nonconforming generalized mixed element model（SIF）

從表1 中不難看出，非協(xié)調(diào)廣義混合元的精度明顯高于非協(xié)調(diào)位移元。

4 結(jié)論

（1）基于廣義H-R 變分原理建立了含參數(shù)的非協(xié)調(diào)廣義混合元的列式?；旌显Ｐ涂梢酝瑫r(shí)引入位移和應(yīng)力邊界條件，模型更加符合實(shí)際情況。另一方面，非協(xié)調(diào)廣義混合元模型中的參數(shù)α的不同取值可以調(diào)節(jié)模型的剛度，參數(shù)α合理的取值可以提高位移結(jié)果的精度。

（2）在非協(xié)調(diào)的廣義混合元模型的基礎(chǔ)上，采用J 積分法對(duì)平板中心水平裂紋和三點(diǎn)彎曲垂直裂紋兩個(gè)典型的實(shí)例進(jìn)行了一類應(yīng)力強(qiáng)度因子分析。數(shù)值結(jié)果表明，在有限元網(wǎng)格模型相同的情況下，非協(xié)調(diào)廣義混合元模型的應(yīng)力強(qiáng)度因子精度明顯高于非協(xié)調(diào)位移元模型。