高文，王生楠

（西北工業(yè)大學(xué)航空學(xué)院，西安 710072）

0 引言

工程結(jié)構(gòu)由于材料夾雜或服役等原因不可避免地存在缺陷或裂紋，這些缺陷或裂紋在疲勞載荷的作用下逐漸擴(kuò)展最終導(dǎo)致結(jié)構(gòu)發(fā)生斷裂。據(jù)統(tǒng)計(jì)，大約60%的結(jié)構(gòu)失效是由疲勞裂紋擴(kuò)展引起的，因此，需要對(duì)含裂紋結(jié)構(gòu)的安全性進(jìn)行量化分析。在線性彈性斷裂力學(xué)（Linear Elastic Fracture Mechanics，簡稱LEFM）范圍內(nèi)，應(yīng)力強(qiáng)度因子（Stress Intensity Factors，簡稱SIFs）用來表征裂紋尖端附近應(yīng)力場強(qiáng)度，是研究裂紋起裂和擴(kuò)展的關(guān)鍵因素。但只有少數(shù)形式簡單的裂紋可以通過試驗(yàn)或解析的方式獲得應(yīng)力強(qiáng)度因子解，對(duì)于形式復(fù)雜的裂紋問題，不可避免地要使用有限差分法、有限元法、邊界元法等數(shù)值方法來求解。在這些方法中，有限元法由于其較高的求解效率和可以模擬任意復(fù)雜幾何形狀而在工程界占據(jù)主導(dǎo)地位。

本文以二次四面體單元VCCM法為基礎(chǔ)，對(duì)應(yīng)用四面體單元進(jìn)行三維裂紋有限元建模方法進(jìn)行研究，對(duì)VCCM法計(jì)算斷裂參數(shù)的影響因素進(jìn)行分析，裂紋建模過程需要繁瑣的手工操作，全自動(dòng)化建模也是本文的研究重點(diǎn)。

1 虛擬裂紋閉合法

1.1 二次四面體單元的VCCM法計(jì)算式

由VCCM法的理論可知，在裂紋的虛擬擴(kuò)展面上需要有單元面片存在以獲得計(jì)算所需的節(jié)點(diǎn)力，且裂紋前緣兩側(cè)的面片寬度需相等。四面體單元的面是三角形，典型的裂紋前緣二次四面體單元分布示意圖如圖1所示。裂紋前緣的三角形單元面片分為兩種：一種是三角面片的一條邊位于裂紋前緣，記為；另一種是三角面片的一個(gè)頂點(diǎn)位于裂紋前緣，記為，Δ為單元寬度，如圖2所示。面片和上的計(jì)算式如式（1）～式（2）所示。

圖1 典型裂紋前緣單元分布示意圖Fig.1 Element face arrangements across the crack front

圖2 裂紋前緣兩種單元面片示意圖Fig.2 Two kinds of arrangements of paring finite element faces

線彈性斷裂力學(xué)范圍內(nèi)，和的關(guān)系如式（3）所示。

由此可得單元面片S（=1，2）對(duì)應(yīng)的SIFs計(jì)算式：

應(yīng)變能釋放率總是假設(shè)為正，注意到式（1）、式（2）和式（4）中和可能為負(fù)，此處的正負(fù)號(hào)表示裂尖局部坐標(biāo)系下裂紋的變形方向。

1.2 節(jié)點(diǎn)力和裂紋張開位移的計(jì)算

裂紋張開位移的計(jì)算較簡單，裂紋面上表面的節(jié)點(diǎn)位移減去下表面對(duì)應(yīng)節(jié)點(diǎn)的位移，可得全局坐標(biāo)系下的裂紋張開位移，再將其轉(zhuǎn)換成裂尖局部坐標(biāo)系下的位移分量，即可得Ⅰ、Ⅱ和Ⅲ型對(duì)應(yīng)的裂紋張開位移。

2 三維裂紋自動(dòng)化建模

