劉培林，王墨偉，冷文浩，金建海

（1.江南大學 數(shù)字媒體學院，江蘇 無錫 214122；2.無錫職業(yè)技術學院，江蘇 無錫 214121；3.中國船舶科學研究中心，江蘇 無錫 214082）

基于SPH流體仿真的船模阻力試驗過程的虛擬展現(xiàn)

劉培林1，2，王墨偉3，冷文浩3，金建海3

（1.江南大學 數(shù)字媒體學院，江蘇 無錫 214122；2.無錫職業(yè)技術學院，江蘇 無錫 214121；3.中國船舶科學研究中心，江蘇 無錫 214082）

通過求解RANS方程計算船體周圍粘性邊界層，同時結(jié)合VOF方法處理自由表面，能夠獲得帶自由表面的船體周圍穩(wěn)態(tài)粘性流場，但是難于比擬拖曳水池模型試驗得到連續(xù)變化的流場，該文提出采用SPH流體仿真方法模擬船模阻力試驗中自由表面及流場的生成與演化，并借助三維交互式圖形顯示技術實現(xiàn)試驗過程的虛擬展現(xiàn)，最后以DTMB 5415船模試驗加以例證。

SPH流體仿真；船模阻力試驗；虛擬展現(xiàn)

0 引 言

隨著全球船舶水動力學CFD（Computational Fluid Dynamics，計算流體力學）應用技術的快速發(fā)展，數(shù)值水池虛擬試驗技術再次成為國際水動力學領域創(chuàng)新研究的核心。中國船舶科學研究中心的科研團隊針對新技術提出數(shù)值水池應具備三大技術特征：基于屬性細分的知識封裝、基準試驗檢驗與大子樣應用驗證及結(jié)果的可信度評估和情景化[1]。其中，情景化的基本內(nèi)涵是比擬物理水池模型試驗，借助虛擬現(xiàn)實技術手段，對試驗對象、環(huán)境、過程和結(jié)果等進行虛擬重建，構(gòu)造并展現(xiàn)高逼真度、高精度的水動力性能試驗“現(xiàn)場”，讓用戶在虛擬實驗室中隨意漫游，并與試驗對象等進行交互，從而能夠從不同位置、不同視角觀察試驗過程中連續(xù)變化的流場，如阻力試驗中的自由表面、分離流動、漩渦等，從而獲得“身臨其境”的用戶體驗。

自由表面模型的數(shù)值模擬方法通常分為兩類：擬合法和捕捉法[2]。在擬合法中，網(wǎng)格產(chǎn)生變形以匹配自由表面形狀，該方法非常適合計算穩(wěn)態(tài)流問題；在捕捉法中，自由表面位于流體域網(wǎng)格內(nèi)部，無需根據(jù)網(wǎng)格變形確定自由表面，該方法適合處理變化較為劇烈的自由表面，常見方法包括VOF（Volume of Fluid，流體體積）法和Level Set（水平集）法。一般說來，通過求解RANS（Reynolds Averaged Navier-Stokes，雷諾平均納維—斯托克斯）方程計算船體周圍粘性邊界層，同時結(jié)合VOF方法處理自由表面，能夠獲得帶自由表面的船體周圍穩(wěn)態(tài)粘性流場[3]，但是難于比擬拖曳水池模型試驗得到瞬時連續(xù)變化的流場。因此，常規(guī)CFD流體仿真方法難以提供船模阻力試驗過程虛擬展現(xiàn)所需的瞬態(tài)流場數(shù)據(jù)。

本文依托江蘇省科技廳資助的研究課題，采用SPH（Smooth Particle Hydrodynamics，光滑粒子流體動力學）流體仿真方法對船模阻力試驗過程中不斷變化的自由表面及流場進行數(shù)值模擬，并借助三維交互式圖形顯示技術實現(xiàn)試驗過程的虛擬展現(xiàn)，嘗試復現(xiàn)阻力試驗過程中自由表面及流場的生成與演化情況，為數(shù)值水池未來提供遠程服務、滿足用戶隨時隨地實施虛擬試驗的需求奠定基礎。

1 SPH流體仿真應用現(xiàn)狀

