基于諧波合成法的大渦模擬脈動風場生成方法研究*

沈 煉 ,韓 艷? ,蔡春聲,2 ,董國朝

(1. 長沙理工大學 土木與建筑學院，湖南 長沙 410114；2. 路易斯安那州立大學 土木與環境工程，美國路易斯安那州，巴吞魯日 70803)

基于諧波合成法的大渦模擬脈動風場生成方法研究*

沈 煉1,韓 艷1?,蔡春聲1,2,董國朝1

(1. 長沙理工大學 土木與建筑學院，湖南 長沙 410114；2. 路易斯安那州立大學 土木與環境工程，美國路易斯安那州，巴吞魯日 70803)

為準確模擬大渦模擬入口處的脈動信息，在充分考慮脈動風場的功率譜、相關性、風剖面等參數前提下，運用諧波合成方法生成了滿足目標風場湍流特性的隨機序列數，通過對FLUENT軟件平臺進行二次開發，將生成的隨機序列數賦給大渦模擬的入口邊界，從而實現了大渦模擬的脈動輸入.基于風洞試驗數據，分別建立了兩種模擬脈動風場的數值模型，一為沒有任何障礙物的空風洞，運用諧波合成方法生成的隨機序列數作為入口邊界來生成脈動信息；二為與真實風洞一致的尖劈粗糙元風洞，采用平均風作為入口邊界，利用尖劈粗糙元對風場的擾動來產生脈動信息.通過對比兩種數值模型發現：基于諧波合成方法生成的脈動風場可作為大渦模擬的入口邊界，可為大渦模擬脈動入口研究提供參考.

大渦模擬；諧波合成；數值模擬；脈動風場

大渦模擬在計算資源與計算效率方面是雷諾時均模擬(RANS)和直接模擬(DNS)的折中，集成了二者優勢，因此在CFD模擬方面得到了廣泛應用.但在實際運用過程中，大渦模擬入口邊界處的脈動信息給定問題還沒有得到完全解決[1].目前，許多學者直接將均勻來流賦給入口邊界，沒有考慮脈動信息的作用．而在實際風場中，來流脈動對結構的脈動風壓和瞬時風荷載起著非常重要的作用[2]，有著不可忽略的地位.因此，非常有必要在大渦模擬過程中考慮脈動信息的影響.

針對大渦模擬脈動入口研究方法，可以將其分為三類[3]，第一類為“預前模擬法”，是運用一個預前模擬區域，將目標風場預先模擬出來，然后賦給主模擬區域的入口邊界.朱偉亮等[2]在預前模擬區域用流場循環的方法生成了脈動風場，并對鈍體平面屋蓋的結構風壓進行了模擬，這種脈動風場生成方法的不足之處是需要生成匹配的數據庫，會消耗巨大的內存和計算時間.序列合成法的另一種儲存形式是用空間序列儲存，Chung[4]利用泰勒方法協調了預前模擬方法兩個計算域之間的離散差異性，對槽道流進行了模擬，運用該方法的缺點是不同時刻的速度不能嚴格的由泰勒級數展開來獲取，因此插值的效果不佳.第二類方法為渦方法，針對渦方法的研究，Mathey等[5]人用一系列的二維隨機漩渦加載到平均風場中去從而得到脈動風場，這種方法生成的脈動風場雖然具有較好的空間相關性和湍動能信息，但不足之處是目標風場很難滿足風場各向異性目標譜的要求.第三類方法為序列合成法，序列合成法是運用傅里葉變換的方法來生成脈動風速時程.2002年Hanna等[6]人提出了一種簡單的各向異性湍流脈動生成方法，生成了隨機正態序列數，它是假設在順風方向和垂直方向的脈動均方根和平均速度有不同的比例，然后通過加權得到脈動信息，這種方法的不足是缺少了雷諾應力信息，同時湍流也會很快衰減.此后，2006年，Kondo[7]等提出了動態湍流輸出方法，但這種方法不能滿足目標譜的需求.2008年，Xie和Castro[8]用雷諾應力的垂向分布構造出了脈動速度分量，得到的脈動量滿足湍流的基本特性并具有雷諾應力信息，但是生成的湍流在入口處并不滿足連續性的需求，需要經過很長時間的發展才能形成真正的湍流.后來，Hemon and Santi[9],Huang等[10-11]也針對大渦模擬的入口脈動邊界做出了相關研究，得到了不錯的計算結果.

