王海鵬，吳永東

(南昌航空大學，南昌 330063)

0 引言

無人機(Unmanned Aerial Vehicle，UAV)是一種由動力驅(qū)動、機上無人駕駛、可重復使用的航空器的簡稱。具有隱蔽性好，生命力強;造價低廉，不懼傷亡;起降簡單，操作靈活等優(yōu)點［1］。機翼各部件中存在的微小裂紋，在飛機的使用過程中會不斷擴展，最終導致機翼的斷裂破壞，嚴重威脅飛行安全。

裂紋擴展屬于典型的強不連續(xù)問題，一直是工程力學分析的熱點、重點和難點問題［2］。在模擬裂紋擴展時，Shephard M等［3］在每次裂紋擴后把裂紋作為邊界條件的一部分，然后重新劃分網(wǎng)格。Klein P等［4］引用可虛擬聯(lián)接單元的概念，用虛擬聯(lián)接模型來解決裂紋擴展問題。王慰軍［5］編寫了基于ABAQUS平臺的裂紋擴展模擬軟件。Ortiz M等［6］發(fā)展了三維有限變形的內(nèi)聚元和一系列不可逆內(nèi)聚力關系，該方法可以近似地模擬出裂尖的軌跡。婁路亮等［7］利用無網(wǎng)格法模擬了裂紋的擴展，得到的裂紋擴展路徑與實驗值吻合的很好。Nicolas Moёs等［8］利用擴展有限元法對裂紋擴展進行模擬，該方法允許裂紋與網(wǎng)格相對獨立，因此在模擬過程中無需重新劃分網(wǎng)格。目前，國內(nèi)外對機翼的裂紋擴展研究較少，對機翼裂紋的擴展過程與形式認識不足。

研究采用Python語言對ABAQUS進行二次開發(fā)，實現(xiàn)模擬裂紋自動擴展過程中有限元網(wǎng)格的自動重劃分功能，采用最大周向應力準則，計算裂紋擴展方向，應力強度因子等斷裂參量。將機翼肋板簡化為二維平板結構，在肋板與梁的連接處放置裂紋，模擬不同初始角度的裂紋的擴展過程，得到相關的力學參量，總結裂紋擴展路徑的總體趨勢。

1 裂紋擴展判斷準則

在實際問題中，裂紋往往處于復合變形狀態(tài)，必須采用復合型斷裂準則，試驗采用最大周向應力準則?；炯僭O是:(1)開裂方向是最大周向應力σθmax所在的方向;(2)該方向的周向應力達到臨界值時，裂紋初始擴展。

裂紋尖端的周向應力為［9］:

式中:r為徑向坐標，θ為角坐標，KⅠ和KⅡ分別為Ⅰ和Ⅱ型的應力強度因子。

裂紋擴展角度θ0可由式(1)對θ的一階偏導得到，注意到當裂紋為純Ⅰ型時，K=0，θ=Ⅱ00，可得

當沿θ0方向的周向應力達到σθC時，裂紋失穩(wěn)擴展，即

臨界值σθC一般由Ⅰ型開裂條件給出，即σθC= σθ(KⅠC，θ，θ0)，KⅠC為材料的斷裂韌度。將 θ0代入式(1)和式(2)，得到裂紋失穩(wěn)擴展條件為:

2 裂紋擴展模擬方法

2.1 裂紋自動擴展程序的設計思路及流程

斷裂力學分析是ABAQUS的一個重要的功能模塊，可以準確地計算出結構中裂紋的相關斷裂力學參數(shù)，但卻不具備模擬在未知路徑下裂紋擴展的能力，裂紋每擴展一步，都需要將草圖創(chuàng)建、裂紋定義、網(wǎng)格劃分、計算提交、計算結果提取等這些步驟重復一遍，這就導致了人工工作量十分巨大。通過ABAQUS提供的二次開發(fā)接口，利用PYTHON語言，編寫裂紋自動擴展的程序模塊，極大的提高了分析效率。

本研究的裂紋自動擴展程序主要有以下部分組成:應力分析及裂紋擴展檢測、計算分析判斷裂紋擴展方向、網(wǎng)格重劃分、計算結果寫入指定數(shù)據(jù)文件。程序的核心思路是使用PYTHON控制seam(裂紋)尺寸，然后移動partition(裂尖區(qū)域劃分)和網(wǎng)格。圖1為程序的主流程圖。

圖1 裂紋自動擴展程序流程圖Fig.1 Flow chart of crack automatic propagation program

帶裂紋的有限元模型的前處理和有限元計算分別由ABAQUS/CAE和ABAQUS/Standard完成。裂紋判斷模塊，將根據(jù)計算結果，利用最大周向應力準則判斷裂紋是否擴展并得出擴展角度。網(wǎng)格重劃分模塊將根據(jù)前一步模型中裂紋尖端的位置，通過給定的擴展增量和擴展角度，計算新裂紋尖端的位置，增加擴展方向上的裂紋，劃分裂尖區(qū)域，并重新劃分整個分析模型的網(wǎng)格，最后提交計算。這一過程會不斷重復，直到分析結果滿足一定的條件，從而實現(xiàn)對裂紋擴展的仿真。

