基于近場動力學(xué)數(shù)值方法的冰-吊艙推進器接觸判斷研究1)

徐 佩 王 超 郭春雨 蘇玉民 葉禮裕

(哈爾濱工程大學(xué)船舶工程學(xué)院,哈爾濱 150001)

引言

極地船舶作為北極冰區(qū)安全通行的重要交通工具,在極地事業(yè)中發(fā)揮著越來越重要的作用.而推進系統(tǒng)作為船舶動力核心部分,其設(shè)計和研究需要得到更多關(guān)注.極地船舶破冰航行過程中開辟出來的碎冰航道很快被新的海冰覆蓋,以及航道周圍的平整冰引起極地船舶轉(zhuǎn)向和調(diào)頭困難等問題,而吊艙推進器具有優(yōu)于常規(guī)的“螺旋槳+舵”推進系統(tǒng)的優(yōu)勢,極大地提高了極地船舶的機動性和操縱性,很好地解決了上述遇到的問題,有利于保證極地船舶在冰區(qū)航行條件下的航行安全,成為了極地船舶廣泛采用的推進形式[1-3].當極地船舶破冰航行時,特別是尾部破冰時,不同操縱狀態(tài)的吊艙推進器與冰層直接發(fā)生切削作用,或船體擠壓產(chǎn)生的碎冰卡在吊艙推進器和船體之間使螺旋槳與冰發(fā)生連續(xù)的切削作用,從而在螺旋槳上產(chǎn)生極端冰載荷,導(dǎo)致吊艙推進器結(jié)構(gòu)損壞、性能下降等問題.因此,冰區(qū)吊艙推進器冰載荷分析成為了吊艙推進器設(shè)計的關(guān)鍵.

然而,由于吊艙推進器幾何結(jié)構(gòu)比較獨特、復(fù)雜,且存在直航、斜航以及操舵等操縱狀態(tài),給冰-吊艙推進器切削過程中的接觸判斷帶來了很大困難.目前,針對冰與海洋結(jié)構(gòu)物的接觸過程,研究人員針對不同的研究方法提出了相應(yīng)的判斷方法.在近場動力學(xué)方法研究方面,薛彥卓等[4-5]采用將船體和冰設(shè)置為近場動力學(xué)方法中不同材質(zhì)物質(zhì)點的形式,利用不同材質(zhì)物質(zhì)點之間相互排斥產(chǎn)生的排斥力作為船體冰載荷開展了破冰船冰載荷的預(yù)報.陸錫奎[6]提出了一種有限元三角形單元與冰物質(zhì)點接觸的預(yù)報方法開展了破冰船冰載荷計算,通過預(yù)判冰物質(zhì)點是否位于以船殼曲面在船內(nèi)側(cè)形心為球心的球體內(nèi)、冰物質(zhì)點與船殼單元曲面的距離是否小于冰物質(zhì)點半徑以及利用重心法判斷冰物質(zhì)點的投影點是否位于船殼單元曲面內(nèi)的方法進行了冰物質(zhì)點與船體的接觸判斷.葉禮裕等[7-10]針對冰區(qū)螺旋槳利用近場動力學(xué)方法和面元法提出了一種連續(xù)接觸判斷算法,通過求解面元法向量、冰物質(zhì)點與面元向量以及平面方程的方式實現(xiàn)了冰槳接觸判斷,但僅適用于螺旋槳與冰的某一航行狀態(tài).在離散元方法研究方面,季順迎等[11]采用將船體離散為一系列圓盤單元和碎冰簡化為三維圓盤的方式開展了船舶在碎冰區(qū)航行時冰載荷預(yù)報,通過圓盤之間面-面、面-弧、弧-弧的接觸判斷以及求解圓盤間作用力的方式計算了船舶冰載荷.李紫麟等[12-13]、季順迎等[14-17]將船體結(jié)構(gòu)簡化為三角形殼單元、碎冰簡化為圓盤,利用面-面、弧-面以及弧-棱的接觸判斷方式開展了破冰船在碎冰區(qū)和平整冰區(qū)的冰載荷預(yù)報.王超等[18]基于STAR-CCM+軟件利用CFD-DEM 耦合方法開展了冰-螺旋槳碰撞狀態(tài)時冰載荷預(yù)報,采用以赫茲-明德林接觸理論為基礎(chǔ)的非線性彈簧阻尼接觸模型直接求解壁面法向力和切向力的方式進行冰載荷計算,但由于模擬固體斷裂方面存在不足,無法進行冰槳切削狀態(tài)的計算.在有限元方法研究方面,丁仕風(fēng)等[19]、鄒早建等[20]基于LS-DYNA 軟件開展了冰槳碰撞和切削狀態(tài)的數(shù)值模擬研究,其中,冰槳之間的接觸判斷采用罰函數(shù)法,即首先檢查冰塊節(jié)點是否穿透槳葉主表面,若未穿透,不做任何處理.若穿透,需要在槳葉主表面和冰從表面之間放置一系列法向彈簧求解冰載荷,但也僅研究了螺旋槳與冰的某一航行狀態(tài).

因此,本文針對復(fù)雜工況下冰-海洋結(jié)構(gòu)物接觸判斷方法研究的不足,借鑒先進的近場動力學(xué)理論和面元法,根據(jù)吊艙推進器獨特的結(jié)構(gòu)和操縱形式,提出吊艙推進器在直航、斜航以及操舵狀態(tài)下與冰的接觸判斷方法,為冰區(qū)吊艙推進器數(shù)值預(yù)報方法的建立及強度優(yōu)化設(shè)計提供技術(shù)支撐.

1 近場動力學(xué)理論

