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

高文，王生楠

（西北工業大學航空學院，西安 710072）

0 引言

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

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

1 虛擬裂紋閉合法

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

由VCCM法的理論可知，在裂紋的虛擬擴展面上需要有單元面片存在以獲得計算所需的節點力，且裂紋前緣兩側的面片寬度需相等。四面體單元的面是三角形，典型的裂紋前緣二次四面體單元分布示意圖如圖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

線彈性斷裂力學范圍內，和的關系如式（3）所示。

由此可得單元面片S（=1，2）對應的SIFs計算式：

應變能釋放率總是假設為正，注意到式（1）、式（2）和式（4）中和可能為負，此處的正負號表示裂尖局部坐標系下裂紋的變形方向。

1.2 節點力和裂紋張開位移的計算

裂紋張開位移的計算較簡單，裂紋面上表面的節點位移減去下表面對應節點的位移，可得全局坐標系下的裂紋張開位移，再將其轉換成裂尖局部坐標系下的位移分量，即可得Ⅰ、Ⅱ和Ⅲ型對應的裂紋張開位移。

2 三維裂紋自動化建模

應用ABAQUS軟件分析裂紋問題的步驟包括：裂紋幾何特征創建，材料屬性賦值，分析步創建，裂紋屬性賦值，創建邊界條件和載荷施加，網格劃分，定義輸出，提交計算，斷裂參數計算及結果保存等，這個過程需要非常繁瑣的手工操作。ABAQUS軟件提供了豐富的二次開發接口，所有在ABAQUS/CAE中實現的前處理和后處理均可通過Python程序實現自動化操作。本文以二次四面體單元VCCM法為基礎，應用Python編程語言，結合開源科學計算庫NumPy和SciPy，采用面向對象的程序設計技術，開發了三維裂紋分析程序包。該程序包由一系列函數、類和模塊組成，其工作流程和在ABAQUS/CAE中分析裂紋問題的工作流程一致。上述裂紋分析步驟通過調用相應的函數或實例化相應的類來完成，各建模步驟需要的參數通過不同的接口文件給定。

2.1 裂紋幾何特征創建

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

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

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

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

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

由此，對橢圓上的離散點序列（，，，…）的每一個點執行同樣的操作，可得位于裂紋面上的點序列（，，，…）和位于裂紋虛擬擴展面上的點序列（，，，…）的坐標，對這三組點進行編號，然后將這三組點按逆時針連接形成裂紋前緣三角面片，同時對面片進行編號并記錄面片和點的拓撲連接關系，然后在ABAQUS中根據前面得到的點和三角面片連接關系及點坐標創建裂紋幾何部件，如圖4所示。需注意，裂紋前緣三角面片應具有良好的形狀以避免在后續的模型網格生成中產生低質量四面體單元。

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

2.2 裂紋插入及裂紋相關數據信息

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

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

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

由此，（，，）為某點的初始坐標，則經旋轉平移后的坐標(，，)可由式（10）計算得到。

2.3 網格劃分

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

2.4 節點集合和單元集合創建

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

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

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

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

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

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

2.5 輸出定義

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

*EL PRINT，ELSET=Set-ElemsF

NFORC

*NODE PRINT，NSET=Set-NodesU

U

2.6 其他步驟

材料屬性賦值：讀取材料屬性定義文件里的材料屬性參數，調用函數完成模型材料屬性賦值。

裂紋屬性賦值：有限元法以無應力面的形式來模擬裂紋，ABAQUS軟件里通過給面片或邊賦Seam屬性來定義裂紋，生成網格時在Seam面片上產生重疊的重復節點，在分析過程中節點分離來模擬裂紋受載。根據存儲的裂紋面相關面片坐標信息，對幾何面片賦Seam屬性。

邊界條件和載荷：根據載荷和邊界條件定義文件里的相關參數，完成載荷施加和邊界條件創建。

提交計算：模型*.inp文件修改完成后，調用函數提交計算。

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

3 數值算例

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

3.1 VCCM法影響因素分析

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

3.1.1 光順面片數量

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

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

圖5 遠端拉伸內埋斜置橢圓裂紋示意圖Fig.5 The problem of embedded elliptical crack

圖6 a/c=0.5時內埋橢圓裂紋有限元網格Fig.6 A typical finite element model for the embedded elliptical crack problem for the case of a/c=0.5

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

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

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

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

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

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

3.1.2 單元寬度

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

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

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

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

3.1.3 單元尺寸

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

表4 復合型內埋橢圓裂紋誤差分析Table 4 Summary of error for mixed mode embedded elliptical crack

從表4可以看出：隨著裂紋前緣網格加密，VCCM法計算的Ⅰ、Ⅱ、Ⅲ型應力強度因子、、和理論解的平均誤差均減小，、和理論解最大誤差同樣減小，和理論解的最大誤差反而有小幅增大。說明隨著網格加密，裂紋前緣上絕大多數位置的數值計算結果和理論解的誤差均減小，個別位置的值和理論解誤差稍有擴大。=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%。表明兩種網格均得到了很好的結果，同時也說明單元尺寸對VCCM法的影響有限。

3.2 Ⅰ型表面裂紋應力強度因子計算

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

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

4 結論

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

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

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

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