陳明達 程細得

(武漢理工大學高性能船舶技術(shù)教育部重點實驗室1) 武漢 430063) (武漢理工大學交通學院2) 武漢 430063)

0 引 言

當船舶駛?cè)肽承┩ê浇ㄖ飪?nèi)的狹窄航道時，船舶的受力情況與在無限水域中有較大的差別.受阻塞效應影響船舶阻力值會增加，如果航道兩側(cè)相對船舶形狀不對稱，左右兩側(cè)流場不對稱會導致船舶受橫向力與首搖力矩的影響，因此，探究船舶在該類航道中航行時的水動力變化趨勢，對通航建筑物的設計、船舶航行的操縱及控制問題具有理論意義和實用價值.

在過去幾十年間，針對限制水域船舶水動力的研究有許多學者做了大量的模型實驗，數(shù)值計算與分析.早年由于計算機能力限制，相關(guān)的研究還是主要以基于勢流理論的細長體理論和三維面元法為主.Hess[1]針對近岸航行船舶的橫向力計算的提出了一套理論模型.熊新民等[2]采用三維Rankine源面元法計算了考慮自由面影響下船舶平行于岸壁航行時的水動力.后來計算機技術(shù)愈發(fā)成熟，計算機性能提高導致基于黏性流的CFD方法在限制水域船舶水動力研究上得到了廣泛應用.桑騰蛟等[3]以KVLCC2船型為對象,利用RANS方法計算其在不對稱岸壁及淺水中斜航水動力,研究不對稱岸壁、水深,以及其聯(lián)合作用對船舶斜航水動力以及水動力導數(shù)的影響.Hoydonck等[4]對比不同的CFD計算條件分析了自由液面、螺旋槳、流體黏性(粘流與勢流理論計算結(jié)果對照)和不同船岸距離與水深對船岸效應的影響.

盡管針對限制航道船舶水動力的研究開始較早，但是多數(shù)學者研究的還是理論上的規(guī)則均勻幾何航道，而通航建筑物內(nèi)的航道形狀往往是不均勻的.對此，采用基于RANS(reynolds average navier-stokes)的CFD方法，針對船舶直航駛?cè)氪l，應用流體分析商業(yè)軟件FLUENT對船舶的黏性繞流場進行數(shù)值模擬，并計算船舶的水動力.通過對比計算結(jié)果與實驗結(jié)果驗證數(shù)值計算的可行性，在此基礎上開展數(shù)值計算，分析船舶在航道的不同位置處水動力的變化趨勢，分析航道形狀對船舶進出船閘限制水域水動力的影響規(guī)律.

1 數(shù)學模型

1.1 控制方程

使用RANS的方法求解NS方程中的連續(xù)性方程和動量守恒方程：

(1)

(2)

用張量形式表達三個方向的動量方程，動量方程中速度的物理量都是時均值，帶有頂線的物理量是雷諾應力項，船舶流場的馬赫數(shù)遠小于臨界值，認為是流體不可壓縮的，故方程中省去了密度的偏導項.

1.2 湍流模型

湍流模型選用的是Menter[5]提出的兩方程SSTk-ω模型，其混合了k-e與k-ω模型的優(yōu)勢，是較為常用的湍流模型.選擇SSTk-ω模型理由有三：①張志榮等[6]認為在常用的湍流模型中，SSTk-ω模型是最適合船舶黏性流場數(shù)值模擬的湍流模型；②Lu等[7]在做了多組湍流模型的數(shù)值計算對比后之后認為SSTk-ω模型在限制水域船舶水動力計算中更具有優(yōu)勢；③由于在計算的中后半段船舶會駛?cè)胍暂^窄的航道，船舶附近的流場較為混亂復雜，故在近壁處理上，需要使用近壁面模型，SSTk-ω模型對于近壁剪切流計算親和性較好.為了保證y+=1 ，參考Frank[8]提供的公式得到第一層邊界層厚度，使船舶面網(wǎng)格附有5層較密的邊界層，邊界層的增長率是1.5，以捕捉船舶附近的流場信息.

