基于POD的大渦模擬入流脈動合成方法研究

鄭德乾， 顧 明， 張愛社

(1.河南工業大學 土木建筑學院，鄭州 450001； 2.同濟大學 土木工程防災國家重點實驗室，上海 200092； 3. 山東建筑大學 土木工程學院，濟南 250101)

基于空間平均的大渦模擬(Large Eddy Simulation, LES)方法可同時獲得結構表面的平均、脈動風荷載，且計算量適中，為應用前景廣闊的數值模擬方法。進行大氣邊界層內高層建筑繞流問題的大渦模擬分析，應采取合適的入流邊界條件生成方法，準確模擬紊流邊界層風場，以保證結構風荷載預測精度[1]。

目前的LES入流脈動合成方法可分為兩類：輔助域法[2-3]與直接合成法[4-9]。輔助域法為先在不含研究對象的計算域中進行LES計算，保存計算域某位置處脈動風速時程作為大渦模擬的入流邊界條件。雖輔助域法生成的入流脈動物理意義清晰，但合成過程較耗時。直接合成法則直接在入流面上直接合成具有指定統計特性及譜特性的脈動風速時程。一般通過在入流面的平均風剖面上疊加由蒙特卡羅法[4-5]或數字濾波法[6]生成的脈動風速時程序列合成LES入流脈動。與輔助域法相比，直接合成法更高效、便捷。本文采用直接合成法獲得LES入流脈動。

在LES計算中，為捕捉流場中盡量多的小尺度渦，致網格數目龐大，需采用并行計算技術提高計算效率。與串行計算不同，在流場并行計算中，不同的計算域網格分塊由不同計算節點負責，限制了多數LES入流脈動直接合成法在并行計算中的應用。解決辦法為在入流面各網格點獨立合成脈動風速時程[7-9]，但該類方法脈動風速合成需通過對三維能量譜離散實現，而風工程中脈動風速譜多為頻譜形式。故需發展既可通過離散頻譜形式風速譜合成脈動風速時程，又適用流體并行計算的LES入流脈動直接合成方法。

本文以Karman譜為目標譜，采用特征正交分解(Proper Orthogonal Decomposition, POD)型譜表示法[10]合成脈動風速時程。針對流體并行計算特點，發展適用于流體并行計算的LES入流脈動合成方法，以提高LES方法進行紊流風場內結構風荷載的預測精度。基于Fluent軟件平臺，分別通過空計算域紊流邊界層風場及寬高比1∶6的單體方形截面高層建筑非定常繞流的LES計算結果與剛性模型測壓試驗[11]結果比較，驗證本文方法的有效性。

1 風洞試驗

某高層建筑寬高比1∶6，見圖1(a)，B=D=0.1 m，H=0.6 m，其剛性模型同步測壓風洞試驗在同濟大學土木工程防災國家重點實驗室風洞試驗室TJ-2大氣邊界層風洞中進行[11]，模型縮尺比1∶500。據該高層建筑所處地面粗糙度類別，按1/500縮尺比模擬B類風場，試驗所測平均風剖面與湍流度剖面見圖1(b)、(c)。試驗時將模型置于轉盤中心，通過轉盤模擬不同風向角。限于篇幅，本文主要對圖1(a)中0°風向角(來流垂直吹向結構)進行大渦模擬計算。風洞試驗中結構表面共布置200測點，此處僅給出用于本文數值模擬與風洞試驗比較的模型表面主要測點(圖1(a))。詳細風洞試驗結果見文獻[11]。

2 LES入流脈動合成

LES計算中入流面網格數量龐大，本文基于POD型譜表示法[10]合成入流面主要網格點的脈動風速時程，采用有限元形函數進行空間插值得出入流面所有網格點的脈動風速時程。據流體并行計算各節點調用關系及數據存儲特點，編程實現流體并行計算時已合成脈動風速時程序列的讀入、賦值。并給出LES入流脈動合成方法的基本原理與實施步驟。

