基于SPH 方法的黃河破冰船冰阻力數(shù)值模擬分析

鄭興，田治宗，謝志剛，張寧波

1 哈爾濱工程大學 船舶工程學院，黑龍江 哈爾濱 150001

2 黃河水利科學研究院，河南 鄭州 450003

0 引 言

黃河由于其所在地理位置緯度較高，冬季嚴寒，因而在冬季會出現(xiàn)結冰現(xiàn)象，特別是在黃河最北端的內(nèi)蒙古流域，在封河初期以及開江解凍期，浮冰在水流的作用下易在河流狹窄彎曲段出現(xiàn)堆積，形成冰排、冰壩，從而堵塞河道，造成凌汛災情，給沿岸人民的生活和財產(chǎn)造成巨大損失[1]。為了應對冰凌災害，需要采取破冰措施進行疏通，目前的破冰措施主要有爆破破冰[2]、破冰船破冰[3]等方法，其中破冰船因具有機動性強、成本低等優(yōu)點，能在部分河段的破冰疏通方面發(fā)揮重要作用。而采用有效的數(shù)值分析方法對黃河破冰船的冰阻力進行估算，對于破冰船的設計和操縱具有重要作用。

近年來，國內(nèi)外一些學者采用數(shù)值模擬方法對冰與結構物的相互作用問題進行了分析，現(xiàn)有的用于模擬冰破壞的數(shù)值方法主要有有限元方法（FEM）、離散元方法（DEM）、近場動力學（PD）方法以及光滑粒子流體動力學（SPH）方法等。 Kolari等[4]提出了一種基于模型更新技術的冰連續(xù)破壞有限元模擬方法，該方法是采用各向異性連續(xù)損傷力學（CDM）模型來預測裂紋的擴展方向，同時也將模型更新技術運用到了有限元模型中裂紋的擴展發(fā)展中。Carney 等[5]建立了冰的現(xiàn)象學破壞模型，并采用顯式動力分析程序LS-DYNA 中的有限元程序?qū)@種冰的塑性敏感性破壞模型予以了求解。Pernas-Sánchez 等[6]建立了基于Drucker-Prager 屈服準則的彈塑性材料模型，得到了冰在高應力率作用下的破碎行為，同時通過采用有限元程序LS-DYNA 中的求解模型，并將拉格朗日、任意拉格朗日?歐拉（ALE）方法和SPH 這3 種方法集成應用到模型研究中，得到了沖擊力對細長冰柱的影響。孔帥等[7]基于高性能DEM 方法，對極區(qū)浮式平臺的冰載荷進行了數(shù)值分析。狄少丞[8]建立了用于模擬海冰與兩單元接觸的離散元模型，研究了冰的單軸壓縮和三點彎曲破壞過程，并與試驗結果進行了對比，顯示吻合較好。薛彥卓等[9]應用近場動力學理論建立了冰材料的彈脆性破壞模型，并針對冰梁的三點彎曲破壞進行了模擬研究，其數(shù)值模擬結果與試驗數(shù)據(jù)的對比顯示其一致性較好。

