復(fù)雜地形下帶自由表面水流的分離渦模擬

張景新

（1.上海交通大學(xué)水動(dòng)力學(xué)教育部重點(diǎn)實(shí)驗(yàn)室，上海 200240；2.上海河口海岸科學(xué)研究中心河口海岸交通行業(yè)重點(diǎn)實(shí)驗(yàn)室，上海 201201）

文章編號(hào)：1001?246X（2015）05?0561?11

復(fù)雜地形下帶自由表面水流的分離渦模擬

張景新1，2

（1.上海交通大學(xué)水動(dòng)力學(xué)教育部重點(diǎn)實(shí)驗(yàn)室，上海 200240；
2.上海河口海岸科學(xué)研究中心河口海岸交通行業(yè)重點(diǎn)實(shí)驗(yàn)室，上海 201201）

將分離渦模型（DES），即一種RANS和LES的混合模型，應(yīng)用于帶自由表面的地表水流運(yùn)動(dòng)，建立一套數(shù)值仿真模型.模型基于有限體積法，水平面內(nèi)采用非結(jié)構(gòu)計(jì)算網(wǎng)格，垂向?yàn)榻Y(jié)構(gòu)化網(wǎng)格，對(duì)流項(xiàng)離散格式采用二階TVD格式，并行基于OpenMP語言庫.算例表明DES模型有助于揭示復(fù)雜地形條件下帶自由表面水流的大渦擬序結(jié)構(gòu).

分離渦；擬序結(jié)構(gòu)；自由表面流動(dòng)；動(dòng)壓模型

0 引言

天然地表水流涉及洋流、河流、湖流等，對(duì)其水動(dòng)力特征的研究有助于水利工程、水環(huán)境工程等的順利開展.基于雷諾平均（RANS）的數(shù)學(xué)模型得到了廣泛的應(yīng)用，該類模型可描繪水流的時(shí)均流動(dòng)特征，且計(jì)算效率較高.高精度數(shù)值模擬技術(shù)（LES、DNS）有助于流動(dòng)細(xì)觀結(jié)構(gòu)的揭示，可描繪出水流運(yùn)動(dòng)的大渦擬序結(jié)構(gòu).但鑒于地表水流的時(shí)空尺度及復(fù)雜的邊界幾何形狀，無論LES還是DNS近期都較難付諸于實(shí)際工程的數(shù)值模擬研究.結(jié)合LES與RANS兩者的混合模型，近年來逐漸被用于工業(yè)流動(dòng)的數(shù)值模擬，將其應(yīng)用于地表水流運(yùn)動(dòng)的數(shù)值模擬，有助于進(jìn)一步研究該類水流運(yùn)動(dòng)的若干機(jī)理問題.分離渦（DES）模型是一種RANS/LES的混合模型，由Spalart等［1］提出，后續(xù)有了長足進(jìn)展［2－5］.本文將分離渦（DES）模型拓展至水動(dòng)力學(xué)數(shù)值模擬，針對(duì)復(fù)雜地形條件下的帶自由表面水流運(yùn)動(dòng)開展研究.

大尺度帶自由表面的地表水流運(yùn)動(dòng)，依據(jù)流動(dòng)的水平尺度遠(yuǎn)大于垂向運(yùn)動(dòng)尺度，通常采用靜壓假定來計(jì)算壓力.但對(duì)于地形陡變、自由表面坡陡較大、局部射流等情況，需要補(bǔ)充動(dòng)壓.Casulli［6－7］建立了完全的動(dòng)壓模型，將壓力分為靜壓和動(dòng)壓.靜壓的求解與常用淺水方程求解方法相似，在靜壓求解的基礎(chǔ)上再通過求解壓力泊松方程計(jì)算動(dòng)壓，可將其視為一種壓力的預(yù)估—校正法.Jankowski［8］詳細(xì)地論述了預(yù)估—校正法，即首先求解靜壓作用下的流場，稱為預(yù)估步；在此基礎(chǔ)上，通過求解動(dòng)壓滿足的泊松（Poisson）方程，進(jìn)一步更新動(dòng)壓作用下的流場.Li［9］將預(yù)估校正法用于表面重力波的模擬.Kocyigit等［10］和Chen等［11］在笛卡爾坐標(biāo)系下求解了三維非靜壓模型.Fringer等［12］建立了大洋潮汐流的非靜壓模型.非靜壓模型正逐漸地被應(yīng)用到大尺度地表水流運(yùn)動(dòng)的數(shù)值模擬中.大尺度地表水流運(yùn)動(dòng)模擬中，網(wǎng)格尺度通常較大，而粗網(wǎng)格不能描述局部地形的陡變，某種意義上而言，地形被數(shù)值平均了.而DES要求較高分辨率的計(jì)算網(wǎng)格，此時(shí)，陡變地形得以刻畫.靜壓假定模型的模擬精度降低，動(dòng)壓的引入是必要的.

