龍卷風維持特性的探索

邸浩宇， 徐晶磊， 高歌

(北京航空航天大學 能源與動力工程學院， 北京 100083)

在自然界中，龍卷風可以產生摧毀性的能量，但龍卷風是如何生成并獲得如此巨大的動能，特別是如何在較長時間維持旋轉一直是個迷。對于龍卷風的研究，國內外學者主要通過現場實測、實驗研究和數值模擬這3種手段。

在早期的現場實測中，由于龍卷風發生的地點和時間的不可預測性，以及測量技術的限制，導致實測研究發展很緩慢。Hoecker[1]針對1957年發生的Dallas龍卷風，通過走訪目擊者以及察看龍卷風的圖片和視頻資料，指出渦核區域存在下沉氣流并給出了龍卷風內部分位置的速度結構。21世紀初，多普勒雷達等的使用大大豐富了數據上的空缺，比較典型的是Alexander和Wurman[2-3]利用移動式多普勒雷達實測了1998年5月30日發生在南達科他州Spencer的龍卷風，并提供了該龍卷風不同高度和時間的速度分布。

相比現場實測，實驗研究和數值模擬更加方便易于收集數據。在龍卷風的實驗研究中，Chang[4]于1971年率先建立了龍卷風發生裝置并測定了龍卷風切向和徑向的速度分布特征。基于Chang[4]的研究，Ward[5]于1972年利用蜂窩板消除龍卷風豎向的渦量，建立了著名的Ward模型。此后，許多科研人員又在Ward模型的基礎上進行實驗并改進建立了諸多模型，其中普渡大學的Church等[6]引入旋轉網格柵改進了Ward模型并進行了實驗，發現龍卷風的渦可分裂呈現出多渦狀態。Jischke和Parang[7]利用Ward模擬器發現了核心半徑這一參數的重要性，并提出將渦流比、高寬比和渦流雷諾數作為龍卷風模擬設計的主要參數。上述實驗均未考慮龍卷風的橫向移動，直到愛荷華州立大學開創了此先河，于2003年建立了可以使龍卷風進行平移運動的模擬發生器。之后Sarkar等[8]利用其對龍卷風進行了風場特性的相關實驗研究，并取得了和實測基本一致的結果。

近年來，隨著計算機技術和計算流體力學(CFD)相關模擬軟件的迅猛發展，對龍卷風的研究更多地轉入了數值模擬方面。Hassenzahl[9]采用限制渦流的方式對龍卷風的渦流特性進行了研究。Maruyama[10-11]于2009年應用大渦模擬研究了龍卷風的渦流特性；并于2011年在龍卷風模型中加入碎片，分析了碎片的運動特性，且利用其驗證了龍卷風的渦核特性。Natarajan等[12-13]則研究了建筑物表面的粗糙度對龍卷風風場的影響，并于2011年通過數值模擬了Ward-type TVC型、WinDEEE Dome型和AVE型這3種龍卷風發生器形成的龍卷風風場的渦流特性，與實驗對比發現Ward-type TVC型與之更相符。徐楓等[14]針對單渦龍卷風的風場特性，分析了切向風速在徑向和高度的速度分布規律，并與Rankin渦模型和參數化氣旋模型相對比，驗證了模擬結果的合理性。

在以往對龍卷風的數值模擬研究中大多都以不可壓縮假設處理，忽略了這種極端大氣現象中氣流可壓縮性的影響。與此同時，忽略了溫度場內氣流之間的熱交換。但龍卷風多發于氣溫很高的夏季時節，而且其內外的溫度差很大，氣體溫度越高其包含的熱能越多，在某些情況下，氣體所儲存的熱能可轉化為機械能的形式而輸出能量做功。北京航空航天大學的高歌教授指出，目前熱機做功基本上都是通過膨脹做功的原理來實現，例如拉瓦爾噴管等。但在自然界中所發生的某些現象雖然沒有相應的做功機構，但仍然可以達到做功的目的，典型的例子就是龍卷風。顯然龍卷風氣體自身所具有的巨大熱能是龍卷風能夠做功的主要動力來源。故本文通過開發計算流體力學軟件CFL3D(cfl3d v5 manual)[15]給定了龍卷風初始速度場和邊界條件，匹配不同的初始溫度場以及龍卷風外圍大氣溫度情況，研究龍卷風風場的發展演化與溫度場的關系。