本文借助于經典的諧波合成方法，在充分考慮湍流的平均風剖面、目標譜、空間相關性等指標下，合成一系列的隨機序列時程數據，把生成的時程數據與大渦模擬的入口建立對應關系，基于對FLUENT軟件平臺進行二次開發，把生成的隨機序列數賦給大渦模擬的入口邊界，從而生成入口脈動信息.最后通過對比數值模型和風洞試驗結果發現：基于諧波合成的脈動信息生成方法能滿足大渦模擬入口邊界的要求，其優勢在于目標譜及相關性參數能根據不同情況的需要進行隨意調整，合成的時程數據可直接賦給入口邊界，生成速度快，模擬時間長短可以隨意控制，模擬的點都相互獨立，可以很好地進行并行計算，大大提高計算效率.

1 數值模擬理論與方法

1.1 大渦模擬

(1)

式中:D為流動區域，x'為實際流動區域里面的空間坐標，x是過濾后的大尺度空間坐標，G(x,x')是過濾函數[13]，表達式為：

(2)

式中：V是控制體積所占幾何空間的大小，過濾后瞬態下的空間N-S方程可表示為：

(3a)

(3b)

上式為動量守恒方程和質量守恒方程，其中帶橫線上標的量表示過濾后的可解尺度量，τij為亞格子尺度應力，具體表達式為：

(4)

為了封閉方程(3)，本文采用Smagorinsky亞格子模型，假定SGS應力形式為：

(5)

在求解方面，本文的N-S方程采用建立在交錯網格上的SIMPLE方法進行求解，對流項和擴散項均采用二階中心差分格式，時間導數采用二階半隱式離散方法，用超松弛方法(SOR)求解壓力Poisson方程，壓力和動量松弛因子分別取0.3和0.7.

1.2 隨機序列合成

(6)

式中：H(ω)為下三角矩陣；對角項為ω的實非負函數，非對角項通常為ω的復函數.

(7)

(8)

脈動風速的合成還需考慮平均風、功率譜、相關性和時間步長的影響.入口風速時程的平均速度采用指數率風剖面公式:

(9)

本文脈動風功率譜采用我國規范采用的Kaimal譜，其表達式為：

(10)

本文空間互相關性采用Davenport給定的經驗公式[16]，如下所示：

(11)

(12)

式中：Sii和Sjj分別為i,j兩點的自譜密度函數.本文對隨機序列數合成、UDF導入和大渦模擬3個過程采用了相同的時間步.在隨機序列數合成過程中，為了防止頻率的失真，時間步ΔT應該滿足：

(13)

在UDF時程數據對接上，要保持計算步的統一，而在大渦模擬過程中，計算步長[17]應滿足庫朗數(CFL)要求，可表示為：

(14)

式中：U為風速；Δx為網格尺寸；ΔT為時間步長，為同時滿足三方面的要求，時間步長取0.002 5 s.

在充分考慮脈動風場的平均速度、功率譜、相關性、時間步長等參數后，基于諧波合成理論，用MATLAB對風速時程數據進行合成，生成滿足目標風場的隨機序列數.然后將生成的隨機序列數與入口網格在時間和空間上建立對應關系，時間上要滿足時間步長的一致，保證諧波合成風速與大渦模擬風速在時間上的同步.空間上要對入口網格坐標進行嚴格控制，將模擬點的速度時程賦給對應的入口網格中心點，整個過程用UDF程序編制，對FLUENT軟件用戶自定義模塊進行加載，重新定義入口速度，從而實現大渦模擬入口脈動信息的輸入.

2 大渦模擬入口脈動風場特性分析