1.3 自由液面

由于孟慶杰等[9]認為對于船舶在極為窄淺水域中航行時，自由面對船舶水動力的影響較大，所以雖然計算模型的船舶弗勞德數(shù)Fr在0.01～0.02，計算時仍考慮自由液面的影響.

VOF的方法就是通過研究單元網(wǎng)格內(nèi)各種填充介質(zhì)占網(wǎng)格總體積的比值F(體積分數(shù))來確定自由面[10].VOF方法的追蹤的是流體自由液面形狀的變化，而不是像勢流理論中的自由液面運動邊界條件一樣追蹤液面上質(zhì)點的運動.

(3)

VOF方法可以追蹤復雜的自由液面形狀、易實施、易拓展到三維空間.但是VOF只適用于壓力基求解器；使用VOF方法的算例的計算域必須充滿介質(zhì)，對于某些高速問題不適用；VOF方法對于在自由液面的捕捉需要高密度的網(wǎng)格支持，變相增加計算壓力.

2 數(shù)值求解

2.1 計算對象

選取12 000 TEU超巴拿馬集裝箱船進行計算，船模幾何模型側(cè)視圖見圖1，x方向相對位置坐標的參考點在首垂線處,主要參數(shù)見表1.

圖1 12 000 TEU幾何模型

表1 船模主要尺度

航道岸壁垂直水底，閘室長1.62Lpp,寬度為0.687 5 m.從閘門處向外航道左側(cè)伸出一道長1.32Lpp引墻，右側(cè)則是逐漸外擴，過渡的岸壁的x方向投影長度為0.44Lpp，引航道的寬度為2.725 m，船舶初始位置首垂線距離閘門1.48Lpp.航道幾何模型見圖2.

圖2 航道幾何模型

與隨船坐標系不同，固定坐標系計算域無來流速度，依靠動網(wǎng)格技術(shù)控制船舶做空間移動，故除TOP面需要設置為對稱面邊界條件(與氣相接觸，是否有切向物理量變化與計算結(jié)果影響不大)，其余壁面均是無滑移固壁(計算域尾段距離船尾初始位置有1倍LPP，理論上計算域尾段對船舶尾流無影響，只需要在該邊界上無質(zhì)量穿透，為了方便計算續(xù)性收斂，故選擇Wall邊界條件).

在航道中心線附近布置運動域，通過interfaces滑移網(wǎng)格與兩側(cè)靜止與進行數(shù)據(jù)交換，運動域分成三個部分，船舶近域網(wǎng)格隨船舶面網(wǎng)格同步運動，首尾方向的網(wǎng)格隨船舶運動漸變，具體分塊見圖3.

圖3 計算域劃分與網(wǎng)格重構(gòu)示意圖

剛體運動域與網(wǎng)格漸變域之間由interior面交界,見圖4，之間不需要產(chǎn)生網(wǎng)格的差值，與重疊網(wǎng)格相比，需要差值的面網(wǎng)格更少，相對的精度與計算穩(wěn)定性更高.

2.2 數(shù)值方法

數(shù)值計算空間離散采用的是有限體積法，空間上采用二階精度的中心差分格式.針對NS方程的求解壓力基求解器，計算過程中對壓力進行修正以滿足連續(xù)性和動量守恒，配合耦合隱式算法Coupled Implicit Solver(適用于非定常計算，全速度范圍求解，收斂性能好，但是占計算資源高)實現(xiàn)計算的時間迭代.

2.3 計算工況

Vantorre等[11]在比利時弗蘭德水動力研究中心的實驗公布的數(shù)據(jù)，針對其中的Test A與Test B船舶進閘工況的數(shù)值模擬，用以驗證CFD計算方法的準確性.由于實驗條件下未能使船模保持勻速，分析航道形狀變化對船舶水動力的影響時需要控制速度變量，故在Test A的基礎上設計一組勻速航行模擬Test C，見表2.

表2 三組不同的工況條件

3 計算結(jié)果分析

3.1 Test A結(jié)果分析