近場動力學(xué)方法由美國Sandia 國家實驗室的Silling 博士提出的一種無網(wǎng)格方法[21-22],他利用積分方式控制方程,而非偏微分方程求解,避免了對裂紋尖端位移場求空間導(dǎo)數(shù),有效地減小了列尖的奇異性,避免了連續(xù)介質(zhì)力學(xué)求解大變形問題的困難,非常適合于求解材料的大尺度變形問題,因此,近場動力學(xué)方法模擬冰的破碎具有非常好的優(yōu)勢.

1.1 基本方程

在近場動力學(xué)方法求解過程中,需要將研究對象離散為一系列的物質(zhì)點,關(guān)注物質(zhì)點在其作用范圍內(nèi)與所有物質(zhì)點相互作用時的物理現(xiàn)象,如圖1所示.在圖中,為了便于研究物質(zhì)點空間位置的變化,對物質(zhì)點的相關(guān)信息進行了定義,物質(zhì)點的位移表示為u,物質(zhì)點之間的相對位置和相對位移分別表示為ξ和η,即ξ=x′-x和η=μ(x′,t)-μ(x,t),ξ+η表示兩個物質(zhì)點之間的相對位置矢量.而物質(zhì)點之間的相互作用通過“鍵”進行表示,與彈簧類似,任意一個物質(zhì)點x與其領(lǐng)域半徑δ 內(nèi)的物質(zhì)點發(fā)生作用,而領(lǐng)域半徑δ 之外的物質(zhì)點與物質(zhì)點x不發(fā)生相互作用,即

圖1 物質(zhì)點x 與物質(zhì)點x′ 之間的相互作用Fig.1 Interaction between material particle x and material particle x′

在參考坐標系中,物質(zhì)點x的近場動力學(xué)運動方程可應(yīng)用虛功原理推導(dǎo)得出[23],即

其中,Hx是物質(zhì)點x及其領(lǐng)域半徑內(nèi)物質(zhì)點構(gòu)成的域,u是物質(zhì)點x的位移,ρ 是物質(zhì)點密度,f是物質(zhì)點x′與x之間的力密度,b是外部載荷.

1.2 本構(gòu)方程

為了模擬冰塊的斷裂、變形以及裂紋形成和延伸的過程,引入鍵伸長率的概念來表示物質(zhì)點之間相互作用力的大小,即[24]

其中,當鍵在拉伸狀態(tài)時,s為正值;在壓縮狀態(tài)時,s為負值.由于標量形式的力密度函數(shù)f不依賴于相對位置ξ,因此,該材料為各向同性的.

此外,為了將材料的破壞準則引入到本構(gòu)方程中,需要定義一個臨界伸長率,當鍵的伸長率大于臨界伸長率時鍵發(fā)生斷裂,斷裂完成后,物質(zhì)點之間的力密度為零.以典型的PMB (prototype micro-elastic brittle)材料為例[25],其力密度函數(shù)表示為

其中,g為線性標量函數(shù),其表達式為:g(s)=cs?s,s為鍵伸長率,c為鍵常數(shù).μ 為與時間相關(guān)的標量函數(shù),μ=0 時,鍵發(fā)生斷裂,μ=1 時,未發(fā)生斷裂.

其中,s0為極限伸長率,是判斷鍵能否發(fā)生斷裂的標準,在具體某個時刻認為是固定不變的,如圖2 所示.對于PMB 材料,極限伸長率s0為

其中,G0是單位斷裂面斷裂時所釋放的鍵能,κ 是體積模量,δ 是物質(zhì)點的領(lǐng)域半徑.

圖2 PMB 材料力和伸長率的關(guān)系[26]Fig.2 Relationship between force and elongation of PMB material[26]

雖然PMB 材料是各向同性的,但是,在某一特定條件下,例如當冰垂向鍵發(fā)生斷裂時冰的近場動力學(xué)物質(zhì)點在該方向上相互作用消失,則其在橫向方向能構(gòu)成承受更大或更小的拉/壓應(yīng)力,此時,冰表現(xiàn)出各向異性特征[27].因此,引入鍵斷裂的概念去清晰地描述單個物質(zhì)點的局部損傷,其定義為

其中,0 ≤φ ≤1,當φ=0 時,材料為初始狀態(tài),當φ=1 時,表明物質(zhì)點周圍的所有鍵發(fā)生斷裂.鍵發(fā)生斷裂后不能繼續(xù)承受拉力,將會導(dǎo)致材料的軟化現(xiàn)象,進而引起損傷的進一步深化,斷裂鍵的累積最終會形成一個裂紋的表面.

1.3 方程離散

近場動力學(xué)是一種無網(wǎng)格方法,其運動方程是一種積分-微分方程,通常不能得出解析解[28].因此,在數(shù)值計算過程中通常將物體離散成物質(zhì)點的形式,然后對每一個物質(zhì)點的運動方程進行積分,方程(1)離散后的每個物質(zhì)點的運動方程可表示為

其中,N為領(lǐng)域范圍內(nèi)子域的數(shù)量,Ne為第e個子域內(nèi)積分點的數(shù)量,x(k)和x(j)分別為第k個和第j個積分點的位置,V(j)為第j個立方體子域的體積.

2 吊艙推進器表面離散及模型簡化

2.1 吊艙推進器模型

本文模擬使用的冰區(qū)吊艙推進器模型由中國船舶重工集團公司第七〇二研究所和課題組聯(lián)合設(shè)計,其中,螺旋槳直徑為0.25 m,葉數(shù)為4,縱傾角為10°,轂徑比為0.368,螺距比(0.7R) 為0.775,縮尺比為1:16.48,艙體長度為0.648 m,艙體最大直徑為0.15 m,艙體前端傾斜角度為12.3°,艙體后端長度為0.22 m,支架高度為0.175 m,支架弦長為0.24 m,支架最大寬度為0.065 m.

2.2 吊艙推進器表面離散