為驗證本文入口處脈動信息輸入的正確性，建立了一空風洞數值模型，里面沒有任何障礙物，具體尺寸為15 m×2.5 m×3.0 m (x×y×z)，總網格數為180萬.入口平均風采用指數率風剖面形式，α值為0.16，取基準風速為20 m/s，地表粗糙高度為0.05 m.入口脈動量由MATLAB程序模擬生成，然后通過UDF程序將這些時程數據賦給入口邊界網格，待計算穩定后提取模型入口0.5 m和1 m高度處的風速時程并與輸入的諧波合成法生成的數據序列對比，如圖1～2所示.

時間/s

時間/s

由圖1～2可以看出，模型入口監測值與輸入的諧波合成法生成的數據序列基本一致，說明了入口邊界風速時程輸入方法的正確性.為了驗證入口邊界功率譜的正確性，將數值模型入口1 m高度處順風方向功率譜模擬值與我國規范進行對比，如圖3所示.

頻率/(rad·s-1)(對數坐標)

從圖3可以看出，功率譜能較好的與我國規范采用的Kaimal譜吻合，滿足脈動風場特性要求.

對流場入口中心0.5 m和1 m高度處的速度時程數據進行相關性分析，定義兩點分別為A和B.圖4給出了入口監測點A和B的互相關性和B點的自相關性檢驗結果，圖4(a)中目標值由公式(11)給出.

時間/s

時間/s

圖4表明了入口處脈動風場相關性呈指數率衰減，模擬值與理論目標值衰減趨勢一致，體現出了模擬風場在入口處具有良好的相關性.

通過對上述平均風速、功率譜、相關性等參數進行分析，發現基于諧波合成法生成的脈動風場滿足湍流的基本特性，可以在大渦模擬的入口進行無縫對接.

3 數值模型

為驗證數值模型內部風場正確性，以TJ-2風洞試驗數據[18]為參照對象，建立了兩個數值模型，一為沒有任何障礙物的空風洞，與第二節所用空風洞驗證模型為同一計算模型，利用上節所提出的諧波合成方法來生成入口脈動信息，簡稱為“空風洞數值模型”；二為與風洞試驗條件一致的“尖劈粗糙元數值模型”.下面分別對兩種模型的尺寸，網格數量，計算參數與邊界條件進行介紹.

3.1 模型尺寸與網格數量

3.1.1 空風洞數值模型

TJ-2風洞具體尺寸為15 m×2.5 m×3.0 m (長×高×寬)，空風洞數值模擬采用與物理風洞一致的模型尺寸，為保證計算精度，采用全六面體網格，并對總網格數為90萬、180萬、400萬3種網格進行了網格無關性測試，權衡計算資源和計算精度后，選取了180萬網格數模型作為最終計算模型，網格在模型底部進行加密，增長因子為1.1，最底層網格高度為0.01 m，具體網格如圖5所示.

圖5 空風洞數值模型網格示意圖

3.1.2 尖劈粗糙元數值模型

尖劈粗糙元數值模型與物理風洞試驗條件保持一致，布置如圖6所示.與上一節類似，為保證精度，對尖劈粗糙元數值模型進行全場六面體網格劃分，并對260萬，450萬，800萬3種網格進行了網格無關性測試，無關性測試結果如表1所示.

通過比較模型內部測量中心處湍流度發現260萬網格數模型較450萬和800萬的偏差較大，但450萬和800萬網格數的計算結果偏差幅度很小，權衡計算精度和計算資源的影響后，最終選取了450萬網格數作為最終計算網格.整個網格在底部進行加密，最底層網格高度為0.005 m，延伸率為1.1，網格如圖7所示.計算過程中，對空風洞數值模型和尖劈粗糙元數值模型的“預定模型中心”位置進行風速監測，具體位置如圖6所示.

表1 網格無關性檢驗

Tab．1 Test of the grid dependence

網格數量260萬450萬800萬湍流度9.9%12.42%12.49%

圖6 尖劈粗糙元模型尺寸示意圖