流體仿真的數(shù)值方法主要有兩類：歐拉網(wǎng)格法（Eulerian Grid-based Method）和拉格朗日粒子法（Lagrangian Particle-based Method）。網(wǎng)格法的理論和算法都比較成熟，但在處理三維復雜問題，尤其是涉及到大變形和不連續(xù)流動時（例如自由表面砰擊、液艙中的晃蕩沖擊、船舶甲板上浪等伴有液面翻轉(zhuǎn)和破碎等強非線性運動），會因為網(wǎng)格發(fā)生嚴重變形、扭曲而嚴重影響計算精度，采用自動重構(gòu)網(wǎng)格和減小時間步長等改進措施則會大幅增加計算量[4]。SPH方法作為一種無網(wǎng)格粒子法，其粒子分布具有自適應性，在處理移動邊界問題，尤其是大變形、強沖擊、自由表面等問題時具有明顯優(yōu)勢[4]。國內(nèi)外學者利用SPH流體仿真方法在潰壩問題、大壩溢流、海岸防護、海冰模擬、水下爆炸沖擊、模擬人類游泳等領域已經(jīng)開展了卓有成效的研究[5-9]。

在船舶水動力學領域，SPH流體仿真方法近年來也逐漸得到了應用。Cartwright等[10]采用SPH方法研究了規(guī)則波中的船體運動問題，結(jié)果表明SPH方法能夠精確地預測船體在規(guī)則波中的運動；Le Touzé等[11]采用 SPH-flow軟件預測了惡劣海況下的甲板上浪問題，二維模型給出了合理的預測結(jié)果；Daniel John Veen[12]采用二維SPH方法研究了規(guī)則海浪中的船艏抨擊問題，數(shù)值結(jié)果與試驗得到的抨擊壓力基本一致；胥飛等[13]采用SPH方法模擬具有代表性的二維外飄型船體物體入水問題，給出了艏外飄型船體的入水抨擊數(shù)值模擬結(jié)果，并與試驗結(jié)果和邊界元方法計算結(jié)果相比較，結(jié)果表明在低速、剖面曲率變化較小的情況下，單相SPH方法可有效地模擬二維外飄型船體入水抨擊過程；Shao等[14]為研究粘性不可壓縮液體晃蕩問題提出了一種改進的SPH方法，包括密度修正及核梯度修正、引入湍流模型、改善固體邊界處理等，并將數(shù)值計算結(jié)果與試驗觀察和其它來源結(jié)果進行比較，結(jié)果吻合很好。

2 經(jīng)典流體仿真方法

在計算機圖形學中，流動現(xiàn)象的模擬是一個重要問題，NS（Navier-Stokes，納維—斯托克斯）方程是流體仿真的理論基礎。目前在流體仿真中應用最為廣泛的是基于歐拉網(wǎng)格法的實現(xiàn)，但大多數(shù)程序只能離線工作，不能解決流體的動態(tài)交互問題[15]。

2.1 NS方程

等溫粘性不可壓縮流體隨時間t運動的經(jīng)典公式服從NS方程[15]：

其中：ρ為質(zhì)量密度，u為速度，p為壓力，μ為流體粘度，f為作用在流體上的外力的總和。▽符號表示梯度算子，▽·符號表示散度算子，·符號表示點積運算。（1）式為流體的動量守恒定律，（2）式為不可壓縮流體的質(zhì)量守恒定律。

2.2 歐拉流體

歐拉流體可看成由規(guī)則排列的流體單元格組成，每個單元格包含若干流體分子或粒子[15]。歐拉法描述了部分流體屬性，例如質(zhì)量密度和壓力場。由于受網(wǎng)格約束，歐拉法難以描述流體的自然流動。

3 SPH流體仿真基礎

3.1 SPH基本原理

SPH是一種采用離散樣本點對連續(xù)場的值和導數(shù)進行近似的插值方法，樣本點是帶有質(zhì)量、位置、速度等屬性的光滑粒子[15]。SPH采用核近似方法將偏微分方程轉(zhuǎn)換成積分方程，然后用粒子近似方法將連續(xù)形式的積分方程轉(zhuǎn)換成離散形式的方程[4]。對于空間中任意一點r處的任意物理量A的值，可通過對其支持域內(nèi)粒子離散求和的方法得到[15]：

