李世隆，楊曉雷，*，袁先旭，郭啟龍

(1.中國科學院 力學研究所 非線性力學國家重點實驗室，北京 100190；2.中國科學院大學 工程科學學院，北京 100049；3.中國空氣動力研究與發(fā)展中心 計算空氣動力研究所，綿陽 621000；4.空氣動力學國家重點實驗室，綿陽 621000)

0 引 言

工程應用中的壁湍流通常發(fā)生在具有一定粗糙度的表面[1-4]，例如大氣邊界層[5-6]、河流[7]、被侵蝕的螺旋槳槳葉繞流[8]。壁湍流對近壁區(qū)的動量和能量輸運起著至關重要的作用，相關研究可以追溯到Colebrook[9]、Nikuradse[10]和 Moody[11]等學者的經典工作。近期的數(shù)值模擬研究采用隨機旋轉的橢球體作為粗糙元生成粗糙壁表面[12-13]。在本研究中，我們嘗試關注橢球體的取向，即橢球粗糙元隨機旋轉與橢球粗糙元取向一致，如何影響粗糙壁湍流的統(tǒng)計特征。

相比于正弦狀[14]和立方體粗糙元[15]生成的粗糙壁，以橢球體作為粗糙元的粗糙壁湍流研究相對較少。本文所采用的橢球體粗糙元最早見于Scotti在2006年發(fā)表的文章[16]。該橢球體的三個半軸長分別為r、1.4r和 2r（r為橢球體的最小半軸長）。粗糙壁表面通過將橢球體放置在 2r×2r的方格中，并進行隨機旋轉獲得[16]。Scotti[16]的研究結果顯示，在過渡粗糙區(qū)，該粗糙壁對應的等效砂粒粗糙度長度ks≈r。針對橢球粗糙壁槽道湍流，Yuan等[17]系統(tǒng)研究了橢球粗糙度對雷諾應力輸運特性的影響。Yuan等[12]著重研究了橢球粗糙壁湍流在粗糙子層中的統(tǒng)計特性。Hantsis和Piomelli研究了橢球粗糙度對被動標量輸運的影響[18]。針對其他橢球粗糙壁湍流，Yuan和Piomelli研究了該粗糙壁對匯流邊界層的積分量、湍動能輸運特性、湍流結構的影響[19]，開展了橢球粗糙壁的加速邊界層直接數(shù)值模擬，發(fā)現(xiàn)粗糙壁可以避免加速導致的湍流層流化[20]。Wu等[21]研究了橢球粗糙壁沖擊射流的流動特性，并與光滑壁的結果進行了對比。研究者也對比了橢球粗糙壁和其他類型粗糙壁的湍流特性。Jouybari等[22]將橢球粗糙壁的湍流結構與立方體粗糙壁及真實葉片表面粗糙壁的湍流結構進行了對比。Mangavelli等[23]研究了槽道湍流中橢球粗糙度與真實葉片表面粗糙度在脈沖式突然加速情況下的作用。橢球粗糙壁的結果也被用于發(fā)展粗糙壁湍流的模型[24-25]及數(shù)據(jù)驅動的粗糙度長度模型[13]。

為了能夠節(jié)省計算量，Chung等[26]根據(jù)Jiménez和Moin關于光滑壁槽道流的最小槽道的研究[27]，提出了粗糙壁湍流的最小槽道理論。在外層相似性理論的支持下，最小槽道可以較為準確地預測粗糙壁湍流的速度虧損，并發(fā)現(xiàn)槽道的展向尺寸為L+z≈350、流向尺寸L+x≈930時，就可以較為準確地預測平均流向速度廓線及流場結構[28]。

橢球粗糙元的取向影響壁面的阻力分布特征和粗糙子層的流動特征。然而，以上使用橢球體作為粗糙元的文獻均采用隨機旋轉的方式生成粗糙壁表面，難以用于研究橢球粗糙元取向的影響。本工作計劃采用三種具有一致取向的橢球體生成粗糙表面，開展粗糙壁湍流直接數(shù)值模擬，將其統(tǒng)計量與隨機旋轉橢球體粗糙壁的結果進行對比，研究粗糙元取向對粗糙壁湍流統(tǒng)計特性的影響。

