膜圓頂及塔架對南極望遠鏡風載影響的仿真分析*

任志敏 溫海焜 李國平 王 圣 趙秀勇

(1 中國科學院南京天文光學技術研究所 南京 210042)

(2 中國科學院天文光學技術重點實驗室(南京天文光學技術研究所) 南京 210042)

(3 中國科學院大學 北京 100049)

(4 國家能源集團科學技術研究院有限公司 南京 210023)

1 引言

南極大陸被證實是地球上最好的天文臺址,其中,位于南極大陸最高點的Dome A地區具有獨特的天文觀測條件: 超長的極夜時間(每年的冬季超過4個月)、低風速(通常低于10 m · s-1)、低溫干燥(冬季溫度在-70°C--60°C之間且夏季溫度低于-40°C;冬季的相對濕度在25%到45%之間)[1]、極好的視寧度條件(夜間視寧度中值為0.31′′,最低能夠達到0.13′′)、極薄的大氣邊界層(邊界層中值為13.9 m)[2-4].上述這些特點使得Dome A成為絕佳的天文觀測臺址.

為了充分利用Dome A的天文觀測優勢,可以利用塔架將望遠鏡安裝在大氣邊界層上方,采用對望遠鏡觀測影響極小的天文圓頂,隔絕熱源對觀測環境的影響并減小望遠鏡附近的風速以避免風載造成的望遠鏡振動與鏡面變形.綜上,我們可以考慮以下方案: (1)將望遠鏡安裝在高于Dome A大氣邊界層的塔架之上;(2)使用輕質膜結構圓頂.由于輕質膜圓頂為完全開放式的結構,而且鋼與膜的組合結構有輕便、吸熱量小的特點,所以圓頂視寧度極佳[5].此外,該類圓頂的整體質量和單個零件尺寸均很小,容易滿足Dome A運輸和安裝的要求;(3)用保溫材料將散熱結構(電機、塔架等)包裹進行隔熱處理[6];(4)采用風屏結構以減小望遠鏡周圍的風速和湍流[7].

計算流體力學(Computational Fluid Dynamics,CFD)在望遠鏡的風載分析中得到了廣泛的應用.De 乙oung等[8]對Gemini望遠鏡圓頂內的氣體流動進行了數值模擬,指出湍流隨風環境和隨動圓頂的開口方向變化而發生相應的改變;三十米望遠鏡(Thirty Meter Telescope,TMT)項目組[9]采用CFD的方法對TMT工作過程受到的外流場進行了模擬仿真,將模擬的風載數據導入望遠鏡集成模型中,從而對望遠鏡的結構進行了優化;中國科學院南京天文光學技術研究所的徐江海等[10]利用CFD分析了經典形式的球形、柱形、蚌殼折疊式圓頂在風載作用下對風速及湍動能的影響;中國科學院光電技術研究所的潘年[11]分析了圓柱形開放式圓頂在風載下的風速及湍動能的分布,同時將CFD分析結果與Gemini望遠鏡的實測結果進行了對比,證明了CFD分析的可靠性.

本文對望遠鏡、圓頂及塔架進行了建模,利用CFD進行數值模擬分析,當圓頂完全打開時,在給定風速和溫度等條件下,分析了不同風向角(α)、方位軸轉動角(β)及鏡筒轉動角(γ)時,望遠鏡周圍的風速、湍動能(Turbulence Kinetic Energy,TKE)、光程差(Pv)等信息,并且研究了風屏對望遠鏡周圍溫度及光程差的影響.旨在用數值方法評價風載對望遠鏡周圍光學傳輸環境的分布特性.

2 數值仿真原理

2.1 湍流模型的數學描述

本文采用標準k-ε湍流模型,這種模型有著計算穩定,好的經濟性和高計算精確度的優點[12].其輸運方程為:

