大射電望遠鏡FAST風環境數值模擬研究*

吳明長，王啟明，郭永衛，趙保慶

(中國科學院國家天文臺，北京 100012)

計劃在中國貴州某喀斯特洼地內建設的500 m口徑球面射電望遠鏡FAST，通過主動控制在觀測方向形成300 m口徑瞬時拋物面。采用光機電一體化的索支撐輕型饋源平臺及平臺內二次調整裝置，在饋源與反射面之間無剛性連接的情況下，實現高精度的指向跟蹤［1－3］。其主動反射基準面為半徑300 m、口徑500 m的球冠。在主動反射面初步設計過程中，風荷載是FAST反射面主要外部載荷［4］，需要進行抗風設計，而現行建筑設計規范GB 50009-2001［5］中并無合適的抗風設計體型系數，因此需要通過數值模擬獲得抗風設計的依據。針對FAST進行的風環境模擬前期研究［6］，給出了FAST反射面風壓系數隨風向的變化，同時給出了FAST索網結構的風振響應初步分析。但隨著項目的進展，該研究結果已不能滿足工程需要，需要進行更深入的研究。

目前，現場測試、風洞試驗、數值模擬是獲得建筑物風壓系數的3種途徑。現場測試最為直觀和真實，但耗費人力物力時間較多。存在如何布置測點，某些關鍵位置難以布置測點及獲得結果種類有限的問題。風洞試驗成本高，而且由于邊界層風洞模擬的局限，結果是否可靠還存在爭議［7］。基于計算流體動力學(Computational Fluid Dynamics，CFD)的數值模擬存在邊界條件難以符合實際情況，理論模型欠完備等缺點。盡管如此，CFD數值模擬成本相對較低，可以給出很多無法測量的結果，對分析流場規律和特點有著難以替代的作用。本文給出了射電望遠鏡FAST風環境數值模擬的主要結果，探討了有關FAST風環境的影響因素和規律。

1 數值模擬方法與模型

1.1 理論模型

大氣邊界層內繞建筑物流動的空氣可以假設為低速、不可壓的粘性牛頓流體，基本控制方程是質量守恒方程和動量守恒方程［8］:式中，Sm是質量源項;ρ是靜壓;τij為應力張量，可以用速度和坐標分量表達;ρgi是重力;Fi是其它體積力或用戶自定義源項。

計算流體動力學軟件FLUENT［9］基于有限體積法，將計算域離散化為一系列控制體積，在控制體上求解質量、動量、能量、組分等通用守恒方程，將偏微分方程組離散化為代數方程組，用數值方法求解偏微分方程組以獲取流場參數的離散數值解。然而，由于實際湍流流場的復雜性，在求解過程中，直接求解這些方程組，尤其是動量方程組，對計算機計算能力要求很高，因此在工程應用中，一般將湍流各變量在時域看作平均和脈動兩部分，將動量方程逐項平均得到不封閉的雷諾方程，并且提出多種湍流封閉模式求解這些雷諾方程［9］。不同的湍流封閉模式有著不同的特點，常用的幾種模型在預測鈍體表面風壓分布時結果差別不大。因此，盡管標準k-ε模型對于分離的預測不夠準確，但由于工程應用精度要求不高，并且由于標準k-ε模型工程上應用廣泛，所以本文有關計算中一般采用標準k-ε模型。

1.2 幾何模型

經過分析試算，確定了將口徑500 m的FAST反射面和水平尺度約3 km×3 km的山體模型放置于XYZ方向尺度分別為40 km×10 km×5.5 km的計算域中。其中，來流方向的計算域尺度為40 km，反射面中心位置距離來流入口20 km，垂直于來流方向的計算域尺度為10 km，計算域底部海拔為500 m，頂部海拔為6000 m。為減小網格劃分的難度并減小計算量，將計算域分為外部和內部兩個。

保持坐標系及內部計算域和計算對象不變，旋轉外部計算域至不同角度，可方便地獲得用于不同風向工況計算的幾何模型，在此基礎上形成網格模型和CFD計算模型。風向約定沿用地學中常用的表達習慣，以正北為0點鐘或12點鐘方向，正東為3點鐘方向。X軸正向指向3點鐘方向，Y軸正向指向12點鐘方向，Z軸正向豎直向上。以3點鐘方向為例，外部計算域和內部計算域的相對位置關系如圖1。其中，中心區域包含了FAST反射面結構、山體表面及可能的擋風墻結構，如圖2。