圖1 風洞試驗測點布置及試驗風剖面

2.1 POD型譜表示法

POD特征正交分解為分析復雜隨機場統計特征的有效工具，該方法通過建立新坐標空間，采用時空分離技術簡化隨機現象描述。文獻[10]將POD法引入橋梁隨機風場模擬中發展了POD型譜表示法。該方法物理意義明確，可通過模態截斷提高計算效率。

對一維N變量零均值的平穩隨機過程，其功率譜密度矩陣可表示為：

(1)

其中：ω為脈動風速圓頻率；Sij(ω)(i,j=1,2,3,…,N)，i=j時Sij(ω)為脈動風速自功率譜密度函數；i≠j時Sij(ω)為脈動風速互功率譜密度函數。

為模擬獲得具有指定譜密度函數形式的脈動風速時程，需對式(1)譜密度矩陣[S(ω)]進行分解。采用原型譜表示法(如諧波合成法)時，對譜密度矩陣[S(ω)]一般采用Cholesky分解，而POD型譜法則采用POD方法分解[10]，即：

(2)

式中：[φ1(ω)], [φ2(ω)],…, [φN(ω)]與η1(ω),η2(ω),…,ηN(ω)分別為譜密度矩陣[S(ω)]的N階特征向量、特征值，特征向量按特征值降序排列，滿足正交化條件[10]：

(3)

式中：δ為Dirac-δ函數。譜矩陣特征值代表對應振型包含振動能量大小，譜振型則由頻域描述隨機過程包含能量的空間分布形式，具有特定物理意義。若前NS階譜振型所含能量足夠大，則可引入模態截斷技術進行模型縮減[10]：

(4)

采用POD型譜表示法時，某點脈動風速時程表達式[10]為：

(5)

式中：θkl為在[0,2π]間均勻分布的獨立隨機相位角；Δω為頻率步長；ωl= (l-1)Δω為頻率采樣點序列。

誤差分析[10]表明，POD型譜表示法總體相對隨機誤差與根方差相對隨機誤差小于原型譜表示法。本文用式(5)的POD型譜表示法合成具有指定頻譜特性的脈動風速時程。

2.2 LES入流脈動合成參數選擇

2.2.1 平均風剖面

用對數律表達式[12]描述平均風速沿高度變化：

(6)

式中：κ=0.42為馮·卡曼常數；u*為摩擦速度；z0為粗糙長度，由圖1(b)試驗平均風剖面擬合獲得，見表1。圖1(b)比較結果顯示平均風剖面擬合值與試驗值吻合較好。

表1 B類1∶500縮尺比風場參數擬合值

2.2.2 湍流度剖面

數值計算時，湍流模型中湍流黏度系數與紊流強度成正比，近地面處湍流強度過大易致湍流黏度系數過高而影響數值計算收斂，因此風洞試驗風場的湍流度剖面采用分段函數進行擬合：

(7)

式中：a，b為常數，由圖1(c)試驗湍流度剖面擬合獲得(表1)。湍流度剖面擬合值與試驗值比較見圖1(c)，由圖1看出，二者一致性較好。

2.2.3 縱向湍流積分尺度

本文缺少相應縱向湍流積分尺度隨高度變化的試驗值，而基于不同湍流積分尺度建議公式值有所不同。縱向湍流積分尺度采用日本AIJ規范[13]建議表達式。為便于同風洞試驗結果比較，本文模擬與風洞試驗相對縮尺比1∶500風場，修正后縱向湍流積分尺度計算式為：

(8)

式中：λL為模型幾何縮尺比，λL=500。

2.2.4 脈動風速譜

Karman譜能反映脈動風速沿高度變化，本文以其為目標風速譜，Karman譜表達式變化形式為：

(9)

2.2.5 脈動風速相干函數

縱向脈動風速豎向、橫向相干函數表達式為[14]：

(10)

