郝嘉凌，劉鎮(zhèn)，胡春也

（河海大學(xué)港口海岸與近海工程學(xué)院，江蘇 南京 210098）

0 引言

近岸水流在大部分時(shí)間段為單向水流，潮流在水流運(yùn)動(dòng)中起著主導(dǎo)作用，潮流的邊界層如明渠單向水流的邊界層一樣可以得到充分的發(fā)展，使邊界層所在水深范圍很大，甚至可以擴(kuò)展到整個(gè)水深[1]。對于該區(qū)域流場，雖然在底床附近仍是低雷諾數(shù)的湍流，但是隨著垂向尺度增大，逐漸發(fā)展的渦旋使湍流變成高雷諾數(shù)的狀態(tài)。因此，簡單地將計(jì)算管流或者水深很小的明渠流的傳統(tǒng)湍流模型應(yīng)用到該區(qū)域，模擬得到的結(jié)果是不精確甚至是不適用的。

在第二屆國際阻力預(yù)測會(huì)議中發(fā)現(xiàn)，湍流模型的選取是解決湍流問題最重要的影響因素，對計(jì)算結(jié)果的影響占到15%[2]。雷諾平均模擬方法（RANS），以雷諾平均運(yùn)動(dòng)方程與脈動(dòng)運(yùn)動(dòng)方程為基礎(chǔ)，通過引入相關(guān)假設(shè)建立起封閉模型，避免了對湍流運(yùn)動(dòng)的直接模擬，降低了大量的計(jì)算成本，同時(shí)又能保證較高的精確度，是目前最流行的湍流模擬方法[3]。雷諾平均模擬方法數(shù)值模型眾多，其中k-ε 模型和Mellor-Yamada 模型在海洋水動(dòng)力模型中應(yīng)用最為廣泛，但Mellor-Yamada模型會(huì)低估充分發(fā)展湍流的垂向混合，導(dǎo)致水流和湍流耗散之間的相位滯后偏小[4]。針對充分發(fā)展的高雷諾數(shù)湍流，以Standard k-ε 模型、RNG k-ε模型、Realizable k-ε 模型應(yīng)用最為廣泛。另外，對于床層附近的低雷諾數(shù)區(qū)域湍流，不同壁面函數(shù)的選取對模擬結(jié)果有很大影響，針對這3 種k-ε模型，標(biāo)準(zhǔn)壁面函數(shù)能表現(xiàn)出對網(wǎng)格和流動(dòng)良好的適應(yīng)性[5]。

潮流的一些重要特征量能反映流場的相關(guān)特性，除了能較為直接地描述水流特性的流速剖面和床層剪切應(yīng)力外，湍動(dòng)黏度的數(shù)值和湍動(dòng)能的衰減可以反映水流的流速梯度及擴(kuò)散程度[6]，雷諾應(yīng)力是水流質(zhì)點(diǎn)脈動(dòng)產(chǎn)生的且能與平均速度梯度建立關(guān)系[7]。本文采用3 種k-ε 模型，用標(biāo)準(zhǔn)壁面函數(shù)法解決在底床附近低雷諾數(shù)湍流的問題，對潮流進(jìn)行了一維垂向數(shù)值模擬，通過對模擬結(jié)果的分析，得出計(jì)算不同的湍流特征值時(shí)適用的模型及其優(yōu)缺點(diǎn)，為進(jìn)一步研究沿岸泥沙的遷移擴(kuò)散、侵蝕沉積過程、懸浮濃度分布等提供參考[8]。

1 控制方程和湍流模型

近岸水流的一維垂向的基本方程為連續(xù)方程和動(dòng)量方程[9]。

連續(xù)方程為：

動(dòng)量方程為：

其中：

式中：u為時(shí)均水平速度；t為時(shí)間；ρ為流體密度；w為時(shí)均垂向流速；z為垂向坐標(biāo)，在床底z=0，在水表面z=h+A，h為平均水深，A為潮流振幅；P為作用在流體微元體上的壓力；μt為湍流運(yùn)動(dòng)黏度；τij為雷諾應(yīng)力；分別為x（平行潮流流速方向）、z方向上的速度時(shí)均值。

