復雜載荷下裂紋SIF計算的子模型技術

2021-04-12 10:01:54黃小平

艦船科學技術 2021年3期

江威，黃小平

(上海交通大學海洋工程國家重點實驗室，上海 200240)

0 引言

疲勞強度校核一直是船舶結構強度中的重要組成部分，疲勞強度評估中，基于斷裂力學方法的壽命預報能夠考慮初始缺陷、載荷次序等因素的影響，是未來發展的趨勢。但將基于斷裂力學的壽命預報方法應用于工程實際結構時，對于重要參數SIF的求解仍存在問題[1]。目前較權威的SIF計算公式[2– 5]僅限于簡單載荷和簡單應力場，尚不能對工程實際結構進行推廣，DNV規范推薦將子模型技術應用于復雜載荷下構件中裂紋SIF求解[6]。子模型技術是在整體分析模型基礎上獲取局部區域精確解的一種高級有限元分析技術，其中子模型是為分析結構細節（孔、角隅、裂紋、接觸面等）而從原模型分割出的一部分區域[7]。通過求解原模型，子模型邊界條件一般可由位移法（DBS）或應力法（SBS）得到[8]，其中位移法因其離散誤差收斂更快被廣泛推薦使用[9–11]，故子模型技術又稱為“切割邊界位移法”或“特定邊界位移法”[12]。子模型技術基于圣維南原理，即實際分布載荷如果被等效載荷代替，應力和應變只會在施加載荷的位置發生變化[13–14]，這表明應力集中效應只會發生在載荷施加部位。

基于Ansys內部的子模型分析技術具有以下優點[15]：

1）減少甚至取消了有限元模型中的復雜傳遞區域；

2）方便對局部區域進行細節分析并得到精確解；

3）可用于驗證原模型的網格劃分是否滿足要求。

子模型分析過程一般包括建立和分析原整體模型、建立子模型、生成切割邊界位移及分析子模型4個步驟，并有以下注意事項[7,16]：

1）子模型和原整體模型相應位置處的單元屬性及位置需一致；

2）切割邊界只能選在殼單元、體單元內或殼體單元連接處；

3）切割邊界應遠離應力集中區域；

4）原整體模型相應位置處網格有必要進行足夠程度的細化。

基于船舶結構的特性及Patran對殼、梁單元的強大支持，行業內一般在Patran中建立全船有限元模型，此時若對相關疲勞熱點進行斷裂力學評估，需考慮結合Patran和Ansys兩套軟件的子模型技術，即通過在Ansys中建立與Patran原整體模型待評估位置對應的子模型進行分析，與單純在Ansys內部實現邊界切割和位移插值的子模型技術相比，除需滿足基本操作要求，還有以下問題：

1）疲勞壽命評估中，對不同尺寸表面裂紋SIF求解需要使用三維單元[17–18]，也就意味著Ansys子模型及原Patran殼模型的對應位置都要使用體單元建模，即刪除原Patran整體模型對應位置處的殼單元，用體單元重建，并用MPC連接重建的體單元和其周圍殼單元，方可保證殼體間變形正確傳遞；

2）對于結合Patran和Ansys的子模型技術的疲勞評估，若疲勞熱點較多且位置各異，僅在Ansys APDL中往復調整子模型位置與原模型一致是不明智的，此時建立子模型應優先考慮方便裂紋建立的做法，然后在Ansys外部實現坐標轉換。由此可見，這種子模型技術的前處理的實現難度較高，各環節的實現低效易出錯。目前較少學者采用子模型技術進行裂紋擴展分析[7]，也很少有人對有限元前處理過程進行有效優化[19]，這在很大程度上限制了子模型技術在船舶疲勞強度校核中的推廣。

