基于非協調廣義混合元的應力強度因子計算

何雨軒 卿光輝

（中國民航大學航空工程學院，天津 300300）

0 引言

裂紋損傷是引發災難事故的主要原因之一，因此裂紋尖端附近的應力、應變和裂紋的擴展規律一直是國內外學者的研究熱點問題。應力強度因子是判斷已有裂紋是否擴展的物理量，高精度的應力強度因子計算結果是裂紋擴展規律分析的重要保證。

在彈性范圍內斷裂力學進展中最重要的成就之一是1968年Rice提出的計算應力強度因子的J積分方法［1-2］。文獻［3］采用位移外推法、J 積分法和虛擬閉合技術計算裂尖應力強度因子，比較了三者的優劣。主要結論是：在有限元網格相同的條件下，J 積分和虛擬閉合技術的強度因子結果精度高。另一方面，J 積分法最大的優點在于解決了位移外推法和虛擬閉合技術高度依賴裂紋尖端區域網格細分的問題，且具有積分路徑的無關性，是當前求解應力強度因子值最常用的計算方法之一。

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

HERMANN［6］最早提出混合有限元法分析板殼問題。混合元模型可以同時引入位移和應力邊界條件［7-8］，較位移元模型增加了已知的應力邊界條件約束，使得模型的剛度更符合實際情況，進而可以提高位移變量和應力變量的精度。最近，Qing 和Tian［9］結合最小勢能原理和H-R 變分原理建立了分析靜力學問題的非協調的混合元。文獻［10］擴展了這種混合單元的應用范疇。非協調混合元模型的結果精度明顯優于非協調位移元。

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

1 非協調四邊形廣義混合元法

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

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

含參數α的廣義混合變分原理：

式中，參數0≤α≤1。

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

非協調的四邊形單元（4 節點）的應力場和應變場可近似表示為［13-14］：

式中，N和Nr分別為單元的協調部分和非協調部分（即單元內部節點）的形函數矩陣。

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

對式（5）中pe,qe,re分別進行變分，消去re后可得非協調的廣義混合元的列式：

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

2 J積分理論與應力強度因子計算

2.1 J積分理論

如圖1 所示，考慮圍繞裂紋尖端的逆時針回路Γ，關于J積分的表達式如下：

圖1 圍繞裂紋尖端的逆時針積分回路Fig.1 Anti clockwise integral circuit around crack tip

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

式中，σij和εij分別為應力張量和應變張量。

Ti為垂直于回路的應力矢量：

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

基于有限元模型的J 積分數值算法參見文獻［3，16］。

2.2 應力強度因子的計算

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

3 數值算例

3.1 算例1:平板中心水平裂紋問題[3]

如圖2所示，給定平面內寬度為2b、長度為2h的均質平板，板厚1.0 mm，有一長度為2a的中心裂縫，受到均勻對稱拉伸作用，平面應力狀態下，材料的彈性模量E=200×103MPa，泊松比為μ=0.3，各向同性、均勻、線彈性，該裂紋板承受均勻應力σ=30 MPa。

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

幾何參數取a=20 mm，b=100 mm，h=200 mm。這是一個典型的Ι 型直裂紋問題，其應力強度因子KΙ的計算公式為［15］

因此，

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

四分之一對稱的有限元模型，如圖3（a）所示。有限元網格劃分數為100×200。

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

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

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

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

3.2 算例2:三點彎曲問題[3]

三點彎曲模型是測試材料斷裂韌度的經典模型。其幾何構型和加載方式如圖5 所示，其中，L、W、a、S、B和P分別表示模型的長度、寬度、裂長、跨度、厚度和載荷。

圖5 三點彎曲模型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。

應力強度因子的表達為［17］

式中，形狀因子函數為

有限元模型的網格劃分數為220×50，非協調位移元和非協調廣義混合元模型的結果與解析解的誤差列于表1。

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

從表1 中不難看出，非協調廣義混合元的精度明顯高于非協調位移元。

4 結論

（1）基于廣義H-R 變分原理建立了含參數的非協調廣義混合元的列式。混合元模型可以同時引入位移和應力邊界條件，模型更加符合實際情況。另一方面，非協調廣義混合元模型中的參數α的不同取值可以調節模型的剛度，參數α合理的取值可以提高位移結果的精度。

（2）在非協調的廣義混合元模型的基礎上，采用J 積分法對平板中心水平裂紋和三點彎曲垂直裂紋兩個典型的實例進行了一類應力強度因子分析。數值結果表明，在有限元網格模型相同的情況下，非協調廣義混合元模型的應力強度因子精度明顯高于非協調位移元模型。