船閘輸水過程閘室船舶系纜力數值模擬

陳 明，梁應辰，宣國祥，陳明棟

（1重慶交通大學 國家內河航道整治工程技術研究中心，重慶400074；2中華人民共和國交通運輸部，北京100736；3南京水利科學研究院，南京 210029）

0 引 言

閘室船舶的停泊安全一直是船閘設計與管理最為關注的問題，而衡量船舶停泊條件的好壞，通常以船舶系纜力的大小予以判斷。研究閘室船舶系纜力，首先需考慮作用于過閘船舶錨泊系統上的水動力荷載。船閘輸水過程中閘室船舶與水體間的相互作用系“浮體—水流耦合動力響應問題，即船閘灌泄水時，閘室船舶錨泊系統將受水流作用，同時船體隨閘室水位的升降而作上下運動，反過來船體的運動狀態必然引起船舶周圍流場、壓力等參數的變化。對于船舶與水體間的耦合動力響應研究，有不少學者做過相關工作，如海上工作平臺在波浪作用下的動力響應[1]、浮體和系泊浮體在波浪作用下的動力響應[2-3]、高速水流沖擊浮體的運動響應[4-5]、錨泊在港的船舶系纜力計算[6]等。以上研究大多集中在波浪環境下浮體與流體間耦合作用的動力響應，這與船閘輸水過程“船舶—水流”間的耦合動力響應差異較大。原因在于，一方面船閘輸水時閘室內水位、流速都隨時間變化而變化，船舶將受非恒定流作用；另一方面，船舶在閘室中的運動屬大位移問題，且運動區域為限制性水域。因此，將已有的關于流固耦合動力響應研究方法直接運用于閘室船舶的系纜力計算存在較大困難。

Jong等[7]研究了閘室內涌浪對船舶的水動力作用及船舶的運動響應，并根據研究成果開發了閘室內船舶縱向力的計算程序，該程序假定了船舶所受縱向力是由三個分力的疊加，即推進波、均勻流和集中灌水形成的射流。Kalkwijk[8]建立了閘室船舶在水平面上的受力方程，但該方程僅適用于某一特定的船閘輸水系統。Natale等[9]針對閘門開孔的輸水系統、環形廊道輸水系統以及閘墻長廊道側支孔輸水系統，得出了閘室船舶含兩個自由度的振動方程。由于該方程考慮了錨泊系統的阻尼系數，但未給出系數的具體取值，Stockstill[10-11]通過物理模型試驗，分析討論了系數的計算公式，并建立了閘室系泊船舶的振動方程。限于閘室船舶受力過程的復雜性，上述振動方程的建立均假設了船舶浮力始終與重力平衡，這與船舶實際運動過程不相符。此外，船舶橫向受力亦是閘室停泊條件的重要研究內容，但上述振動方程只考慮了船舶的縱向受力。

本文基于“船舶—水流”間的耦合動力響應，建立船舶運動控制方程及縱向和橫向受力方程。利用FLUENT軟件，借助用戶自定義函數（UDF）編制閘室船舶系纜力并行計算程序，通過建立船閘整體輸水系統三維數學模型進行閘室船舶系纜力數值模擬研究，為快速、方便獲取閘室船舶的系纜力提供了可能。

1 模型試驗條件與計算條件

選取的船閘輸水型式為帶格柵消能室的短廊道集中輸水系統（圖1），該船閘閘室有效尺寸為140 m×14 m×2.5 m（長×寬×檻上最小水深），左右側廊道對稱布設，閥門最大工作水頭為9 m，閥門雙邊按7 min勻速開啟。結合船閘基本尺度，選擇停靠于閘室中央的300 t級單體船舶作為研究對象。鑒于實際船舶的體型繁多，研究時將船體進行一定的概化，其基本尺度為35.0 m×9.2 m×1.4 m（總長×型寬×設計吃水）。

圖1 船閘輸水系統Fig.1 Shiplock filling system

2 物理模型