應(yīng)用ABAQUS軟件分析裂紋問題的步驟包括：裂紋幾何特征創(chuàng)建，材料屬性賦值，分析步創(chuàng)建，裂紋屬性賦值，創(chuàng)建邊界條件和載荷施加，網(wǎng)格劃分，定義輸出，提交計(jì)算，斷裂參數(shù)計(jì)算及結(jié)果保存等，這個(gè)過程需要非常繁瑣的手工操作。ABAQUS軟件提供了豐富的二次開發(fā)接口，所有在ABAQUS/CAE中實(shí)現(xiàn)的前處理和后處理均可通過Python程序?qū)崿F(xiàn)自動(dòng)化操作。本文以二次四面體單元VCCM法為基礎(chǔ)，應(yīng)用Python編程語言，結(jié)合開源科學(xué)計(jì)算庫NumPy和SciPy，采用面向?qū)ο蟮某绦蛟O(shè)計(jì)技術(shù)，開發(fā)了三維裂紋分析程序包。該程序包由一系列函數(shù)、類和模塊組成，其工作流程和在ABAQUS/CAE中分析裂紋問題的工作流程一致。上述裂紋分析步驟通過調(diào)用相應(yīng)的函數(shù)或?qū)嵗鄳?yīng)的類來完成，各建模步驟需要的參數(shù)通過不同的接口文件給定。

2.1 裂紋幾何特征創(chuàng)建

FRANC3D軟件通過其特有的“裂紋模板”功能來完成裂紋前緣單元布置和裂紋快速插入。ZENCRACK軟件通過“Crack-Block”技術(shù)引入不 同 形 式 的 裂 紋。J.C.Sobotka等提 出 在ABAQUS軟件上通過Python語言二次開發(fā)生成橢圓裂紋前緣網(wǎng)格和裂紋快速插入的方法。上述技術(shù)和方法的目的有兩個(gè)：一是在裂紋前緣生成高質(zhì)量單元以提高斷裂參數(shù)計(jì)算精度；二是在待分析結(jié)構(gòu)上快速引入裂紋。本文通過顯式創(chuàng)建裂紋前緣單元面片、布種控制單元分布來完成基于二次四面體單元VCCM法的三維裂紋建模，以此為基礎(chǔ)開發(fā)參數(shù)化裂紋模板庫，用于不同形式、不同尺寸的裂紋幾何特征的快速引入。

待分析結(jié)構(gòu)中引入裂紋幾何特征可分為三步：①讀取裂紋定義文件里的裂紋形式、尺寸、裂紋前緣單元數(shù)量、單元寬度、裂紋插入位置和方位等參數(shù)；②根據(jù)第①步中的裂紋參數(shù)創(chuàng)建裂紋幾何部件；③根據(jù)第①步中的裂紋位置參數(shù)通過平移、旋轉(zhuǎn)等操作，將裂紋部件移動(dòng)至指定位置和待分析結(jié)構(gòu)進(jìn)行切割、合并等幾何布爾運(yùn)算，完成裂紋幾何特征的引入。橢圓裂紋是一種非常重要的裂紋形式，在飛機(jī)結(jié)構(gòu)的損傷容限評(píng)定中，對(duì)于結(jié)構(gòu)內(nèi)部、表面、孔邊等部位的初始損傷通常用橢圓/圓形裂紋、二分之一橢圓/圓形裂紋、四分之一橢圓/圓形角裂紋來表示。本文以橢圓裂紋為例來說明參數(shù)化裂紋模板的實(shí)現(xiàn)過程。

設(shè)為全局坐標(biāo)系，坐標(biāo)系原點(diǎn)坐標(biāo)為（0，0，0），基底向量為=(1，0，0)，=(0，1，0)，=(0，0，1)，裂紋初始位置位于全局坐標(biāo)系平面上，橢圓中心位于坐標(biāo)系原點(diǎn)，裂紋面坐標(biāo)系記為，如圖3（a）所示。裂紋前緣上一點(diǎn)的坐標(biāo)（，，）可由式（7）得到：

式中：，為橢圓的短半軸和長半軸尺寸；為點(diǎn)橢圓離心角，∈(-π，π]。