圖7 尖劈粗糙元數值模型網格示意圖

3.2 邊界條件設置

空風洞數值模型的入口條件是基于對FLUENT軟件平臺進行二次開發，將生成的隨機序列數賦給入口邊界；而尖劈粗糙元數值模型則采用平均風作為入口邊界，詳細邊界條件與計算參數如表2所示.

4 計算結果

空風洞數值模型采用超線程12核工作站進行計算，模擬150 s時長(60 000步)需耗時60 h，其計算結果如圖8(a)所示.尖劈粗糙元模型同樣采用超線程工作站進行計算，模擬同樣時長需耗時124 h，計算結果如圖8(b)所示.

表2 計算參數與邊界條件

Tab．2 Computational parameters and boundary conditions

計算參數空風洞模型尖劈粗糙元模型湍流模型大渦模擬大渦模擬網格類型六面體結構網格六面體結構網格入口邊界UDF導入的隨機序列數均勻來流,U=20m/s計算域側面對稱邊界對稱邊界計算域頂面自由滑移自由滑移計算域底面無滑移固壁無滑移固壁尖劈粗糙元表面無滑移固壁無滑移固壁

(a) 空風洞數值模型計算云圖

(b)尖劈粗糙元數值模型計算云圖

通過對 “預定模型中心”不同高度處的風速監測，得到監測中心處的平均風剖面與湍流度信息，對其無量綱處理并與風洞試驗結果進行對比，其結果如下圖所示.通過比較圖9和圖10發現，在風剖面和湍流度方面，兩種數值模型相對風洞試驗值的偏差均較小，能滿足工程精度需要，也再一次證明了諧波合成方法生成脈動風場的正確性.用最小二乘法對空風洞模型風剖面進行指數率擬合，得到α指數為0.161，非常接近風洞試驗的擬合值0.162，也接近于我國規范給出的B類風場特性.

相對風速U/Ur

湍流度Iu

對“預定模型中心處”1 m高度處的時程數據進行分析，得到空風洞模型和尖劈粗糙元模型的順風方向功率譜，將其進行擬合并與Kaimal譜進行對比，如圖11所示.

頻率/Hz

通過對兩種方法同一監測點順風方向功率譜與Kaimal譜進行對比發現，數值模擬結果與風洞試驗的風速功率譜在低頻段(頻率<1)吻合較好，但在高頻段(頻率>1)出現陡降，而工程領域最為關心的頻率出現在慣性子區(0.2<頻率<1)，本文兩種數值模擬方法均能較好的捕捉該區域的風譜，因此所生成的脈動風能較好的滿足工程需要.目前，數值模擬的風速時程普遍存在風譜高頻段衰減現象，是由于高頻段代表風場中小渦的能量貢獻，而LES中小尺度渦是用亞格子模型進行封閉，不能直接對小尺度渦進行求解，因此數值模擬不能捕捉到脈動風場的高頻段風譜，通過改進湍流模型和加密網格數量可以改善風譜高頻衰減問題.

在計算效率方面，將兩種數值模型所需的網格數和耗時量作出對比，如表3所示.

表3 計算效率對比

Tab．3 Efficiency of the computational

模型網格數/萬計算時間/h空風洞模型18060尖劈粗糙元模型450124

由表3可以發現，在達到允許精度要求下，空風洞模型所需的網格數不到尖劈粗糙元模型的一半，計算耗時量方面也遠遠優于尖劈粗糙元模型.

此外，空風洞數值模型具有更好的適應性，特別是對不同脈動風場的模擬，空風洞模型只需對目標風場的特性參數進行調整，操作簡單方便，可調性強.而尖劈粗糙元數值模型只能通過不斷調整尖劈粗糙元的位置和數量來對目標風場進行試算，計算工作量大，消耗時間長.

5 結 論

本文基于諧波合成方法生成了隨機序列數，通過對商業軟件FLUENT進行二次開發，利用UDF程序把時程數據轉化成大渦模擬的入口脈動信息，同時建立了與物理風洞試驗條件一致的尖劈粗糙元數值模型，通過對比研究得到了以下成果和結論.

