基于DEM-SPH 流固耦合的冰區船舶快速性預報

吳 捷 劉 璐 田于逵 季順迎

（1. 大連理工大學 工業裝備結構分析優化與CAE 軟件全國重點實驗室 大連 116024； 2. 大連理工大學 船舶工程學院大連 116024； 3.中國船舶科學研究中心 無錫 214082； 4. 深海技術科學太湖實驗室 無錫 214082）

0 引 言

極地船舶在冰區破冰航行過程中需要考慮船舶快速性能，而極地海冰和海水組成的流固混合場對船體結構的快速性能具有重要影響[1]。船舶的快速性一般研究航行阻力和推進力平衡狀態下的主機功率、螺旋槳轉速和扭矩等參數，同時考慮船舶航行過程中的船體結構、推進裝置與流體的流固耦合作用[2]。對于冰區航行的極地船舶，還需考慮海冰與結構的相互作用，主要包括海冰與船體結構，以及在船體周圍伴流場作用下，破碎海冰同螺旋槳和推進器相關結構產生的碰撞作用等[3]。因此，極地船舶的快速性研究需要充分考慮船體和推進器結構在冰水混合場中受到的作用力，對船舶在特定航速下的阻力和推力等力學參數進行預報分析，從而獲得極地船舶快速性的相關參數。

結構在浮碎冰區航行過程中主要與碎冰發生碰撞、摩擦等作用，會受到碎冰密集度、形狀與尺寸等參數的影響[4]。由于結構在浮碎冰區的運動會引起自由液面興波，因此當碎冰密集度較低時，結構周圍碎冰會在波浪作用下提前被排開，導致結構與碎冰的作用次數減少，進而降低了結構的航行阻力；當碎冰密集度較高時，結構首部處的碎冰會在周圍碎冰的持續作用下與結構同步運動，形成類似附連水效應，從而增大了結構的航行阻力。

在破冰船開辟的冰區水道中，一般存在尺寸較小的海冰，其會在螺旋槳附近流場的作用下，同槳葉與局部結構發生高速碰撞、切削作用[5]。然而，針對碎冰受到流場擾動并與結構相互作用的物理過程的研究目前并不充分，該場景下局部結構的動力響應及運動受到的海冰影響也不明確。此外，高速運動的結構會對海水形成一定沖擊作用，所造成的強烈流場波動會提前導致海冰的裂紋擴展、斷裂與破碎，從而影響海冰與結構的相互作用過程，進而受到海水的流固耦合作用。目前還缺少流場影響下的海冰與結構相互作用試驗研究，相關的數值模型研究也極其匱乏。

數值模擬方法是研究極地船舶與海洋工程結構同海冰相互作用、冰載荷、冰激結構響應的重要手段。近年來，有限元[6]、離散元[7]、近場動力學[8]等方法在海冰破壞過程模擬中應用廣泛。離散元方法（discrete element method，DEM）對脆性材料破壞過程的模擬具有較高的計算效率，其對材料破壞過程中的裂紋擴展、破碎特性等均能進行出色的模擬[9-10]，故在海冰與結構相互作用的數值模擬中受到了廣泛關注。結合結構的有限元方法還可建立海冰與結構相互作用的離散元-有限元耦合模型，在船舶、風機等結構的冰激響應和振動特性分析中取得了良好效果[11-12]。另外，在船體結構與海冰的相互作用數值模型中，還可進一步結合海水的計算流體力學方法，分析流場影響下的海冰與結構相互作用過程[13-14]。

然而，各方法之間由于收斂時間步長不同，會導致計算時間尺度差異問題，在不同時間步長上的參數傳遞效率較低；不同形態海冰單元與流體存在復雜形狀的流固界面，傳統方法計算海冰與流體間耦合作用的效率較低。因此，海冰與結構的相互作用研究中需要考慮海水與海冰、結構的流固耦合作用，建立結構在冰水混合場中航行運動的流固耦合模型，為冰載荷和冰激結構響應計算提供重要的方法基礎。