圖3 橢圓裂紋前緣面片頂點(diǎn)Fig.3 Vertices of element faces across the crack front

由此，對(duì)橢圓上的離散點(diǎn)序列（，，，…）的每一個(gè)點(diǎn)執(zhí)行同樣的操作，可得位于裂紋面上的點(diǎn)序列（，，，…）和位于裂紋虛擬擴(kuò)展面上的點(diǎn)序列（，，，…）的坐標(biāo)，對(duì)這三組點(diǎn)進(jìn)行編號(hào)，然后將這三組點(diǎn)按逆時(shí)針連接形成裂紋前緣三角面片，同時(shí)對(duì)面片進(jìn)行編號(hào)并記錄面片和點(diǎn)的拓?fù)溥B接關(guān)系，然后在ABAQUS中根據(jù)前面得到的點(diǎn)和三角面片連接關(guān)系及點(diǎn)坐標(biāo)創(chuàng)建裂紋幾何部件，如圖4所示。需注意，裂紋前緣三角面片應(yīng)具有良好的形狀以避免在后續(xù)的模型網(wǎng)格生成中產(chǎn)生低質(zhì)量四面體單元。

圖4 橢圓裂紋模板Fig.4 Elliptical crack front template

2.2 裂紋插入及裂紋相關(guān)數(shù)據(jù)信息

裂紋幾何部件創(chuàng)建完成后，根據(jù)裂紋定義接口文件里的裂紋位置、方位等信息，在ABAQUS軟件的Assembly模塊里，通過平移、旋轉(zhuǎn)等操作將裂紋部件放置到指定位置，同時(shí)計(jì)算旋轉(zhuǎn)、平移操作后的裂紋前緣三角面片頂點(diǎn)的坐標(biāo)，然后和待分析結(jié)構(gòu)件進(jìn)行切割、合并等幾何布爾運(yùn)算完成裂紋幾何特征的創(chuàng)建。對(duì)于表面裂紋或穿透裂紋，由于裂紋前緣三角面片是預(yù)先產(chǎn)生的，布爾運(yùn)算后可能會(huì)在自由表面附近產(chǎn)生不規(guī)則的、影響后續(xù)網(wǎng)格劃分的小邊和小面，需要進(jìn)一步做小面刪除、短邊合并、新面添加等操作，完成后對(duì)裂紋前緣三角面片頂點(diǎn)坐標(biāo)、面片編號(hào)、面片面積等數(shù)據(jù)信息進(jìn)行更新。

裂紋平面初始位置位于全局坐標(biāo)系的平面上，裂紋面局部坐標(biāo)系和全局坐標(biāo)系重合。一個(gè)空間坐標(biāo)系可以用另一個(gè)參考坐標(biāo)系的三次空間旋轉(zhuǎn)來表達(dá)，由此可知，空間任意位置的裂紋可通過將初始裂紋最多旋轉(zhuǎn)三次來實(shí)現(xiàn)。假設(shè)裂紋插入位置的坐標(biāo)為（t，t，t），繞裂紋面局部坐標(biāo)系軸、軸、軸依次旋轉(zhuǎn)、和，則對(duì)應(yīng)的坐標(biāo)轉(zhuǎn)換矩陣為

其中，R、R、R、T的定義如下：

由此，（，，）為某點(diǎn)的初始坐標(biāo)，則經(jīng)旋轉(zhuǎn)平移后的坐標(biāo)(，，)可由式（10）計(jì)算得到。

2.3 網(wǎng)格劃分

