基于WMLES方法的復式斷面明渠三維紊流數(shù)值模擬

曾 誠，丁少偉，周 婕，王玲玲，陳 辰

(1.河海大學水利水電學院，江蘇 南京 210098； 2.河海大學力學與材料學院，江蘇 南京 210098)

在平原地區(qū)，天然河道和人工渠道斷面大多為復式結構，包括主槽和灘地。平時水流主要集中在主槽中流動，汛期水流漫過灘地，主槽和灘地共同過流，形成復式斷面流動。由于主槽和灘地上水流之間較大的流速差異，導致在灘槽交界處出現(xiàn)較大的動量交換，使得內(nèi)部水流結構體現(xiàn)出三維流動特征。在垂直于主流的橫斷面上，出現(xiàn)由紊流驅(qū)動的二次流，引起橫向的物質(zhì)輸運、岸坡沖刷等問題。因此，研究復式斷面明渠紊流的三維流動特征，對于河道行洪斷面設計和加固具有重要意義。

許多學者對復式斷面明渠灘槽交界處二次流和動量交換機理進行了試驗和理論研究。Rajaratnam等[1]、Myers[2]、Wormleaton等[3]、Tominaga等[4]、Shiono等[5]、Sanjou等[6]先后對復式斷面明渠紊流進行了試驗測量，對主槽和灘地間的動量交換機理進行了研究。他們的研究結果表明灘槽交界面上的切應力比固體邊界上的平均剪切應力大許多倍，證明復式斷面明渠流動中橫向動量交換強于垂向動量交換，在流動中占據(jù)主導地位。

在復式斷面明渠流的數(shù)值模擬研究方面，大量成果均基于雷諾時均方程的渦黏性紊流模型。Sofialidis等[7]發(fā)現(xiàn)當相對水深比和雷諾數(shù)較小時，非線性k-ω模型比非線性k-ε模型有更好的表現(xiàn)。林斌良等[8]采用非線性k-ε模型和雷諾應力模型(Reynolds stress model，RSM)對復式斷面明渠紊流進行了模擬研究，結果表明兩種模型均能較好地模擬復式斷面明渠三維紊流結構，RSM得出的主流速分布更接近試驗結果，兩種模型得出的雷諾應力分布趨勢與試驗結果一致，但在數(shù)值上偏大較多。梁愛國等[9-11]采用RSM對復式斷面彎道明渠三維流場結構進行了模擬研究，揭示了水深比和彎轉(zhuǎn)角度對流場中漩渦結構的影響。Cater等[12-13]運用大渦模擬(LES)方法對非對稱矩形復式斷面明渠進行了模擬，得到了詳細的流場信息，且精度較高。許棟等[14]采用LES方法對復式斷面明渠水流進行了數(shù)值模擬研究，揭示了水深比和雷諾數(shù)對其紊流結構的影響。

前人的研究結果[12-14]表明，LES方法對于復式斷面明渠流動具有較好的模擬精度，但是在邊界層內(nèi)需要足夠細密的網(wǎng)格來保證最小尺度的漩渦結構能夠被捕捉到。同時，對于復式斷面明渠流問題而言，由于垂直方向上存在兩個邊界層區(qū)域(灘地底面和主槽底面)，導致利用LES方法對該類問題進行模擬研究的計算成本很高。Shur等[15]提出了一種壁面建模的大渦模擬(wall-modeled LES，WMLES)方法，克服了傳統(tǒng)LES方法在模擬高雷諾數(shù)下近壁面紊流時要求非常精細網(wǎng)格尺度的問題。與傳統(tǒng)LES方法相比，WMLES無須加密邊界層中的網(wǎng)格，可節(jié)約計算成本。目前，這種方法已經(jīng)越來越多地應用于數(shù)值模擬中[16-17]。本文采用WMLES方法對復式斷面明渠流動進行模擬研究，并將計算得到的主流速、二次流、床面切應力、紊動強度、紊動能及雷諾應力與已有的試驗測量值[4]和LES模擬結果[12-13]進行對比分析，以驗證WMLES方法對于此類流動問題的適用性和可靠性。

1 數(shù)學模型

1.1 控制方程

連續(xù)方程和動量方程的過濾形式如下：

(1)

(2)

由濾波造成的亞格子應力是未知的，須通過建模進行求解。本文采用Boussinesq假設對亞格子應力進行計算：

(3)

νt=min[(κdw)2,(CsmagΔ)2]·

|S-Ω|{1-exp[-(y+/25)3]}

(4)

式中：dw為計算點與壁面間距；S為應變速率；Ω為渦量；參數(shù)κ=0.41，Csmag=0.2；Δ為濾波尺寸；y+為垂直壁面的無量綱距離。傳統(tǒng)LES方法中，濾波尺寸通常通過下式進行計算：

