陳立家 劉錠坤 田延飛 黃 巍

(武漢理工大學(xué)航運學(xué)院1) 武漢 430063) (內(nèi)河航運技術(shù)湖北省重點實驗室2) 武漢 430063) (長江引航中心江陰引航站3) 江陰 214400)

0 引 言

水流模擬的方法包括基于紋理貼圖、幾何形態(tài)、粒子系統(tǒng)、流體動力學(xué)、海浪波譜、柏林噪聲,以及元胞自動機等仿真方法[1-3].對于不同類型水體的仿真，由于對物理真實感和渲染真實感以及實時性、交互性的要求不同，選用的模擬方法各有不同，使用過程中需結(jié)合實際情況進(jìn)行必要的設(shè)置或改進(jìn).

在河流仿真與流體繪制的研究中，冶運濤等[4]采用基于紋理的流場動態(tài)可視化方法來模擬河道水流，綜合了幾何形狀映射和顏色映射方法的優(yōu)點.劉晨龍[5]基于Perlin噪聲和海浪譜的方法繪制水面波浪效果，結(jié)合紋理坐標(biāo)變換和圖像變換的紋理動畫技術(shù)實現(xiàn)河流流動，但該方法生成的河流物理特性較弱.余偉等[6]利用斷面水位線結(jié)合DEM地形構(gòu)建多層水面疊加模型，采取流速驅(qū)動紋理塊自適應(yīng)于河道移動的河流渲染方法，該方法能良好表達(dá)水體的自然特性，但在支流處無法智能實現(xiàn)斷面處理.翟小明等[7]基于真實河道岸界數(shù)據(jù)求解河流速度場，通過瓦片式向量圖實現(xiàn)河流表面的實時繪制，該方法計算簡便，繪制過程基于GLSL著色器并在GPU中執(zhí)行，體現(xiàn)出真實的河流流動方向.

在河流水體的三維模擬研究中，康蘇海[8]基于非靜壓假設(shè)下的三維水動力數(shù)學(xué)模型，結(jié)合立體成像原理實現(xiàn)了三維流場可視化，并將其應(yīng)用于長江復(fù)雜水域的三維水流求解，計算結(jié)果較為精確，但未注重水體渲染；馮小香等[9]基于平面曲線坐標(biāo)系和垂向坐標(biāo)系建立了三維水流數(shù)學(xué)模型，對彎曲河段和船閘引航道口門區(qū)的水流進(jìn)行建模，模擬精度能滿足通航水流條件研究的要求；揚程等[10]用非結(jié)構(gòu)化網(wǎng)格對黃河沙坡頭的連續(xù)彎道處水流運動進(jìn)行研究，垂直方向上進(jìn)行等距分層，采用有限體積法進(jìn)行離散，較好的模擬了天然河道的水流運動.目前，在大規(guī)模場景的河流繪制研究上并未有較好的方法既體現(xiàn)水體的真實物理特性又兼顧實時、細(xì)致的渲染效果.

本文主要研究船舶操縱模擬器中對內(nèi)河特殊航段的三維水流建模與仿真.為增強河流繪制的物理特性，選用三維水動力模型進(jìn)行流場解算，以葛洲壩下游水域為研究對象，水平方向上采用結(jié)構(gòu)化網(wǎng)格，垂向上采用σ坐標(biāo)系進(jìn)行非等距分層，運用多重紋理映射技術(shù)進(jìn)行河流流體渲染，結(jié)合流場的三維特點對水面高度圖及法向貼圖、流動紋理貼圖進(jìn)行制作，并通過折射、反射紋理實現(xiàn)水面光照效應(yīng)，整個過程在OpenGL平臺中利用GLSL著色器進(jìn)行流體渲染，利用GPU并行算法提高繪制質(zhì)量和效率.

1 水動力三維數(shù)學(xué)模型

1.1 控制方程

控制方程基于水平邊界擬合曲線坐標(biāo)系和σ坐標(biāo)系，轉(zhuǎn)換關(guān)系為

(1)

笛卡爾坐標(biāo)系下河流的三維紊流控制方程包括動量方程、連續(xù)方程、溫鹽方程和體積分?jǐn)?shù)方程.

(2)

(3)

(4)

(5)

(6)

(7)

式中：x，y，z分別為笛卡爾直角坐標(biāo)系的東、北和垂直方向的坐標(biāo)；u，v，w分別為笛卡爾坐標(biāo)下x，y，z方向上的流速分量；ρ為水溫T和鹽度S下的水密度；ρ0為水密度平均值，取999.7 kg/m3；g為重力加速度；ζ為自由表面水位；h為水深；H為實際水深，H=h+ζ；pa為大氣壓力；fc=2ωsinφ，為科里奧利系數(shù)；K，vT為水平和垂直方向上的渦流黏性系數(shù)；ωs為單位時間從外界流入的水量；c為水體積分?jǐn)?shù)；w0為浮力引起的垂直方向上的速度；Ts，Ss，cs分別為水源的溫度、鹽度和體積分?jǐn)?shù)；λ為非守恒系數(shù)；Kc，vc為水平和垂直方向上的擴散系數(shù).

