飛機薄壁結構面外載荷下三維疲勞裂紋擴展仿真分析

曾葦鵬，魯嵩嵩，劉斌超，劉漢海，隋福成，鮑蕊

（1.北京航空航天大學 航空科學與工程學院，北京 100083）

（2.中國航空工業(yè)集團有限公司 沈陽飛機設計研究所，沈陽 110035）

0 引言

薄壁結構是飛機的基本結構形式之一，是飛機設計時用來承受面內的拉伸和剪切載荷的重要結構。然而，部分薄板結構在實際應用中，不可避免地要承受面外彎曲載荷，由于薄板結構的抗彎截面剛度很小，面外彎曲載荷會引發(fā)極大的彎曲應力，進而影響此類結構的使用壽命。例如，某型飛機服役過程中，來自鎖鉤和鎖環(huán)的集中載荷，將導致座艙蓋相關側型材薄板結構發(fā)生面外彎曲，從而承受較高水平的面外彎曲載荷。由于飛機服役中的載荷環(huán)境具有交變的特征，此類面外彎曲載荷將導致側型材薄板結構存在較大的疲勞失效風險，進而影響飛機的使用安全。因此，有必要對面外彎曲載荷下的薄板結構開展疲勞/損傷容限評定。

隨著計算機軟硬件水平的提升，數(shù)值模擬成為分析裂紋問題的有力手段。I.S.Raju等基于有限元計算結果，給出了拉伸或彎曲載荷下半橢圓表面裂紋應力強度因子隨參數(shù)角度、裂紋深度、裂紋長度、板厚和板寬變化的經驗公式，并分析了板寬對裂紋前緣應力強度因子變化的影響；X.B.Lin等采用有限元方法計算應力強度因子，模擬了裂紋擴展，并討論了敏感性和網格非正交對應力強度因子的影響。

目前，以有限元為代表的各類數(shù)值仿真方法已經基本解決了各類裂紋的應力強度因子計算問題。部分構型與承載簡單的結構的裂紋擴展模擬問題也已基本解決，但復雜裂紋的裂紋擴展模擬研究仍有待繼續(xù)深入。為了模擬復雜裂紋的裂紋擴展，研究人員開發(fā)了三維裂紋擴展分析軟件FRANC3D，且已有許多研究者利用該軟件進行了不同結構在部分承載形式下的裂紋擴展模擬研究，并驗證了FRANC3D計算的可靠性和準確性。但這些研究大多是針對單向或雙向拉伸載荷進行的，針對彎曲載荷下的薄板結構的研究較少。

本文針對面外彎曲載荷下的裂紋擴展分析方法開展研究。采用將彎矩等效成拉應力的方法，計算面外彎曲載荷下疲勞裂紋的應力強度因子，并分析該等效方法的準確性；采用FRANC3DABAQUS聯(lián)合仿真的方法，模擬疲勞裂紋擴展，得到面外彎曲載荷下板的裂紋擴展特征；在此基礎上分析影響應力強度因子計算結果的因素，討論仿真方法的適用性。

1 基于彎矩等效拉伸應力方法的薄板應力強度因子計算

圖1 NASGRO中的等效拉伸載荷方法［16］Fig.1 Equivalent tensile load method in NASGRO［16］

將此等效方法用于中心帶孔板（長230 mm、寬300 mm、厚3 mm、圓孔直徑6 mm、初始裂紋總長12 mm，如圖2所示）的應力強度因子計算。由于此板僅受彎曲載荷作用，在采用此方法計算等效拉伸應力時，僅考慮部分。對應的幾何修正因子的數(shù)值可在應力強度因子手冊中查表得到，其值接近1。

圖2 中心裂紋板尺寸及受載情況Fig.2 Size of central crack plate and the applied load

等效方法的應力強度因子計算結果與文獻［1］中的試驗結果對比如圖3所示，可以看出：等效方法的應力強度因子計算結果明顯偏大（大于100%）。因此，該等效方法不適用于彎曲應力占主導地位的應力強度因子計算。