Δ=(ΔxΔyΔz)1/3

(5)

式中：Δx、Δy、Δz分別為沿x、y、z軸方向的網(wǎng)格尺寸。這種計算方法適用于各項同性或者近似各向同性網(wǎng)格，導致計算網(wǎng)格在單元數(shù)量上過于密集，進而導致計算成本提高。與傳統(tǒng)的LES方法相比，WMLES方法在計算濾波尺寸時引入了壁面阻尼函數(shù)項(如式(6)所示)，使得該方法適用于各項異性網(wǎng)格，從而降低網(wǎng)格數(shù)量，提高計算效率。

Δ=min[max(Cwdw,Cwhmax,hwn),hmax]

(6)

式中：hwn為沿壁面法向的網(wǎng)格尺度；hmax為壁面網(wǎng)格的最大尺度；Cw為常數(shù)，取為0.15。

1.2 數(shù)值方法與邊界條件

為保證流場充分發(fā)展，選取SSTk-ω計算的穩(wěn)態(tài)場作為WMLES計算的初始流場。離散格式選取具有三階精度的QUICK格式，壓力場和速度場的解耦計算采用PISO算法。

在入口和出口邊界上采用周期性邊界條件，自由水面邊界采用剛蓋假定，固體壁面為無滑移邊界。

2 計算域及網(wǎng)格劃分

選取Tominaga等[4]的物理模型試驗測量數(shù)據(jù)進行模型驗證和分析。試驗在長為12.5 m，主槽寬度b=0.2 m，河道總寬度B=0.4 m的復式斷面水槽中進行，采用激光多普勒測速儀測量流速。選取水深比hr=0.5的工況進行模擬(hr=h/H，其中h為灘地水深，H為主槽水深)，主槽水深和灘地水深分別為0.08 m和0.04 m，平均流速Um為0.349 m/s，雷諾數(shù)為54 500，摩阻流速為0.016 4 m/s。圖1為數(shù)值模擬計算區(qū)域示意圖，x、y、z分別為流向、展向和垂向坐標。數(shù)學模型的計算域大小為10H×5H×H(x×y×z，下同)。計算中為了便于結果處理與分析，采用主槽水深H對坐標系無量綱化。水槽沿x向坡度S0=6.397×10-4。

圖1 計算區(qū)域和尺寸

采用六面體網(wǎng)格對計算區(qū)域進行剖分，在x、y、z方向均采用均勻結構化正交網(wǎng)格，未對邊界層進行加密。網(wǎng)格劃分與數(shù)值模擬的精度及計算結果的收斂性密切相關，因此本文選取3組不同精度的網(wǎng)格進行網(wǎng)格無關性驗證，工況C1、C2、C3網(wǎng)格數(shù)量分別為128×100×64、160×200×80、 200×400×100。3組不同網(wǎng)格精度工況的流向時均流速U沿橫斷面分布如圖2所示，圖中以平均速度Um對U進行無量綱化。由圖2可知，工況C1由于網(wǎng)格精度不足，與另外兩組工況計算結果差距較大。工況C2和C3流向時均流速U的分布基本一致，誤差在5%范圍以內(nèi)。說明更密的網(wǎng)格對模型的計算結果沒有顯著的影響，但是工況C3的計算時間遠遠超過工況C2的計算時間。為了節(jié)約計算成本，最終選取工況C2的網(wǎng)格劃分方案(160×200×80)作為本次計算的網(wǎng)格方案。Cater等[12]，Xie等[13]也運用LES方法對Tominaga等[4]的物理模型試驗工況進行驗證分析，他們采用的網(wǎng)格數(shù)量分別為500×320×128和192×384×96，分別是本次模擬網(wǎng)格數(shù)的8倍和2.7倍。

圖2 不同網(wǎng)格精度下流向(x方向)時均流速分布

3 計算結果與分析

計算約60個流動周期(一個流動周期T=10H/Um)后流場達到充分發(fā)展狀態(tài)，開始統(tǒng)計流場的時均流動特性，統(tǒng)計時間為50個流動周期。U、V、W分別為x、y、z方向的時均流速。為對WMLES方法在復式斷面明渠流動模擬中的表現(xiàn)進行全面的評價，將模擬計算結果與Tominaga等[4]的試驗測量數(shù)據(jù)，Xie等[13]基于動態(tài)亞格子模型的大渦模擬結果和Cater等[12]基于Smagorinsky模型的大渦模擬結果進行比較。

3.1 主流速等值線