1 數值模擬設置

1.1 數值模擬幾何模型

本文龍卷風的初場是給定的，且主要研究的是龍卷風的發展演化，因此其發生區域給定在一尺寸為1 000 m×1 000 m×5 m的長方體內。模型的網格是采用軟件Pointwise生成的結構網格，因為參數變化主要集中在核心半徑區域，所以在中心加密，從內到外網格由密變疏。網格維數為193×193×5，網格俯視截面圖如圖1所示。

1.2 控制方程及計算方法

本文選用的控制方程為連續性方程、Navier-Stokes方程以及能量方程。三維非定常流體的連續性方程為

圖1 龍卷風數值模擬的網格劃分Fig.1 Mesh generation for numerical simulation of tornado

本次模擬采用可壓縮流動模擬，故應用于本文的Navier-Stokes方程和能量方程分別為

式中：ρR為單位體積流體所受的質量力，R為質量力；p為黏性流體的動壓；τij為黏性應力張量；e為單位質量流體的能量；λ為導熱系數；T為流體溫度；q為由于熱輻射或其他原因在單位時間內傳入的熱量。

計算方法方面，無黏項空間離散格式采用三階Roe格式+MUSCL插值，黏性通量空間離散格式為中心差分，時間推進方式采用近似因子分解算法,雙時間步推進實現精確地實時計算。

1.3 初始速度場及邊界條件的設置

龍卷風的初始速度場是根據文獻[16]中所提供的某種情況下龍卷風的速度剖面圖并查找龍卷風相關風強級別后加以修改的，如圖2所示。rc為初始速度場的渦核半徑(速度最大處的半徑)。在本文的模型中，初始速度場的半徑r為140 m，渦核半徑rc為50 m，所取最大速度為60 m/s，屬于F2級龍卷風速度范圍。

圖2 初始速度場的速度分布Fig.2 Velocity distribution of initial velocity field

本文忽略了龍卷風在水平方向上的移動，輸入文件中每一步的時間步長dt為0.1 s，每個步長內迭代5次，計算步數NTSTEP為100、1 000、2 000、4 000和8 000。在數值模擬中邊界條件的選擇直接關乎到模擬的可信度甚至成敗，由于本文研究的是龍卷風內部某一截平面上的變化，相當于是水平截取龍卷風中的一部分，故設置6個面的邊界條件均為無黏面(inviscid surface)。

1.4 不同初始溫度場的設置

1) 算例1初始溫度場為均一的溫度場，即場域中初始溫度均一樣，為288 K。

2) 算例2的初始溫度場采用類似于算例1相對穩定后的形式。

2 模擬結果與分析

在龍卷風的發展過程中，氣流的旋轉運動必然會改變局部氣流的參數變化乃至影響到全場。在這些參數中，本文著重關注速度場、溫度場、壓力場和密度場這4個參數場，通過對比分析來解析龍卷風的發展過程。速度場直觀地展現了場中機械動能的大致分布，顯示出了龍卷風破壞力最強的地方；溫度場的變化則是明確了內能轉化為動能的方向以及全場能量變化的趨勢；壓力場最大的作用是提供足夠的壓力來抵消氣流旋轉產生的離心力來維持龍卷風旋轉，而且壓力主要依托溫度和密度存在；密度場的變化主要體現了氣流的運動方向，是向某一方向匯聚還是發散，可以說明該狀態下龍卷風是被加強還是被削弱。

2.1 算 例 1

