螺旋槳在冰槳銑削下的強度計算分析

常欣, 王錫棟, 王超, 孫盛夏

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

螺旋槳在冰槳銑削下的強度計算分析

常欣, 王錫棟, 王超, 孫盛夏

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

冰區(qū)航行時，相比于水動力，螺旋槳冰銑削工況下的冰載荷對螺旋槳的危害更大，故冰區(qū)航行對螺旋槳的強度要求較高。本文以顯示算法理論為基礎，利用ANSYS/LS-DYNA為工具，在模擬球形冰和剛體的碰撞與相關實驗對比驗證選用的計算方法的正確性后，冰槳銑削工況下的螺旋槳載荷響應進行了模擬研究，計算監(jiān)視了冰槳銑削過程中槳葉的等效應力、位移響應以及冰的破壞過程并得到分析計算結果。計算結果表明：冰槳銑削過程螺旋槳受到的等效應力主要集中在葉梢及導邊區(qū)域，應力值在冰槳接觸后迅速達到峰值；在銑削過程中葉梢受到的的彎曲變形最為嚴重，在冰區(qū)槳設計環(huán)節(jié)應著重關注葉梢部分。

冰； 螺旋槳； 顯示動力分析； 銑削； 強度計算； 數值模擬

隨著北極航道的通航和極地油氣資源的開發(fā)利用，冰區(qū)船舶的重要性與日俱增，而對于冰區(qū)船舶來說，海冰的存在是危害船舶航行安全的一項重要隱患。冰載荷和水動力載荷是螺旋槳與冰相互作用時螺旋槳受到的兩種主要載荷，螺旋槳所受到的水動力載荷遠遠小于冰載荷，槳葉的強度主要取決于冰載荷的大小。因此，研究海冰影響下船舶螺旋槳強度問題對確保船舶冰區(qū)航行條件下的可靠性與安全性有著重要而深遠的意義。國外對冰槳接觸研究有很長的歷史，Veitch在Belyashov和Shpakov實驗模型的基礎上提出了冰槳接觸模型的基礎概念[1]。Veitch的模型可估算任意切削角時做用于葉片上的冰載荷。雖然Veitch假定了一個合適的螺旋槳與冰相互作用的情況，沒有考慮螺旋槳槳葉間的相互作用；而且只有球形冰的形狀被考慮在內，并忽略了在螺旋槳與冰相互作用時冰的形狀和質量的變化。Akinturk 等做了一系列冰池螺旋槳模型試，選用不同的冰型研究了冰阻對螺旋槳水動力性能的影響，并不斷改進數值方法分析槳-冰相互作用，得到的結果與試驗結果吻合較好，研究表明槳-冰相互作用載荷受螺旋槳的幾何參數、進速系數、攻角和槳的銑削深度等因素的影響較明顯[2-6]。為了探討固體冰的破壞，Soininen對莫爾—庫侖破壞準則的滑移線理論進行了研究，碎冰被認為是使用粘性顆粒擠出的壓力分布模型，作用在葉片的總負荷是由每一段的有效載荷相加而得，繼而提出并驗證了計算有效載荷冰槳接觸模型[7]。Wang Jungyong開展了冰槳相互作用下的混合載荷模型試驗研究，數值結果和實驗結果吻合較好[8]。T.J. Huisman等進行了冰槳接觸下的泡沫塑料模型冰模型試驗，測定冰槳接觸過程的銑削力[9]。

國內對冰區(qū)螺旋槳的研究長期以來都處于空白狀態(tài)。何菲菲對使用導管槳和全回轉推進器的內河破冰船進行了推進性能計算[10]。胡志寬等對螺旋槳在冰槳碰撞條件下的強度問題進行了動力學分析計算[11]，研究了冰槳接觸條件下螺旋槳應力的影響因素。郭春雨等就冰級槳的水動力性能搜集整理了國外現有的試驗和數值研究方法及相應的進展情況[12]。

本文運用非線性顯示有限元方法，進行冰槳接觸狀態(tài)計算，獲取冰槳流一體下冰槳銑削總載荷，模擬螺旋槳與海冰的銑削過程，對銑削過程中的接觸力和最大等效應力的變化及槳葉的最大變形進行了研究。

1 數學模型

LS-DYNA為一個以顯式為主，兼顧隱式非線性動力有限元分析的計算程序[13]。其主要用于求解三維非彈性結構的高速碰撞，爆炸沖擊下的大變形動力響應，在航空、汽車領域廣泛應用。對于動力學問題，LS-DYNA的顯示求解方法采用中心差分法求解在時間t時的加速度：

