考慮波浪形底面影響的邊界層風(fēng)場大渦模擬

孫麗明，曹曙陽，李 明，楊 青，張恩臻

（同濟(jì)大學(xué) 橋梁工程系土木工程防災(zāi)減災(zāi)國家重點(diǎn)實(shí)驗(yàn)室，上海 200092）

考慮波浪形底面影響的邊界層風(fēng)場大渦模擬

孫麗明，曹曙陽，李 明，楊 青，張恩臻

（同濟(jì)大學(xué) 橋梁工程系土木工程防災(zāi)減災(zāi)國家重點(diǎn)實(shí)驗(yàn)室，上海 200092）

為了研究強(qiáng)風(fēng)作用下海面上的平均風(fēng)剖面，采用大渦模擬方法對帶有余弦波形狀底面、自由滑移頂面的渠道流進(jìn)行了模擬，波幅、波高比2a／λ為0.1，基于平均速度Ub和渠道高度的雷諾數(shù)Reb為6760。統(tǒng)計(jì)結(jié)果表明，波浪形下表面對上方風(fēng)場分布有著顯著的影響，y／H=0.3（y+≈200）是內(nèi)區(qū)、外區(qū)的分界線，內(nèi)區(qū)受壁面影響顯著，外區(qū)受壁面影響較小。流動的分離點(diǎn)位于x／λ=0.14，而再附點(diǎn)位于x／λ=0.65。展向速度脈動峰值出現(xiàn)在上坡處，且超過豎向脈動。表面壓力要比表面摩擦力大一個數(shù)量級，是阻力的主要來源。隨著波幅的增大，回流區(qū)面積增大，速度峰值和展向脈動峰值也會增大，展向脈動峰值甚至?xí)^流向脈動峰值。

大渦模擬；波浪形渠道流；平均風(fēng)剖面；渠道流；摩擦阻力；正弦波

0 引 言

自然界中的流體存在著各種各樣的流動形態(tài)，下墊面的形狀以及表面的粗糙度分別會以形狀阻力（form drag）和表面摩擦阻力（skin friction）的形式對上層的流體流動造成影響，根據(jù)下墊面的形狀以及粗糙度的不同，兩者對流動影響的貢獻(xiàn)也不同。而流體的流動有很大一部分是發(fā)生在下墊面不是水平光滑的情況下的，例如，丘陵地帶的風(fēng)場正是氣體在連續(xù)的小山包上方的流動，海水對海底沙石的沖刷和搬運(yùn)也是發(fā)生在高低起伏的海床上方的，更復(fù)雜的情況還有海面上的風(fēng)場分布，在海氣接觸面上，由于風(fēng)以及重力等的作用，海水無法再保持水平光滑，海面波濤起伏，海氣作用相互耦合，使得海上的風(fēng)場分布發(fā)生改變，從而使得按照傳統(tǒng)規(guī)范給出的風(fēng)剖面設(shè)計(jì)的海上構(gòu)筑物如橋梁、海上鉆井平臺、海上風(fēng)力發(fā)電機(jī)等面臨失效的風(fēng)險，甚至導(dǎo)致重大的財產(chǎn)損失和人員傷亡。因此，研究海上風(fēng)場分布有著重要的現(xiàn)實(shí)意義，而研究海上風(fēng)場分布就需要建立兩項(xiàng)流海氣相互作用模型，作為基礎(chǔ)研究，靜止的曲面渠道流中正弦波（或余弦波）形狀可以近似模擬波浪的形狀，對上層氣體流動形成連續(xù)的阻礙，研究這種流場分布可以為更深入的研究奠定基礎(chǔ)。

風(fēng)洞實(shí)驗(yàn)作為研究流動的重要手段之一，一直受到廣大研究者的重視。關(guān)于波浪形渠道流的研究，Hanratty和他的同事們做了大量的風(fēng)洞實(shí)驗(yàn)，例如，Zilker，Cook ＆Hanratty 1977［1］；Zilker＆ Hanratty 1979［2］；Buckles，Hanratty＆ Adrian 1984［3］，1986［4］（簡稱 BHA）；Frederick ＆Hanratty 1988［5］；Kuzan，Hanratty＆Adrian 1989［6］；Hudson，Dykhno＆ Hanratty 1996［7］，Gong et al.1996［8］。其中 BHA 實(shí)驗(yàn)給出了大波幅波浪渠道流平均風(fēng)速和脈動風(fēng)速的大量測量數(shù)據(jù)，而本文中考慮大波幅和中等波幅兩種情況。

流動的理論分析（Sykes 1980［9］；Hunt，Leibovich＆ Richards 1988［10］；Belcher，Newley ＆ Hunt 1993［11］）表明，地形引起的湍流結(jié)構(gòu)改變在壁面附近變異性非常大，傳統(tǒng)對數(shù)率只在離壁面非常近的范圍內(nèi)成立，粘性底層的尺度要比波浪形渠道的波長小很多，通常為波長的5%，使直接捕捉這部分渦尺度變得非常困難，從而影響求解表面摩擦力的準(zhǔn)確性，因?yàn)楸砻婺Σ亮εc表面流場有著密切的聯(lián)系。

