基于ALE的礦粉貨物液化晃蕩問題并行數值模擬

丁峻宏，金允龍，王 惠（上海超級計算中心，上海003；上海船舶運輸科學研究所，上海0035）

基于ALE的礦粉貨物液化晃蕩問題并行數值模擬

丁峻宏1，金允龍2，王 惠1
（1上海超級計算中心，上海201203；2上海船舶運輸科學研究所，上海200135）

含水量較高的礦粉貨物在海上運輸過程中易出現液化，形成自由表面并使晃蕩現象加劇，嚴重威脅船舶運輸安全。針對船載液化礦粉晃動和艙壁沖擊問題，采用ALE有限元方法對其進行了細致建模和計算模擬，從三維角度考察了在船艙一定裝載率和運動狀態下液化礦粉的晃蕩現象和特性；同時，借助不同仿真軟件，對計算結果的合理性和準確性進行了相互比對和分析。模型求解借助了高性能計算資源，以解決問題求解時間長和多組計算工況帶來的大規模計算需求；結合所建計算模型特點和流固耦合特性，研究了多核環境下兩種不同區域分解策略和實現方式，通過并行計算性能數據比較分析，以探求更為合理的并行加速策略。

礦粉液化；ALE；晃蕩；流固耦合；并行計算；區域分解

0 引 言

隨著中國經濟社會的快速發展，對礦產資源的需求一直處于高位，許多重要能源和礦物原材料持續表現出高峰需求。近幾年來，從國外進口的鐵礦和有色金屬精礦供應日趨增多，以鎳礦為例，一方面有限的國內資源供應難以滿足市場需求，另一方面鎳又是高品質不銹鋼生產過程中的重要添加劑，工業需求催生了國外鎳礦的大量進口。然而，鎳礦貨物在運輸過程中，時常由于種種原因導致含水量超標，從而在遇到搖擺和振動時發生水分析出礦體表面的現象（亦稱礦粉"液化流動"），形成自由液面，導致船體穩性損失及局部強度超載，嚴重威脅船舶運輸安全[1]。

僅以2010年底到2011年底為例，在近一年時間內，在中國周邊海域就有四艘船舶相繼由于礦粉貨物液化造成船體穩定性不足最終導致船體傾覆沉沒，包括“建富星”號、“南遠鉆石”號、“宏偉”號和“皇后”號，事故的高度頻發引發業界對這種貨物的運輸安全性日趨重視，各國行業組織[2]也開始致力于研究制定更加統一實用的安全規范。

不同礦源地和不同含水量條件下的礦粉物性存在差別，可以在實驗室條件下開展相應的物性研究和小比例液化礦粉晃蕩試驗模擬，但是實驗模型力學行為和實際船舶運輸條件下相比仍有一定的區別，因此借助數值仿真手段進行液化礦粉晃蕩和安全性研究具有其獨特和重要的作用。

由于船艙晃蕩引起液化礦粉或礦砂的流動，涉及到非線性、大變形和流固耦合等復雜物理特性。過去有不少針對LNG液艙晃蕩開展的研究，可以從技術手段上提供借鑒，比如說有的研究采用歐拉法結合多相流來模擬，并選用流體體積法（VOF法）跟蹤自由表面[3]；有的研究采用任意拉格朗日—歐拉法（ALE）來求解含自由液面的液體晃蕩問題[4]；有的采用無網格的拉格朗日法來模擬，對大變形區域的模擬和自由表面的捕捉也有其獨到之處[5]。但是，從現有文獻綜合來看，目前國內外使用數值仿真方法來開展粘稠泥沙或礦砂流動模擬的先例不多[6-7]；雖然國內近年來也開始有一些從事液化礦物運輸安全性方面的研究出現，但多專注于船舶傾覆機理分析和預防監控應急措施[8-9]，鮮見開展液化礦粉船艙晃蕩數值模擬的研究案例。

此外，以往類似研究都基于二維分析或者小比例研究模型，如若擬實建立三維船艙模型和細致的空間網格，則可以對液化礦粉船艙晃蕩過程進行更為立體真實的仿真和考察。但大規模的流固耦合計算和深入比較分析，卻并非一般計算機的計算能力能夠企及。