1.1 Standard k-ε 模型

Launder 和Spaldingt 提出了基于湍動(dòng)能k和湍動(dòng)耗散率ε的Standard k-ε 模型，對于不可壓流體，與之相應(yīng)的輸運(yùn)方程分別為：

式中：湍動(dòng)黏度μt=ρCuk2/ε，Cu為黏度系數(shù)，取Cu=0.09；μ為流體動(dòng)力黏度；P為湍流動(dòng)能的產(chǎn)生項(xiàng)，P=μt（?u/?z）；σk、σε分別為湍動(dòng)能k、湍動(dòng)耗散ε對應(yīng)的普朗特?cái)?shù)，取σk=1.0，σε=1.3；C1、C2為模型常數(shù)，取C1=1.44，C2=1.92。

1.2 RNG k-ε 模型

RNG k-ε 模型的湍動(dòng)能k和湍動(dòng)耗散率ε分別對應(yīng)的輸運(yùn)方程如下：

式中：μeff為修正湍動(dòng)黏度；Cu為黏度系數(shù)，取Cu=0.085；、η為模型系數(shù)；Eij為時(shí)均應(yīng)變率；模型常數(shù)η0=4.377，β=0.012，αk=αε=1.39，C1ε=1.42，C2ε=1.68。

1.3 Realizable k-ε 模型

Realizable k-ε 模型的湍動(dòng)能k和湍動(dòng)耗散率ε分別對應(yīng)的輸運(yùn)方程如下：

其中：

式中：模型常數(shù)σk=1.0，σε=1.2，=1.9；、A0、As為模型系數(shù)，A0=4.0；φ、W、U*為通用變量；為時(shí)均轉(zhuǎn)動(dòng)速率張量，對于無旋流場為0；ωk為角速度。

2 3 種k-ε 模型在近岸潮流計(jì)算中的應(yīng)用

2.1 計(jì)算區(qū)域

英國Menai 海峽連通Caernarfon 灣和北愛爾蘭海，全長約20 km，平均寬度約為800 m。海峽受以半日潮M2 為主的潮流影響，潮流振幅為2.25 m，大潮時(shí)，海峽西南入口的水流速度可達(dá)2.5 m/s[10]。當(dāng)水流流向?yàn)闁|北方向時(shí)，海平面處于上升狀態(tài)，該過程為漲潮（定義為x正方向）。當(dāng)水流流向?yàn)槲髂戏较驎r(shí)為落潮（定義為x負(fù)方向），水道見圖1。

圖1 測點(diǎn)位置圖Fig.1 Location map of measuring points

采用1 200 kHz 的ADCP 布置在海峽截面中部的底床上，對該水道的垂向流速進(jìn)行了測量。測點(diǎn)平均水深12.5 m 左右，床層以上1 m 以及水面下1 m 區(qū)域是沒有數(shù)據(jù)的，垂向數(shù)據(jù)點(diǎn)間隔距離為0.5 m。本文模擬時(shí)長遠(yuǎn)大于實(shí)測時(shí)長，對一個(gè)潮周期的模擬值和實(shí)測數(shù)據(jù)進(jìn)行比較分析。

2.2 數(shù)值方法

Menai 海峽受強(qiáng)潮影響，含沙量較小，淡水輸入量小，且湍流發(fā)展充分，無分層現(xiàn)象[10-11]。因此針對計(jì)算區(qū)域的水流，不考慮鹽度分層及懸沙分層的影響，對水流進(jìn)行一維垂向數(shù)值模擬。

在以上3 個(gè)模型中，對k 方程和ε 方程求解時(shí)都采用理論上無條件穩(wěn)定的Crank-Nicolson 時(shí)間積分方案。由于潮流水深隨時(shí)間在變化，網(wǎng)格數(shù)設(shè)置為nz-1不變，網(wǎng)格尺寸Δz為變化值，為了避免造成穩(wěn)定性問題，選定足夠小的時(shí)間步長Δt=0.5 s。變量u，ρ，μt，k和ε設(shè)置在網(wǎng)格中心，只有w設(shè)置在網(wǎng)格底部。網(wǎng)格設(shè)置如圖2 所示。