近年來，隨著計(jì)算機(jī)技術(shù)的發(fā)展，數(shù)值模擬的方法被廣泛應(yīng)用于流體計(jì)算。波浪形渠道流的 DNS （Maa?＆Schumann1994［12］，1996［13］；De Angelis，Lombardi ＆ Banerjee 1997［14］；Cherukat et al.1998［15］）成功的再現(xiàn)了中等雷諾數(shù)下實(shí)驗(yàn)室 觀測到的現(xiàn)象，并且提供了更加詳細(xì)的湍流結(jié)構(gòu)數(shù)據(jù)。其中De Angelis et al.和Cherukat et al.均發(fā)現(xiàn)了上坡處展向速度脈動增大；De Angelis et al.采用了偽譜方法，他們在波谷的下游處（上坡處）發(fā)現(xiàn)了高剪應(yīng)力區(qū)和條帶結(jié)構(gòu)的存在，同時指出準(zhǔn)流線渦也在此處產(chǎn)生。大渦模擬方法也越來越多的被應(yīng)用于各種復(fù)雜幾何外形的流動模擬。最早利用LES技術(shù)對非水平均勻渠道流進(jìn)行研究的有Krettenauer＆Schumann （1992）［16］，Walko，Cotton＆ Pielke（1992）［17］，以及D?rnbrack＆Schumann（1993）［18］，他們對正弦波上方的湍流對流流動進(jìn)行了模擬，對波浪形渠道流進(jìn)行了初步的研究。Henn and Sykes（1999）［19］采用了LES研究了不同波幅的波浪形渠道流，他們也發(fā)現(xiàn)了上坡處展向速度脈動增大，并指出在波峰背風(fēng)面存在的分離剪切層導(dǎo)致了波谷上方產(chǎn)生強(qiáng)湍流，他們在模擬中采用了湍動能輸運(yùn)方程、二階閉合模型（Lewellen 1977［20］），并利用坐標(biāo)轉(zhuǎn)換將曲線坐標(biāo)系，轉(zhuǎn)換成直角坐標(biāo)系，在此過程中，增加了數(shù)值離散的難度，加大了計(jì)算量。Calhoun and Street（2001）［21］同樣采用了LES方法，不同的是他們采用的動態(tài)混合亞格子模型（dynamic mixed subgrid-scale model，Zang at al.1993［22］and Salvetti et al.1997［23］）。但是以上研究都是基于頂部是無滑移條件的，而在風(fēng)工程應(yīng)用中，波浪上方僅僅是大氣，沒有頂部壁面，因此該邊界條件不能適用于風(fēng)工程。

本文采用大渦模擬方法，非正交網(wǎng)格，對散度離散進(jìn)行修正，亞格子模型采用動態(tài)的Smagorinsky模型（Lilly，1992［24］），考慮到達(dá)到一定高度后，上層的流動趨于層流化，因而頂部采用自由滑移條件。

本文安排如下：首先介紹數(shù)值模型與方法，包括各個算例的網(wǎng)格與計(jì)算參數(shù)等；作為對數(shù)值模型的驗(yàn)證，接著對直渠道流進(jìn)行了計(jì)算；然后詳細(xì)討論了中等波幅（2a／λ=0.1）波浪形渠道流的計(jì)算結(jié)果，包括統(tǒng)計(jì)后的平均流場和瞬態(tài)流場；最后給出了大波幅（2a／λ= 0.2）波浪形渠道流的部分計(jì)算結(jié)果，并作了對比。

1 數(shù)值模型與方法

1.1 LES數(shù)值模型

采用直角坐標(biāo)系，x、y、z分別為流向、豎向和展向，如圖1所示。底面的高程可以表示成x、y的函數(shù)（實(shí)際只是x的函數(shù)），即h（x，y）=a cos（2πx／λ），其中a為波幅，λ為波長。過濾后的速度矢量為u=（u，v，w），直角坐標(biāo)為x=（x，y，z）。過濾后的Navier-Stokes方程可以表示為：

圖1 波浪形渠道流示意圖Fig.1 Sketch of channel with a periodic wavy wall

連續(xù)性方程為：

其中，動量方程（1）中的6項(xiàng)分別為瞬態(tài)項(xiàng)、對流項(xiàng)、壓力源項(xiàng)、耗散項(xiàng)、亞格子應(yīng)力項(xiàng)和驅(qū)動力項(xiàng)，P= （Px，0，0），Px是x 方向的驅(qū) 動力，計(jì)算過程中根 據(jù)平均速度Ub來調(diào)整Px，以保證Ub=1為定值。采用Boussinesq渦粘假設(shè) （Boussinesq eddy viscosity assumption），令

因此，在程序中定義

方程 （1）可以寫成

其中νEff稱為有效粘性系數(shù)，νSGS為亞格 子湍動粘 度，本文選用動態(tài)的Smagorinsky亞格子模型。模型采用有限體積法，對流項(xiàng)和擴(kuò)散項(xiàng)采用二階中心差分，對梯度項(xiàng)進(jìn)行網(wǎng)格非正交修正，時間導(dǎo)數(shù)項(xiàng)采用二階向后差分，網(wǎng)格采用結(jié)構(gòu)化網(wǎng)格，壓力速度耦合采用PISO算法。

1.2 網(wǎng)格與計(jì)算參數(shù)

本文首先計(jì)算了兩個直渠道流的算例，分別為FC1和FC2（FC=Flat Channel），渠道尺寸均為（L，W，H）=（4π，2π，2），見圖1，對直渠道流，波長λ=∞，基于平均流速和半槽寬的雷諾數(shù)Reb=0.5HUb／ν= 3250，其中，平均流速

u是流向速度。FC1流向、豎向和展向網(wǎng)格數(shù)均為64，流向和展向網(wǎng)格均勻分布，豎向網(wǎng)格按余弦函數(shù)分布，即

FC2流向、豎向和展向網(wǎng)格數(shù)均為96，流向和展向網(wǎng)格均勻分布，豎向網(wǎng)格按雙曲正切函數(shù)分布。直渠道流的邊界條件設(shè)置為：下底面和上頂面均為無滑移條件，四個側(cè)面均為周期條件。亞格子模型采用動態(tài)的Smagorinsky模型，空間平均采用沿流向和展向平均，時間平均取200L／Ub。