其中：mj為粒子j的質(zhì)量，ρj為粒子j的質(zhì)量密度為一個徑向?qū)ΨQ的光滑函數(shù)（也稱作核函數(shù)，可簡寫為Wij），h為光滑長度（也稱作光滑核半徑，它決定了支持域的大小）。

3.2 SPH粒子近似在拉格朗日粒子法中，用粒子代替網(wǎng)格可顯著簡化NS方程的形式。假定粒子總數(shù)不變，每個粒子的質(zhì)量保持不變，（2）式可忽略。粒子完全定義了流體，隨流體一起移動。拉格朗日形式的NS方程定義如下[15]：

其中：右邊第一項為壓力項，第二項為粘性力項，二者共同構(gòu)成流體內(nèi)力。第三項為外力，用于平衡內(nèi)力，包括重力、浮力等。各物理量的近似定義如下[15]：

（1）質(zhì)量密度：粒子質(zhì)量是用戶定義的常量，質(zhì)量密度的計算僅依賴粒子質(zhì)量：

（2）壓力：流體內(nèi)部粒子的壓力由理想氣體定律決定，同時考慮牛頓第三定律：

（3）粘性力：流體內(nèi)部分子摩擦產(chǎn)生粘性力，降低流體動能，μ為粘性系數(shù)：

（4）重力：重力均勻地作用在所有流體粒子上，g為向下的重力加速度：

（5）浮力：浮力由溫度擴散引起。對于等溫流體，人工浮力可定義為：

其中：b為人工浮力擴散系數(shù)。當質(zhì)量密度小于閾值時，粒子產(chǎn)生浮力。

3.3 SPH形式運動方程

采用SPH方法模擬流體運動時，須滿足拉格朗日形式的流體動力學方程[5]。

3.3.1 動量方程

連續(xù)場的動量守恒方程定義為：

其中：Γ為耗散項。針對耗散影響，可以采用人工粘性、層流粘性、層流粘性和SPS湍流幾種方式進行處理，人工粘性算法簡單且健壯性好，應用較為廣泛。本文采用人工粘性，（10）式可近似為：

3.3.2 連續(xù)性方程

連續(xù)場的質(zhì)量守恒方程定義為：

通過求解質(zhì)量守恒方程得到密度變化。（13）式可近似為：

對于弱可壓縮流體，為了減小密度和壓力波動，可參考delta-SPH模型[16]引入一個數(shù)值耗散項[17]（δ為delta-SPH系數(shù)，通常取0.1）：

3.3.3 狀態(tài)方程

為了根據(jù)粒子密度確定壓力，須將流體看作是弱可壓縮介質(zhì)。壓力計算公式定義為：

4 SPH流體仿真數(shù)值算法

SPHysics是一個由歐美研究人員開發(fā)的可求解自由表面流動問題的開源SPH數(shù)值計算程序，Du-alSPHysics是它的開源并行版本[17]。目前國外已將DualSPHysics應用于潰壩、晃蕩、造波、主動消波、防波提、浮體運動、振蕩水柱（一種波浪發(fā)電裝置）、系泊錨鏈等復雜問題的數(shù)值模擬，并開展了模型試驗驗證[18]。

本文采用DualSPHysics計算程序，通過建立自由表面流動的數(shù)值計算模型，對繞船體自由面粘性流問題進行數(shù)值模擬，生成展現(xiàn)精細流場所需的隨時間變化的流場數(shù)據(jù)。

4.1 光滑核函數(shù)選擇

光滑核函數(shù)應具備歸一性、緊支性、非負性、單調(diào)遞減性、Dirac（狄拉克）δ函數(shù)性質(zhì)、對稱性、光滑性等特性[8]。針對特定問題選擇合適的光滑函數(shù)，對算法的效率和精度有著非常重要的影響。常用的光滑核函數(shù)有高斯型、三次樣條型和五次樣條型[5]，高斯核函數(shù)足夠光滑，對無序粒子來說非常穩(wěn)定精確，但在理論上不滿足緊支性條件，且計算量較大[8]。三次樣條核函數(shù)形似高斯核函數(shù)，計算量相對較小，目前應用最為廣泛。五次樣條核函數(shù)在計算精度和時間成本方面能夠達到最佳折衷。DualSPHysics提供了Cubic Spline三次樣條核函數(shù)和Quintic五次樣條核函數(shù)[17]。三次樣條核函數(shù)能滿足本文仿真要求，定義如下：

