獨立C型罐LNG運輸船鞍座結構的數值模擬方法研究

婁本強，張洪軍，孫 立，董國忠

(1.江蘇科技大學 船舶與海洋工程學院, 鎮江 212100) (2.上海海威斯特保溫工程有限公司, 上海 201208) (3.中遠海運重工有限公司, 揚州 225211)

隨著我國LNG需求開始增加，在儲氣設施集約建設思想的指導下，我國正引領“碳中和LNG”(carbon neutral LNG)的行業實踐[1]，在完善沿海省份LNG 產業鏈布局的同時，逐步增加 LNG 儲運體系的建設力度.LNG的物流可以通過專用汽車或內河沿海船舶運輸到末端用戶，日本、挪威等國已有使用內河和近海船舶將遠洋燃氣資源由港口轉運、分銷至沿線LNG接收終端的先例[2].借鑒國外相關經驗，我國亦可以采用中小型LNG運輸船為紐帶，借助長江航運便利，進行支線運輸，供內陸城市使用[3].獨立式C型儲罐LNG運輸船具有裝載靈活、投資回收周期短等優點，適用于我國東南沿海及長江流域內河的LNG運輸，因此，對該船型的研究具有現實工程意義和經濟價值.

C型儲罐與船體結構之間相互獨立，鞍座結構支撐和限制儲罐移動起到保護的作用.在對C型罐LNG運輸船進行結構強度計算時，除計入波浪載荷外，還應考慮因加速、減速、轉彎和升沉等引起的慣性載荷項.《國際散裝運輸液化氣體船舶構造與設備規則》(簡稱IGC規則)[4]在慣性載荷合成方法上作了重要修改，用三維加速度橢球代替二維橢圓圖，雖然增加了計算難度，卻更接近實際，因此，如何合理地準確施加IGC規則中的載荷，會影響支承結構的數值仿真結果，以及對安全評估合理性的判斷.

文中以某22 000 m3獨立雙耳筒C型儲罐LNG運輸船作為研究對象，利用自編軟件工具AWQLoading[5]，將靜水載荷、航行中的等效波浪載荷及慣性載荷等加載到典型工況的全船有限元模型上，進行直接計算.文中摒棄線性迭代、轉向非線性有限元計算，討論接觸非線性邊界以及動/靜摩擦系數設置對鞍座結構的仿真效果，比較不同模型在模擬層壓木時對層壓木支反力與鞍座結構應力分布的影響.在現有的文獻中，鮮有文獻對如何計入接觸面間動/靜摩擦力的介紹，甚至直接忽略摩擦影響.文中的工作是對現有C型罐LNG運輸船結構強度計算的細化與補充，具有一定的參考價值.

1 LNG運輸船及C型儲罐支撐系統

1.1 LNG儲罐

國際海事組織(IMO)對儲罐進行了分類，包括整體型、內部絕熱型、薄膜型、半薄膜型、獨立型等.其中獨立型為自持式儲罐，采用全壓式或半冷半壓式，與船體結構非直接相連.通常將此類獨立儲罐劃分為A型、MOSS型(屬B型)和C型3種型式，如圖1.

圖1 常見的獨立型LNG儲罐及其所對應的全船有限元模型

C型儲罐為臥式貨艙容器，儲罐形狀特殊(圖2)，獨立放置于船體艙內，儲罐材料采用5%～9%的鎳鋼，使用聚氨酯泡沫材料包覆，確保其在運輸過程中，液貨保持-163℃的載運溫度.雙耳筒型儲罐較圓柱單筒型儲罐能提供較大的倉容，可以更充分地利用內部空間；較三耳型儲罐，雙耳筒型儲罐的技術和工藝相對成熟，因此，雙耳型儲罐的應用更加廣泛.

圖2 3種常見的C型儲罐型式

1.2 目標船型的艙容與結構布置