式中k為湍動能、ε為湍流動能耗散率、ρ為流體密度、t為時間、μ為動力粘度、μt為湍流粘度、ui為流速、xi和xj為相互垂直的兩坐標軸、σk和σε為湍流動能k和耗散率ε對應的普朗特數(分別取為1和1.3)、Gk為由平均速度梯度引起的k的產生項、C1ε和為經驗常數(分別取為1.44和1.92).

2.2 湍流對光程差的影響

大氣溫度和速度場的不均勻形成了大氣湍流,湍流會造成流場的密度及折射系數隨時間的不均勻變化.從光波波前的角度分析大氣湍流對望遠鏡成像質量的影響可以分為以下兩種: 第1種是像運動與像模糊,主要表現為波前傾斜造成的像點圍繞著“平均位置”不規則的抖動,以及位相畸變帶來的不同像點位移的距離不同;第2種是像的閃爍,是對數振幅漲落造成的成像強度的變化[13].

湍動能是湍流速度脈動的度量,其大小直接關系到流體內部動量、熱量和水汽的輸運.湍動能定義為:

湍動能從大尺度渦旋向小尺度渦旋傳遞,部分湍動能在各級旋渦轉為流體熱能,進而導致大氣溫度起伏.折射率隨介質的密度和成分而變化,大氣折射率場的結構是由溫度場、濕度場和壓力場的結構決定的,在空氣中可以用柯西-勞倫斯公式表示為[14]:

式中N為空氣折射率;T為絕對溫度,單位為K;λ為光的波長,單位為μm;P為干燥的大氣壓強,e為水汽的壓強,單位均為hPa.

光波通過湍流介質的光程公式為:

式中L為光在介質中傳播的幾何路程,u、v、w為笛卡爾坐標系的3個方向.由(3)式可知,折射率的變化取決于環境溫度,大氣壓強和水汽含量.Dome A的空氣濕度非常低,同時壓力的波動對折射率的影響也可以忽略不計,折射率主要受湍流動能帶來的溫度波動的影響[14].由(4)式可知,折射率的變化會導致光束的光程差異變大,使來自天體的星光波前畸變,即光強和相位發生改變,導致成像圖像發生模糊、抖動、平移等現象[15].圖1表示了湍動能對視寧度的影響途徑.

圖1 湍動能對視寧度的影響(θdome為圓頂視寧度,ΔTc-a為圓頂內外的溫度差)Fig.1 Effect of TKE on seeing (θdome is the seeing degree of the dome,ΔTc-a is the temperature difference inside and outside the dome)

3 數值仿真

3.1 計算模型

望遠鏡的研制需要經歷概念設計、詳細設計、加工與檢測、外場建設、集成和調試等階段.對于方案制定的概念設計階段,望遠鏡的詳細結構以及臺址周邊環境的具體情況尚未明確.鑒于這些不確定性,在研究風載對望遠鏡的影響時,本文對望遠鏡、圓頂、塔架以及外部環境做出了適當的簡化.仿真望遠鏡為口徑2.5 m的地平式望遠鏡.輕質圓頂的直徑為12 m,兩側為驅動電機,圓頂保持完全打開的狀態時高度為0.42 m.塔架包含圓頂塔架和望遠鏡塔架兩部分,其中內側的望遠鏡塔架為傾角7°向上收縮的正六邊形結構,分為3層結構,每層的側面都由交錯的斜桿構成,塔架底邊邊長為4.3 m,高度為10.7 m;外側的圓頂塔架為傾角10°向上收縮的四邊形結構,分為兩層結構,每層的側面都由交錯的斜桿構成,塔架底邊邊長為17.2 m,高度為15 m;兩個塔架每個桿件的截面尺寸均為0.2 m×0.2 m.仿真模型的整體結構如圖2所示.通過改變風向角(α)-方位軸轉動角(β)-鏡筒轉動角(γ)的組合(為表述方便,文中形如0-0-120代表模型處于風向角0°、方位角0°、高度角120°的位置)實現模型在不同運行位置時受風載變化的仿真模擬.