其中：q=r/h。

4.2 時間步算法選擇

粒子相互作用使得每個時間步中速度、密度、位置等物理量均發(fā)生變化。將SPH運動方程中動量、密度、位置的導數(shù)分別寫成如下形式：

DualSPHysics提供了兩種時間積分算法，其中Verlet算法計算簡單，Symplectic算法穩(wěn)定性好但計算量較大[17]。本文主要使用了Verlet算法。該算法基于兩個泰勒展開式。算法被拆成兩部分，通常情況下按照下式計算：

其中：n表示時間步，t=nΔt。每隔Ns個時間步（通常取50），則按照下式計算：

算法的第二部分用于防止積分值隨時間推移產(chǎn)生發(fā)散。當出現(xiàn)數(shù)值穩(wěn)定性問題時，可適當減小Ns。當Ns小于10時，表明Verlet算法不再適用，此時應采用Symplectic算法。

在忽略摩擦或粘性影響的情況下，Symplectic算法具有時間可逆性，因此能提高長時間求解的精度。顯式二階Symplectic算法包括預估和校正兩個階段。

在預估階段，即時間步的中間，計算密度、速度和位置的值：

4.3 實體邊界處理

當粒子位于實體邊界附近時，支持域被邊界截斷，（3）式不再成立[4]。一般采用虛粒子法來解決邊界問題，將邊界看作是一組獨立于流體粒子的粒子集[5]。DualSPHysics提供實體不可穿透性和周期性開放邊界，以及根據(jù)固定力作用而移動的邊界[17]：

（1）動態(tài)邊界條件

這類邊界粒子與流體粒子滿足相同的方程，但不隨作用其上的力一起移動，相反它們要么是位置固定的，要么是按照一個施加的運動函數(shù)移動。當流體粒子靠近邊界時，邊界粒子的密度增加，導致壓力增加，從而產(chǎn)生一個作用于流體粒子上的排斥力。

（2）周期性開放邊界條件

這類邊界條件是通過允許靠近開放側(cè)邊界的邊界粒子與在計算域的另一側(cè)的靠近互補開放側(cè)邊界的流體粒子的相互作用來實現(xiàn)的。粒子的支持域被鄰近的開放邊界截斷，支持域的剩余部分應用于互補的開放邊界。

（3）預施加的邊界條件

通過預定義的運動函數(shù)或提供時間軸的輸入文件來定義邊界粒子的運動。這類邊界粒子具有與動態(tài)邊界條件相似的行為，但是它們是不固定的，獨立于作用其上的力而移動。邊界粒子運動時會影響流體粒子，提供了定義復雜仿真行為的能力（例如造波槳）。

4.4 SPH數(shù)值修正

為了避免流體粒子局部集中或速度太大導致穿透實體邊界等問題，需要對SPH流體仿真數(shù)值算法進行一些工程化修正，以使其能夠較為真實地反映流動情況。

4.4.1 速度修正算法

為了減少速度波動，粒子移動的速度被修正為相鄰粒子的速度平均值[17]：

其中：ε為常數(shù)，取值范圍為0～1.0。通過計入粒子速度以及所有相互作用的粒子的平均速度，重新計算粒子速度。速度修正可使粒子分布更有條理，避免高速粒子穿透邊界。

4.4.2 密度修正算法

采用常規(guī)SPH仿真方法，粒子的壓力場會呈現(xiàn)不穩(wěn)定性。解決該問題的最直接且計算量最小的方法是密度再初始化，即對粒子密度進行濾波和賦值。當粒子靠近邊界時，支持域內(nèi)的粒子數(shù)不足，需要使用修正的核函數(shù)初始化粒子的密度[5]。DualSPHysics提供了Shepard濾波器，可對密度場進行修正。每隔Ns個時間步（通常取30），采用下式進行修正[17]：

5 流場虛擬展現(xiàn)數(shù)值算法

為了獲得較為精細的流場數(shù)據(jù)，需要定義較小的粒子間距和時間步長。在計算域大小不變的情況下，粒子間距越小，粒子數(shù)目就會越多。復雜的三維SPH流體仿真模型往往包含了千萬以上的流體粒子，很難直接用于渲染顯示。因此就需要對流體粒子進行采樣，以降低流體粒子的空間分辨率，實現(xiàn)流場數(shù)據(jù)壓縮，從而方便進行實時渲染，以實現(xiàn)流體的透明和景深效果。最后還應參照自然界中水的真實效果加入反射、折射、焦散等光照效果。

