基于斷裂力學的船型網箱疲勞損傷評估研究

劉安民，謝凌俊，蘇一鳴，丁洪濤，何文濤

（中國海洋大學 工程學院，山東 青島 266000）

利用簡化計算方法進行海洋結構物疲勞壽命評估時，波浪載荷與結構響應都是通過經驗公式計算的。這一方面會限制簡化計算方法在缺乏數據樣本的新型結構物上的應用，另一方面也會帶來計算精度上的問題[1-3]。直接計算方法是公認較為準確的海洋結構物疲勞強度計算方法，與簡化計算方法相比，直接計算方法中結構的疲勞載荷并不是通過經驗公式計算的，而是通過波浪載荷軟件，直接運用勢流理論和莫里森方程對有限元模型進行計算得到的，并且通過有限元結構分析得到了結構的應力。這不僅能更好地反映全船疲勞載荷分布，也能夠更好地反映出結構細節之處的載荷特點[4]。

譜分析方法是直接計算方法的一個重要方法，譜分析方法物理意義明確，過程合理，考慮因素全[5]，被認為是一種準確的方法，也得到了較為廣泛的研究。賴明雁[6]認為，譜分析方法是建立在真實的海況與真實裝載工況下的直接計算方法，計算精度高。徐帥[7]利用譜分析方法開展了LNG 船疲勞可靠性的分析，對譜分析方法的計算精度表示了肯定。宋憲倉[8]展開了半潛平臺關鍵節點疲勞壽命評估的研究，并對結果進行了寬帶修正，研究表明，相比于確定性方法，譜分析方法能夠考慮波浪能量的分布，結果更為精確。曾琪[9]利用挪威船級社開發的海洋結構分析軟件SESAM 對江海直達船進行了基于S-N曲線法的疲勞譜分析研究，并且將結果與簡化計算方法進行了比較，認為譜分析相比于簡化計算方法能得到更加準確的結果。白曉東[10]開展了深海網箱結構時域方法和頻域方法的疲勞分析，并探討了包括Rayleigh 分布、Gamma 分布和GEV 分布在內的多種分布模型對疲勞壽命結果的影響，表明對于非高斯過程，利用Rayleigh 分布對應力范圍概率密度進行描述可能會產生較大的誤差。由于譜分析方法是建立在應力響應窄帶高斯過程的基礎上的，因此不少學者對非高斯過程[11]、雙峰或者三峰高斯過程[12]、高低頻載荷共同作用[13]、窄帶修正[14]等方面也做了不少工作，以提高譜分析方法的精度。

為了加深對基于斷裂力學方法的海洋結構物疲勞裂紋擴展分析與壽命預測研究的認識和了解，形成更加成熟完善的研究體系，學者們從不同層面展開了廣泛的研究。Tanaka 等[15-17]對管節點結構開展了疲勞裂紋擴展分析研究，豐富了海洋工程中常見節點結構的疲勞性能研究。王小松[18]以在役FPSO 為對象，根據英國標準BS 7910 進行了基于斷裂力學的疲勞研究，并進行了多種規范之間的比較，研究結果對于營運安全具有一定的指導意義。

本文所采用的海洋結構物疲勞壽命評估方法將譜分析方法和斷裂力學方法相結合，利用有限元計算的方式得到結構物在真實海況下的載荷條件，從疲勞載荷的角度保證疲勞評估的精度[19]，然后通過斷裂力學進行疲勞評估。本文將以船型網箱為例，進行海洋結構物的運動特性以及結構應力分析，并利用子模型技術進行局部應力分析。之后，運用多尺度子模型實現殼單元到實體單元的轉變，并在實體單元部分進行裂紋擴展，得到疲勞壽命。

1 疲勞分析方法和有限元模型

1.1 譜分析方法

譜分析方法的主要分析流程：

1）對海洋結構物進行水動力分析與應力響應分析，并提取疲勞熱點的應力傳遞函數Hstress。