1)實現了隨機序列數在大渦模擬入口處的無縫對接，開發了大渦模擬入口脈動速度輸入模塊，可以很好地進行并行計算，解決了大渦模擬入口處湍流信息給定問題.

2)通過比較平均風剖面、湍流度和功率譜等結果發現：基于諧波合成方法生成的大渦模擬脈動風場具有較高的精度，能滿足工程需要，但風速功率譜在高頻段會出現衰減現象.

3)通過對比空風洞數值模型和尖劈粗糙元模型的計算效率發現，空風洞數值模型在網格數量和計算耗時方面均要遠遠優于尖劈粗糙元數值模型.

4)相比尖劈粗糙元數值模型，基于諧波合成脈動生成方法對不同風場模擬只需調整風場特性參數，操作方便，模擬速度快，是一種極具前景的脈動生成方法，也可為大渦模擬入口處脈動生成方法精細化研究提供參考.

基于諧波合成法的脈動風場生成方法是譜合成方法的一種拓展，擁有許多優點的同時也存在一些不足，如生成的隨機風速時程不滿足N-S方程的連續性，風譜高頻衰減等問題還有待進一步解決.

[1] JIANG Wei-mei，MIAO Shi-guang. 30 years review and perspective of the research on the large eddy simulation and atmospheric boundary layer[J]. Advances in Nature Sciences, 2004, 14: 11-19.

[2] 朱偉亮，楊慶山. 基于 LES 模型的近地脈動風場數值模擬[J]. 工程力學，2010，27(9)：17-21.

ZHU Wei-liang, YANG Qing-shan. Large eddy simulation of near ground turbulent wind field [J]. Engineering Mechanics, 2010, 27 (9): 17-21. (In Chinese)

[3] BABA-AHMADI M H, TABOR G. Inlet conditions for LES using mapping and feedback control[J]. Computer & Fluid, 2009, 38(6): 1229-1311.

[4] CHUNG Y M, SUNG H J. Comparative study of inflow conditions for spatially evolving simulation[J]. AIAA Journal, 1997, 35(2): 269-274.

[5] MATHEY F, COKLJAT D, BERTOGLIO J P,etal. Assessment of the vortex method for large eddy simulation inlet conditions[J]. Progress in Computational Fluid Dynamics, 2006, 6: 1-3.

[6] HANNA S R, TEHRANIAN S, CARISSIMO B,etal. Comparisons of model simulations with observations of mean flow and turbulence within simple obstacle arrays[J]. Atmospheric Environment, 2002, 36: 5067-5079.

[7] KONDO H, ASAHI K, TOMIZUKA T,etal. Numerical analysis of diffusion around a suspended expressway by a multi-scale CFD model[J]. Atmospheric Environment, 2006, 40: 2852-2859.

[8] XIE Z T, CASTRO I P. Efficient generation of inflow conditions for large eddy simulation of street-scale flows[J]. Flow Turbulence and Combustion, 2008, 81: 449-470.

[9] HEMON P, SANTI F. Simulation of a spatially correlated turbulent velocity field using biorthogonal decomposition[J]. Journal of Wind Engineering and Industrial Aerodynamics, 2007, 95 (1): 21-29.

[10]HUANG S H, LI Q S, WU J R. A general inflow turbulence generator for large eddy simulation[J]. Journal of Wind Engineering and Industrial, 2010, 98: 600-617.

[11]盧春玲，李秋勝，黃生洪，等. 大跨度屋蓋風荷載的大渦模擬研究[J]. 湖南大學學報：自然科學版，2010，37(10):7-12.

LU Chun-ling, LI Qiu-sheng, HUANG Sheng-hong,etal. Large eddy simulation of wind loads on long-span roofs[J]. Journal of Hunan University : Natural Sciences, 2010,37 (10): 7-12. (In Chinese)

[12]DEARDOFF J W. The use of sub grid transport equation in a three dimensional model of atmospheric turbulence[J]. ASME Journal of Fluid Engineering, 1973, 95: 429-438.

