基于船舶阻力計算方法確定浮泥的適航密度值

高志亮，龐啟秀，張瑞波

（交通運輸部天津水運工程科學(xué)研究所工程泥沙交通行業(yè)重點實驗室，天津300456）

基于船舶阻力計算方法確定浮泥的適航密度值

高志亮，龐啟秀，張瑞波

（交通運輸部天津水運工程科學(xué)研究所工程泥沙交通行業(yè)重點實驗室，天津300456）

效仿通過船模阻力試驗來界定浮泥適航密度值的方法，提出采用船舶阻力計算方法確定浮泥的適航密度值?；谇蠼釴S方程的CFD方法，建立了模擬浮泥中船舶運動流場的數(shù)學(xué)模型，模型中浮泥的本構(gòu)關(guān)系用改進的Herschel?Bulkley模型來描述。以連云港航道現(xiàn)場泥樣和超大型油船KVLCC2為例，采用所建立的數(shù)學(xué)模型計算了船舶以不同速度在不同密度浮泥中航行受到的阻力。通過數(shù)學(xué)手段分析了船舶阻力隨浮泥密度變化的趨勢，確定該港口航道浮泥的適航密度值為1 210 kg/m3，該值與采用流變試驗方法和船模阻力試驗方法得到的適航密度值接近，證明了文章提出方法的有效性。

船舶阻力；計算流體動力學(xué)；浮泥；適航密度

在淤泥質(zhì)河口和海岸地區(qū)，浮泥的發(fā)育常常導(dǎo)致航道的淤積，給船舶的正常通航帶來困難，為此需要定期對航道進行疏浚維護。為了減少航道的疏浚量和疏浚次數(shù)，適航水深技術(shù)被提出，即利用航道底部密度值較小的浮泥作為可航行水深。如何確定浮泥的適航密度，保證浮泥與船體接觸不會造成船舶結(jié)構(gòu)的損壞，或者嚴(yán)重影響船舶的操縱和控制，是應(yīng)用適航水深技術(shù)的關(guān)鍵。

長期以來，浮泥的流變試驗是確定其適航密度的實用化方法，被國內(nèi)外多個港口所采用［1-7］。該方法通過試驗測定不同密度浮泥的應(yīng)力和應(yīng)變關(guān)系，采用賓漢流體模型對試驗數(shù)據(jù)進行外推插值，得到零剪切速率對應(yīng)的應(yīng)力，作為浮泥的屈服剪切應(yīng)力，取密度與屈服剪切應(yīng)力關(guān)系曲線上曲率最大的點所對應(yīng)的密度值作為浮泥的適航密度值。然而對于密度值低于1 200 kg/m3的浮泥，其抵抗剪切的能力較弱，剪切應(yīng)力隨剪切速率變化的拐點出現(xiàn)得比較晚，采用賓漢模型外推方法得到的屈服剪切應(yīng)力存在較大的不確定性，從而影響浮泥適航密度值的準(zhǔn)確性。

為了彌補流變試驗方法的不足，我國頒布的《海港總體設(shè)計規(guī)范》規(guī)定浮泥適航密度值還應(yīng)通過船模阻力試驗的方法來確定，即在不同密度的浮泥中進行船模拖曳試驗［8-10］，測量船體以不同速度航行時所受的阻力，取浮泥密度與船模阻力關(guān)系曲線上曲率變化比較明顯的起始點所對應(yīng)的密度值作為浮泥的適航密度值。船模阻力試驗方法原理簡單，浮泥的流變特性由船舶的受力特點綜合反映，但該方法在實際應(yīng)用過程中存在著一些缺陷。首先，試驗需要有專門的場地和儀器設(shè)備，并進行費時費力的系列試驗，因而成本較高，效率較低；其次，泥樣配制過程繁瑣，且使用的浮泥量較大，因而很難保證每組試驗浮泥密度均勻一致，使得試驗結(jié)果的可靠性下降；再次，試驗中需要鋪設(shè)特定密度和厚度的浮泥，因此試驗一般在尺寸較小的水槽中進行，而水槽的邊壁效應(yīng)會對船舶阻力產(chǎn)生額外的影響，使得試驗結(jié)果不能完全反映浮泥特性的變化。