基于拉格朗日粒子的光滑粒子流體動力學（smoothed particle hydrodynamics，SPH）方法在自由液面、大變形模擬中具有明顯優勢，在海洋工程、化工攪拌等領域中應用較為廣泛[15]。在DEM與SPH 的耦合模型中，不僅需重點考慮固體單元截斷流體粒子支持域造成相關物理量計算不準確的問題，還需考慮不同時間步長引起的計算參數傳遞問題[16]。固體粒子模型將邊界抽象為由若干層固定不動的粒子所組成的粒子墻，是SPH 中使用最為廣泛的邊界構造方法[17-18]。但由于固體粒子需消耗大量的計算資源，故采用粒子構造邊界的方法難以適應不同形態固體顆粒的模擬，而通過光滑函數的積分正則化修正粒子密度和速度，可有效構造復雜的固體邊界并計算其與流體粒子的相互作用[19]。采用若干層固體粒子構造固體邊界，并假定固體粒子的密度等物理量與流體粒子相同，繼而將流體粒子的光滑函數導數擬合為與固體邊界距離的函數，在實際的SPH 模擬中即可直接采用該擬合函數計算流體粒子受到的邊界作用[20]。以上SPH 方法中的流固邊界計算方法為流體粒子與固體邊界的快速計算研究奠定了良好的基礎，據此可進一步發展海冰與海水的DEM-SPH 耦合模型。

本文分別采用離散元方法和光滑粒子流體動力學方法模擬海冰和海水，并建立冰水耦合的DEM-SPH 模型，對剛性船體和螺旋槳結構進行簡化建模，建立合適的流體便捷模型，模擬并分析船舶航行過程中所收到的總阻力與螺旋槳推力情況，反推出自航所需的螺旋槳轉速，進而分析船體在冰區航行過程中的快速性。

1 基于DEM-SPH 的冰水耦合模型

本文采用基于平行黏結模型的離散元方法模擬海冰與結構作用過程中的斷裂和破碎，并考慮海冰破裂后碎冰塊之間的接觸、碰撞和摩擦作用。采用SPH 方法模擬海水，同時考慮固定計算域邊界處的速度和壓力邊界條件，實現固定流體域上的海水流動模擬。采用SPH 粒子與固定粒子邊界相對運動的擬合項，直接計算固體與流體之間的相互作用力，建立海冰與海水耦合的DEM-SPH 模型。

1.1 海冰離散元方法

海冰離散元方法基于海水離散元模型（圖1），其中采用密排六方形式的球體單元表征海冰，球體單元之間可傳遞黏結力、接觸力等，如圖1（a）所示[21]。黏結力表示連續海冰內部的作用情況，通過對黏結力施加一定的失效準則，連續海冰可發生斷裂、破碎等動力學過程；接觸力則表示海冰發生破壞后由接觸、碰撞產生的相互作用。

圖1 海冰離散元模型

平行黏結模型是離散元方法中實現材料破壞模擬的經典模型，其將2 個球體膠結在一起，不僅可傳遞力，還可傳遞力矩。如圖1（b）所示，平行黏結模型將相鄰單元之間的黏結作用簡化為梁或黏結圓盤[22]，采用梁模型計算單元之間的法向和切向作用力，如式（1）所示：

式中：Fb為單元之間的黏結力，N；Mb為單元之間的黏結力矩，N·m；Fbn和Fbs分別為單元之間的法向和切向黏結力，N；Mbn和Mbs分別為單元之間的法向和切向黏結力矩，N·m。

此處可采用增量的形式計算黏結力和力矩，如式（2）所示：

式中：ΔUn和ΔUs分別為單元間相對位移的法向和切向分量，m；Δθn和Δθs分別為相對轉動位移的法向和切向分量，m。各參數表達式見式（3）：