本文采用DirectX應用程序接口編寫流場虛擬展現(xiàn)程序SphVizKit，嘗試虛擬重建試驗對象和試驗環(huán)境，并以三維實時交互方式展現(xiàn)船模阻力試驗過程中的自由表面及流場的生成與演化情況。

5.1 流場數(shù)據(jù)壓縮

將計算域沿X、Y、Z方向等分為Cx，Cy，Cz份。根據(jù)每個時間步中流體粒子的位置，建立流體粒子與格子的多對一映射關系：

其中：dp為粒子間距。將格子內(nèi)部的流體粒子的位置、密度、速度、壓力等數(shù)值進行求和、取平均，使得每個格子僅對應一個流體粒子。

5.2 流體立體渲染

SPH流體仿真既非網(wǎng)格法，也非水平集法，本身不產(chǎn)生自由面，只能得到粒子的位置、密度、速度、壓力等數(shù)值，因此渲染基于粒子的流體是比較困難的[19]。流體粒子渲染的常規(guī)做法是采用Marching Cube（行進立方體）算法[20]，即根據(jù)三維密度場提取多邊形網(wǎng)格等值面，這種表面重建方法非常耗時，僅適用于非實時渲染場合。本文采用Screen Space Fluid Rendering（屏幕空間流體渲染）方法[21]，可直接在二維屏幕空間中渲染基于粒子的流體，無需產(chǎn)生多邊形網(wǎng)格，僅處理視錐體以內(nèi)的流體粒子，其余不可見部分則被剔除掉。該方法完全基于現(xiàn)代圖形硬件，因此非常適合實時交互式渲染。算法描述如下：

第一步，生成深度圖像：將流體粒子渲染成球狀的Point Sprite（點精靈，一種簡單高效的粒子系統(tǒng)渲染技術，常用于模擬爆炸、煙霧、火焰等效果），得到包含流體表面像素的深度圖像。

（1）在Vertex Shader（頂點著色器）中將三維空間中的流體粒子轉(zhuǎn)換為二維空間中的粒子球，其坐標定義為四元數(shù)：

其中：MWorld-View-Proj為4×4轉(zhuǎn)換矩陣。

（2）在Pixel Shader（像素著色器）中計算每個粒子球的實際位置及深度值：

其中：MView-Proj為4×4轉(zhuǎn)換矩陣。

第二步，光滑深度圖像：為了得到光滑的流體表面，需要對深度圖像進行模糊化處理。采用Bilateral Filter（雙邊濾波器，一種保邊去噪的濾波器）算法[22]，對平滑變化區(qū)域進行模糊處理，同時保留輪廓邊界。雙邊濾波器包含兩個權重函數(shù)，一個由空間距離決定，另一個由色彩距離決定：

其中：σ為權重參數(shù)，x表示空間距離或色彩距離。

第三步，計算表面法向：為了渲染流體表面，需要得到光滑的表面法向。根據(jù)視點與相鄰像素點的深度距離的偏差計算法向量。

第四步，生成不透明度圖像：為了讓流體看起來更真實，在根據(jù)深度圖像生成流體表面的同時，還需要計算流體的不透明度，并根據(jù)不透明度對顏色進行衰減，從而渲染出半透明的流體表面。顏色衰減的實現(xiàn)參考Bouguer-Lambert-Beer（布格—朗伯—比爾）光吸收基本定律：

其中：I0為入射光的強度，k為吸收系數(shù)，d為介質(zhì)厚度。

最后，生成最終圖像：合并光滑后的深度圖像、不透明度圖像和背景圖像，渲染出最終的圖像。渲染流程如圖1所示。

圖1 屏幕空間流體粒子渲染流程Fig.1 The diagram of screen space fluid particle rendering

5.3 水體光照效果