LNG運輸船是布置型船舶，艙容要求高，主尺度的選擇要從儲罐總布置及船舶性能等方面來考慮，結合吃水高、內河船舶航行視野寬闊等要求，選擇儲罐型式、儲罐尺寸、壓力要求、重量及鞍座結構等內容，與主尺度相互影響[6].文中目標船采用球鼻首型船艏，機艙設置在尾部.機艙前設置4個C型獨立液貨儲罐.貨艙臥式容器采用2個鞍座支承，其中1個是固定鞍座，另外1個是滑動鞍座，鞍座與容器間墊有層壓木[7].2～4號儲罐是常見的標準雙筒圓柱形，艙容均為5 570 m3；為配合船艏線型，1號儲罐通常采用的雙筒錐形，儲罐艙容約為5 300 m3.在滿載情況下，雙筒圓柱儲罐的液貨總質量約是3 066.5 t，雙筒錐形儲罐的液貨總質量約為2 680.4 t.在中拱和中垂工況下，船體梁將產生較大位移變形.C型獨立儲罐長度和體積都較大，距離舯剖面較近的儲罐的鞍座系統,將承受來自液貨和儲罐的慣性力，引起局部應力升高[8]，文中選取2號儲罐作為目標.

1.3 基于IGC規則的附加慣性加速度

雙耳筒型儲罐的支撐鞍座呈“碗”狀布置(圖3)，當面臨環境載荷時，由于慣性，儲罐會產生運動滯后，船體與儲罐之間發生耦合作用.在已發表的研究中有使用線性疊加載荷系數的子工況[9]或者用近似正余弦作用力分布[7]代替整體載荷.

圖3 C形獨立儲罐

按照IGC規則和相關壓力容器規范要求，設計中不僅要考慮低溫環境，同時還應考慮船體運動對貨罐的影響.在更新的IGC規則中，使用加速度橢球代替之前的二維橢圓.慣性加速度沿縱向、橫向和垂向的最大分量Ax，Ay，Az是根據北大西洋中在超越概率水平為10-8波浪載荷下得到：

(1)

(2)

(3)

式中：x為液艙重心至船舯的縱向水平距離；z為實際水線面至液艙中心的垂向距離；B為是船舯最大寬度；HGM為橫傾高度;Cb為方形系數;V為船速;K為系數，K=13×HGM/B，當K≥1.0時，取計算值，當K<1.0時取1.0.

對液艙內某參考點而言，內部壓力Pgd計算為：

(4)

式中：aβ為動載荷引起的加速度；Zβ為沿β方向至液貨艙的液柱高度，m；ρ為貨物密度，kg/m3.

之前，計算該點所對應的最大壓力是采用迭代法，間隔一定的傾斜角進行計算[10]，迭代法對雙耳圓筒型C型罐而言稍顯臃繁.文中采用解析法，用球坐標(φ,θ)表示aβ，見式(5)，利用矢量數量積理論，得到aβ·Zβ極值的解析解——式(6) (文獻[11]中應用于薄膜型LNG液艙).解析解具有精度好、運算效率高等優點.

aβ=(0,0,1)+(Axcosφ,Aysinφcosθ,Azsinφsinθ)

(5)

aβ·Zβ=(zA-zi)+

(6)

值得注意的是縱艙壁通常設置為水密，左右兩側的壓力需要分開計算.使用PTC.MathCAD中編程實現解析解的圖形化結果，輸入2號儲罐右側外型，將數據和計算集成在一個工作表中，得到橢球無因次加速度矢量與重力夾角β= 35.6°，并得到可視化結果：液貨高度Zβ=13.74 m，對應的加速度aβ=1.50 m/s2(圖4).

圖4 第2號液貨艙右艙的PTC.MathCAD軟件的可視化結果

1.4 基于設計波法的全船有限元直接計算

采用Patran對全船進行有限元建模，使用Nastran求解.總體坐標系原點設置于船體的形心，根據X、Y、Z左手坐標系統建立模型，其中X軸沿船長方向船艏為正；Y軸沿船寬方向右舷為正；Z軸由基線開始，以沿型深方向向上為正.船體的殼板、船體強框架腹板、底部肋板和液貨艙內強框架腹板等結構使用板殼單元(shell)；強框架面板、扶強材、縱骨及船底、舷側上的加強筋和液貨艙支撐構件使用梁單元(beam)模擬.采用等效設計波方法將波浪附加彎矩和靜水彎矩施加在船體上的，其中波浪載荷是準靜態載荷，采用水動力計算軟件PRECAL?計算目標船在規則波作用下的船體運動響應、波浪載荷，得到以船舯彎矩為控制參數的等效設計波，波幅為垂向波浪彎矩最大值與典型規則波中船舯彎矩的傳遞函數比值.