圖1 3點鐘風向下外部計算域和內部計算域Fig.1 External and internal calculation domains with the wind direction of 3 o’clock

圖2 山體表面、FAST反射面及擋風墻結構的幾何模型(a)3 km×3 km山體表面;(b)中心350 m范圍開挖后地形;(c)FAST反射面在洼地內的相對位置;(d)不同高度、不同徑向位置的擋風墻Fig.2 Geometric model of the hill surface，the reflector of FAST，and the windbreaks(a)Hill surface of 3km×3km;(b)The excavated terrain within 350m from the center;(c)The relative position of the reflector in the depression;(d)Windbreaks of different heights and radial positons

1.3 網格模型

山體表面形狀本身是不規則的，反射面臨近地域在開挖后帶來更多片狀不規則區域，這種不規則性，給幾何建模、網格劃分及生成有效的CFD計算模型帶來困難。采用分區劃分及非結構化網格劃分的方法，獲得了多種計算工況下的網格模型，典型網格劃分結果的部分表面網格見圖3。相應的網格設置如表1，不同范圍的邊界之間，通過網格劃分軟件的邊界網格一致功能實現。

圖3 地形表面的網格(a)山體表面網格，由外到內逐漸加密;(b)距離反射面中心350 m范圍內的地表網格Fig.3 The CFD mesh grid at the surface of the terrain(a)Mesh grid at the hill surface(appearing denser in the inner area);(b)The mesh grid at the terrain surface within the radial distance of 350m from the center

表1 不同地形范圍內的網格大小設置Table 1 The mesh sizes for the different regions of the terrain surface

1.4 邊界條件

計算域入口邊界條件采用剪切流風剖面，根據建筑物結構荷載規范GB 50009-2001，近地面風速沿高度方向近似按指數分布，如下式:

對于山地和叢林，指數α取0.16，速度基準高度Vb一般以地面以上Zb=10 m計，Vb取規范中要求的極限風速14 m/s。

來流湍流強度采用澳大利亞規范中第2類地貌的取值，

顯然，來流速度和湍流特性與坐標相關，這種變化的邊界條件在FLUENT軟件中可以通過用戶自定義函數(User Defined Function，UDF)實現。

為保證植被光照需求，所設計的FAST反射面其透孔率約為35.4%。反射面各單元是多孔薄面板形式，相對于總體尺度，1.2 mm板厚和5 mm孔直徑非常小，而且孔的數量又很多，試圖直接將孔的模型反映在計算模型中是不現實的，直接忽略也缺乏依據。為了考察反射面板透孔率對流場的影響，在后文中嘗試使用FLUENT軟件中多孔介質跳躍的邊界條件模擬透孔反射面的空氣動力學特性。具體做法是首先通過數值模擬計算空氣流過透孔面板后壓降隨速度的變化關系，然后根據這種速度壓降關系，計算得到多孔介質跳躍的邊界條件所要求的面滲透性α和慣性阻力因子C2。出口采用完全發展邊界條件，即流場任意物理量沿出口法向梯度為0。計算域頂部和底部取自由滑移的運動壁面條件，地面和物體表面取靜止壁面條件。

考慮到山體表面巖石、樹木等不平整因素的影響，將山地表面分為兩個區域，距離反射面中心350 m之外和之內兩部分。因為350 m之內大部分需要開挖和施工，基本沒有高大樹木等，取較小的粗糙高度，而對350 m之外的部分，則取較大的粗糙高度。

2 計算方法和結果

計算中采用基于壓力的求解器，穩態計算，湍流模型采用標準k-ε模型，采用標準壁面函數近似處理近壁網格劃分不夠細致的問題。計算過程中根據實際迭代情況，調整松弛因子，一般收斂準則為殘差小于1×10－3，同時以關鍵物理量是否穩定作為收斂與否的判據。

首先通過分析和試算確定了入口速度計算的基準高度，然后利用數值模擬獲得了FAST反射面透孔面板的速度壓降特性，并以此為基礎計算了是否考慮面板透孔率兩種情形的反射面風壓系數分布特點及其與風向的關系，最后給出了擋風墻效果數值模擬研究的部分結果。

2.1 基準高度對計算結果的影響

