基于LES和DES的定日鏡結構風致響應分析

盧春玲， 陳建通， 陳錦焜， 王 強

(1. 桂林理工大學 土木與建筑工程學院， 桂林 541004； 2. 廣西建筑新能源與節能重點實驗室， 桂林 541004；3. 桂林理工大學 有色金屬礦產勘查與資源高效利用協同創新中心， 桂林 541004)

塔式太陽能發電系統需要在開闊場地內大范圍安裝定日鏡，定日鏡是其主要設備。當定日鏡離中央集熱塔較遠時，微小的風致振動也會降低定日鏡的聚光效率，造成效率損失。此外，在開闊的野外場地上，結構高度低，易受近地面復雜風場影響。并且定日鏡受風面積大、剛度低，在風力動力荷載下容易引起定日鏡振動、損壞，甚至倒塌。所以，對定日鏡的風致響應研究具有重要意義。

針對定日鏡的抗風問題，國內外學者進行了大量的研究。Vasquez-Arango等[1]結合實測與有限元分析，對定日鏡的模態進行了驗證，為定日鏡有限元模擬提供了依據。Terrés-Nícoli等[2]通過風洞試驗，對定日鏡結構在風荷載作用下的應力分布進行了分析研究。Wolmarans等[3]采用單向流固耦合方法，研究了中型定日鏡周圍渦流脫落過程和瞬態風荷載。Gong等[4]在風洞試驗和數值模擬的基礎上研究了定日鏡的脈動風壓特征和風致動態響應。王延忠等[5]采用流固耦合數值模擬方法，研究了定日鏡在不同仰角下的風致響應。馮煜等[6]采用AR模型模擬風場時程，對三維風場下定日鏡風致動態響應進行了分析。王鶯歌[7]將理論分析、風洞試驗及數值模擬相結合，對定日鏡結構的抗風問題進行了綜合性的研究。相較于數值模擬，現場實測和風洞試驗研究獲取的流場信息相對較少，且數據后處理復雜、試驗費用昂貴。對定日鏡的數值模擬研究，雷諾平均法(Reynolds-averaged Navier-Stokes, RANS)較為常用，但RANS模型對脈動過程進行時均化處理，這對風致響應分析非常不利。而LES和DES避免了RANS抹平細節的缺陷，可以較準確的模擬脈動風壓及風荷載。但其對網格精度及入口湍流的準確度要求較高，利用LES和DES對定日鏡結構進行結構風荷載和風致響應的研究還較少。

本文選用LES和DES，結合離散再合成的隨機湍流生成法對定日鏡在0度風向角下0°、30°、60°鏡面仰角(鏡面法線與水平夾角)的多種典型工況[8]下進行數值模擬，得到定日鏡正反面風荷載時程數據。建立有限元模型，將表面風荷載時程施加在定日鏡有限元模型上，對其進行風致響應分析。將計算結果與風洞試驗的相應數據進行對比，以驗證LES和DES的準確性。

1 風洞試驗

定日鏡風洞測力試驗在湖南大學風洞實驗室的HD-2單回流邊界層風洞中的低速試驗段內進行。該風洞低速試驗段尺寸為寬5.5 m，高4.4 m，長15 m。采用有機玻璃制成定日鏡剛性模型，縮比為1/10，高1.2 m，寬0.9 m，模擬出鏡面、轉動軸、支撐臂組件及立柱等定日鏡主要組成構件，使模型與實體結構在幾何外形上保持相似，實體結構及模型如圖1。

(a) 定日鏡實體結構

風洞試驗根據定日鏡實體結構所處地區的地貌特征，采用格柵、尖梯、擋板、粗糙元裝置模擬出我國建筑結構荷載規范(GB 50009—2001)規定的B類風場，采用IFA300熱膜風速儀對大氣邊界層模擬風場進行調試和測定。風洞試驗參考點選在 1.0 m高度處，對應于實際中 10 m 高度，采用皮托管測控流場的參考風速，參考風速為8 m /s。定日鏡模型的基座處固定有六分力天平，從而獲得模型六個方向上的測力數據。在定日鏡正面和背面各布置144個測點，測點布置如圖2，通過64通道DTCnet壓力掃描閥、ESP-64HD型傳感器及INV306大容量自動數據采集和信號處理系統測量各風壓時程。采樣頻率330 Hz；每個測點記錄20 000個數據長度，總采樣時間為60.6 s。