波浪形渠道流兩個算例分別為WC1、WC2（WC =Wavy Channel），WC1的計(jì)算域尺寸為（L，W，H）=（4，2，1），如圖1所示，H=1為平均槽寬，波高a= 0.05，波長λ=1，波陡為ka=2πa／λ≈0.314，與 Maa? ＆Schumann（1996）DNS相同，豎向最小、最大無量綱網(wǎng)格高度分別為0.001和0.05，無量綱化采用如下公式

其中，h為底面高程，h（x，y）=a cos（2πx／λ），豎向網(wǎng)格的伸展率約為1.05。基于平均流速和平均槽寬的雷諾數(shù)Reb=HUb／ν=6760。WC2的計(jì)算域尺寸為（L，W，H′）=（4.286，1.072，1.072），如圖1，平均槽寬H= H′-a=1.072-0.01072=0.9648，H′為最大槽寬，波高a=0.1072，波長λ=1.072，波陡為ka=2πa／λ≈0.628，與Henn and Sykes（1999）LES的BHA1L算例及實(shí)驗(yàn)BHA相同，豎向最小、最大無量綱網(wǎng)格高度分別為0.001和0.05，Reb=H′Ub／ν=6760。兩個算例的網(wǎng)格數(shù)均為（Nx，Ny，Nz）=（192，96，96），網(wǎng)格示意圖如圖2所示。圖2中，由于網(wǎng)格太密，網(wǎng)格線采用每四根網(wǎng)格線一顯示。兩個算例底面均為四個波長的余弦波，如圖1所示，流向長度為L，展向長度和波長分別為W 和λ。四個側(cè)面均采用周期條件，底面采用無滑移條件，而為了風(fēng)工程研究的需要，頂面應(yīng)該是與大氣接觸的，不存在固定壁面，所以應(yīng)采用自由滑移條件，這是與以往文獻(xiàn)中的邊界條件均不同的（Nakayama 2002除外），頂面邊界條件的不同使得平均風(fēng)剖面的結(jié)果在0.5H 以上的部分與以往的數(shù)值模擬和風(fēng)洞模擬有所不同。

計(jì)算結(jié)果的平均先對時間進(jìn)行，平均采用的時間長度為100H／Ub，此處的平均流速定義為這里均控制Ub=1，再對所有展向取平均，最后對四個周期相同相位的數(shù)值取平均。

圖2 波浪形渠道流網(wǎng)格示意圖Fig.2 Grids of the wavy channel（Every four grid lines shown）

2 結(jié)果分析

2.1 直流渠道

直渠道流的研究比較廣泛，結(jié)果也比較成熟，有大量的DNS數(shù)據(jù)和實(shí)驗(yàn)數(shù)據(jù)可供參照對比，是檢驗(yàn)數(shù)值算法的很好算例。為了驗(yàn)證數(shù)值方法的可靠性，本文首先采用相同的LES模型（包括離散方法，湍流模型）對直渠道流進(jìn)行了模擬。無量綱的流向平均速度剖面如圖3（a）所示，與KMM （Kim，Moin and Moser 1987［25］）DNS對比，顯示了很好的符合性，在粘性地層y+＜10的范圍內(nèi)，LES結(jié)果與 DNS結(jié)果差別很小，基本滿足u+=y+的壁面率；在y+＞20，LES的速度略大于KMM，與Eckelmann（1974）風(fēng)洞實(shí)驗(yàn)數(shù)據(jù)吻合得很好，F(xiàn)C2在對數(shù)層的速度值略低于FC1，差別不大，說明了網(wǎng)格的精度足夠，數(shù)值計(jì)算結(jié)果是收斂的。其中，y+=yuτ／ν，u+=u／uτ。

圖3（b、c）顯示了速度脈動沿槽道高度方向的分布，速度脈動均采用壁面剪切速度uτ無量綱化，與KMM 相比，剖面的整體形狀符合的比較好，兩個算例的最大流向速度脈動均出現(xiàn)在y+=12處，與產(chǎn)生最大湍動能的位置相符。三個方向脈動速度的峰值位置均與KMM 相同，在壁面的地方吻合較好，但是在y+≈5以外 LES算得的流向速度脈動值要略大于KMM的DNS值，F(xiàn)C2與Eckelmann實(shí)驗(yàn)值較吻合，而展向和豎向速度脈動則要略小于DNS值，在壁面附近數(shù)值模擬的值都要明顯小于實(shí)驗(yàn)值，實(shí)驗(yàn)測得的豎向、展向脈動在壁面附近增長得更快，展向脈動速度的峰值更加明顯。圖4雷諾應(yīng)力也顯示了很好的吻合性。其中，

圖3 近壁速度分布Fig.3 Velocity distribution near the wall

圖4 雷諾應(yīng)力Fig.4 Reynolds stress

2.2 中等波幅波浪形渠道流

2.2.1 平均流場

由于波浪形下表面沿展向是均勻不變的，沿流向的四個波長的對應(yīng)位置又呈現(xiàn)出周期性變化的特點(diǎn)，所以統(tǒng)計(jì)特性在展向和同相位處是統(tǒng)計(jì)均勻的，故所有的統(tǒng)計(jì)量取均值的方法都是首先對時間平均，再沿展向平均，最后對同相位平均。所有的圖形都是統(tǒng)計(jì)后的量在豎直平面內(nèi)繪出的。