文章第1部分介紹自由表面水流的DES模型及數(shù)值求解方法；第2部分為模型算例，驗(yàn)證模型精度，討論DES的相關(guān)技術(shù)問題.

1 數(shù)學(xué)模型

1.1 控制方程描述

基本控制方程為考慮地球自轉(zhuǎn)效應(yīng)（柯氏力）的三維自由表面水流運(yùn)動(dòng)方程，其中壓力分解為靜壓ph＝ρg（ζ－z）和動(dòng)壓pn，連續(xù)性方程和動(dòng)量方程如下其中g(shù)為重力加速度，柯氏力f＝2ωsin?，?為地球緯度，ω為地球旋轉(zhuǎn)角速度，ζ為水位值．若忽略動(dòng)壓作用，則上述控制方程退化為常用的靜壓方程，只包含連續(xù)性方程和水平動(dòng)量方程．流場變量 qx＝Du，qy＝Dv，qz＝Dw，qσ＝Dω~，σ坐標(biāo)系下的垂向速度

上述控制方程中采用了垂向σ坐標(biāo)［13］變換，可擬合不規(guī)則床面及波動(dòng)水面.

1.2 湍流模型

湍流模型采用S?A［14］一方程模型，計(jì)算變量υ~通過如下輸運(yùn)方程獲得：

渦粘性系數(shù)υt通過υt＝的關(guān)系式計(jì)算得到．其中fv1＝χ3/（χ3＋），χ≡/υ．υ為分子運(yùn)動(dòng)學(xué)粘性系數(shù)．為渦量值，而＝＋（υ~/k2 d2）fv2，其中fv2＝1－χ/（1＋χfv1）．方程（6）中的耗散項(xiàng)系數(shù)fw由fw＝g［（1 ＋）/（g6＋）］1/6計(jì)算得到，其中各參數(shù)計(jì)算公式為g＝r＋cw2（r6－r），r≡/（）．S?A模型中的各計(jì)算參數(shù)分別為cb1＝0.135 5，＝2/3，cb2＝0.622，κ＝0.41，cw1＝cb1/κ2＋（1＋cb2），cw2＝0.3，cw3＝2.0，cv1＝7.1．

S?A模型控制方程中耗散項(xiàng)中的特征長度d為距固壁的距離，表征了受固壁限制湍流運(yùn)動(dòng)的特征長度. DES模型通過修改該變量實(shí)現(xiàn)RANS至LES的轉(zhuǎn)換，以替代方程中的d，即d~＝min（d，CDESΔ）．表達(dá)式中Δ＝max（4A/π，Δz），A為水平網(wǎng)格面積，系數(shù)CDES取0.65［15］.

求解變量υ~的邊界條件涉及固壁、自由表面.采用S?A模型，固壁邊界條件滿足ν~t＝0．自由表面的值可設(shè)定為常數(shù)，某些參考值可見相應(yīng)的研究成果，如設(shè)定為3ν～5ν［14，16］，3ν［17－18］，也有研究取值范圍為［0，0.1ν］［2］．以上的研究成果多針對(duì)氣動(dòng)力學(xué).Yue等［19］針對(duì)氣水兩相流，采用LES模擬了明渠流動(dòng)，氣水交界面基于VOF法，結(jié)果表明自由表面處的渦粘性系數(shù)為幾倍的分子粘性系數(shù).

1.3 DES的“灰區(qū)”

DES模型由RANS至LES的轉(zhuǎn)換過程中，存在一個(gè)“灰區(qū)”.其實(shí)質(zhì)是模擬的雷諾應(yīng)力不足 （Model?Stress Depletion，MSD），根本原因是RANS模型不能激發(fā)充分的湍流脈動(dòng)量.Spalart等［3－4］為克服DES的這一缺陷，先后提出了DDES（Delayed DES）和IDDES（Improved Delayed DES）.另一種解決辦法是在RANS/LES界面處通過數(shù)值方法產(chǎn)生湍流脈動(dòng)量［20］，也可明顯改善MSD問題.區(qū)域分離渦（Zonal?DES，ZDES）模型采用了另一種思路.與DES模型比較，DES模型由RANS至LES的轉(zhuǎn)換是逐漸過渡的，即方程中的湍流特征尺度的變化是連續(xù)的.ZDES模型中的特征尺度當(dāng)計(jì)算域由RANS轉(zhuǎn)換至LES計(jì)算域時(shí)，直接以LES模型網(wǎng)格尺度代替，即Δ＝（AΔz）1/3（A為本模型中水平網(wǎng)格面積）．與此同時(shí)，相關(guān)的計(jì)算參數(shù)fv1，fv2和fw修改為

