江 威，黃小平

(上海交通大學(xué) 海洋工程國家重點(diǎn)實(shí)驗(yàn)室，上海 200240)

0 引 言

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

基于Ansys內(nèi)部的子模型分析技術(shù)具有以下優(yōu)點(diǎn)[15]：

1）減少甚至取消了有限元模型中的復(fù)雜傳遞區(qū)域；

2）方便對局部區(qū)域進(jìn)行細(xì)節(jié)分析并得到精確解；

3）可用于驗(yàn)證原模型的網(wǎng)格劃分是否滿足要求。

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

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

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

3）切割邊界應(yīng)遠(yuǎn)離應(yīng)力集中區(qū)域；

4）原整體模型相應(yīng)位置處網(wǎng)格有必要進(jìn)行足夠程度的細(xì)化。

基于船舶結(jié)構(gòu)的特性及Patran對殼、梁單元的強(qiáng)大支持，行業(yè)內(nèi)一般在Patran中建立全船有限元模型，此時若對相關(guān)疲勞熱點(diǎn)進(jìn)行斷裂力學(xué)評估，需考慮結(jié)合Patran和Ansys兩套軟件的子模型技術(shù)，即通過在Ansys中建立與Patran原整體模型待評估位置對應(yīng)的子模型進(jìn)行分析，與單純在Ansys內(nèi)部實(shí)現(xiàn)邊界切割和位移插值的子模型技術(shù)相比，除需滿足基本操作要求，還有以下問題：

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

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

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

1 子模型技術(shù)簡介

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

圖 1 應(yīng)力強(qiáng)度因子計(jì)算流程Fig. 1 The procedure of SIF calculation

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

1）在Patran中采用shell單元建立厚度為t的平板模型，對板左端施加六自由度約束，對板右端施加垂直于板厚方向的拉應(yīng)力。對板中心區(qū)域網(wǎng)格進(jìn)行t×t細(xì)化，細(xì)化區(qū)域不少于10t范圍，為實(shí)現(xiàn)對三維表面裂紋SIF的求解，需將細(xì)化區(qū)域的板單元拉伸為實(shí)體單元，并在殼體邊界使用RSSCON Surf-Vol型MPC關(guān)聯(lián)對應(yīng)節(jié)點(diǎn)自由度，然后提交計(jì)算。

2）在Ansys中建立與Patran體單元部分對應(yīng)的子模型，并對裂紋進(jìn)行網(wǎng)格劃分。然后輸出Patran中實(shí)體單元邊界節(jié)點(diǎn)的位置信息，經(jīng)坐標(biāo)轉(zhuǎn)換及插值施加為Ansys子模型的邊界條件，即可實(shí)現(xiàn)SIF的求解。

2 殼體邊界MPC的創(chuàng)建

如圖2所示，Patran模型中體單元和殼單元的邊界節(jié)點(diǎn)需要使用RSSCON Surf-Vol型MPC進(jìn)行關(guān)聯(lián)。此MPC需將上下層邊界節(jié)點(diǎn)（即體單元邊界節(jié)點(diǎn)）設(shè)為獨(dú)立點(diǎn)，中間層邊界點(diǎn)（即殼體聯(lián)合邊界節(jié)點(diǎn)）設(shè)為非獨(dú)立點(diǎn)。

圖2 邊界MPC的創(chuàng)建Fig.2 Creation of boundary MPC

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

2.1 平面逐周分層法

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

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

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

通過Patran輸出切割邊界節(jié)點(diǎn)的位置信息，并使用插件進(jìn)行逐周分層處理，即可快速實(shí)現(xiàn)各層邊界點(diǎn)的重新排序，并匯總輸出準(zhǔn)確對應(yīng)的各層節(jié)點(diǎn)號。相對GUI方式，視邊界點(diǎn)節(jié)數(shù)量通?？蓪PC創(chuàng)建效率至少提高百倍。

2.2 曲面逐周分層法

除常見的平板結(jié)構(gòu)，船舶也存在諸多典型曲面結(jié)構(gòu)物，本節(jié)針對部分應(yīng)用廣泛的典型曲面結(jié)構(gòu)使用逐周分層法。

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

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

除圓柱殼模型應(yīng)用廣泛，橢圓柱殼由于其在裝備布置和后屈曲性能方面的優(yōu)勢，也逐漸被應(yīng)用于潛艇設(shè)計(jì)[23]，如941型戰(zhàn)略核潛艇。圖5為實(shí)際厚度為t的橢柱形殼體混合結(jié)構(gòu)。在橢柱坐標(biāo)系中，外層橢柱面，中層橢柱面及內(nèi)層橢柱面分別對應(yīng)3個底橢圓面，其焦點(diǎn)分別為A1，B1，A2，B2和A3，B3，其中中層橢柱面的底橢圓長軸為a，將底橢圓面中心O取為橢柱坐標(biāo)系原點(diǎn)，以上參數(shù)視作已知。

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

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

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

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

蘇聯(lián)部分Y級潛艇艇首采用橢球殼結(jié)構(gòu)，圖7為實(shí)際厚度為t的橢球殼體混合模型。在橢球坐標(biāo)系中，外層橢圓截面，中層橢圓截面及內(nèi)層橢圓截面的焦點(diǎn)分別為A1，B1，A2，B2及A3，B3，中層橢圓截面長軸為a，原點(diǎn)O取為橢球中心，以上參數(shù)視為已知。

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

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

2.3 注意事項(xiàng)

針對本文所提逐周分層法，尚有幾點(diǎn)說明：

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

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

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

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

3 基于矩陣法的坐標(biāo)變換

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

1）兩坐標(biāo)系各平面間的相對轉(zhuǎn)角其實(shí)并不直觀，即便先通過計(jì)算方向向量來求得各向轉(zhuǎn)角，也費(fèi)諸多周折。

2）體單元中或存在成千上萬的節(jié)點(diǎn)，若每次求解SIF均通過Ansys APDL處理坐標(biāo)變換，在批量求解各尺寸裂紋的SIF時，也會增加前處理的負(fù)擔(dān)。

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

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

4 基于映射網(wǎng)格的裂紋建模

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

5 子模型技術(shù)的驗(yàn)證

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

圖8 疲勞熱點(diǎn)及子模型Fig.8 Fatigue hotspot and sub-model

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

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

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

6 結(jié)語

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

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

圖9 橫浪工況Fig.9 Beam sea load case

圖10 裂紋位置及構(gòu)件尺寸Fig.10 Crack location and component size

圖11 I型裂紋SIFFig.11 Mode-ISIFof crack

2）針對不同坐標(biāo)系間節(jié)點(diǎn)位置轉(zhuǎn)換困難問題，提出更加直觀的矩陣法并編寫插件FEM_coor_transfer_V1.0進(jìn)行實(shí)現(xiàn)，且成功應(yīng)用于艙體疲勞熱點(diǎn)位置處的節(jié)點(diǎn)坐標(biāo)轉(zhuǎn)換。

3）針對非常見構(gòu)件求解裂紋SIF時自由劃分方式的局限性，提出更具優(yōu)勢的映射劃分方式并編寫插件Crack_mapper_V1.0進(jìn)行實(shí)現(xiàn)，且成功應(yīng)用于艙體疲勞熱點(diǎn)處各尺寸表面裂紋的SIF求解。