確定計算域入口速度剖面時，由于FAST臺址洼地地表起伏很大，不能作為CFD分析模型中Z向坐標的起點。需要引入基準高度Hb，才能恰當描述計算域入口速度剖面，如下式:

結合山地表面海拔高度，以及計算域頂部和底部的海拔高度，對若干不同Hb值進行了試算，得到反射面下表面和上表面風壓系數之差與Hb的關系，如圖4。

首先，基準高度應小于所涉及的地形邊界的最低海拔高度，對本文所涉及的地形，這個高度約為810 m。另外，由圖4可以看出，在810 m上下40 m的變化范圍內，反射面風壓系數的極值和均值對參考高度是不敏感的，因此后續計算中，一系列分析計算中基準高度Hb的取值均為810 m。

2.2 透孔面板的速度壓降特性

圖4 反射面下表面、上表面風壓系數差值與基準高度的關系(其中，Diff表示下表面減去上表面的差值，Diff Min表示在整個反射面范圍內各計算節點上這個差值的最小值，Diff Max表示相應的最大值，Diff Aver表示這個差值的平均值，Diff RMS表示這個差值的均方根值，下同)Fig.4 The relation between the wind pressure coefficient difference from the upper to lower surfaces of the reflector and the reference altitude(The“Diff”denotes the wind pressure coefficient of the lower surface subtracted by the one of the upper surface.The“Diff Min”，“Diff Aver”，and “Diff RMS”denote the minimum，the average，and the rms values of the coefficient differences between the reflector surfaces)

FLUENT軟件提供了多孔介質跳躍邊界條件，用于模擬速度或壓降特性均為已知的薄膜，是多孔介質模型的一維簡化，可用于流體通過篩子、過濾器、過濾紙和多孔板的壓降模擬。這種邊界類型需要設置面滲透性數值α、介質厚度t、壓力變化系數C23個參數。

對于簡單各項同性多孔介質，流體流過多孔介質的壓降Si由粘性損失和慣性損失兩部分組成［9］(實際應用中應考慮介質厚度t)，除物性參數μ和ρ外，面滲透性α和慣性阻力因子C2是需要求解的。

另外，物理實驗表明，流體穿過多孔板或絲網的壓降與速度成二次函數關系:

通過物理實驗獲得慣性阻力項系數a1和粘性阻力項系數a2后，就可以通過對比獲得面滲透性α和慣性阻力因子C2。

因此，為了以多孔介質跳躍邊界條件表示FAST反射面，首先需要獲取描述速度和壓降特性的系數a1和a2。另外，既然物理實驗可以獲得這樣的壓降關系，在條件有限的情況下，也可以通過數值模擬近似獲得a1和a2的值。對于FAST反射面板的某局部孔來說，周圍的孔允許空氣的流出，所以在建立空氣流過孔板的CFD模型時，應允許空氣自由地從孔板周邊流過。圖5給出了若干不同流速下，空氣流過FAST孔板(透孔率35.4%)的壓降。

根據速度壓降特性，就可以通過二次曲線擬合得到慣性阻力項系數a1和粘性阻力項系數a2，并通過對比獲得面滲透性α和慣性阻力因子C2。

氣動系數K值可以通過速度壓降特性擬合得出，所得K值可以與Dryden H L和Schubauen G B于1947年的實驗結果作對比［10］(圖6)。由圖6可以看出，35.4%透孔率的FAST反射面面板的數值模擬結果與實驗結果誤差小于3%。

圖5 FAST反射面板的速度壓降特性曲線Fig.5 The simulated relation between the wind-pressure drop and the wind speed for the FAST reflector panel under consideration

圖6 FAST反射面板氣動系數值與實驗結果對比Fig.6 The comparison of the simulated aerodynamic coefficients of the FAST reflector panel to experimental results from references

由于FAST反射面板的厚度和孔徑相比較小，并且不同來流方向等效透孔率近似相等，所以平面孔板的數值模擬和實驗結果用于模擬空氣流過FAST球冠形透孔反射面的流動特性時具有一定的合理性，空氣是否垂直流向面板可以作為次要因素忽略。

2.3 面板透孔率對反射面風壓分布的影響