在吊艙推進器表面離散過程中,由于其結(jié)構(gòu)、組成比較復(fù)雜,且螺旋槳的表面曲率較大,給冰-吊艙推進器的接觸判斷問題帶來很大困難.為了解決冰-吊艙推進器接觸判斷的問題,本文借助面元法面元劃分的思想,將吊艙推進器表面進行網(wǎng)格劃分來逼近吊艙推進器表面.根據(jù)文獻[7]中對螺旋槳槳葉網(wǎng)格劃分的討論,螺旋槳槳葉徑向和弦向均采用余弦網(wǎng)格劃分,其中,余弦網(wǎng)格劃分的節(jié)點表示為

其中,R0為螺旋槳半徑;Rh為槳轂半徑;rj和sj分別為徑向和弦向網(wǎng)格劃分點的位置;bj為rj處葉剖面的弦長;Nr和Nc分別為徑向和弦向網(wǎng)格數(shù).

圖3 展示出了吊艙推進器表面網(wǎng)格劃分,其中,吊艙推進器表面離散為若干個四邊形單元,在網(wǎng)格劃分過程中,將吊艙推進器劃分為3 段,第一段為艙體前端到支架前端,第二段為與支架相接的部分,第三段為支架后端到艙體尾部.第一段、第三段軸向和周向均采用平均網(wǎng)格劃分,艙體第二段軸向劃分位置與支架相同,支架軸向采用余弦網(wǎng)格劃分,高度上采用平均網(wǎng)格劃分,周向采用平均網(wǎng)格劃分.

圖3 吊艙推進器網(wǎng)格劃分Fig.3 Grid division of podded propulsor

2.3 冰?吊艙推進器模型簡化

根據(jù)冰-吊艙推進器的實際接觸場景,數(shù)值模擬冰-吊艙推進器的切削過程存在兩個問題:(1) 根據(jù)極地海冰實際場景,冰塊通常是不規(guī)則的,這給數(shù)值建模帶來很多問題,且很難將冰塊離散為物質(zhì)點的形式;(2)若不考慮冰-冰、冰-吊艙推進器-船體之間的擠壓作用,只建立冰-吊艙推進器接觸模型,一旦冰塊和螺旋槳接觸,冰塊會由于不受約束而遠離螺旋槳,難以形成連續(xù)的切削過程.

此外,由于數(shù)值計算方法和計算能力的限制,目前建立冰-流-吊艙推進器耦合計算模型模擬冰區(qū)吊艙推進器的實際切削過程有很大難度.因而,為了便于開展冰-吊艙推進器切削狀態(tài)的研究,本文建立的冰-吊艙推進器切削數(shù)值模型忽略了流體的影響.雖然這種數(shù)值模型模擬結(jié)果與實際情況存在差異,無法體現(xiàn)出冰的存在對吊艙推進器水動力性能的影響,但是能夠比較真實地再現(xiàn)切削過程中冰塊的破碎和瞬態(tài)冰載荷的變化.根據(jù)Wang 等[29-30]開展的系列實驗測量數(shù)據(jù)可以看出,冰槳切削過程中螺旋槳受到的冰載荷要比水動力載荷高一個量級以上,水動力載荷對于螺旋槳總載荷的作用幾乎可忽略不計.事實上,有些學(xué)者為了便于分析冰槳切削狀態(tài)時冰的破碎模式和冰載荷大小,在空氣中開展了冰槳切削實驗,Bach 等[31]認為研究冰槳接觸冰載荷可在空氣中開展冰槳銑削實驗,將水動力載荷從總載荷中區(qū)分出來,有利于冰槳銑削過程的基礎(chǔ)研究.

為此,分別針對吊艙推進器的不同操縱狀態(tài)對冰-吊艙推進器模型及運動模式進行簡化,在此之前,為了便于說明螺旋槳旋轉(zhuǎn)、吊艙推進器操縱和冰塊運動,對吊艙推進器坐標系進行定義,如圖4 所示.螺旋槳局部坐標系為OP-XP,YP,ZP,OP位于槳盤面中心,XP正方向是從螺旋槳指向艙體的方向,YP正方向垂直向上,ZP正向根據(jù)右手定則確定.為了分析吊艙推進器與船體連接處扭矩的變化,建立了吊艙推進器局部坐標系OU-XU,YU,ZU,OU位于支架最上端翼型弦長的中點處,XU,YU和ZU方向與XP,YP和ZP的方向相同.吊艙推進器整體坐標系為O-X,Y,Z,O位于艙體軸線與支架弦長中點垂線的交點處,X,Y和Z方向與XP,YP和ZP的方向相同.

圖4 吊艙推進器的坐標系Fig.4 Coordinate system of podded propulsor

當?shù)跖撏七M器在直航狀態(tài)與冰切削時,假設(shè)吊艙推進器處于靜止狀態(tài),螺旋槳繞自身坐標系旋轉(zhuǎn).同時,為了獲得相對穩(wěn)定和周期性的冰載荷,假設(shè)冰塊由大量的物質(zhì)點構(gòu)成,且以恒定的速度逐漸向吊艙推進器靠近,冰塊的形狀可根據(jù)需要將其簡化為規(guī)則形狀,例如長方體等.在數(shù)值模擬過程中,為了使冰塊的運動狀態(tài)更加接近實際情況,在冰塊后端進行位移約束和速度設(shè)置,其中,對冰塊后端幾層的物質(zhì)點進行位移約束,同時,進行速度設(shè)置,從而保證冰塊后端的物質(zhì)點以恒定的速度推動整塊冰運動,其中,模型簡化如圖5 所示.

圖5 冰-吊艙推進器切削工況模型簡化(直航狀態(tài))Fig.5 Simplification of model for milling conditions of ice-podded propulsor(straight ahead)