針對破冰船的破冰阻力預測，可以采用SPH方法進行計算分析。SPH 方法屬于無網(wǎng)格拉格朗日粒子算法，與傳統(tǒng)的網(wǎng)格方法相比，無網(wǎng)格粒子算法的求解不依賴于網(wǎng)格，而是通過一系列的節(jié)點（粒子）來求解整個系統(tǒng)，且這些粒子點之間不需要網(wǎng)格進行連接，可以自由運動，特別適合分析大變形和不連續(xù)的自由表面問題。另外，SPH 方法在模擬固體力學方面因為粒子的拉格朗日特性，粒子在發(fā)生屈服后會自然生成裂紋，因此SPH 方法也非常適合模擬固體大變形以及破壞問題。隨著SPH 方法在流體力學和固體力學領域的快速發(fā)展，該方法也被用于模擬冰動力學、冰破壞以及冰與結構物的相互作用方面。例如，Das[10]在LS-DYNA 軟件中使用SPH 模型對冰梁的四點彎曲破壞進行了模擬，并采用馮米塞斯屈服準則針對冰粒子的破壞失效予以了判別，結果顯示采用該方法一旦冰粒子失效，偏應力分量將直接降為0。Zhang 等[11]使用一種改進的SPH方法，并結合D-P 屈服準則和黏聚力軟化模型的彈塑性本構模型，對冰的彎曲和壓縮破壞特性進行了模擬。Gutfraind 等[12]和Oger 等[13]在SPH 方法中應用基于Mohr-Coulomb 屈服準則的流變學，對在水面上漂浮和在風力作用下移動的碎冰場進行了模擬。Shen 等[14]提出了一個二維河面中冰的漂移和阻塞數(shù)值模型，并采用SPH 方法針對河面冰的漂移運動進行了求解與模擬。Ji 等[15]提出了一種新的海冰動力學黏彈塑性本構模型，其中SPH 方法被用于模擬矩形海盆中冰的運動。此外，Ji 等[16]還發(fā)展了一種混合拉格朗日?歐拉（HLE）海冰動力學方法，其中海冰覆蓋層用具有自身厚度和密集度的SPH 模型表示。Pan 等[17]提出了一個新的SPH 非牛頓模型，用于研究冰蓋和冰架的耦合動力學問題。喬岳[18]應用SPH 方法建立了船?冰?水耦合作用數(shù)值模型，并對冰阻力預報方法予以了分析，不過其主要關注的還是冰阻力的預報方法，對于冰?水耦合模型也只是采用了簡單的排斥力模型，未考慮碎冰區(qū)以及有波浪作用下船?冰?波的相互作用情況。隨后，Zhang 等[19-20]應用SPH 方法建立了船?冰?波耦合數(shù)值模型，模擬研究了冰的破壞過程和波浪?冰的相互作用過程，并預測了船體所受的冰載荷。以上一系列研究表明，使用SPH 方法研究破冰阻力可行且有效。

根據(jù)SPH 方法的無網(wǎng)格特性以及其在流體問題、固體材料大變形及破壞問題，以及流?固耦合問題上的適用性與優(yōu)勢，本文將基于SPH 方法建立船?冰?水耦合數(shù)值模型，并對黃河破冰船的冰阻力進行預報。該方法的構建包括冰的彈塑性本構方程、船?冰?水數(shù)值模型、破碎冰環(huán)境以及與黃河河道中流場的作用。通過與數(shù)值計算結果的對比驗證，對黃河破冰船冰阻力相關因素的影響規(guī)律進行分析，以為破冰船的設計和操縱提供指導。

1 數(shù)學模型的建立

1.1 SPH 的基本方程及相關數(shù)值方法

流體動力學的基本控制方程包括質(zhì)量守恒方程和動量守恒方程，具體如下：

式中：D/Dt為物質(zhì)導數(shù)；α，β 為x，y，z方向上的笛卡爾分量；ρ 為粒子密度；ν 為粒子速度；g為重力加速度；σαβ為應力張量，xβ為位置坐標在β 方向上分量；t為時間。本文中的連續(xù)性方程引入了Antuono 等[21]提出的一種人工擴散項，主要用于消除流體模擬中壓力場不合理的高頻振蕩。為了提高數(shù)值計算的穩(wěn)定性，減小不穩(wěn)定振蕩，引入了Monagham 型人工黏度耗散項。固體相的質(zhì)量守恒方程粒子可近似表達式為

式中： ρi為具有速度分量vi的粒子i的密度； ρj，mj分別為具有速度分量vj的粒子j的密度和質(zhì)量；Wij為SPH 方法中的光滑核函數(shù)。冰粒子動量方程的離散形式為

1.2 冰的破壞模型

為了模擬冰的破壞過程，本文將彈塑性本構模型應用到了SPH 方法中。結合Drucker-Prager屈服準則，得到具有非關聯(lián)流動規(guī)律的冰模型應力應變方程的應力應變關系為：

在破壞過程中，還需要對應力張量進行塑性修正，需要采用黏聚力軟化方法降低塑性流動階段冰的黏聚力，最終實現(xiàn)對冰破壞過程的模擬。由于冰的破壞模型不是本文的研究重點，故在此不予展開，具體的過程可參見文獻[20]。

