淺水中低速斜航船舶水動力預報及驗證與確認分析

鄒璐

（上海交通大學 船舶海洋與建筑工程學院 高新船舶與深海開發裝備協同創新中心，上海 200240）

淺水中低速斜航船舶水動力預報及驗證與確認分析

鄒璐

（上海交通大學 船舶海洋與建筑工程學院 高新船舶與深海開發裝備協同創新中心，上海 200240）

淺水中的斜航船舶受到淺水阻塞效應和不對稱流的綜合影響。為預報該運動中的船舶水動力，文章采用基于定常雷諾平均納維—斯托克斯方程的計算流體動力學方法，對淺水中做斜航運動的船舶粘性繞流場進行數值模擬。考慮低航速運動的特點，忽略航速影響下的自由面興波，由數值計算得到水動力系數在漂角影響下的變化規律。針對計算精度問題，在數值模擬中從驗證和確認角度分析和評估計算結果：通過網格收斂性分析分析數值誤差與不確定度；結合試驗數據考察計算模型的誤差。此外，從計算區域尺度、湍流模型、邊界條件、船體下沉和縱傾作用方面對模型誤差的影響因素進行探討，可為改進計算模型、提高數值模擬精度提供參考依據。

淺水；斜航運動；CFD；水動力系數；誤差和不確定度；驗證和確認

0 引 言

如何準確預報船舶操縱水動力一直以來是船舶操縱性研究領域的熱點問題，因為這與船舶運動模擬和操縱控制密切相關。目前主要的預報方法包括：經驗公式估算、約束模試驗以及基于計算流體動力學（Computational Fluid Dynamics,CFD）的數值模擬。近年來，基于CFD的數值計算方法得到越來越廣泛的應用，已成為與約束模試驗并駕齊驅、相輔相成的一種水動力預報方法。然而，該方法不可避免地存在精度問題。對數值計算中的誤差和不確定度的準確評估以保證數值模擬的可靠性是CFD數值預報的重要環節。目前驗證和確認（Verification and Validation）是能有效評估數值和模型誤差的方法，國際上有許多相關研究進展[1-5]。

實際工程中船舶常出現斜航運動的情況。此時，船體的縱軸線與航行方向成一定夾角，即所謂的漂角。船體因受到不對稱流的影響，船尾出現較直航情況下愈加明顯的流體分離，船尾伴流變大。在淺水中斜航的船舶由于受到阻塞效應和不對稱流的綜合作用，船體受到更復雜的水動力影響。淺水斜航船舶水動力預報在國際上是比較研究的代表案例。“船舶操縱預報方法驗證與確認專題討論會”SIMMAN2008[6]對此設置了相關基準算例。意大利船模水池INSEAN（Italian Ship Model Basin）對該算例進行了平面運動機構約束模試驗（Planar Motion Mechanism，PMM），為數值模擬提供了重要的對比研究用數據。國外的Simonsen等[7]、Toxopeus等[8]和國內的田喜民等[9]、錢永峰[10]采用基于雷諾平均納維—斯托克斯方程（Reynolds Averaged Navier-Stokes，RANS）的CFD方法均對淺水斜航問題做過相關研究，但這些研究并沒有系統討論計算中的誤差和不確定度。作者前期工作[11]中已對淺水斜航有一定研究，本文主要是在驗證和確認基礎上更深入地探討數值計算的精度水平、量化評估誤差和不確定度，以及分析可能的誤差源。

1 物理問題的描述

本文參考SIMMAN2008專題討論會中的基準算例，選取淺水中在固定漂角下低速斜航運動的超大型油船船模KVLCC2（KRISO Very Large Crude-oil Carrier）為研究對象。圖1表示的是KVLCC2斜航運動的示意圖，β為漂角。本文著重考察的是縱向力X，橫向力Y和轉首力矩N，計算結果以無量綱形式表示為：