英國勞氏船級社 (LR)的指南[12]要求對規定吃水下滿載和隔艙裝載等狀態進行直接計算，滿載狀態是實際營運中的主要狀態，通常為進一步確保結構強度安全，需要校核隔艙裝載狀態.全船有限元模型的加載過程復雜，手工操作耗時費力，自編AWQLoading輔助加載程序，通過編輯Patran輸出*.bdf文件，實現調整空船質量、加載液貨艙和壓載水質量等功能，設置二次開發接口，將等效設計波載荷以法向力的形式加載在典型工況下目標船的船殼濕表面上，與靜水浮力疊加.和船級社的商業程序(如NAUTICUS?，POSEIDON?和SHIPRIGHT?)進行對比，性能相近.

在AWQLoading中設置參數交互控件，讀入合成加速度，將慣性載荷以場的形式加載，得到目標船在“準靜態隔艙裝載動橫搖”的結果(圖5)，該等效設計波主導目標參數，包括波浪垂向中拱彎矩和指定點最大壓力對應的橫搖加速度.

圖5 準靜態隔艙裝載動橫搖工況

2 C型儲罐支承鞍座結構的模型化

2.1 C型儲罐的支承系統

儲罐通過支承系統與船體非直接相連，貨物及儲罐的重量及其動載荷作用在鞍座上.為防止熱脹冷縮對貨罐產生的附加作用，通常在船底設置雙聯鞍座形式支撐(圖6)，固定鞍座位于艉部，滑動鞍座于艏部.隨著C型儲罐設計容積的增加，也有工程師提出采用多聯鞍座支撐的方案[6].支承系統沿C型罐的周向布置，形成150°左右的包角，提供有效支撐的同時，還隔斷了儲罐和船體結構之間的熱量交換.層壓木通常選用北歐地域生長的木材，浸漬樹脂，削成薄片，經高溫和高壓加工成型.成型后的層壓木質地堅實、不易開裂[13]，還具有隔熱效果好、老化速度慢等優點.

圖6 雙耳C型LNG船雙聯鞍座布置

固定鞍座的面板上設有兩道擋板形成的槽，由伸入層壓木凹槽內的止移鍵板固定層壓木，在層壓木與鞍座面板之間和液貨艙殼體之間填充環氧樹脂膠泥，起到調平和固定儲罐的作用；滑動鞍座處層壓木分為上下兩部分，由位于鞍座面板及儲罐筒體上的小擋板固定(圖6).在儲罐外側沿鞍座周圍布置不銹鋼片，緩解剪切力.

2.2 鞍座結構的剛度系數

計算鞍座結構的整體等效剛度以兩根串聯的彈簧為例(圖7(a))，彈簧A的剛度為K1，彈簧B的為K2，串聯后彈簧組的總剛度為K.在外力F作用下，彈簧伸長量分別是x1、x2，彈簧組的總伸長量為x.假設不考慮彈簧組自重，變形均在彈性限度內，彈簧A與彈簧B受力相等，即F=K1·x1，F=K2·x2，那么K1x1=K2x2，則彈簧組的總伸長量x=x1+x2，其中x1=F/K1，x2=F/K2，所以，x=F/K1+F/K2，用公式表示為：

K=K1·K2/(K1+K2)

(7)

圖7 串聯彈簧與固定鞍座結構示意圖

結合上述串聯彈簧關系推廣至多個彈簧串聯，如包含不銹鋼片、層壓木、環氧樹脂等的鞍座結構(圖7(b))的等效剛度：

(8)

式中：H抗壓方向的高度；A構件橫截面積；E材料彈性模量;下標1，2，3分別代表不銹鋼片、層壓木、環氧樹脂[14].

2.3 鞍座面板與層壓木的面間摩擦力