圖5 Test A計算結(jié)果與實驗數(shù)據(jù)對比圖

對比發(fā)現(xiàn)阻力的數(shù)值計算結(jié)果與實驗數(shù)據(jù)雖有一定差值，但是相差較小且總體趨勢一致，而橫向力與首搖力矩的計算結(jié)果差值較大.但是可以發(fā)現(xiàn)船舶橫向受力的計算結(jié)果與實驗數(shù)據(jù)初始分流位置明顯(x>-0.3)，此時船模大部分區(qū)域駛?cè)胍龎Ψ秶鷥?nèi)，有理由懷疑是船尾的某個變量與引墻相互作用干擾的數(shù)值計算的結(jié)果，由于數(shù)值計算未考慮螺旋槳，而實驗條件中船尾配有一定轉(zhuǎn)速的螺旋槳，可能是螺旋槳產(chǎn)生的尾流受引墻反射影響到船體進而影響船舶橫向受力.

Kaidi[13]在研究螺旋槳對船模的岸壁效應影響做了多組對照仿真，發(fā)現(xiàn)順時針旋轉(zhuǎn)的螺旋槳對單右側(cè)限制航道船舶水動力的影響較大.當量綱一的量的sbd(ship-bank dinstance)等于0.25時，螺旋槳的轉(zhuǎn)動嚴重影響船舶去流段與岸壁之間的流場壓力分布，使其出現(xiàn)一個明顯的高壓區(qū)，導致“岸推”現(xiàn)象加重.Test A引墻與船舶左舷的sbd為0.2，順時針旋轉(zhuǎn)的螺旋槳可能使得船舶去流段橫向受到一個額外向左的吸力，導致實驗橫向力結(jié)果偏小.為了進一步驗證計算方案的準確性，也為了證明當前的誤差分析的合理性，故針對Vantorre的Test B做一組數(shù)值仿真.

3.2 Test B結(jié)果分析

Test B與Test A的主要區(qū)別在于Test B的水深只有0.209 m，且Test B進行過程中大部分時間螺旋槳是停止轉(zhuǎn)動的，0轉(zhuǎn)速的螺旋槳相對于船舶主體來說只相當于一個小構(gòu)件附體，對水動力的影響較小，故Test B的數(shù)據(jù)結(jié)果能代表裸船體在進閘程中的水動力變化.

Test B數(shù)值計算結(jié)果與實驗數(shù)據(jù)吻合較好，見圖6.當x>0.6時船模速度劇烈震蕩導致阻力計算結(jié)果與實驗數(shù)據(jù)有一定出入，考慮到interface數(shù)據(jù)交換、網(wǎng)格拓撲能力以及幾何模型質(zhì)量等等不可避免的誤差，可以認為該方案的計算結(jié)果是符合物理規(guī)律的.

Test B實驗數(shù)據(jù)與計算結(jié)果的對比與在驗證了CFD計算方案的準確性的同時，也在一定程度上能證明Test A的計算誤差來源于實驗條件中螺旋槳轉(zhuǎn)動擾動流場帶來的影響.

圖6 Test B計算結(jié)果與實驗數(shù)據(jù)對比圖

3.3 Test C結(jié)果分析

Test A與Test B兩組計算證明了CFD仿真方案的準確性，Test C是勻速工況，船模速度取Fr相似條件下的服務航速Vs=0.11 m/s.

3.3.1阻力分析

Test C船模阻力分布圖見圖7.由圖7可知，當船舶位置x<0時阻力變化幅度較小，可以認為小阻塞比時單側(cè)岸壁對裸船體的阻力影響不大(側(cè)壁效應)，從船首部開始進入船閘時(x>0)船舶阻力才開始逐步增加.許立汀中間速度理論的船舶在限制航道內(nèi)受到的阻力增加，是因為收縮的岸壁使得船舶航道的阻塞率上升，船舶的回流速度增加等于變相增加船舶航速，但是觀察阻力計算數(shù)據(jù)發(fā)現(xiàn)，當航道收縮到一定程度，船舶阻力反而減小，說明阻塞率與船舶阻力并不是完全正相關(guān)，或者說還存在其他因素影響該航道內(nèi)的船舶阻力.