2）基于海洋結構物為線性系統的假設，根據海況i的波浪譜，求得應力響應譜密度函數。

3）基于應力響應是窄帶高斯過程的假設，求得應力范圍的概率密度函數f(ΔS)。

4）根據線性累積損傷理論，可以得到海況i下每年的疲勞損傷度，見式（1）。

式中：pi為海況概率；f0i為海況i對應的平均跨零率；ni為循環次數；Ni為循環破壞的壽命；S表示應力范圍；A和m為疲勞實驗參數。

5）循環步驟1）、2）、3），在對所有短期海況進行遍歷后即可得到1 a 內的總損傷度，見式（2）。

1.2 裂紋擴展增量計算

對于裂紋穩定擴展階段，最熟悉的擴展模型莫過于Paris 模型[20]：

式中：C和m為材料參數；為裂紋擴展速率；ΔK為應力強度因子范圍；ΔKth為應力強度因子范圍門檻值。當應力強度因子范圍ΔK小于門檻值ΔKth時，裂紋將不能夠擴展。

單一短期海況下的裂紋擴展增量E(da)可以根據式（4）–（6）進行計算：

式中：ΔNi為第i個短期海況下的循環次數；ΔTi為對應短期海況下的持續時間；fi為對應海況下的平均過零率；fiΔK為應力強度因子范圍的概率密度函數；m0和m2為應力強度因子響應譜0 階和2 階譜矩。

由此可得，單一短期海況下的裂紋擴展增量為：

1.3 有限元模型

1.3.1 模型參數和坐標系

船型網箱主要由船艏艙室、船艉浮箱以及船艏船艉之間的箱型梁組成。船型網箱主要參數見表1。

表1 船型網箱主要參數Tab.1 Main parameters of ship shaped cages

坐標系原點在船底中線上，x軸正方向由船艏指向船艉，y軸由左舷指向右舷，z軸由船底指向上，如圖1 所示。同時，波浪入射角θ的定義方式也在圖1 中進行了展示，波浪入射角θ即波浪入射方向與x軸正方向之間的夾角。

圖1 網箱參考坐標系定義Fig.1 Definition of the cage reference coordinate system

1.3.2 網箱有限元模型

由于海洋結構物主尺度較大，結構復雜，同時結構的厚度效應也會在水動力分析中有明顯的體現，因此對于海洋結構物，一般利用殼單元以及梁單元進行有限元建模。網格主要為四邊形單元，在結構不連續、曲率變化大處等將利用三角形單元進行離散。對于網箱上大量的T 形材、角鋼等細長結構將利用梁單元進行建模，如圖2a 所示。

圖2 船型網箱有限元模型Fig.2 Finite element model of ship shaped cage:a) element type;b) constraint settings

結構應保持自由動態平衡狀態，需要對其自由度進行約束[21]。DNV 等船級社做出了相應的規范要求[22]，采用三點約束法對網箱進行約束，如圖2b 所示。約束點1 位于船艏底板中線縱艙壁下，約束其x、y、z三個自由度上的平動位移；約束點2 位于船艉右浮箱橫縱型材的交點下，約束其y、z兩個自由度上的平動位移；約束點3 位于船艉左浮箱橫縱型材的交點下，約束其z軸方向上的自由度。

1.3.3 服役環境與載荷

海洋結構物在服役期間，同時遭受風浪流等載荷的共同作用，由于風和流的作用較弱，因此本文只考慮波浪載荷的作用。本文船型網箱服役于北大西洋海域，采用北大西洋海域波浪散布圖。進行頻域內分析時，由于結構的對稱性，單位波幅規則波入射方向設為0°～180°，具體為0°、30°、60°、90°、120°、150°、180°。同時，波浪頻率設置為 0.2～2.4 rad/s，間隔0.1 rad/s。

1.4 子模型