四面體單元可以方便的自動(dòng)化生成，在有限元軟件ABAQUS中對(duì)模型進(jìn)行網(wǎng)格劃分的步驟主要包括布種、網(wǎng)格生成技術(shù)和算法指定、單元形狀和類型選擇、單元尺寸增長比選擇等。布種包括全模型布種和局部邊或面布種，兩者分別從總體和局部描述了模型網(wǎng)格密度。讀取網(wǎng)格設(shè)置接口文件里的全模型布種尺寸和局部邊布種尺寸及邊的查找坐標(biāo)，完成全局布種和局部邊的布種，裂紋前緣三角面片的每條邊上強(qiáng)制布置一個(gè)單元，單元類型設(shè)定為二次四面體單元（C3D10），網(wǎng)格生成技術(shù)為自由網(wǎng)格劃分，算法選擇默認(rèn)網(wǎng)格生成算法。在關(guān)心部位的網(wǎng)格密度足夠的情況下，單元尺寸增長比可以有效控制模型總體單元規(guī)模。根據(jù)經(jīng)驗(yàn)，其他網(wǎng)格參數(shù)保持不變的情況下，單元尺寸增長比取1.0時(shí)生成的模型四面體單元數(shù)量是取1.5時(shí)生成單元數(shù)量的3～4倍。單元尺寸增長比設(shè)置過小會(huì)導(dǎo)致網(wǎng)格規(guī)模偏大，過大則可能會(huì)產(chǎn)生低質(zhì)量單元進(jìn)而影響計(jì)算結(jié)果。為了兼顧網(wǎng)格質(zhì)量和計(jì)算效率，通過迭代的方式?jīng)Q定單元尺寸增長比的大小。具體思路如下：網(wǎng)格劃分其余設(shè)置不變，單元尺寸增長初值設(shè)為1.5，每次迭代初值減0.1直至1.0停止循環(huán)，每次循環(huán)對(duì)模型進(jìn)行一次網(wǎng)格劃分，劃分完成后對(duì)網(wǎng)格質(zhì)量進(jìn)行評(píng)估，如果該次循環(huán)沒有警告單元，跳出循環(huán)，模型網(wǎng)格劃分完成；如果迭代至循環(huán)結(jié)束每一次循環(huán)均有警告單元，取警告單元數(shù)量最少的網(wǎng)格為最終的模型網(wǎng)格。

2.4 節(jié)點(diǎn)集合和單元集合創(chuàng)建

由上述分析可知，計(jì)算所需的節(jié)點(diǎn)力需要用裂紋面上方的單元來計(jì)算，裂紋張開位移需要裂紋上表面位移節(jié)點(diǎn)和下表面對(duì)應(yīng)節(jié)點(diǎn)的位移差來計(jì)算，此過程均需判斷四面體單元和面片的位置關(guān)系。和節(jié)點(diǎn)力節(jié)點(diǎn)（圖2中面片上的節(jié)點(diǎn)1，2，3，4，5和面片上的節(jié)點(diǎn)1，2，3）關(guān)聯(lián)的四面體單元可分為3種情況：①四面體單元一個(gè)頂點(diǎn)和三角面片的頂點(diǎn)重合；②四面體單元的一條單元邊和三角面片的邊重合；③四面體單元的一個(gè)單元面和三角面片重合。如果根據(jù)四面體單元的頂點(diǎn)坐標(biāo)去判斷四面體單元和三角面片的位置關(guān)系，需要進(jìn)行較為繁瑣的情況甄別和頂點(diǎn)選取，而四面體單元的內(nèi)切球球心必定和四面體單元位于面片的同一側(cè)，本文通過判斷四面體內(nèi)切球心和三角面片的位置來確定四面體單元和面片的位置關(guān)系。

點(diǎn)和面的位置關(guān)系可以通過判斷該點(diǎn)和面內(nèi)三點(diǎn)坐標(biāo)組成的行列式的符號(hào)來判定。設(shè)三維空間點(diǎn)，，，，坐標(biāo)分別為(a，a，a)，(b，b，b)，(c，c，c)，(d，d，d)，定義行列式：

式（11）的幾何意義為，，三點(diǎn)逆時(shí)針排列，當(dāng)det(，，，)＜0時(shí)，點(diǎn)位于過、三點(diǎn)平面的上方；當(dāng)det(，，，)＞0時(shí)，點(diǎn)位于過、三點(diǎn)平面的下方。

