低速風洞非定常風特性模型及驗證

徐讓書，葛 寧，王治敏

（沈陽航空航天大學 航空航天工程學部，遼寧 沈陽 110136）

1 引言

風洞是以人工的方式產生并且控制氣流，用來模擬飛行器或實體周圍氣體的流動情況，并可量度氣流對實體的作用效果以及觀察物理現象的一種管道狀實驗設備，它是進行空氣動力實驗最常用、最有效的工具之一。低速風洞是指試驗段氣流速度小于馬赫數小于0.3的風洞。從20世紀四十年代開始，國外開始建造低速風洞，迄今低速風洞已經得到了很大的發展[1]。

在大氣邊界層的紊流風能夠激起低頻運動發生共振現象，造成建筑結構破壞。故為研究此破壞程度，需要在風洞中模擬出建筑所在環境，即非定常的大氣邊界層風洞。采用回流式風洞，回流式低速風洞試驗中，風速往往在試驗前通過標定的手段來確定[2]，但只能實現產生穩定的風速，而不能實現對隨時間變化風的模擬。為了研究由于風速脈動導致的建筑結構破壞現象，從而達到建筑風洞對非定常風的相關要求（實現給定陣風譜在時域上的風速），需要實現在風洞試驗段中產生的功率譜滿足某一給定函數，經過反傅里葉變換后在時域上的值是隨時間變化的氣流，來模擬出給定陣風譜的風速。

文獻[3]系統的研究了泵與風機，給出風扇產生壓力的對應關系；文獻[4]對內燃機中的一維非定常流動進行了系統的分析，但沒有考慮動力裝置提供的壓力；文獻[5]應用計算機計算直流伺服電動機的動態特性曲線計算及其實驗測定，給出了通過計算機模擬電機的計算方法；1996年，文獻[6]研究作用于海洋工程結結構物的陣風及其譜模擬實現了在系泊平臺模型試驗中進行了風速的模擬；文獻[7]系統地介紹了國內外低速風洞動態試驗技術的新進展，而其他文獻中鮮有論述。

指在通過改變風洞動力裝置的輸入功率來實現對非定常風的模擬。風實際流動的物理現象是十分復雜的，根據實際情況需對其建立合理簡化的物理模型。將物理模型轉化數學模型，這就既需要以風場中的流體為研究對象建立起相關的運動方程，又要對風洞的動力裝置建立守恒方程，并對封閉的方程組耦合求解。方程組是非線性的微分方程組，無法直接得到解析解，所以通過數值方法得到數值解，從而實現通過風洞對非定長風的模擬，并進行試驗驗證。

2 風洞系統物理模型

研究的回流式低速風洞結構，如圖1所示。

圖1 風洞輪廓圖Fig.1 Wind Tunnel Contour

為控制方便，該風洞的動力裝置采用他勵式直流電機，只需控制電機的輸入電壓即可改變對氣流的作用力。

根據風洞的實際輪廓可知，由于風洞設計要求湍流度較小，在非拐角處其風速方向主要是沿著壁面的切線方向，而其他方向的風速幾乎為0。回流式風洞有4個拐角，但經過導流片導流后流動方向依然是壁面切線方向。風洞中沿軸向每一截面的物理量可視作相同，故可將其流動簡化為一維流動。在數值計算中，沿流動方向合理分段后，每一段可作為一個控制體。

控制體截面的面積和周長是影響風洞流動的重要因素，其值隨風洞位置變化，根據實際風洞的輪廓，將原光滑輪廓曲線采用梯形簡化，可較精確地模擬出流場的信息。電機風扇產生推力推動風洞中的流體運動，經過蜂窩器、阻尼網、導流片等幾個元件以及壁面產生摩擦力阻礙空氣流動，均可視作阻力元件，僅影響空氣流速。可分別在該元件所在控制體單元賦予相應的阻力系數實現對阻力元件的模擬。由于洞體較大，內部摩擦力較小，故其溫度幾乎與外界環境相同。

3 風洞系統數學模型

3.1 一維流動數學模型