為節(jié)省運算時間，本次模型不考慮溫度、鹽度和濃度的影響，采用如下所示的忽略二階水平對流項的標(biāo)準(zhǔn)k-ε模型求得垂直渦流黏性系數(shù)，其中：

(8)

(9)

(10)

(11)

(12)

式中：k為紊動動能；ε為紊動耗散率；Pr為普朗特數(shù)，Pr=0.96；σε=1.3;cμ，c1ε，c2ε，c3ε為經(jīng)驗常數(shù)，cμ=0.09，c1ε=1.44，c2ε=1.92，c3ε在穩(wěn)定條件下，vTN>0，c3ε=0.001，在非穩(wěn)定條件下，vTN<0，c3ε=2.

模型采用有限差分方法進(jìn)行數(shù)值離散，流場在水平方向上采用矩形網(wǎng)格進(jìn)行劃分，垂直方向采用σ坐標(biāo)系進(jìn)行非等比例分層.

1.2 邊界條件

1) 自由表面與底部邊界條件[11].

在自由表面z=ζ(x,y,t)處：

(13)

(14)

式中：τs(x),τs(y)為表面風(fēng)應(yīng)力在x和y方向上的分量；ρa為空氣密度，取1.225 kg/m3；W=(wx,wy)為水面上10 m高度處的風(fēng)速；Cd為風(fēng)應(yīng)力拖曳系數(shù)，Cd=(0.63+0.066|W|)10-3.

在底部z=-h(x,y)處：

(15)

(16)

fb=[k/ln(zb/z0)]2

(17)

式中：τb(x),τb(y)為底部應(yīng)力在x和y方向上的分量；fb為底部拖曳系數(shù)；κ為卡曼常數(shù),κ=0.4；zb為無量綱層厚度；zo為底部粗糙度參數(shù).

2) 固邊界采用無滑移邊界條件，無質(zhì)量和熱量交換，固邊界處法向和切向速度為零.

3) 開邊界條件包括水位、水量邊界條件和流速邊界條件.采用實測水位、水量和流速資料作為開邊界條件.

4) 出水淺灘采用動邊界處理，即對計算區(qū)域內(nèi)灘地干濕過程采用水位判別法處理，設(shè)置判別水深為干水深δ1和濕水深δ2.當(dāng)某單元水深小于δ1時，該處流速為零，灘地干出，不參與計算；單元水深在δ1和δ2之間時，只進(jìn)行連續(xù)方程的計算；單元水深大于δ2時參與所有控制方程的計算.

2 水面渲染方案

本文基于GLSL著色器的多重紋理映射技術(shù)進(jìn)行水面渲染，最終形成具有真實感的河流流體，渲染過程在GPU中執(zhí)行，使用可編程渲染管線處理數(shù)據(jù)，以減小CPU運算量，提高計算效率，增強渲染效果.渲染流程見圖1，首先將水動力模型計算出的水位信息結(jié)合三維水波算法[12]制作水流高度圖，并利用凹凸紋理映射技術(shù)將水流高度圖轉(zhuǎn)化為法向貼圖，經(jīng)著色器渲染生成水波波紋；添加流動紋理貼圖控制波紋的偏移，從而控制河流流向和大小，實現(xiàn)河流流動效果；最后結(jié)合菲涅爾效應(yīng)生成的混合紋理實現(xiàn)水面光照效應(yīng).

圖1 基于多重紋理映射技術(shù)的流動水面渲染

2.1 法向貼圖

使用凹凸映射技術(shù)將高度圖中儲存的高度值進(jìn)行轉(zhuǎn)化，生成法向貼圖，將法向貼圖任一點周圍的四個法向量混合疊加形成混合法向量，再經(jīng)過著色器渲染即可形成水面波紋.

首先對水動力方程計算出的水位值插值得出水面上任一點高度值，生成初始水面高度；再根據(jù)三維水波算法，生成具備起伏變換效果的水波，并將水波高度與初始水面高度疊加，最終得出可實時變換的動態(tài)水面高度圖.在頂點著色器中通過紋理坐標(biāo)讀取高度信息，并將其置換到水面模型中，可獲得凹凸不平的波浪效果.

法向貼圖中各點處光線的法向量值可以將相鄰點高度值作差后叉乘得到，空間坐標(biāo)系中法線由X，Y，Z軸上三個分量組成，法線信息儲存在貼圖的R，G，B三個通道中.