由此可以得出四面體單元和三角面片的位置關(guān)系判斷方法：①根據(jù)四面體單元的頂點(diǎn)坐標(biāo)計(jì)算四面體單元的內(nèi)切球球心坐標(biāo)；②將三角面片的三個(gè)頂點(diǎn)坐標(biāo)和四面體球心坐標(biāo)帶入式（11），根據(jù)行列式的符號(hào)即可得四面體單元和面片的位置關(guān)系。

以面片為例來說明載荷面片相關(guān)數(shù)據(jù)信息的建立過程：①根據(jù)三角面片的頂點(diǎn)坐標(biāo)獲得載荷節(jié)點(diǎn)（圖2中面片上1，2，3，4，5節(jié)點(diǎn)）的坐標(biāo)；②根據(jù)坐標(biāo)查找網(wǎng)格節(jié)點(diǎn)，獲得網(wǎng)格節(jié)點(diǎn)編號(hào)，同時(shí)計(jì)算該節(jié)點(diǎn)的載荷分配系數(shù)；③根據(jù)單元和節(jié)點(diǎn)連接關(guān)系獲得包含該節(jié)點(diǎn)的所有四面體單元；④根據(jù)四面體單元和三角面片位置判斷方法對(duì)第③步中獲得的四面體單元進(jìn)行過濾，保留位于面片上方的單元。

對(duì)每一個(gè)位移面片和載荷面片按照上述方法建立計(jì)算所需的數(shù)據(jù)信息，所有的位移面片上的位移節(jié)點(diǎn)取并集得到節(jié)點(diǎn)集合Set-NodesU，所有載荷面片上節(jié)點(diǎn)力計(jì)算相關(guān)的單元取并集得到單元集合Set-ElemsF，這兩個(gè)集合將用于定義輸出位移和節(jié)點(diǎn)力。

2.5 輸出定義

模型創(chuàng)建完成后，調(diào)用ABAQUS/Standard生成模型的*.inp文件，通過修改*.inp文件定義輸出計(jì)算所需的相關(guān)節(jié)點(diǎn)的位移和節(jié)點(diǎn)力至*.dat文件。定義輸出節(jié)點(diǎn)力和位移的命令如下：

*EL PRINT，ELSET=Set-ElemsF

NFORC

*NODE PRINT，NSET=Set-NodesU

U

2.6 其他步驟

材料屬性賦值：讀取材料屬性定義文件里的材料屬性參數(shù)，調(diào)用函數(shù)完成模型材料屬性賦值。

裂紋屬性賦值：有限元法以無應(yīng)力面的形式來模擬裂紋，ABAQUS軟件里通過給面片或邊賦Seam屬性來定義裂紋，生成網(wǎng)格時(shí)在Seam面片上產(chǎn)生重疊的重復(fù)節(jié)點(diǎn)，在分析過程中節(jié)點(diǎn)分離來模擬裂紋受載。根據(jù)存儲(chǔ)的裂紋面相關(guān)面片坐標(biāo)信息，對(duì)幾何面片賦Seam屬性。

邊界條件和載荷：根據(jù)載荷和邊界條件定義文件里的相關(guān)參數(shù)，完成載荷施加和邊界條件創(chuàng)建。

提交計(jì)算：模型*.inp文件修改完成后，調(diào)用函數(shù)提交計(jì)算。

斷裂參數(shù)計(jì)算及結(jié)果保存：有限元模型求解完成后，讀取輸出文件*.dat里相關(guān)節(jié)點(diǎn)的位移和節(jié)點(diǎn)力，根據(jù)前述的VCCM法理論并應(yīng)用光順技術(shù)計(jì)算裂紋前緣SIFs，計(jì)算完成后將計(jì)算結(jié)果和對(duì)應(yīng)位置的坐標(biāo)、局部坐標(biāo)系等相關(guān)信息保存至文件以供下一步操作。

3 數(shù)值算例

通過數(shù)值算例對(duì)VCCM法計(jì)算斷裂參數(shù)的影響因素進(jìn)行分析，在此基礎(chǔ)上，對(duì)表面裂紋的SIFs進(jìn)行計(jì)算。

3.1 VCCM法影響因素分析