龍卷風的旋轉外形是漏斗狀的，其風速包含水平運動速度、豎直方向速度、徑向速度和切向速度。本文忽略龍卷風水平運動速度，又由于發展演化過程中主要是內外氣流的混合過程，所以主要考慮切向速度和徑向速度的合速度的變化，即龍卷風在平面上的速度變化。圖3給出了算例1隨時間推移在平面上的速度變化曲線。由于在均勻溫度場情況下，場內的變量變化并不劇烈，主要集中在初始速度型區域，所以只給出半徑為3rc內的變量曲線。由圖3可以看出，速度型隨時間雖然是呈整體降低的形式，但降低的程度很小，而且擴散不厲害，也就是說必然有額外的能量提供來抵消動能的耗散。在整個場域中，氣流除了有以動能為主的機械能外，就只有氣流本身所儲存的內能。而氣體的內能主要是以溫度的形式來體現的。

圖3 不同計算步數下的速度變化曲線(算例1)Fig.3 Velocity variation curves under different calculation steps (Case 1)

圖4給出了算例1在不同計算步數下的溫度變化曲線。從圖中可以直觀地看出，隨著計算步數的增加，即隨著時間的推移，中心溫度開始降低最終形成中心低溫區，而速度型邊緣氣流溫度的變化則是先升高降低到保持為288 K。整個溫度場從均一溫度場變為有一定梯度的不均一溫度場，顯然是能量的轉化導致了溫度的變化，氣流的內能在溫度降低時將部分能量釋放出來轉化成了動能。這是氣流以一定速度旋轉的結果，可以用漩渦管的Ranque-Hilsch效應來解釋，即當氣流在向心旋流器中，通過向心螺旋通道送入旋渦管產生旋流時，會發生溫度分離現象，旋渦管中外圍的氣流溫度比中心的氣流溫度高且該溫差隨初始氣流溫度的升高而升高。

圖5的壓力變化曲線和圖6的密度變化曲線應該結合圖4的溫度變化曲線分析，因為這3個參數必須滿足氣體狀態方程P=ρRT。從圖5和圖6中參數的變化情況來看，3個參數是符合該狀態方程的。從圖5中可以看出，，外圍高壓區和中心低壓區的形成向氣流提供了足夠的向心力來阻止擴散。從圖6中可以看出，中心密度是先急劇下降之后又回升，是因為剛開始是壓力場提供的向心力并不能夠抵消氣流旋轉所產生的離心力，因而中心區域氣流會先向外擴散，通過和外圍氣流擠壓換熱產生足夠強的壓力并回流之后，再形成比較穩定的密度場。

圖4 不同計算步數下的溫度變化曲線(算例1)Fig.4 Temperature variation curves under different calculation steps (Case 1)

圖5 不同計算步數下的壓力變化曲線(算例1)Fig.5 Pressure variation curves under different calculation steps (Case 1)

圖6 不同計算步數下的密度變化曲線(算例1)Fig.6 Density variation curves under different calculation steps (Case 1)

最終，氣流流場內的變量經過一系列的變化形成了一個各個參數互相匹配的新變量場，由此得到的新的旋渦場所具有的各個參數場的形式即為所有穩定的旋渦場所具有的共有特性。

2.2 算 例 2

由2.1節可知，算例1的數值模擬的理論基礎是Ranque-Hilsch效應，通過驗證可知漩渦場的旋轉會自發的產生一個與之匹配的溫度場來維持。而龍卷風也是一個比較特殊的漩渦場，故算例2在算例1的基礎上將均一的溫度場改為類似于算例1相對穩定后的溫度場形式，將其放入原初始速度場中計算并分析結果。這樣做的目的是觀察加大溫度梯度之后，龍卷風風場的變化。

算例2初始溫度場為中心低溫253 K，到初始速度場邊緣處線性升高為303 K，初始速度場外氣流溫度也為303 K。溫度場函數關系式為

在有溫度梯度的溫度場情況中，場內的變量變化會延伸到部分外圍氣流，所以給出的都是半徑為5rc內的變量曲線圖。

圖7給出了算例2隨時間推移的速度變化曲線，可以看出在2 000步的時候速度達到最大值，約為62.7 m/s。圖8將算例1中最大速度曲線同算例2中的速度曲線對比，可以發現算例2中速度曲線的各個速度最大值均大于算例1中出現的最大值，說明在算例2中速度在被加強。