圖2 仿真模型結構圖(圖中的風向角-方位軸轉動角-鏡筒轉動角為0-0-30,風向初始角度與塔架側面垂直,方位軸初始角度與風向初始角度同向,鏡筒轉動角為望遠鏡觀測方向與水平面的夾角,圖中的箭頭表示風向角、方位軸轉動方向和鏡筒轉動方向)Fig.2 Three-dimensional (3D) model of telescope and accessory equipment (In the figure,the wind direction angle-azimuth axis rotation angle-lens tube rotation angle is 0-0-30.The initial angle of the wind direction is perpendicular to the side of the tower,and the initial angle of the azimuth axis is in the same direction as the initial angle of the wind direction.The rotation angle of the lens tube is the angle between the observation direction of the telescope and the horizontal plane.The arrows in the figure indicate the wind direction angle,the rotation direction of the azimuth axis and the rotation direction of the lens tube)

3.2 邊界條件

通常,望遠鏡的工作上限風速為10 m·s-1,所以設置地面15 m高處為10 m·s-1的入口風速以得到極限工作情況下2.5 m望遠鏡周圍的流場分布情況,其他不同高度的入口風速值按照Dome A風速分布(5)式[16]設置;同時取Dome A的平均溫度213.15 K及大氣壓強0.5358 atm作為仿真流體域的環境溫度及壓強;設置空氣為不可壓縮理想氣體;出口為自由流出邊界;圓頂、望遠鏡和塔架設定為粗糙度0.5、無滑移的壁面;流場的其他壁面設置為可滑移的壁面.

3.3 網格劃分

為了提高計算效率、計算收斂性和計算速度,在模擬中采用了分區域劃分網格的方法,在遠離分析目標的流體域采用平均網格尺寸為100 mm的六面體網格,在分析目標周圍流場參數梯度變化較大的區域將網格加密為平均尺寸為50 mm的四面體網格,對分析目標表面進一步加密為平均尺寸25 mm的四面體網格,整個計算區域的網格數在各仿真條件下均不少于340萬,網格的劃分情況如圖3所示.

圖3 部分流體域的網格劃分Fig.3 Meshing of partial fluid domains

3.4 計算結果

3.4.1 風速分析

對風向角(α)分別為0°、30°、60°、90°(120°、150°、180°與60°、30°、0°風向角完全對稱)時,不同β和γ對應的2.5 m望遠鏡周圍穩態風速分布進行分析.由于望遠鏡鏡筒附近以及觀測方向的風速分布會影響到望遠鏡的觀測,所以本文在進行分析時主要考慮鏡筒觀測方向10 m內以及鏡筒周圍5 m內的風速分布.

首先分析α對穩態風速分布的影響,如圖4所示,當γ相同時,只要任意兩組的α-β值相同,則兩組間的速度云圖分布也相同,由此可認為塔架和膜圓頂周圍的風速分布對α的改變不敏感,α=30°、60°、90°時的鏡前風速分布規律均可以由α=0°時的分布結果來表示.為了直觀地表示鏡筒觀測方向的風速分布,圖中云圖平面均垂直于地面和望遠鏡鏡面且過鏡面的中心點.

圖4 α=0°時2.5 m望遠鏡周圍穩態風速分布云圖Fig.4 Cloud map of steady-state wind speed distribution around 2.5 m telescope at α=0°

對圖4做橫向對比分析.當β=0°、30°、60°時,隨著γ逐漸增大鏡筒觀測方向越來越背離來流方向,所以鏡面前方的風速有顯著的減小,鏡面所受風載也會減小[17];但γ的逐漸增大導致鏡面前方的風速分布出現明顯的分層結構,風速的分層會導致空氣結構的分層,空氣的折射率在風速的突變層間產生較大的變化,從而影響光線的穩定傳播.當β=90°時,由于鏡面與風向的角度始終為90°,所以鏡筒觀測方向的風速分布在各γ均一致,均是在鏡面處有明顯的風速分層結構.