本文提出“逐周分層法”、“轉換矩陣法”及“映射劃分法”，著眼于解決從Patran整體模型到Ansys子模型實現過程中殼體單元間MPC創建低效、不同坐標系間節點位置信息轉換困難及裂紋自由劃分局限性的問題，并基于VBA和APDL語言編寫插件M PC_arranger_V 1.0，FEM_coo r_transfer_V 1.0及C rack_m apper_V 1.0進行實現。同時選取某B型LNG燃料艙關鍵疲勞節點，參考DNV及ABS相關規范[20]，對子模型邊界位移進行驗證，并進行疲勞熱點處裂紋的SIF計算。

1 子模型技術簡介

基于斷裂力學進行疲勞評估需要求解SIF，目前Ansys自帶SIF求解命令KCALC，如圖1所示。

圖 1 應力強度因子計算流程Fig. 1 The procedure of SIF calculation

針對一種結合Patran和Ansys的子模型技術[1]，以承受簡單拉伸載荷的平板為例介紹其實現流程：

1）在Patran中采用shell單元建立厚度為t的平板模型，對板左端施加六自由度約束，對板右端施加垂直于板厚方向的拉應力。對板中心區域網格進行t×t細化，細化區域不少于10t范圍，為實現對三維表面裂紋SIF的求解，需將細化區域的板單元拉伸為實體單元，并在殼體邊界使用RSSCON Surf-Vol型MPC關聯對應節點自由度，然后提交計算。

2）在Ansys中建立與Patran體單元部分對應的子模型，并對裂紋進行網格劃分。然后輸出Patran中實體單元邊界節點的位置信息，經坐標轉換及插值施加為Ansys子模型的邊界條件，即可實現SIF的求解。

2 殼體邊界MPC的創建

如圖2所示，Patran模型中體單元和殼單元的邊界節點需要使用RSSCON Surf-Vol型MPC進行關聯。此MPC需將上下層邊界節點（即體單元邊界節點）設為獨立點，中間層邊界點（即殼體聯合邊界節點）設為非獨立點。

圖2 邊界MPC的創建Fig.2 Creation of boundary MPC

為解決MPC創建工作量大，效率低同時易出錯的問題，本文提出“逐周分層”思想，分別以平板模型和部分典型曲面模型為例解釋，并基于VBA語言編寫插件MPC_arranger_V1.0進行高效實現，另外，本文所提體單元皆為8節點6面體solid單元。

2.1 平面逐周分層法

在船體結構的有限元模型中，以平板單元最為常見。如圖3所示，以某簡單平板為例，先對其中央區域進行網格細化，然后拖拽板單元建立實際厚度為t的殼體混合模型。以局部區域中心為原點O建立柱坐標系，以下皆以殼單元邊界上任一待建立MPC的節點P為例，對逐周分層法做出解釋。其中，O點和P點在原Patran總體坐標系中的位置分別為(xO,yO,zO)和(xP,yP,zP)，可視作已知。

柱坐標系中P點相對O點的位置 (ρ,θ,z)，可通過下列方程求得：

圖3 平面逐周分層法Fig.3 Circum ferentially-layering in plane

通過Patran輸出切割邊界節點的位置信息，并使用插件進行逐周分層處理，即可快速實現各層邊界點的重新排序，并匯總輸出準確對應的各層節點號。相對GUI方式，視邊界點節數量通?？蓪PC創建效率至少提高百倍。

2.2 曲面逐周分層法

除常見的平板結構，船舶也存在諸多典型曲面結構物，本節針對部分應用廣泛的典型曲面結構使用逐周分層法。

大多數潛艇艇體[21]及海洋平臺結構[1,22]一般采用圓柱殼模型，圖4為實際厚度為t的圓柱殼體混合模型。在圓柱坐標系中，外層節點所在圓柱面，中層節點所在圓柱面及內層節點所在圓柱面分別對應底部3個同心圓面，其中中層圓柱面的底圓面半徑為r，將底圓面圓心O取為圓柱坐標系原點，以上參數視作已知。P點的圓柱系坐標 (ρ,θ,z)同式（1）。