圖2 垂向一維水動(dòng)力模型中設(shè)置的網(wǎng)格和尺寸Fig.2 Mesh and dimensions set in vertical one-dimensional hydrodynamic model

模型采用無滑移不可入固底邊界條件，即z=0 處u=0，在無風(fēng)情況下的自由表面z=h+η處，流速梯度du/dz=0，us=ks=εs=0。為了計(jì)算精確，將垂向水流分為黏性層、過渡層和對數(shù)層，由于黏性層近似層流，k-ε 模型對于該區(qū)域是不適用的。對黏性層的k和ε值進(jìn)行線性處理，在黏性層與過渡層交界，采用標(biāo)準(zhǔn)壁面函數(shù)令該處的湍動(dòng)能，湍動(dòng)耗散率，其中：u*為摩阻流速；κ為卡門常數(shù)取為0.41；Zb為黏性底層厚度，調(diào)試后可取Zb=0.1 m[8，12]。

2.3 計(jì)算結(jié)果與分析

2.3.1 流速剖面

模型對恒定流情況進(jìn)行了模擬，選取Gonzalez等人用1 200 kHz ADCP 在Chicago Sanitary and Ship Canal 測量的穩(wěn)定明渠流數(shù)據(jù)進(jìn)行驗(yàn)證[13]。比較3 種湍流模型在恒定流和潮流情況下計(jì)算流速剖面的差異，見圖3、圖4。

圖3 恒定流條件下3 種模型計(jì)算的流速剖面Fig.3 Velocity profiles calculated by three models under steady flow conditions

圖4 潮流條件下3 種模型計(jì)算的流速剖面Fig.4 Velocity profiles calculated by three models under tidal current conditions

從圖3 模擬結(jié)果與實(shí)測數(shù)據(jù)可以看出，對于恒定流而言，3 種模型計(jì)算的流速剖面結(jié)果很相似，與實(shí)測數(shù)據(jù)都很接近。然而近岸水流以潮流場為主，且漲、落潮峰值和潮流時(shí)域特性對港工建筑、工程整治措施影響顯著[14-15]。潮流條件下3種模型計(jì)算的流速剖面結(jié)果有較為明顯的差異，結(jié)合表1 中流速方差值可知，RNG k-ε 模型和Realizable k-ε 模型相比于Standard k-ε 模型對潮流流速剖面模擬更為精確，尤其是在落急時(shí)刻Standard k-ε 模型誤差較大。對于恒定流和潮流模擬，模擬結(jié)果與實(shí)測數(shù)據(jù)的誤差隨著水深增加而變大，且模擬值較實(shí)測值偏大，這可能是由于模型中未考慮風(fēng)應(yīng)力條件導(dǎo)致的。

表1 潮流流速剖面與剪切應(yīng)力的模擬結(jié)果與實(shí)測值的方差Table 1 Square deviation of simulated and measured values of tidal current velocity profile and shear stress

確定垂向流速為對數(shù)模式u=，通過對實(shí)測數(shù)據(jù)的擬合可得到摩阻流速u*，水流底部剪切應(yīng)力，因此流速剖面和底部剪切應(yīng)力可以相互驗(yàn)證。對一個(gè)周期內(nèi)的潮流底部剪切應(yīng)力進(jìn)行了模擬，結(jié)果見圖5。結(jié)合表1 中潮流底部剪切應(yīng)力的方差值可以看出，在漲潮過程，3種模型的模擬結(jié)果與實(shí)測值都很接近，而落潮過程的模擬結(jié)果在相位和數(shù)值大小上與實(shí)測值有一定的誤差。整個(gè)漲落潮過程，Realizable k-ε 模型模擬最為精確，RNG k-ε 模型次之。

圖5 潮流底部剪切應(yīng)力Fig.5 Shear stress at the bottom of the tidal current

2.3.2 雷諾應(yīng)力

在渦流模型方法中，不直接處理雷諾應(yīng)力項(xiàng)，而是引入湍動(dòng)黏度，把湍流應(yīng)力表示成湍動(dòng)黏度的函數(shù)[7]。渦黏假定建立了雷諾應(yīng)力與平均速度梯度的關(guān)系，即：