式中：u為平動速度，m/s；ω為轉動速度，rad/s；Δt為時間步長，s。

根據三維梁理論，單元之間的最大拉伸σmax和剪切應力τmax見式（4）：

式中：A為黏結圓盤的橫截面積（A=πR2），m2；J為黏結圓盤的極慣性矩（J=πR4/2），m4；I為黏結圓盤的慣性矩（I=πR4/4），m4。

在單元黏結作用過程中，若黏結圓盤的最大拉應力超出其法向黏結強度，或最大剪應力超出其切向黏結強度，即σmax＞或τmax＞，則該黏結作用失效，從而可模擬海冰材料內部的斷裂和破壞過程。和滿足Mohr-Coulomb 準則，如圖2 所示。同時，在基于平行黏結模型的離散元模擬中，和受到單元尺寸的影響，故需通過材料的典型力學試驗進行校核，如三點彎曲和單軸壓縮試驗。

圖2 法向和切向黏結強度滿足Mohr-Coulomb 準則

通過海冰材料的三點彎曲和單軸壓縮數值試驗，可建立海冰離散元試樣的材料強度與黏結強度的關系式[23]，見式（5）：

式中：σf為海冰的彎曲強度，MPa；hi為海冰厚度，mm；R為球體單元的半徑，mm。

海冰發生破壞后，單元之間的黏結作用轉變為接觸摩擦作用，這里采用球體單元的線性接觸模型計算接觸力。該接觸力可分為由剛度導致的彈性力Fnce和阻尼導致的黏性力Fncv。法向的接觸力Fnc可表示為式（6）：

式中：δn為法向重疊量，m；Cn為法向阻尼系數；δn為單元間接觸點處的法向相對速度，m/s；Kn為兩相鄰單元間的法向剛度（單位是N/m），表達式見（7）：

式中：E為海冰的彈性模量，Pa；i和j表示兩相鄰單元的編號。

法向阻尼系數Cn可表示為式（8）：

式中：mij為等效質量，mij=mimj/(mi+mj)，kg；ξn為無量綱阻尼系數，為恢復系數。

切向接觸力Fcs可表示為式（9）：

式中：Ks為兩相鄰單元間的切向剛度；Cs為切向阻尼系數；t和■分別為切向接觸重疊量和相對切向速度的單位向量；δs為切向的接觸重疊量，可根據相對切向速度和逐步積分獲得；sδ˙ 為兩單元接觸點處的相對切向速度，m/s；μ為滑動摩擦系數。Ks和Cs分別按照Kn和Cn的一定比例計算，見式（10）：

式中：α和β均為0.1[24]。

單元間的接觸力Fc和由切向接觸力引起的接觸力矩Mc，可表示為式（11）：

式中：n為接觸的單位法向量；s為切向上的單位矢量，

考慮每個單元受到的黏結和接觸作用，單元的運動方程可表示為式（12）：

式中：Is為球體的轉動慣量，Is= 0.4mR2；Fh為流體給海冰單元的作用力，N；Mh為流體給海冰單元的作用力矩，N·m。

1.2 基于SPH的計算流體動力學方法

相比于基于網格的計算流體力學方法，SPH 方法更易于處理自由表面問題以及冰-水或船-水間的固液交界面。因為SPH 方法本身是一種無網格的粒子法，其不需要像傳統的網格方法一般引入體積分數來表征自由液面或固液交界面，在計算船冰水耦合問題上相比于傳統計算流體力學動力學方法具有一定的優勢。同時，相對于傳統的弱可壓縮格式以及隱式不可壓縮格式的SPH 方法，本文采用顯式不可壓縮格式的SPH 方法求解流體壓力，在計算精度與計算效率上找到了較好的平衡點，更易于進行三維大規模數值模擬[25-26]。

SPH 方法可以將連續的場函數進行拉格朗日格式離散。如圖3 所示，粒子i處的任意場函數A（ri）見式（13）：