將高度圖轉(zhuǎn)化后的法線貼圖進(jìn)行法向量疊加，為

1) 已知點pa，將其坐標(biāo)分別在x、y和沿45°角方向上進(jìn)行位移，得到pb，pc，pd點的坐標(biāo).Normalpa、Normalpb、Normalpc、Normalpd分別為坐標(biāo)點對應(yīng)的法向量.

2) 將pa點與pb點法向量在x方向上進(jìn)行混合生成新的法向量Normalpapb，將pc點與pd點法向量在x方向上混合生成法向量Normalpcpd，再將Normalpapb與Normalpcpd在y方向上混合得到混合法向量Normal，經(jīng)過著色器渲染即可形成水面波紋.

3) 引入二維偏移變量，該偏移值大小為該任一點處x,y方向上的速度-時間函數(shù)，用該值去影響pa,pb,pc,pd坐標(biāo)值，可得到隨時間變化的波紋.

因高度圖信息基于流場的水位信息，偏移變量與水面實際速度相關(guān)，因此水面波紋的效果與實際波動狀態(tài)相符.

2.2 流動紋理

要實現(xiàn)水面的流動效果，還需要采用動態(tài)紋理技術(shù)，通過流動紋理貼圖控制紋理的移動體現(xiàn)出實際水流的方向和流速.流動紋理的R,G,B,A各通道用于保存流體信息，R通道保存x方向的速度信息，G通道保存y方向的速度信息，通過流動紋理的R值和G值對法向紋理坐標(biāo)的方向進(jìn)行旋轉(zhuǎn)，改變紋理映射的方向，從而改變波紋的流動方向；B通道保存河流的速度大小，通過B的值對紋理坐標(biāo)進(jìn)行縮放，可改變法向貼圖在著色器中的映射數(shù)量，從而控制波紋大??；A通道為各點處混合法向量參數(shù)，通過A的值控制河道外混合法向量為零向量，對河道外部分不進(jìn)行渲染，以控制河道邊界.

2.3 菲涅爾效應(yīng)混合紋理

在現(xiàn)實中，對水面進(jìn)行觀察時，觀察者水面角度越小，反射效果越明顯；角度越大，折射效果越明顯，這種現(xiàn)象稱為菲涅爾效應(yīng).為提高水面模擬的真實感，考慮到菲涅爾效應(yīng)的影響，分別獲取水面反射紋理與折射紋理，并將兩者進(jìn)行混合，使用頂點著色器和片元著色器對水面進(jìn)行渲染，最后得到與實際水面光照效應(yīng)接近的渲染效果.采用以下方式進(jìn)行簡單計算.

FresnelRatio=F+(1-F)×

pow((1-dot(eveDir,normal)),5)

(18)

finalColor=waterColor+reflectionColor×

FresnelRatio+refractionColor×

(1-FresnelRatio)

(19)

式中：F=0.200 37，為空氣到水面的反射系數(shù)；eyeDir為視線方向；normal為光線與水面交點處的法線方向；finalColor為菲涅爾效應(yīng)渲染的水面顏色；waterColor為初始水的頂點顏色；reflectionColor為水面的反射顏色；reflractionColor為水面的折射顏色.

3 案例研究

3.1 計算范圍

模型計算區(qū)域為葛洲壩下游水域，包括主航道與三江航道，及其交匯水域，該處河段水流狀況復(fù)雜，三江下引航道口門與大江來流夾角過大，存在的較大橫流與回流使得船舶通航難度加大.因此，有必要對該類水文條件復(fù)雜的通航水域進(jìn)行三維流場模擬，并開展水域內(nèi)的船舶操縱模擬研究，從而提高內(nèi)河航行安全.研究所用資料為長江勘測規(guī)劃設(shè)計研究院在葛洲壩河段的河道水文數(shù)據(jù)，計算域范圍為30°40′29.9N″～30°43′31.9″N，111°15′26.6″E～111°18′34.2″E，水域面積約6 km2，零水位處水深分布見圖2.

圖2 模型計算區(qū)域水深

3.2 計算條件

根據(jù)《三峽(正常運行期)—葛洲壩水利樞紐梯級調(diào)度規(guī)程》，葛洲壩大江1號船閘最大通航流量為35 000 m3/s，葛洲壩上游運行水位為66.0 m.模擬方案控制下游大江航道下泄流量為35 000 m3/s，同時二號船閘充水，控制三江航道流量為1 600 m3/s.試驗采用Windows7操作系統(tǒng)，以O(shè)SG為視景渲染引擎，使用GLSL著色器語言進(jìn)行河流流體渲染，OpenGL為4.3版.

計算網(wǎng)格平面上采用矩形網(wǎng)格，分別采用方案一30 m網(wǎng)格與方案二10 m網(wǎng)格尺度計算，見圖3，其中CS1、CS2分別為兩種網(wǎng)格在同一位置處的河流橫截面，緯度為30°42′15.3″N.