根據(jù)吊艙推進器在斜航狀態(tài)與冰的切削場景[33-34],對冰-吊艙推進器計算模型進行簡化,如圖6 所示.假設(shè)吊艙推進器繞整體坐標系中的Y軸進行偏轉(zhuǎn),螺旋槳變?yōu)榱嗽谛焙綘顟B(tài)局部坐標系

此外,考慮到吊艙推進器在操舵狀態(tài)與冰層或者冰脊的接觸場景,假設(shè)吊艙推進器在整體坐標系中位置不變,繞整體坐標系中的Y軸進行操舵,而螺旋槳除繞整體坐標系中的Y軸轉(zhuǎn)動外也圍繞自身的局部坐標系進行旋轉(zhuǎn).其中,對冰塊的形狀和運動狀態(tài)的假設(shè)與吊艙推進器直航狀態(tài)時相同,吊艙推進器操舵狀態(tài)與冰的接觸場景和模型簡化如圖7 所示.

圖6 冰-吊艙推進器切削工況模型簡化Fig.6 Simplification of model for milling conditions of ice-podded propulsor(azimuthing condition)

圖7 冰-吊艙推進器切削工況模型簡化(操舵狀態(tài))Fig.7 Simplification of model for milling conditions of ice-podded propulsor(dynamic azimuthing condition)

3 冰?吊艙推進器不同作用方式下接觸判斷方法

3.1 吊艙推進器直航狀態(tài)時接觸判斷方法

為了減少冰-螺旋槳接觸檢測的計算量和復(fù)雜的運算過程,同時,保證接觸判斷的精度和效率,根據(jù)螺旋槳與冰的接觸算法進行吊艙推進器直航狀態(tài)下螺旋槳與冰物質(zhì)點的判斷[7].首先,將冰物質(zhì)點轉(zhuǎn)換為極坐標的形式,檢測物質(zhì)點是否位于槳葉葉根和葉梢所包圍的范圍內(nèi),如果在此范圍以外說明冰槳未接觸,否則進行下一步判斷,如圖8(a)所示.如果在葉根和葉梢所包圍的范圍內(nèi),再判斷物質(zhì)點是否位于槳葉導(dǎo)邊和隨邊之間的角度方位內(nèi),若不在該方位角內(nèi),說明沒有接觸,反之將其定義為可能與螺旋槳接觸的物質(zhì)點,如圖8(b)所示.

圖8 冰與槳葉接觸的區(qū)域判別法(直航狀態(tài))Fig.8 Regional discrimination method of ice contacting with blades(straight ahead)

3.2 吊艙推進器斜航狀態(tài)時接觸判斷方法

當開展吊艙推進器斜航狀態(tài)冰載荷預(yù)報時,首先實現(xiàn)吊艙推進器斜航狀態(tài)的調(diào)整,根據(jù)螺旋槳和支架、艙體的不同運動特征,分別設(shè)置了不同的偏轉(zhuǎn)方程,螺旋槳偏轉(zhuǎn)方程為

其中,X,Y和Z分別為直航狀態(tài)時槳葉面元3 個方向的坐標,th1 為偏轉(zhuǎn)角度,β 為吊艙推進器的斜航角度,X′,Y′,Z′分別為斜航狀態(tài)時槳葉面元3 個方向的坐標.

支架和艙體偏轉(zhuǎn)方程為

其中,X,Y和Z分別為直航狀態(tài)時面元3 個方向的坐標,X′′,Y′′和Z′′為中間變量,X′,Y′和Z′分別為斜航狀態(tài)時面元3 個方向的坐標.

當?shù)跖撏七M器處于斜航狀態(tài)時螺旋槳也呈現(xiàn)出了斜航狀態(tài),此時,螺旋槳的局部坐標與冰物質(zhì)點進行極坐標轉(zhuǎn)化的坐標系OP-XP,YP,ZP出現(xiàn)了一定的夾角,如圖9 所示.

在圖9 中,冰物質(zhì)點無法直接使用螺旋槳的局部坐標系進行極坐標轉(zhuǎn)化,故吊艙推進器在斜航狀態(tài)時冰物質(zhì)點與螺旋槳的接觸判別方法與直航狀態(tài)時出現(xiàn)差異.因此,為了繼續(xù)使用冰物質(zhì)點與螺旋槳接觸的徑向和周向判別法,提出了一種吊艙推進器斜航狀態(tài)時冰槳接觸判斷方法,即在冰物質(zhì)點與槳葉接觸判別之前加入冰物質(zhì)點的三維偏轉(zhuǎn)方程,如式(18) 所示,該方程可將冰物質(zhì)點轉(zhuǎn)化到螺旋槳的局部坐標系中,使冰物質(zhì)點在螺旋槳斜航狀態(tài)時的局部坐標系中進行極坐標變換.其中,冰物質(zhì)點的偏轉(zhuǎn)方程為

圖9 冰與槳葉接觸的區(qū)域判別(斜航狀態(tài))Fig.9 Regional discrimination method of ice contacting with blades(azimuthing condition)

其中,為冰物質(zhì)點在螺旋槳斜航狀態(tài)時局部坐標系坐標,XP,YP,ZP為直航狀態(tài)時螺旋槳局部坐標系坐標.

3.3 吊艙推進器操舵狀態(tài)時接觸判斷方法