圖1 斜航運動示意圖Fig.1 Illustration of the drift motion

圖2 KVLCC2船型及淺水域幾何示意圖Fig.2 Schematic of the flow field(h:water depth)

算例中船后不帶附體（槳/舵）；船模縮尺比為45.714，對應主尺度：垂線間長LPP=7.0 m，船寬B= 1.269 m，吃水T=0.455 m。船型幾何如圖2所示。根據INSEAN的試驗設置，計算中水深吃水比h/T= 1.2；低船速U=0.533 m/s，對應的雷諾數Re=U·LPP/ν=3.697×106，佛汝德數（ν為水的運動粘度，g為重力加速度）。此外，選取四個漂角β=0°，2°，4°，6°進行系列計算。由于航速較低，假定自由面興波及船體自由下沉和縱傾對船舶水動力的影響可忽略不計，本文采用疊模方法模擬靜水面并保持船體正浮。

2 數值方法

假定繞斜航KVLCC2船模的粘性流是不可壓縮的，相應的流動問題可由如下控制方程描述：

本文采用船舶水動力學專業軟件SHIPFLOW中的RANS求解器XCHAP，對物理問題展開數值計算。XCHAP基于有限體積法，用Roe格式[12]離散對流項，中心差分離散擴散通量。離散的方程由交替方向隱格式（alternating direction implicit scheme）[13]求解。壓力和速度通過人工壓縮格式（artificial compressibility scheme）[14]全耦合。為封閉控制方程，XCHAP中主要有兩個湍流模型可用：剪切應力輸運模型（Shear Stress Transport k-ω model，k-ω SST）[15]和顯式代數應力模型（Explicit Algebraic Stress Model，EASM）[16]。

3 計算設置和網格生成

數值計算中的坐標系定義為固定在船上的右手笛卡爾直角坐標系：坐標原點置于船中剖面、中縱剖面和水面的交點，如圖3所示。x，y，z軸分別指向船首、右舷、船底。由于船體周圍流場不對稱，數值求解選取的計算區域包含整個船體，由七個邊界面組成：船體表面；靜水面（z=0）；淺水底的位置根據計算條件位于水面下1.2T處；計算域入口在船首柱（F.P）前的一倍船長LPP處；出口在船尾柱（A.P）后的LPP處；兩邊的計算域邊界面（即岸壁處）距船中心線均為1.8LPP。在計算域的各個邊界面上分別設置邊界條件：船體表面滿足無滑移條件（No-slip），無壁面函數；入口/出口邊界面上設置入流/出流條件（Inflow/Outflow），即入口處流場速度等于船速，出口處除壓力外的所有的流動參數法向梯度為零；水面處由于不考慮興波影響，采用滑移條件（Slip）模擬為對稱面；淺水底及岸壁同樣設置滑移條件以簡化計算。計算中首先采用EASM模型進行湍流模擬。

為保證數值求解中網格離散的質量，本文采用結構化重疊網格技術對計算域進行網格離散。由圖4可見，重疊網格由兩部分組成。首先對船體及其附近的區域用半圓柱型的H-O型貼體網格離散，在船首尾附近流場變化比較突出的區域加密網格以提高求解的精度。但是很明顯這種網格形式不適合離散船身下方平直的淺水底，因而在H-O網格外部用嵌套的H-H網格來離散平直幾何邊界（特別是水面、岸壁和水底）附近的區域。

圖3 坐標系及邊界條件Fig.3 Manoeuvring coordinate system and boundary conditions

圖4 網格離散示意圖Fig.4 Schematic of grid distribution

4 驗證和確認（Verification and Validation）

4.1 驗證：數值誤差和不確定度分析