模型按重力相似準則設計，比尺為1:36。模型設計總長約為10 m，包括上游水庫、上游引航道、輸水系統、閘室以及船舶。上引航道水位采用溢流式平水槽控制，輸水閥門采用可無級調速的步進電機驅動啟閉機控制，閘室水位變化過程利用電阻式點壓力傳感器測定。閘室船舶系纜力的量測，采用南京水利科學研究院研制的全環電阻式測試裝置。該裝置是先通過在船舶軸線上的首尾處分別選擇相應測點，然后在船首的測點處安裝縱向和前橫向全環電阻式測力儀，在船尾的測點處安裝后橫向全環電阻式測力儀。測力儀是由測力環、矩形鋼圈及滾輪組成，其中測力環上貼有應變片。測力儀安裝完畢后，將船舶安插在首尾的兩根垂向圓形立柱上，立柱與滾輪直接接觸并起到導向作用。當船舶受水流作用時滾輪首先受力，同時經矩形鋼圈將作用力傳至測力環上，隨之便引起測力環的變形，最后由采集系統對應變信號進行采集處理。物理模型總體布置及閘室船舶系纜力測試裝置如圖2所示。

圖2 物理模型布置及船舶系纜力測試裝置Fig.2 Layout of laboratory model and experiment facility for hawser forces

3 數學模型

3.1 水流控制方程

對船閘輸水過程產生的三維非恒定流的模擬選用RNG k-ε紊流模型，該模型在粘性水動力的研究中得到了廣泛應用[12-13]，其方程組在此不再贅述。

3.2 船舶縱向和橫向受力方程

在不受任何邊界約束的情況下，船舶在自由面上的運動為六自由度問題。而停靠在閘室中的船舶，為防在水流作用下發生撞墻、翻滾等海損事故，需通過纜繩將其系于浮式系船柱上以束縛船舶在水平面上的平動和沿豎向的轉動。船舶受到纜繩的約束后僅剩下三個自由度，即縱搖、橫搖以及沿豎向的平動。一般而言，閘室水面波動在輸水系統設計時需考慮一定的控制，因此縱搖和橫搖并不明顯，而參與最為顯著的運動為沿水面升降的平動過程。為簡化計算，建立船舶的受力方程時，只考慮船舶沿豎向平動過程的“船舶—水流”耦合動力響應。

船閘灌水時閘室船舶將受到水流和邊界約束的共同作用。其中，水流作用力可通過水壓力（含靜水壓力和動水壓力）和粘滯力來描述。為更好地實現對數值計算結果的驗證，本文采用的船舶約束邊界條件與物理模型試驗裝置一致。邊界約束力包括縱向約束力、前橫向和后橫向約束力。在此需說明的是，由于這三種約束力與船舶受纜繩的拉力作用效果相當，因而在模型實試驗中通常將縱向約束力、前橫向和后橫向約束力分別視為縱向系纜力、前橫向和后橫向系纜力。由此，便可得出船舶在水平面（XZ平面）上的受力情況（圖3）：船舶受到縱向系纜力FL、前橫向系纜力FT1和后橫向系纜力FT2、縱向水壓力Fpx和橫向水壓力Fpz、縱向粘滯力Fvx和橫向粘滯力Fvz。若將水流作用力均移至船舶質心時，則同時產生相應的力矩Mpy和Mvy。水流作用力的計算如下：

式中：px、pz分別為船舶任意表面上的水壓力在 x、z方向上的分量；As為船舶表面；τx、τz分別為船舶任意表面上的粘滯應力在x、z方向上的分量；（x0，y0，z0）為船舶質心坐標；（x，y，z）為船舶表面上任意一點坐標。

根據船舶在XZ平面上的受力平衡，可建立以下表達式：

在上式中，水流作用力可根據剖分的船體網格采用數值積分求出，則僅剩三個未知數，即縱向系纜力FL、前橫向系纜力FT1和后橫向系纜力FT2。因此，通過方程組便可求出船舶所受的縱、橫向系纜力。

圖3 船舶在水平面及豎直面上的受力圖Fig.3 All forces on horizontal and longitudinal section of vessel

3.3 船舶運動控制方程