[13]祝志文，陳魏，向澤，等. 大跨度斜拉橋主梁氣動力特性的大渦模擬[J]. 湖南大學學報：自然科學版，2013，40(11): 26-33.

ZHU Zhi-wen, CHEN Wei, XIANG Ze,etal. Large eddy simulation of aerodynamics on main girder of long span cable-stayed bridges[J]. Journal of Hunan University: Nature Sciences, 2013, 40 (11): 26-33. (In Chinese)

[14]王福軍. 計算流體動力學分析[M].北京：清華大學出版社，2004.

WANG Fu-jun. Computational fluid dynamics analysis [M]. Beijing: Tsinghua University Press, 2004. (In Chinese)

[15]JCR H, WRAY A, MOIN P. Eddies, streams, and convergence zones in turbulent flows[R]. Center for turbulence research report CTR-S88, 1988: 193-208.

[16]韓艷. 橋梁結構復氣動導納函數與抖振精細化研究[D]. 長沙：湖南大學土木工程學院，2007.

HAN Yan. Study on complex aerodynamic admittance functions and refined analysis of buffeting response of bridge[D]. Changsha: College of Civil Engineering, Hunan Univ, 2007. (In Chinese)

[17]Fluent Inc. Fluent user's guide, version：ansys12.0 [M]. USA：Fluent Inc Lebanon, 2009.

[18]龐加斌，林志興. 關于風洞中用尖劈和粗糙元模擬大氣邊界層的討論[J]. 流體力學實驗與測量，2004，18(2)：32-37.

PANG Jia-bin, LIN Zhi-xing. Discussion on the simulation of atmospheric boundary layer with spires and roughness elements in wind tunnels[J]. Experiments and Measurements in Fluid Mechanics, 2004, 18(2): 32-37. (In Chinese)

Research on Generating Method of Fluctuating Wind Field of LES Based on WAWS

SHEN Lian1，HAN Yan1?, CAI Chun-sheng1, 2, DONG Guo-chao1

(1. School of Civil Engineering and Architecture，Changsha Univ of Science Technology，Changsha,Hunan 410114，China; 2.Department of Civil and Environmental Engineering，Louisiana State Univ，Baton Rouge，USA，LA 70803)

In order to accurately simulate the inlet boundary conditions of turbulence information of Large Eddy Simulation (LES), the fluctuating time-history data satisfying the target wind field were simulated with the weighted amplitude wave superposition (WAWS) considering the integral scale, power spectrum, turbulence intensity and wind velocity profile. The fluctuating time-history data were given to the inlet boundary of LES by the secondary development of the commercial software FLUENT. Consequently, the fluctuating information of the Large Eddy Simulation was obtained. The two numerical models were established on the basis of the experiment data in the wind tunnel. The first was the flow field without any obstacles in which the fluctuating time-history data simulated by WAWS method were the inlet boundary condition. The second was the flow field with spires and roughness element model in accordance with those in the wind tunnel, in which the mean wind velocity was given to the inlet boundary and the wind turbulences were generated by the spires and roughness elements. The results show that fluctuating wind field based on the WAWS method can be the inlet boundary for the LES and provides a reference for the investigation of the inlet boundary conditions of LES.

Large Eddy Simulation；weighted amplitude wave superposition；numerical simulation；fluctuating wind field

2014-11-11

國家重點基礎研究計劃(973計劃)(2015CB057700)；國家自然科學基金資助項目(51408061，51178066，51278069)，National Natural Science Foundation of China(51408061，51178066，51278069) ；交通部西部交通科技項目重大科技專項資助項目(2011318824140)；長沙理工大學橋梁工程安全控制省部共建教育部重點實驗室開放基金資助項目(12KB01)；湖南省研究生創新項目(CX2014B375)

沈 煉(1988-)，男，湖南岳陽人，長沙理工大學博士研究生

?通訊聯系人，E-mail:ce_hanyan@163.com

1674-2974(2015)11-0064-08

O35

A