在裂紋尖端，應力出現(xiàn)奇異點，為了確保得到的結果真實可靠，必須在裂尖放置奇異性單元。在ABAQUS中進行斷裂計算時，要在裂紋的尖端布置一圈1/4奇異性單元，用來模擬應力應變的奇異性。引入的奇異性單元，可以看作是把四邊形八結點等參單元的一條邊上的3個結點重合，使該邊退化為一個點，將與裂紋尖端相鄰的兩條邊上的中點移動到距離裂紋尖端1/4處而形成［10］。理論和實際都表明，這種奇異單元具有的奇異性，可以滿足裂紋尖端應力和應變奇異性的要求。因此，可在裂尖布置相應的1/4奇異性單元，如圖2a所示，裂尖網(wǎng)格采用圖2b所示的方式劃分。

圖2 裂尖網(wǎng)格劃分示意圖Fig.2 Diagram of crack tip mesh

2.2 程序驗證

計算模型如圖3所示，矩形板的2個圓孔的圓周上分別有一條斜裂紋，對這2條裂紋在遠場拉伸應力作用下的擴展情況進行模擬。文獻［8］利用擴展有限元法對上述問題進行了分析，本節(jié)將把分析結果與文獻［8］中的結果進行對比，以此來驗證程序的可靠性與準確性。

該裂紋板是邊長為127mm的正方形板，兩圓孔位于板的中心位置，兩圓心相距25.4mm，圓孔半徑為2mm，兩條裂紋的初始長度均為2.54mm，裂紋與水平線的初始角度均為θ=45°，裂紋板的上下兩表面均施加大小為13790 Pa的拉應力。材料的彈性模量和泊松比分別為E=689.5 MPa和 ν=0.3。文獻［8］中裂紋的擴展增量為Δa=1.27mm，本研究分別模擬了 Δa=1.02、1.27、1.52mm情況下裂紋的擴展過程以及對應的斷裂力學參數(shù)。為了便于與文獻給出的結果進行對比，均采用左側裂紋的值。

圖3 雙孔邊裂紋板Fig.3 Plate with cracks emanating from two hole

圖4是利用裂紋自動擴展程序模擬得到的裂紋擴展路徑及與文獻［8］的對比?？梢钥吹?，利用本研究的方法得到的裂紋擴展路徑與文獻給出的結果基本吻合，表明了程序的正確性。

圖5是應力強度因子KⅠ和KⅡ隨裂紋長度變化的規(guī)律。對于KⅠ，除了中間3步的結果與參考值有較大變化外，其它擴展步得到的結果都與文獻結果吻合得很好，KⅡ在裂紋的整個擴展過程中與文獻給出的值都基本符合，誤差在可接受范圍內(nèi)。

圖6是裂紋在不同步長下的擴展軌跡，由裂紋在每一步擴展過程中裂尖的坐標得出。可以看出，只要裂紋擴展的步長值在合理的范圍之內(nèi)，裂紋的擴展軌跡就不會發(fā)生較大偏差。因此，在分析實際問題時，當結構的尺寸相對于裂紋較大時，可以適當增加裂紋擴展的步長，這樣能在不影響最終結果準確性的前提下減少運算量，節(jié)省運算時間。

圖4 雙孔邊裂紋板裂紋擴展路徑模擬結果Fig.4 Simulation result of a plate with cracks emanating from two holes

圖5 應力強度因子變化規(guī)律Fig.5 Variation of stress intensity factor

圖6 不同步長下的裂紋擴展軌跡Fig.6 Path of crack propagation under different step

3 無人機機翼肋板裂紋擴展模擬

3.1 機翼肋板結構分析建模與載荷確定

翼肋的形式多種多樣，有構架式、腹板式、整體式等。不論使用何種形式，都是用來保持機翼外形、承受局部氣動載荷。翼肋的厚度遠小于其長度和寬度，而且與機翼上安置翼肋處的翼型在外形上基本一致。因此，可將翼肋簡化為二維翼型平板結構進行分析［11］。采用的翼型為NACA－2412，肋板中間的孔為翼梁所在位置。

根據(jù)文獻［12］的氣動分析結果，可將翼肋所受的載荷分段施加，這樣更符合實際情況(圖7)。具體的載荷為:在翼肋上表面從左至右分5段施加均布載荷，大小依次為 0.12、0.096、0.06、0.096、0.12 MPa;翼肋下表面載荷同樣處理，載荷大小依次為0.18、0.156、0.12、0.156、0.18MPa。計算時翼梁所在的位置采用固定約束。肋板材料為鋁合金，參數(shù)為:E=72 GPa，ρ=2780 kg/m3，μ =0.33。

裂紋取在肋板與梁的連接處，如圖8所示，其中θ為裂紋與x軸 (與肋板弦線平行)的夾角。本文模擬了初始夾角θ不同時裂紋的擴展形態(tài)，計算了每擴展一步對應的應力強度因子和裂紋擴展角。計算時設定裂紋的初始長度為3mm，裂紋擴展增量為0.2mm。