從圖9中可以看出，龍卷風在旋轉過程中，中心溫度升高減小了內外氣流之間的溫度差，在2 000步時溫差達到最小，同時速度也達到了最大值。在整個變化過程中，溫度型都基本保持了相似的形態，僅僅在溫差上有變化，說明溫度型的選擇是正確的。在內外氣流之間匯流換熱的時候，外部高溫氣流的熱能輸運到內部低溫氣流使其溫度升高。此時，內部氣流的旋轉速度增大則表明外部高溫氣流提供給內部低溫氣流的熱能并非全部被用來增加內能，而是有部分熱能轉化為了動能加強了氣流的旋轉。這是由匯流換熱效應所引起的變化。在匯流換熱中，根據總壓靜壓的關系：

圖7 不同計算步數下的速度變化曲線(算例2)Fig.7 Velocity variation curves under different calculation steps (Case 2)

圖8 速度變化曲線對比Fig.8 Comparison of velocity variation curve

圖9 不同計算步數下的溫度變化曲線(算例2)Fig.9 Temperature variation curves under different calculation steps (Case 2)

(1)

由式(1)以及動量方程

dP=-ρVdV

可推導出總溫和總壓的關系式為

(2)

式中：Ma為氣流馬赫數；k為絕熱指數；P*和T*分別為總壓和總溫，總溫代表了熱能的大小，總壓代表的做功能力的大小。

由式(2)可知，總溫的降低會增大總壓，而且在此基礎上氣流馬赫數(速度)越高，其溫度變化帶來的總壓變化越大，由熱能儲存成的機械能增多。高速氣流降溫、低速氣流升溫和降溫帶來的總壓增強遠大于升溫帶來的總壓損失。顯然，圖8中速度的增大與圖9中溫度的變化驗證了這一點。

因此可以認為，大氣中的溫差及內外對流是龍卷旋渦形成的前提條件，而使龍卷旋渦得以持續甚至增強的能量則來源于周圍的熱能。龍卷風外圍的熱氣流與旋渦中心的冷氣流形成的有一定溫度梯度的溫度場，成為把大氣熱能轉化為旋渦流動動能的主要因素。自然界中，龍卷風自轉的同時還發生沿地表的平移，那么龍卷風就源源不斷地獲得環境的熱能，從而維持自身旋轉。

從圖10可以看出，此時算例2中提供的壓力要大于算例1中所匹配出來的壓力。從圖11中可以更加直觀的看出這一結果，對比圖11與圖6

圖10 不同計算步數下的壓力變化曲線(算例2)Fig.10 Pressure variation curves under different calculation steps (Case 2)

圖11 不同計算步數下的密度變化曲線(算例2)Fig.11 Density variation curves under different calculation steps (Case 2)

發現在相同半徑內，算例1中的密度變化曲線所對應的值大都在1以下，而算例2的密度變化曲線所對應的值則大都在1以上，顯然在相同半徑內算例2的密度積分要遠大于算例1的密度積分，也就是說在這一半徑內，算例2中包含有更多質量的氣流，且這些氣流的速度也很高，需要更大的壓力來抵制離心力。這些集中在一定范圍內的高速高質量的氣流所具有的高動能表明該龍卷風的強度被大大加強，進一步說明了算例2的初始溫度場加強了龍卷風的破壞能力。由此可以認為以這一溫度型為基礎，溫度梯度越大對龍卷風風場的加強作用越大。

3 結 論

1) 均一溫度場的旋渦在旋轉過程中會發生Ranque-Hilsch效應，產生一個不均一的溫度場，通過這個溫度場產生過程中氣流內能所轉化的動能來維持旋轉，與此同時引起壓力場和密度場的變化，最終這些變量互相匹配成一個新的旋渦變量場。由于新旋渦變量場是自適應生成的，所以其具有的變量場形式可代表所有穩定旋渦的變量場。