開展吊艙推進器操舵狀態(tài)冰載荷預(yù)報時,首先實現(xiàn)吊艙推進器中螺旋槳、艙體和支架的操舵運動,根據(jù)冰-吊艙推進器接觸載荷預(yù)報程序的特點和優(yōu)勢,選擇以迭代時間步長為時間變量,設(shè)置了以迭代時間步長為變量的操舵方程.在吊艙推進器操舵過程中,將操舵過程分為3 個階段(以俯視圖為例),第一階段為直航狀態(tài)逆時針操舵到左側(cè)最大斜航角度,第二階段為左側(cè)最大斜航角度順時針操舵到右側(cè)最大斜航角度,第三階段為右側(cè)最大斜航角度操舵到直航狀態(tài),后續(xù)操舵過程依次進行.在吊艙推進器操舵方程設(shè)置的過程中,由于螺旋槳既要實現(xiàn)自身的旋轉(zhuǎn)又要進行操舵,而艙體和支架只需要進行操舵,為此,分別針對螺旋槳、艙體和支架設(shè)置單獨的操舵方程,其中,吊艙推進器操舵角度變化方程如式(19)～式(20)所示.求解完成后將操舵角度方程代入吊艙推進器偏轉(zhuǎn)方程(13)～式(17)中即可實現(xiàn)吊艙推進器的操舵運動.

螺旋槳操舵角度變化方程為

其中,th1 為某一時刻吊艙推進器的偏轉(zhuǎn)角度,β 為操舵過程中最大斜航角度,tt1 為迭代步數(shù),ttt為螺旋槳偏轉(zhuǎn)次數(shù),ω 為操舵速率,dt為時間間隔.

支架和艙體操舵角度變化方程為

其中,tttt為艙體和支架的偏轉(zhuǎn)次數(shù),其他定義同式(9).

其中,X,Y和Z為冰物質(zhì)點在直航狀態(tài)時坐標系坐標,X′,Y′和Z′為冰物質(zhì)點在操舵過程中不同斜航角度時的坐標,其他參數(shù)的定義如上文所示.

圖10 冰與槳葉接觸的區(qū)域判別(操舵狀態(tài))Fig.10 Regional discrimination method of ice contacting with blades(dynamic azimuthing condition)

在確定了可能與螺旋槳發(fā)生接觸的物質(zhì)點后,需要通過幾何關(guān)系來檢測冰物質(zhì)點與槳葉是否發(fā)生接觸.由于螺旋槳槳葉具有葉面和葉背兩個面,因此,根據(jù)葉面和葉背的面元以及冰物質(zhì)點之間的關(guān)系,求解兩個面元的中點、法向量以及物質(zhì)點與面元的向量,求解得到的葉面和葉背的中心點坐標分別為(Pcf(1),Pcf(2),Pcf(3)) 和(Pcb(1),Pcb(2),Pcb(3)),葉面和葉背的法向量分別為nf和nb,物質(zhì)點與面元的向量分別為a和b,其中,向量a=對于每一個可能與槳葉接觸的物質(zhì)點,當t時刻物質(zhì)點位于葉面之外,而t+Δt時刻物質(zhì)點坐標與面元法向量和面元中心點坐標滿足式(24)時,說明物質(zhì)點和槳葉發(fā)生接觸,相關(guān)判斷過程如圖11 所示.

圖11 物質(zhì)點與槳葉面元接觸的判斷過程Fig.11 Judgment process of material point contacting with blade surface

3.4 冰?吊艙推進器接觸載荷計算

在冰-吊艙推進器切削過程中,槳葉與冰物質(zhì)點將發(fā)生接觸,接觸完成后,物質(zhì)點將會進入槳葉內(nèi)部,但為了真實反映冰槳切削的物理過程,將進入到槳葉內(nèi)部的物質(zhì)點進行位置的重新分配,即將該物質(zhì)點分配到槳葉附近,其中,滲入到槳葉內(nèi)的物質(zhì)點和新分配物質(zhì)點的距離為

其中,nb和nf分別為螺旋槳葉背和葉面中面元的法向量;a和b分別為冰物質(zhì)點的坐標與葉背和葉面面元中心的向量.

新分配的冰物質(zhì)點坐標計算公式如下

其中,k為第k個冰物質(zhì)點,V0為冰物質(zhì)點的速度,n為接觸點面元的法向量.

此外,在冰與槳葉接觸后的冰物質(zhì)點x(k)在t+Δt時刻的速度可表示為

通過冰-吊艙推進器接觸判斷和物質(zhì)點接觸力的計算,可判斷出計算域內(nèi)任意一個冰物質(zhì)點是否與螺旋槳接觸,若發(fā)生接觸,則判斷出物質(zhì)點與螺旋槳中那個槳葉的第幾個面元發(fā)生接觸以及對應(yīng)的接觸力大小.在t+Δt時刻螺旋槳表面任意一個面元的接觸力可通過疊加與面元接觸的所有物質(zhì)點的接觸力,如式(29)所示

從而可以計算出t+Δt時刻螺旋槳表面上任意一個面元j的接觸壓力,如式(31)所示

在吊艙單元冰載荷預(yù)報過程中,由于本文模擬未涉及流體的影響,所以,假設(shè)吊艙單元整體軸向力和螺旋槳受到的軸向力相等,如式(34)所示,而吊艙單元扭矩的計算公式,如式(35)所示

其中,LX和LY分別為支架最上端翼型弦長中點與槳盤面中心之間的距離在X軸和Y軸的分量,LX和LY對應(yīng)于吊艙推進器直航狀態(tài).

由于吊艙推進器在斜航狀態(tài)和操舵過程中槳盤面到支架弦長中點的力臂在X軸和Z軸的分量不斷變化,故提出了吊艙推進器操舵過程中吊艙單元扭矩的求解公式,如式(36)所示

其中,th1 為某一時刻吊艙推進器的偏轉(zhuǎn)角度,當th1 ≥0 時,利用式(36a) 進行吊艙單元扭矩的求解,反之利用式(36b)進行求解.當偏轉(zhuǎn)角度th1 達到最大值時為吊艙推進器斜航狀態(tài)時扭矩的計算公式.