式中：r為空間兩點距離；y1，y2，z1，z2分別為空間兩點橫、豎坐標；U(z1)，U(z2)分別為高度z1，z2處平均風速；Cy，Cz分別為橫、豎向相干指數衰減系數，與平均風速、離地高度、地面粗糙度有關，取Cy=16，Cz=10。

2.3 入流面網格點脈動風速合成步驟

為盡量捕捉到較小渦運動，LES計算中需較密網格布置，使入流面網格點數增大，若同時生成入流面所有網格點脈動風速時程耗時耗力。故采取入流面網格點分組循環[4]，或分若干塊并假設各分塊內網格點風速完全相關[5]等措施生成入流面脈動風速。前者實施較麻煩，后者會造成所得入流面上速度場不連續。本文采用有限元形函數，空間插值出入流面所有網格點的LES入流脈動，盡量保證速度場的連續性。考慮入流面網格多為規則矩形，故采用矩形等參單元進行形函數計算、插值，實施步驟為：

(1)在LES入流面上生成一定數目的點坐標，且應含入流面邊界上4個頂點，并在近地面及模型高度附近區域較密分布，見圖2(a)中“●”；

(2)據點空間位置關系，在步驟(1)生成點分別組成的矩形有限單元內，選所含相應的LES入流面網格點(圖2(b))，分別計算網格點在所屬矩形單元的形函數；

(3)用式(5)POD型譜表示法同時生成步驟(1)點的脈動風速時程，即獲得步驟(2)中各矩形單元4頂點的脈動風速時程；

(4)據步驟(2)所得入流面網格點在所屬單元的形函數及步驟(3)生成矩形單元4頂點的脈動風速時程，依次插值出最終LES入流面所有網格點的脈動風速時程。

圖2 LES入流面網格點脈動風速空間插值示意圖

除空間插值問題外，LES計算中較小時間步長亦會帶來脈動風速時程時間插值問題。本文通過用POD型譜表示法直接生成與LES計算相同時間步長的脈動風速時間序列避免該問題。

2.4 并行計算中入流面網格點脈動風速時程讀入

與串行計算相比，流體采用并行計算技術時，需對計算域網格進行分區，不同網格分區由不同計算節點負責，并通過不同區間的重疊網格實現計算節點間數據傳輸，見圖3(a)。并行計算存在主進程及若干計算節點等多個進程，見圖3(b)，其中主進程不含任何網格信息，主要負責外部數據存儲文件的開、關及進行命令解釋與傳遞；0#計算節點負責命令分發及數據收集[15]。因此，LES并行計算時，需針對并行計算特點采用相應措施確保合成的脈動風速時程準確讀入與賦值。本文先由主進程讀取外部存儲的脈動風速時程數據文件，傳遞給0#計算節點；0#計算節點將脈動風速時程數據分發給所有計算節點；各計算節點將各自脈動風速數據分別賦值給負責的入流面分區網格點。通過合成的脈動風速時程數據分時刻存儲及LES計算中分時刻讀入網格點風速值方法提高LES計算效率；通過將LES非定常迭代時間步引入風速數據存儲矩陣辦法解決LES計算因故中斷退出后繼續計算時入流面網格點風速值的正確讀入問題。以上處理方法，本文用Fluent UDF編程實現。

圖3 Fluent并行計算網格分區及命令傳遞示意圖[15]

3 單體方形截面高層建筑大渦模擬

基于Fluent軟件平臺，采用與風洞試驗[11]相同比例縮尺比，進行B類1:500縮尺風場方形截面高層建筑的LES非定常繞流計算。計算域94D×28D× 36D，見圖4(a)。采用非均勻結構化網格對計算域進行離散，近壁面處網格加密，見圖4(b)、(c)，最小網格尺度1/2000D，對應壁面y+值小于5.0。采用網格總數分別為151萬(G1)、90萬(G2)兩種劃分方法。采用速度入口邊界條件，為考察本文方法的有效性及精度，LES入流脈動合成分別采用本文方法及Fluent內置改進合成法[9]。計算域其它位置邊界條件設置見圖4(a)。