圖3 等效公式和試驗［1］反推的應力強度因子變程Fig.3 Stress intensity factor variation from simplified formula and experiment［1］

2 仿真與分析

2.1 FRANC3D與ABAQUS聯(lián) 合 仿 真 流 程

FRANC3D與ABAQUS聯(lián)合仿真的具體流程為：在FRANC3D中導入無裂紋的有限元模型，使用FRANC3D提供的功能插入指定位置和形狀的裂紋，F(xiàn)RANC3D對插入裂紋的有限元模型重新劃分網格；再調用ABAQUS求解有限元模型，F(xiàn)RANC3D利用ABAQUS的計算結果，計算裂尖裂紋擴展驅動力，即應力強度因子，判斷是否達到斷裂韌度，如果達到了則停止裂紋擴展，如果沒達到，則根據(jù)應力強度因子計算裂紋擴展后的前緣；重復從重新劃分網格到擴展裂紋前緣的過程，使裂紋不斷擴展。仿真流程圖如圖4所示。

圖4 FRANC3D與ABAQUS聯(lián)合仿真流程Fig.4 FRANC3D and ABAQUS joint simulation process

FRANC3D的主要作用是重新劃分裂紋前緣周圍的網格、計算應力強度因子以及計算裂紋前緣擴展。FRANC3D在計算應力強度因子時，提供了-積分、虛擬裂紋閉合技術等方法。在獲取裂紋前緣各節(jié)點的應力強度因子后，F(xiàn)RANC3D采用最大周向應力準則或其他準則判斷裂紋的擴展方向，并采用Paris公式或其他裂紋擴展速率公式，計算一定循環(huán)次數(shù)后裂紋前緣各節(jié)點的擴展方向和擴展距離，并對擴展后的節(jié)點進行曲線擬合，得到擴展后的裂紋前緣。

2.2 基于有限元仿真的應力強度因子計算

中心帶孔板的有限元模型如圖5所示，模型尺寸同圖2。

圖5 中心帶孔板有限元模型Fig.5 Finite element model of plate with a central hole

為了提高計算效率，僅建立1/2模型，并在對稱面施加對稱邊界條件、兩端施加固支邊界條件。模型中板表面施加合力大小1.8 kN、外徑10 mm的環(huán)狀均勻壓強，以模擬薄板受面外法向載荷的情況。模型為線彈性本構，材料為某2系鋁合金，彈性模量70 GPa，泊松比0.30，單元類型為八節(jié)點六面體線性減縮積分單元C3D8R，網格尺寸在0.04～6.00 mm的范圍內。

在FRANC3D中，插入總長2=12 mm的初始直線前緣穿透裂紋，重新劃分網格后裂尖附近網格如圖6所示（裂尖局部區(qū)域網格較細）。裂紋擴展模型中設置應力比=0.1，擴展偏角準則選擇最大拉應力準則，裂紋擴展速率公式選用Paris公式d/d=(Δ)，其中，和為材料參數(shù)，依 據(jù) 試 驗 中 的 數(shù) 據(jù)，為2.136×10，為3.133。

圖6 裂尖附近有限元網格Fig.6 Finite element mesh around crack tip

采用FRANC3D-ABAQUS聯(lián)合方法，計算穿透裂紋的裂紋前緣在受拉面一點不同裂紋長度下的應力強度因子變程Δ，計算結果如圖7中黑色線條所示，可以看出：Δ在裂紋長度超過8 mm后，隨著受拉面裂紋長度的增加而震蕩上升，但上升速度緩慢。計算結果震蕩的原因在于，裂紋前緣不斷變化的形狀，影響了應力強度因子的計算。裂紋擴展速率上升緩慢的原因在于，模擬過程中受拉面裂紋長度大于受壓面，導致裂紋前緣與板表面夾角減小，減緩了受拉面應力強度因子的增加。