4 不同接觸判斷方法在冰?吊艙推進器切削狀態(tài)中的應(yīng)用及分析

4.1 直航狀態(tài)下冰?吊艙推進器接觸判斷方法的應(yīng)用及分析

在數(shù)值模擬過程中,以螺旋槳模型網(wǎng)格劃分數(shù)目22×22,冰的密度ρ=800 kg/m3、泊松比μ=0.25、彈性模量E=280 MPa、物質(zhì)點大小dx=0.6D/45、領(lǐng)域半徑δ=3.015dx[35]、極限伸長率scr=0.06、時間步長dt=3.6×10-6為參數(shù)進行設(shè)置,并與Wang等[30]在冰水池內(nèi)開展的冰-吊艙推進器切削試驗進行對比和分析,對比結(jié)果如圖12 所示.

圖12 給出了數(shù)值模擬過程中單槳葉X軸方向上軸向力和扭矩隨時間的變化曲線與試驗值的對比,由于冰塊長度設(shè)置有限,這里僅展示出了600 ms 時間范圍內(nèi)槳葉冰載荷的變化曲線.由圖12 可知,兩條曲線波峰附近區(qū)域?qū)儆诒鶚邢鲄^(qū)域,從整體來看數(shù)值計算值和試驗值的冰載荷數(shù)量級、冰載荷變化趨勢以及冰載荷峰值基本一致,兩條曲線的吻合性較好,但是由于冰-螺旋槳切削作用的隨機性,不可能使兩條曲線完全一致.此外,在數(shù)值計算時,由于未考慮流體的影響,當槳葉與冰未發(fā)生切削作用時槳葉受到的載荷為零.而在冰水池試驗時,由于流體的存在,當槳葉與冰未發(fā)生切削作用時槳葉受到了水動力載荷的作用,且可以觀察出力也發(fā)生了波動,這可能是由于在冰槳相互作用下槳葉附近的水流比較復(fù)雜,伴流場不均勻引起的,故槳葉與冰未發(fā)生切削作用時數(shù)值計算值和試驗值載荷的變化存在一定差異,因此,可以證明在冰載荷數(shù)值預(yù)報過程中參數(shù)的設(shè)置具有合理性和準確性.

圖12 軸向力和扭矩的數(shù)值計算值與試驗值對比Fig.12 Comparison between numerical calculation value and test value of axial force and torque

在直航狀態(tài)模擬過程中,以冰塊推送速度V=0.2 m/s,螺旋槳轉(zhuǎn)速n=3 rad/s,冰槳切削深度h=35 mm 為例進行分析,其中,冰-吊艙推進器切削過程以及冰載荷變化曲線如圖13 和圖14 所示.

圖13 給出了吊艙推進器直航狀態(tài)與冰切削過程中不同時刻時冰-吊艙推進器的切削狀態(tài),圖中冰塊的顏色代表其內(nèi)部鍵的破壞狀態(tài),藍色為鍵未發(fā)生斷裂,紅色為鍵發(fā)生斷裂.圖13(a)為冰-吊艙推進器開始切削的初始階段,此時冰槳未接觸,冰塊未發(fā)生破碎,冰物質(zhì)點全部為藍色,與冰的實際運動狀態(tài)相符.隨著冰塊的運動和螺旋槳的旋轉(zhuǎn),當t=104 ms時,螺旋槳槳葉與冰塊開始接觸,冰塊前端的物質(zhì)點顏色變?yōu)榧t色,表明此時冰物質(zhì)點發(fā)生了破碎現(xiàn)象.隨著冰塊向吊艙推進器逐漸靠近,槳葉在冰塊內(nèi)部的切削深度逐漸增加,冰物質(zhì)點破碎的區(qū)域也逐漸增加,切削后的碎冰體積也相對較大,如圖13(c) 所示.在圖13(d) 中,槳葉在冰塊內(nèi)部形成的凹槽深度和寬度達到最大,槳葉切削后的碎冰形狀也類似,冰-螺旋槳的切削過程達到穩(wěn)定,因此,隨著冰塊的繼續(xù)運動,冰槳開始進入穩(wěn)定性和周期性的切削過程.

圖13 冰-吊艙推進器的切削過程(直航狀態(tài))Fig.13 Milling process of ice-podded propulsor(straight ahead)

圖13 冰-吊艙推進器的切削過程(直航狀態(tài))(續(xù))Fig.13 Milling process of ice-podded propulsor(straight ahead)(continued)

由圖14 可知,在t=100 ms 左右時吊艙推進器冰載荷開始增大,而冰載荷在t=100～170 s 范圍內(nèi)的峰值明顯小于其他時間范圍,主要原因是該時間段內(nèi)冰槳初步接觸,槳葉在冰塊內(nèi)的切削深度較小,該現(xiàn)象與圖13 中冰-吊艙推進器的切削過程相一致.隨著冰與螺旋槳的切削,吊艙推進器冰載荷在t=170～250 s 范圍的峰值與在t=100～170 s 的冰載荷的峰值相比有所增大,但仍然小于其他時間段.當?shù)跖撏七M器冰載荷的第3 個峰值出現(xiàn)后,冰載荷開始出現(xiàn)周期性分布,表明槳葉與冰進入了穩(wěn)定的切削過程,與圖13 中冰槳接觸時冰物質(zhì)點的破碎現(xiàn)象相對應(yīng),從而證明了吊艙推進器直航狀態(tài)下與冰切削過程中冰槳接觸判斷方法的正確性.將螺旋槳不同方向的軸向力或扭矩進行對比可知,該工況下螺旋槳X軸方向上的力明顯大于Y軸和Z軸,且X軸方向上的力為負值,表明槳葉葉面與冰物質(zhì)點的接觸面積大于葉背與冰物質(zhì)點的接觸面積,從而產(chǎn)生了負方向的軸向力.Y軸方向上的力為負正交替出現(xiàn)的現(xiàn)象,該現(xiàn)象與槳葉的切削過程以及坐標系相關(guān),即槳葉進入冰塊時槳葉由下向上切削冰,槳葉受到的冰載荷向下為負,而槳葉旋出冰塊時槳葉由上向下切削冰,槳葉受到的冰載荷向上為正.Z軸方向上的力始終為正,即與螺旋槳的旋轉(zhuǎn)方向和Z軸的方向相關(guān).但X軸方向上的扭矩明顯小于Y軸和Z軸.吊艙單元整體扭矩隨著時間的變化也呈現(xiàn)出了由小到大再到周期性變化的現(xiàn)象,與實際情況相符.通過對比3 個方向上吊艙單元扭矩的變化可知,Z軸方向上的扭矩大于Y軸也大于X軸,主要原因與X軸方向上力較大以及支架的長度相關(guān).此外,由圖14 可知,在冰-吊艙推進器切削過程中吊艙推進器冰載荷出現(xiàn)了劇烈的波動,將可能導(dǎo)致艙體內(nèi)的軸系、推進電機、軸承、密封裝置等部件和子系統(tǒng)失效,對冰區(qū)航行船舶的安全航行產(chǎn)生非常大的影響.