基于計算流體動力學(xué)（CFD）方法的數(shù)值船池技術(shù)［11-13］具有成本低、效率高、重復(fù)性好和無外界干擾等優(yōu)點。在清水中的船舶阻力預(yù)報領(lǐng)域，數(shù)值船池技術(shù)已經(jīng)跟模型試驗方法并駕齊驅(qū)［14-16］，將該技術(shù)代替船模阻力試驗方法來確定浮泥的適航密度值已經(jīng)成為可能。本文以連云港現(xiàn)場泥樣和超大型油船KVLCC2為例，應(yīng)用CFD方法計算船體在不同密度的浮泥中航行所受的阻力，并建立了船舶阻力與浮泥密度的數(shù)學(xué)關(guān)系，通過分析阻力隨密度變化的趨勢，確定該港浮泥的適航密度值。

1 數(shù)學(xué)模型

1.1控制方程

浮泥視為不可壓縮粘性流體，不考慮浮泥-空氣交界面波動和湍流影響，流體的運動由連續(xù)性方程和動量方程控制。在笛卡爾坐標(biāo)系下，控制方程描述如下

式中：腳標(biāo)i和j采用愛因斯坦求和約定，即相同腳標(biāo)表示從1至3求和，1、2、3分別代表空間方向的3個分量；x為空間坐標(biāo)；u為流體速度；ρ為流體密度；P為壓力；μ為流體表觀粘度。

浮泥為非牛頓流體，其表觀粘度隨剪切速率變化而變化。本研究采用改進型Herschel-Bulkley（HB）模型［17］描述浮泥的本構(gòu)關(guān)系，該模型表達如下

式中：τ為流體剪切應(yīng)力；τc為流體屈服剪切應(yīng)力；γ˙為流體剪切速率；γ˙c為流體臨界剪切速率；μd為流體動力粘度；n為冪指數(shù)。參數(shù)τc、γ˙c、μd和n通過浮泥流變試驗確定。

由于研究中的船速為中低速，故船體阻力成分將忽略興波阻力，而只計及摩擦阻力和粘壓阻力。當(dāng)流場中的流體速度和壓力求得后，船體所受阻力由公式（4）計算得到。

式中：R為船體阻力；SB為船體表面；n為船體表面法向量；腳標(biāo)i采用愛因斯坦求和約定；腳標(biāo)1代表空間沿船長方向的分量。

1.2數(shù)值方法

研究中采用通用CFD求解器FLUENT［17］對控制方程（1）和（2）進行求解。通過有限體積法在同位網(wǎng)格上對方程進行離散，其中對流項采用二階迎風(fēng)格式，耗散項和壓力項均采用中心差分格式。速度場和壓力場的耦合通過SIMPLE算法實現(xiàn)。

2 浮泥適航密度值確定

本節(jié)將敘述如何基于上述數(shù)學(xué)模型，采用船舶阻力計算方法代替船模阻力試驗方法來確定浮泥的適航密度值。研究中采用的浮泥取自連云港航道，計算船型為超大型油船KVLCC2。目前的計算機性能仍然難以對全尺度（縮尺比1:1）的船舶流場進行模擬，研究中船體的縮尺比為1:45.7。表1和圖1分別給出了該船的信息和型線。

首先，將現(xiàn)場取得的浮泥配制成不同密度（1 075～1 406 kg/m3）的泥樣，采用R/S-CC流變儀對配制好的泥樣進行流變測試，得到剪切應(yīng)力和剪切速率的關(guān)系。對試驗數(shù)據(jù)進行最小二乘擬合，求得HB模型（式（3））中的4個待定參數(shù)τc、γ˙c、μd和n，浮泥的本構(gòu)關(guān)系隨之確立。圖2和圖3給出通過流變試驗測得浮泥在不同密度下的應(yīng)力-應(yīng)變關(guān)系以及采用HB模型對試驗數(shù)據(jù)進行擬合后的結(jié)果，而模型參數(shù)列于表2?？梢钥闯?，應(yīng)力-應(yīng)變關(guān)系曲線上存在一個拐點，該點對應(yīng)的剪切速率γ˙即為臨界剪切速率γ˙c，γ˙c隨著浮泥密度增大而逐漸減小。當(dāng)γ˙＜γ˙c時，浮泥的剪切應(yīng)力隨其剪切速率的增大而急劇增大；當(dāng)γ˙＞γ˙c時，剪切應(yīng)力隨剪切速率變化變得平緩。從流變特性來看，浮泥隨著其密度增大從牛頓流體向賓漢流體轉(zhuǎn)變。以上所述的不同密度浮泥的流變特點均能被改進的HB模型很好地描述。

