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

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

0 引言

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

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

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

1 控制方程和湍流模型

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

連續(xù)方程為：

動量方程為：

其中：

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

1.1 Standard k-ε 模型

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

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

1.2 RNG k-ε 模型

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

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

1.3 Realizable k-ε 模型

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

其中：

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

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

2.1 計算區(qū)域

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

圖1 測點位置圖Fig.1 Location map of measuring points

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

2.2 數(shù)值方法

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

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

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

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

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

2.3.1 流速剖面

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

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

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

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

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

確定垂向流速為對數(shù)模式u=，通過對實測數(shù)據(jù)的擬合可得到摩阻流速u*，水流底部剪切應(yīng)力，因此流速剖面和底部剪切應(yīng)力可以相互驗證。對一個周期內(nèi)的潮流底部剪切應(yīng)力進行了模擬，結(jié)果見圖5。結(jié)合表1 中潮流底部剪切應(yīng)力的方差值可以看出，在漲潮過程，3種模型的模擬結(jié)果與實測值都很接近，而落潮過程的模擬結(jié)果在相位和數(shù)值大小上與實測值有一定的誤差。整個漲落潮過程，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)力項，而是引入湍動黏度，把湍流應(yīng)力表示成湍動黏度的函數(shù)[7]。渦黏假定建立了雷諾應(yīng)力與平均速度梯度的關(guān)系，即：

式中：δij為“Kronecker delta”符號（當i=j時，δij=1；當i≠j時，δij=0）

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

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

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

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

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

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

3 結(jié)語

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

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

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

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