(a) 迎風面測點布置

(b) 背風面測點布置圖2 風洞試驗模型風壓測點布置Fig.2 Distributions of pressure measuring points in the model of wind tunnel test

2 數值模擬方法

2.1 計算模型及網格劃分

定日鏡結構的計算模型為全尺寸，定日鏡整體高度為(H)11.93 m，鏡面寬度(L)9 m，轉動軸和立柱為直徑1 m的圓柱，鏡面與支架之間留有空隙。計算域z方向高為5H,Y方向寬度為8L,鏡面法向(X)方向長度為20L,流域入口距模型鏡面5L，其結構阻塞度為2.6%，計算域尺寸如圖3。

圖3 計算域尺寸Fig.3 Computational domain size

流場網格設置采用混合網格。圖4為定日鏡表面網格劃分方案，圖5為流場網格劃分方案。定日鏡結構處于X方向長3L，Y方向寬2L，Z方向高1.5H的網格加密區中。加密區采用非結構化網格，網格尺寸為0.01～0.5 m，在定日鏡結構表面的網格最小，并以5%的增長率逐漸增大。加密區嵌套于非加密區之中，非加密區采用結構化網格，總網格數約為510萬。這樣既能保證計算效率又能保證定日鏡周圍的流場計算精度。

圖4 定日鏡表面網格劃分Fig.4 Mesh for heliostat surface

圖5 流場整體網格劃分Fig.5 Global mesh generation for flow field

2.2 湍流模型

本文采用LES和DES對定日鏡結構進行數值模擬計算。對于LES，FLUENT提供了四種亞格子模型，分別是Smagorinsky-Lilly模型、動態Smagorinsky-Lilly模型、壁面自適應局部渦流黏性模型和動態動能亞網格模型。動態動能亞網格模型通過亞網格尺度湍流動能的傳輸來描述亞網格尺度湍流。彌補了Smagorinsky-Lilly模型在高雷諾數下，能量耗散難符合預期的缺點。因此，采用了動態亞網格動能模型。其應力尺度為

(1)

其中，ksgs通過求解其輸運方程(2)得到

(2)

DES中，采用分離渦的一種優化模型——延遲分離渦(delayed detached eddy simulation, DDES)，分離渦模型是Spalart等[9]提出的一種混合數值模型，該方法結合了RANS和LES，通過對網格尺寸相關的長度比例l的判斷，自動在兩種方法之間轉換。如式(3)、(4)所示

l=min(d,CDESΔDES)

(3)

ΔDES=max(δx，δy，δz)

(4)

式中:CDES=0.65，δx、δy、δz表示三個方向的網格比例，d表示網格到壁面的距離。當l=d時，啟用RANS計算；當網格從壁面開始增長時，當l=CDESΔDES時，轉換為LES進行計算。某些情況下，網格縱橫比在近壁面不夠大，即使網格厚度很小，LES也會被啟用。為了保證計算精度和效率，DDES對d進行修正如式(5)、(6)

d≡d-fdmax(0,d-CDESΔ)

(5)

fd=1-tanh[(8rd)3]

(6)

(7)

式中：rd是湍流尺度與網格到壁面的距離之比，其中vt是動力黏度，v是分子黏度，Uij是速度梯度，k是Von Karman常數。當網格靠近壁面時fd=0，啟用RANS計算。當網格遠離壁面時rd=1，即fd=1，湍流模型轉化為LES。通過對d的修正，DDES延遲了LES的發生，提高模擬的準確度和計算效率。DES亞網格模型沿用動態動能亞網格模型，RANS模型設定為Realizablek-ε兩方程模型。

2.3 求解設置

在數值求解中，采用不可壓縮流。在LES和DES計算前，導入瞬態化處理的RANS模型計算結果作為初始流場，以加快計算收斂速度。數值計算采用SIMPLE算法，時間離散采用二階隱式格式，空間離散采用二階中心格式，計算收斂殘差標準取為1×10-4。進行流體計算時，采用Fluent提供的monitor模塊對定日鏡結構上的每個鏡面進行風壓監測并對其積分，得到風荷載時程數據。監測面布置如圖6。

圖6 數值模擬模型監測面布置Fig.6 Pressure monitoring surface on the model for the numerical simulation