表1KVLCC2船型信息Tab.1Main particulars of KVLCC2

圖1KVLCC2船體型線Fig.1Body plan of KVLCC2

接下來采用所建立的數(shù)學(xué)模型計算KVLCC2油船在浮泥中航行受到的阻力。計算中船速取3種，分別為0.41 m/s、0.83 m/s和1.24 m/s，對應(yīng)的弗汝德數(shù)（Fn）分別為0.05、0.1和0.15。浮泥密度取12種，從1 075 kg/m3增加至1 406 kg/m3，相鄰密度差約為30 kg/m3。每次取1種船速和1種密度的浮泥進行計算，共36個計算組次。因為流場關(guān)于船舶中縱剖面對稱，故只取中縱剖面一側(cè)的流場作為計算區(qū)域，圖4給出了流場計算域范圍及網(wǎng)格劃分情況。計算域向上游、下游、側(cè)方和底部分別延伸1倍船長、2倍船長、1倍船長和10倍設(shè)計吃水的距離。計算時船體固定在流場中，在上游邊界施加速度入口條件，流速與船速大小相同方向相反；在下游邊界施加壓力出口條件；在上部邊界和中縱邊界上施加對稱條件；其余邊界均施加不可滑移固壁條件。計算網(wǎng)格劃分時，除了在船首球鼻和船尾槳軸附近區(qū)域采用非結(jié)構(gòu)化網(wǎng)格外，其余區(qū)域均采用結(jié)構(gòu)化的網(wǎng)格布置形式。船體表面附近網(wǎng)格尺寸為3 mm，遠離船體區(qū)域網(wǎng)格尺寸逐漸增大，流場網(wǎng)格總數(shù)約為71萬。

圖2低密度浮泥的應(yīng)力-應(yīng)變關(guān)系（ρ＜1 200 kg/m3）Fig.2Rheological behaviour of fluid mud with lower densities（ρ＜1 200 kg/m3）

圖3不同密度浮泥的應(yīng)力-應(yīng)變關(guān)系（ρ=1 100～1 400 kg/m3）Fig.3Rheological behaviour of fluid mud with different densities（ρ=1 100～1 400 kg/m3）

表2不同浮泥密度下HB模型參數(shù)Tab.2Parameters of HB model for fluid mud with different densities

圖4KVLCC2流場計算域及網(wǎng)格布置Fig.4Computational domain and mesh arrangement for KVLCC2

圖5不同航行條件下船舶阻力比較Fig.5Comparison of ship resistance in different navigation conditions

圖5給出了船體在不同速度和浮泥密度下航行的阻力計算值和采用六次多項式對數(shù)據(jù)進行擬合的結(jié)果?？梢钥闯觯白枇﹄S著航速和浮泥密度的增加而增大。在不同航速條件下，根據(jù)阻力隨浮泥密度變化的趨勢，可以將密度值分為（1 050～1 200 kg/m3）、（1 200～1 300 kg/m3）和（1 300～1 400 kg/m3）3個區(qū)間。當(dāng)密度小于1 200 kg/m3時，阻力隨密度變化比較平緩；當(dāng)密度值介于1 200 kg/m3和1 300 kg/m3之間時，阻力隨密度增長的比例開始變得明顯；當(dāng)密度值超過1 300 kg/m3后，阻力隨密度增大而迅速增加。