風洞系統中的氣體流動可視為考慮摩擦和可壓縮的一維非定常無黏氣體動力學問題，它滿足相應的質量、動量和能量守恒方程。能量轉化的過程是電機產生的推力對空氣做功，轉化為空氣的動能以及空氣與阻力元件和壁面的耗散，沒有其它形式的能量轉化。風洞內部空氣可視為理想氣體，用理想氣體狀態方程化簡能量方程。在直角坐標系下的一維流動方程組的微分形式為：

式中：u—方程所要計算的物理量即壓強和速度；f—對應沒有時間偏導數的項。根據一位非定常流動的守恒原理可列出方程簡化整理后，得：

式中：G—源項，m2/s；ρ—空氣密度，kg/m3；v—流動速度，m/s；A—風洞截面積，m2；x—沿洞體切向方向坐標位置，m；p—壓強，Pa。

將風扇產生的壓力和阻力元件產生的摩擦均加在方程組的源項中：

式中：f—摩擦阻力系數；P—風扇作用于流體的壓強，Pa；c—風洞截面周長，m。

第一項為阻力源項，不同控制體的摩擦阻力系數根據所在位置的阻力元件和壁面的摩擦可得；第二項為動力源項，風扇產生的壓強傳遞到氣流中，推動風扇所在的一個控制體單元做功，其他位置的P為0。

3.2 電機風扇數學模型

風洞中的動力裝置通過改變輸入電壓來控制電機轉速從而是風扇拍擊空氣產生推力。根據電機風扇的動態特性應用基爾霍夫定律、達朗貝爾原理的耦合關系列出方程組，其形式如下：

式中：Ia—電機電樞電流，A；U—電機兩端電壓，V；Ra—電樞電阻，Ω；Gaf—電樞的運動電勢系數；If—勵磁繞組電阻，Ω；ω—電樞機械角速度，rad/s；La—電樞電感；Rω—摩擦系數；J—風扇轉動慣量，kg·m2。

本試驗采用直流他勵式電機，其中，Ra、Gaf、If、La、Rω不隨時間變化，通過改變輸入電壓來計算出電樞機械角速度，并換算出轉速。再通過泵與風機中風扇轉速與產生速度之間的關系式，如圖2所示。得到風扇處瞬時提供的壓力值。

圖2 風扇轉速關系圖Fig.2 Fan Speed Diagram

圖中：H—揚程，可根據物理關系轉化成壓強，m；Q—空氣流量，可計算出流體速度，m3/s；η—風扇轉速，與電機中ω值相同，rad/s。

進過分析，上述方程組封閉可以求解。先通過計算一維流動模型得到達到該風速在風扇處所需的壓強，再通過風扇壓強與轉速的關系得到電機轉速，通過電機模型用該轉速最終計算得到該時刻輸入的功率，實現通過控制輸入功率的方式得到大氣邊界層所需的非定常風。

4 風洞系統模型求解

4.1 一維流動模型求解

根據上述方程組可知，其形式比較復雜，無法得到精確的解析解，因此采用數值迭代方法耦合求得數值解。由于以流體為研究對象所列的方程組形式是三個未知數的偏微分方程，與一維Euler方程組類似，但形式不完全相同。根據分析可借鑒數值計算方法中具有二階精度的MacCormack兩步差分迭代法[8]，但由于方程組中存在源項不能直接應用。根據所求方程組的形式對該計算方法進行進一步改進，將原差分格式中的分母項還原中既有并迭代求解。求解流體運動方程的兩步差分迭代格式：

式中：f—對應沒有時間偏導數的項；t—時間步長；下角標j—該控制體所在位置；j-1—相鄰的左側控制體；j+1—相鄰的右側控制體；上標n—當前時刻預測值作為中間量；n+1—矯正值即下一時刻的物理量。

通過上述計算模型可將原偏微分方程組經過數值方法直觀地從當前時刻物理量計算出下一時刻的物理量。這樣，可應用計算機編程實現迭代計算每一時刻的物理量并畫出相應曲線。

4.2 電機風扇模型求解

求解風扇模型的守恒方程的4階龍格-庫塔公式解方程組形式：

代入可得：