2) 特征如算例2的不均一的溫度場更加貼近于實際龍卷風的溫度場，其通過將氣流匯聚到核心區域的加大旋轉氣流質量的同時加強氣流旋轉速度，來強化龍卷風旋渦并加劇其破壞力，其是由匯流換熱原理所引起的。

3) 通過算例1的引申和算例2的驗證，得到龍卷風的生成和強化主要依托于其溫度場的變化，溫度場梯度的大小在一定程度上決定了該龍卷風的破壞能力。因此，為使龍卷風的強度迅速降低甚至消亡，破壞龍卷風的溫度場也許是一種行之有效的方法。

參考文獻 (References)

[1] HOECKER W H.Wind speed and air flow patterns in the Dallas tornado of April 2,1957[J].Monthly Weather Review,1960,88(5):167-180.

[2] ALEXANDER C R,WURMAN J.The 30 May 1998 Spencer,South Dakota,storm.Part Ⅰ:The structural evolution and environment of the tornadoes[J].Monthly Weather Review,2005,133(1):72-96.

[3] WURMAN J,ALEXANDER C R.The 30 May 1998 Spencer,South Dakota,storm.Part Ⅱ:Comparison of observed damage and radar-derived winds in the tornadoes[J].Monthly Weather Review,2005,133(1):97-119.

[4] CHANG C C.Tornado wind effects on buildings and structures with laboratory simulation[C]∥Proceedings of the Third International Conference on Wind Effects on Buildings and Structures,1971:231-240.

[5] WARD N B.The exploration of certain features of tornado dynamics using a laboratory model[J].Journal of the Atmospheric Sciences,1972,29(6):1194-1204.

[6] CHURCH C R,SNOW J T,AGEE E M.Tornado vortex simulation at Purdue University[J].Bulletin of the American Meteorological Society,1977,58(9):900-908.

[7] JISCHKE M C,PARANG M.Properties of simulated tornado-like vortices[J].Journal of the Atmospheric Sciences,1974,31(2):506-512.

[8] SARKAR P P,HAAN F L,GALLUS JR W A,et al.A laboratory tornado simulator:Comparison of laboratory,numerical and full-scale measurements[C]∥10th Americas Conference on Wind Engineering.Baton Rouge,LA：American Association for Wind Engineering,2005.

[9] HASSENZAHL H C.Numerical investigations of a tornado vortex using vorticity confinement[D].Madison:University of Wisconsin Madison,2007.

[10] MARUYAMA T.A numerically generated tornado-like vortex by large eddy simulation[C]∥Proceedings of 7th Asia Pacific Conference on Wind Engineering.Taipei:Wind Engineering,2009.

[11] MARUYAMA T.Simulation of flying debris using a numerically generated tornado-like vortex[J].Journal of Wind Engineering and Industrial Aerodynamics,2011,99(4):249-256.

[12] NATARAJAN D.Numerical simulation of tornado-like vortices[D].London:The University of Western Ontario,2011.

[13] NATARAJAN D,HANGAN H.Numerical study on the effects of surface roughness on tornado-like flows[C]∥11th Americas Conference on Wind Engineering(11ACWE).Baton Rouge,LA:American Association for Wind Engineering,2009.

[14] 徐楓,肖儀清,李波,等.龍卷風風場特性的 CFD 數值模擬[J].空氣動力學學報,2013,31(3):350-356.

XU F,XIAO Y Q,LI B,et al.Study of the adaptive optimal design method based on design variables space[J].Acta Aerodynamica Sinica,2013,31(3):350-356(in Chinese).

[15] KRIST S L,BIEDRON R T,RUMSEY C L.CFL3D user’s manual(version5.0)：NASA-TM-1998-208444[R].Washington,D.C.:NASA,1998.

[16] 王錦,周強,曹曙陽,等.龍卷風風場的試驗模擬[J].同濟大學學報(自然科學版),2014,42(11):1654-1659.

WANG J,ZHOU Q,CAO S Y,et al.Physical study on tornado-like flow based on tornado vortex simulator[J].Journal of Tongji University(Natural Science),2014,42(11):1654-1659(in Chinese).