驗證（Verification）是用來評估數值誤差和不確定度。數值求解過程中的誤差由計算機的舍入誤差、計算迭代誤差和對數學方程的離散誤差組成，其中占主導的部分為離散誤差。離散誤差的產生是源自數值求解需要對數學方程進行空間和時間上的離散。理論上隨著離散點的增加，離散誤差應該趨于零，然而實際數值計算中的離散誤差很難避免，需要通過所謂的“收斂性分析（Convergence study）”來評估離散問題帶來的誤差和不確定度。對于定常問題，則主要在于網格收斂性分析。本文采用E?a等[4]的最小二乘（Least Squares Root，LSR）理論對一系列不同密度的網格進行分析和評估。

為具體考察本研究中的數值誤差和不確定度，以h/T=1.2，β=4°的斜航運動為例進行網格收斂性分析。根據E?a等[4]的方法，網格收斂性分析需要通過加密網格來考察數值誤差和不確定度。首先選定網格尺度比值（hi+1和hi為兩相鄰網格的步長），按照這個網格加密比例，生成6個密度不同的網格。表1列舉了這6個網格密度下的網格節點數，其中以hi/h1表示網格尺度比以區別網格的疏密程度，h1對應的是最密的網格。然后，對這些系統加密的網格進行數值計算得到相應的水動力系數預報結果。

表1 網格尺度Tab.1 Dimensions of grid series

圖5 X′，Y′，N′的網格收斂性Fig.5 Grid convergence of X′,Y′,N′

應用LSR方法對數值結果進行曲線擬合，得到網格收斂趨勢圖。網格5給出了各網格密度下的X′，Y′，N′數值計算結果和曲線擬合得到的網格收斂趨勢及其相關參數（外推值S0，精度p），其中也給出了理論精度pth（在數值方法中已知，本文為二階精度）對應的擬合結果用作比較。需要說明的是，網格6由于網格過于稀疏給出的結果不精確，因此其數值解沒有包含在曲線擬合當中。由圖5可見，擬合得到的精度p都大于理論的pth=2.0，特別是Y′，原因很可能是其各個數值解呈現一種分散（scatter）的形式，使得Y′值的變化與網格疏密程度無明顯相關性。而X′和N′的結果圖顯示，隨著網格的加密，二者的數值解都趨于收斂，相鄰網格數值解之差也逐漸變小。

表2列出了由LSR方法得到的不確定度USN，其中選取最后加密的三個網格（網格3，網格2，網格1）的結果以作比較。由網格3到網格1，三個水動力系數的不確定度并不明顯隨網格密度變化，結合圖5也可看到三個網格密度對應的數值解較為接近。然而表2中數值不確定度的值仍然較大，一方面是由于本研究中用到的重疊網格某種程度上將網格離散復雜化，因而使其數值求解精度受限；另一方面對本文中的p?2.0情況，LSR方法的不確定度評估公式的準確度還有待考察。由網格收斂趨勢圖（圖5）可知，高密度網格能給出較為精確的數值解，但同時也需要更多的計算資源。在結合數值不確定度結果并權衡數值精度和計算時間的基礎上，下文選擇對網格3的結果展開進一步討論，并用該網格進行系統的淺水斜航計算以考察漂角對船舶水動力系數的影響。

表2 X′，Y′，N′的數值不確定度Tab.2 Numerical uncertainties of X′,Y′,N′

表3 X′，Y′，N′的確認分析結果Tab.3 Validation results of X′,Y′,N′

4.2 確認：模型誤差的評估