ZDES模型中，流動(dòng)由RANS轉(zhuǎn)換至LES時(shí)，特征尺度的變化及上述各參數(shù)的設(shè)定，使得模型迅速轉(zhuǎn)化為LES，可激發(fā)出流動(dòng)的脈動(dòng)量，從而獲得更多的湍流脈動(dòng).Breuer等［21］以ZDES模擬了繞平板的強(qiáng)分離流動(dòng)，獲得了與LES模擬吻合較好的結(jié)果.Deck［22－24］通過對(duì)粘性系數(shù)等的模擬結(jié)果比較，驗(yàn)證了ZDES在克服MSD問題方面的可行性.同時(shí)，該模型實(shí)現(xiàn)簡便，并未增加額外的計(jì)算量.

1.4 數(shù)值方法

1.4.1 坐標(biāo)變換

對(duì)于正交網(wǎng)格，控制體各面上的法向?qū)?shù)可由該面兩側(cè)單元形心處的物理量沿兩點(diǎn)連線直接求導(dǎo)［7，12］.對(duì)于非正交網(wǎng)格系統(tǒng)，控制面上的導(dǎo)數(shù)計(jì)算可借助局部坐標(biāo)系求解.本文模型，在控制面上引入局部坐標(biāo)系（ξ，η），笛卡兒坐標(biāo)系下的導(dǎo)數(shù)計(jì)算通過鏈導(dǎo)法則，轉(zhuǎn)化為局部坐標(biāo)系下的相應(yīng)計(jì)算：

其中J＝xξyη－xηyξ為坐標(biāo)系轉(zhuǎn)換系數(shù)的雅可比行列式．局部坐標(biāo)ξ由單元形心指向鄰單元形心，η沿兩單元界面，方向?yàn)棣畏较蚰鏁r(shí)針旋轉(zhuǎn)，該坐標(biāo)系見圖1.

圖1 局部坐標(biāo)系Fig.1 Local coordinates on a control cell surface

1.4.2 靜壓模型求解

數(shù)值模型的時(shí)間積分采用半隱格式，通過參數(shù)θ實(shí)現(xiàn)［25］．參數(shù)θ的取值范圍已有相關(guān)報(bào)道［10－11，25］．本文模型計(jì)算中θ取值0.5.

連續(xù)性方程的數(shù)值離散為

動(dòng)量方程分別離散為

上述離散方程中F為顯式方式離散項(xiàng)，包括對(duì)流項(xiàng)、水平粘性擴(kuò)散項(xiàng)及柯氏力等.

首先，離散方程中僅保留靜壓項(xiàng)ph，引入臨時(shí)變量，，，，η?，單獨(dú)靜壓作用下連續(xù)性方程的離散形式為

動(dòng)量方程離散形式為

采用矩陣表達(dá)式，將上述離散方程改寫為緊湊的表達(dá)形式，

上述表達(dá)式中的矩陣分別如下

其中，Ai是三對(duì)角形式的系數(shù)矩陣.

將方程（18）、（19）代入方程（17），可得到關(guān)于待求水位變量ζ?的控制方程

將水位的控制方程（21）在計(jì)算單元水平面內(nèi)積分，利用高斯積分定理，得到離散方程

方程（22）中，符號(hào)f代表水平單元的邊，ΔSi代表水平單元面積，Δlis為第s條邊的邊長，NS為各水平單元的總邊數(shù)，（cosαis，sinαis）為第s條邊的法向單位向量．單元各邊建立局部坐標(biāo)系（圖1），表達(dá)式（22）中變量在笛卡爾坐標(biāo)系（x，y）下的導(dǎo)數(shù)計(jì)算可借助局部坐標(biāo)系（ξ，η）下的相應(yīng)計(jì)算獲得，具體表達(dá)式

將表達(dá)式（23）寫成如下的緊湊形式其中各系數(shù)為

代數(shù)方程組（24）采用雙共軛梯度法（Bi?CGSTAB）求解，可獲得新時(shí)刻的水位變量ζ?；新時(shí)刻的流速變量，通過方程（18）和（19）進(jìn)一步求解得到.如果模型僅做靜壓計(jì)算，上述的求解變量即為n＋1時(shí)步的最終變量，一個(gè)計(jì)算循環(huán)完成；若模型設(shè)定為非靜壓模式，需要進(jìn)一步求解動(dòng)壓項(xiàng).

1.4.3 動(dòng)壓模型求解

動(dòng)壓模型在靜壓求解的基礎(chǔ)上進(jìn)行，動(dòng)壓pnn＋1作用下的控制方程

n＋1時(shí)刻更新的流場須滿足連續(xù)性條件，即流速變量滿足連續(xù)性方程（1）.方程（1）中不顯含變量qz，首先需要將方程（1）改寫為將方程（25）、（26）和（27）代入方程（28），得到如下關(guān)于動(dòng)壓pn＋1n的泊松方程上述控制方程中關(guān)于動(dòng)壓pn＋1

將方程（25）、（26）和（27）代入方程（28），得到如下關(guān)于動(dòng)壓pn＋1n的泊松方程

其中各系數(shù)的計(jì)算表達(dá)式

1.4.4 TVD離散格式

本模型采用二階TVD數(shù)值格式離散動(dòng)量方程的對(duì)流項(xiàng).有限體積法中，界面處的任一物理量〈φ〉f通過如下插值法獲得

