陳 明，梁應(yīng)辰，宣國祥，陳明棟

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

0 引 言

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

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

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

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

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

圖1 船閘輸水系統(tǒng)Fig.1 Shiplock filling system

2 物理模型

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

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

3 數(shù)學模型

3.1 水流控制方程

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

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

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

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

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

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

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

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

3.3 船舶運動控制方程

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

3.4 數(shù)值計算

3.4.1 數(shù)值方法

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

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

3.4.2 網(wǎng)格剖分及邊界條件

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

圖4 船體表面網(wǎng)格Fig.4 Computation grid on the ship surface

圖5 網(wǎng)格剖分及邊界條件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 數(shù)值計算與模型試驗結(jié)果的比較

4.1 船舶運動過程

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

4.2 船舶系纜力

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

圖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上下振蕩，該變化特性與實測結(jié)果基本一致。在量值上，在個別時段內(nèi)實測值大于計算值，其余時段兩者均吻合良好。

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

5 算 例

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

圖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可知，布設(shè)輔助消能工前后縱向系纜力的過程線總體上未出現(xiàn)明顯的差異。從橫向系纜力過程線的對比情況（圖12）來看，布設(shè)輔助消能工前后前橫向系纜力和后橫向系纜力均基本沿FT=0軸上下振蕩，但兩者幅值的變化存在一定的差異，其中在輔助消能工布設(shè)前，前橫向系纜力的變化范圍為-1.82～1.63 kN，后橫向系纜力的變化范圍為-1.81～1.72 kN；輔助消能工布設(shè)后，前橫向系纜力的變化范圍為-1.23～1.25 kN，后橫向系纜力的變化范圍為-1.24～1.24 kN。由此可知，輔助消能工布設(shè)后較布設(shè)前降低了28%左右。以上計算結(jié)果表明，輔助消能工的功能主要體現(xiàn)在對閘室橫向流速均勻分布的調(diào)整上，為閘室船舶提供了較為均勻的流場環(huán)境，減小了船閘左右側(cè)的動水壓力差和局部水流力，從而達到了降低橫向系纜力的目的。

6 結(jié) 論

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

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

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

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

[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]劉應(yīng)中,繆國平.船舶在波浪上的運動理論[M].上海:上海交通大學出版社,1987.Liu Yingzhong,Liao Guoping.Theory of ship motion on the wave[M].Shanghai:Shanghai Jiaotong University Press,1987.

[3]肖 越.系泊系統(tǒng)時域非線性計算分析[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]鄒志利,張日向,張寧川等.風浪流作用下系泊船系纜力和碰撞力的數(shù)值模擬[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]王化明.限制水域操縱運動船舶粘性流場及水動力數(shù)值研究[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.