圖4 圓柱面逐周分層法Fig.4 Circum ferentially-layering on cylinder

除圓柱殼模型應用廣泛，橢圓柱殼由于其在裝備布置和后屈曲性能方面的優勢，也逐漸被應用于潛艇設計[23]，如941型戰略核潛艇。圖5為實際厚度為t的橢柱形殼體混合結構。在橢柱坐標系中，外層橢柱面，中層橢柱面及內層橢柱面分別對應3個底橢圓面，其焦點分別為A1，B1，A2，B2和A3，B3，其中中層橢柱面的底橢圓長軸為a，將底橢圓面中心O取為橢柱坐標系原點，以上參數視作已知。

其中，P點的橢柱系坐標 (ρA2+ρB2,θ,z) 的分量θ和z見式（1）， ρA2為P點與焦點A2點間的距離, ρB2為P點與焦點B2點間的距離，可通過下列方程求得：

圖 5橢柱面逐周分層法Fig.5 Circum ferentially-layering on elliptic cylinder

目前眾多深水潛器[24]、潛艇艇首[25]、MOSS型LNG艙[26]都采用圓球殼模型，圖6為實際厚度為t的圓球殼體混合模型。在圓球坐標系中，中層圓球面半徑為r，原點O取在圓球中心，以上參數視作已知。

其中，P點的圓球系坐標 (R,θ,?) 中分量θ表達式見式（1），R和 ?可通過下列方程求得：

蘇聯部分Y級潛艇艇首采用橢球殼結構，圖7為實際厚度為t的橢球殼體混合模型。在橢球坐標系中，外層橢圓截面，中層橢圓截面及內層橢圓截面的焦點分別為A1，B1，A2，B2及A3，B3，中層橢圓截面長軸為a，原點O取為橢球中心，以上參數視為已知。

圖6 圓球面逐周分層法Fig.6 Circum ferentially-layering on sphere

其中，P點的橢球系坐標 (RA2+RB2,θ,?) 中分量θ和?參見式（1）和式（3），RA2及RB2可由以下方程求得：

2.3 注意事項

針對本文所提逐周分層法，尚有幾點說明：

1）逐周分層是一種方便MPC創建的思想。簡言之，即先按邊界點周向位置進行排序，再按層向位置歸類，反之亦可，本文暫時采用先逐周后分層的手段。此外，在進行周向排序和層向歸類時，判據并不唯一，本文僅取其中易于實現的一類判據進行實現。

圖 7橢球面逐周分層法Fig.7 Circum ferentially-layering on ellipsoid

2）使用逐周分層法具有一定的前提條件。各船級社疲勞評估規范均要求將局部區域的網格細化的較為規整，一是方便熱點應力的插值，二是網格過渡和單元拖拽都比較簡單，三是方便Ansys中子模型的建立。所以本文所提逐周分層法暫時考慮較為規整的邊界。簡言之，若整體模型對應位置處網格細化較為合理，這對于大多數復雜結構并不難實現，則推薦使用本文方法，若前期網格細化過于隨意，會使切割邊界及子模型的建立變得復雜，也會影響本方法的使用。

3）在實際操作中，各種誤差均需考慮，如曲面網格光順程度的影響，可直接在插件中進行設置。本文僅在理論層面進行說明，暫不涉及誤差分析。

3 基于矩陣法的坐標變換

在Ansys中求解SIF，需要對體單元節點在原Patran坐標系中的位置進行坐標轉換，并經插值才可施加為Ansys子模型的邊界條件。閆小順[1]自編patran2ansys宏插件完成此部分工作，但插件中需要輸入不同坐標系間的相對平動及轉角位移，才能實現坐標轉換。實際上，這僅可推廣于部分較為簡單的平動坐標變換，對于轉動坐標變換，存在以下問題：

1）兩坐標系各平面間的相對轉角其實并不直觀，即便先通過計算方向向量來求得各向轉角，也費諸多周折。