根據(jù)計(jì)算結(jié)果，波浪形下表面對上層流體流動平均特性的影響在0.3倍渠道高度以下（稱為內(nèi)區(qū)，inner region）影響明顯，此高度往上（稱為外區(qū)，outer region），流體受到的影響明顯減弱，平均速度剖面趨于一致，如圖6、圖7 所示（y／H=0.3對應(yīng)y+≈200），故云圖只繪出了內(nèi)區(qū)，即0.3H 以下部分。流向平均速度云圖如圖5所示，從圖中可以看出，分離點(diǎn)（separation point）在x／λ=0.14，再附點(diǎn)（reattachment point）位于x／λ=0.65，在分離點(diǎn)與再附點(diǎn)之間存在著一個回流區(qū)，形狀呈扁長的卵狀，回流區(qū)中心位于比波谷略靠上游的位置，相應(yīng)的無量綱坐標(biāo)為（0.48，-0.037），對應(yīng)的平均速度最小值約為-0.08。整個回流區(qū)的長度約為0.5λ，回流區(qū)覆蓋波谷上游（波峰的背風(fēng)面，或稱下坡面）0.36λ，覆蓋波谷下游（上坡面）0.15λ，最大高度約為0.04H?；亓鲄^(qū)的形狀由平均流線圖6也可以清楚地看出。波谷上方速度的豎向梯度較小，變化較緩慢，而在波峰處以及鄰近的兩側(cè)等值線密集，速度的豎向梯度很大，尤其是0～0.1λ和0.7λ～1.0λ范圍內(nèi)，0.1H 往上部分，速度變化勻緩。將0.7λ～1.1λ區(qū)間貼近壁面的部分稱為邊界層區(qū)（boundary layer region），邊界層區(qū)的顯著特點(diǎn)是速度的展向變化大、展向脈動大，這在后面的瞬態(tài)速度向量圖和展向脈動圖中可以看出。

圖5 平均流向速度云圖Fig.5 Mean streamwise velocity contour

圖6 平均流線圖Fig.6 Mean streamlines

流向平均速度剖面如圖7所示，圖中藍(lán)色帶星號的是當(dāng)前的LES結(jié)果，黑色帶圓圈的是Maa?等人的DNS 結(jié)果，粉色實(shí)線是 Nakayama［26］的 LES 結(jié)果，由圖可見速度剖面在0.5倍高度以下吻合得很好，由于 Maa?頂部采用的是無滑移條件，頂部速度等于0，而本文為模擬大氣邊界層，頂部采用的是自由滑移條件，頂部速度最大，約為1.19，Maa? 的DNS中，由于上下部都受到壁面摩擦阻力作用，剖面在渠道中部更趨于外凸。圖中每隔約0.1λ繪制一個剖面，為避免相鄰的剖面交錯在一起，圖中速度的值縮小了12倍，可以看出剖面的變化趨勢與云圖的結(jié)論一致，波峰處的剖面形狀更接近于直渠道流，而在波峰背風(fēng)面的分離點(diǎn)之后，由于負(fù)壓力梯度的影響，速度剖面開始出現(xiàn)反彎，在接近壁面處的速度為負(fù)值，一直到再附點(diǎn)，壁面附近的負(fù)速度消失，即回流區(qū)部分。同樣，剖面形狀也顯示出波峰處由于迎風(fēng)截面受到壓縮，速度加速，底部速度梯度很大。波谷處的速度剖面變化則相對平緩。

圖7 平均流向速度剖面Fig.7 Mean streamwise velocity over the first wavelength

無量綱化以后的流向速度剖面如圖8所示，每隔0.1λ繪制一個剖面，為便于對比，相隔0.5λ的剖面繪在一起，圖中帶圓圈的是前半個波長上的剖面，帶三角的是后半個波長上的剖面，藍(lán)色實(shí)線是平均剖面，無量綱化時采用沿下表面平均后的uτ。Maa?的DNS只繪出了0.5H 以下的部分，可以看出渠道中部，其速度要比本文偏大，這與上邊界條件設(shè)置是一致的。相隔半個波長的兩個剖面顯示出很強(qiáng)的規(guī)律性，即他們基本關(guān)于平均速度剖面對稱，分布在其上下兩側(cè)。在x=0.1、0.7、0.8、0.9、1.0處的5個剖面落在回流區(qū)之外，速度均不出現(xiàn)負(fù)值；而其余5個剖面均落在回流區(qū)，近壁速度先負(fù)后正。外區(qū)（y＞0.3，y+＞200）部分，速度剖面基本重合，受下表面的影響可以忽略。由外區(qū)的直線關(guān)系可以看出，外區(qū)的速度剖面還是滿足對數(shù)率的，只是此時由于采用的不是當(dāng)?shù)啬Σ了俣?，而是沿整個下表面平均后的摩擦速度，對數(shù)率的系數(shù)發(fā)生了變化。

平均豎向速度云圖如圖9所示，可以分為三個區(qū)域，正速度區(qū)被負(fù)速度區(qū)分隔成了兩塊，第一塊位于負(fù)速度區(qū)底部的回流區(qū)內(nèi)，區(qū)間范圍為0.14λ～0.5λ，起點(diǎn)與分離點(diǎn)重合，終點(diǎn)位于波谷；第二塊位于邊界層區(qū)內(nèi)，范圍為0.65λ～1.0λ，起點(diǎn)與再附點(diǎn)重合，由于受到上坡的擠壓，平均豎向 速度變正，豎向速度迅速的從0增加到峰值0.108，位置位于波峰迎風(fēng)側(cè)，無量綱位置坐標(biāo)為（0.86，0.05），高度與波峰齊平，之后又向波峰處迅速減小，在峰頂附近減為0；然后豎向速度沿下坡緩慢的變?yōu)橄蛳碌乃俣龋俣冉^對值逐漸增大，在波谷上方的下游出現(xiàn)最值-0.042，無量綱坐標(biāo)為（0.53，0.05），高度與波峰高程齊平，之后又逐漸由負(fù)變正，負(fù)值區(qū)的范圍為0～0.65λ，負(fù)速度區(qū)的覆蓋范圍要大于正速度區(qū)，正負(fù)最值的比約為2.5。