1.3 船?冰?水耦合模型

在冰?水的SPH 耦合計算模型中，將破冰船當作流體相和冰相的固壁邊界來處理，采用虛粒子邊界方法。該方法是先對破冰船邊界建模，然后劃分網(wǎng)格，隨后在這些網(wǎng)格節(jié)點上布置SPH 粒子以實現(xiàn)對固壁邊界的布置。這些船體邊界粒子作為邊界虛粒子，將分別參與到冰相和流體相的求解中。在求解冰相控制方程時，冰粒子計算域中的船體固壁邊界粒子可以作為虛粒子來施加邊界條件。這些船體固壁邊界粒子將參與冰相的連續(xù)性和動量方程求解，這樣，船體邊界附近的冰粒子將滿足應力和速度的連續(xù)性。在求解流體控制方程時，流體粒子計算域中的船體固壁粒子也將充當流體相的邊界虛粒子。這些固壁粒子同樣將被計入流體相的連續(xù)性和動量守恒方程求解中，用來滿足流體相的壓力和速度邊界條件。

由于船體邊界粒子的速度與船體速度相同，當這些船體邊界粒子作為虛粒子來近似船舶邊界與冰之間的耦合作用時，需要從其鄰域內(nèi)的冰粒子中插值計算出這些邊界虛粒子的應力，用以準確近似邊界附近冰粒子的應力梯度。然后，再從流體鄰域內(nèi)的體粒子中插值計算得到這些固壁邊界虛粒子的壓力，用以近似船體邊界附近流體粒子的壓力梯度。該計算過程可參見圖1。圖中：i，f分別為第i號和第f號粒子；kh為粒子支持域的影響范圍，其中k= 2.0，取1 倍粒子尺寸，h為光滑長度。

圖1 船?冰?水耦合模型示意圖Fig. 1 Schematic diagram of ship-ice-water coupling model

當冰層發(fā)生破壞時，破壞位置處冰粒子的密度值變化較大，為了保證密度的連續(xù)性，在船?冰耦合求解中，船體邊界粒子的密度根據(jù)其領域內(nèi)冰粒子的密度插值計算得到。同理，在船?水耦合求解中，船體邊界粒子的密度也可以由其鄰域內(nèi)流體粒子的密度插值得到。此外，還可通過忽略船體邊界虛粒子與周圍鄰域內(nèi)冰粒子和流體粒子的黏性作用，來施加可滑移邊界條件。通過上述方法，可實現(xiàn)流體相、冰相與固壁邊界的耦合。這種邊界處理方法的主要優(yōu)點是計算過程簡單，不需要計算船體邊界和耦合界面復雜的幾何信息。船體表面的邊界粒子應足夠密，以防止冰粒子和流體粒子穿透船體壁面邊界[19]。需要指出的是，由于直接采用SPH 方法進行大尺度數(shù)值模擬計算效率較低，為了提高計算效率，在上述耦合模型中只考慮施加了流體對冰的作用力，忽略了冰對水的作用。此外，數(shù)值模型中也未考慮船體的運動響應或是力學響應。

另外需要指出的是，在數(shù)值模型中均是采用虛擬粒子法來處理計算域的邊界[20]。例如，在后文的數(shù)值模擬中，對于平整冰的模擬計算，需沿著平整冰面布置相應的邊界虛粒子；而在黃河河道碎冰條件的數(shù)值模擬中，則是沿著河道邊界布置對應的邊界虛粒子。

2 平整冰情況下的數(shù)值模擬與驗證

本節(jié)將采用SPH 耦合計算模型，針對使用區(qū)域黃河上冰面的特點，先給出平整冰情況下的冰阻力計算結果，用以為實際河道中破冰船冰阻力的預測提供參考。

2.1 船模模型建立

所研究的輕型黃河破冰船的總布置圖如圖2（a）所示，其三維模型如圖2（b）所示，主要參數(shù)如表1 所示。

圖2 28 m 黃河破冰船模型示意圖Fig. 2 Schematic diagram of a 28 m ice breaker model operating on the Yellow River