進行FAST反射面附近空氣流場特性數值模擬時，究竟是否需要考慮反射面面板的透孔率，需要通過比較才能給出結論。圖7給出了5點半風向下(該風向為通過不同風向比較計算確定的不利風向)，不考慮反射面透孔率和35.4%透孔率兩種情形的反射面風壓系數的分布特點，可以看出孔的存在總體效果是均化了反射面風壓系數分布。

FAST所處地域山地地表海拔高度差異很大，不同風向下反射面的風壓系數分布應當有所差異，圖8給出了是否考慮反射面透孔率兩種情形下，反射面風壓系數隨風向的變化。

從圖8可以看出，其一，考慮反射面透孔率時，反射面風壓系數的極值明顯減小;其二，不考慮反射面透孔率時，不利風向約為5點半方向，考慮反射面透孔率時，不利風向約為5點半方向和9點半方向。

圖9給出了是否考慮反射面透孔率兩種情況下，各風向下反射面風壓系數沿半徑方向的統計。

圖7 是否考慮面板透孔率兩種情形下反射面風壓系數分布(a)不考慮反射面透孔率;(b)反射面透孔率35.4%Fig.7 The wind pressure coefficients of the reflector panel for the two cases of perforation rate.(a)null perforation，and(b)35.4%perforation rate

圖8 是否考慮面板透孔率兩種情形下反射面風壓系數隨風向的變化(a)不考慮反射面透孔率;(b)反射面透孔率35.4%Fig.8 The wind pressure coefficients of the reflector panel for different cases of perforation rate and different wind directions.(a)null perforation，and(b)35.4%perforation rate

圖9 是否考慮面板透孔率兩種情形下各風向反射面壓力系數沿半徑方向的統計(a)不考慮反射面透孔率;(b)考慮反射面透孔率35.4%Fig.9 Statistics of the wind pressure coefficients at various radii of the reflector panel for different cases of perforation rate under and different wind directions.(a)null perforation，and(b)35.4%perforation rate

可以看出，與不考慮面板透孔率的情況相比，35.4%透孔率的反射面面板的風壓系數從遠離球心的壓力為主(負值)轉為指向球心的壓力(正值)為主。

2.4 擋風墻問題

一般認為在反射面周邊適當位置設置擋風墻能對某些山口的風速起到阻擋作用，從而改善反射面風壓環境。實際擋風墻效果如何，需要通過計算或試驗才能確定。

以距離中心275 m位置為例，在高度方向，根據擋風墻上沿與反射面邊沿的相對位置分為:高出邊沿10 m和5 m，與邊沿平齊(0 m)，低于反射面邊沿5 m、10 m、15 m、20 m，以及低于反射面邊沿40 m(無擋風墻)共8種情形，如圖10。另外，還可以考察同樣的擋風墻高度，擋風墻距反射面中心不同位置處的效果，如圖11。

圖10 不同高度擋風墻對反射面下表面、上表面風壓系數差值的影響(距離中心275 m)Fig.10 Results of the wind pressure coefficients for different heights of the windbreaks.(All for the position of275m from the center of the reflector)

圖11 不同位置擋風墻對反射面下表面、上表面風壓系數差值的影響(擋風墻高出反射面邊沿10 m)Fig.11 Results of the wind pressure coefficients for different positions of the windbreaks.(All for a height of windbreak of 10m above the edge of the reflector)

從圖10和圖11可以看出，擋風墻高度變化對反射面風壓系數極值改變的效果大于擋風墻徑向位置改變的影響。這說明建設擋風墻時，擋風墻距中心位置的距離選擇要求稍低。

此外，擋風墻效果如何，還要看其對反射面風壓系數分布的影響。圖12給出了反射面透孔率為35.4%時無擋風墻和有擋風墻的反射面風壓系數分布特點對比。除了少數局部位置有所變化外，總體范圍差異并不明顯。

圖12 考慮反射面透孔率35.4%時，無擋風墻和有擋風墻的效果(a)無擋風墻;(b)距離275 m處高出邊沿10 m的擋風墻Fig.12 Contours of wind pressure coefficients for the panel of 35.4%perforation rate with and without windbreaks(a)without windbreaks，and(b)with a windbreak at a distance of275m from the center and with a height of 10m above the edge of the reflector

3 結果分析和討論