本文針對某型運輸船舶中一定裝載率下的液化礦粉，建立了大規模全三維流固耦合有限元模型，計算則基于ALE有限元顯式算法，并借助高性能計算機和多核并行求解工具LS-DYNA軟件，研究了液化礦粉在晃動過程中的力學行為和周期性變化規律，并與其他仿真結果進行了對比。

1 原理與方法

1.1 ALE算法和VOF法

對于呈液化流動狀的泥沙、礦砂或礦粉計算模型來說，拉格朗日算法將網格處理為隨材料一起運動，界面跟蹤和邊界設置比較容易，但很難處理由此引起的網格大變形和大位移；歐拉算法將材料處理為從固定網格中穿越，避免了網格畸變的同時卻不利于追蹤界面和邊界，難點還有多相流的模擬[10]。ALE算法系統計算所有單元的質量、動量和能量輸運量、單元密度和速度等參量，并隨時間步進行更新，最后采用更新后的密度和狀態方程中單位內能來計算單元中的壓力值。

ALE算法在拉格朗日和歐拉兩種坐標系基礎上，又引入任意參考坐標系，使得物質導數描述為：

式中：Xi是拉格朗日坐標，xi是歐拉坐標，wi是物質速度v和網格速度u之間形成的相對速度差。

公式（2）至公式（4）表述了ALE算法的守恒控制方程：

a.質量守恒方程

b.動量守恒方程

c.動量守恒方程

式中：ρ為流體密度，σi j為應力張量，b為單位質量的體積力，E為單位質量的總能量。

有限差分法被用于實現每個節點任意時刻的速度u和位移x更新：

式中：M為對角質量矩陣，Fint為內力矢量，Fext為與體積力和邊界條件等有關的外力矢量。

VOF法通過研究網格單元中流體和網格體積比函數Vr來確定自由面，追蹤流體的變化，被廣泛用于求解流體和固體力學中的很多問題。在晃蕩計算模型求解時，Vr取值為1或0時分別代表全部為指定相流體所占據單元和無指定相流體單元兩種狀態，部分填充單元（0＜Vr＜1）構成了自由液面，界面追蹤算法在網格更新和物質輸送之前執行。

1.2 多CPU/核并行求解

流固耦合工程問題非線性特點強，計算規模大而耗時久，最好能借助并行計算手段[11]。目前有限元程序并行求解時，多采用大規模并行處理（MPP）模式。如圖1所示，一般通過區域分解方法（Domain decomposition method，簡稱DDM）來處理分析對象，即把整個模擬問題區域劃分成很多相對小的求解子區域,把每個子區域分配給不同的CPU處理器（或核心，Core）分別進行求解，區域分解允許每個核單獨求解一部分問題，而核之間通過交互機制進行數據交換，最后將所有子區域的解匯總得到整個區域的全局解。

圖1 并行求解的模型分區和CPU布局Fig.1 Model of domain decomposition for parallel computing and CPU distributing

在對區域分解并行求解的計算效果進行比較和評價時，常見的考核性能評價指標包括模型計算總耗費時間TC，加速比SN和并行效率EN。其中SN和EN分別用來考量并行加速處理效果和并行計算時整體系統的資源利用率，即：

2 應用實例

2.1 計算模型

本文開展液化礦粉船艙晃蕩數值模擬研究的基本計算模型見圖2（隱去部分艙殼單元），圖3顯示的是后續為開展壓力時程對比而選定的艙壁中下部監測點位置。整體模型基本尺寸為29.5 m（長）× 32.3 m（寬）×18.4 m（高），所裝載的礦粉高度距內底板8.0 m，礦粉液化后假定上層還出現了深0.3 m的稀泥漿，船艙重心距內底板5.3 m。礦粉、泥漿和空氣各對象有限元模型均采用ALE六面體體網格劃分，其中泥漿網格還進行了局部加密；船艙采用四面體殼單元來模擬，整個計算模型單元和節點總數均在24萬左右。

圖2 礦粉貨物晃蕩計算有限元模型Fig.2 FEM model for sloshing of ore fines cargo

圖3 艙壁壓力監測單元位置示意圖Fig.3 Diagram of pressure monitoring elements at bulkhead

各對象的基本物理屬性見表1。其中船艙艙殼用MAT_RIGID剛性材料模型來模擬，而礦粉、泥漿和空氣三者都采用MAT_NULL材料模型，同時使用EOS_GRUNEISEN狀態方程[12]定義可壓縮流體材料的壓力：