式中：δij為“Kronecker delta”符號（當(dāng)i=j時(shí)，δij=1；當(dāng)i≠j時(shí)，δij=0）

雷諾應(yīng)力模擬結(jié)果見圖6。從圖6 可以看出，RNG k-ε 模型和Realizable k-ε 模型能對漲潮情況下的潮流雷諾應(yīng)力精準(zhǔn)模擬，Standard k-ε 模型模擬結(jié)果誤差較大。對于雷諾應(yīng)力分布比較分散的落潮情況，3 種模型模擬的結(jié)果與實(shí)測值都有較大的誤差。在潮流底部，模擬值急劇減小，是由于模型對于垂向上進(jìn)行了分層，從對數(shù)層計(jì)算到黏性底層時(shí)，水流接近層流，流速梯度驟減，這與實(shí)際的湍流情況是相符的。

圖6 雷諾應(yīng)力模擬結(jié)果Fig.6 Reynolds stress simulation results

2.3.3 湍動(dòng)能與湍動(dòng)黏度的無量綱數(shù)

湍動(dòng)能與湍動(dòng)黏度的無量綱數(shù)模擬結(jié)果見圖7、圖8。湍動(dòng)能的無量綱數(shù)、湍動(dòng)黏度的無量綱數(shù)N′=μt/κu*h，其中，湍動(dòng)能的產(chǎn)生項(xiàng)P=。從圖7 可以看出，對于湍動(dòng)能的無量綱數(shù)的模擬，在漲急時(shí)刻，RNG k-ε 模型的模擬結(jié)果要明顯優(yōu)于其他兩種；在落急時(shí)刻，3種模型模擬的結(jié)果與實(shí)測數(shù)據(jù)趨勢一致，RNG k-ε模型模擬結(jié)果較優(yōu)，但仍有一定誤差。從圖8 可知，對于湍動(dòng)黏度的無量綱數(shù)的模擬，在漲急時(shí)刻，Realizable k-ε 模型模擬結(jié)果與實(shí)測數(shù)據(jù)最接近，其他2 種模型的計(jì)算值相對實(shí)測值較??；在落急時(shí)刻，3 種模型模擬結(jié)果與實(shí)測數(shù)據(jù)趨勢一致，湍動(dòng)黏度的分布規(guī)律符合一維垂向湍動(dòng)黏度的標(biāo)準(zhǔn)拋物線形[12]，但存在一定誤差，表現(xiàn)為在底部和表面位置模擬值偏大，而在中部位置模擬值偏小。

圖7 湍動(dòng)能的無量綱數(shù)模擬結(jié)果Fig.7 Dimensionless number simulation results of turbulent kinetic energy

圖8 湍動(dòng)黏度的無量綱數(shù)模擬結(jié)果Fig.8 Dimensionless number simulation results of turbulent viscosity

3 結(jié)語

根據(jù)以上模擬結(jié)果與分析，可得如下結(jié)論：

1）3 種模型在漲急時(shí)刻的特征值模擬結(jié)果比在落急時(shí)刻的精確。

2）對于恒定流流速剖面的模擬，3 種模型的模擬值幾乎沒有差異且均適用；對于潮流流速剖面模擬，RNG k-ε 模型和Realizable k-ε 模型模擬都比較精確。

3）在漲潮過程中，對于底部剪切應(yīng)力、雷諾應(yīng)力、湍動(dòng)黏度的無量綱數(shù)，選用Realizable k-ε模型計(jì)算最精確，RNG k-ε 模型模擬湍動(dòng)能的無量綱數(shù)最精確；對于落潮過程，3 種模型能正確地模擬特征值的分布規(guī)律，但有一定誤差，原因有兩方面：①實(shí)測切應(yīng)力數(shù)據(jù)垂向分布發(fā)生局部突變時(shí)，模型沒有考慮此效應(yīng)，落潮切應(yīng)力垂向分布發(fā)生突變主要受漲落潮不對稱和地形變化影響；②由于該區(qū)域海峽呈收縮狀，漲潮水流收束，落潮水流分散，模型只計(jì)算了x方向特征值而忽略橫向水流的影響。