圖3 SPH 核函數示意圖

式中：r為粒子的坐標；N為粒子相鄰粒子的總數；m為粒子的質量，kg；ρj為該粒子的密度，kg/m3；rij為粒子間的相對坐標，可寫為rij=ri-rj；h為光滑長度；W(rij,h)為核函數，其是與粒子間距離有關的偶函數[27]。

這里采用三次樣條曲線方程作為核函數，見式（14）所示：

式中：q為正則化粒子間距離，q=rij/h。

粒子i處的任意場函數梯度 ?·A(ri)可表示為式（15）：

式中：rij為粒子間的距離（rij=|rij|=|ri-rj|），m；?iWij為核函數梯度。

在不可壓縮格式下的流體拉格朗日型Navier-Stokes 方程可表示為式（17）和下頁式（18）：

式中：u為流體的速度，m/s；P為流體的壓力，Pa；v為流體運動黏度，對于海水，可以取v=10-6m2/s；g為重力加速度，m/s2。

由于不可壓縮流體密度保持恒定（即 ?·u= 0），Navier-Stokes 方程可以近似表示為式（19）和式（20）：

式中：ρi=ρj=ρ0，ρ0為不可壓縮流體的密度，kg/m3；黏度項 (v?2u) 可表示為式（21）[28]：

式中：uij為粒子間的相對速度（uij=ui-uj），m/s；η為保證分母不為0 的較小值，η=0.1h。

采用兩步預測修正的方式可獲得SPH 粒子的壓力泊松方程[29]。在預測步，計算只受黏度項及重力項影響下的粒子速度，粒子經過預測步后的速度u*可表示為式（22）：

式中：Δt為運動積分的時間步長，s；un為粒子在當前時間步的速度，m/s；粒子下一時間步的速度un+1可由預測步狀態進行壓力梯度修正得到，見式（23）：

基于速度散度為0 的條件（即 ?·un+1= 0），可將上式轉化為式（24）所示壓力泊松方程：

為了提高計算效率，采用可顯式計算的近似解計算壓力泊松方程[30]，粒子下一時間步的壓力可表示為式（25）：

式中：Aij可表達為式（26），Bi可表達為式（27）：

式中：u*ij為2 個SPH 粒子在預測步中的相對速度（u*ij=u*i-u*j），m/s。

使用上式計算所得的粒子壓力Pni+1帶入控制方程中，可以得到下一時間步粒子的速度，下一時間步粒子的位置rn+1可表示為式（28）：

為施加固定流體域上的邊界條件，采用設置虛粒子層的方式構建緩沖區，從而可設置流體的出入口邊界[31-32]。流體的出入口邊界模型大致分為壓力出入口邊界、速度出入口邊界、開放邊界。3 種邊界均采用虛粒子方法在邊界處形成緩沖區，避免流體計算域內的粒子核函數計算域缺失，虛粒子緩沖區如圖4 所示。

圖4 邊界緩沖區示意圖

在本模型中，出入口邊界的幾何法向均指向真實計算域內，并在負法向上添加虛粒子以構建緩沖區。緩沖區內虛粒子初始間距與真實粒子一致，且緩沖區寬度需大于2 倍光滑長度[33]，即最少也需2層虛粒子。通過緩沖區內的虛粒子與真實粒子的位置和物理量交換，實現計算域邊界處速度、壓力的邊界條件。

1.3 DEM-SPH流固耦合模型

為得到擴展多面體顆粒和SPH 粒子之間的耦合作用，本文將單元視為一種剛性固體邊界，并通過擬合的方法得到兩者之間的耦合作用。對于SPH 粒子核函數計算域被固體邊界截斷的情況，在邊界外緊密排列了足夠多的固定邊界粒子，進而模擬邊界對SPH 粒子的作用[20,34]，如下頁圖5 所示。固定粒子的間距與SPH 初始粒子間距相同，固體粒子質量和密度也與SPH 粒子相同，速度為0。邊界粒子的壓力見式（29）：