式中：p為介質壓力，ρ和ρ0分別為當前和初始密度，μ=ρ/ρ0-1，C為沖擊波速，γ0為GRUNEISEN常數，S1、S2和S3是vs-vp曲線的斜率系數，a是壓縮時γ0的一階體積修正，E為內能。以本文液化礦粉為例：可取波速C=1 647 m/s，S1=1.921，S2=-0.096，S3= 0，a=0，γ0=0.35。

在裝載液化礦粉船艙晃蕩的模擬過程中，設定船艙繞重心橫向往復擺動，最大搖擺角度為10°，取晃蕩激勵周期T為10 s，每個工況的計算時間取為3T。

表1 各建模對象的基本材料參數表Tab.1 Basic material parameters of each modeling object

2.2 并行求解

本文建立的計算模型規模不小，每個工況問題物理時間久，而且ALE算法求解時間也比拉格朗日算法要長得多。如果使用普通工作站（4CPU Cores；內存8～16 GB），計算一個工況至少需要近7天，如果計算工況多的話，則對計算效率、磁盤空間直至研究的順利開展都提出了挑戰，因此必須要借助高性能計算資源。上海超級計算中心搭建有專門針對工業和工程領域應用需求的“蜂鳥”集群，該集群包括65臺HS23刀片計算節點，每個節點含兩個Intel八核Intel處理器（主頻2.6 GHz,共享64 G內存）。

利用LS-DYNA中ALE算法進行問題求解時，采用坐標遞歸對分方法（Recursive Coordinate Bisection，RCB）來進行計算模型的多核區域分解處理。該方法將網格沿垂直于模型最長坐標方向將模型對分，將相鄰邊界節點復制分布到兩側區域，持續遞歸分割直至滿足分割收斂條件。通常采用程序推薦的基于ALE單元均分策略來組織進行計算區域的剖分方案A（見圖4（1））；進一步考慮到船艙尺寸特征、晃蕩方向性特點以及流固耦合界面的均衡分配需求，本文在此基礎上又提出了新的子區域剖分策略方案B（見圖4（2）），并就這兩種不同區域分解方案對應并行加速效果進行測試對比。

圖4 不同的區域分解策略及分區結果Fig.4 Different domain decomposition strategies and results

3 計算分析

3.1 船載液化礦粉晃蕩分析

一般搖擺船艙內液體的晃動可以認為是駐波、行波、水躍和渦旋等形式的組合，如何從三維的角度去再現艙內裝載礦粉液化后隨船搖擺下的晃蕩現象非常重要，由于篇幅關系，本文在此僅給出了幾個特定時刻艙內裝載礦粉晃蕩時自由表面變化和密度分布云圖（見圖5）。可以發現，與一般艙載均質低密度低粘度液體晃蕩現象和特點有所不同，絕大部分液化礦粉的晃動形式偏向駐波，而其上層部分礦粉和表層稀泥漿在兩側艙壁之間往復流動現象更似行波，從三維視角來看自由表面流動平緩且不規則（見圖6）。另外，液化礦粉在艙內晃蕩流動時，與船舶本身橫搖運動之間有一定相位差，該現象在船舶從搖擺至兩側最大傾角繼而向平衡位置運動時最為明顯，若在海上持續受到不利外載作用并導致礦粉逐漸向一側傾斜聚集，終將影響船舶的運輸安全。

圖5 不同時刻礦粉貨物的自由表面形狀Fig.5 Free surface shapes of ore fines cargo at different time

圖6 礦粉貨物的流體密度云圖（t=25 s）Fig.6 Fluid density contour plot of ore fines cargo(t=25 s)

3.2 仿真數據對比和應急措施驗證

船載礦粉發生液化產生晃蕩時多發生在海上運輸途中，一般很難直接獲取來自現場的數據和資料。因此仿真結果數據的對比，只能借助于相似比例模型試驗或者其他仿真算法。將本文利用LS-DYNA軟件計算得到的艙壁監測點（見圖3）壓力時程曲線，與利用CFD軟件FLUENT計算得到的結果進行了對比，見圖7。兩組結果數據初始時刻壓力值有所不同是因為FLUENT模擬時考慮了靜壓，但曲線后續變化趨勢表明隨著船舶開始晃蕩后靜壓的影響逐漸降低，曲線整體周期性規律和最大峰值基本趨于一致，但相位和幅值略有差異，這是受不同流體模型參數設置以及不同軟件具體算法的影響。初步計算表明兩種軟件所得結果還是比較接近的，這可為借助不同方法開展晃蕩仿真數值模擬時的合理性和準確性評判提供相互驗證。