針對海洋結構物這類尺寸能夠達到百米，并且結構復雜、細節豐富的超大型結構，在建模時保留全部細節將會導致建模工作量的顯著增加，同時也會大大降低計算效率，還可能帶來應力奇異現象等不利影響。因此，一般利用簡化模型、粗網格模型對結構進行整體建模與分析，而采用子模型技術進行局部區域的建模與應力分析。

子模型技術的主要分析流程如圖3 所示，主要步驟如下[23-24]：

圖3 SESAM 子模型Fig.3 SESAM submodel

1）全局模型的建模與分析。對全局模型進行簡化建模，采用較粗的網格進行單元劃分，然后進行全局模型的分析。根據全局模型的分析結果，確定需要進行重點分析的局部區域（子模型區域）。

2）子模型的建模。根據實際尺寸、結構以及分析目標等，對子模型進行精細化建模，添加整體模型中未考慮的細節結構，并利用精細網格進行單元劃分。

3）子模型邊界繼承。根據子模型的切割邊界，從全局模型的分析結果中讀取計算當前子模型的邊界條件。

4）子模型分析。利用子模型進行局部結構的應力分析等，得到局部結構更精確的應力分析結果等。

在對大型結構進行有限元分析時，為了兼顧計算效率與計算精度，通常不會采用單一尺度（如殼單元、實體單元、梁單元）的單元進行結構離散，而是利用多尺度子模型技術建立多尺度耦合模型，在保證計算精度的同時提高計算速率，殼-實體耦合模型便是其中的重要形式[25]。在本文海洋結構物裂紋擴展分析中，殼-實體耦合模型的采用主要基于以下2 點考慮：1）載荷計算是在SESAM 中以殼、梁單元的形式完成的，因此殼單元的保留有利于更加靈活便捷地實現載荷的傳遞；2）實體單元的采用主要為了體現結構細節，得到更加精準的局部應力，實現裂紋在實體中的擴展。

1.5 多尺度耦合模型

殼-實體耦合模型中，一部分單元為殼單元，一部分單元為實體單元。由于殼單元節點具有六自由度，而實體單元節點僅僅具有3 個平動自由度，因此殼單元與實體單元連接處將會出現自由度不一致的情況，需要進行內部分布耦合約束以實現節點之間的協調[26]。

殼-實體耦合在 ABAQUS 中的實現（通過shell-to-solid coupling 命令）如圖4 所示。為了減少后續有限元計算的時間，殼-實體耦合模型的范圍只是SESAM 中子模型的局部，相當于從子模型中再次提取最為關心的疲勞熱點部分作為子模型，這樣處理的好處是既能更加突出重點區域的應力分布，也能縮小模型，提高計算效率。從SESAM 中子模型導出分析區域（二級子模型），并將其中實體區域內結構以實體單元重新建模，最后通過shell-to-solid coupling實現殼體單元與實體單元之間的連接，形成殼-實體耦合子模型。

圖4 殼-實體耦合子模型Fig.4 Shell-solid coupling submodel

此外，殼-實體耦合子模型的載荷是通過節點對節點的形式從上一級子模型中繼承的，并賦予到圖4中殼單元最外層節點上。這部分工作主要是將從SESAM 軟件中導出的邊界上所有節點的六自由度運動逐一導入到ABAQUS 殼-實體耦合子模型中對應的節點上，從而實現載荷的傳遞，如圖5 所示。

圖5 載荷傳遞示意圖Fig.5 Schematic diagram of load transfer

2 應力分析結果及疲勞分析

2.1 運動特性分析

在SESAM 軟件中，根據三維勢流理論，可以求得波浪場的速度勢，并可以進一步求解出單位波幅規則波作用下結構的六自由度運動響應。本文船型網箱在0°、30°、60°、90°、120°、150°、180°入射方向下的六自由度（垂蕩、縱蕩、橫蕩、橫搖、縱搖、艏搖）運動響應傳遞函數如圖6 所示。