做縱向對比分析.在γ分別為30°、60°、90°、120°、150°時,2.5 m望遠鏡周圍的風速分布均隨著β的增大而增大,鏡面前方的風速分層越來越緊密且逐漸貼近鏡面.風向來流與鏡面的角度逐漸增大,鏡面所受的風作用力變小,鏡面的變形會減小.

為了定量分析光路上的風速變化情況,本文按照圖5的排布方式取垂直于鏡面的25條長度為10 m的光路,每條光路上均勻取1000個采樣點進行采樣,圖6為風速分布采樣圖.

圖5 垂直于鏡面采樣線的位置及編號Fig.5 Position and number of the sample line perpendicular to the mirror

圖6 α=0°時垂直于鏡面方向的風速分布Fig.6 The wind speed distribution perpendicular to the mirror in the case of α=0°

對圖6分析得到,當鏡面處于迎風(γ＜90°)狀態時,同一組內25條取樣線上發生風速波動的空間范圍一致,但鏡面邊緣取樣線(6、8、10、12、14-25)和鏡面內側取樣線(1、2、3、4、5、7、9、11、13)分別呈現兩種不同的變化趨勢.各取樣線上的風速值均在離鏡面約2 m的范圍內發生劇烈的變化,并且在2 m后逐漸趨于穩定風速值10 m·s-1.鏡面邊緣取樣線(6、8、10、12、14-25)的波動變化趨勢為由鏡面的0 m · s-1迅速增大,然后在離鏡面0.3 m遠處出現一個明顯的風速拐點,此后風速的變化變得緩和;鏡面內側取樣線(1、2、3、4、5、7、9、11、13)的波動變化趨勢為由鏡面的0 m·s-1較平穩地增加到10 m·s-1.此外,γ越大風速值的變化范圍越大,在0-0-30時,風速值的變化范圍僅為0-10.5 m·s-1,在0-90-90時,風速的變化范圍可達到0-13.9 m·s-1.

當鏡面處于背風(γ＞90°)狀態時,同一組內25條取樣線上的風速波動范圍和變化趨勢一致,每條取樣線上的風速值均是從鏡面處的0 m·s-1迅速增大到2-4 m·s-1,隨后下降到0.5-1.5 m·s-1,之后迅速增大到最大風速值后緩慢歸于一個穩定的風速值.但各取樣線上發生風速波動的空間范圍則隨著鏡面的背風程度(γ)的增大而增大,在0-60-120時,25條取樣線上的速度波動發生在離鏡面3 m的范圍內,在0-30-150時,各取樣線上的速度波動則直到離鏡面6 m處才逐漸平緩.

此外,在所有來流風與鏡面平行的情況下,鏡面前方的氣流變化規律均相似,由鏡面處的0 m·s-1增加到離鏡面0.3 m處的14 m · s-1,然后緩慢變化到穩定風速10 m · s-1.總體而言,在迎風狀態下,鏡面前方的空氣流速變化更加穩定,發生變化的空間范圍也更小,所以鏡面前方空氣的狀態更為均勻、光線的傳輸更加穩定.背風狀態下,鏡面前方的空氣流速變化較不穩定,空氣密度的均勻性也必然較差,從而影響到光程差的大小.本文嘗試使用風屏結構緩和這種劇烈的風速變化,改善望遠鏡周圍的風速分布環境.

3.4.2 湍動能分析

對α分別為0°、30°、60°、90°時,不同β和γ對應的2.5 m望遠鏡周圍的TKE分布進行分析.由于望遠鏡鏡筒附近以及觀測方向的湍動能分布會影響到望遠鏡的觀測,所以本文在進行分析時主要考慮鏡筒觀測方向10 m內以及鏡筒周圍5 m內的湍動能分布.