1 數(shù)值模擬方法與算例設置

1.1 數(shù)值模擬方法

采用計算流體力學程序VFS-Wind 開展本研究相關的模擬。該程序已成功應用于模擬具有復雜幾何邊界的湍流[29-32]。近期工作中，我們進一步發(fā)展了結合清晰界面浸沒邊界方法和光滑界面浸沒邊界方法的混合方法，使得VFS-Wind程序可開展復雜邊界下顆粒湍流的顆粒解析模擬[33]。需要指出的是，當前工作的模擬只采用了清晰界面浸沒邊界方法，沒有采用上面提及的混合方法。流動的控制方程為不可壓Navier-Stokes方程：

其中：i,j=1,2,3；xi和ui代表三個方向的坐標和速度分量；p為壓強；υ為運動黏性系數(shù)。

曲線浸沒邊界方法（curvilinear immersed boundary method, CURVIB）[34]用于模擬粗糙壁表面的復雜幾何。該方法基于曲線（或笛卡爾）網格模擬流動，采用三角形網格離散復雜幾何邊界。為了在流動模擬中施加邊界條件，將背景網格點標記為流體點和固體點，將至少有一個固體點鄰居的流體點進一步標記為浸沒邊界（immersed boundary, IB）點（如圖1黑實點）。進而，在壁法向通過壁面邊界條件（如圖1中a點）和流體點速度（如圖1中c點）構建IB點速度（如圖1中b點）ub=ua+uchb/hc，作為邊界條件施加給外流，c點的速度由周圍的流體點插值得到，其中ua、ub、uc分別是點a、b、c處的速度。

圖1 曲線浸沒邊界方法處理粗糙壁面邊界的二維示意圖（修改自文獻[32]）Fig.1 Treatment of rough surfaces for the CURVIB method(adapted from[32])

1.2 算例設置

以橢球作為粗糙元生成粗糙壁的表面幾何，共考慮12種粗糙壁表面。其中包括四種不同的粗糙元取向（隨機取向、垂直放置、朝下游傾斜45°、朝上游傾斜45°）。對應每種取向，選取三種不同的流向和展向粗糙元間距（l）。

表1中列出了所有的模擬算例。表中， θz為沿橢球第二長軸（展向方向）旋轉的方向，其取值為0°和45°時分別對應圖2（a）和圖2（b）所示的粗糙元取向，其沒有取值時對應圖2（c）所示的隨機取向。表中，Reτ為基于摩擦速度、槽道高度h和黏性系數(shù)的雷諾數(shù)。在模擬中，采用定流量的方式驅動槽道湍流，所得雷諾數(shù)不完全相同，但基本在1 000附近。

表1 粗糙壁湍流模擬的算例設置Table 1 Numerical setups for simulations of rough-wall turbulence

以下簡要說明算例的命名規(guī)則。每個算例名稱的第一個字母代表粗糙元的形狀，即字母“E”代表橢球體。該橢球體半軸分別為r、1.4r、 2r（對于所有模擬算例，r/h=0.07），其中心位于壁面下方z=-0.5r處。第一個字母后面的數(shù)字表示粗糙元的間距，例如，E20D0Ry 中的 20 表示l= 2.0r。字母“D”表示粗糙元的取向。字母組合“Rn”（0°）代表粗糙元的初始取向，如圖2（a）長軸和短軸分別沿著壁法向方向和流向方向。字母數(shù)字組合“p45”和“n45”分別代表橢球粗糙元沿第二長軸（沿展向方向）正向和負向旋轉45°，如圖2（b）所示。字母組合“Ry”代表橢球的隨機旋轉，如圖2（c）所示。如果橢球較近，則相鄰橢球之間允許重合。圖3顯示了所有模擬粗糙壁表面幾何的局部示意圖。