正如文獻［18］指出，通過觀察船舶阻力與浮泥密度關(guān)系曲線變化趨勢來確定浮泥適航密度值存在較大的人為因素；針對傳統(tǒng)方法在理論上的不足，本研究在如何利用船舶阻力與浮泥密度關(guān)系上提出了更為可靠的處理方法。前文已經(jīng)采用曲線擬合的方法得到阻力隨密度變化關(guān)系的數(shù)學(xué)表達式，因而可以通過數(shù)學(xué)計算的方法來得到阻力-密度曲線的曲率，然后取最大曲率所對應(yīng)的密度值作為浮泥的適航密度值。圖6給出了阻力-密度擬合曲線的曲率隨浮泥密度的變化，其中數(shù)值的正、負號分別表示曲線朝上或朝下彎曲?？梢钥闯觯?種不同航速下阻力-密度曲線的最大曲率所對應(yīng)的密度值非常接近，分別為1 219 kg/m3、1 210 kg/m3和1 216 kg/m3。取上述三者中的最小值1 210 kg/m3作為該流域浮泥的適航密度值，換算成重度值為11.9 kN/m3，該值與文獻［10］中通過浮泥流變試驗和船模阻力試驗綜合分析確定的適航重值12.2 kN/m3比較接近，證明基于船舶阻力計算方法確定浮泥的適航密度值是有效的。

圖6不同航速下阻力-密度曲線的曲率隨密度變化Fig.6Curvature of resistance?density curve changing in function of density for different forward speeds

3 結(jié)語

本文基于CFD方法建立了模擬浮泥中船舶運動流場的數(shù)學(xué)模型，并提出采用數(shù)值模擬手段代替物理模型試驗手段對船舶在浮泥中航行的阻力進行預(yù)報，將數(shù)值計算得到的船舶阻力用于浮泥適航密度值的確定。應(yīng)用本文方法對連云港航道浮泥的適航密度值進行了界定，研究結(jié)果表明：（1）數(shù)學(xué)模型中浮泥本構(gòu)關(guān)系的建立可以通過改進的HB模型并結(jié)合流變試驗方法實現(xiàn)；（2）采用數(shù)學(xué)手段分析船舶阻力與浮泥密度關(guān)系曲線可以避免傳統(tǒng)方法利用該曲線確定浮泥適航密度值的人為主觀性；（3）基于船舶阻力計算方法與基于流變試驗方法和船模阻力試驗方法確定的浮泥適航密度值相吻合；（4）船舶阻力計算方法具有廉價、高效、可靠等優(yōu)點，可作為現(xiàn)行浮泥適航密度確定方法的一種有益補充。本文僅采用了超大型油船作為計算船型，在下一階段的工作中將會對其他船型進行計算，以分析不同船型的選用對確定浮泥適航密度值的影響。

［1］陳學(xué)良.連云港浮泥測試及“適航深度”的確定［J］.水運工程,1998（8）:28-31.

［2］陳曉峰,黃建維,劉建軍.汕頭港外航道適航水深的研究和運用［J］.海洋工程,2002,20（4）:26-31. CHEN X F,HUANG J W,LIU J J.Investigation and utilization of the nautical depth of the outer navigation channel of Shantou Harbour［J］.The Ocean Engineering,2002,20（4）:6-31.

［3］錢平生.天津港適航水深下界面淤泥容重值確定與適航水深測量［J］.海洋測繪,2002,22（3）:35-39.

［4］梅劍云,孫月,李占元.國華臺電煤港適航水深應(yīng)用研究［J］.水道港口,2008,29（5）:314-317. MEI J Y,SUN Y,LI Z Y.Nautical depth applied to port of Guohua Taishan power plant［J］.Journal of Waterway and Harbor, 2008,29（5）:314-317.

［5］Wurpts R.15 years experience with fluid mud:definition of the nautical bottom with rheological parameters［J］.Terra et Aqua, 2005,99:22-32.

［6］Fontein W F,van der Wal J.Assessing nautical depth efficiently in terms of rheological characteristics［C］//Special Publication of the Hydrographic Society.Proceedings of the 15th International Congress of the International Federation of Hydrographic Societ?ies.London:Hydrographic Society,2006:149-152.

［7］Fontein W F,Byrd R W.The nautical depth approach,a review for implementation［C］//Randell R E.Proceedings of the World Dredging Congress XVIII Annual Seminar.Spotsylvania:Western Dredging Association,2007:767-772.