表1 28 m 黃河破冰船主要參數(shù)Table 1 Main parameters of a 28 m ice breaker operating onthe Yellow River

為了采用SPH 方法進行模擬，首先在船體表面建立均勻分布的固壁粒子來表示船體邊界的作用，然后通過求解冰粒子對船體固壁邊界粒子的作用力，得到船?冰耦合過程中船體所受的冰載荷。當船體邊界粒子作為冰粒子的邊界虛粒子時，可以利用相鄰冰粒子的應力張量的體積積分來估算冰粒子對船體邊界粒子的作用力。冰粒子對船體邊界粒子的作用力可以表示為

選取2015年5月～2017年5月我院ICU接受機械通氣治療的新生兒140例作為研究對象，其中，男79例，女61例，平均胎齡（35.36±2.77）d，平均出生體重（2014.05±160.77）g，呼吸窘迫綜合征57例，病理性黃疸49例，高膽紅素血癥21例，缺氧缺血性腦病13例。將其隨機分為研究組與對照組，各70例，兩組一般資料比較，差異無統(tǒng)計學意義（P＞0.05）。

2.2 平整冰面冰阻力的模擬與驗證

在給出SPH 方法的計算結果之前，將引入Lindqvist 經(jīng)驗公式對所計算的冰阻力進行驗證?；诒枇εc速度呈線性關系的假設，Lindqvist公式將冰阻力分為3 個部分進行計算，包括擠壓阻力、浸沒阻力和彎曲破壞阻力[23]，相關公式可參見文獻[23]。Lindqvist 計算公式的主要參數(shù)如表2 所示。

表2 Lindqvist 計算公式的主要參數(shù)Table 2 Main parameters of Lindqvist formula

表3 V = 4 kn，hi = 0.3 m 時不同彎曲強度情況下的破冰阻力預報精度比較Table 3 Comparison of prediction accuracy of ice breaking resistance with different bending strength when V = 4 kn and hi = 0.3 m

圖3 V = 4 kn，hi = 0.3 m 時不同彎曲強度下破冰總阻力的時間歷程比較Fig. 3 Time histories comparison of ice breaking resistance with different bending strengths when V = 4 kn and hi = 0.3 m

為了研究不同冰厚情況下的破冰阻力，圖4給出SPH 方法計算的不同冰厚和不同彎曲強度下冰厚對冰阻力的影響。由圖4 可以看出，冰阻力的變化與彎曲強度大致成線性增長，且冰阻力隨冰厚的變化也較為明顯，當冰厚增大1 倍，其他條件不變時，冰阻力將增大2.5 倍左右，這充分說明冰厚是影響船舶破冰阻力的最重要因素。圖5（圖中， εˉP為應變率幅值）給出了不同冰厚情況下當t= 20 s 時平整冰情況下的破冰形狀?？傮w來看，破冰船在破冰過程中在船艏前會形成周向裂紋，隨后，冰層會遵循相同的破裂模式，這與圖3中冰阻力呈周期性變化的特點是相對應的。此外，在船體的左舷和右舷也發(fā)生了一些局部冰被擠壓破碎的現(xiàn)象，這與實際破冰船在平整冰中破冰過程中的類似現(xiàn)象相符[24]。另外，在船后形成的水道中也存在一些碎浮冰。結果表明，破冰過程并不總是兩側對稱的。由圖5 可以發(fā)現(xiàn)，當冰厚較小時，船艏附近冰的裂紋擴展不明顯，破碎后的塊狀浮冰邊界不清晰；但隨著冰厚的增加，船艏附近冰的破碎程度增加，在出現(xiàn)斷裂后有清晰的塊狀浮冰。這說明采用SPH 方法不僅可以預報冰阻力，還可以預報破碎后的冰裂紋情況。

圖4 不同彎曲強度下冰厚對冰阻力的影響Fig. 4 Influence of ice thickness on ice breaking resistance at different bending strength

圖5 不同冰厚情況下t = 20.0 s 時的破冰形狀Fig. 5 Ice breaking shapes at different ice thickness when t = 20 s

3 河道中冰阻力模擬與分析