圖8 無量綱平均流向速度剖面Fig.8 Normalized mean streamwise velocity profiles

圖9 平均豎向速度云圖Fig.9 Mean vertical velocity contour

從平均豎向速度剖面圖10中，也可以看出速度的3個區(qū)域，在回流區(qū)底部豎向速度有正值，回流區(qū)上方速度變?yōu)樨?fù)值。豎向速度在沖擊帶達(dá)到最大值，沖擊帶的概念見下文展向速度脈動部分。

圖10 平均豎向速度剖面圖Fig.10 Mean vertical velocity profiles

相比平均流向速度和豎向速度，平均展向速度盡管在波谷上方以及邊界層區(qū)也會出現(xiàn)峰值，平均速度剖面的波動性很大，其他規(guī)律不是那么明顯。

展向速度沿展向的變化較大，呈現(xiàn)正負(fù)交替的現(xiàn)象，故沿展向作平均后的結(jié)果會使得一些信息喪失，從而研究展向速度的平均值意義不大，本文將在下文的瞬態(tài)速度場中研究展向速度的分布規(guī)律。

與平均風(fēng)速一樣，平均后的脈動特性也顯示了很強(qiáng)的規(guī)律性。圖11（a）為平均流向脈動速度平方的云圖，流向脈動速度平方在波谷上方出現(xiàn)一個最大值區(qū)域，起始點(diǎn)位置基本與分離點(diǎn)相同，無量綱坐標(biāo)為（0.145，0.068），終點(diǎn)為（0.82，0.039），峰值為0.0545，峰值點(diǎn)位于波谷上方，約0.8倍波幅高度處，坐標(biāo)為（0.48，0.041）。脈動值在壁面附近迅速減小，尤其是邊界層區(qū)0.65λ～1.1λ范圍內(nèi)，脈動速度梯度很大，回流區(qū)梯度相對較小，壁面處減為0。從剖面角度看，流向各個位置處，出現(xiàn)峰值的高度大約都在0.1H 以下，邊界層區(qū)的峰值位置高度非常接近壁面，小于0.01H，波峰處的峰值高度約為0.01H，波谷處約為0.04H。

圖11 平均脈動速度平方云圖Fig.11 Mean velocity fluctuation

與流向、豎向脈動不同，平均展向速度脈動的最大值出現(xiàn)在上坡處，如圖11（c）所示，最大值為0.0272，是流向最大值的0.5倍，是豎向最大值的1.8倍。最大 值 出現(xiàn) 的位 置坐 標(biāo)為 （0.73，0.013），離壁面非常近?？梢孕蜗蟮貙⒃俑近c(diǎn)到波峰之前的這段上坡區(qū)域叫做沖擊帶 （impact zone），當(dāng)波幅增大，2a／λ=0.2時，展向脈動甚至超過了流向脈動，這在Henn and Skyes，De Angelis，Cherukat等人的研究中均有體現(xiàn)。這意味著在沖擊帶存在著與波陡有關(guān)的局部能量產(chǎn)生機(jī)理。

平均展向渦量ωz云圖如圖12所示，在0.8H 以下渦量明顯，壁面附近渦量出現(xiàn)正負(fù)峰值，0.8H 以上，渦量接近于0。壁面附近，渦量以回流區(qū)為界，分成正負(fù)渦量區(qū)，回流區(qū)以內(nèi)，渦量為正值，即渦逆時針旋轉(zhuǎn)，無滑移的壁面壁面對流體施加的剪力作用沿壁面切線方向向左，最大值約為21，位于波谷處。渦量負(fù)值區(qū)則分布在回流區(qū)兩側(cè)，并且在分離點(diǎn)處向上抬升，拖出一條很寬的尾巴，一直蔓延到波谷上方，最小值則以0.9λ為中心，相應(yīng)的最小值約為-88，此處位于沖擊帶，負(fù)渦量的絕對值約為正渦量的4倍。

圖12 平均展向渦量ωz云圖Fig.12 Mean spanwise vortictyωz

統(tǒng)計(jì)量在豎直面內(nèi)的分布與下表面的幾何形狀（波長／波幅）有著密切的聯(lián)系，幾何形狀在很大程度上決定了回流區(qū)的形狀與位置，下表面相對波幅減小甚至?xí)?dǎo)致回流區(qū)消失，極限情況就是直渠道流（2a／λ=0），Henn and Skyes的研究中，2a／λ=0.031時就未出現(xiàn)回流區(qū)。

從脈動速度的云圖可以看出，流向速度脈動是占主導(dǎo)地位的，所以湍動能的分布更接近于流向速度脈動，如圖13所示，峰值出現(xiàn)在波谷上方，高度比波峰高度略高，最大值為0.029。