4.2 斜航狀態(tài)下冰?吊艙推進器接觸判斷方法的應(yīng)用及分析

在斜航狀態(tài)模擬過程中,吊艙推進器沿整體坐標系中的Y軸進行轉(zhuǎn)動,以吊艙推進器斜航角度β=30°,60°,-30°和-60°(順時針為負、逆時針為正)、螺旋槳轉(zhuǎn)速n=3 r/s,冰塊運動速度V=0.2 m/s、冰槳切削深度h=35 mm 為例進行分析,其中,冰-吊艙推進器切削過程以及冰載荷變化曲線如圖15 和圖16 所示.

圖15 不同斜航角度時冰-吊艙推進器的切削過程Fig.15 Milling process of the ice-podded propulsor at different angles of azimuthing

由圖15 可知,當?shù)跖撏七M器斜航狀態(tài)與冰切削時,與槳葉接觸的冰物質(zhì)點顏色發(fā)生了變化,未接觸的冰物質(zhì)點顏色仍然為藍色,表明接觸區(qū)域內(nèi)冰塊發(fā)生了破碎,即槳葉與冰物質(zhì)點能夠正常判斷.當?shù)跖撏七M器斜航角度β=30°偏轉(zhuǎn)到β=60°時,槳葉拋出的碎冰體積越小,主要原因是隨著吊艙推進器斜航角度的增加,槳葉與冰塊的接觸面積越小,故冰塊內(nèi)部的破碎區(qū)域越小.同時,在圖15(a)中可以觀察出,槳葉拋出后的碎冰更多地碰撞到艙體上,將引起吊艙推進器的冰激振動或者瞬時高應(yīng)力.因而,在冰區(qū)吊艙推進器設(shè)計過程中對艙體前端的厚度進行適當調(diào)整.當斜航角度β=-30°偏轉(zhuǎn)到β=-60°時,槳葉拋出后的碎冰體積也越小,與β=30°偏轉(zhuǎn)到β=60°時的現(xiàn)象相同.以斜航角度β=30°和β=-30°為例,對比了吊艙推進器逆時針和順時針偏轉(zhuǎn)相同角度時冰槳的切削過程,由圖15(a)和圖15(c)可知,吊艙推進器順時針轉(zhuǎn)動時槳葉切削后的碎冰體積更大,即對冰塊的破壞程度越大,主要原因是β=30°時,槳葉的旋轉(zhuǎn)方向與冰塊的運動方向相反,槳葉導(dǎo)邊逐漸切削冰塊.而當β=-30°時,槳葉的旋轉(zhuǎn)方向與冰塊的運動方向相同,且槳葉存在后傾等原因,使槳葉對冰塊的破壞程度更大.因此,在吊艙推進器斜航狀態(tài)破冰過程中,需要綜合考慮吊艙推進器的斜航角度、螺旋槳的旋轉(zhuǎn)方向以及船體的運動方向等因素.

此外,為了進一步對比吊艙推進器不同斜航角度時冰載荷的變化以及證明接觸判斷方法的正確性,展示出了螺旋槳中某一槳葉在X軸方向的軸向力和Y軸方向的扭矩以及吊艙單元整體在Y軸和Z軸方向的扭矩,如圖16 所示.

圖16 吊艙推進器冰載荷隨斜航角度的變化Fig.16 Change of ice load of the podded propulsor with the angle of azimuthing

由圖16(a) 和圖16(b) 可知,吊艙推進器的斜航角度越小,螺旋槳軸向力和扭矩越大,槳葉與冰接觸的時間越長.當?shù)跖撏七M器逆時針和順時針偏轉(zhuǎn)相同角度時,順時針偏轉(zhuǎn)時螺旋槳受到的冰載荷較大,但槳葉與冰的接觸時間較短,與圖15 中吊艙推進器順時針偏轉(zhuǎn)時槳葉切削后的碎冰體積越大相一致,證明了本文提出的吊艙推進器斜航狀態(tài)與冰接觸判斷方法的正確性.吊艙單元整體扭矩的變化如圖16(c)和圖16(d) 所示,由圖16(c) 可知,吊艙推進器逆時針轉(zhuǎn)動時,斜航角度越大,吊艙單元沿Y軸的操舵扭矩越小,而順時針轉(zhuǎn)動時,斜航角度越大,吊艙單元沿Y軸的操舵扭矩越大,且逆時針和順時針偏轉(zhuǎn)相同角度時,順時針偏轉(zhuǎn)時吊艙單元沿Y軸的操舵扭矩較大,主要原因與順時針偏轉(zhuǎn)時螺旋槳受到X軸方向的軸向力越大相關(guān).吊艙單元整體扭矩沿Z軸方向的變化如圖16(d)所示,圖中展示出吊艙推進器斜航角度越大,吊艙單元沿Z軸的扭矩越小.吊艙推進器逆時針和順時針偏轉(zhuǎn)相同角度時,吊艙單元沿Z軸的扭矩差異較小.因此,可以得出吊艙推進器順時針偏轉(zhuǎn)相同角度與冰層切削時對冰的破壞效果較好,但螺旋槳冰載荷和吊艙單元操舵扭矩較大,對吊艙推進器結(jié)構(gòu)強度的要求也越高.