應(yīng)用四面體單元計(jì)算得到的斷裂參數(shù)沿裂紋前緣 會(huì)出現(xiàn) 震蕩，H.Okada等采用 光順技術(shù)來解決這一問題。光順技術(shù)的核心思想是以一小段裂紋前緣上的均值來取代單個(gè)面片的值。本文以Ⅰ型和復(fù)合型內(nèi)埋橢圓裂紋為例對(duì)光順面片數(shù)量、單元寬度Δ以及單元尺寸對(duì)VCCM法計(jì)算斷裂參數(shù)的影響進(jìn)行分析。

3.1.1 光順面片數(shù)量

H.Okada等的研究表明，隨著變大，沿裂紋前緣計(jì)算結(jié)果波動(dòng)變小，但當(dāng)進(jìn)一步取至21時(shí)，裂紋前緣曲率變化大的位置的計(jì)算結(jié)果和理論解的相對(duì)誤差擴(kuò)大。原因是值越大參與光順的面片覆蓋的裂紋段就越長，裂紋前緣曲率變化大的裂紋段上變化比較大，而光順?biāo)惴ǖ暮诵氖侨【担纱藢?dǎo)致數(shù)值計(jì)算結(jié)果和真實(shí)值出現(xiàn)偏差。因此，值的選擇傾向于相對(duì)小的值以避免光順范圍內(nèi)的裂紋段上結(jié)果有明顯波動(dòng)，文獻(xiàn)［11］給出的建議是取5～15。然而，當(dāng)裂紋前緣某段曲率變化很大時(shí)，即使取相對(duì)小的值也不能解決上述問題。以文獻(xiàn)［11］中的Ⅰ型內(nèi)埋橢圓裂紋算例為例，取9，當(dāng)取較大值時(shí)，VCCM法計(jì)算結(jié)果和理論解最大相對(duì)誤差均小于2%；當(dāng)取0.2時(shí)，VCCM法的計(jì)算結(jié)果和理論解的最大相對(duì)誤差擴(kuò)大至3.9%。文獻(xiàn)［11］沒有給出取更小值時(shí)VCCM法計(jì)算結(jié)果和理論解的對(duì)比，但有理由相信兩者的誤差會(huì)更大。該問題其中一個(gè)改善方法是沿曲線裂紋前緣非等尺寸布置單元，曲率大的地方單元尺寸適當(dāng)減小，如圖4所示。由此在特定的值下，裂紋前緣曲率大的地方參與光順計(jì)算的面片覆蓋的裂紋段減小，該裂紋段上的波動(dòng)范圍也會(huì)相對(duì)減小。

承受遠(yuǎn)端拉伸載荷內(nèi)埋斜置橢圓裂紋示意圖如圖5所示，當(dāng)=0時(shí)，裂紋尖端承受Ⅰ型載荷。模型尺寸=10，可視為無限大體；取0，裂紋前緣單元數(shù)取360，單元寬度Δ取裂紋前緣 平 均 單 元 尺 寸/，為 裂 紋 前 緣 長度，/分別取0.1，0.2，0.4，0.5，0.6，0.8，1.0，模型單元數(shù)為113 590～173 198，節(jié)點(diǎn)數(shù)為159 272～239 448，/0.5時(shí)內(nèi)埋橢圓裂紋有限元網(wǎng)格如圖6所示。需注意，大多數(shù)情況下，Δ取都能獲得形狀良好的裂紋前緣三角面片，但當(dāng)/取很小的值（比如0.1）時(shí)，有可能會(huì)產(chǎn)生自相交而導(dǎo)致建模失敗，此時(shí)，可通過改變裂紋前緣單元數(shù)或減小單元寬度Δ來解決，這一過程由程序自動(dòng)完成。

圖5 遠(yuǎn)端拉伸內(nèi)埋斜置橢圓裂紋示意圖Fig.5 The problem of embedded elliptical crack

圖6 a/c=0.5時(shí)內(nèi)埋橢圓裂紋有限元網(wǎng)格Fig.6 A typical finite element model for the embedded elliptical crack problem for the case of a/c=0.5