表面的摩阻力以及表面壓力也是重要的空氣動力學(xué)特性指標(biāo)之一，許多領(lǐng)域都很關(guān)心與流體接觸面的受力，海洋學(xué)中研究風(fēng)浪相互作用時海面的拖拽系數(shù)，再如風(fēng)工程領(lǐng)域關(guān)心的大氣邊界層粗糙度類別，實(shí)際就是表面摩阻力的工程化算法。故本文對波浪形下表面的表面摩阻力τw和表面壓力pw作了研究，根據(jù)牛頓液體的應(yīng)力應(yīng)變關(guān)系τ0=μδˉu／δy，其中δˉu／δy是第一層網(wǎng)格處的平均速度剪切率，這里采用τw=τ0／（0.5ρU2b）進(jìn) 行 無量 綱 化，Ub=1為 平 均 來 流速度。τw沿流向的分布如圖14（a）所示，在x／λ≈0.14～0.66之間，τw為負(fù)值，此區(qū)間與回流區(qū)區(qū)間簡本重合，再次說明了回流區(qū)流體是沿順時針方向運(yùn)動的，流體対壁面的作用力是向左的，從分離點(diǎn)開始，負(fù)的摩阻力逐漸增大，在0.2λ～0.35λ之間增長緩慢，幾乎不變，負(fù)的最大摩阻力出現(xiàn)在波谷x／λ=0.5處，值為-0.0047，之后一直減小，大約在再附點(diǎn)處減到0。進(jìn)入沖擊帶之后，摩阻力變正，由于此處速度梯度大，摩阻力也大，且增長迅速，最大值0.0223出現(xiàn)在x／λ=0.917處，之后一直減小到分離點(diǎn)的0值。

圖13 平均湍動能k云圖Fig.13 Mean kinetic energy k

圖14 壁面平均壓力平均摩阻力τw與pw沿流向的分布Fig.14 Mean surface drag and pressure distribution along x-direction

壁面平均壓力pw的展向分布如圖14（b）所示，這里pw=（p-pref）／0.5ρU2b，正值 區(qū) 間 為 0.419λ～0.872λ，其余為負(fù)值區(qū)間，最大正壓 0.154，位于0.71λ處，最大負(fù)壓-0.14，位于波峰處。在靠近分離點(diǎn)處，pw曲線的斜率明顯降低。

pw與形狀阻力對應(yīng)，τw與表面摩擦阻力對應(yīng)，由圖可見，形狀阻力要比表面摩擦阻力大一個數(shù)量級。2.2.2 瞬態(tài)流場

這里采用λ2識別法來識別渦，其識別標(biāo)準(zhǔn)是S2+Ω2的特征值有兩個為負(fù)值。圖15是采用λ2識別法來識別近壁渦結(jié)構(gòu)的可視化圖形，圖中顯示了λ2=-20的等值面。由圖可知，流向渦占主導(dǎo)地位，是主要的渦結(jié)構(gòu)形式，這些流向渦或多或少與x軸存在一定的交角，故可稱之為準(zhǔn)流向渦，他們絕大多數(shù)起始于上坡面（沖擊帶），跨過波峰向下游的波谷延伸。圖中也可見到展向渦，展向渦的尺度相對較小。有些地方甚至可以見到發(fā)卡渦，這種擬序結(jié)構(gòu)（coherent structure）在直渠道流中很常見。

圖15 λ2等值面圖Fig.15 Iso-surface ofλ2=-20

圖16為瞬態(tài)速度場垂直于流向的切片，切片位于x／λ=1.75，在沖擊帶上方，圖中顯示的高度為0.3H 以下。圖中展向的運(yùn)動速度明顯，流場中出現(xiàn)一系列類似渦形態(tài)的展向運(yùn)動結(jié)構(gòu)，貼近壁面處的展向“類渦”運(yùn)動最為明顯。這種流動可以形象地稱作次生流動，主要出現(xiàn)在沖擊帶上方。由于壁面的阻擋和反射作用，使得近壁面處的豎向脈動衰減，而展向沒有壁面阻擋，這部分湍動能會注入到展向脈動中，即展向脈動會相應(yīng)的增加，這也可能是壁面附近展向速度|w′|要比豎向速度|v′|大的原因。

圖16 瞬態(tài)速度場x／λ=1.75，z=0.5λ～1.5λFig.16 Slice of instantaneous velocity vectors at x／λ=1.75

2.3 大波幅波浪形渠道流

當(dāng)波幅增大后，流場的流動形態(tài)會發(fā)生變化，為了對比，以下給出了 WC2（2a／λ=0.2）的部分流場情況。

由圖17（a）和圖17（b）可以看出，由于波幅的增大，回流區(qū)更加飽滿，覆蓋了兩個波峰之間的絕大部分，分離點(diǎn)由0.14前移到0.079，再附點(diǎn)從0.65延后到0.766；波幅增大后，正負(fù)峰值速度的絕對值均顯著增大，最大負(fù)速度由-0.08變?yōu)?0.218，頂部最大速度也由1.19增大到1.41；流線圖中，旋渦更加豐滿，渦頂部幾乎保持水平，且位置與波峰齊平。

圖17 大波幅平均統(tǒng)計(jì)結(jié)果Fig.17 Averaged results of the large amplitude wavy channel

平均豎向速度依然分為三個區(qū)域，如圖17（c）所示，最大正負(fù)值速度分別由-0.042和0.108變?yōu)?0.09和0.16，底部的正速度區(qū)面積增大，從變長的卵形變?yōu)榻频闹苯侨切?，相?yīng)的負(fù)值區(qū)域有所減小；正峰值位置變化不大，負(fù)峰值后移，且更貼近壁面。圖17（d）平均展向脈動的峰值仍然位于沖擊帶上方，位置略微后移，峰值由0.0272變?yōu)?.072，且超過了流向脈動的峰值0.065和豎向脈動峰值0.028，為三個方向脈動最大的方向，壁面對三個方向脈動能量分配的影響更顯著。

3 結(jié)論與展望

本文采用LES方法對帶有余弦波形狀底面、自由滑移頂面的渠道流進(jìn)行了模擬，得到的結(jié)論總結(jié)如下：