4.3 操舵狀態(tài)下冰?吊艙推進器接觸判斷方法的應(yīng)用及分析

在數(shù)值模擬過程中,以吊艙推進器操舵角度范圍β=-30°～30°(順時針為負、逆時針為正)、操舵速率ω=10(°)/s、螺旋槳轉(zhuǎn)速n=3 r/s、冰層運動速度V=0.02 m/s、冰槳切削深度h=35 mm 為例進行分析,冰-吊艙推進器切削過程以及冰載荷變化曲線如圖17 和圖18 所示.

圖17 冰-吊艙推進器的切削過程(操舵狀態(tài),仰視圖)Fig.17 Milling process of ice-podded propeller(dynamic azimuthing condition,upward view)

由圖17 可知,當t=0 ms 時,吊艙推進器位于冰層中部,此時,吊艙推進器為直航狀態(tài),冰物質(zhì)點全部呈現(xiàn)為藍色,表明冰槳未發(fā)生切削.由圖17(a) 到圖17(b) 可知,在吊艙推進器逆時針操舵過程中,冰物質(zhì)點始終為藍色,表明螺旋槳與冰一直未發(fā)生接觸,主要原因是冰層逐漸靠近吊艙推進器,吊艙推進器逆時針操舵,使螺旋槳逐漸遠離冰層,導(dǎo)致冰層和螺旋槳之間始終保持一定的距離.當?shù)跖撏七M器達到最大斜航角度時,如圖17(b) 所示,隨著時間的增加,吊艙推進器開始順時針轉(zhuǎn)動,此時,螺旋槳逐漸靠近平層,因此,兩者在很短的時間內(nèi)發(fā)生切削作用,槳葉開展進入冰層,如圖17(c)所示,此時,冰物質(zhì)點變?yōu)榧t色,表明海冰發(fā)生了破碎,即螺旋槳與冰發(fā)生了切削作用.隨著吊艙推進器的操舵和冰層的運動,螺旋槳逐漸進入冰層,切削后的碎冰從槳葉后端拋出.隨著吊艙推進器操舵角度的變化,被槳葉高速拋出的碎冰與吊艙推進器的支架發(fā)生碰撞,引起吊艙推進器的冰激振動和瞬時高應(yīng)力,對吊艙推進器結(jié)構(gòu)的損害非常大,在冰區(qū)吊艙推進器設(shè)計過程中需特別注意.此外,在圖17(d)中可以清晰地觀察出吊艙推進器操舵完成后在冰層內(nèi)部留出的凹槽.

由圖18 可知,當操舵時間t<3000 ms 時,螺旋槳和吊艙單元的冰載荷都為0,說明此時間范圍內(nèi)冰和螺旋槳未發(fā)生接觸,與圖17(a) 到圖17(b) 的過程相一致.當3000 ms

圖18 吊艙推進器操舵過程冰載荷變化曲線Fig.18 Ice load curve of podded propulsor during dynamic azimuthing condition

5 結(jié)論

本文采用近場動力學(xué)方法和面元法耦合開展了不同操縱狀態(tài)的吊艙推進器與冰切削時接觸判斷方法的研究,通過對比吊艙推進器在直航、斜航以及操舵狀態(tài)與冰的切削過程,得出的主要結(jié)論如下:

(1)吊艙推進器在直航狀態(tài)與冰切削時,通過對海冰破碎過程、冰載荷變化曲線以及試驗值的對比和描述,證明了本文提出的冰-吊艙推進器接觸判斷方法的可靠性、模型建立的正確性以及計算方法的有效性.

(2)吊艙推進器在斜航狀態(tài)與冰切削時,通過對比吊艙推進器在不同斜航角度與冰切削時海冰的破碎過程、冰載荷的變化規(guī)律證明了吊艙推進器斜航狀態(tài)時接觸判斷方法的正確性;在吊艙推進器斜航狀態(tài)與冰切削時,建議吊艙推進器由直航狀態(tài)逆時針偏轉(zhuǎn)處于斜航狀態(tài),此時,與冰切削時吊艙推進器受到的冰載荷較小,對其結(jié)構(gòu)損壞程度也較小.

(3)吊艙推進器在操舵狀態(tài)與冰切削時,通過分析冰塊的破碎和冰載荷的時程曲線證明了吊艙推進器操舵狀態(tài)時接觸判斷方法的正確性;通過對海冰破碎過程的模擬發(fā)現(xiàn),槳葉高速拋出的碎冰更多地碰撞到艙體和支架上,將引起吊艙推進器的冰激振動或者瞬時高應(yīng)力,在冰區(qū)吊艙推進器設(shè)計過程中需特別關(guān)注.

本文重點分析了吊艙推進器處于不同操縱狀態(tài)時螺旋槳與冰的接觸判斷方法,未涉及到冰-艙體、冰-支架接觸判斷方法的研究,未來針對吊艙推進器的不同操縱狀態(tài)開展冰-艙體、冰-支架碰撞狀態(tài)時接觸判斷方法的研究.