模擬計算在高性能工作站上進行，工作站搭載18核36線程的Intel Xeon Gold 6140處理器，運行內存192 GB。LES與DES計算時間步長取0.05 s，每個時間步長最大迭代次數為20次，每個工況進行了10 000個時間步長的計算。LES模擬每個時間步長計算耗時約40 s，每個工況耗時約110 h；DES模擬每個時間步長計算耗時約35 s，每個工況耗時約95 h。

2.4 邊界條件

本文LES和DES模擬入口湍流生成方法采用DSRFG法[10]，相對于以往的方法，DSRFG法具有以下優點：該方法基于嚴格的理論推導，能產生滿足任意形式功率譜及各向異性的湍流脈動風速場，能夠嚴格滿足流體連續性條件，從而保證了大渦模擬計算的穩定性，易于并行化處理等。采用DSRFG方法生成的入口瞬時風速分布云圖，如圖7所示。

圖7 DSRFG法產生的入口瞬時速度場Fig.7 Instantaneous velocity field at the inflow boundary generated by DSRFG method

計算域的入口處邊界條件設定為速度入口，按照定日鏡實際工作環境，大氣邊界層平均風速度剖面的模擬采用如式(8)的指數率形式

(8)

(a) 平均風速剖面圖

(b) 湍流強度剖面圖圖8 數值模擬與風洞試驗入流條件對比圖Fig.8 Inflow boundary condition of numerical simulation and wind tunnel test

湍流強度參考日本規范[11]中的第II類地貌，取值如式(9)

(9)

式中：I0為由規范確定的湍流強度；zg為梯度風高度，m；B類風場下，I0=0.23，zg=350 m。數值模擬與風洞試驗湍流強度剖面圖如圖8(b)。

2.5 有限元模型

通過ANSYS軟件選用不同的材料單元模擬了定日鏡主要組成構件，建立了定日鏡有限元模型。采用shell63單元模擬鏡面和轉動軸、立柱；采用beam188單元模擬出不同型號空心鋼管構成的后部支撐鋼架；利用耦合的方法，模擬出支撐鋼架與鏡面和鏡面與鏡面之間的駁接抓連接。通過命令流，可快速實現有限元模型不同工況的轉換。圖9給出了鏡面仰角為0°的有限元模型示意圖。實測頻率由業主方提供，實測試驗在24塊鏡面中點、立柱中點和頂端各布置一個加速度傳感器，采用環境激勵，通過模態分析軟件將采集數據進行處理分析，得到定日鏡實體結構的模態參數。對有限元模型進行模態分析，將有限元模型前三階自振頻率與實測頻率進行對比，結果如表1。將Fluent流體計算得到的定日鏡表面風荷載時程加載到有限元模型上進行瞬態動力分析，得到定日鏡各測點的位移響應，位移測點的布置如圖10。

圖9 有限元模型Fig.9 Finite element model of heliostat

圖10 位移測點布置Fig.10 Displacement monitoring points on the FEM

表1 不同鏡面仰角下有限元模型與實體結構的自振頻率Tab.1 Natural frequency of heliostat finite element model and Entity structure different elevation angle Hz

3 結果分析

3.1 流場分析

風洞試驗得到的流場信息十分有限，而數值模擬方法可以得到豐富的可視化流場信息，可以更好探究流場作用機理。

圖11為兩種模擬方法的各工況下在Y=0、Z=0.15H和Z=0.85H平面上的平均風速流線圖。從圖中可以看出，兩種模擬方法都表現出相似的流場特性，這反映了定日鏡結構在不同鏡面仰角下流場的一般規律。氣流在定日鏡邊緣發生流動分離，渦旋作用于定日鏡背面形成負壓區。隨著仰角增大，定日鏡后方的兩側尾流區域變長且集中。這在定日鏡群中，可能會對周邊定日鏡產生不利影響。

(a) LES-00

(b) DES-00

(c) LES-30

(d) DES-30

(e) LES-60

(f) DES-60圖11 各工況下平均風速流線圖Fig.11 Mean veloicty contours around heliostat with different elevation angle