文獻［1］中的試驗反推得到的應力強度因子變程如圖7中紅色點所示，對比模擬結果與試驗結果可以看出：受拉面裂紋長度較短時（＜8 mm），F(xiàn)RANC3D計算結果明顯大于試驗結果，這是因為模擬中初始裂紋在受拉面與受壓面的裂紋長度相同，而試驗中受拉面裂紋長度遠小于此受拉面的裂紋長度，導致試驗中受拉面應力強度因子計算值受受壓面的影響而降低；受拉面裂紋長度較長時（≥8 mm），應力強度因子變程Δ的計算結果比試驗結果稍高（30%～40%），可能是三維裂紋中的裂尖約束效應導致了這一結果。

圖7 FRANC3D計算和試驗［1］反推的受拉面應力強度因子變程Fig.7 Stress intensity factor range on the tension face calculated by FRANC3D and from experiment［1］

總體而言，應力強度因子隨裂紋長度變化的計算結果比較符合試驗測量的情況，表明FRANC3D-ABAQUS聯(lián)合方法能夠用于面外彎曲載荷下的薄板裂紋擴展的初步分析。

2.3 裂紋前緣形狀特點

裂紋前緣形狀的數(shù)值模擬結果如圖8所示，可以看出：裂紋前緣呈1/4橢圓形狀，與相關文獻中描述的彎曲載荷下裂紋前緣形狀相符。然而，本文中，沿厚度方向的裂紋長度達到約3/4板厚時，裂紋向受壓面擴展極為緩慢，裂紋在整個模擬過程中未擴展到受壓面（受壓面裂紋長度為0 mm），這與試驗觀測的結果不符。在前期的試驗研究中，雖然受壓面裂紋長度始終小于受拉面裂紋長度，受壓面裂紋仍會持續(xù)擴展。

圖8 裂紋前緣變化Fig.8 Propagation of crack front

在FRANC3D-ABAQUS聯(lián)合仿真中，無法擴展的部分會阻礙受拉面裂紋的擴展，從而導致受拉面的應力強度因子計算結果偏小。不僅如此，在受拉面裂紋長度達到一定值后，F(xiàn)RANC3D自動劃分的網格將發(fā)生嚴重畸變，裂紋擴展模擬將終止。

為了探究上述無法擴展的受壓面裂紋對應力強度因子計算結果的影響，本文仍然采用FRANC3D和ABAQUS聯(lián)合仿真方法和2.2節(jié)所述的有限元模型，但不模擬裂紋擴展，而是據(jù)試驗中測量的受拉面與受壓面裂紋長度設置模型中裂紋的長度，并插入1/4橢圓形狀裂紋直接計算應力強度因子，并與之前FRANC3D-ABAQUS聯(lián)合仿真中相同受拉面裂紋長度下的應力強度因子計算結果對比。根據(jù)試驗結果修正受壓面裂紋長度與未修正前的對比結果如表1所示，為試驗中受拉面裂紋長度，為試驗中受壓面裂紋長度，可以看出：FRANC3D-ABAQUS聯(lián)合仿真的應力強度因子計算結果低于修正長度差后的應力強度因子計算結果，表明受壓面裂紋不能擴展這一因素將導致應力強度因子的計算結果偏小10%～15%。

表1 應力強度因子對比Table 1 Stress intensity factor comparison

3 影響疲勞裂紋擴展模擬結果的因素

3.1 邊界條件

薄板結構多采用分布式的鉚接或螺接形式與其他結構連接。在裂紋擴展分析中，常將這種分布式連接的邊界條件簡化為簡單的邊界約束條件。然而，這種邊界條件的約束強度，實際上強于鉸接卻稍弱于完全固支。因此，有必要研究邊界約束強度對應力強度因子計算的影響。

計算三種不同邊界條件（如圖9所示）下的應力強度因子。三種邊界條件下的應力強度因子的計算結果如圖10所示。

圖9 邊界條件示意圖Fig.9 Schematic images of boundary conditions

圖10 不同邊界條件下的應力強度因子Fig.10 Stress intensity factor under different boundary