當鞍座面板與層壓木相互接觸并發生相對運動或有運動趨勢時，接觸面上會產生阻礙運動的摩擦力f，按摩擦條件分為靜摩擦(static friction)與動摩擦(sliding friction)兩種.當只存在運動趨勢而未發生相對運動時，面間存在與切向力互為反力的靜摩擦力；當靜摩擦力達到與法向力W和靜摩擦系數μs的乘積相當程度時，會發生相對運動，摩擦力將“由靜轉動”，大小是法向力W與動摩擦系數μd的乘積，F為剪力，且μs>μd.

(9)

3 數值模擬方法

雖然是否考慮支撐系統內的接觸和摩擦非線性對船體總縱強度的影響不大[12]，然而儲罐自身及液貨重量通過鞍座集中傳遞給船體梁，會在鞍座位置處產生較大剪力和集中載荷.在校核鞍座局部強度時，其對應力分布的影響還是值得探討的.無論IGC規則還是各船級社并沒有對層壓木的模擬給出強制性要求，文中將對數值仿真方案進行介紹和討論.

3.1 臨時桿單元法

臨時桿單元法是模擬層壓木的傳統模型方法，該方法將層壓木等效為彈簧或桿單元[15].使用臨時桿單元法進行線性迭代：計算并判斷單元的拉/壓狀態，修改模型，刪除受拉的層壓木單元，更新模型并再次計算直至無層壓木單元受拉為止.迭代次數視工況而定，通常經過2～5輪即可得到收斂結果，若計算模型太復雜將導致迭代步驟增多，人工操作容易出錯.推薦借助軟件工具，設置自動迭代功能，對Patran底層數據文件進行二次開發：讀取*.bdf、判斷拉/壓狀態、修改模型、調用Nastran求解器、提交任務等步驟[16].

3.2 線性間隙單元法

MSC?還提供一種可以在求解器內部完成迭代計算的單元——線性間隙單元[17].對該單元的工作原理進行如下簡述：從隨機選取一組向量開始，一些節點處于接觸狀態，保證接觸部位處于受壓狀態，且沒有侵徹；若檢測到受限循環，即迭代返回到前一個狀態，重啟計算直至收斂.在文獻[16]中詳細介紹該單元的理論、控制方程及工作原理，并驗證了線性間隙單元的計算與使用臨時桿單元法的計算結果，發現二者等效.線性間隙單元仿真的優點在于不需要人工介入，設置簡單；問題在于對某些復雜工況會產生計算不收斂的現象.

3.3 計及摩擦力的非線性間隙單元法

最理想的方案是既能保證計算的收斂性，又可以考慮層壓木的“開/關”狀態以及端面的動/靜摩擦效果，非線性間隙單元(GAP)正好滿足上述要求[16].在單元卡片CGAP中定義(圖8)幾何信息，初始間隙用U0來表示，載荷通過層壓木傳遞，端點A、B的位移分別是UA、UB，剛度取決于端點的相對偏移量(即UA-UB).單元處于自由狀態時(UA-UB≥U0)，分離剛度為KB(一般取KB= 10-14KA)；當偏移量小于初始間隙時，剛度恢復到接觸剛度值KA[18]，在單元屬性卡片PGAP中除了賦值接觸時的法向剛度KA、切向剛度KT和分離時的剛度KB，還可以定義動/靜摩擦系數(μs和μd).

圖8 GAP單元模擬層壓木

(10)

在NLPARM控制卡中對收斂條件的參數進行設置，通過“更新剛度”和“子步內迭代次數”控制計算效率，避免計算發散或歧義.提交計算后，求解器會在不刪除單元的前提下，自動判斷，更新剛陣，得到收斂的計算結果.

4 計算結果與分析

4.1 層壓木的支反力結果