圖12為兩種模擬方法的各工況下Y=3處渦量圖。0°和30°時，氣流沖擊定日鏡后，氣流從流場駐點處分別向定日鏡頂端和底端流去，部分氣流在鏡面頂端發生分離，渦旋向上發展，再附著現象不明顯；定日鏡底部也相似的現象。隨著仰角增大，60°時，氣流在底部分離，發生強烈的渦旋脫落，形成眾多強度較大的小渦后再附著在鏡面的中上部。與DES相比，LES在結構近壁面形成的渦旋尺度更加細碎。

(a) LES-00

(b) DES-00

(c) LES-30

(d) DES-30

(e) LES-60

(f) DES-60圖12 各工況下Y=3 m平面渦量圖Fig.12 Vorticity magnitude contours on X-Z plane at Y=3 m with different elevation angle

3.2 風壓相關系數

在對建筑結構風致響應的考慮中，脈動風壓的空間相關性起這重要的作用。兩點之間的脈動風壓相關系數計算式為

(10)

式中：Cov(X,Y)表示兩點脈動風壓協方差；Var[X]和Var[Y]分別表示兩點脈動風壓方差。

以定日鏡中間B列測點為代表進行比較，選取了8個測點來研究脈動風壓隨高度變化的規律，分析定日鏡表面脈動風壓豎向相關性。不同工況的脈動風壓相關系數如圖13。

圖13 B列各測點豎向相關性Fig.13 Vertical correlation of column B

迎風面脈動風壓的豎向相關系數隨著測點間距離的增大而減小。背風面的相關性要強于迎風面，背風面測點的相關性均在0.9以上。隨著鏡面仰角的增大，迎風面相關系數逐漸增大，背風面相關系數逐漸減小。對比發現，DES結果的相關性均強于LES的結果。

3.3 風荷載

風荷載時程具有隨機性，對風荷載時程進行譜分析可以較直觀的看出風荷載的特性。圖14為定日鏡0°工況鏡面總體、頂層和底層的歸一化風荷載功率譜的對比圖。定日鏡屬于低矮結構，其風荷載譜受特征湍流影響較大，即定日鏡表面渦旋脫落的影響。從圖中可以看出，在低頻段時，風荷載功率譜與卡門譜有吻合度較好，而風荷載功率譜峰值出現0.7 Hz， 說明定日鏡頂部與底部的渦旋脫離頻率集中在0.7 Hz附近的低頻段。風荷載功率譜與卡門在1 Hz附近出現分離，并且衰減較快。這在Yan等[12]和Daniels等[13]的研究也有相似的情況，這是由于數值模擬方法有限的網格分辨率造成的[14]。

圖14 順風向荷載時程及功率譜Fig.14 Time history and power spectrum of along-wind load

3.4 等效靜力風荷載

為驗證數值風洞模擬的準確性，對模擬結果進行時域分析和頻域分析計算出等效靜力風荷載，并與風洞試驗測力結果對比，以驗證數值模擬結果的準確性。等效靜風荷載的計算方法有多種，本文分別采用慣性風荷載法、慣性力-荷載響應相關法[15]對定日鏡的等效靜力風荷載進行時域分析和頻域分析。

慣性風荷載方法中等效靜力風荷載可表示為

(11)

慣性力-荷載響應相關法采用定日鏡前三階振型對應的自振頻率進行計算，通過編寫Matlab程序，對等效靜力風荷載進行求解計算。

表2為各工況的等效風荷載計算結果。時域法和頻域法結果比較接近，但和風洞結果相比都較小。這是由于數值模擬測點設置有限，只考慮了作用于定日鏡鏡面上的風荷載，忽略了支撐臂、轉動軸和立柱的風荷載導致的。頻域分析結果小于時域分析，這是由于頻域法只考慮了三階振型。LES和DES的結果差距較小，表面這兩種湍流模型對等效風荷載的預測有較好的一致性。

表2 定日鏡順風向等效風荷載數值模擬與風洞試驗對比Tab.2 Comparison between numerical simulation and wind tunnel test of equivalent along-wind load of heliostat N

3.5 位移響應分析

各工況不同結點的位移時程趨勢基本相同，圖15為LES0°鏡面仰角下的各測點的位移時程。從圖中可以看出，定日鏡底部的鏡面位移最小，頂部鏡面的位移最大。

圖15 各測點位移時程Fig.15 Displacement time history of each monitoring point