圖2 橢球粗糙元取向的示意圖Fig.2 Sketch of ellipsoid orientation

圖3 粗糙壁幾何的局部示意圖（所顯示區(qū)域在流向和展向的尺寸分別為h和 0.5h）Fig.3 Sketch of rough surfaces in an area with size h × 0.5h

采用流向（x）、壁法向（y）和展向（z）尺寸為 3h、h、h的全尺寸槽道，該槽道尺寸小于通常采用的槽道尺寸。但計算域大小仍約為Chung等測試的最小槽道的9倍，在保證解析近壁流場結構和容納足夠多粗糙元的同時，降低了計算量，使得大量不同粗糙壁湍流的直接數(shù)值模擬成為可能。小槽道與全尺寸槽道結果的對比見附錄。所模擬算例的計算域在流向（x）、壁法向（y）和展向（z）對應的網格節(jié)點數(shù)分別為513、300和172。離開壁面第一個網格的寬度為其中uτ為摩擦速度。在其他兩個方向上，無量綱的網格寬度約為 Δx+=Δz+=6。時間步長為 2 ×10-3h/Ub，其中Ub為流向截面的平均速度。在粗糙壁表面，通過CURVIB方法施加無滑移邊界條件，在上邊界施加自由滑移邊界條件，在水平方向施加周期邊界條件。通過在流向施加保證定流量的平均壓力梯度驅動流動。在所有模擬中，首先將流場發(fā)展至充分發(fā)展狀態(tài)，再求解 6 0h/Ub平均值，計算湍流統(tǒng)計量。

2 模擬結果

本節(jié)介紹不同粗糙壁的模擬結果，包括粗糙壁附近的流場結構、平均速度、雷諾應力、分散應力、零平面位移、等效砂粒粗糙度長度等。分散應力描述了時間平均流場在空間上的非均勻程度，通過將瞬時速度ui進行如下分解得到：其 中上 加橫 線ˉ·表 示 時 間 平均，表示水平方向的空間平均。近壁處的空間平均包括固體和流體。

圖4顯示了不同粗糙壁附近的瞬時流動結構?？梢钥吹剑鲃咏Y構的復雜性與橢球體取向的隨機性密切相關。對于具有相同取向的粗糙元生成的粗糙表面，通?？梢杂^察更為規(guī)則的流動結構，例如E20D0Rn、E20Dp45Rn。而對于隨機取向橢球體生成的粗糙壁表面，即E20D0Ry、E28D0Ry和E35D0Ry，由于粗糙元周圍流動的差異以及相鄰粗糙元間的相互作用，可以觀察到更為復雜的流動結構。如圖4所示，粗糙元周圍流動的一個重要特征是流動沿粗糙元迎風面的向上運動。一個有趣的現(xiàn)象是，當橢球體朝上游或朝下游傾斜時，其上方的流動結構具有一定的相似性。當間距較大（l=2.8r，3.5r），橢球體朝下游旋轉時（E28Dp45Rn、E35Dp45Rn）的流動結構相較于朝上游旋轉時（E28Dn45Rn、E35Dn45Rn）的具有更高幅值的流向速度。

圖4 使用 Lambda2 準則（λ2 = -0.5）可視化的近壁流動結構（等值面使用時間平均的流向速度著色，圖像中僅顯示了 一部分區(qū)域）Fig.4 Near-wall flow structures visualized by the λ2-criteria (λ2 = -0.5).Iso-surfaces are colored by the time-averaged streamwise velocity.Only part of the flow filed is shown

為了展示粗糙元尾跡的三維幾何特征，圖5顯示了時間平均流向速度u=0的等值面?？梢钥吹?，隨機旋轉橢球元素尾跡的形狀與相同取向橢球體的尾跡形狀有顯著不同。另一方面，粗糙元間距是影響尾跡形狀的關鍵因素。對于本文所考慮的粗糙元間距，流向速度小于0的尾跡區(qū)域都直接作用于下游粗糙單元。在粗糙元間距較小時（l=2.0r，2.8r），相鄰行的尾跡區(qū)域相互連接；在粗糙元間距較大時（l=3.5r），相鄰行的尾跡各自獨立。