（1）統(tǒng)計(jì)結(jié)果表明波浪形下表面通過其波幅和波長影響著流場在空間上的分布形態(tài)，根據(jù)受壁面形狀影響的大小，在豎向可將流場分成內(nèi)區(qū)和外區(qū)，他們以y／H=0.3（y+≈200）為界，內(nèi)區(qū)受壁面形狀影響很大；外區(qū)基本不受壁面形狀的影響，外區(qū)在流向不同位置的平均速度剖面趨于一致。

（2）分離點(diǎn)位于x／λ=0.14，再附點(diǎn)位于x／λ= 0.65，分離點(diǎn)與再附點(diǎn)之間為回流區(qū)，長度約為半個波長，位置偏上游；流向平均速度剖面在回流區(qū)會有反彎，回流區(qū)的平均豎向速度基本為正，回流區(qū)上方的平均豎向速度為負(fù)。

（3）再附點(diǎn)之后的上坡面到波峰偏下游為沖擊帶，此處存在著一些特殊的物理現(xiàn)象：流向速度在壁面處的梯度很大，過流面積的壓縮使流體加速；沖擊帶上方存在一個正的豎向速度區(qū)，且在（0.86，0.05）達(dá)到峰值；更顯著的特征是展向脈動速度在沖擊帶上方達(dá)到峰值，且超過豎向脈動，約為流向脈動的一半；負(fù)值渦量（渦沿順時針旋轉(zhuǎn)）在沖擊帶上方有最值，并且是最大正值渦量（波谷處）的4倍。

（4）在內(nèi)區(qū)，相隔半個波長的無量綱化平均速度剖面基本關(guān)于對整個場平均的平均速度剖面對稱，且不再符合傳統(tǒng)的x+=y+；在外區(qū)（y+＞200），平均速度剖面依然符合對數(shù)率；內(nèi)區(qū)的速度在波峰后減速、變負(fù)，然后再加速變正。（5）波谷上方是平均流向脈動豎向脈動以及湍動能k 的最值區(qū)（最大值）。

（6）表面摩阻力τw沿流向的分布規(guī)律也比較顯著，負(fù)摩阻力區(qū)間基本與回流區(qū)一致，最大正摩阻力位于沖擊帶，最大負(fù)摩阻力位于波谷處，且正負(fù)最值比約為5；表面壓力pw比表面摩阻力τw大一個數(shù)量級，對阻力的貢獻(xiàn)占主導(dǎo)地位。

（7）瞬態(tài)流場中，壁面附近是渦的主要存在區(qū)域，從λ2等值面圖中可以看到準(zhǔn)流向的擬序結(jié)構(gòu)；沖擊帶上方存在著類似渦形態(tài)的次生展向流動，近壁展向速度在展向?qū)诱?fù)交替分布，特征長度約為1／3λ。

（8）波幅增大后，流動分離更加明顯，回流區(qū)面積增大，形狀更加飽滿，正負(fù)峰值速度的絕對值均顯著增大；展向脈動顯著增加，并超過了流向脈動，壁面對三個方向脈動能量分配的影響更加顯著。

本文是為研究風(fēng)浪耦合作用下的海上風(fēng)剖面而做的基礎(chǔ)研究，進(jìn)一步深入的研究包括：

（1）移動下邊界波浪形渠道流的研究，考慮下表面作強(qiáng)制平移運(yùn)動或者波動對上層風(fēng)場的影響（后者為相位傳播，質(zhì)點(diǎn)不“隨波逐流”，更接近真實(shí)海浪的運(yùn)動），可以研究波齡對拖拽力的影響，然而強(qiáng)迫運(yùn)動也是一種理想化的模型，不能考慮風(fēng)浪的耦合作用，即風(fēng)-浪之間動量的互相傳遞，也不能考慮波浪的復(fù)雜形態(tài)，僅僅把浪當(dāng)作單一的簡諧波考慮。

（2）建立氣液兩相流模型，考慮風(fēng)浪的相互耦合，甚至考慮波浪的復(fù)雜形態(tài)，甚至考慮在強(qiáng)風(fēng)作用下波浪的破碎對拖拽力的影響。

［1］ ZILKER D P，COOK G W，HANRATTY T J.Influence of the amplitude of a solid wavy wall on a turbulent flow.Part 1.Nonseparated flows［J］.J.Fluid Mech.，1977，82：29-51.

［2］ ZILKER D P，HANRATTY T J.Influence of the amplitude of a solid wavy wall on a turbulent flow.Part 2.Separated flows ［J］.J.Fluid Mech.，1979，90：257-271.

［3］ BUCKLES J，HANRATTY T J，ADRIAN R J.Turbulent flow over large-amplitude wavy surfaces［J］.J.Fluid Mech.，1984，140：27-44.

［4］ BUCKLES J，HANRATTY T J，ADRIAN R J.Separated turbulent flow over a small amplitude wave.Laser Anemometry in Fluid Mechanics［M］｛II（ed.R.J.Adrian，D.F.G.Durao，F(xiàn).Durst，H.Mishina ＆J.H.Whitelaw），1986：347-357.Ladoan，Lisbon.

［5］ FREDERICK K A，HANRATTY T J.Velocity measurements for a turbulent nonseparated flow over solid waves［J］.Exps. Fluids，1988，6：477-486.

［6］ KUZAN J D，HANRATTY T J，ADRIAN R J.Turbulent flows with incipient separation over solid waves［J］.Exps.Fluids，1989，7：88-98.

［7］ HUDSON J D，DYKHNO L，HANRATTY T J.Turbulence production in flow over a wavy wall［J］.Exps.Fluids，1996，20：257-265.

［8］ GONG W，TAYLOR P A，D?RNBRACK A.Turbulent boundarylayer flow over fixed aerodynamically rough two-dimensional sinusoidal waves［J］.J.Fluid Mech.，1996，312：1-37.