應用4階經典龍格-庫塔可對該一階全微分方程組進行求解，得到下一時刻的角速度和電樞電流。角速度轉化的轉速通過其與壓強的關系得到風扇對流體產生的壓強，既流體源項中下一時刻P的值，并繼續對下一時刻其他物理量進行迭代計算。當達到設定時間時，停止計算。應用Matlab軟件進行計算，輸入初始各個位置測定的靜壓，溫度，空氣密度，阻力元件摩擦系數，通過計算得到時間步長和空間步長；將上述一維流動模型和電機風扇模型通過編程形式在程序中實現，即可通過初始值計算下一時刻每個位置物理量。迭代計算到給定時間后結束計算。該計算模型流程圖，如圖3所示。

圖3 計算模型流程圖Fig.3 Flow Chart of Calculation Model

5 結果與分析

本次試驗在北京三十五中學風洞完成，其風洞主要參數：實驗段高度1m，寬度1.2m，水力直徑1.140m，當地壓力為標準大氣壓力，溫度為25℃，電機額定功率為75kW，額定轉速為750r/min，效率為85.1%，動力段直徑為2m，面積為3.1m2，風扇葉片展長0.445m，4個拐角摩擦阻力系數0.18，蜂窩器阻力系數0.336，阻尼網阻力系數1.642。本次試驗現場，如圖4所示。

圖4 北京35中學風洞實驗室Fig.4 Wind Tunnel Laboratory of Beijing 35 Middle School

試驗段的檢測方法為標定法。在試驗開始前，在試驗段坐標支架上固定測速儀，使測速儀的方向平行于洞壁；為保證試驗段不受干擾，在風洞收縮段前的洞壁上分別打四個孔，四個孔的位置在同一截面上，并通過軟管連接到測壓儀，其作用是在收縮段前測量平均總壓。標定時先用測壓儀測得當地靜壓，啟動電機，讀取測速儀測得不同風速下對應的收縮段前總壓。根據伯努利方程可知動壓與風速的關系，但由于本次試驗測得壓力與速度的位置不同，所以必須乘上標定系數才能使等式成立。于是在不同的風速下測得對應的總壓以及標定系數存入計算機中。試驗時，去掉測速儀，啟動電機后不斷采集總壓值，并提取對應的標定系數，根據公式可計算出對應的風速。

通過上述測速方法對試驗的風洞進行控制，設定穩定風速為20m/s，當風速通過負反饋調節PID控制達到穩定后，停止PID控制功能，輸入在原電壓基礎上增加周期為25s振幅為0.3當前測得的輸入電壓的正弦電壓，結果，如圖5所示。

圖5試驗風速變化曲線Fig.5 Test of Wind Speed Change Curve

圖5 結果表明電機經PID控制在48s左右達到穩定，穩定后實驗段風速以正弦波形式隨時間改變。

應用Matlab編程計算，其中PID控制部分不是研究重點，用線性代替，當速度達到設定值時輸入與試驗相同的信號，得到結果，如圖6所示。

圖6 模擬風速變化曲線Fig.6 Simulation of Wind Speed Change Curve

從圖6中可以看出當輸入正弦信號后，風速的響應也是正弦曲線，與輸入信號有1s左右延遲。將試驗的速度曲線在啟動50s后經快速傅里葉變換后得到頻域曲線，如圖7所示。

圖7 試驗結果頻域圖Fig.7 Test Results in Frequency Domain

在圖7中表明，試驗的非定常風速的頻率與輸入正弦電壓頻率一致，振幅與風洞對應輸入功率所標定的風速一致；模擬結果做傅里葉變換后與該試驗結果相符，經計算其相對誤差約為2.4%。風速與輸入功率頻域上其主頻率相同，風速頻譜上主頻率附近有微小擾動應為測量誤差不影響實驗結果，可忽略不計。

6 結論

以上給出的通過計算機計算方法較好的模擬出低速風洞非定常風，模擬結果與試驗對照表明：

（1）將在回流式風洞中采用控制輸入功率的方式對陣風譜的模擬應用一維非定常的物理模型進行簡化，可大量減少計算量也能保證計算準確性。

（2）采用MacCormack二步差分的改進方法適用于一維流動偏微分方程組的求解，與電機風扇模型耦合求解計算出滿足給定陣風譜的風速曲線，與實驗結果基本吻合，驗證了計算的準確性。

（3）通過改變輸入功率產生的非定常風與該功率在頻譜上的分布基本相同，可通過標定的方法改變功率在風洞中實現非定常風。