=0.5，取不 同 值 時(shí)VCCM法 計(jì) 算結(jié) 果和理論解的對(duì)比如表1所示，為了更清晰地展示取不同對(duì)VCCM法的影響，表1中的誤差分析保留三位小數(shù)，后續(xù)的誤差分析均保留兩位小數(shù)。

表1 a/c=0.5時(shí)Ne取不同值VCCM法計(jì)算結(jié)果誤差分析Table 1 Summary of error for modeⅠembedded elliptical crack with different Ne computing SIFs

從表1可以看出：=2時(shí)，計(jì)算結(jié)果和理論解的平均誤差和最大誤差分別為0.507%和1.639%，但波動(dòng)明顯，如圖7（a）所示；隨著進(jìn)一步變大，計(jì)算結(jié)果波動(dòng)性變小，和理論解的相對(duì)誤差變小，當(dāng)取至30時(shí)，計(jì)算結(jié)果和理論解依然吻合非常好，平均誤差為0.145%，最大相對(duì)誤差為0.360%，如圖7（b）所示，并沒有出現(xiàn)文獻(xiàn)［11］中誤差明顯擴(kuò)大的現(xiàn)象，說明在當(dāng)前的裂紋前緣網(wǎng)格布置下（如圖4和圖6所示），VCCM法計(jì)算斷裂參數(shù)的穩(wěn)定性有提高。綜合各方面因素，本文后續(xù)的分析中，取10。

圖7 Ne取不同值VCCM法計(jì)算結(jié)果和理論解的對(duì)比Fig.7 The distributions of the stress intensity factor along the crack front for different Ne

/取不同值時(shí)VCCM法計(jì)算結(jié)果和理論解的誤差分析如表2所示，可以 看出：當(dāng)/≥0.4時(shí)，VCCM法計(jì)算結(jié)果和理論解吻合非常好，平均誤差小于0.24%，最大誤差小于0.71%；隨著/變小，誤差稍有擴(kuò)大，但依然和理論解吻合很好，當(dāng)=0.2和=0.1時(shí)，平均誤差分別為0.53%和1.64%，最大相對(duì)誤差分別為1.02%和3.03%。誤差變大的原因是小的值靠近長軸的裂紋前緣單元尺寸小，單元寬度取定值導(dǎo)致的單元形狀不夠優(yōu)良。

表2 Ⅰ型內(nèi)埋橢圓裂紋誤差分析Table 2 Summary of error for modeⅠembedded elliptical crack with fixed-width element faces

3.1.2 單元寬度

由前面的分析可以看出，當(dāng)取很小的值，比如0.1時(shí)，VCCM法的計(jì)算結(jié)果和理論解的誤差稍有擴(kuò)大，原因是裂紋前緣的單元形狀不夠優(yōu)良。本節(jié)裂紋前緣面片采用變寬度，單元寬度等于單元尺寸，如圖8所示。模型尺寸=10，裂紋前緣單元數(shù)取360，/分別取0.1，0.2，0.4，0.5，0.6，0.8，1.0，模 型 單 元 數(shù) 為94 414～131 068，節(jié)點(diǎn)數(shù)為133 788～183 942。VCCM法計(jì)算結(jié)果和理論解的誤差分析如表3所示，可以看出：/=0.1時(shí)，相較于裂紋前緣單元寬度取定值，計(jì)算結(jié)果和理論解的平均誤差和最大誤差均有改善；/=0.2時(shí)，平均誤差有小幅改善，但最大相對(duì)誤差稍有擴(kuò)大；/≥0.4時(shí)，兩種裂紋前緣單元布置的數(shù)值計(jì)算結(jié)果很接近。

兩種裂紋前緣單元布置的計(jì)算結(jié)果說明單元寬度對(duì)VCCM法的影響很小，裂紋前緣單元具有良好的形狀即可獲得非常好的結(jié)果。實(shí)際上，對(duì)于真實(shí)裂紋，小的值對(duì)應(yīng)的裂紋形狀是不穩(wěn)定的，在疲勞載荷的作用下，裂紋前緣陡峭的地方隨著裂紋的擴(kuò)展會(huì)快速變得平緩。