圖5 使用流向速度u = 0 等值面標記的尾跡形狀（為清晰起見，未顯示粗糙度單元表面的等值面）Fig.5 Wake shapes illustrated by u = 0.The iso-surface u = 0 at the surface of the roughness element is not shown for clarity

圖6顯示了不同粗糙壁算例的平均流向速度廓線?？梢钥吹剑兴憷季哂幸粋€清晰的對數(shù)區(qū)域。在該區(qū)域，粗糙壁的流向速度低于光滑壁，該速度差被稱為粗糙度函數(shù)ΔU。圖6顯示，當粗糙度單元間距為 2.0r時，隨機旋轉橢球粗糙壁（E20D0Ry）的ΔU高于其他算例；而對于高一些的粗糙元間距（l=2.8r，3.5r），垂直放置橢球粗糙壁（E28D0Rn、E35D0Rn）的ΔU高于其他算例。當以相同的方式朝上游或下游傾斜橢球體時，粗糙度函數(shù)比垂直放置的有所降低（例如，E20Dp45Rn或E20Dn45Rn與E20D0Rn相比）。一個有趣的現(xiàn)象是，橢球朝上游傾斜（E20Dn45Rn）時的 ΔU與朝下游傾斜（E20Dp45Rn）時的ΔU基本相同。對于所有粗糙元取向，ΔU均隨粗糙元間距增加而增加。

圖6 不同粗糙表面的平均流向速度廓線（虛線為光滑壁的線性律和對數(shù)律，不同粗糙壁參數(shù)見表1）Fig.6 Vertical profiles of the mean streamwise velocity computed from cases with different rough surfaces.The dashed lines show the linear profile and the logarithmic law for smooth walls.Parameters of rough walls are shown in Table 1

圖7和圖8定量考察了砂粒粗糙度ks和零平面位移d隨粗糙元間距l(xiāng)的變化規(guī)律。零平面位移是指平均流向速度為零的位置，代表粗糙壁面對邊界層的抬升作用，實際中可能并不準確是速度為零的位置[6]。這里，基于對數(shù)律計算d和ks的值，即，8.5，其中U+為通過壁面摩擦速度無量綱化的平均流向速度， κ =0.4 為卡門常數(shù)。可將對數(shù)率化為y=kseκ(U+-8.5)+d， 這 是 關 于ks和d的 線 性 式 。 將0.15＜y＜0.3之間約20個數(shù)據(jù)點的y與U+作為觀測數(shù)據(jù)，將ks與d作為優(yōu)化參數(shù)，可以得到與對數(shù)率較為符合的ks與d。壁面摩擦速度通過驅動流動的流向平均壓力梯度計算獲得。在圖7中可以看到，由于橢球粗糙表面迎著正流向速度的面積增加，ks隨著單元間距的增加而增加。橢球垂直放置粗糙壁的增長率最高，而橢球隨機取向壁面的增長率最低。另一方面，朝上游和朝下游傾斜橢球壁面的ks和增長率大致相同。間距為 2.8r或 3.5r時，橢球垂直放置粗糙壁的ks值高于其他橢球取向粗糙壁。

圖7 砂粒粗糙度 ks 隨粗糙元間距的變化Fig.7 Variations of sandgrain roughnessks with roughness element spacing

圖8 零平面位移 d 隨粗糙元間距的變化Fig.8 Variations of zero-plane displacementd with element spacing

不同橢球粗糙壁的零平面位移d隨粗糙元間距l(xiāng)的變化如圖8所示?？梢钥吹?，d值隨l的增加而減小。在所有橢球粗糙壁中，垂直取向橢球粗糙壁的d值最大。一個有趣的現(xiàn)象是，當粗糙度間距為 2.8r或3.5r時，橢球朝上游傾斜粗糙壁的d值高于橢球朝下游傾斜粗糙壁，而ks值基本相同。