圖5 SPH 剛性邊界模型

式中：PB為邊界粒子的壓力，Pa；d為SPH 粒子到邊界的距離，m；nB為邊界面法向。

此狀態下，邊界粒子對SPH 粒子的作用力可表示為式（30）：

式中：rib為SPH 粒子與邊界間距離，m；uib為SPH 粒子與邊界間的相對速度，m/s。

式（30）右側第1 項可簡化為邊界對SPH 粒子的法向作用力，見式（31）：

式中：F?W可表示為式（32）：

式（30）右側第2 項可簡化為邊界對SPH 粒子的切向作用力，見式（33）：

式中：F?W、F?W/r只與h/Δx及SPH 粒子到邊界的正則化距離q有關。此處：q=d/h；h/Δx=1.25。F?W、F?W/r的關系曲線如圖6 所示。

圖6 F? W和 F? W/r 與q 的關系曲線

對F?W、F?W/r和q之間的關系進行多項式擬合[35]，見式（35）、式（36）：

SPH 方法中的粒子運動均為顯式積分，因此計算過程中可看作流固界面上存在一對排斥力分別作用于流體粒子和固體界面，固體粒子所受SPH 粒子作用力之和即視為固體單元受到的流體作用力。

2 冰區船舶的快速性預報分析

由于冰區船舶的快速性需要同時考慮海冰和海水的耦合作用，采用流固耦合算法實現難度較大，目前尚未有相關的研究報道。為此，本文基于海冰離散元方法建立的DEM-SPH 耦合方法的優勢和特點，形成了冰區船舶快速性預報的方法流程。該模型分別對船舶在冰區的航行阻力和推進力進行模擬：在船舶航行阻力的計算中，模擬船舶在一定航速下的航行阻力；在推進力的計算中，模擬螺旋槳在不同轉速下的推進力。通過擬合的方式匹配航行阻力和推進力，從而預報船舶在特定航速下實現自航所需的螺旋槳轉速。

2.1 船舶快速性分析的數值計算模型

為分析船體結構在不同航行工況下的自航螺旋槳轉速，本文在模擬中使船體保持不動，約束其六自由度運動，并給定流場速度以模擬不同螺旋槳轉速下船體結構所受的強制力，進而擬合出強制力為0 時的螺旋槳轉速，該轉速即為所需自航條件下的螺旋槳轉速。為了考慮船冰間的相互作用，文中采用DEM-SPH 耦合的方法模擬船體結構-海冰-海水三者間的相互作用。此時，若仍同時考慮船體冰阻力和螺旋槳推力的計算，則所需SPH 計算精度較高、計算規模過大，目前較難實現，故將船體快速性預報分為2 個模擬算例。

（1）船體航行阻力模擬

首先計算船體冰區航行過程中的總阻力，其中包含了冰阻力和水阻力。此時流場包裹全船，但計算精度相對較低，在計算過程中螺旋槳并不旋轉，具體計算設置如圖7 所示。

圖7 阻力計算模型示意圖

計算過程中，船體結構保持不動，通過前部的速度入口及尾部的壓力出口得到穩定的流場，其中，速度入口處流體速度等于該工況下的船速；海冰從流體入口處進入計算區域，在進入計算區域之前，其速度保持不變并與流體入口速度相同。

（2）螺旋槳推力模擬

在計算得到船體航行總阻力后，將SPH 計算域縮小，只考慮船體尾部螺旋槳附近的流場計算，其余位置的流場使用恒定流場代替，如圖8 所示。在SPH 計算域外部，海冰受恒定流場的拖曳力和浮力作用，進入SPH 計算域后才改為計算流體與海冰單元間的相互作用。同時，考慮速度入口與壓力出口并不對稱，將底部流場邊界設為開放邊界，兩側邊界設為速度入口邊界，其速度大小及方向與該工況下的船速相同。