圖6 運動響應傳遞函數RAOFig.6 Motion response transfer function RAO:a) heave;b) roll;c) surge;d) pitch;e) sway;f) yaw

由圖6 可見，對于船型網箱，六自由度運動中垂蕩運動是其中最激烈的運動模式。隨著波浪頻率逐漸增大，垂蕩運動幅值先增大、后快速減小到0 附近。這可能是由于在波浪頻率等于1.1 rad/s 附近時，網箱出現共振現象，導致垂蕩運動的迅速增大。共振區的左側區域，頻率較小，波周期較大，結構所承受的浮力發生變化，出現了明顯的回復力。對于縱蕩以及橫蕩運動，基本都呈現出隨著波浪頻率的減小、波浪周期的增大，運動幅值不斷增大的趨勢，這主要是由于水動力分析時未考慮系泊力的影響，無法提供回復力以抵消結構運動導致的。同時可以看出，對于縱蕩運動，波浪角度為0°以及180°時，運動幅值最大。因為此時波浪方向與船型網箱的長度方向一致。橫蕩的運動幅值則在90°時最大，因為此時波浪方向正好與船型網箱長度方向垂直，承受著最大的橫向波浪力。

橫搖與縱搖的運動幅值隨著波浪頻率變化趨勢基本一致，均是先增大到一定程度后逐漸減小。這是因為橫搖與縱搖均會受到回復力的影響，在超過共振區之后，運動幅值受到回復力的作用逐漸減小。同時，縱搖與縱蕩一致，運動幅值在0°及180°時最大，而在90°時最小。橫搖則相反，運動幅值在90°時最大，而在0°和180°時最小。艏搖運動是六自由度運動最為微弱的運動，并且運動幅值小，基本不發生運動。

2.2 結構應力分析

網箱在服役過程中，始終遭受著復雜的靜水壓力與波動載荷的作用，并最終反映到結構應力上。對結構進行應力分析有利于了解網箱整體的受力情況，判斷結構整體強度，評估最危險工況與結構上最危險的節點等。同時，局部節點疲勞強度分析等工作的進行也需要應力分析的支持，特別是局部應力分析。因此，本節對網箱模型在波浪載荷作用下的應力進行分析研究。

2.2.1 整體應力分析

網箱在靜水壓力與單位波幅波浪力共同作用下的整體應力分析結果如圖7 所示。由于結構的對稱性，僅對0°、30°、60°、90°波浪入射方向進行展示。從圖7 中可以看出，網箱整體上應力分布并不均勻，呈現出縱向箱型梁與船艏底部應力高，而艉浮箱與橫向箱型梁的應力偏低的趨勢。同時，盡管網箱整體應力并不高，但是在各波浪入射角下，均出現了明顯的局部節點應力集中現象。如船中角隅處、船艏與兩側箱型梁連接處、上下箱型梁與之間斜撐的連接點、橫縱箱型梁連接點等局部區域都是應力集中明顯的區域，這些區域都是潛在的裂紋源或者疲勞損傷嚴重的區域。

圖7 網箱整體應力響應Fig.7 Overall stress response of cages

2.2.2 局部應力分析

利用子模型技術得到的網箱局部應力分析結果如圖8 所示，其中節點D1位于船艏浮箱與橫箱型梁相貫處，節點D2和D3位于船中角隅處。從圖8 中可以看出，通過子模型技術得到的局部應力與通過整體分析得到的相同位置處的節點應力在應力分布與應力值上大體一致，這說明子模型技術能夠比較真實完整地繼承整體模型的載荷條件。同時，由于子模型技術采用了更加精細的網格進行單元劃分，在各拐角處的應力集中現象均更加顯著。由此可見，利用子模型技術進行網箱上局部結構的應力分析能夠更好地反映局部結構的真實應力分布，同時能夠更好地突出應力集中現象，對疲勞熱點的疲勞強度分析、疲勞壽命預測具有重要的意義。