［8］張文忠.不同容重浮泥對船型阻力影響的物理模型試驗研究［D］.天津：天津大學(xué),2002.

［9］蔡南樹,龐啟秀,楊樹森,等.廣州港南沙港區(qū)港池適航水深綜合論證研究［J］.水道港口,2009,30（4）:253-256. CAI N S,PANG Q X,YANG S S,et al.Research of nautical depth in harbor basin of Nansha district of Guangzhou port［J］.Jour?nal of Waterway and Harbor,2009,30（4）:253-256.

［10］龐啟秀,楊樹森,楊華,等.淤泥質(zhì)港口適航水深技術(shù)研究與應(yīng)用［J］.水利水運工程學(xué)報,2010（3）:33-39. PANG Q X,YANG S S,YANG H,et al.Research and application of the technique of nautical depth in muddy harbors［J］.Hydro?Science and Engineering,2010（3）:33-39.

［11］劉應(yīng)中,張懷新,李誼樂,等.21世紀(jì)的船舶性能計算和RANS方程［J］.船舶力學(xué),2001,5（5）:66-84. LIU Y Z,ZHANG H X,LI Y L,et al.Calculations of the ship performances and solving of RANS equations in the 21st century［J］.Journal of Ship Mechanics,2001,5（5）:66-84.

［12］朱仁傳,繆國平,向紅貴,等.數(shù)值水池及其在船舶與海洋工程中的應(yīng)用［J］.上海造船,2007（4）:21-23.

［13］Stern F,Yang J,Wang Z,et al.Computational ship hydrodynamics:Nowadays and way forward［J］.International Shipbuilding Progress,2013,60（1-4）:3-105.

［14］Larsson L,Stern F,Bertram V.Benchmarking of computational fluid dynamics for ship flows:the Gothenburg 2000 workshop［J］. Journal of Ship Research,2003,47（1）:63-81.

［15］Carrica P M,Wilson R V,Stern F.An unsteady single?phase level set method for viscous free surface flows［J］.International Jour?nal for Numerical Methods in Fluids,2007,53（2）:229-256.

［16］倪崇本.基于CFD的船舶阻力性能綜合研究［D］.上海：上海交通大學(xué),2011.

［17］ANSYS Inc.ANSYS Fluent User′s Guide［R］.USA:Canonsburge,2010.

［18］金镠,虞志英,何青,等.淤泥質(zhì)港口航道適航密度確定方法的改進［J］.水運工程,2013（2）:91-95. JIN L,YU Z Y,HE Q,et al.Improvement of method for determining nautical density in muddy harbor［J］.Port&Waterway Engi?neering,2013（2）:91-95.

Determination of nautical density of fluid mud based on ship resistance calculation

GAO Zhi?liang,PANG Qi?xiu,ZHANG Rui?bo
（Tianjin Research Institute for Water Transport Engineering，Key Laboratory of Engineering Sediment,Ministry of Transport,Tianjin 300456,China）

Following the model test for ship resistance measurement,the method based on ship resistance cal?culation was proposed to determine the nautical density of fluid mud.The CFD method based on solving the NS equations was applied to calculate the resistance of ship advancing in fluid mud.The modified Herschel?Bulkley model was employed to describe the constitutive relation of fluid mud.To test its performance,the proposed method was applied to calculate the resistance of a KVLCC2 tanker in fluid mud dredged from Lianyun port.Based on the mathematical analysis of the changing trend between computed ship resistance and mud density,the nautical densi?ty for this port was defined as 1 210 kg/m3,which agreed well with the values defined by the rheology test and ship model test.

ship resistance;CFD;fluid mud;nautical density

U661.31+1

A

1005-8443（2015）04-0285-08

2014-10-11；

2014-11-24

中央級公益性科研院所基本科研業(yè)務(wù)費專項資金(TKS130201,TKS140201);留學(xué)人員科技活動項目擇優(yōu)資助經(jīng)費

高志亮(1982-),男,廣東省人,博士,助理研究員,主要從事計算船舶流體力學(xué)研究。

Biography：GAO Zhi?liang（1982-）,male,assistant professor.