首先分析α對湍動能分布的影響,如圖7所示,與速度分布規律類似,當γ相同時,只要任意兩組間的α-β值相同,則兩組的湍動能云圖分布相同,由此可以認為塔架和圓頂周圍的湍動能分布對來流風向(α)的改變不敏感,α=0°、60°、90°時鏡前湍動能分布規律均可以由α=0°時的分布結果來表示.為了直觀地表示鏡筒觀測方向的湍動能分布,圖中云圖平面均垂直于地面和望遠鏡鏡面且過鏡面的中心點.

然后對圖7做橫向對比分析.當β=0°、30°、60°時,隨著γ逐漸增大鏡筒觀測方向的湍動能數值和分布范圍逐漸增大,導致環境溫度變得更不均勻、空氣折射率的不均勻使光程差變大.當β=90°時,由于鏡筒與風向的角度始終為90°,所以鏡筒觀測方向的湍動能分布在各γ均一致,均是在鏡面表面有薄薄一層湍動能稍大的區域.

最后對圖7做縱向對比分析.當γ=30°、60°時,不同β下鏡面前方的湍動能分布差異不大,均是鏡面表面有薄薄一層湍動能值較小的區域;當γ=90°時,隨著β的增大,湍動能的數值有所減小,但是湍動能均集中分布在鏡面附近沒有變化;當γ=120°、150°時,隨著β的增大,望遠鏡觀測方向的湍動能分布越來越緊密且逐漸貼近鏡面.

同樣,對光路方向的湍動能做定量分析,每條取樣線上的湍動能分布如圖8所示.

圖8 α=0°時垂直于鏡面方向的湍動能分布Fig.8 TKE distribution perpendicular to the mirror direction at α=0°

對圖8分析得到,當鏡面處于迎風(γ＜90°)的狀態時,不同γ取值下各組25條取樣線上發生TKE波動的范圍與趨勢均一致,每條取樣線上的湍動能值均在離鏡面約2 m的范圍內發生劇烈的變化,并且在2 m后逐漸趨于0;波動變化趨勢均為從鏡面表面到離鏡面0.2 m的范圍內由接近0值迅速增大到最大值,然后在離鏡面2 m的范圍內迅速降到很小的值,之后則會緩慢趨于穩定.迎風狀態下各組的湍動能值均比較小,在0-0-60湍動能的最大值僅為約0.2 J·kg-1.

當鏡面處于背風(γ＞90°)狀態時,每組25條取樣線上的湍動能波動趨勢在不同γ取值下均一致,由鏡面處大于0的值逐漸增大到最大值,之后隨著遠離鏡面逐漸歸于0值;但各取樣線發生湍動能變化的空間范圍則隨著鏡面的背風程度(γ)的增大而增大.在0-0-120時,各取樣線上的湍動能波動發生在離鏡面4 m的范圍內,在0-0-150時,各取樣線上的湍動能波動則直到離鏡面8.5 m處才逐漸平緩.此外,湍動能值波動范圍隨著背風程度的增大而略有減小.背風狀態下各組的湍動能值均比較大,在0-60-120時,湍動能的最大值達到13.5 J·kg-1.

3.4.3 光程差分析

由于光程的大小受空氣折射率的影響,空氣折射率又由湍動能造成的溫度變化決定,不同的湍動能分布就代表著不同的光程差,由前文分析得到α=30°、60°、90°的湍動能分布情況可由α=0°來表示,所以在此只分析α=0°時,最大光程差Pv隨β,γ值的變化關系.

由(3)-(4)式可得到不同β、γ時各組取樣光線所經過的實際光程S1、S2、...、S25,然后可以求出最大光程差Pv=Smax-Smin.圖9為β=0°、30°、60°、90°時Pv隨γ值變化的情況.

圖9 β=0°、30°、60°、90°時Pv隨γ的變化Fig.9 Change of Pv with γ at β=0°,30°,60°,90°

β取不同值時曲線的變化趨勢相同,均是隨著γ的增大Pv先緩慢減小,在γ=90°時達到最小值,然后再逐漸增大,并在γ=150°時達到最大值.整體來看,β值越小光程差越大,最大的光程差為0-0-150時的229.4 nm;β=90°時光程差的變化幅度不大且數值較小,最小為0-90-60時的24 nm.