圖9比較了不同橢球粗糙壁的雷諾正應力，不同粗糙壁參數(shù)見表1?？梢钥吹?，雷諾正應力的三個分量按大小從高到低依次為的最大值位于靠近粗糙壁的位置，而的最大值位置離粗糙壁稍遠。對于不同粗糙元間距，不同粗糙壁的和的最大值相近，而的最大值隨粗糙元間距增加而減小。與隨機取向橢球粗糙壁相比，相同取向橢球粗糙壁的雷諾正應力的峰值更高。朝上游和下游傾斜橢球粗糙壁的峰值基本相同，且高于其他粗糙壁的。圖10顯示了不同粗糙壁面的流向分散應力。流向分散應力的峰值幅值比其他兩個分量高約一個數(shù)量級。分散應力衡量了時間平均速度場的空間不均勻程度?？梢钥吹剑煌植诒诘姆稚τ酗@著差異。在粗糙元間距為 2.0r、 3.5r時，與橢球體垂直放置的粗糙壁相比（E20D0Ry、E20D0Rn），隨機旋轉橢球體粗糙壁的的最大值更高；在粗糙度單元間距為2.8r時，其最大值基本相同。由圖9可見，橢球體朝上游和下游傾斜45°粗糙壁的雷諾應力大致相同（例如E20Dp45Rn 與 E20Dn45Rn）。另一方面，圖10顯示，其對分散應力的影響是不同的：當朝下游傾斜45°時，的最大值增加（E20Dp45Rn）；而在朝上游傾斜45°時，的最大值降低（E20Dn45Rn）。在增加粗糙元間距時，朝上游傾斜45°時的的最大值幅值逐步增大，而對于其他橢球取向粗糙壁，的最大值幅值增大或減小的趨勢并不一致。其原因可能在于，朝上游或朝下游傾斜45°的影響主要局限于粗糙子區(qū)，使得分散應力有較大不同。而對于外區(qū)流動，兩種粗糙壁的平均高度相同，正面密實度和水平密實度也相同，使得水平方向平均的湍流統(tǒng)計量較為接近。另一方面，對于所有橢球粗糙壁，具有較大幅值的范圍均隨粗糙元間距的增大而逐步增大。

圖10 不同粗糙表面的流向分散應力的垂向分布Fig.10 Vertical profiles of the streamwise dispersive stress computed from cases with different rough surfaces

圖11 不同粗糙表面的垂向和橫向分散應力的垂向分布Fig.11 Vertical profiles of the vertical and spanwise dispersive stresses computed from the cases with different rough surfaces

3 結 論

開展了以橢球體作為粗糙元的粗糙壁槽道湍流直接數(shù)值模擬，研究了橢球體取向對粗糙壁槽道湍流統(tǒng)計特性的影響。研究得到以下結論：

1）直接模擬結果顯示粗糙元取向對湍流統(tǒng)計量有顯著影響。對于小粗糙元間距（l=2.0r，r為橢球體的最小半軸長），隨機取向粗糙壁的砂粒粗糙度長度ks高于其他相同取向的粗糙壁；對于大粗糙元間距（l=2.8r，3.5r），垂直放置粗糙壁的ks高于其他粗糙壁。

2）將粗糙元朝上游/下游傾斜45°時的ks基本相同，并且這兩種粗糙壁的雷諾正應力（以壁面摩擦速度無量綱化）也基本相同，且高于其他粗糙壁。

3）粗糙元朝上游/下游傾斜45°時的粗糙壁的流向分散應力顯著不同，其中，粗糙元朝下游傾斜45°時的流向分散應力更高，且遠高于其他粗糙壁。

4）對于零平面位移，垂直放置橢球粗糙壁的值明顯高于其他粗糙壁。增加粗糙元間距，粗糙壁的流向雷諾正應力的最大值逐步降低，而其他方向的雷諾正應力沒有明顯的一致性趨勢。