(1)

(2)

節(jié)點的速度與位移由下式得到

Vt+Δt/2=Vt-Δt/2+atΔtt

(3)

ut+Δt=Vt+Vt+Δt/2Δtt+Δt/2

(4)

其中

Δtt+Δt/2=0.5(Δtt+Δtt+Δt)

Δtt-Δt/2=0.5(Δtt-Vtt+Δt)

新的幾何構型有初始位置與位移增量疊加得到

xt+Δt=x0+ut+Δt

(5)

對于顯式方法，求解非線性問題時為了保證計算的收斂，通常需要采用較小的時間步長，其步長的設置必須滿足方程：

(6)

式中：ωmax為系統(tǒng)的最高固有振動頻率，由系統(tǒng)中最小單元的特征值方程得到

|Ke-ω2Me|=0

(7)

為確保計算能夠收斂，LS-DYNA采用了變步長積分方法。該方法會選擇當前網格中尺度最小的單元來確定步長。

2 計算模型

2.1冰的材料模型

在冰槳切削的模擬過程中，冰材料模型的準確選取直接影響結果的準確性。與水相比，海冰的固有性質非常復雜，故往往在數值模擬中采用一種典型的各向同性彈性斷裂材料，對應LS-DYNA材料庫里面的MAT13，本文取冰的相關參數如下：密度900 kg/m3、體積模量5.26 GPa、剪切模量2.20 GPa、塑性失效應變0.35%、屈服應力2.12 MPa、階段壓力-4.00 MPa、塑性硬化模量4.26 GPa。

為了驗證本文冰的數值模型和碰撞數值模擬計算方法的可行性，本文參照Kim等在2000年做的球型冰沖擊試驗[14]，進行對比驗證。其氮氣炮實驗裝置如圖1所示，Kim等實驗的球型冰直徑為42.7 mm，速度為73.5 m/s。

圖1 氮氣炮實驗裝置Fig.1 Nitrogen gas cannon experimental setup

在建模過程中，剛體半徑為0.1 m，厚度0.02 m。剛體采用拉格朗日網格，并對其進行全約束，密度為4 400 kg/m3，彈性模量為109 GPa，泊松比0.34。冰直徑為42.7 mm，網格數為18.9萬。最終計算模型如圖2所示。

圖2 沖擊計算模型Fig.2 Impact calculation model

圖3為球型冰沖擊鋼板接觸力時程曲線，通過對比實驗結果發(fā)現，力的時間分布和曲線的形狀基本相同，但力峰值相差較大。在 0～0.1 ms 的瞬間，本文數值模擬結果與Kim的實驗結果吻合較好。同時也驗證了冰的數值模型和碰撞計算方法的可行性，為接下來的冰與槳銑削碰撞研究提供了基礎。

圖3 接觸力時程曲線Fig.3 Contact force curves

2.2模型建立

本文計算的螺旋槳模型采用的是以加拿大海岸警衛(wèi)R級破冰船上裝載的四葉1200系列R-class冰級槳為原型。選冰級為PC7，即破冰船破冰厚度為1.5 m，其實槳直徑為4.1 m，本文分析對象為模型槳，直徑為200 mm。槳葉的有限元模型及網格劃分采用ANSYS前處理模塊處理，由于歐拉-拉格朗日耦合算法中結構網格與流體網格相互獨立，可以重疊，所以僅需對螺旋槳和流場分別進行建模和網格劃分。由于槳轂處與槳葉相連的部分區(qū)域結構表面不規(guī)則，因此采用了四面體網格進行劃分，而對于其他區(qū)域的槳轂則按照圓柱體進行了結構化網格劃分，劃分完成的螺旋槳網格總數為40 683。另外，為了避免槳葉與海冰接觸時的沙漏現象，將槳葉與海冰的單元設置為全積分單元。海冰和流域的均為規(guī)則幾何體，故直接在LS-PrePost劃分網格。為了減少海冰的邊界效應對結果的影響，又能充分利用計算資源的前提下，本文所選用的海冰尺寸取為200 mm×88 mm×50 mm，并且海冰與螺旋槳不發(fā)生接觸的四個面上施加徑向和周向的位移約束，劃分網格均用六面體網格，總數共計64 000。劃分完的螺旋槳和海冰網格如圖4所示。

圖4 螺旋槳和海冰網格Fig.4 Grid of propeller and ice