圖3 交匯河段計算網(wǎng)格

垂向網(wǎng)格劃分為6層，其中表層與底層厚度較小，垂向網(wǎng)格劃分方法見圖4.

圖4 垂向網(wǎng)格劃分方法

3.3 參數(shù)設(shè)置及調(diào)整

參照河段地貌與水流阻力特征，河段模型糙率系數(shù)取值為全區(qū)域一致0.025；垂直渦流黏度系數(shù)最小值取分子黏度1.2×10-6；水平渦流黏度系數(shù)為參數(shù)化亞尺度湍流渦旋中能量耗散的經(jīng)驗系數(shù)，為有效防止數(shù)值震蕩，該值大小為方案一1.1，方案二1.0.

考慮到對流項對水動力計算過程中的影響，時間步長取值受到條件限制為

(20)

式中：Δt為時間步長；Δx為網(wǎng)格步數(shù)；U為水流流速.根據(jù)該河段實際水流特征，時間步長取值為方案一1.5 s，方案二0.5 s.

3.4 結(jié)果分析

圖5～6為主河段和三江河段交匯處的表層、中層、底層流場及回流區(qū)表層流場分布狀況.垂向流速從表層到底層逐漸減小，河口交匯處存在斜流與回流.三江河段與主河段流速、流量相差較大，結(jié)合圖7截面處流速的橫向分布，可明顯看出長江主流相對于三江航道中心線而言的斜流流態(tài).沿三江流向的水流在交匯區(qū)域靠左岸處流速降低，并在河口處形成回流，回流以逆時針方向旋轉(zhuǎn).受網(wǎng)格精度影響，30 m網(wǎng)格中無明顯回流流態(tài)，10 m網(wǎng)格中回流明顯，匯流口流線清晰，匯流角明顯.從實驗結(jié)果來看，三維流場模擬精度與實際吻合度較高.

圖5 河流交匯處流場分布(30 m網(wǎng)格)

圖6 河流交匯處流場分布(10 m網(wǎng)格)

圖7 CS2截面處橫向流速分布

圖8為交匯河口水面流動效果圖，方案二湍流特性較方案一好，該處部分區(qū)域流速較大，水面波動狀況復(fù)雜，方案一受網(wǎng)格精度影響，渲染無回流流態(tài).由于回流區(qū)流速相較于主流流速太小，因而回流流動并不明顯，回流區(qū)渲染效果體現(xiàn)為回流流速慢，流態(tài)相對平靜.

圖8 河流交匯處水面流動效果

本文所用方法的河流流體繪制過程在GPU中執(zhí)行，10 m網(wǎng)格和30 m網(wǎng)格的河流流體渲染幀率平均分別在42幀和45幀左右，可實現(xiàn)場景的流暢運行.總體上水流渲染結(jié)合流場的物理特性，水體流動變換自然，體現(xiàn)出河流的真實流向，河口匯流處波動現(xiàn)象與實際符合度較好.

4 結(jié) 束 語

本文以葛洲壩下游水域為例，建立了內(nèi)河河流三維數(shù)值流場，對壩后彎曲河段及主、支流交匯河口的復(fù)雜流態(tài)進(jìn)行模擬.使用GLSL著色器進(jìn)行流體渲染，以計算出的水位信息為基礎(chǔ)生成水流高度圖，以此制作了符合水面實際光影效果的法向貼圖，并根據(jù)計算出的速度場信息制作了流動紋理貼圖，最后結(jié)合菲涅爾效應(yīng)混合紋理實現(xiàn)了河流流體繪制.

實驗表明該方法計算出的水流仿真精度較高，能還原出斜流、回流等復(fù)雜流態(tài)，結(jié)構(gòu)化網(wǎng)格的數(shù)據(jù)結(jié)構(gòu)簡單，模型整體的計算效率高.流體渲染通過GPU并行計算完成，可編程渲染管線能靈活的進(jìn)行紋理映射，渲染過程結(jié)合流場的物理特性，生成的河面起伏波動效果較好，通過河流的速度場信息驅(qū)動紋理移動，能體現(xiàn)出實際的流速流向.總體來說，基于GLSL的多重紋理映射技術(shù)能較為真實的繪制出河流水體，生成與實際相符的水面光照效應(yīng)和流動效果.

本文方法建立的三維水流模型有效提高了內(nèi)河航道水流仿真精度，渲染方案能夠?qū)崿F(xiàn)水流模擬的真實感.對特定航段水流仿真的成果可應(yīng)用于航道水流條件分析、內(nèi)河船舶操縱模擬器系統(tǒng)建設(shè)、船舶安全操縱決策等方面.