圖8 變寬度單元面片的橢圓裂紋模板Fig.8 Elliptical crack front template with changing element width

表3 Ⅰ型內(nèi)埋橢圓裂紋誤差分析Table 3 Summary of error for modeⅠembedded elliptical crack with changing-width element faces

3.1.3 單元尺寸

本節(jié)對(duì)復(fù)合型內(nèi)埋橢圓裂紋的SIFs進(jìn)行計(jì)算并和理論解進(jìn)行對(duì)比，研究單元尺寸對(duì)VCCM法計(jì)算斷裂參數(shù)的影響。模型尺寸=10，模型尺寸大于10倍裂紋尺寸，可視為無限大體，取π/4，裂 紋 前 緣 單 元 數(shù)分 別 取180和360，分別取0.4，0.6，0.8，裂紋前緣單元寬度Δ取定值。取180時(shí)，模型單元數(shù)為56 644～65 957，節(jié)點(diǎn)數(shù)為82 847～95 730；取360時(shí)，模型單元數(shù)為87 619～128 285，節(jié)點(diǎn)數(shù)為125 930～180 278。復(fù)合型內(nèi)埋橢圓裂紋誤差分析如表4所示。

表4 復(fù)合型內(nèi)埋橢圓裂紋誤差分析Table 4 Summary of error for mixed mode embedded elliptical crack

從表4可以看出：隨著裂紋前緣網(wǎng)格加密，VCCM法計(jì)算的Ⅰ、Ⅱ、Ⅲ型應(yīng)力強(qiáng)度因子、、和理論解的平均誤差均減小，、和理論解最大誤差同樣減小，和理論解的最大誤差反而有小幅增大。說明隨著網(wǎng)格加密，裂紋前緣上絕大多數(shù)位置的數(shù)值計(jì)算結(jié)果和理論解的誤差均減小，個(gè)別位置的值和理論解誤差稍有擴(kuò)大。=180時(shí)，VCCM法計(jì)算的、、和理論解的平均誤差為0.24%～0.35%、0.56%～1.21%、0.47%～0.68%；=360時(shí)，VCCM法計(jì)算的、、和理論解平均誤差為0.23%～0.27%、0.36%～0.57%、0.27%～0.32%。表明兩種網(wǎng)格均得到了很好的結(jié)果，同時(shí)也說明單元尺寸對(duì)VCCM法的影響有限。

3.2 Ⅰ型表面裂紋應(yīng)力強(qiáng)度因子計(jì)算

圖9 有限寬板表面裂紋Fig.9 The problem of semi-elliptical surface crack

圖10 a/c分別取0.6和0.8時(shí)裂紋前緣應(yīng)力強(qiáng)度因子分布Fig.10 The distributions of the SIFs for a/c being 0.6 and 0.8

4 結(jié)論

（1）通過顯式創(chuàng)建裂紋前緣單元面片結(jié)合布種控制單元分布實(shí)現(xiàn)了基于二次四面體單元VCCM法的三維裂紋建模。全模型采用常規(guī)四面體單元顯著簡化了建模流程，提高了三維裂紋自動(dòng)化建模的穩(wěn)定性和魯棒性。

（2）本文開發(fā)了基于二次四面體單元VCCM法的參數(shù)化三維裂紋分析程序包，實(shí)現(xiàn)了三維裂紋從建模到斷裂參數(shù)計(jì)算的全流程自動(dòng)化，有效提升了三維裂紋問題的分析效率。

（3）單元寬度和單元尺寸對(duì)二次四面體單元VCCM法影響有限，數(shù)值算例結(jié)果和文獻(xiàn)結(jié)果一致性良好。

（4）本文的建模方法同樣可用于其他有限元軟件的三維裂紋全自動(dòng)化建模。此外，本文的程序包可用于三維裂紋疲勞擴(kuò)展的全自動(dòng)化數(shù)值模擬。