為了研究實際河道中破冰船的冰阻力，本節(jié)將對河道中碎冰區(qū)的船舶冰阻力進行分析。首先，根據(jù)黃河某段河道的衛(wèi)星圖片（圖6（a））繪制出河道的邊緣模型（圖6（b））；隨后，為了建立三維河道的剖面圖形，先采用Solidworks 軟件建立河流出口與入口這兩處河道的斷面，然后通過河道兩側地形進行誘導拉伸，得到每個不同位置處的剖面形狀。在模擬河道中冰的堆積與運動情況時，需要有關河道流場流速方面的分布信息。由于直接采用SPH 方法對大尺度的河道流速進行求解有困難，本文采用OpenFOAM 中的icoFoam求解器予以求解。用該求解器求解非穩(wěn)態(tài)不可壓縮N-S 方程，離散格式采用有限體積法。河道網(wǎng)格采用blockMesh 工具生成六面體結構網(wǎng)格，使用snappyHexMesh 工具讀取河道表面stl 文件后畫出河道輪廓并細化網(wǎng)格。隨后，在河道進口處給定流速U=0.5，0.7，1.0 m/s，即可獲得給定流速對應的河道流場分布。

圖6 河道模型的建立Fig. 6 Setting up of the river channel model

建立好河道邊界后，需在河道區(qū)域內(nèi)建立浮冰的SPH 粒子模型。該碎冰模型可采用Voronoi圖來建立，其建模原理如下：首先，根據(jù)碎冰航道的面積以及碎冰密集度和浮冰的平均尺寸，計算出該碎冰區(qū)域內(nèi)浮冰塊的數(shù)量，其中碎冰密集度根據(jù)實際的碎冰覆蓋面積與相應冰區(qū)總面積之比得到，浮冰平均尺寸為碎冰覆蓋總面積與碎冰數(shù)量之比；然后，在所需計算的河道模型區(qū)域內(nèi)生成相同數(shù)量的隨機樣點，根據(jù)隨機點生成Voronoi多邊形，以保證多邊形內(nèi)任何位置點到該多邊形隨機樣點的距離最近；最后，根據(jù)選定的碎冰密集度對Voronoi 多邊形進行收縮，以確定SPH 粒子的填充區(qū)域。根據(jù)該方法生成的河道碎冰區(qū)計算模型如圖7 所示。

圖7 碎冰建模示意圖Fig. 7 Setting up of broken ice model

圖8 所示為碎冰密集度C= 70%時破冰船在河道中作業(yè)時的模擬結果。碎冰計算模型中，浮冰的平均尺寸為12.1 m2。該船從河道入口出發(fā)，行駛軌跡與河道主航道重合，航速保持V= 4 kn不變。由圖可見，采用SPH 方法能夠?qū)λ楸诤拥乐械钠七M行有效模擬，而且在破冰船的作用下，船行駛后尾部有明顯的開闊水道。為了對不同碎冰密集度下的破冰阻力進行分析，圖9 給出了航速V= 4 kn、冰厚hi= 0.3 m 時不同碎冰密集度下破冰船冰阻力的時間歷程比較。由圖9 可以看出，當碎冰密集度較小（C= 50%）時，浮冰對船體破冰阻力的作用比較小，當碎冰密集度增大到C= 90%時，浮冰對船體破冰阻力的作用明顯增大。另由圖9 所示的時歷曲線還可以發(fā)現(xiàn)，當碎冰密集度較小時，冰阻力的作用時間比較短，且每次作用時冰阻力的最大幅值也較小。表4 給出了不同碎冰密集度下的冰阻力平均值?？梢婋S著碎冰密集度的增加，冰阻力的作用時間逐漸變得連續(xù)，冰阻力的平均幅值也逐漸增大。碎冰密集度對破冰船破冰阻力作用時間的分布有著明顯的影響。

圖8 碎冰密集度為70%時破冰船在河道中作業(yè)的模擬結果Fig. 8 Numerical results of icebreaker working in the river when the concentration of broken ice C = 70%

圖9 V = 4 kn，hi = 0.3 m 時不同碎冰密集度下破冰船總阻力的時間歷程比較Fig. 9 Time histories comparison of breaking ice resistance with different concentrations of broken ice when V = 4 kn and hi = 0.3 m