圖7 Test C船模阻力分布圖

由于航道與船舶的相對形狀是變化的，采用平均阻塞比的方法定義航道的阻塞效應，平均阻塞率為

(4)

式中：B為船寬；T為船舶吃水；ΔB為船舶首尾垂線平行范圍內(nèi)船岸平均距離；H為航道深度.

為了分析船舶在航行中阻力隨阻塞率的增大反而減小的原因，將船舶阻力按成分分離，研究航道形狀與阻塞率對船舶阻力的影響，見圖8.

圖8 航道不同位置阻塞率分布圖與Test C船模不同成分阻力分布圖

由于船舶的Fr(0.010時，黏壓阻力開始急劇上升；當x>0.7時，黏壓阻力逐漸下降.影響船舶黏壓阻力主要因素有中間速度和形狀因子，摩擦阻力隨船舶位置移動是呈穩(wěn)步增加的，說明船舶的實際速度是持續(xù)增長的，那么可以推出導致船舶黏壓阻力乃至船舶總阻力下降的關(guān)鍵因素就是船舶的形狀系數(shù).

限制航道船舶的形狀系數(shù)不僅與船舶自身形狀有關(guān)，還與航道相對船舶形狀有關(guān).黏壓阻力來源于船首尾的壓力差，在Test C的劇烈變化區(qū)域(x>0)，船舶前半段附近的航道形狀與后半段航道形狀有明顯差別.為了具體分析航道形狀對船舶阻力影響，在CFD-POST 中故將船舶在船中處將船舶分割，分別計算船舶前、后半段段船舶阻力的變化趨勢，見圖9.

圖9 Test C前、后半段船模x方向受力大小分布圖

由圖9可知船舶前半段船模x方向受力變化幅度大于后半段x方向受力變化幅度，即船舶總阻力主要受船舶前半段附近航道形狀控制.當x>0，船首開始進入船閘室時，船舶前半段受力開始增加，當x>0.5時，此時船舶前半段完全進入船閘，在之后的航行中船舶前半段左右兩側(cè)的航道形狀不在發(fā)生變化，而此時船模前半段x方向受力開始減小.需要注意的是，由于船閘室的長度(1.622Lpp)是有限的，當船舶逐漸駛?cè)氪l室時，會不可避免的受下級閘門(不可穿透邊界條件)的影響.

對于垂直船舶航行方向的壁面可能使船舶阻力的減小原因提一種解釋.從能量的角度的看黏壓阻力就是船舶單位時間將附進的流體向前推所需要做的功，當x越大，船舶離下級閘門越近，遠前端受船舶作用的流體的質(zhì)量就越小，則船舶對前方的流體單位時間內(nèi)需要做的功就越少，船舶所受的黏壓阻力就會減小.

在0.50.709后，船舶前半段x方向受力降速加快，而后半段x方向受力增速放緩，船舶阻力開始出現(xiàn)降低；當x>0.89后，雖然船舶前半段x方向受力開始趨向平穩(wěn)，此時而后半段x方向受力出現(xiàn)上升，使得船舶阻力仍是持續(xù)降降低.

3.3.2橫向受力分析

船模橫向力初始穩(wěn)定在-0.5處，當x>-1時，橫向力開始增加;當x>-0.75時，橫向力開始減??；當船舶完全進入引墻范圍時(x>-0.25)，船舶橫向力進入第一個波谷；之后橫向力開始出現(xiàn)急速的上升；船舶前半段完全駛?cè)腴l室時(x>0.5)，橫向力組件下降直至調(diào)轉(zhuǎn)方向，最后當船舶完全駛?cè)腴l室時，橫向力趨近于0，見圖10.

圖10 Test C船模橫向力分布圖

關(guān)于船舶橫向力變化趨勢分析，在后處理中分別對x=0.26，x=0.36與x=0.46處船模做iso surface分割處理，將船模表面按站位分割成20段，分別求解各個分段的橫向力數(shù)值，將結(jié)果匯總見圖11.