圖7 不同模擬軟件計算得到的壓力時程對比Fig.7 Comparison of pressure time history due to different simulation softwares

圖8 安裝隔板后礦粉貨物的流體密度云圖（t=25 s）Fig.8 Fluid density contour plot of l ore fines cargo after installing middle baffle(t=25 s)

在易液化礦粉裝船以及航行途中發生橫傾時，現場都有一些規范和應急舉措，基于現有三維仿真模型，還可以從中選取一些建模方便的方案開展驗證和評估，充分發揮數值計算模型的優勢。如圖8所示，有種措施是在艙內中部用木板縱向埋于礦粉表面并與艙前后兩側固定，與圖6相比，計算表明該措施的確可以降低晃蕩幅度，從而減小自由液面的影響。

3.3 并行求解性能分析

以船載液化礦粉晃蕩計算模型計算三個搖擺周期（30 s）為測試對象，從圖9（1）來看，1核和2核由于兩種區域分區結果無異，因此計算耗時也一致。4核時方案B開始略有不同，使用8核和16核時該方案節省計算時間現象最為明顯，到32核和64核時能繼續保持其比較優勢，表明其仍具有一定的可擴展性。從計算耗時具體數值來看，單核計算需要25天，雙核也要13天，不只是計算時間上令人難以接受，更重要的是對軟硬件計算環境工作穩定性也提出了很高要求，一般很難絕對保證問題能在預期時間順利算完；使用16核后，計算耗時可顯著降低到3天之內，給多工況求解和深入研究提供了時間上的可能性。

圖9 不同區域分解策略下并行性能數據對比和分析Fig.9 Parallel performance data comparison and analysis with different decomposition strategies

進一步從圖9（2）分析來看，兩種分區方案對應的加速比效果均隨核數增加呈非線性上升趨勢，其中方案B自4核之后始終保持領先并略呈擴大趨勢，表明其負載均衡得到了進一步改善。但從并行效率數據來看，32核之后兩種方案并行效率都已開始回落至50%甚至更低，從計算資源的有效配置和使用角度來看，不建議再使用過多的計算核心數。

因此對于本文以及類似同等規模計算模型而言，求解時使用16～32核應該是比較合理的選擇，單工況求解時所耗時間將會大大降低，2～3天就可完成計算；而在此基礎上若采用方案B執行計算模型的區域分解，則可將單工況總計算耗時在常用方案A的基礎上再進一步縮減3～10小時。

4 結 論

本文詳細闡述了如何利用高性能計算軟硬件資源，與流固耦合算法以及并行計算技術緊密結合，就船舶遠洋運輸過程中礦粉液化后產生的晃動沖擊問題開展深入研究，給出了具體實現方法和相關參數設置，得到了一些有意義的計算結果并對此進行了深入分析，最后有以下認識：

（1）結合ALE算法和VOF法的三維有限元方法，可有效捕捉自由液面的變化，模擬液化礦粉受船舶搖擺引起的往復晃蕩和艙壁沖擊等力學行為，驗證減小晃蕩的實際措施，并為艙壁結構承壓分析打下基礎；不同數值模擬軟件計算結果可以相互進行比對，以間接驗證該問題仿真方法的可行性和正確性。

（2）以高性能計算資源為支撐，可對大規模船載液化礦粉晃蕩沖擊問題進行細致建模和深入分析，并將問題求解時間的影響控制在有利范圍內；開展并行計算時，應該根據求解問題類型和模型力學特點選用合理的區域分解方案和合適的多核并行加速策略，以實現計算資源的優化配置和最佳利用。

[1]雷 海.紅土鎳礦粉的安全運輸[J].航海技術,2011(1):27-28. Lei Hai.Transprotation safety of laterite nickle ore fines[J].Marine technology,2011(1):27-28.

[2]Class N K.Guidelines for the safe carriage of nickel ore(Second edition)[M].Tokyo:Nippon Kaiji Kyokai,2012.