［9］ SYKES R I.An asymptotic theory of incompressible turbulent boundary-layer flow over a small hump［J］.J.Fluid Mech.，1980，101：647-670.

［10］HUNT J C R，LEIBOVICH S，RICHARDS K J.Turbulent shear flow over low hills［J］.Q.J.R.Met.Soc.，1988，114：1435-1470.

［11］BELCHER S E，NEWLEY T M J，HUNT J C R.The drag on an undulating surface induced by the flow of a turbulent boundary layer［J］.J.Fluid Mech.，1993，249：557-596.

［12］MAA?C，SCHUMANN U.Numerical simulation of turbulent flow over a wavy boundary.Direct and Large Eddy Simulation-I［M］.（ed.P.R.Voke，L.Kleiser＆J.P.Chollet）.Kluwer，1994.

［13］MAA?C，SCHUMANN U.Direct numerical simulation of separated turbulent flow over a wavy boundary.Flow Simulation with Higher Performance Computers-II（ed.E.H.Hirschel）［M］.1996，227-241.

［14］DE ANGELIS V，LOMBARDI P，BANERJEE S.Direct numerical simulation of turbulent flow over a wavy wall［J］.Phys.Fluids，1997，9：2429-2442.

［15］CHERUKAT P，NA Y，HANRATTY T J，et al.Direct numerical simulation of a fully developed turbulent flow over a wavy wall［J］.Theoret.Comput.Fluid Dyn.，1998，11：109-134.

［16］KRETTENAUER K，SCHUMANN U.Numerical solution of turbulent convection over wavy terrain［J］.J.Fluid Mech.，1992，237：261-300.

［17］WALKO R L，COTTON W R，PIELKE R A.Large-eddy simulations of the effects of hilly terrain on the convective boundary layer［J］.Boundary-Layer Met.，1992，58：133-150.

［18］D?RNBRACK A，SCHUMANN U.Numerical simulation of turbulent convective flow over wavy terrain［J］.Boundary-Layer Met.，1993，65：323-355.

［19］DOUGLAS S.HENN and R.IAN SKYES.Large-eddy simulation of flow over wavy surfaces［J］.J.Fluid Mech，1999，383：75-112.

［20］LEWELLEN W.Use of invariant modeling：Handbook of Turbulence，Plenum［M］.New York，1977.

［21］RONALD JCalhoun，ROBERT L Street.Turbulent flow over a wavy surface：Neutral case［J］.Journal of Geophysical Research，2001，106（C5）：9277-9293.

［22］ZANG Y，STREET R，KOSEFF J.A dynamic mixed subgridscale model and its application to turbulent recirculating flow ［J］.Phys.Fluids，1993，5：3186-3196.

［23］SALVETTI M，ZANG Y，STREET R，et al.Large eddy simulation of free-surface decaying turbulence with dynamic subgridscale models［J］.Phys.Fluids，1997，9：2405-2419.

［24］LILLY D K.A proposed modification of the Germano subgridscale closure method［J］.Phys.Fluids A，1992，4：633-635.

［25］JOHN KIM，PARVIZ MOIN，ROBERT MOSER.Turbulent statistics in fully developed channel flow at low Reynolds number［J］.J.Fluid Mech.，1987，177：133-166.

［26］NAKAYAMA A，SAKIO K.Simulation of flows over wavy rough boundaries［R］.Center for Turbulence Research Annual Research Briefs，2002.

Large-eddy simulation of fully developed turbulent flow over a wavy surface

SUN Liming，CAO Shuyang，LI Ming，YANG Qing，ZHANG Enzhen
（State Key Laboratory of Disaster Reduction in Civil Engineering，Tongji University，Shanghai 200092，China）

Large eddy simulation（LES）is used to study flow over a sinusoidal bottom wall，with amplitude-wavelength ratio 2a／λ=0.1 and bulk velocity based Reynolds number Reb=6760.Different from a conventional flat channel，flow over a wavy lower surface is affected by the bottom wall in terms of form drag，leading to the change of both mean and instantaneous fields.Wall bounded flow turns into a separation flow with the increasing of Reynolds number，changing the way in which momentum flux transports and mixes.Statistical properties，transient flow fields as well as surface drag are studied.Averaged fields and turbulent structures are different from those of a flat channel.Distribution and shape of streamwise vortex are closely related to the configuration of the lower surface.Comparison of averaged fields with 2a／λ=0.2 shows the dependence of flow characteristics on the wave steepness.

LES；wavy channel；mean wind profile；channel flow；drag force；sinusoidal wave

O357.5

Adoi：10.7638／kqdlxxb-2012.0172

0258-1825（2014）04-0534-10

2012-10-24；

2012-12-27

國家自然科學(xué)基金 （50978202）；國家自然科學(xué)基金中日科研國際合作項(xiàng)目 （51021140005）

孫麗明（1986），男，江蘇泰州人，碩士，研究方向：橋梁結(jié)構(gòu)抗風(fēng)，強(qiáng)風(fēng)作用下海上風(fēng)剖面數(shù)值模擬.E-mail：sunliming59045＠126.com

曹曙陽（1966），男，江蘇姜堰人，博士、同濟(jì)大學(xué)特聘教授，國際風(fēng)工程協(xié)會秘書長.E-mail：shuyang＠#edu.cn

孫麗明，曹曙陽，李明，等.考慮波浪形底面影響的邊界層風(fēng)場大渦模擬［J］.空氣動力學(xué)學(xué)報，2014，32（4）：534-543.

10.7638／kqdlxxb-2012.0172. SUN L M，CAO S Y，LI M，et al.Large-eddy simulation of fully developed turbulent flow over a wavy surface［J］.ACTA Aerodynamica Sinica，2014，32（4）：534-543.