其中符號(hào)f代表界面，φD和φC分別為該界面順風(fēng)向和迎風(fēng)向單元形心處的相應(yīng)變量．通過通量限制器ψ（rf），實(shí)現(xiàn)二階的TVD格式．ψ（rf）是變量rf的相關(guān)函數(shù)，相關(guān)的計(jì)算可參考文獻(xiàn)Darwish等［26］.本模型引入了如下的通量限制器：

SUPERBEE ψ（rf）＝max（0，min（1，2rf），min（2，rf））；

MINMOD ψ（rf）＝max（0，min（1，rf））；

OSHER ψ（rf）＝max（0，min（2，rf））；

MUSCL ψ（rf）＝（rf＋ rf）/（1＋ rf）；

VAN LEER ψ（rf）＝（rf＋ rf）/（1＋rf）；

SWEBY ψ（rf）＝max（0，min（1，1．5rf），min（1.5，rf））；

QUICK ψ（rf）＝max（0．0，min（2rf，（3.0＋rf）/4．0，2．0））；

UMIST ψ（rf）＝max（0．0，min（2rf，（1．0＋3rf）/4．0，（3．0＋rf）/4．0，2．0）；

MC ψ（rf）＝max（0．0，min（（1．0＋rf）/2．0，2．0，2．0rf））．

模型中的上述TVD格式均做了嚴(yán)格驗(yàn)證，下文計(jì)算采用OSHER格式.

1.4.5 湍流S?A模型的離散

控制方程（6）采用有限體積法（FVM）離散求解，具體離散格式如下

其中F為顯式方式離散項(xiàng)，包括對(duì)流項(xiàng)、水平粘性擴(kuò)散項(xiàng)等；生成項(xiàng)與耗散項(xiàng)采用隱格式求解.上述離散方程的具體求解過程及所形成的三對(duì)角方程組的求解同動(dòng)量方程的相應(yīng)求解方法.

2 模型應(yīng)用

將上述數(shù)值模型應(yīng)用于自由表面水流運(yùn)動(dòng)的模擬，重點(diǎn)考察小尺度渦結(jié)構(gòu)，獲得流動(dòng)的細(xì)觀特征.本文針對(duì)系列沙丘地形下的水流運(yùn)動(dòng)開展小尺度渦結(jié)構(gòu)的模擬研究，模型計(jì)算細(xì)節(jié)及結(jié)果分析如下文所述.

沙丘地形下帶自由表面水流運(yùn)動(dòng)

單一沙丘的幾何形狀與Balachandar等［27］的實(shí)驗(yàn)相同.該實(shí)驗(yàn)共布設(shè)了22個(gè)相同的沙丘，詳細(xì)測量了第17個(gè)沙丘范圍內(nèi)的流場.圖2描繪了單個(gè)沙丘的幾何形狀及測量點(diǎn)位置.本文模擬了5個(gè)沙丘組成的系列地形條件下的水流運(yùn)動(dòng)，重點(diǎn)考察局部流動(dòng)的渦結(jié)構(gòu)，討論了流動(dòng)的沿程演化等問題.覆蓋沙丘地形的計(jì)算域網(wǎng)格尺度設(shè)計(jì)滿足LES的要求，向上下游開邊界網(wǎng)格尺度逐漸增加，直至完全處于RANS計(jì)算范疇.依據(jù)網(wǎng)格尺度設(shè)計(jì)的計(jì)算域分區(qū)見圖3，其中RANS計(jì)算域不僅覆蓋近底床區(qū)域，還延伸至上下游開邊界處的全水深區(qū)域.對(duì)于LES，開邊界的入流條件，即來流的湍流脈動(dòng)信息的充足與否對(duì)計(jì)算影響明顯.而本算例中開邊界設(shè)置為RANS計(jì)算域，沒有提供流場脈動(dòng)量信息，勢必對(duì)湍流小尺度渦結(jié)構(gòu)的模擬產(chǎn)生影響.本算例設(shè)計(jì)了5個(gè)連續(xù)的沙丘地形，預(yù)期前部沙丘地形誘導(dǎo)湍流脈動(dòng)，為后續(xù)流場提供湍流脈動(dòng)成分.與LES模擬中采用數(shù)值生成脈動(dòng)量的方法相比較，可視其為湍流脈動(dòng)產(chǎn)生的物理方法.若僅研究沙丘地形下的流動(dòng)特征，布置5個(gè)相同的沙丘地形，顯然效率不高.天然水流（如河流），空間尺度不允許全計(jì)算域采用LES，而只能局部關(guān)注區(qū)域采用LES，其它區(qū)域采用RANS.天然河流等無論地形，還是岸線均是不規(guī)則的，這些不規(guī)則的幾何形狀影響流動(dòng)，可激發(fā)出湍流的脈動(dòng)量.正如本算例所設(shè)置的地形，可作為脈動(dòng)量生成的物理因素.上游來流中脈動(dòng)量的重要性，通過以下模擬結(jié)果進(jìn)行分析.