當光線照射到透明材質(zhì)的交界面時，一部分光被反射掉，另一部分光發(fā)生折射，這種現(xiàn)象被稱為Fresnel Effect（菲涅爾效應）。在流場渲染時考慮菲涅爾效應可使水體顯得更加真實。通過線光源對水面進行照射，計算漫反射光、鏡面反射光和環(huán)境光，采用菲涅爾方法解決特殊觀察角度下的問題。算法描述如下[21]：

其中：F為菲涅爾函數(shù)，a為折射的流體顏色，b為反射的場景顏色，ks和α為鏡面反射常數(shù)，n為表面法向量，v為視向量，h為視向量與光向量之間的半角向量。在實時圖形渲染時，可將理論物理中的菲涅爾公式簡化如下：

其中：bias為偏移系數(shù)（取值范圍為0～1.0），scale為縮放系數(shù)（取值范圍為0～1.0），power為菲涅爾指數(shù)（通常取1.0），l為入射光向量。

6 DTMB 5415 船模虛擬試驗應用案例

本文以CFD國際標模DTMB 5415為例，比擬拖曳水池中開展的約束模阻力試驗，對在有限尺寸槽形水池中作勻速航行船體周圍的自由面粘性流進行了數(shù)值模擬，展現(xiàn)了船模阻力試驗中自由表面及流場的生成與演化過程。

6.1 數(shù)值模型建立

首先建立船體幾何模型，然后生成船體網(wǎng)格模型，并輸出為DualSPHysics支持的網(wǎng)格文件格式，接下來編寫XML格式的SPH求解輸入文件，設定合適的模型參數(shù)、運動參數(shù)和求解控制參數(shù)，最后對SPH模型進行預處理，生成SPH仿真粒子模型。

6.1.1 建立船體幾何模型

根據(jù)船模阻力試驗相似準則，要求船模與實船保持幾何相似，同時參照CFD數(shù)值試驗結(jié)果，設定縮尺比為1/24.8，船模主尺度參數(shù)如表1所示。

在CAD軟件中建立DTMB 5415船體幾何模型，在模型上方的開口區(qū)域加入曲面，使其成為一個封閉體，導出為IGES中性幾何文件。

6.1.2 生成船體網(wǎng)格模型

將IGES文件導入Patran軟件，修復模型中的縫隙，劃分三角形面網(wǎng)格。通過二次開發(fā)程序，得到VTK（Visualization Toolkit，一種開源可視化工具函數(shù)庫）格式的網(wǎng)格模型。

6.1.3 SPH輸入文件準備

為了減小壁面效應，在計算時間和計算能力允許的情況下，應適當?shù)卦龃笏氐某叽?。在初始位置處，下游邊界距離船艉0.5個船長，側(cè)向0.35個船長，深度方向0.25個船長。槽形水池主尺度參數(shù)如表2所示。

初始條件下，船模處于靜止狀態(tài)。零時刻起，船模在1 s內(nèi)加速到2.096 m/s，然后保持恒定速度和方向航行6 s，最后在1 s內(nèi)減速到0 m/s。按照0.04 s的輸出時間間隔，迭代計算200個時間步內(nèi)的流場變化情況。參照DualSPHysics提供的模板，手工編寫XML格式的求解輸入文件。

表1 DTMB 5415船模主尺度參數(shù)Tab.1 Main scale parameters of the DTMB 5415 ship model

表2 槽形水池主尺度參數(shù)Tab.2 Main scale parameters of the trough shaped water basin

圖2 SPH流體仿真計算域模型Fig.2 Computational domain model of SPH fluid simulation

6.1.4 SPH模型預處理

將水池作為固定邊界，船體作為運動邊界，按照一定的密度均勻填充流體粒子，建立SPH流體仿真的數(shù)值模型。調(diào)用DualSPHysics預處理程序，生成水池邊界粒子、船體邊界粒子和均勻分布的流體粒子如圖3所示。

圖3 SPH流體仿真粒子模型Fig.3 Particle model of SPH fluid simulation

6.2 控制參數(shù)選取

為了兼顧SPH流體仿真的計算效率和精度，需要選取合適的求解控制參數(shù)，本文選擇粒子間距dp值為0.02，時間步長為1e-4。

6.3 試驗過程展現(xiàn)

首先對DualSPHysics輸出的200個時間步的流場數(shù)據(jù)進行壓縮。然后采用SphVizKit程序進行實時交互式立體渲染和表面渲染，虛擬展現(xiàn)DTMB 5415船模在水池中加速、勻速、減速航行過程中自由面及流場的生成與演化情況，如圖4和圖5所示。

