基于諧波合成法的大渦模擬脈動風(fēng)場生成方法研究

沈煉 韓艷 蔡春聲 董國朝

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

關(guān)鍵詞：大渦模擬；諧波合成；數(shù)值模擬；脈動風(fēng)場

中圖分類號：O35 文獻(xiàn)標(biāo)識碼：A

文章編號：1674-2974（2015）11-0064-08

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

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

本文借助于經(jīng)典的諧波合成方法，在充分考慮湍流的平均風(fēng)剖面、目標(biāo)譜、空間相關(guān)性等指標(biāo)下，合成一系列的隨機(jī)序列時(shí)程數(shù)據(jù)，把生成的時(shí)程數(shù)據(jù)與大渦模擬的入口建立對應(yīng)關(guān)系，基于對FLUENT軟件平臺進(jìn)行二次開發(fā)，把生成的隨機(jī)序列數(shù)賦給大渦模擬的入口邊界，從而生成入口脈動信息.最后通過對比數(shù)值模型和風(fēng)洞試驗(yàn)結(jié)果發(fā)現(xiàn)：

基于諧波合成的脈動信息生成方法能滿足大渦模擬入口邊界的要求，其優(yōu)勢在于目標(biāo)譜及相關(guān)性參數(shù)能根據(jù)不同情況的需要進(jìn)行隨意調(diào)整，合成的時(shí)程數(shù)據(jù)可直接賦給入口邊界，生成速度快，模擬時(shí)間長短可以隨意控制，模擬的點(diǎn)都相互獨(dú)立，可以很好地進(jìn)行并行計(jì)算，大大提高計(jì)算效率.

1 數(shù)值模擬理論與方法

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

2 大渦模擬入口脈動風(fēng)場特性分析

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

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

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

對流場入口中心0.5 m和1 m高度處的速度時(shí)程數(shù)據(jù)進(jìn)行相關(guān)性分析，定義兩點(diǎn)分別為A和B.圖4給出了入口監(jiān)測點(diǎn)A和B的互相關(guān)性和B點(diǎn)的自相關(guān)性檢驗(yàn)結(jié)果，圖4（a）中目標(biāo)值由公式（11）給出.

圖4表明了入口處脈動風(fēng)場相關(guān)性呈指數(shù)率衰減，模擬值與理論目標(biāo)值衰減趨勢一致，體現(xiàn)出了模擬風(fēng)場在入口處具有良好的相關(guān)性.

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

3 數(shù)值模型

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

3.1 模型尺寸與網(wǎng)格數(shù)量

3.1.1 空風(fēng)洞數(shù)值模型

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

3.1.2 尖劈粗糙元數(shù)值模型

尖劈粗糙元數(shù)值模型與物理風(fēng)洞試驗(yàn)條件保持一致，布置如圖6所示.與上一節(jié)類似，為保證精度，對尖劈粗糙元數(shù)值模型進(jìn)行全場六面體網(wǎng)格劃分，并對260萬，450萬，800萬3種網(wǎng)格進(jìn)行了網(wǎng)格無關(guān)性測試，無關(guān)性測試結(jié)果如表1所示.

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

3.2 邊界條件設(shè)置

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

4 計(jì)算結(jié)果

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

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

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

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

在計(jì)算效率方面，將兩種數(shù)值模型所需的網(wǎng)格數(shù)和耗時(shí)量作出對比，如表3所示.

由表3可以發(fā)現(xiàn)，在達(dá)到允許精度要求下，空風(fēng)洞模型所需的網(wǎng)格數(shù)不到尖劈粗糙元模型的一半，計(jì)算耗時(shí)量方面也遠(yuǎn)遠(yuǎn)優(yōu)于尖劈粗糙元模型.

此外，空風(fēng)洞數(shù)值模型具有更好的適應(yīng)性，特別是對不同脈動風(fēng)場的模擬，空風(fēng)洞模型只需對目標(biāo)風(fēng)場的特性參數(shù)進(jìn)行調(diào)整，操作簡單方便，可調(diào)性強(qiáng).而尖劈粗糙元數(shù)值模型只能通過不斷調(diào)整尖劈粗糙元的位置和數(shù)量來對目標(biāo)風(fēng)場進(jìn)行試算，計(jì)算工作量大，消耗時(shí)間長.

5 結(jié) 論

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

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

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

3）通過對比空風(fēng)洞數(shù)值模型和尖劈粗糙元模型的計(jì)算效率發(fā)現(xiàn)，空風(fēng)洞數(shù)值模型在網(wǎng)格數(shù)量和計(jì)算耗時(shí)方面均要遠(yuǎn)遠(yuǎn)優(yōu)于尖劈粗糙元數(shù)值模型.

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

基于諧波合成法的脈動風(fēng)場生成方法是譜合成方法的一種拓展，擁有許多優(yōu)點(diǎn)的同時(shí)也存在一些不足，如生成的隨機(jī)風(fēng)速時(shí)程不滿足N-S方程的連續(xù)性，風(fēng)譜高頻衰減等問題還有待進(jìn)一步解決.

參考文獻(xiàn)

[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 模型的近地脈動風(fēng)場數(shù)值模擬[J]. 工程力學(xué)，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， et al. 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， et al. 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， et al. 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]盧春玲，李秋勝，黃生洪，等. 大跨度屋蓋風(fēng)荷載的大渦模擬研究[J]. 湖南大學(xué)學(xué)報(bào)：自然科學(xué)版，2010，37（10）：7-12.

LU Chun-ling， LI Qiu-sheng， HUANG Sheng-hong， et al. 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]. 湖南大學(xué)學(xué)報(bào)：自然科學(xué)版，2013，40（11）： 26-33.