表2為本文LES計算工況，壓力速度耦合均采用SIMPLEC算法，空間離散采用具有二階精度的bounded central differencing格式；時間離散采用二階全隱格式，時間步長0.0005 s；亞格子模型采用Dynamic Smagorinsky-Lilly模型[16]。

表2 單體方形截面高層建筑LES計算工況

圖4 LES計算域、邊界條件設置及離散網格示意圖

4 結果與討論

圖5為本文考慮工況下，數值模擬所得距離模型前方5D處的模擬風剖面、速度監測點脈動風速譜，與試驗[11]風剖面及Karman譜的比較。由圖5看出，①所有工況下，數值模擬所得平均風剖面與風洞試驗結果基本一致；②數值模擬所得紊流度剖面在遠離地面處與試驗結果吻合相對較好，而在近地面處均有所衰減，相同網格密度下與結果[9](G1_LES1工況)相比，采用本文方法合成LES入流脈動(G1_LES2工況)，數值模擬所得紊流度剖面近地面衰減有所改善；③數值模擬所得脈動風速譜在低頻處均稍高于Karman譜而高頻處產生衰減，可能因LES計算網格分辨率不足以捕捉到更小尺度的渦(高頻渦)所致。采用較密網格布置及本文方法合成LES入流脈動(G1_LES2工況)時，脈動風速譜在高頻處的衰減有所減弱，與Karman譜一致性更好。本文LES數值模擬計算所得脈動風速譜基本能捕捉到無量綱頻率在0.2～1.0范圍內工程中所需高頻脈動風速譜段即慣性子區段。

圖5 LES模擬風場與試驗結果的比較

圖6為結構表面2/3H高度處測點平均和根方差風壓系數的LES計算結果與風洞試驗[11]結果比較，圖中CP,mean，CP,rms分別為測點平均與根方差風壓系數，此處及下文中結構表面風壓系數均用模型高度處來流風速無量綱化。由圖6看出，①所有工況下，LES計算所得平均風壓系數均與風洞試驗結果基本吻合，但在結構背風面LES計算結果均較高估計風吸力，此因結構背風面處于流動分離后的尾流區，流動相對復雜，氣流分離會使風吸力較大；②根方差風壓系數雖與風洞試驗結果有偏差，但趨勢基本一致，相同網格密度下，與結果[9](G1_LES1工況)相比，用本文方法合成LES入流脈動(G1_LES2工況)結果稍高地估計結構表面風壓脈動；采用較稀疏網格布置(G2_LES2工況)過高估計側面風壓脈動值，此因較稀疏網格布置不足以準確捕捉結構側面的流動分離現象。

圖6 LES模擬所得高層建筑表面2/3H高度處風壓系數統計值與試驗結果比較

圖7 LES模擬所得高層建筑表面2/3H高度處風壓系數自譜與試驗結果比較

圖8 LES模擬所得高層建筑表面2/3H高度處風壓系數水平相干性與試驗結果比較

圖7、圖8分別為數值模擬所得結構表面2/3H高度測點風壓系數自譜和水平相干性與風洞試驗[11]結果比較。由圖7看出，①在結構迎風面，LES計算所得風壓系數自譜在高頻處衰減較快，此因結構迎風面風壓直接受來流紊流影響，模擬脈動風速譜均在高頻處衰減所致；②在結構側、背風面，LES計算所得風壓系數自譜與風洞試驗結果差別較小，結構側面LES計算及風洞試驗所得橫風向旋渦脫落無量綱頻率均為0.1，LES計算結果在高頻處所含頻率成分較少，此因網格過濾效應作用所致。由圖8看出，LES計算所得結構測點風壓系數的水平相干性隨頻率變化較平緩光滑，在低頻處與風洞試驗結果吻合較好，而在高頻處與風洞試驗結果均值線較接近，LES結果基本反映出測點風壓系數水平相干性隨頻率變化規律。