參考拖曳水池的內(nèi)部照片，建立水池的三維數(shù)字模型。將試驗對象、試驗環(huán)境置于統(tǒng)一的場景當中，調(diào)整試驗對象和各種試驗件的位置和顯示比例，虛擬重建試驗對象和試驗環(huán)境，允許從岸上、船艏、船艉等不同方位和視角觀察船模及流場，以獲得直觀、全面的流動細節(jié)，如圖6和圖7所示。

圖4 DTMB 5415船體周圍瞬時變化流場Fig.4 Transient varying fluid field around the DTMB 5415 ship hull

圖5 DTMB 5415船體周圍自由面速度云圖Fig.5 Velocity nephogram of free surface around the DTMB 5415 ship hull

圖6 從船艏視角觀察船模及流場Fig.6 Observing ship model and fluid field from the perspective of stem

圖7 從船艉視角觀察船模及流場Fig.7 Observing ship model and fluid field from the perspective of stern

7 結(jié) 語

針對阻力虛擬試驗難以獲得連續(xù)變化流場數(shù)據(jù)的問題，本文提出采用SPH流體仿真方法，對船模阻力試驗中自由表面及流場的生成與演化開展粗粒度的數(shù)值模擬，同時借助三維交互式圖形顯示技術實現(xiàn)試驗過程的虛擬展現(xiàn)，最后以DTMB 5415船模試驗加以例證。

由于知識能力所限，本文的研究工作難免失之偏頗，仍有一些需要改進的地方。例如，為了減小SPH流體粒子的規(guī)模，水池與船模的尺寸比例遠未達到物理試驗中的比例，可采用可變分辨率的SPH仿真算法[23]加以改進；未考慮航行過程中船模的姿態(tài)變化，可采用流體驅(qū)動對象這類邊界條件[17]加以改進。鑒于本文研究工作主要是針對獲取虛擬展現(xiàn)所需的瞬態(tài)流場數(shù)據(jù)而展開，不是以替代CFD方法為目標，因此未與CFD計算結(jié)果進行比較。

[1]趙 峰,吳乘勝,張志榮,等.實現(xiàn)數(shù)值水池的關鍵技術初步分析[J].船舶力學,2015,19(10):1210-1220. Zhao Feng,Wu Chengsheng,Zhang Zhirong,et al.Preliminary analysis of key issues in the development of numerical tank[J].Journal of Ship Mechanics,2015,19(10):1210-1220.

[2]Wackers J,Koren B,Raven H C,et al.Free-surface viscous flow solution methods for ship hydrodynamics[J].Arch Comput Methods Eng,2011,18.

[3]Zhang Zhirong,Zhao Feng,Li Baiqi.Numerical calculation of viscous free-surface flow about ship hull[J].Journal of Ship Mechanics,2002,6(6):10-11.

[4]鄭 興,段文洋.光滑質(zhì)點流體動力學（SPH）及其算法特性[J].船舶力學,2008,12(4):550-551. Zheng Xing,Duan Wenyang.Smoothed particle hydrodynamics(SPH)and its numerical behavior[J].Journal of Ship Mechanics,2008,12(4):550-551.

[5]Alejandro Jacobo Cabrera Crespo.Application of the smoothed particle hydrodynamics model SPHysics to free-surface hydrodynamics[D]:Ph.D.Thesis.University of Manchester,2008.

[6]Altomare C,Rogers B D,Domínguez J M,et al.Numerical modeling of coastal defences using Smoothed Particle Hydrodynamic method[C].4th International Conference on the Application of Physical Modelling to Port and Coastal Protection, 2012.