以上討論的是在網格離散過程中的數值誤差和不確定度，為了更深入地探討誤差問題，還需對數值建模過程中的模型精度作進一步分析，即評估模型誤差和不確定度。本文采用美國機械工程師協會（American Society of Mechanical Engineers，ASME）發表的V&V 20-2009標準[3]中的方法，結合INSEAN的試驗數據（包括試驗不確定度），針對同一工況h/T=1.2，β=4°進行確認分析。在ASME標準的確認分析方法中，引入了確認不確定度和比較誤差兩個參數。S為數值解，USN為其對應的數值不確定度；D為試驗數據，UD為對應的試驗測量不確定度。若，說明模型誤差隱含在由數值不確定度和試驗不確定度帶來的“噪聲（noise）”里，因而無法對模型誤差做出評估；而若，則表示建模過程存在模型誤差，需要結合比較誤差E的大小和符號來改進模型。表3列出了本文研究中水動力系數X′，Y′，N′的確認不確定度和比較誤差。X′，Y′，N′三者的平均誤差在11%左右，與SIMMAN2008專題討論會[6]上其他機構研究結果的精度水平一致。比較兩個參數可以看到，X′和N′的確認不確定度都遠大于比較誤差，意味著模型誤差并不明顯；Y′的十分接近，比較誤差稍大，所以還需要進一步的考察來確定模型可能存在的問題。

4.3 模型誤差的來源分析

對淺水斜航KVLCC2船模的水動力預報中，可能存在的模型誤差源有以下幾個方面：

a.計算區域縱向尺度的影響；

b.湍流模型的影響；

c.水底和岸壁的簡化邊界條件的影響；

d.忽略船體自由下沉和縱傾的影響。

為檢查每一誤差源，本文逐一進行模型改進和新的數值計算。首先為考察計算區域縱向尺度的影響，將x方向尺度增大，把入口和出口邊界面都設置在船前和船后2LPP處。而y和z方向由于要保證和模型試驗的淺水、岸壁條件一致，所以保持不變。然后，用Menter的k-ω SST模型代替之前采用的EASM模型來判別湍流模型的影響。此外，原計算中水底和岸壁的邊界條件采用滑移條件簡化計算，而實際上斜航船體對水底/岸壁存在相對運動，在水底面（受船體擾動影響更大）上會形成邊界層，因此在底/岸上采用移動壁面（無滑移）條件。最后，低速下的自由面興波往往很小，對水動力的影響不大，所以將船體浮態的改變和自由面效應一并忽略。但是在淺水條件下，由于其阻塞作用，對船體浮態（下沉和縱傾）可能有一定影響，而下沉和縱傾的改變也會引起水動力的變化。再者，模型試驗中船體是可以自由下沉和縱傾的。為此，將初始計算得到的下沉和縱傾看作流場中船舶調整后的姿態，將其加入新的計算當中來模擬船體自由下沉和縱傾對水動力的影響。

表4 各誤差源對應的Uval和Tab.4 Uvalandof different sources of modelling error

表4 各誤差源對應的Uval和Tab.4 Uvalandof different sources of modelling error

X′ Y′ N′22.53 1.63計算區域縱向尺度2LPP?原結果Uval%D E%D 26.31 12.29 19.78 20.27 22.50 1.81湍流模型（k-ω SST）Uval%D E%D 25.60 9.29 20.00 21.59 Uval%D E%D 25.76 9.89 18.57 12.39 23.13 1.08底/岸邊界條件初始下沉和縱傾誤差源總和Uval%D E%D Uval%D E%D Uval%D E%D 26.76 14.26 26.60 13.55 25.73 9.78 20.54 25.16 20.49 24.89 20.27 23.45 22.78 0.47 24.22 6.03 25.15 10.25

對以上四個誤差源分別進行調整和重新計算后，得到確認所需參數，如表4所列。與原結果對比可見，將計算區域的縱向尺度加大后，縱向力系數X′的比較誤差降低了3%左右，而橫向力和轉首力矩的誤差變化不大；三個水動力系數的對比結果并無改變。湍流模型k-ω SST對X′和N′結果的影響與計算區域尺度類似，但是將Y′的誤差降低了近8%并使其比較誤差小于確認不確定度，在某種程度上降低了模型誤差的影響。此外，底/岸的移動岸壁邊界條件增大了X′和Y′的誤差，卻稍微降低了N′的誤差，原因還有待查證。而加入初始下沉和縱傾量后，也使得三個水動力系數的比較誤差較原計算有所增大。最后，在計算中綜合四個誤差源改進后，可以看到本質上并沒有改變原確認結果，但是Y′和N′的確認不確定度和比較誤差都有所增加，尤其是N′。

