基于四面體單元的三維裂紋自動化建模

高文，王生楠

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

0 引言

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

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

1 虛擬裂紋閉合法

1.1 二次四面體單元的VCCM法計算式

由VCCM法的理論可知，在裂紋的虛擬擴展面上需要有單元面片存在以獲得計算所需的節(jié)點力，且裂紋前緣兩側(cè)的面片寬度需相等。四面體單元的面是三角形，典型的裂紋前緣二次四面體單元分布示意圖如圖1所示。裂紋前緣的三角形單元面片分為兩種：一種是三角面片的一條邊位于裂紋前緣，記為；另一種是三角面片的一個頂點位于裂紋前緣，記為，Δ為單元寬度，如圖2所示。面片和上的計算式如式（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）對應(yīng)的SIFs計算式：

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

1.2 節(jié)點力和裂紋張開位移的計算

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

2 三維裂紋自動化建模

應(yīng)用ABAQUS軟件分析裂紋問題的步驟包括：裂紋幾何特征創(chuàng)建，材料屬性賦值，分析步創(chuàng)建，裂紋屬性賦值，創(chuàng)建邊界條件和載荷施加，網(wǎng)格劃分，定義輸出，提交計算，斷裂參數(shù)計算及結(jié)果保存等，這個過程需要非常繁瑣的手工操作。ABAQUS軟件提供了豐富的二次開發(fā)接口，所有在ABAQUS/CAE中實現(xiàn)的前處理和后處理均可通過Python程序?qū)崿F(xiàn)自動化操作。本文以二次四面體單元VCCM法為基礎(chǔ)，應(yīng)用Python編程語言，結(jié)合開源科學(xué)計算庫NumPy和SciPy，采用面向?qū)ο蟮某绦蛟O(shè)計技術(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ù)和方法的目的有兩個：一是在裂紋前緣生成高質(zhì)量單元以提高斷裂參數(shù)計算精度；二是在待分析結(jié)構(gòu)上快速引入裂紋。本文通過顯式創(chuàng)建裂紋前緣單元面片、布種控制單元分布來完成基于二次四面體單元VCCM法的三維裂紋建模，以此為基礎(chǔ)開發(fā)參數(shù)化裂紋模板庫，用于不同形式、不同尺寸的裂紋幾何特征的快速引入。

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

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

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

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

由此，對橢圓上的離散點序列（，，，…）的每一個點執(zhí)行同樣的操作，可得位于裂紋面上的點序列（，，，…）和位于裂紋虛擬擴展面上的點序列（，，，…）的坐標(biāo)，對這三組點進行編號，然后將這三組點按逆時針連接形成裂紋前緣三角面片，同時對面片進行編號并記錄面片和點的拓?fù)溥B接關(guān)系，然后在ABAQUS中根據(jù)前面得到的點和三角面片連接關(guā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)等操作將裂紋部件放置到指定位置，同時計算旋轉(zhuǎn)、平移操作后的裂紋前緣三角面片頂點的坐標(biāo)，然后和待分析結(jié)構(gòu)件進行切割、合并等幾何布爾運算完成裂紋幾何特征的創(chuàng)建。對于表面裂紋或穿透裂紋，由于裂紋前緣三角面片是預(yù)先產(chǎn)生的，布爾運算后可能會在自由表面附近產(chǎn)生不規(guī)則的、影響后續(xù)網(wǎng)格劃分的小邊和小面，需要進一步做小面刪除、短邊合并、新面添加等操作，完成后對裂紋前緣三角面片頂點坐標(biāo)、面片編號、面片面積等數(shù)據(jù)信息進行更新。

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

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

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

2.3 網(wǎng)格劃分

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

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

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

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

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

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

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

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

2.5 輸出定義

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

*EL PRINT，ELSET=Set-ElemsF

NFORC

*NODE PRINT，NSET=Set-NodesU

U

2.6 其他步驟

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

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

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

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

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

3 數(shù)值算例

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

3.1 VCCM法影響因素分析

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

3.1.1 光順面片數(shù)量

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

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

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

圖6 a/c=0.5時內(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，取不 同 值 時VCCM法 計 算結(jié) 果和理論解的對比如表1所示，為了更清晰地展示取不同對VCCM法的影響，表1中的誤差分析保留三位小數(shù)，后續(xù)的誤差分析均保留兩位小數(shù)。

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

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

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

/取不同值時VCCM法計算結(jié)果和理論解的誤差分析如表2所示，可以 看出：當(dāng)/≥0.4時，VCCM法計算結(jié)果和理論解吻合非常好，平均誤差小于0.24%，最大誤差小于0.71%；隨著/變小，誤差稍有擴大，但依然和理論解吻合很好，當(dāng)=0.2和=0.1時，平均誤差分別為0.53%和1.64%，最大相對誤差分別為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時，VCCM法的計算結(jié)果和理論解的誤差稍有擴大，原因是裂紋前緣的單元形狀不夠優(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é)點數(shù)為133 788～183 942。VCCM法計算結(jié)果和理論解的誤差分析如表3所示，可以看出：/=0.1時，相較于裂紋前緣單元寬度取定值，計算結(jié)果和理論解的平均誤差和最大誤差均有改善；/=0.2時，平均誤差有小幅改善，但最大相對誤差稍有擴大；/≥0.4時，兩種裂紋前緣單元布置的數(shù)值計算結(jié)果很接近。

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

圖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é)對復(fù)合型內(nèi)埋橢圓裂紋的SIFs進行計算并和理論解進行對比，研究單元尺寸對VCCM法計算斷裂參數(shù)的影響。模型尺寸=10，模型尺寸大于10倍裂紋尺寸，可視為無限大體，取π/4，裂 紋 前 緣 單 元 數(shù)分 別 取180和360，分別取0.4，0.6，0.8，裂紋前緣單元寬度Δ取定值。取180時，模型單元數(shù)為56 644～65 957，節(jié)點數(shù)為82 847～95 730；取360時，模型單元數(shù)為87 619～128 285，節(jié)點數(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法計算的Ⅰ、Ⅱ、Ⅲ型應(yīng)力強度因子、、和理論解的平均誤差均減小，、和理論解最大誤差同樣減小，和理論解的最大誤差反而有小幅增大。說明隨著網(wǎng)格加密，裂紋前緣上絕大多數(shù)位置的數(shù)值計算結(jié)果和理論解的誤差均減小，個別位置的值和理論解誤差稍有擴大。=180時，VCCM法計算的、、和理論解的平均誤差為0.24%～0.35%、0.56%～1.21%、0.47%～0.68%；=360時，VCCM法計算的、、和理論解平均誤差為0.23%～0.27%、0.36%～0.57%、0.27%～0.32%。表明兩種網(wǎng)格均得到了很好的結(jié)果，同時也說明單元尺寸對VCCM法的影響有限。

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

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

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

4 結(jié)論

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

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

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

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