對于考慮反射面透孔率的情況，采用多孔介質跳躍邊界條件模擬反射面，FLUENT計算結果給出的是上下表面之差，風壓系數規定指向上部(球心)為正值。而當不考慮反射面透孔率時，就需要計算反射面兩個側面風壓系數差值，同樣以指向上部(球心)為風壓系數正值。這樣，反射面某部分承受正風壓時，說明該部分承受的是指向球心的風壓，反之，反射面某部分承受負風壓時，說明該部分承受的是遠離球心的風壓。

對于CFD計算模型，邊界條件和物理屬性是建立在簡化和假設之上的，這決定了它獲得的總是依賴于模型的離散近似解。首先，計算依據的地形數據經過空間域采樣和插值等必要處理后存在誤差，計算中采用的幾何模型和實際對象的差異會帶來幾何模型的誤差。其次，對幾何模型進行網格劃分，有限的網格是對連續的幾何模型的一種近似，這是離散誤差。再次，計算中涉及的邊界條件，如來流特性和壁面屬性等，均與實際情形有差異。最后，數值計算過程中，數據截斷和舍入誤差，控制終止條件而導致的迭代誤差等均是CFD數值模擬中的誤差來源。

如果CFD數值模擬對象較為簡單，則有可能結合試驗或理論結果進行一定的誤差估計，但是對于本報告涉及的FAST風環境模擬，這兩種方法均不可行，試驗也僅僅能給出一部分的規律性結果，難以給出精確的結果。

因此，本文給出的模擬計算的結果其意義在于分析FAST臨近地域風環境的規律和特點，為后續抗風設計做準備，體現在以下幾個方面:

(1)通過分析和計算確定了810 m為計算基準高度。

(2)35.4%透孔率反射面板的速度壓降特性數值模擬結果與前人實驗結果一致性很好。是否考慮反射面透孔率，結果有較大差異。

(3)不同風向下反射面風壓系數差異較大，將不利風向5點半方向作為各種分析的基本風向。

(4)建設擋風墻時，距中心位置距離的選擇要求較低。

這些分析和計算結果為后續反射面受風荷載的變形計算以及反射面高精度測控打下了基礎。

4 結論

本文針對FAST反射面及其周邊山地地形，建立了幾何模型、網格模型和CFD分析計算模型。通過分析和試算確定了入口風速剖面的基準高度為810 m。并通過數值模擬獲得了35.4%透孔率的FAST反射面板的速度壓降關系，氣動系數結果與前人實驗結果誤差小于3%。模擬計算了不同風向下反射面風壓系數分布特征，得出的不利風向為5點半方向。研究還表明，擋風墻高度選擇的效果優于位置選擇的效果。

致謝:感謝Richard David先生提供文獻數據來源。感謝計算機網絡信息中心超級計算中心提供軟硬件計算資源。

［1］南仁東.500 m球反射面射電望遠鏡FAST［J］.中國科學G輯，2005，35(5):449－466.

［2］FAST項目組.FAST初設總冊 ［R］.北京:中國科學院國家天文臺，2008.

［3］Nan Rendong，Li Di，Jin Chengjin，et al.The Five-Hundred-Meter Aperture Spherical Radio Telescope(FAST)Project［J］.International Journal of Modern Physics D，2011，20(6):989.

［4］FAST項目組.500 m口徑球面射電望遠鏡(FAST)項目初步設計(主動反射面分冊)［R］.北京:中國科學院國家天文臺，2008.

［5］中國建筑科學研究院.建筑結構荷載規范(GB50009-2001)［S］.北京:中國建筑工業出版社，2006.

［6］林斌，范峰，錢宏亮.大射電望遠鏡FAST主動反射面在喀斯特地貌下的繞流數值模擬及其風振響應分析［C］//第四屆全國現代結構工程學術研討會論文集.寧波，2004.

［7］魏慶鼎，陳凱，盧占斌，等.風洞模擬大氣邊界層的若干問題 ［C］//第11屆全國結構風工程學術會議論文集.三亞，2004.

［8］賀德馨.風工程與工業空氣動力學 ［M］.北京:國防工業出版社，2006.

［9］Fluent6.3 User's Guide ［EB/OL］.http://www.ansys.com/products/fluid%2Ddynamics/fluent/.

［10］Dryden H L，Schubauen G B.The Use of Damping Screens for the Reduction of Wind Tunnel Turbulence ［J］.Journal of the Aeronautical Sciences，1947，14(2):221－226.