2）體單元中或存在成千上萬的節點，若每次求解SIF均通過Ansys APDL處理坐標變換，在批量求解各尺寸裂紋的SIF時，也會增加前處理的負擔。

為解決以上問題，本文保留patran2ansys插值部分功能，并使用VBA語言編寫插件FEM_coor_transfer_V 1.0，基于轉換矩陣法對其坐標變換部分進行改進，原理如下：

為保證轉換矩陣可逆，式中(p1,q1,r1),(p2,q2,r2),(p3,q3,r3)需取為Patran體模型中任意不共線的3個節點坐標，這里推薦選取易于識別的體單元邊界點。(a1,b1,c1),(a2,b2,c2),(a3,b3,c3)為Ansys子模型中與上述3點對應的節點坐標，這種對應是指幾何形狀的對應。 (xp,yp,zp)為 Patran中某待轉換節點坐標；(xa,ya,za)為Ansys中待求節點坐標； (xO,yO,zO)為Patran坐標系原點坐標。簡言之，本文所提矩陣變換法更關注模型的宏觀幾何形狀，并不考慮坐標系間的相對位置關系。根據公式（5）及Patran中可輸出的原節點坐標，經FEM_coor_transfer_V1.0處理后可直接輸出指定格式的節點坐標信息，方便實現子模型邊界信息交互。

4 基于映射網格的裂紋建模

求解復雜載荷下構件中的裂紋擴展問題需要計算不同尺寸裂紋的SIF，若通過部分裂紋的SIF推得所有可能的SIF，一般采用擬合公式或插值的方法。擬合公式工作量較大，只有對常見的結構，如縱骨端部[27]，才有必要進行擬合，對于不常見或沒有經驗公式可用的結構，插值也是一種推薦的手段。但插值中存在外插準確性不高的問題，故應盡可能地采用內插值方式，前提是必須最大限度地獲取所有可能出現的裂紋的SIF。

5 子模型技術的驗證

如圖8所示，本節以某B型LNG燃料艙艙頂橫向支座處弧形肘板趾端處的疲勞熱點為例，按照規范對熱點進行不少于10t范圍的t×t細化，并對熱點所在底板在板厚方向設置4層體單元，采用改進的子模型技術及所編插件，完成熱點處各殼體單元間MPC的創建、節點坐標轉換及裂紋SIF的求解。目標艙段的疲勞載荷根據ABS規范中的疲勞工況施加，包括迎浪、橫浪及斜浪工況。

圖8 疲勞熱點及子模型Fig.8 Fatigue hotspot and sub-model

對復雜載荷或復雜結構使用子模型技術，參考DNV規范，需對邊界條件進行驗證。如圖9所示，本文選取對橫向肘板趾端影響較大的橫浪工況進行驗證，經對比，可見Patran中體單元和Ansys子模型的邊界及整體位移吻合度較高，可證明本文所提子模型技術的合理性。

上述疲勞熱點處構件尺寸及裂紋位置如圖10所示，其中T1=30mm，L=25mm，B=9.5mm，T2=12mm。

裂紋SIF大小受復雜載荷成分及焊趾的影響。此處焊縫2L較長且板厚B較小，參考ABS規范，并考慮裂紋實際擴展情況，a取0.5～9 mm，a/c取0.15和0.2，I型SIF計算結果如圖11所示。

6 結語

本文針對使用子模型技術求解復雜載荷下裂紋SIF時各環節存在的問題，結合Patran和Ansys兩款有限元軟件，基于VBA及APDL語言，進行如下改進：

1）針對殼體單元間MPC的創建低效問題，提出更加高效的逐周分層法并編寫插件M PC_arranger_V1.0進行實現。給出對平板及各典型曲面結構進行逐周分層的相關公式及邏輯流程，且成功應用于艙體疲勞熱點處MPC的快速建立。相對GUI操作，一般視邊界點數量可將MPC創建效率提高數百倍，具有一定推廣意義。