流域的建模采用了圓柱形的結構，由于流固耦合分析對計算資源占用較高，耗時很長，為了降低計算成本，本文中設置流域直徑為2倍螺旋槳直徑，流場進出口距螺旋槳距離為1倍螺旋槳直徑，劃分完成后的流場網格如圖5所示，網格總數為11 040。

圖5 流場網格Fig.5 Grid of flow field

建好模型后通過修改K文件，設置螺旋槳單元和海冰單元為拉格朗日網格，流場網格為歐拉網格，螺旋槳與流場之間及海冰與流場之間分別定義*CONSTRAINED_LAGRANGE_IN_SOLID關鍵字，設置螺旋槳與流場的耦合方式為罰函數約束，海冰與螺旋槳的耦合方式為考慮了結構侵蝕的罰函數約束。

3 銑削過程螺旋槳受力分析

槳切削過程中的螺旋槳受力與二者間的切削深度，即海冰距離槳軸的徑向位置有著密切的關系。選擇切削深度h=20 mm，螺旋槳進速設置為1.2 m/s，進速系數J=0.6。為便于分析，本文在結果處理時選用單個槳葉與海冰銑削的過程進行分析，主槳葉初始位置、冰槳相對位置如圖6所示，圖中圓形虛線為螺旋槳的運動軌跡，槳葉形虛線為螺旋槳開始切削與結束切削的槳葉位置，h為螺旋槳切削深度，α為螺旋槳切削度數。

圖6 螺旋槳與海冰相對位置Fig.6 Relative position of Propeller and ice

3.1槳葉應力分析

計算得到不同階段槳葉壓力面、吸力面應力云圖如圖7、8所示。數值計算結果表明，銑削工況下槳葉最大應力區(qū)域多處在與海冰發(fā)生銑削的陰影區(qū)域內。通過圖9中槳葉壓力面應力分布隨時間的變化趨勢可以看到，冰槳開始接觸后葉梢靠近導邊一側，即在接觸點附近出現應力集中。隨著時間的推移，槳葉不斷旋轉，槳葉與海冰中的接觸區(qū)域也有局部變形產的接觸區(qū)域由導邊一側逐漸移向隨邊一側，槳葉最大應力分布區(qū)域亦隨之發(fā)生移動并保持在螺旋槳與海冰發(fā)生接觸的陰影區(qū)域附近。最后，隨著槳葉的進一步旋轉，螺旋槳開始脫離海冰，二者間的相互作用結束，并且下一槳葉開始與海冰接觸。

圖7 壓力面應力云圖Fig.7 Stress counter in pressure side

圖8 吸力面應力云圖Fig.8 Stress counter in suction side

單片槳葉與海冰接觸過程中的槳葉最大應力曲線和冰槳接觸力曲線如圖9、10所示。由圖9的槳葉最大等效應力曲線可見，在冰槳切削過程中槳葉的最大等效應力也可以分為三個階段：

主槳葉切削過程中海冰的破壞形式如圖9。

第一階段，海冰與螺旋槳開始接觸(0.012～0.016 8 s)，此時主要是螺旋槳的導邊切削海冰，冰槳之間的接觸面積較小，最大等效應力的值在接觸開始后迅速增大至1.62 GPa(t=0.013 4 s)左右，之后又迅速減小，這一階段海冰的運動狀態(tài)主要表現為完整冰開始出現裂紋并逐漸開裂。從圖8也可以看出，冰槳接觸力在這一階段也是迅速上升,在t=0.016 8 s時，接觸力達到27.5 kN。

第二階段，螺旋槳進一步切入海冰(0.016 8～0.031 2 s)，此時主要是螺旋槳葉梢中部切削海冰，海冰已經部分碎裂，故強度遠不及第一階段，并且開始有部分碎冰與螺旋槳接觸，接觸面積也比第一階段大，導致槳葉的最大等效應力迅速下降，隨后再穩(wěn)步上升，這一階段應力峰值為847 MPa。冰槳接觸力在這一階段也有所下降，接觸力最大為16 kN。

第三階段也是螺旋槳導邊切削海冰(0.031 2～0.036 0 s)，并伴有大量的碎冰與螺旋槳接觸，但此時的海冰是有裂紋的冰，故產生的應力遠沒有第二階段那么大，應力峰值約為第二階段的70%，大小為507 MPa， 冰槳間的接觸力主要是碎冰與螺旋槳的相互作用，并在這一階段持續(xù)上升，最后達到16.6 kN。槳葉上方海冰已基本碎裂，切削過程結束后，應力又迅速下降，圖9中最后最大等效應力又開始迅速上升主要是因為后一片槳葉的導邊開始切削海冰了，導致應力的驟升。