4.4 漂角對水動力系數的影響

在以上分析的基礎上，對其他漂角情況進行類似的計算（分別采用原計算設置和綜合改進四個誤差源的設置），得到相關水動力系數，如圖6所示。

計算得到的X′，Y′，N′均隨β的增加而變大，與INSEAN的試驗數據相比可以看到Y′，N′的變化趨勢與試驗一致，但是X′結果和試驗數據的趨勢不同。在前期深水情況的數值計算中[11]，X′的計算結果與INSEAN試驗值的對比也存在類似問題，但計算結果與其他比較試驗結果一致。另外，本文X′計算結果與Toxopeus等人[8]的也相符，所以考慮X′變化趨勢的誤差與試驗中測量的誤差和不確定度有關。Y′，N′變化趨勢的數值預報結果雖然較好，但是數值上的精度仍然不夠高。另一方面，圖6顯示綜合改進前文所述的四個誤差源后，得到的結果與原計算結果的趨勢一致。數值上X′和N′較原結果都有所增大；Y′值在β=0～4°上的改變并不明顯，但在β=6°處比原值要大。

圖6 X′，Y′，N′隨漂角的變化Fig.6 X′,Y′,N′against drift angle

5 結 語

本文采用基于定常RANS方程的數值計算方法預報了在淺水中低速斜航運動船舶的水動力系數。為了把握數值計算中的誤差和不確定度以便為提高計算精度提供必要的參考信息，文中采用系統的驗證和確認方法，評估和分析計算中的數值誤差以及模型誤差。研究結果顯示本文預報的淺水中水動力系數變化趨勢與試驗數據一致，但計算的精度問題仍然存在，平均誤差在11%左右。在確認分析中，縱向力和轉首力矩與試驗的比較誤差比數值和試驗的不確定度之和要小，所以不易判別模型誤差；而橫向力系數的比較誤差與確認不確定度非常接近甚至前者稍大，所以存在某種模型誤差。文中考慮了四種可能的誤差源：計算區域縱向尺度、湍流模型、底/岸邊界條件、初始下沉和縱傾，分析結果表明這四個誤差源對水動力系數的影響各不相同，唯一將橫向力誤差降低的是湍流模型k-ω SST；而綜合所有的誤差源并沒有直接提高橫向力和轉首力矩的精度，可見確定模型誤差源的復雜性，在今后的研究中還需更多的探討。

[1]Roache P J.Verification and validation in computational science and engineering[M].Hermosa Publishers,Albuquerque, 1998.

[2]ITTC,7.5-03-01-01.CFD general:Uncertainty analysis in CFD-verification and validation methodology and procedures [K].Quality Manual,2002.

[3]ASME V&V 20-2009.Standard for verification and validation in computational fluid dynamics and heat transfer[S].American Society of Mechanical Engineers,New York,UK,2009.

[4]E?a L,Vaz G,Hoekstra M.A verification and validation exercise for the flow over a backward facing step[C]//Proceedings of the V European Conference on Computational Fluid Dynamics,ECCOMAS CFD 2010.Eds.Pereira J C F,Sequeira A,Lisbon,Portugal,2010.

[5]Zou L,Larsson L.A V&V study based on resistance submissions to the Gothenburg 2010 workshop on numerical ship hydrodynamics(Chapter 5,Book name:Numerical Ship Hydrodynamics-An Assessment of the Gothenburg 2010 Workshop. Editors:Larsson L,Stern F,and Visonneau M)[C].Springer,ISBN:978-94-007-7188-8,2013.

[6]SIMMAN 2008 Workshop[C].Available at:http://www.simman 2008.dk,Copenhagen,Denmark,2008.