圖8 推力計算模型示意圖

在該模型計算中，需計算得到不同螺旋槳轉速下螺旋槳受到的推力及推力減額情況。其中，推力減額被認為是由于船體受流場阻力的增加引起的，因此需計算零轉速下船體所受的流場作用力，并與高轉速下所得計算值進行對比，得到對應轉速下的推力減額情況。

船體模型采用DTMB 5415 單體船型，配備雙軸槳推進，螺旋槳旋轉方向如圖9 所示。

圖9 船體模型示意圖

計算工況按照縮尺比1∶30 設計，船長4.733 m、寬0.634 m、吃水0.205 m，海冰及流場計算參數列于表1 中，相關計算參數均為模型尺度下的參數值。

表1 船體冰阻力計算參數

2.2 浮冰區船舶自航下的螺旋槳轉速分析

基于上述DEM-SPH耦合模型及邊界條件模型，建立了船體結構冰區航行阻力計算模型，其中浮冰厚度為0.033 m、尺寸設為1 m × 1 m、密集度為80%，與模型試驗保持一致；航速設置為0.282 m/s，并設置了3 個螺旋槳轉速，分別為100 r/min、200 r/min、400 r/min。船體破冰結果如圖10 所示。

圖10 船體在浮冰區航行過程模擬

兩側浮冰基本沒有發生破碎的情況，中間的浮冰因與船體相互作用大多從中間斷裂，并向船體兩側滑開。由于受到兩側浮冰的限制，船后水道只比船寬稍寬，但其中的碎冰較少。不過，隨著船體結構遠離，兩側破裂的大塊浮冰會再次占據航道，船后遠端的航道則基本重新閉合。

船體在破冰過程中受到的阻力情況見圖11。其中，10 s 之前由于海冰未同船體結構碰撞，故其阻力主要為流體阻力。可以看出，流體阻力遠小于冰阻力，這也與其他學者的模擬結果一致。在進入冰區后，冰載荷呈現較大波動，但總體保持穩定（平均阻力為28.7 N），這說明此時冰阻力主要為船體破冰產生的作用力，破壞后的碎冰對船體的浸沒阻力與摩擦力在該工況中占比較小。

圖11 船體浮冰區航行的阻力時程

不同螺旋槳轉速下的尾部流場模擬結果參見圖12。從圖中可以看出：隨著螺旋槳轉速的增加，在螺旋槳和舵槳之間產生高速流場，經過舵槳后的流速明顯降低；并且由于浮冰破碎后向兩邊排開，因此船體底部沒有海冰下潛，破碎的海冰與螺旋槳之間不會發生碰撞。

圖12 不同螺旋槳轉速ωp 下的尾部流場模擬結果

圖13 所示為2 個螺旋槳所受的總推力以及各個螺旋槳所受的力矩情況。從圖中可以發現：當轉速由200 r/min 提高到400 r/min 后，螺旋槳所受推力與力矩都明顯增加；對比相應轉速下計算得到的流場情況，發現400 r/min 下的螺旋槳后部流場速度也明顯提高。計算得到的流場現象與推力統計現象兩者一致。

圖13 螺旋槳所受推力及扭矩時程

圖14 所示為尾部流場對船體結構的作用力以及計算得到的推力減額。從圖中可以看出，當轉速由200 r/min 增至400 r/min 后，推力減額數值明顯提升。從前文分析可知，螺旋槳附近的流場流速增加最為明顯，由此可以認為該方法得到的推力減額數值，就是由于流場流速增加所導致的船體阻力增加部分。

圖14 不同轉速下的推力減額情況

將計算結果統計后列于表2 中，計算出各個轉速下的強制力情況；再以強制力為x軸、螺旋槳轉速為y軸，對所得數據進行二次擬合，得到的曲線截距即為自航時所需的螺旋槳轉速；最后便可得出該工況下自航所需的螺旋槳轉速為281.52 r/min，如圖15 所示。