圖9 槳葉最大等效應力曲線Fig.9 Blade maximum equivalent stress curve

圖10 冰槳接觸力曲線Fig.10 Contact force curve of ice and propeller

圖11 海冰的破壞過程Fig.11 Ice failure process

圖12 六個單元和兩個節(jié)點在槳葉上的位置Fig.12 Six elements and two nodes locations on the blade

為了更好地分析局部單元應力情況，選取導邊附近單元5168、5844、6588及葉梢附近單元5040、5630、6350。單元位置如圖12所示。圖13給出了切削深度為20 mm下的六個單元的等效應力曲線與槳葉最大等效應力曲線的對比圖。

從圖13可以看出，在槳葉切削的第一階段，單元5844的應力分布更接近槳葉所有單元的最大等效應力曲線，說明第一階段，槳葉的最大應力區(qū)域破壞較為嚴重的區(qū)域是導邊靠下緣部分第二階段，單元5040的應力分布曲線更接近槳葉所有單元的最大等效應力曲線，說明在這個時間段，槳葉應力集中在單元5040附近的葉梢區(qū)域。隨著切削時間的深入，槳葉應力集中區(qū)域從導邊不斷向葉梢靠攏，從圖9的不同時刻的云圖也可以看出來，應力集中區(qū)域隨著時間變化而變化。

圖13 各單元應力曲線Fig.13 Each element stress curves

3.2槳葉變形分析

槳在與海冰的切削過程中，槳葉產生的變形主要是軸向的，因此本文著重分析槳葉在軸向位移，圖14為槳葉在不同時刻的變形云圖。對比切削過程中海冰的破壞情況，海冰在切削過程中產生大變形，其變形量遠大于槳葉的變形，因為螺旋槳的彈性模量和密度遠大于海冰。圖15中，槳葉的最大變形主要集中在葉梢附近，故選取葉梢和導邊附近處的節(jié)點6269、6822進行分析，如圖15所示。

圖14 槳葉變形云圖Fig.14 Counter of blade deformation

兩個節(jié)點的軸向位移隨時間的變化如圖15所示。兩節(jié)點位移曲線基本一致，峰值均出現在螺旋槳切削的第一階段末尾，且位于導邊附近的節(jié)點6822產生的變形更大，最大可達0.002 81 m。說明在切削過程，冰對螺旋槳葉梢及導邊有很大的破壞作用。在螺旋槳設計環(huán)節(jié)，葉梢是最薄的，剛度也是最低的，因此在冰區(qū)槳設計時需要對這一部分重點關注。

圖15 節(jié)點的軸向位移曲線Fig.15 Axial displacement curves

4 結論

1)冰槳切削過程中的螺旋槳槳葉應力主要集中于葉梢及導邊附近區(qū)域，切削過程主要分為三個階段，第一階段槳葉應力最大，這與實際符合。隨著槳葉的不斷轉動，應力集中區(qū)域由葉梢的導邊一側逐漸向隨邊一側移動，整個接觸過程中最大應力的覆蓋區(qū)域主要集中在葉梢附近。

2)螺旋槳在銑削時受到的接觸力與接觸面積有關，隨著槳葉切入海冰的深度的增加，接觸面積也不斷增加，螺旋槳受到的冰槳接觸力隨之增加，但槳葉最大應力的變化卻并非隨槳葉轉動而逐漸增大。槳葉應力在冰槳接觸后即迅速達到較高水平，隨后由于海冰碎裂，碎冰與槳相互作用，海冰強度下降，導致槳葉應力峰值有所下降，最后隨著槳葉脫離海冰的作用范圍，槳葉應力又迅速下降。

針對其他變參數工況，比如改變螺旋槳轉速、海冰速度、海冰銑削深度等，還需進一步研究，以確定相關規(guī)律。

[1] VEITCH B. Predictions of ice contact forces on a marine screw propeller during the propeller-ice cutting process[J]. Acta polytech scand mech eng ser, 1995, 118: 1-110.

[2] AKINTURK A, JONES S J, MOORES C. Testing propulsion systems for performance in ice [J]. Mari-tech, 2003: 28-30.

[3] WANG J. Prediction of propeller performance on amodel podded propulsor in ice (propeller ice interaction) [D]. City of Saint John：Memorial University, 2007: 1-37, 100-140.

[4] WANG J, AKINTURK A, JONES S J, et al. Ice loads on a model podded propeller blade in milling conditions[C]//ASME 2005 24th International Conference on Offshore Mechanics and Arctic Engineering. American Society of Mechanical Engineers, 2005: 931-936.