圖2 沙丘地形及測點(diǎn)位置Fig.2 Geometry of dunes

圖3 計(jì)算域及RANS/LES區(qū)域劃分Fig.3 Schematic of computational zones

本算例以上游水深L和自由表面來流速度U0為特征量，控制參數(shù)為Re＝5．7×104，且Fr＝0.44．水平網(wǎng)格最高分辨率5mm，垂向網(wǎng)格設(shè)計(jì)滿足近壁面第一層網(wǎng)格中心點(diǎn)z＋≈1，繼而垂直向上以1.15的伸展率逐漸增加，直至達(dá)到5mm.明渠流動(dòng)的大渦尺度受水深這一特征尺度限制，相關(guān)研究［28］給出了滿足LES計(jì)算需要的網(wǎng)格尺度與水深的關(guān)系，可作為明渠水流LES模擬網(wǎng)格設(shè)計(jì)的參考.本文網(wǎng)格尺度與水深關(guān)系為L/24（L為上游水深值）.DES計(jì)算中RANS和LES的區(qū)域劃分以網(wǎng)格尺度控制，本算例中5mm的計(jì)算網(wǎng)格控制了RANS/LES的分區(qū)界面位于O（10z＋），而該區(qū)域是近壁湍流發(fā)展的活躍區(qū)域［20］.

模型首先以靜壓模式運(yùn)行至流動(dòng)穩(wěn)定狀態(tài)，再切換至動(dòng)壓模式計(jì)算至穩(wěn)定狀態(tài)，繼而由RANS模式轉(zhuǎn)換至DES模式.在DES模型運(yùn)行過程中，記錄若干點(diǎn)流動(dòng)變量的計(jì)算值，當(dāng)達(dá)到統(tǒng)計(jì)意義上的穩(wěn)定狀態(tài)后，再繼續(xù)運(yùn)行15個(gè)大渦周期（L/uτ），uτ為上游流場的平均摩阻流速.該15個(gè)大渦周期時(shí)段的流場模擬信息用于流動(dòng)分析，數(shù)據(jù)采樣頻率100Hz.

時(shí)均流動(dòng)分析

對(duì)于本文系列沙丘地形下的流動(dòng)，RANS模型獲得穩(wěn)定的定常流動(dòng)，然而DES模型預(yù)報(bào)了小尺度渦的脈動(dòng)，流動(dòng)是非定常的，僅存在統(tǒng)計(jì)意義上的穩(wěn)定狀態(tài).將流場的時(shí)間序列模擬結(jié)果作時(shí)均分析，可獲得時(shí)均流場信息.

實(shí)驗(yàn)測量了圖2所示6條垂線位置處的流速時(shí)間序列，并給出了相應(yīng)的流場特征.6個(gè)測點(diǎn)相對(duì)于沙丘峰值點(diǎn)的位置分別為x/h＝2，4，5，6，12和18.本文算例分別記錄了第一、第三和第五個(gè)沙丘相應(yīng)的6個(gè)測點(diǎn)的流場信息，用于模型驗(yàn)證.圖4給出了流向速度的對(duì)比，其中符號(hào)為實(shí)測值，線分別對(duì)應(yīng)于三個(gè)沙丘的計(jì)算值.分析6個(gè)測量點(diǎn)的模擬結(jié)果與實(shí)驗(yàn)結(jié)果，第五個(gè)沙丘范圍內(nèi)的相應(yīng)計(jì)算值與實(shí)測值吻合最好，精度最高.第一個(gè)沙丘范圍內(nèi)的相應(yīng)計(jì)算值精度最低.而RANS模型的結(jié)果顯示幾個(gè)沙丘范圍內(nèi)的計(jì)算值基本一致，并未顯示出沿流向流場的明顯變化.分析本文DES模型的網(wǎng)格設(shè)計(jì)，入口為RANS計(jì)算域，對(duì)于其后的LES計(jì)算域而言，輸入的湍流脈動(dòng)成分不足，導(dǎo)致湍流場發(fā)展不充分.當(dāng)水流流經(jīng)沙丘時(shí)，該局部地形變化將誘發(fā)流場脈動(dòng)，該脈動(dòng)量可視為其后流場的輸入條件，故模擬結(jié)果顯示湍流場沿程逐漸發(fā)展，至第五個(gè)局部沙丘范圍時(shí)，模擬結(jié)果與實(shí)測值吻合精度已明顯提高.

圖4 流向速度驗(yàn)證Fig.4 Stream velocity at different sites