[7]王 剛.中小尺度海冰動力學的粘彈塑性本構(gòu)模型及SPH數(shù)值模擬[D].大連:大連理工大學,2007. Wang Gang.Viscoclastic-plastic constitutive model for sea ice dynamics under meso-small scale and numerical simulation with SPH method[D].Ph.D.Thesis.Dalian University of`Technology,2007.

[8]Liu G R,Liu M B.Smoothed Particle Hydrodynamics-A meshfree particle method[M].World Scientific Publishing,2003.

[9]Cohen R C Z,Cleary P W,Mason B.Simulations of human swimming using smoothed particle hydrodynamics[C].Proceeding of the Seventh International Conference on CFD in the Minerals and Process Industries,2009.

[10]Cartwright B K,Eng B,Sc M.The study of ship motions in regular waves using a mesh-free numerical method[D].Master Thesis.University of Tasmania,2012.

[11]Le Touzéa D,Marsha A,Ogerb G,et al.SPH simulation of green water and ship flooding scenarios[J].Journal of Hydrodynamics,2010,22(5).

[12]Daniel John Veen.A Smoothed Particle Hydrodynamics study of ship bow slamming in ocean waves[D].Ph.D.Thesis. Curtin University of Technology,2010.

[13]胥 飛,劉 樺.艏外飄型船體入水砰擊的二維SPH模擬[J].水動力學研究與進展,2013,28(5):585-589. Xu Fei,Liu Hua.Numerical simulation of 2-D bow flare slamming using SPH method[J].Journal of Hydrodynamics,2013, 28(5):585-589.

[14]Shao J R,Li H Q,Liu G R,et al.An improved SPH method for modeling liquid sloshing dynamics[J].Computers and Structures,2012:100-101

[15]Micky Kelager.Lagrangian fluid dynamics using Smoothed Particle Hydrodynamics[Z].Department of Computer Science, University of Copenhagen,2006.

[16]Marrone S,Antuono M,Colagrossi A,et al.δ-SPH model for simulating violent impact flows[J].Comput.Methods Appl. Mech.Engrg,2011,200:1526-1542.

[17]Crespo A J C,Domínguez J M,Rogers B D,et al.DualSPHysics:Open-source parallel CFD solver based on Smoothed Particle Hydrodynamics(SPH)[C].Computer Physics Communications,2015.

[18]http://dual.sphysics.org/index.php/validation/(Access Date 2016.4.8)[M/OL].2016.

[19]Green S.Screen space fluid rendering for games(Presentation)[C].Game Developers Conference,2010.

[20]Lorensen W E,Cline H E.Marching cubes:A high resolution 3d surface construction algorithm[J].Computer Graphics, 1987,21(4):163-169.

[21]van der Laan W J,Green S,Sainz M.Screen space fluid rendering with curvature flow[C].SIGGRAPH 2009,Symposium on Interactive 3D Graphics and Games,2009.

[22]Paris S,Kornprobst P,Tumblin J,et al.Bilateral filtering:Theory and applications[C].Foundations and Trends in Computer Graphics and Vision,01/2009.

[23]Vacondio R,Rogers B D,Stansby P K,et al.Variable resolution for SPH:A dynamic particle coalescing and splitting scheme[J].Comput.Methods Appl.Mech.Engrg.,2013,256:

[24]Monaghan J J.Smoothed Particle Hydrodynamics[J].Reports on Progress in Physics,2005,68:1703-1759.

Virtual presentation of the procedure of ship model resistance test using SPH-based fluid simulation

LIU Pei-lin1,2,WANG Mo-wei3,LENG Wen-hao3,JIN Jian-hai3
(1.School of Digital Media,Jiangnan University,Wuxi 214122,China;2.Wuxi Institute of Technology, Wuxi 214121,China;3.China Ship Scientific Research Center,Wuxi 214082,China)

Steady viscous free-surface water flow around ship hull can be acquired by the combination of computing the viscous boundary layer around ship hull by solving the RANS equation and processing the free-surface with VOF method.But it is difficult to obtain continuous varying flow field as analogous to model test in towing tank.Simulating the generation and evolution of free-surface and flow field in ship model resistance test using SPH-based fluid simulation method and making a virtual presentation of the test procedure using three dimensional interactive graphics techniques are presented in this paper.Finally, DTMB 5415 ship model test is taken as an example.

SPH-based fluid simulation;ship model resistance test;virtual presentation

U661.33

A

10.3969/j.issn.1007-7294.2016.08.002

1007-7294（2016）08-0935-12

2016-05-26

江蘇省科技支撐計劃項目—船舶水動力數(shù)值試驗中的虛擬現(xiàn)實技術（BE2012028）

劉培林（1972-），女，博士研究生，教授，E-mail:lpl6666@163.com；王墨偉（1980-），男，高級工程師。