從圖10可以看出：兩端鉸支下的計算結果明顯高于兩端固支與四端固支，這與彎曲理論計算的兩端鉸支梁最大應力是兩端固支梁最大應力的兩倍的結論相符合；四端固支與兩端固支下的值幾乎相同，表明自由端上的邊界條件對的計算結果影響很小，約束端上的邊界條件對的計算結果影響較大，并且約束端上的約束越強，計算的值越小。約束端上的邊界條件影響明顯是因為約束端的邊界條件對垂直于裂紋方向上的正應力影響較大，而該方向的應力是裂紋擴展的主要驅動力。

3.2 初始裂紋前緣形狀

裂紋萌生的位置和形狀，會影響后續(xù)的裂紋擴展過程。針對直線前緣穿透裂紋、角裂紋和表面裂紋三種不同初始裂紋形式，通過FRANC3DABAQUS聯(lián)合仿真，研究初始裂紋前緣形狀對裂紋擴展的影響。本文所研究的初始直線前緣穿透裂紋長度為6 mm，角裂紋形狀與表面裂紋形狀與尺寸如圖11所示，分別為半徑0.1 mm的四分之一圓與半徑0.1 mm的半圓。

圖11 初始裂紋形狀Fig.11 Initial crack shapes

裂紋擴展后在受拉面（如圖8所示）的應力強度因子計算結果如圖12所示，可以看出：初始角裂紋在剛擴展時應力強度因子較小；在裂紋長度約6.5 mm時，初始角裂紋的應力強度因子超過初始直線前緣裂紋；在裂紋長度大于6.5 mm時，初始角裂紋的應力強度因子隨裂紋長度的變化規(guī)律與初始直線前緣裂紋相似，但數(shù)值略高（約20%）。

圖12 不同初始裂紋的應力強度因子Fig.12 Stress intensity factors under different initial cracks

初始裂紋為角裂紋和表面裂紋的裂紋擴展預測結果如圖13所示。由于表面裂紋在擴展一定長度后，就會與圓孔邊界面連通，從而產生近似角裂紋的情況，只比較連通前的角裂紋與表面裂紋的裂紋擴展情況（受拉面上的裂紋長度與循環(huán)數(shù)的關系）。可以看出：角裂紋在裂紋擴展初期的擴展速率較表面裂紋更快。

圖13 不同初始裂紋下擴展長度與循環(huán)次數(shù)的關系Fig.13 Relationship between propagation length and cycles under different initial cracks

4 仿真方法的適應性

4.1 FRANC3D計算流程的局限性

FRANC3D裂紋擴展模擬方法無法考慮受壓面的裂紋擴展，本文從FRANC3D的計算流程出發(fā)，分析該問題出現(xiàn)的原因。

FRANC3D計算流程為：①根據(jù)有限元計算結果，計算最大應力強度因子；②根據(jù)應力比和計算時設置的其他參數(shù)，按照計算流程（如圖14所示），計算出Δ；③根據(jù)選用的裂紋擴展速率公式，計算一定循環(huán)周次后裂紋前緣上每個節(jié)點的擴展量。

圖14 FRANC3D計算ΔK的流程Fig.14 Process ofΔK calculation in FRANC3D

本文模擬的承受彎曲載荷的模型，受壓部分的應力始終為壓縮應力，故FRANC3D計算的受壓部分的應力強度因子近似為0，應力強度因子變程近似為0，基于裂紋擴展公式計算得到的裂紋擴展長度d也近似為0。而裂紋從受拉面向受壓面的擴展，主要是靠裂紋前緣擴展的節(jié)點“拉動”裂紋前緣的擬合曲線，從而“帶動”部分未擴展的節(jié)點，如圖15所示。而在靠近受壓面處，裂紋擴展長度始終接近0，因此擴展的節(jié)點也無法“帶動”未擴展的節(jié)點，導致受壓面的裂紋在仿真中始終無法擴展。