計(jì)算時(shí)均Reynolds應(yīng)力〈－u′w′〉，并與實(shí)測值比較，圖5給出了6個(gè)測量點(diǎn)處兩者的比較.圖中符號(hào)為實(shí)驗(yàn)測量值，線分別為第一、第三和第五個(gè)沙丘范圍內(nèi)相應(yīng)位置的計(jì)算值.Reynolds應(yīng)力〈－u′w′〉的極值大致出現(xiàn)在z/L＝0，數(shù)值模擬結(jié)果與實(shí)測值吻合較準(zhǔn)確.第一個(gè)沙丘范圍內(nèi)的計(jì)算值明顯低于實(shí)測值，第三和第五個(gè)沙丘范圍內(nèi)的相應(yīng)計(jì)算值與實(shí)測值符合的較好.該結(jié)果再次驗(yàn)證了入流湍流脈動(dòng)量對(duì)計(jì)算結(jié)果有明顯的影響.系列沙丘地形條件下，水流流經(jīng)某個(gè)局部沙丘，將激發(fā)湍流脈動(dòng)量.對(duì)后續(xù)流場，該脈動(dòng)量可視為入流條件，故模擬精度沿流向逐漸提高.圖5顯示第三個(gè)沙丘地形模擬結(jié)果較之第五個(gè)沙丘范圍內(nèi)的相應(yīng)計(jì)算結(jié)果精度更高，與湍流場逐漸充分發(fā)展的預(yù)期存在偏差，這可能在于流動(dòng)經(jīng)過這一系列沙丘后，尚未達(dá)到統(tǒng)計(jì)意義上的穩(wěn)定狀態(tài).相應(yīng)的實(shí)驗(yàn)測量均取在第17個(gè)沙丘范圍內(nèi)完成，也是考慮到湍流場的充分發(fā)展需要一定的沿程距離.

圖5 Reynolds應(yīng)力模擬與比較Fig.5 Predicted mean Reynolds stress and experimental data

瞬時(shí)流動(dòng)分析

圖6 Q的瞬時(shí)等值云圖（灰度由渦量渲染）Fig.6 Snapshot of criterion Q isosurfaces contoured by vorticity

較之RANS模型，DES可模擬小尺度的渦結(jié)構(gòu)，給出流動(dòng)的細(xì)觀特征.湍流場的復(fù)雜渦結(jié)構(gòu)可以用速度梯度張量的第二不變量描述，即變量Q ＝（ΩijΩij－SijSij）/2＞0，其中Sij和Ωij分別是速度梯度張量ΔV的對(duì)稱和反對(duì)稱部分．圖6給出了Q＝5的瞬時(shí)云圖，并以渦量值渲染.小尺度渦沿流向逐漸發(fā)展，模擬結(jié)果逐漸顯示出大渦模擬的內(nèi)容.上游局部沙丘激發(fā)出湍流脈動(dòng)，從而促使湍流場充分發(fā)展.這一算例的設(shè)計(jì)，入流邊界條件并未包含湍流脈動(dòng)信息，但流動(dòng)的發(fā)展，特別是地形變化引起的水流分離運(yùn)動(dòng)，極大地激發(fā)出了脈動(dòng)量，使得小尺度渦結(jié)構(gòu)向下游越發(fā)清晰，湍流運(yùn)動(dòng)逐漸充分演化.

考察湍流場的沿程演化，分別在第一、第三和第五個(gè)沙丘范圍內(nèi)取一橫斷面，各斷面相對(duì)沙丘峰值點(diǎn)位置為0.2λ（λ為單個(gè)沙丘全長）.圖7描述了各個(gè)斷面內(nèi)的瞬時(shí)流速矢量分布及渦量分布.瞬時(shí)流速矢量及渦量分布顯示了不同斷面模擬結(jié)果中小尺度渦的發(fā)展情況，即湍流場充分發(fā)展的程度.數(shù)值模擬結(jié)果所包含的渦結(jié)構(gòu)是否充分，是大渦模擬成功與否的一個(gè)標(biāo)志.結(jié)果顯示來流所含湍流場脈動(dòng)信息對(duì)湍流運(yùn)動(dòng)模擬、大渦結(jié)構(gòu)的捕獲非常重要.

圖7 不同斷面內(nèi)的瞬時(shí)流場（x/λ＝0.2）Fig.7 Instantaneous flows at different cross?sections（x/λ＝0.2，withλthe dune length）

通過時(shí)均處理，可獲得時(shí)均流場信息，進(jìn)一步提取瞬時(shí)流場脈動(dòng)量，同樣可分析湍流場的沿程發(fā)展.圖8給出了中垂面（y＝0，x?z）內(nèi)瞬時(shí)脈動(dòng)速度（u′，w′）的分布，分別對(duì)應(yīng)于第一、第三和第五個(gè)沙丘范圍.與前述分析類似，沿流向，脈動(dòng)速度逐漸增大，湍流運(yùn)動(dòng)逐漸發(fā)展.大渦模擬中通常采用數(shù)值方法在來流中生成脈動(dòng)速度，本算例的入流開邊界處于RANS計(jì)算域，缺少脈動(dòng)量的輸入.而局部地形變化激發(fā)的湍流脈動(dòng)量，可視為其后流場模擬的入流條件，故極大地提高了模擬精度.若針對(duì)天然河流采用該計(jì)算方法，僅在重點(diǎn)關(guān)注區(qū)域采用LES，開邊界局部仍采用RANS，將簡化邊界條件的設(shè)定.同時(shí)，天然河道的不規(guī)則幾何特征將極大地激發(fā)湍流脈動(dòng)，從而為LES模擬區(qū)域提供充分的脈動(dòng)信息.