3.4.4 風屏對湍流的改善作用研究

在望遠鏡周圍安裝風屏可以有效減少望遠鏡附近風速突變及湍動能,但同時也會造成一定程度的溫升[7].本文研究了風屏對2.5 m望遠鏡周圍風速、湍動能及溫度的影響.由于Dome A地區沒有固定的主風向,所以將風屏模型簡化成360°圍繞圓頂塔架頂端外側均勻排布的矩形長桿,每根桿的尺寸均為0.1 m×0.1 m×4.2 m,桿與桿的間距為0.05 m,桿的高度略低于望遠鏡觀測時鏡筒的最大高度,如圖10所示.

圖10 望遠鏡周圍風屏的結構及溫度采樣點Fig.10 Structure of the wind screen around the telescope and temperature sampling points

圖11和圖12分別為α=0°時,不同β、γ對應的湍動能和風速分布情況,為了直觀地表示鏡筒觀測方向的湍動能和風速分布,圖中云圖平面均垂直于地面且平行于望遠鏡觀測方向.由圖11分析得到在不同β、γ時,2.5 m望遠鏡觀測方向與遠處自由穩定大氣的湍動能值差別不大,均為略大于0的數值.從圖12可以得到在不同β、γ時,高速風在通過風屏后,風速有大幅度減小,同時由于膜圓頂對風的進一步阻擋作用,望遠鏡周圍的風速值＜1 m·s-1,望遠鏡所受風載顯著降低.

圖11 加裝風屏后α=0°時望遠鏡和圓頂及塔架周圍的湍動能云圖Fig.11 TKE cloud map around the telescope,dome and tower at α=0° after installing the wind screen

圖12 加裝風屏后α=0°時望遠鏡和圓頂及塔架周圍的風速云圖Fig.12 Cloud map of wind speed around the telescope,dome and tower at α=0° after installing the wind screen

為了定量研究風屏帶來的溫升效應,對望遠鏡周圍16個點的溫度值進行采樣分析,圖10為采樣點的分布情況,表1為采樣點的溫度.

表1中的數據并沒有顯示出望遠鏡周圍的溫度有特別的分布趨勢與規律性.但是,16個采樣點的最高溫度為213.196 K,比環境溫度213.150 K高0.046 K;平均溫度TA為213.194 K,較環境溫度高0.044 K.仿真關于風速、湍動能、溫度的結果與文獻[7]中得到的10 m·s-1的風速在通過風屏后最大溫升不大于0.050 K,風速在離風屏1 m遠、3 m高的位置處降為來流風速的1/3至1/4的結論一致.

同時對受風屏影響的最大光程差進行了分析,雖然湍動能的數值較無風屏的情況有很大的減小,但風屏帶來的溫升效應會直接影響到空氣折射率的變化.圖13為風屏作用下不同β時光程差隨γ的變化,與圖9對比可得到,受到風屏溫升效應的影響雖然Pv的變化趨勢一致,但不同β時的Pv值均有顯著的增大.可同時考慮到鏡面受到的風載壓力也會變小,從而由鏡面的面形變化導致的成像變化也會減小[17],這需要今后進一步的計算和實驗.

圖13 加裝風屏后,β=0°,30°,60°,90°時Pv隨γ的變化Fig.13 After installing the wind screen,change of Pv with γ at β=0°,30°,60°,90°

4 結論

本文借助CFD和有限元方法,分析了2.5 m望遠鏡及膜圓頂和塔架結構,在給定風速和環境溫度的條件下,得到了在不同的風向角-方位軸轉動角-鏡筒轉動角度時望遠鏡周圍空氣流速、湍流強度等的變化情況以及對光程差的影響;同時還分析了風屏對望遠鏡周圍的空氣流速、湍流強度及溫度的影響.對后續天文望遠鏡塔架結構設計,膜圓頂的結構優化以及風屏結構的應用提供了參考.