以最大主流速Umax對主流速進行無量綱化，計算得到的主流速等值線分布如圖3(d)所示。計算得到的最大主流速出現(xiàn)在主槽中心z/H=0.67附近。由圖3可以看出WMLES計算結果與Tominaga等[4]的實測值吻合。在主槽與灘地交界處，由于二次流的影響，主流速等值線向主槽側(cè)的自由水面凸起。同時可以觀察到主流速在接近自由水面的位置速度減小的現(xiàn)象。在斷面上由于強烈的二次流動，渦體將動量較高的水體向水槽兩邊輸運，導致主流速等值線在拐角附近發(fā)生較大程度的彎曲。在自由水面附近靠近壁面處受壁面的影響，主流速等值線也出現(xiàn)較大程度的彎曲。與Cater等[12-13]的模擬結果相比，WMLES模擬的斷面主流速在灘槽交界面處與實測值更接近。

圖3 橫斷面主流速U/Umax分布

3.2 二次流

圖4(d)是通過WMLES計算得到的斷面流速矢量圖，圖中V、W已通過最大主流速Umax進行無量綱化。Prandtl根據(jù)形成原因?qū)⒍瘟鞣譃閮深悾旱谝活怭randtl二次流常產(chǎn)生在彎曲河道中，離心力的作用使水流偏離主流方向，此類二次流強度較大；第二類Prandtl二次流常發(fā)生在復式斷面河道邊角附近，由于雷諾應力的各向異性導致水流偏離主流方向，此類二次流強度相對較小。與試驗中觀察到的現(xiàn)象一致，圖4中可以看出在灘地與主槽交界處產(chǎn)生較強的二次流，存在一對旋轉(zhuǎn)方向相反的渦延伸至自由水面，靠近主槽的渦稱為深槽渦，靠近灘地的渦稱為淺灘渦。另外，計算結果顯示在主槽靠近壁面處有一漩渦，稱為自由水面渦。自由水面渦的影響范圍超過主槽寬度的一半，它將緩慢流動的流體從側(cè)壁附近輸送到水槽的中心部分，這也是主流速最大值出現(xiàn)在自由水面以下的原因之一。與試驗測量結果[4]和前人計算結果[12-13]相比，WMLES方法模擬計算得到的橫斷面漩渦結構在位置和方向上都較為一致。

圖4 橫斷面流速矢量

3.3 床面切應力

圖5 橫斷面床面切應力分布

3.4 紊動強度

水流流速的脈動強弱程度通常可以由紊動強度來反映，紊動強度用脈動流速的均方根來表示：

(7)

(8)

式中：urms、vrms、wrms分別為流向、展向和垂向紊動強度。

圖6 流向紊動強度urms/uτ分布

圖7 展向紊動強度vrms/uτ分布

圖8 垂向紊動強度wrms/uτ分布

圖9 紊動能分布

3.5 雷諾應力

圖10 雷諾應力斷面分布

圖11 雷諾應力斷面分布

由圖10可知，-〈u′v′〉的值與速度梯度?U/?y有關，代表了灘槽交界面處動量交換的強度。在主槽內(nèi)-〈u′v′〉的絕對值出現(xiàn)最大值，同時可以發(fā)現(xiàn)在灘槽交界的拐角處y/H>2.5的部分區(qū)域-〈u′v′〉的值為正值。表明動量交換在主槽內(nèi)最強烈，水體在灘地靠近床面附近動量從主槽向灘地輸運。在主槽靠近壁面?zhèn)?〈u′v′〉的值逐漸增大，這與該處存在自由表面渦有關。

由圖11(b)可知，在主槽內(nèi)靠近床面-〈u′w′〉的值逐漸增大，主槽與灘地交界處的左側(cè)出現(xiàn)負值，這與試驗結果是一致的。在灘地內(nèi)-〈u′w′〉模擬值與試驗測量值相比偏小，兩者的變化規(guī)律一致。

4 結 語

本文利用WMLES方法對水深比為0.5的復式斷面明渠水流進行三維大渦模擬研究。與傳統(tǒng)LES方法相比，WMLES無須在邊界層內(nèi)對網(wǎng)格加密，很大程度上降低了計算成本。就本文算例而言，與Cater等[12]和Xie等[13]的LES數(shù)值模型相比，利用WMLES方法使得計算網(wǎng)格數(shù)量分別減少到其網(wǎng)格數(shù)量的1/3和1/8。WMLES計算結果與試驗測量值[4]和LES計算結果[12-13]的比較表明,WMLES能夠準確模擬平均流速、床面切應力、紊動強度、紊動能和雷諾應力分布，并且能夠準確模擬出斷面內(nèi)的二次漩渦結構。本文驗證了WMLES方法對于復式斷面明渠流動模擬問題的適用性和可靠性，為后續(xù)研究打下基礎。