圖8 中垂面（y＝0，x?z）內(nèi)瞬時(shí)脈動(dòng)速度場（u′，w′）的沿程演化Fig.8 Instantaneous velocity vectors（u′，w′）in the central x?z plane

3 結(jié)論

將靜壓假定條件下的自由表面淺水動(dòng)力學(xué)模型拓展至考慮動(dòng)壓的數(shù)學(xué)模型，極大地拓寬了模型的應(yīng)用范圍.通過修改一方程湍流模型（S?A模型），建立了自由表面水流運(yùn)動(dòng)的DES數(shù)值模型.建立的數(shù)值模型基于FVM，采用非結(jié)構(gòu)網(wǎng)格，并基于OpenMP技術(shù)實(shí)現(xiàn)了軟件的并行化.DES模型較之RANS模型，在小尺度渦結(jié)構(gòu)的捕獲方面具有一定的優(yōu)勢，可用于流動(dòng)細(xì)觀特征的研究.通過RANS與DES模型算例的對(duì)比，可知DES模型能夠提供更豐富的小尺度旋渦信息，而這些在RANS模型中均被時(shí)均化處理方法所掩蓋.從流動(dòng)細(xì)觀入手，通過適當(dāng)?shù)慕y(tǒng)計(jì)方法，進(jìn)而獲得流動(dòng)的宏觀特征，較之直接的時(shí)均化方法研究，可獲得更加豐富的流場信息.

DES模型局部采用RANS，強(qiáng)分離區(qū)采用LES，對(duì)于高雷諾數(shù)、復(fù)雜邊界幾何條件的流動(dòng)模擬是兼顧了計(jì)算效率和計(jì)算精度的一種可行辦法.但也存在一些問題，最突出的一個(gè)就是“灰區(qū)”問題，即計(jì)算模型由RANS轉(zhuǎn)換至LES時(shí)，所輸入的湍流脈動(dòng)量不足.對(duì)于開邊界處于RANS區(qū)的情況，入流脈動(dòng)量不足也將導(dǎo)致計(jì)算結(jié)果精度降低.

天然大尺度的自由表面水流運(yùn)動(dòng)，多為高雷諾數(shù)、復(fù)雜幾何邊界條件，LES的實(shí)現(xiàn)尚有一定的難度.而采用LES/RANS混合模型，既可實(shí)現(xiàn)局部關(guān)注區(qū)域的高分辨率模擬，又可大大降低計(jì)算量，是一種可推廣應(yīng)用的數(shù)值方法.

［1］ Spalart PR，Jou W H，Strelets M，Allmaras SR.Comments on the feasibility of LES forwings，and on a hybrid RANS/LES approach［M］∥Liu C，Liu Z.Advances in DNS/LES.Greyden Press，1997.

［2］ Spalart PR，Allmaras SR.A one?equation turbulencemodel for aerodynamic flows［J］.La Recherche Aérospatiale，1994，1：5－21.

［3］ Spalart PR，Deck S，Shur M L，Squires K D，Strelets M K，Travin A.A new version of detached?eddy simulation，resistant to ambiguous grid densities［J］.Theor Comput Fluid Dyn，2006，20：181－195.

［4］ Shur K L，Spalart PR，Strelets M K，Travin A K.A hybrid RANS?LES approach with delayed?DES and wall?modelled LES capabilities［J］.International Journal of Heat and Fluid Flow，2008，29：1638－1649.

［5］ Spalart PR.Detached?eddy simulation［J］.Annu Rev Fluid Mech，2009，41：181－202.

［6］ Casulli V.A semi?implicit finite difference method for non?hydrostatic，free?surface flows［J］.International Journal for Numerical Methods in Fluids，1999，30：425－440.

［7］ Casulli V，Zanolli P.Semi?implicit numericalmodeling of nonhydrostatic free?surface flows for environmental problems［J］. Mathematical and Computer Modeling，2002，36：1131－1149.

［8］ Jankowski JA.A non?hydrostatic model for free surface flows［D］.Germany：University of Hannover，1999.

［9］ Li B，F(xiàn)leming C A.Three?dimensionalmodel of Navier?Stokes equations for water waves［J］.Journal of Waterway，Port，Coastal and Ocean Engineering，2001，127（1）：16－25.

［10］ Kocyigit M B，F(xiàn)alconer R A，Lin B.Three?dimensional numericalmodeling of free surface flowswith non?hydrostatic pressure ［J］.International Journal for Numerical Methods in Fluids，2002，40：1145－1162.