表4 不同碎冰密集度情況下SPH 計算的冰阻力平均值Table 4 Average ice resistance by SPH with difference concentrations of broken ice

為了研究碎冰情況下不同冰厚對破冰船冰阻力的影響。圖10 給出V= 4 kn、碎冰冰密集度C= 70%時，采用SPH 方法計算得到的不同冰厚下破冰船冰阻力的時間歷程比較。圖中數(shù)值結果對應的模擬時間歷程為120 s。由圖可見，冰厚對冰阻力的曲線形狀變化影響不明顯，但對冰阻力的幅值變化影響顯著。當hi= 0.3 m 時，冰阻力的最大幅值為2.73×105N；當hi= 0.4 m 時，冰阻力的最大幅值為3.11×105N；當hi= 0.5 m 時，冰阻力的最大幅值為4.93×105N，當hi= 0.6 m 時，冰阻力的最大幅值增大為5.81×105N。由圖可知最大冰阻力的增幅是冰層厚度增幅的2.7 倍，說明冰厚是影響破冰阻力的主要因素。表5 示出了不同冰厚對冰阻力的影響。當hi= 0.3 m 時，冰阻力的平均值為31.88 kN；當hi= 0.4 m 時，冰阻力的平均值為33.53 kN；當hi= 0.5 m 時，冰阻力的平均值為48.48 kN；當hi= 0.6 m 時，冰阻力的平均值為61.03 kN?？梢姳枇Φ钠骄迪鄬τ谧畲蠓到档土? 個數(shù)量級，此時，平均冰阻力的增幅是冰層厚度增幅的1.91 倍。

表5 不同冰厚情況下SPH 計算冰阻力平均值Table 5 Average ice resistance by SPH with difference ice thicknesses

圖10 V = 4 kn，C = 70%時不同冰厚下破冰船總阻力的時間歷程比較Fig. 10 Time histories comparison of breaking ice resistance with different ice thicknesses when V =4 kn and C = 70%

為了驗證SPH 結果的可靠性，采用文獻[25]中有關碎冰中的船舶冰阻力進行了驗證，計算公式如下：

表5 關于不同冰厚條件下SPH 計算結果與經(jīng)驗公式的比較說明，采用SPH 方法模擬河道中碎冰的冰阻力具有一定的合理性。

4 結 語

本文首先對SPH 方法進行了拓展，通過結合冰的彈塑性本構方程和河道的建模方式，完成了船?冰?水耦合數(shù)值模型的建立；然后，采用該數(shù)值模型對平整冰面中黃河破冰船的冰阻力進行預報，并與Lindqvist 經(jīng)驗公式進行比較驗證了其有效性；接著，將該數(shù)值模型拓展應用到了黃河河道中破冰船的冰阻力預報中；最后，通過對不同情況下冰阻力計算結果的比較，分析了黃河破冰船冰阻力的規(guī)律和特性。主要得到如下結論：

1） 采用SPH 方法能夠較好地預報黃河破冰船在平整冰面情況下的破冰阻力，通過與Lindqvist經(jīng)驗公式的對比，顯示采用SPH 方法得到的冰阻力誤差結果小于17.6%。

2） 可以采用Voronoi 圖的方法建立破碎冰面的SPH 粒子模型，該模型可以通過對碎冰密集度的調(diào)整來實現(xiàn)對不同碎冰冰面初始時刻的劃分，能為河道中碎冰的運動及其對船舶的碰撞作用提供研究基礎。

3） 借助icoFoam 求解器求解的河道流場，采用SPH 方法能夠?qū)崿F(xiàn)對真實河道中碎冰運動和冰阻力的模擬計算。通過對不同碎冰密集度情況下冰阻力的對比研究，發(fā)現(xiàn)在該河道中碎冰條件下冰阻力的作用方式與平整冰面不同，河道中破冰阻力的作用時間歷程變化非常劇烈，且作用時間較短。此外，冰厚也是影響冰阻力最重要的一個因素。

以上規(guī)律和特性可為黃河破冰船的實際操作提供參考。