[5] WANG J, AKINTURK A, JONES S J, et al. Ice loads acting on a model podded propeller blade [J]. Journal of offshore mechanics and Arctic engineering, 2007, 129(3): 236-244.

[6] WANG J, AKINTURK A, BOSE N, et al. Experimental study on a model azimuthing podded propulsor in ice [J]. Journal of marine science and technology, 2008, 13(3): 244-255.

[7] SOININEN H. A propeller-ice contact model [D]. Helsinki：Helsinki University of Technology, 1998: 20-70.

[8] WANG J, AKINTURK A, BOSE N. Numerical prediction of propeller performance during propeller-ice interaction[J]. Marine technology, 2009, 46(3): 123-139.

[9] HUISMAN T J , BOS R W, BROUWER J, et al. Interaction between warm model ice and a propeller[C] //Proceedings of the ASME 2014 33rd International Conference on Ocean, Offshore and Arctic Engineering, San Francisco, 2014.

[10] 何菲菲.破冰船破冰載荷與破冰能力計算方法研究[D]. 哈爾濱：哈爾濱工程大學，2010： 28-30.

HE Feifei. A computational method for ice-breaking loads and capability of icebreaker [D]. Harbin: Harbin Engineering University, 2010: 28-30.

[11] 胡志寬, 桂洪斌, 夏鵬鵬, 等. 冰載荷下船舶螺旋槳強度的有限元分析[J]. 船舶工程, 2013, 5: 6.

HU Zhikuan, GUI Hongbin, XIA Pengpeng, et al. Finite element analysis of ship propeller strength under ice loads[J]. Ship engineering, 2013 ,5: 6.

[12] 郭春雨，謝暢，趙大剛. 冰級槳水動力性能研究綜述[J].船海工程， 2014， 43(4): 1-8.

GUO Chunyu, XIE Chang, ZHAO Dagang. An overview of hydrodynamic performance research for ice-class propeller[J]. Ship & ocean engineering， 2014, 43(4)： 1-8.

[13] 郝好山，胡仁喜，康士廷, 等. ANSYS 12.0 LS-DYNA非線性有限元分析從入門到精通[M]. 北京: 上海交通大學出版社, 2010: 9.

HAO Haoshan, HU Renxi, KANG Shiting, et al. ANSYS 12 LS-DYNA nonlinear finite element analysis from entry to mastery [M]. Beijing: Shanghai Jiao Tong University Press, 2010: 9.

[14] KIM H, KEDWARD K T. Modeling hail ice impacts and predicting impact damage initiation in composite structures[J]. AIAA journal, 2000, 38(7): 1278-1288.

本文引用格式：

常欣，王錫棟，王超，等. 螺旋槳在冰槳銑削下的強度計算分析[J]. 哈爾濱工程大學學報， 2017， 38(11): 1702-1708.

CHANG Xin, WANG Xidong, WANG Chao, et al. Calculation and analysis of propeller strength under ice milling[J]. Journal of Harbin Engineering University， 2017， 38(11): 1702-1708.

Calculationandanalysisofpropellerstrengthundericemilling

CHANG Xin, WANG Xidong, WANG Chao, SUN Shengxia

(College of Shipbuilding Engineering, Harbin Engineering University, Harbin 150001, China)

When sailing in ice region, the ice-load caused by ice milling is more harmful to propeller compared to the hydrodynamic, so ships sailing in ice region need high strength requirements. After analyzing and sorting out relevant research at home and abroad, and based on explicit algorithm, this paper first proved the correctness of the calculation method by comparing related experiment with numerical simulation, then simulated the load response of propeller under the ice-milling working condition and got the result of equivalent stress of blade, displacement of blade and the process of the ice breaking. By analysis of the result, this paper concludes that the equivalent stress of ice milling process is mainly concentrated in the tip region and leading edge, the equivalent stress value rapidly will reach the peak once the ice come into contact with propeller, in the milling process the deformation of the tip region is most serious, so designer should focus on the tip part in ice propeller design process.

ice; propeller; explicit dynamics analysis; milling; strength calculating; numerical simulation

10.11990/jheu.201609008

http://www.cnki.net/kcms/detail/23.1390.u.20170427.1407.034.html

U661.1

A

1006-7043(2017)11-1702-07

2016-09-05.

網絡出版日期：2017-04-27.

國家自然科學基金項目(51379040，51679052).

常欣(1978-), 男, 副教授,博士；

王超(1981-), 男, 副教授,博士后.

王超，E-mail:wangchao0104@hrbeu.edu.cn.