[7]Simonsen C D,Stern F,Agdrup K.CFD with PMM test validation for manoeuvring VLCC 2 tanker in deep and shallow water[C]//Proceedings of the International Conference on Marine Simulation and Ship Manoeuvrability,MARSIM 2006. Terschelling,Netherlands,M-4-1,2006a.

[8]Toxopeus S L,Simonsen C D,Guilmineau E,Visonneau E G,Xing T,Stern F.Investigation of water depth and basin wall effects on KVLCC2 in manoeuvring motion using viscous-flow calculations[J].Journal of Marine Science and Technology, 2013,DOI 10.1007/s00773-013-0221-6.

[9]田喜民,鄒早建,王化明.KVLCC2船模斜航運動粘性流場及水動力數值計算[J].船舶力學,2010,14(8):834-840. Tian Ximin,Zou Zaojian,Wang Huaming.Computation of the viscous flow and hydrodynamic forces on a KVLCC2 model in oblique motion[J].Journal of Ship Mechanics,2010,14(8):834-840.

[10]錢永峰.淺水中作斜航運動船體粘性繞流計算[D].華中科技大學,2007.

[11]Zou L,Larsson L,Orych M.Verification and validation of CFD predictions for a manoeuvring tanker[J].Journal of Hydrodynamics,Ser.B,2010,22(5):438-445.

[12]Roe P L.Approximate riemann solvers,parameter vectors,and difference schemes[J].Journal of Computational Physics, 1981,43:357-372.

[13]Chakravarthy S,Osher S.A new class of high accuracy TVD schemes for hyperbolic conservation laws[J].AIAA Paper, No.85-0363,1985.

[14]Chorin A.A numerical method for solving incompressible viscous flow problems[J].Journal of Computational Physics, 1967,2:12-26.

[15]Menter F R.Zonal two equation k-ω turbulence models for aerodynamic flows[C]//Proceedings of the 24th Fluid Dynamics Conference.Orlando,AIAA paper 93-2906,1993.

[16]Gatski T B,Speziale C G.On explicit algebraic stress models for complex turbulent flows[J].Journal of Fluid Mechanics, 1993,254:59-78.

Numerical predictions of hydrodynamic forces on a ship during a low-speed drift motion in shallow water including verification and validation

ZOU Lu
(School of Naval Architecture,Ocean and Civil Engineering,Shanghai Jiao Tong University,Collaborative Innovation Center for Advanced Ship and Deep-Sea Exploration,Shanghai 200240,China)

Travelling at a drift angle in shallow water,the ship is affected by both the asymmetry flow and the blockage from shallow water.To predict the hydrodynamic forces on the ship,this paper applies a Computational Fluid Dynamics method based on the steady state Reynolds-Averaged Navier-Stokes equations to solve the viscous flow around the ship during a low-speed drift motion in shallow water.Considering the ship moves at a low speed,the effects of free surface are assumed negligible.Tendencies of the hydrodynamic forces against the drift angle are simulated from numerical computations.For the purpose of evaluating the degree of accuracy in the numerical results,verification and validation analyses are performed.The numerical error and uncertainty are estimated from a grid convergence study,while the modelling errors are investigated along with the experimental data.In addition,the contributions of computational domain size, turbulence model,boundary condition,as well as the influence of sinkage and trim to the modelling errors are discussed.This study offers a useful reference for the improvement of the numerical model and the ac-curacy in numerical predictions.

shallow water;drift motion;CFD;hydrodynamic forces;error and uncertainty; verification and validation

U661.1

：Adoi:10.3969/j.issn.1007-7294.2016.07.007

1007-7294（2016）07-0841-08

2016-01-20

國家自然科學基金項目（51309152）；高性能船舶技術教育部重點實驗室開放基金課題資助項目（2013033101）

鄒 璐（1983-），女，博士，助理研究員，E-mail：luzou@sjtu.edu.cn。