位移時程數據具有隨機性，對定日鏡各工況選取對應測點的順風向位移響應譜進行分析，可以更好地揭示其結構動力反映機理。各工況對應結點歸一化位移響應譜如圖16，位移響應譜共振峰值能量如表3。可以看出，第一階振動模態對定日鏡的風致響應起主要作用。對于定日鏡上部RT測點對應一階模態的共振峰值較為明顯，隨著仰角增大，定日鏡上部RT測點的一階共振峰值表現出更高的能量。0°和30°時，定日鏡下部存在兩個振動峰值，第一個振動峰值對應一階模態，第二個振動峰值對應三階模態。隨著仰角增大，定日鏡底部第二個峰值逐漸減弱。從圖11可以看出，由于仰角較小時，氣流流向定日鏡表面后，較多的氣流從定日鏡底端繞流，在底部與地面之間形成豐富的旋渦，發生劇烈的旋渦脫落現象，引起多階共振響應，動力響應較為強烈。而隨著仰角增大，氣流在定日鏡底端分離，脫落漩在鏡面中上部再附著，上部的共振峰值表現出更高的能量。兩種模擬方法對功率譜能量峰值頻率的預測有較好的一致性，而LES比DES在定日鏡上部和下部的共振峰值預測出更高的能量。如圖17所示，這是由于LES比DES在定日鏡近壁面形成更多尺度細碎的小渦，并再附著在定日鏡背風面，使LES預測出更高的脈動能量造成的。

圖16 定日鏡各工況歸一化位移功率譜Fig.16 Normalized displacement power spectrum of heliostat with different elevation angle

(a) LES-00

(b) DES-00圖17 Y=3 m的X-Z剖面上LES和DES渦量分布圖Fig.17 Vorticity magnitude contours on X-Z plane at Y=3 m by LES and DES

表3 各工況歸一化位移功率譜共振峰值Tab.3 Peak value of normalized displacement power spectrum with different elevation angle

3.6 風振系數

進行結構風致響應分析可以得到設計人員關注的風振系數，其反映的是脈動風對結構的影響。各測點位移風振系數可由下式求出

(12)

定日鏡各部位風振系數隨鏡面仰角而變化，對定日鏡各工況選取對應測點的風振系數進行分析，得到定日鏡不同測點各工況風振系數如圖18。鏡面仰角為0°時，定日鏡底端有豐富的渦旋脫落現象，脈動響應影響較大，此時下部的風振系數最大，約為3.1。隨著鏡面仰角增大到60°，在鏡面底端分離的渦旋再附在定日鏡中上部，中部的風振系數增大到2.0，上部的風振系數增大到3.4。

圖18 定日鏡不同測點各工況風振系數Fig.18 Wind vibration coefficient of heliostat with different elevation angle

4 結 論

本文采用LES和DES，對0°風向角下0°、30°、60°鏡面仰角的多種定日鏡典型工況進行數值風洞模擬計算，利用ANSYS軟件建立有限模型，得到定日鏡結構的動力特性，將模擬結果與風洞結果進行對比分析，主要結論如下：

(1) 運用LES和DES可以較準確地預測定日鏡順風向等效風荷載。LES的結果更接近風洞試驗結果，但兩種模擬方法差距不大，其結果差距在2%以內。

(2) 兩種方法模擬的流場表現出相同的規律。隨著仰角的增大，定日鏡底部的渦旋出現脫離再附著現象，且尾跡區發展得更狹長且集中。這在定日鏡群中，可能會對后方定日鏡產生不利影響。

(3) 隨著仰角增大，流場對定日鏡影響發生變化，定日鏡下部的共振峰值能量逐漸減小，上部的共振峰值能量逐漸增大。兩種方法的位移響應譜都表現出相同的趨勢，但LES在結構近壁面的渦旋尺度更加細碎，對定日鏡的共振峰值預測出更高的能量。

(4) 定日鏡風振系數各部位最不利工況隨鏡面仰角變化。鏡面仰角為0°時，下部的風振系數最大為3.1。中上部的最不利工況發生在鏡面仰角為60°時，風振系數分別為2.0、3.4。

(5) 對于豎向相關系數，兩種方法的模擬結果都有相同的變化趨勢，DES結果的相關性均強于LES的結果。

(6) 綜合本文的模擬結果，結合DSRFG與LES、DES的模擬方法，可以有效獲得定日鏡各測點的風荷載時程。結合有限元分析，可為定日鏡的結構抗風設計提供初步的參考。