船閘灌水時，船舶必然隨閘室水位的上升而上升，一旦閘室水位上升，則船舶浸水面積加大，勢必引起船體所受浮力增加，從而打破原有平衡狀態或改變原有加速度。因此，船舶總體上沿豎向（y方向）做變加速的上升運動。由于船舶與立柱間的摩擦力較小，研究中可忽略不計，由此可作出船舶在豎直方向上的受力圖（圖3）。根據牛頓第二運動定律，可建立相應的運動控制方程：

3.4 數值計算

3.4.1 數值方法

本文應用FLUENT大型流體力學計算平臺，對閘室船舶系纜力進行數值模擬。針對閥門開啟和船舶運動采用動網格模塊；對RNG k-ε紊流模型方程組采用控制體積法進行離散，壓力和速度的耦合求解利用SIMPLEC算法；對固壁邊界采用壁面函數法；閥門井及閘室自由液面采用水汽兩相流的VOF模型處理；對船舶的運動控制及系纜力的計算程序采用用戶自定義函數（UDF）進行編寫。為避免因船舶運動速度的過快而導致網格在重構時出現負體積，計算過程中對船舶運動速度及其運動方向上的水流作用力進行時均化處理，即將整個計算時間劃分為若干個運動時段，在每個運動時段內均作勻速運動，當前時段內的速度大小為上一時段內的速度值與作時均處理后的速度變化值之代數和，其計算公式見（6）式。考慮不影響計算結果的精度，本文選取每20個計算步為一個運動時段，則運動時段總數ns=tt/20，其中tt為總時間步數。

實現船閘輸水全過程的閘室船舶系纜力數值模擬，不僅需建立整體輸水系統三維數學模型，且要考慮工作閥門和船舶在同一模型中的多體運動。由此則帶來了計算網格數量龐大，同時還牽涉到大量的網格更新，普通PC機尚難于完成計算任務。為此，本文采用多處理器的工作站同時調用16個CPU進行并行計算，從而提高了計算效率，縮短了計算時間。

3.4.2 網格剖分及邊界條件

為兼顧計算精度和減小計算量，網格剖分時將閘室區域劃分成多塊區域，其中在閥門段、廊道進出口區域以及船舶動區域采用非結構網格劃分并進行網格加密處理，其余較規則的區域采用結構化網格進行劃分。船體表面網格見圖4。整個計算區域剖分的網格單元總數約為124萬個，節點總數約為36萬個，網格剖分見圖5。進口邊界施加沿水深分布的靜水壓力，兩側閥門井水體與大氣相通，采用空氣壓力進口，閘室出流同樣與大氣相通，采用空氣壓力出口，兩側工作閥門的開啟過程采用動邊界處理。

圖4 船體表面網格Fig.4 Computation grid on the ship surface

圖5 網格剖分及邊界條件Fig.5 Grid generation and boundary conditions

圖6 船舶運動速度變化Fig.6 Computed time series of vessel velocity

圖7 船舶位移與閘室水位比較Fig.7 Time series of curves of vessel displacement and water surface elevation

4 數值計算與模型試驗結果的比較

4.1 船舶運動過程

圖6為計算程序記錄的船舶運動速度隨時間的變化情況，在整個灌水過程中船舶的運動速度基本上是處于增減的交替變化狀態，即做變加速的上升運動。盡管速度的變化存在瞬態性，但從整體趨勢上看是呈先增加后減小的變化特性，在船閘灌水初期和末期上升速度較慢，在灌水中期上升速度較快，其中最大值達0.036 4 m/s。從船舶質心豎向位移曲線與閘室水位實測變化過程線間的比較（圖7）可以看出，兩條曲線整體上處于平行狀態，由此表明船舶上升過程與閘室水位的變化基本同步，這與船舶實際過閘過程一致。因此分析認為，通過采用動網格技術，基于“浮體—水流”耦合動力響應所提出的船舶運動數值方法是合理可行的。

4.2 船舶系纜力