[3]劉楨兵.基于VOF法的液艙晃蕩數值模擬及載荷計算[J].重慶理工大學學報（自然科學）,2011,25(3):24-29. Liu Zhenbing.Numerical simulation of liquid tank sloshing and pressure calculation based on volume of fluid method[J]. Journal of Chongqing University of Technology(Natural Science),2011,25(3):24-29.

[4]管延敏,葉恒奎,陳慶任.基于ALE有限元法二維液體晃蕩的數值模擬[J].船舶力學,2010,14(10):1094-1099. Guan Yanmin,Ye Hengkui,Chen Qingren.Numerical simulation of 2D sloshing with Arbitrary Lagrangian-Eulerian finite element method[J].Journal of Ship Mechanics,2010,14(10):1094-1099.

[5]陳正云,朱仁慶,祁江濤.基于SPH法的二維液體大幅晃蕩數值模擬[J].船海工程,2008,37(2):44-47. Chen Zhengyun,Zhu Renqing,Qi Jiangtao.Numerical simulation of sloshing in two dimensional liquid tank based on SPH method[J].Ship&Ocean Engineering,2008,37(2):44-47.

[6]Yim S C.3-D Wave-structure interaction with coastal sediments-A multi-physics/multi-solution-techniques approach[R]. Corvallis:Oregon State University,2008:1-10.

[7]Karmakar S,Kushwaha R L.Simulation of soil deformation around a tillage tool using computational fluid dynamics[J]. Transactions of the ASAE,2005,48(3):923-932.

[8]王 威.淺談海運精礦粉的安全運輸[J].南通航運職業技術學院學報,2009,8(2):63-65. Wang Wei.On safe transport of ore concentrate powder[J].Journal of Nantong Vocational&Technical Shipping College, 2009,8(2):63-65.

[9]趙月林,孟淑嫻.易流態化貨物安全運輸研究[J].大連海事大學學報,2012,38(4):15-18. Zhao Yuelin,Meng Shuxian.Safety transportation of easily fluidized cargoes[J].Journal of Dalian Maritime University,2012, 38(4):15-18.

[10]Aquelet N,Souli M,Gabrys J,et al.A new ALE formulation for sloshing analysis[J].Structural Engineering and Mechanics,2003,16(4):423-440.

[11]李 政,金先龍,亓文果.流體—結構耦合問題的有限元并行計算研究[J].計算力學學報,2007,24(6):727-732. Li Zheng,Jin Xianlong,Qi Wenguo.Study on the parallel algorithm for finite element analysis of fluid-structure interaction[J].Chinese Journal of Computational Mechanics,2007,24(6):727-732.

[12]LSTC.LS-DYNA keyword user’s manual[M].California:Livermore Software Technology Corporation,2012.

ALE-based parallel numerical simulation for sloshing problem of liquefied ore fines cargo

DING Jun-hong1,JIN Yun-long2,WANG Hui1
(1 Shanghai Supercomputer Center,Shanghai 201203,China;2 Shanghai Ship and Shipping Research Institute,Shanghai 200135,China)

Ore fines cargo with high moisture content is liable to liquefaction,making free surface emerging and intensifying sloshing behavior,and finally pose a threat to marine transportation security.With the benefit of Arbitrary Lagrangian-Eulerian(ALE)finite element method,detailed modeling and refined simulation were implemented to deal with shipping liquefied ore fines sloshing problem,and to investigate three-dimensional sloshing phenomena and characteristics with certain charging ratio and motion state. Two different numerical nodes were used to make comparison and analysis for rationality and accuracy of results.High-performance computing(HPC)resource was utilized to meet the challenge of large scale computing requirement due to time-consuming solving and load cases.In view of model size and fluid-structure intersection property,two domain decomposition plans were proposed to better improve load balance for multicore processors,and compared with parallel performance data to determine more reasonable parallel acceleration strategy.

liquefaction of ore fines;ALE;sloshing;fluid-structure interaction;parallel computing; domain decomposition

TP391.9 U661.3 U674.13+4.1

A

10.3969/j.issn.1007-7294.2015.08.006

1007-7294（2015）08-0927-07

2014-12-18

國家高技術研究發展計劃（863計劃）課題（2012AA01A308）；國家重點基礎研究發展計劃資助課題（2012CB723804）

丁峻宏（1977-），男，工學博士，高級工程師，E-mail:Jhding@ssc.net.cn；金允龍（1964-），男，研究員。