圖7 簡化翼肋的有限元模型和載荷模型Fig.7 Simplified finite element model of the wing rib and load model

圖8 裂紋位置Fig.8 Location of crack

3.2 結果分析

分別對θ=－70°～70°的裂紋擴展情況進行模擬，圖9是不同初始角度下裂紋的擴展形態(tài)?？梢钥吹?，裂紋的擴展軌跡可以分為兩種情況。當裂紋初始角度θ=－70°～40°時，裂紋最終向著翼肋的下表面擴展;當裂紋初始角度θ=－50°～70°時，裂紋沿著初始方向向前擴展，不發(fā)生方向的改變。

圖10與圖11是裂紋初始角度θ=－70°～40°時裂紋斷裂力學參量的變化規(guī)律。對于應力強度因子KⅠ，除θ=30°和θ=40°外，KⅠ值最終都以大致相同的速率增長，而對于 θ=30°和 θ=40°，KⅠ在經(jīng)過較長的擴展步之后，其值也呈現(xiàn)與其他情況相同的變化趨勢，且θ越大，對應的KⅠ值也越大。對于應力強度因子KⅡ，除了KⅡ值從第16步開始基本不再起伏變化，而其他的裂紋，從第34擴展步開始，其KⅡ值也在0上下小幅度變化。對于裂紋擴展角，θ=40°的裂紋擴展角從第34步開始穩(wěn)定在0°，裂紋擴展方向基本不再變化，對于其他的裂紋，從第16步開始裂紋擴展的方向就不再變化。

圖9 不同初始角度裂紋的擴展形態(tài)Fig.9 Crack propagation morphology of different initial crack angle

圖10 應力強度因子變化規(guī)律 (θ=－70°～40°)Fig.10 Variation of stress intensity factor(θ=－70°～40°)

圖11 裂紋擴展角變化規(guī)律 (θ=－70°～40°)Fig.11 Variation of crack propagation angle(θ=－70°～40°)

圖12是裂紋初始角度時θ=50°～70°裂紋應力強度因子的變化規(guī)律。隨著裂紋的增長，KⅠ的代數(shù)值呈下降趨勢，KⅡ的值不斷增加。對于θ=50°和θ=60°的裂紋，其應力強度因子的變化趨勢都較平緩，值的增長或減小速率大致相同;而對于θ=70°的裂紋，應力強度因子的值則呈指數(shù)變化，迅速地增長或減小。

4 結論

1)通過ABAQUS提供的二次開發(fā)接口，利用PYTHON語言，編寫裂紋自動擴展的程序，實現(xiàn)了對裂紋擴展的模擬，并與相關文獻結果進行比較，驗證了該模擬方法的有效性和準確性。

2)裂紋的擴展形態(tài)、擴展過程中的應力強度因子和裂紋擴展角的變化規(guī)律，與裂紋的初始角度有密切的關系，總體可分為兩種情況:一種是裂紋沿著初始方向擴展，不發(fā)生偏轉;另一種是裂紋在擴展的初始階段發(fā)生較大的偏轉，隨后擴展方向基本不再變化。

圖12 應力強度因子變化規(guī)律(θ=50°～70°)Fig.12 Variation of stress intensity factor(θ=50°～70°)

3)研究結果對研究機翼肋板裂紋的止裂以及裂紋擴展速率等有一定的參考價值。

［1］秦博，王蕾.無人機發(fā)展綜述［J］.飛航導彈，2002(8):4－10.

［2］吳圣川，吳玉程.斷裂分析及CAE軟件的現(xiàn)狀與發(fā)展［J］.計算機輔助工程，2011，20(1):1－2.

［3］Shephard M，Yehia N，Burd G，et al.Automatic crack propagation tracking［J］.Computers and Structures，1985，20(1):211－223.

［4］Klein P，Gao H.Crack nucleation and growth as strain localization in a virtual－bond continuum ［J］. Engineering Fracture Mechanics，1998，61(1):21－48.

［5］王慰軍.基于ABAQUS的裂紋擴展仿真軟件及應用［D］.杭州:浙江大學，2006:17－19.

［6］Ortiz M，Pandolfi A.用于三維裂紋擴展分析的有限變形不可逆內(nèi)聚元［J］.力學進展，2008，38(5):630－640.

［7］婁路亮，曾攀，聶蕾.裂紋擴展的無網(wǎng)格數(shù)值模擬方法［J］.航空材料學報，2001，21(3):51－56.

［8］Moёs N，Dolbow J，Belytschko T.A finite element method for crack growth without remeshing［J］.International Journal for Numerical Methods in Engineering，1999，46:131－150.

［9］程靳，趙樹山.斷裂力學［M］.北京:科學出版社，1996:30－31.

［10］姜薇.基于ABAQUS的船舶典型結構裂紋擴展模擬及分析研究［D］.武漢:華中科技大學，2011:14－15.

［11］季武強.飛機機翼結構拓撲優(yōu)化方法研究及其實現(xiàn)［D］.沈陽:沈陽航空航天大學，2013:39－40.

［12］郭遠帥.某無人機機翼斷裂力學分析［D］.南昌:南昌航空大學，2013:37－41.