［11］ Chen X J.A fully hydrodynamic model for three?dimensional，free?surface flows［J］.International Journal for Numerical Methods in Fluids，2003，42：929－952.

［12］ Fringer O B，Gerritsen M，Street R L.An unstructured?grid，finite?volume，nonhydrostatic，parallel coastal ocean simulator ［J］.Ocean Modeling，2006，14：139－173.

［13］ Phillips N A.A coordinate system having some special advantages for numerical forecasting［J］.Journal of Meteorology，1957，14：184－185.

［14］ Spalart PR.Strategies for turbulencemodelling and simulations［J］.International Journal of Heat Fluid Flow，2000，21：252 －263.

［15］ Shur M，Spalart P R，Strelets M，Travin A.Detached?eddy simulation of an airfoil at high angle of attack［C］∥Fourth International Symposium on Engineering Turbulence Modelling and Experiments，Rodi W，Laurence D，eds.Corsica：Elsevier，New York，24－26 May，1999.

［16］ Spalart PR，Rumsey C L.Effective inflow conditions for turbulencemodels in aerodynamic calculations［J］.AIAA Journal，2007，45（10）：2544－2553.

［17］ Aupoix B，Spalart PR.Extensions of the Spalart?Allmaras turbulencemodel to account for wall roughness［J］.International Journal of Heat and Fluid Flow，2003，24：454－462.

［18］ Eca L，Hoekstra M，Hay A，Pelletier D.A manufactured solution for a two?dimensional steady wall?bounded incompressible turbulent flow［J］.International Journal of Computational Fluid Dynamics，2007，21（3/4）：175－188.

［19］ YueW，Lin C L，Patel V C.Large eddy simulation of turbulent open?channel flow with free surface simulated by level set method［J］.Physics of Fluids，2005，17：1－12.

［20］ Keating A，Piomelli U.A dynamic stochastic forcingmethod as a wall?layermodel for large?eddy simulation［J］.Journal of Turbulence，2006，7（2）：1－24.

［21］ Breuer M，Jovicˇic＇N，Mazaev K.Comparison of DES，RANS and LES for the separated flow around a flat plate at high incidence［J］.International Journal for Numerical Method in Fluids，2003，41：357－388.

［22］ Deck S.Numerical simulation of transonic buffet over a supercritical airfoil［J］.AIAA J，2005，43（7）：1556－1566.

［23］ Deck S.Zonal?detached eddy simulation of the flow around a high?lift configuration［J］.AIAA J，2005，43（11）：2372－2384.

［24］ Deck S，Weiss P，Pamiès M，Garnier E.Zonal detached eddy simulation of a spatially developing plate turbulent boundary layer［J］.Computers＆Fluids，2011，48：1－15.

［25］ Casulli V，Cattani E.Stability，accuracy and efficiency of a semi?implicitmethod for three?dimensional shallow water flow ［J］.Computers and Mathematics with Applications，1994，27（4）：99－112.

［26］ Darwish MS，Moukalled F.TVD schemes for unstructured grids［J］.International Journal of Heat and Mass Transfer，2003，46：599－611.

［27］ Balachandar R，Polatel C，Hyun B?S，Yu K，Lin C?L，YueW，Patel V C.LDV，PIV，and LES investigation of flow over a fixed dune［C］∥Proc，Symp Held in Monte Verta：Sedientation and Sediment Transport，Monte Verita，Switzerland，Kluwer Academic，Dordrecht，2002：171－178.

［28］ Hinterberger C，F(xiàn)rohlich J，RodiW.Three?dimensional and depth?averaged large?eddy simulations of some shallow water flows ［J］.Journal of Hydraulic Engineering，2007，133（8）：857－872.

A Detached Eddy Simulation Model for Free Surface Flows w ith Uneven Bottom

ZHANG Jingxin1，2
（1.MOE Key Laboratory ofHydrodynamics，Shanghai Jiao Tong University，Shanghai 200240，China；
2.Key Laboratory of Estuarine＆Coastal Engineering，Ministry of Transport，Shanghai 201201，China）

A detached?eddy simulation（DES）model is proposed based on a fully hydrodynamic pressuremodel instead of hydrostatic model.The numerical scheme is based on finite volumemethod（FVM）on unstructured grids in the horizontal plane，andσcoordinate in vertical direction to fix free surface and uneven bottom.The in?house codes are paralleled using OpenMP.The proposed model is shown particularly effective in prediction of small?scale vortical structures.

detached eddy simulation；coherent structure；free surface flow；hydrodynamic model

O352

A

2014－12－11；

2015－01－13

國家重點(diǎn)基礎(chǔ)研究發(fā)展計(jì)劃（973）（2014CB046200）及水利部公益性行業(yè)科研專項(xiàng)經(jīng)費(fèi)（201401027）資助項(xiàng)目

張景新（1975－），男，副教授，從事環(huán)境流體力學(xué)研究，E?mail：zhangjingxin＠sjtu.edu.cn

Received date： 2014－12－11；Revised date： 2015－01－13