圖8比較了船舶縱向系纜力數值計算與實測值的過程線，由圖可知，在灌水初期（0～170 s），兩者的變化特性和變化幅值吻合良好，均出現了3次明顯的周期性變化，且周期逐漸減小，其中在T=32 s時出現了第一次波谷，計算值為-4.9 kN，實測值為-4.0 kN；在灌水中期（170～395 s），縱向系纜力沿FL=0軸上下振蕩，仍然呈周期性變化且波形較為復雜；在灌水后期（395 s之后），實測值主要偏向于正向值，在該時段內，盡管兩者過程線的變化特性存在一定差異，但均具有逐步向零值靠攏的趨勢。縱向系纜力出現上述變化特性，原因在于灌水初期閘室水位較低，水流進入閘室時容易引起閘室水面的傾斜，從而激發水面周期性的長波運動，由此引起的船舶受力主要為正向（x正向）波浪力，則船舶縱向系纜力方向自然為x負向。在灌水中期船舶受水流局部力與波浪力的共同作用，此時存在正向波和反射波的多重疊加，從而導致了船舶受到了復雜的水流作用。之后，隨著輸水流量的下降，閘室斷面流速逐漸減小，閘室水面在縱向上以擺動為主，無明顯的水力坡降，船舶受力以水流局部力為主，且正向和負向的水流作用相當。到達灌水末期，水面波動趨于平穩，水流作用力逐步向零值衰減。

圖8 縱向系纜力計算與實測過程線Fig.8 Computed and observed time series of longitudinal hawser force

圖9 船舶橫向系纜力計算與實測過程線Fig.9 Computed and observed time series of transverse hawser forces

圖9比較了船舶橫向系纜力的計算值與實測值的過程線，由圖可知，計算得出的船舶前橫向系纜力和后橫向系纜力的變化曲線均基本在軸FT=0上下振蕩，該變化特性與實測結果基本一致。在量值上，在個別時段內實測值大于計算值，其余時段兩者均吻合良好。

由以上對比分析可知，縱向和橫向系纜力的實測值在個別時段內大于計算值。分析其原因，一方面是由于試驗過程中啟閉機運行時偶爾發出較強的電磁信號干擾了應變信號的采集質量，另一方面是由于建立船舶受力方程時只考慮了船舶沿豎向平動的單自由度運動所致。但綜合上述對比成果，計算結果在整體上達到了較好的精度要求，尤其是在船閘灌水初期，船舶縱向系纜力的吻合程度較高。因此，本文提出的數值模擬方法是合理可行的，可用于閘室船舶系纜力研究。

5 算 例

基于上述計算條件，針對格柵消能室內有無設置輔助消能工進行閘室船舶系纜力數值模擬研究。輔助消能工的布設方法是在消能室內對稱布置縱向消力檻和階梯式消力梁，同時在每個頂面格柵孔下方設置“T”形橫截面的縱向消力擋板，見圖10。由于停泊位置對船舶的受力影響較大，本次計算采用的停泊位置為閘室上半區域且靠右閘墻。

圖10 輔助消能工三維立體圖Fig.10 3D stereogram of auxiliary energy dissipaters

圖11 船舶縱向系纜力過程線比較Fig.11 Comprison of computed and observed time series of longitudinal hawser force

圖12 船舶橫向系纜力過程線比較Fig.12 Comprison of computed and observed time series of transverse hawser forces

對于集中輸水而言，船舶的主要受力階段在船閘灌水初期，圖11～12給出了船閘灌水前200 s的船舶系纜力過程線。由圖11可知，布設輔助消能工前后縱向系纜力的過程線總體上未出現明顯的差異。從橫向系纜力過程線的對比情況（圖12）來看，布設輔助消能工前后前橫向系纜力和后橫向系纜力均基本沿FT=0軸上下振蕩，但兩者幅值的變化存在一定的差異，其中在輔助消能工布設前，前橫向系纜力的變化范圍為-1.82～1.63 kN，后橫向系纜力的變化范圍為-1.81～1.72 kN；輔助消能工布設后，前橫向系纜力的變化范圍為-1.23～1.25 kN，后橫向系纜力的變化范圍為-1.24～1.24 kN。由此可知，輔助消能工布設后較布設前降低了28%左右。以上計算結果表明，輔助消能工的功能主要體現在對閘室橫向流速均勻分布的調整上，為閘室船舶提供了較為均勻的流場環境，減小了船閘左右側的動水壓力差和局部水流力，從而達到了降低橫向系纜力的目的。

6 結 論