5 結 論

針對流體并行計算特點，基于POD型譜表示法，發展適于流體并行計算的LES入流脈動合成方法。基于Fluent軟件平臺，進行寬高比1∶6的單體方形截面高層建筑非定常繞流LES計算，通過數值模擬所得風場特性、風壓系數統計值及譜特性，與風洞試驗及文獻數值模擬結果比較，驗證本文方法的有效性，結論如下：

(1)對紊流風場，本文方法可較好模擬平均風剖面及紊流度剖面，一定程度上能抑制近地面處紊流度及模擬脈動風速譜高頻處的衰減現象。

(2)對結構表面風壓預測，采用較密網格布置及本文方法合成LES入流脈動，不僅能較好進行結構表面風壓分布預測，且能較準確給出風壓系數自譜與水平相干性結果。

參 考 文 獻

[1]Tamura T. Towards practical use of LES in wind engineering [J].Journal of Wind Engineering and Industrial Aerodynamics, 2008, 96(10-11): 1451-1471.

[2]Lund T S, Wu X, Squires K D. Generation of turbulent inflow data for spatially-developing boundary layer simulations[J].Journal of Computational Physics,1998, 140(2): 233-258.

[3]朱偉亮, 楊慶山. 湍流邊界層中低矮建筑繞流大渦模擬[J]. 建筑結構學報, 2010, 31(10): 41-47.

ZHU Wei-liang, YANG Qing-shan. Large eddy simulation of flow around a low-rise building immersed in turbulent boundary layer[J]. Journal of Building Structures,2010, 31(10): 41-47.

[4]Kondo K, Murakami S, Mochida A. Generation of velocity fluctuations for inflow boundary condition of LES[J]. Journal of Wind Engineering and Industrial Aerodynamics, 1997, 67-68: 51-64.

[5]李清雅, 葉繼紅. 三維空間曲面結構風荷載的數值模擬[J]. 振動與沖擊, 2009, 28(4): 121-126.

LI Qing-ya, YE Ji-hong. Numerical simulation of wind pressure on spatial structures[J]. Journal of Vabraion and Shock, 2009, 28(4): 121-126.

[6]Klein M, Sadiki A, Janicka J. A digital filter based generation of inflow data for spatially developing direct numerical or large eddy simulations[J].Journal of Computational Physics, 2003, 186(2): 652-665.

[7]Smirnov A, Shi S, Celik I. Random flow generation technique for large eddy simulations and particle-dynamics modeling [J]. Journal of Fluids Engineering, 2001, 123(2): 359-371.

[8]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 Aerodynamics, 2010, 98(10-11): 600-617.

[9]ZHENG De-qian, ZHANG Ai-she, GU Ming. Improvement of inflow boundary condition in Large Eddy Simulation of flow around a tall building [J]. Engineering Applications of Computational Fluid Mechanics, 2012, 6(4): 633-647.

[10]胡 亮. 基于特征正交分解的橋梁風場隨機模擬[D]. 武漢：華中科技大學, 2007.

[11]唐 意.高層建筑彎扭耦合風致振動及靜力等效風荷載研究 [D].上海：同濟大學, 2006.

[12]Richards P J, Hoxey R P Appropriate boundary conditions for computational wind engineering models using the k-ε turbulence model[J]. Journal of Wind Engineering and Industrial Aerodynamics, 1993, 46-47: 145-153.

[13]Architectural Institute of Japan (AIJ), Recommendations for loads on buildings [S]. Japan, 2004.

[14]Davenport A G. The dependence of wind load upon meteorological parameters[C]. Proceedings of the international research seminar on wind effects on building and structures, University of Toronto Press, Toronto, 1968: 19-82.

[15]Fluent 6 documentation [M]. Fluent Inc, 2005.

[16]Lilly D K. A proposed modification of the Germano subgrid-scale closure model [J]. Physics of Fluids,1992, 4: 633-635.