圖15 裂紋前緣節(jié)點擴展Fig.15 Propagation of nodes on crack front

本文嘗試在FRANC3D的框架下，通過調整部分設置，實現(xiàn)受壓面上的裂紋擴展。首先，通過調整應力強度因子計算方法，改善受壓面上應力強度因子計算結果為0的情況。然而，F(xiàn)RANC3D提供的-積分、虛擬裂紋閉合技術、裂紋張開位移三種計算方法，計算結果相似，無論做出何種選擇，受壓面的應力強度因子計算結果始終約為0；其次，更改到Δ計算流程中的選項對應力強度因子變程計算結果的影響。然而，任何調整均無法從接近0的得到能使受壓面裂紋擴展的Δ；最后，通過調整裂紋擴展速率公式使得裂紋面擴展的嘗試也以失敗告終，原因在于，基于Δ的裂紋擴展速率公式無法在Δ=0的情況下，計算出裂紋擴展速率。

綜上所述，在計算受彎曲載荷下的薄板的裂紋擴展時，由于FRANC3D基于值的計算框架，必將導致總受壓縮應力的那一面的裂紋無法在模擬中擴展。

4.2 仿真方法適用的裂紋范圍

為了進一步研究FRANC3D方法的適用性，計算面外彎曲載荷下其他類型裂紋的應力強度因子，包括表面裂紋和直線前緣穿透裂紋，這些裂紋也有可能在飛機結構中產生并受面外彎曲載荷作用。

對于兩端受彎曲載荷作用下的有限板（=150 mm，=130 mm，=3 mm，如圖16所示），采用FRANC3D-ABAQUS聯(lián)合方法，計算中心表面半橢圓裂紋擴展過程中的應力強度因子的變化，計算結果如表2所示，表2中同時給出應力強度因子擬合公式的計算結果。

圖16 有限板中的表面裂紋［5］Fig.16 Surface crack in finite plate［5］

表2 有限板表面裂紋的應力強度因子Table 2 Stress intensity factor of surface crack in finite plate

從表2可以看出：擬合公式的計算結果與FRANC3D-ABAQUS聯(lián)合方法僅在=2.5 mm、=π/2處，差異較大；其他情況下，差異很小。

對于受彎曲載荷的無限大板（如圖17所示），采用FRANC3D-ABAQUS聯(lián)合方法，計算中心直線穿透裂紋在裂紋擴展過程中應力強度因子的變化，計算結果如表3所示，表3中同時給出文獻［3］不考慮裂紋面接觸時的理論應力強度因子計算結果。

圖17 無限大板中心的直線穿透裂紋［3］Fig.17 Linear penetrating crack in the center of infinite plate［3］

表3 無限大板直線穿透裂紋的應力強度因子Table 3 Stress intensity factor of linear penetrating crack in infinite plate

從表3可以看出：兩方法的計算結果在裂紋較短（＜10 mm）時差異較小，在裂紋擴展超過10 mm后，差異較大。

綜上所述，針對受彎曲載荷的薄板，F(xiàn)RANC3 D-ABAQUS聯(lián)合方法能夠很好地計算受拉面表面裂紋或穿透裂紋擴展中的應力強度因子，而在分析裂紋前緣形狀復雜且板厚度較大的情況時，計算結果不理想。

5 結論

（1）對于面外彎曲載荷下的薄板，在采用FRANC3D-ABAQUS聯(lián)合方法進行疲勞裂紋分析時，施加較強的邊界約束會導致應力強度因子的計算結果增大。

（2）在承受面外彎曲載荷的薄板中，初始非孔邊表面裂紋比初始孔邊角裂紋在裂紋擴展初期的擴展速率低。

（3）FRANC3D-ABAQUS聯(lián)合方法在分析面外彎曲載荷作用下薄板孔邊裂紋擴展時，存在局限性。現(xiàn)有方法中僅受拉面的裂紋擴展，由于FRANC3D方法對受壓面的裂紋應力強度因子變程計算結果近似為0，受壓面裂紋無法擴展。