圖11 Test C船模不同位置處橫向力占比示意圖

由圖11可知，船模表面的橫向力會集中在相鄰的3～4個站上，且隨著船模進入船閘(x增加)，橫向力集中位置會相應的后移，且移動規(guī)律與船模運動規(guī)律反向且一致，說明船模橫向力的集中與船模相對航道的某一個固定位置，分析橫向力集中位置發(fā)現(xiàn)橫向力的集中區(qū)域出現(xiàn)在閘門與船模交界處，該位置可能是引起船模橫向力劇烈變化的關(guān)鍵.

通過分析x=0.46處z=0.185液面x方向速度云圖分布(圖12)，發(fā)現(xiàn)船舶右舷附近流場存在一處高回流速度區(qū)域.參考伯努利定理中對速度對壓力的影響，切向速度越大法相壓強越小，回流集中區(qū)域正是導致船舶右舷存在低壓區(qū)的原因.

圖12 x=0.46處z=0.185液面流場速度u云圖

船舶駛?cè)腴l室的同時閘室內(nèi)的水需要流出，當船模航行至處于，當船舶航行至x=0.46，閘門處斷面系數(shù)最大.閘門處左右舷斷面積相等，隨著流體質(zhì)點后移右舷斷面積逐漸大于左舷，由于流體的連續(xù)性，閘室內(nèi)的流體會傾向于從右舷流出，使得在閘門處右舷單位時間流過的流體體積大于左舷，即右舷的回流速度大于左舷，見圖13.

圖13 x=0.46處閘門位置剖面速度u云圖 (面向x正方向觀測)

船舶右舷的低壓是由船舶自身形狀與航道形狀共同作用，當船舶開始駛?cè)腴l室時(x>0)，時船舶局部橫剖面積變大，閘門處斷面系數(shù)增大，左右舷回流速度差增大，右舷低壓區(qū)更加明顯，導致船舶向右的橫向力增大，當船舶去流段開始于閘門位置接觸時，船舶局部橫剖面積減小，船舶橫向力減小.圖14為Test C船模首搖力矩分布圖

圖14 Test C船模首搖力矩分布圖

低壓區(qū)固定在于閘門處相近的船舶右舷，這意味著船舶的橫向力合力位置會隨著船舶航行慢慢從船舶前半段向后移動，這就可以解釋當x>0.18時船舶首搖力矩開始減小甚至出現(xiàn)反向.

4 結(jié) 論

1) 基于黏性CFD理論，研究了某船舶進入巴拿馬船閘的航行過程中的水動力進行了研究，對不同的物理條件進行了數(shù)值計算，并與實驗結(jié)果進行比較，驗證了本文計算方法合理性，為限制航道船舶操縱提供一定理論依據(jù).

2) 船舶在船閘內(nèi)航行時，會受左、右、前與底部岸壁的影響，使得船舶水動力發(fā)生變化，影響船舶航行的安全。對船舶而言，由于受航道阻塞率影響，船舶在駛?cè)氪l過程中會出現(xiàn)明顯的阻力增加現(xiàn)象，但是當船舶逐漸駛?cè)氪l時，下級閘門的存在會使導致船舶阻力降低；對船舶橫向受力而言，單側(cè)的岸壁對船舶橫向力的影響較小，但是航道形狀的劇烈變化會導致船舶橫向力的快速增加，且橫向力的受力中心會隨著船舶航行而后移，從而使得船舶在駛?cè)氪l過程中所受首搖力矩也會出現(xiàn)明顯的變化。

3) 仿真是基于文獻[14]的基礎上，針對其提出的展望，作出的一些計算上的調(diào)整，包括對于船模的速度的控制及考慮自由液面對水動力計算的影響，所以在計算精度上也有所提升，但是仿真仍未考慮船舶航行時的浮態(tài)變化.限制航道中船舶浮態(tài)變化大，對水動力結(jié)果有一定的影響，是下一步的研究方向.