圖8 網箱局部應力響應Fig.8 Local stress response of cages

殼-實體耦合子模型與殼模型在相同載荷條件下的位移云圖對比如圖9 所示?？梢钥闯觯瑲?實體耦合模型的計算結果與殼模型計算結果的一致性較好，證明了利用多尺度子模型技術建立殼-實體耦合子模型的可行性。

圖9 多尺度子模型計算效果Fig.9 Result of multi-scale submodel:a) shell-solid coupling model;b) shell model

2.3 裂紋擴展和疲勞壽命

表面裂紋一般用半橢圓描述，因此本文也利用雙參數半橢圓裂紋進行裂紋的描述。裂紋長度為2a，裂紋深度為c，裂紋長深比為a/c?；诖耍谄诹鸭y擴展分析中，僅需要針對裂紋深度與裂紋長度進行應力強度因子傳遞函數與裂紋增量的計算。本文將初始裂紋定義為a=c=0.5 mm 的圓形裂紋。

基于斷裂力學的疲勞裂紋擴展主要流程：

1）初始化裂紋半長、深度、循環次數與時間。

2）在殼-實體耦合模型中，逐次施加某一浪向下所有不同頻率的單位波幅規則波對應的邊界條件，插入當前尺寸裂紋，并計算對應的SIF，以此構建該裂紋尺寸、該浪向下的應力強度因子傳遞函數。

3）根據應力強度因子傳遞函數，結合各短期海況的波浪譜，進行各短期海況下裂紋擴展增量的計算，并進行疊加。

4）當該浪向下所有短期海況均已循環完畢或者達到裂紋更新閾值時，進行裂紋尺寸更新，并重復步驟2）—4）。

D1、D2、D3這3 個節點在裂紋擴展過程中裂紋半長及裂紋深度隨時間變化的曲線如圖10 所示?？梢钥闯觯? 個節點的裂紋擴展速率存在著明顯的差別，節點D1的擴展速率最高，其次分別為節點D2、D3。對應地，按照裂紋深度達到80%板厚時認為結構達到疲勞失效的判定標準，節點D1的疲勞壽命最短，而節點D3的疲勞壽命最長，三者疲勞壽命分別為11.3、17.3、33.1 a。

3 結論

本文以船型網箱為例，進行了海洋結構物頻域內的響應分析，介紹了網箱的主要參數、參考坐標系和有限元模型建模與邊界條件，同時開展了結構的運動特性與應力分析，并對疲勞熱點開展疲勞分析?？傻玫揭韵轮饕Y論：

1）通過水動力分析得到了該船體網箱的運動特性，為疲勞壽命評估提供了更準確的波浪載荷與結構響應。六自由度的運動幅值均與浪向和波頻有著明顯的關聯，對于有明顯共振頻率的運動，其幅值往往在共振頻率附近達到最大值，而對于沒有共振頻率的運動，其幅值往往隨著波頻的減小而增大。垂蕩是六自由度運動中最激烈的運動模式，在0°與180°波浪入射角、1.1 rad/s 的頻率下運動響應最激烈?？v搖和縱蕩的運動幅值也在0°以及180°時最大。

2）網箱整體應力符合強度要求，但在不同浪向下均可觀察到明顯的應力集中現象，應力集中區域主要為船中角隅處、船艏與兩側箱型梁連接處、上下箱型梁與之間斜撐的連接點、橫縱箱型梁連接點等。

3）子模型技術能夠完整準確地繼承局部區域的載荷條件并展開應力分析，能夠更準確地反映局部區域的應力集中現象，為疲勞強度分析提供更準確的載荷條件。

4）實現了基于子模型技術的海洋結構物強度分析與壽命預測。運用子模型進行局部應力計算，并將多尺度子模型技術應用到其中，解決了大尺度模型與小尺寸裂紋之間的耦合、模型與載荷的傳遞、計算效率與計算精度的平衡等問題。