本文針對帶格柵消能室的船閘集中輸水系統，以300 t級單船為研究對象，開展了灌水過程的閘室船舶系纜力數值模擬研究，得出以下主要結論：

（1）基于“船舶—水流”耦合動力響應，建立了船舶縱向和橫向受力方程及其運動控制方程。通過充分利用FLUENT軟件提供的用戶自定義函數功能，實現了閘室船舶系纜力的并行計算。將計算結果與模型試驗結果對比發現，本文提出的數值計算方法具有良好的精度，可用于閘室船舶系纜力研究，為下一步考慮更為全面和更為復雜的船舶三自由度運動的系纜力數值模擬打下了基礎；

（2）由船舶系纜力過程線的變化規律可知，船舶縱向受力主要發生在灌水初期，此階段的水流作用力主要由坡降力組成，船舶橫向上的水流作用力主要由動水壓力差和局部水流力組成。因此，對于集中輸水系統而言，應采取工程措施控制好灌水初期的水面坡降和調整閘室流場的均勻分布，以達到減小船舶系纜力的目的；

（3）通過對消能室內有無設置輔助消能工情形下的閘室船舶系纜力數值模擬，研究表明在消能室內采用消力檻、消力梁和消力擋板的組合式布置方式，可有效降低船舶橫向系纜力。研究成果對類似工程的設計和建設具有重要的指導意義。

[1]Korsmeyer F T,Lee C H,Newann J N,et al.The analysis wave interactions with tension leg platforms[C]//Houston:OMAE Conference,1988.

[2]劉應中,繆國平.船舶在波浪上的運動理論[M].上海:上海交通大學出版社,1987.Liu Yingzhong,Liao Guoping.Theory of ship motion on the wave[M].Shanghai:Shanghai Jiaotong University Press,1987.

[3]肖 越.系泊系統時域非線性計算分析[D].大連:大連理工大學,2005.Xiao Yue.Non-linear time-domain anaylsis and calculation for moored systems[D].Dalian:Dalian University of Technology,2005.

[4]Sigalotti L Di G,Lopez H,Trujillo L.An adaptive SPH method or strong shocks[J].Journal of Computational Physics,2009(228):5888-5907.

[5]肖 瀟,蔣昌波,程永舟.水流對浮體作用的SPH方法模擬[J].船舶力學,2011,15(8):861-866.Xiao Xiao,Jiang Changbo,Cheng Yongzhou．Simulation of flow-induced floating-body motion with SPH method[J].Journal of Ship Mechanics,2011,15(8):861-866.

[6]鄒志利,張日向,張寧川等.風浪流作用下系泊船系纜力和碰撞力的數值模擬[J].中國海洋平臺,2002,17(2):22-27.Zou Zhili,Zhang Rixiang,Zhang Ningchuan,et al.The numerical simulation of mooring force&impact force of ship moored to offshore platform[J].China Offshore Platform,2002,17(2):22-27.

[7]de Jong R J,Vrijer A.Mathematical and hydraulic model investigation of longitudinal forces on ships inlocks with door filling system[R].Delft Hydraulic Laboratory,1980.

[8]Kalkwijk J P T.Hydrodynamic forces and ship motions induced by surges in a navigation lock[R].Dept.of Civil Engineering,Delft Univ.of Technology,1973:73-1.

[9]Natale L,Savi F.Minimization of filling and emptying time for navigation locks[J].Journal of Waterway,Port,Coast.,and Oc.Engrg,2000(126):274-280.

[10]Stockstill R L.Mooring model coefficients for barge tows in a navigation lock[J].Journal of Waterway,Port,Coastal and Ocean Engineering,2003:233-236.

[11]Stockstill R L.Modeling hydrodynamic forces on vessels during navigation lock operations[J].ASCE,http://www.ascelibrary.org,2004.

[12]Yakhot V,Orzag S A.Renormalizationn group analysis of turbulence:Basic theory[J].J Scient Comput,1986,1:3-11.

[13]王化明.限制水域操縱運動船舶粘性流場及水動力數值研究[D].上海:上海交通大學,2009.Wang Huaming.Numerical study on the viscous flow and hydrodynamic forces on a manoeuvring ship in restricted waters[D].Shanghai:Shanghai Jiaotong University,2009.