5）對于垂直取向粗糙壁，流向分散應力隨著粗糙元的間距增加而增加；對于其他粗糙壁，沒有發(fā)現(xiàn)隨粗糙元間距變化的單調變化趨勢。

6）垂向和橫向分散應力小于流向分散應力近一個數(shù)量級。隨機取向橢球粗糙壁的垂向和橫向分散應力顯著高于其他由相同取向橢球體形成的粗糙壁。

本文研究了粗糙元取向對粗糙壁槽道湍流統(tǒng)計量的影響。更為深入的機理研究及基于這些機理研究發(fā)展粗糙壁湍流的工程模型將在未來工作中開展，并將進一步探索不同粗糙壁是否存在統(tǒng)一的標度形式。

附錄A 槽道尺寸驗證

本附錄給出E20D0Ry小槽道模擬結果與全尺寸槽道結果的對比。圖A1為小槽道和全尺寸槽道的粗糙壁表面示意圖。全尺寸槽道的流向、壁法向、展向尺寸分別為6h、h、3h，對應的網格數(shù)分別為1 025、300、513。

圖 A1 全尺寸槽道與小槽道的粗糙壁表面示意圖（使用粗糙壁面的高度染色）Fig.A1 Sketch of the rough surfaces (colored by height)for the full-size channel and minimal channel

圖A2對比了得到的平均速度廓線、雷諾應力和分散應力?？梢钥吹?，平均速度廓線在對數(shù)區(qū)吻合較好，而在接近槽道中心處，小槽道預測平均速度廓線略高。在近壁面，平均流向速度廓線、雷諾應力和分散應力有一定的差別。

以下，嘗試分析產生上述差別的原因。對于槽道中心處平均速度廓線的差別，小槽道在水平方向的尺寸有較大影響。對于近壁面的平均速度和湍流統(tǒng)計量差別，粗糙壁表面統(tǒng)計量的差別可能是主要原因。E20D0Ry粗糙壁面由隨機旋轉橢球組成，而小槽道中包含的粗糙元個數(shù)較少，使得其統(tǒng)計量難以與全尺寸槽道粗糙表面一致。以粗糙壁面平均高度kˉ為例，全尺寸槽道kˉ =0.047h，而小槽道kˉ =0.052h。為了定量評估kˉ的差異的影響，檢查了等效砂粒粗糙度長度和最大流向雷諾正應力隨kˉ的變化。其中，所采用的Yuan和Piomelli[17]的結果來自于以同樣方式生成的粗糙壁表面的槽道湍流直接數(shù)值模擬。在他們的工作中，壁面的復雜幾何形狀采用VOF（volume of fluid）方法處理。從圖A3可以看到，等效沙粒粗糙度長度ks和最大流向雷諾正應力基本符合全尺寸槽道結果的變化趨勢。從而證明，kˉ是造成上述差別的主要原因。對于其他兩個方向的雷諾正應力和雷諾剪切應力，小尺寸槽道和全尺寸槽道結果的差別相對較小。對于分散正應力，流向和展向分量的差異相對較大。差異主要體現(xiàn)在小槽道的峰值較小，位置離壁面較遠，這與小槽道kˉ 較大是一致的。綜上，小尺寸槽道可以較準確地預測近壁面和對數(shù)區(qū)的流動，可以用于本文的相關研究。

圖 A2 小槽道與全尺寸槽道湍流的結果對比Fig.A2 Comparisons between minimal and full-size channels

圖 A3 砂粒粗糙度和流向雷諾正應力最大值隨粗糙壁面平均高度的變化（黑色圓點為Yuan和Piomelli[17]不同橢球半徑的結果，紅色方塊為全尺寸槽道的結果，藍色三角為小槽道的結果）Fig.A3 Variations of equivalent sandgrain roughness length and the maximum value of the streamwise Reynolds normal stress with the average height of rough surfaces.Black circles sare the results from Yuan and Piomelli[17] of different semi-axis lengths,red squares and blue triangles are respectively the results from full-size and minimal channels