表2 浮冰工況中不同螺旋槳轉速的計算結果表

圖15 通過強制力擬合獲得自航轉速

通過更改航速計算3 個不同航速工況，將結果匯總于表3 中；同時也統計了不同航速下所需的螺旋槳轉速，如圖16 所示。

表3 不同航速工況下的計算結果

圖16 不同航速下所需的螺旋槳轉速

從圖中可以發現，隨著航速增加，所需的螺旋槳轉速也出現明顯增加，這主要是因航速增加后，阻力也快速增加所導致。該現象具有一定合理性，在沒有試驗對比的條件下，可初步認為所用方案可行。

2.3 層冰區船舶自航下的螺旋槳轉速分析

本文基于浮冰計算模型建立了船體平整冰區航行阻力計算模型，平整冰厚度與浮冰一致，平整冰區域寬1.33 m、航速設置為0.188 m/s，設置3 個螺旋槳轉速，分別為300 r/min、500 r/min、700 r/min。船體阻力計算結果參見圖17，其中破冰結果如圖17（a）所示。船體破冰后碎冰會向兩側冰層下方運動，因此船尾部航道中浮冰較少。船體航行破冰阻力時程如圖17（b）所示，取60 ～80 s 段計算破冰的平均阻力，該工況下船體平均阻力為75.56 N。

圖17 船體阻力計算結果

圖18 所示為不同螺旋槳轉速下的流場模擬結果。從圖中可以看出，由于冰層破碎后，碎冰大部分被擠壓到兩側平整冰層下方，下潛到船體底部的海冰較少，因此破碎的海冰與螺旋槳之間不會頻繁發生碰撞。

圖18 不同螺旋槳轉速下流場模擬結果

圖19 所示為2 個螺旋槳所受的總推力及所受的力矩情況。從圖中可以發現，該工況中3 個轉速下的推力呈現線性增加趨勢，且對應扭矩的增加也符合線性規律。

圖19 螺旋槳所受推力及扭矩時程

圖20 所示為尾部流場對船體結構的作用力以及計算得到的推力減額值。從圖中可以看出，隨著螺旋槳轉速的增加，該工況中推力減額的變化情況基本呈線性增加趨勢，與其推力的線性增加趨勢一致。

圖20 不同轉速下的推力減額情況

將計算結果統計后列于表4 中，計算各個轉速下的強制力，并采用浮冰工況類似的擬合方式獲得自航時所需的螺旋槳轉速（該工況下所需螺旋槳轉速為501.8 r/min），計算航速分別為0.282 m/s 和0.188 m/s，而后將結果匯總于表5。

表4 層冰工況中不同螺旋槳轉速的計算結果

表5 計算結果匯總表

3 結 語

由于冰區船舶的快速性需要同時考慮海冰和海水的耦合作用，目前尚未有采用冰-水-船三者耦合開展冰區船舶快速性預報的研究報道。本文采用海冰離散元方法和海水光滑粒子流體動力學方法建立了冰水耦合的DEM-SPH 模型，結合DEM-SPH 耦合方法的優勢和特點，形成了冰區船舶快速性預報的方法流程，并將船體和螺旋槳視為剛體模型，形成了阻力-推力預報模型。文中采用DTMB 5415 船模模擬了船舶在浮冰區和層冰區的航行過程，分析了不同航速下船體阻力和螺旋槳推力，將尾部流場對船體的阻力增額作為推力減額值，擬合了船模強制力隨螺旋槳轉速的變化曲線，從而獲得了船舶自航下的螺旋槳轉速。基于上述模型分別針對浮冰區和層冰區進行了自航螺旋槳轉速計算，所得計算結果具有一定合理性。由于目前未能開展相關模型試驗研究，因此后期還需要進一步結合模型試驗結果驗證預報結果的準確性，并優化計算分析流程，從而使本模型可以為極地航行船舶的船機槳匹配工作提供計算參考。