結合加載完成的C型LNG運輸船的全船模型，共對5種層壓木的數值仿真方案進行對比：只考慮層壓木接觸開/關，使用方案①“傳統臨時桿單元法”和方案②“非線性間隙單元法模擬層壓木”；方案③“在非線性間隙單元中賦值動/靜摩擦系數，μs=0.5和μd=0.25[9]”；方案④”假設無窮大靜摩擦系數 (μs=∞)”；方案⑤“假設較大靜摩擦系數(μs=50)的情況，忽略動摩擦阻力(μd=0)”.將滿載垂蕩工況下，支撐第2號儲罐的滑動與固定鞍座層壓木單元的支反力F結果繪制成圖表，見圖9，其中x軸為層壓木沿周向布置距船舯的距離d，單位為m.

圖9 滿載動垂蕩工況，鞍座支反力結果

同理，隔艙動橫搖工況，2號儲罐下的滑動與固定鞍座層壓木的支反力結果見圖10.

圖10 隔艙裝載動橫搖工況，鞍座支反力結果

觀察不同數值模擬方案下層壓木支反力的異同，發現方案①和②的結果相近，可知當僅考慮接觸非線性時，使用GAP仿真與傳統方法的結果是一致的，驗證了GAP單元的可靠性.方案③是最接近工程實際的模擬，最大支反力小于方案①和②，且分布緩和.方案④、⑤是假想模型：方案④與不計入摩擦力方案的結果相近，究其原因可能是由于無窮大的靜摩擦力限制了節點的位移，使計算回歸到僅考慮接觸非線性的模型；當設置較大靜摩擦系數(方案⑤，μs=50)時，接觸面間將只存在相對運動的趨勢，不存在相對滑動，可以看出周向靜摩擦力將分擔很大一部分的徑向力，此假設模型的支反力最小.因此，對所研究的C型LNG船而言，鞍座接觸面之間的摩擦對保護儲罐和船體結構是有益的，靜摩擦力抑制儲罐移動，動摩擦力將減緩儲罐的相對滑動.

4.2 儲罐限位裝置應力結果

鞍座結構的應力大小及其分布是反映局部強度的重要指標，文中選取方案②、③，觀察應力結果及其分布.兩種方案的區別之處在于是否在定義單元屬性中設置動/靜摩擦系數.

滿載垂蕩工況下，鞍座結構Von mises應力云圖如圖11，可知最大應力出現在固定鞍座上龍骨位置的肘板與支座面板相交處.若不計入摩擦效果，鞍座結構的最大應力為227.4 MPa；若計入動/靜摩擦力，鞍座結構的最大應力將有所緩解，降至213.2 MPa.

圖11 滿載動垂蕩工況，鞍座的應力云圖

另外還發現兩個模型中的最大應力出現在同一位置，其他較大的應力出現在鍵板、鍵板肘板、鞍座面板相交角處.當鞍座結構應力超過許用應力設計值時，可考慮在鞍座面板下設置支撐加強，改善應力分布，提高局部強度.

隔艙裝載動橫搖工況的鞍座受力云圖如圖12所示：若不計入摩擦力，最大應力出現在固定鞍座的基座結構處，值是120.0 MPa；若考慮動/靜摩擦力，最大應力轉移到滑動鞍座的端部，值是131.8 MPa.滑動鞍座結構的應力峰值比固定鞍座結構處應力峰值大.

圖12 隔艙動橫搖工況，鞍座的應力云圖

觀察有限元的應力結果可知動/靜摩擦力對局部強度校核是有影響的.支撐端部相較其他位置更易發生損壞，應該在此處加強支撐或者優化結構形式，防止出現應力過大.此處列舉的只是計算的一個方面，還有加強環、重磅板、真空環、縱艙壁和中部殼體及前后封頭等各位置的應力比較，就不一一贅述.

5 結論

(1)使用一維桿單元或者彈簧單元模擬層壓木，計算結果偏于安全，但存在一定誤差.采用本文提出的對Patran的數據文件進行編輯的后處理方法可取；利用編程實現循環迭代，縮短人工操作時間，提高計算效率.

(2)在C型LNG運輸船鞍座結構的仿真中，引入非線性間隙單元(GAP)，并計入接觸非線性和動/靜摩擦力.雖然對計算機硬件要求提高，也對操作者有限元計算的經驗和熟練程度的要求提高，但可以較真實地描述實際工程情況，降低引入人為誤差的風險，是較理想的模擬方案，可以作為安全評價或改進結構設計的依據.