徐讓書，葛 寧，王治敏

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

1 引言

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

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

文獻(xiàn)[3]系統(tǒng)的研究了泵與風(fēng)機(jī)，給出風(fēng)扇產(chǎn)生壓力的對應(yīng)關(guān)系；文獻(xiàn)[4]對內(nèi)燃機(jī)中的一維非定常流動(dòng)進(jìn)行了系統(tǒng)的分析，但沒有考慮動(dòng)力裝置提供的壓力；文獻(xiàn)[5]應(yīng)用計(jì)算機(jī)計(jì)算直流伺服電動(dòng)機(jī)的動(dòng)態(tài)特性曲線計(jì)算及其實(shí)驗(yàn)測定，給出了通過計(jì)算機(jī)模擬電機(jī)的計(jì)算方法；1996年，文獻(xiàn)[6]研究作用于海洋工程結(jié)結(jié)構(gòu)物的陣風(fēng)及其譜模擬實(shí)現(xiàn)了在系泊平臺模型試驗(yàn)中進(jìn)行了風(fēng)速的模擬；文獻(xiàn)[7]系統(tǒng)地介紹了國內(nèi)外低速風(fēng)洞動(dòng)態(tài)試驗(yàn)技術(shù)的新進(jìn)展，而其他文獻(xiàn)中鮮有論述。

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

2 風(fēng)洞系統(tǒng)物理模型

研究的回流式低速風(fēng)洞結(jié)構(gòu)，如圖1所示。

圖1 風(fēng)洞輪廓圖Fig.1 Wind Tunnel Contour

為控制方便，該風(fēng)洞的動(dòng)力裝置采用他勵(lì)式直流電機(jī)，只需控制電機(jī)的輸入電壓即可改變對氣流的作用力。

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

控制體截面的面積和周長是影響風(fēng)洞流動(dòng)的重要因素，其值隨風(fēng)洞位置變化，根據(jù)實(shí)際風(fēng)洞的輪廓，將原光滑輪廓曲線采用梯形簡化，可較精確地模擬出流場的信息。電機(jī)風(fēng)扇產(chǎn)生推力推動(dòng)風(fēng)洞中的流體運(yùn)動(dòng)，經(jīng)過蜂窩器、阻尼網(wǎng)、導(dǎo)流片等幾個(gè)元件以及壁面產(chǎn)生摩擦力阻礙空氣流動(dòng)，均可視作阻力元件，僅影響空氣流速?？煞謩e在該元件所在控制體單元賦予相應(yīng)的阻力系數(shù)實(shí)現(xiàn)對阻力元件的模擬。由于洞體較大，內(nèi)部摩擦力較小，故其溫度幾乎與外界環(huán)境相同。

3 風(fēng)洞系統(tǒng)數(shù)學(xué)模型

3.1 一維流動(dòng)數(shù)學(xué)模型

風(fēng)洞系統(tǒng)中的氣體流動(dòng)可視為考慮摩擦和可壓縮的一維非定常無黏氣體動(dòng)力學(xué)問題，它滿足相應(yīng)的質(zhì)量、動(dòng)量和能量守恒方程。能量轉(zhuǎn)化的過程是電機(jī)產(chǎn)生的推力對空氣做功，轉(zhuǎn)化為空氣的動(dòng)能以及空氣與阻力元件和壁面的耗散，沒有其它形式的能量轉(zhuǎn)化。風(fēng)洞內(nèi)部空氣可視為理想氣體，用理想氣體狀態(tài)方程化簡能量方程。在直角坐標(biāo)系下的一維流動(dòng)方程組的微分形式為：

式中：u—方程所要計(jì)算的物理量即壓強(qiáng)和速度；f—對應(yīng)沒有時(shí)間偏導(dǎo)數(shù)的項(xiàng)。根據(jù)一位非定常流動(dòng)的守恒原理可列出方程簡化整理后，得：

式中：G—源項(xiàng)，m2/s；ρ—空氣密度，kg/m3；v—流動(dòng)速度，m/s；A—風(fēng)洞截面積，m2；x—沿洞體切向方向坐標(biāo)位置，m；p—壓強(qiáng)，Pa。

將風(fēng)扇產(chǎn)生的壓力和阻力元件產(chǎn)生的摩擦均加在方程組的源項(xiàng)中：

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

第一項(xiàng)為阻力源項(xiàng)，不同控制體的摩擦阻力系數(shù)根據(jù)所在位置的阻力元件和壁面的摩擦可得；第二項(xiàng)為動(dòng)力源項(xiàng)，風(fēng)扇產(chǎn)生的壓強(qiáng)傳遞到氣流中，推動(dòng)風(fēng)扇所在的一個(gè)控制體單元做功，其他位置的P為0。

3.2 電機(jī)風(fēng)扇數(shù)學(xué)模型

風(fēng)洞中的動(dòng)力裝置通過改變輸入電壓來控制電機(jī)轉(zhuǎn)速從而是風(fēng)扇拍擊空氣產(chǎn)生推力。根據(jù)電機(jī)風(fēng)扇的動(dòng)態(tài)特性應(yīng)用基爾霍夫定律、達(dá)朗貝爾原理的耦合關(guān)系列出方程組，其形式如下：

式中：Ia—電機(jī)電樞電流，A；U—電機(jī)兩端電壓，V；Ra—電樞電阻，Ω；Gaf—電樞的運(yùn)動(dòng)電勢系數(shù)；If—?jiǎng)?lì)磁繞組電阻，Ω；ω—電樞機(jī)械角速度，rad/s；La—電樞電感；Rω—摩擦系數(shù)；J—風(fēng)扇轉(zhuǎn)動(dòng)慣量，kg·m2。

本試驗(yàn)采用直流他勵(lì)式電機(jī)，其中，Ra、Gaf、If、La、Rω不隨時(shí)間變化，通過改變輸入電壓來計(jì)算出電樞機(jī)械角速度，并換算出轉(zhuǎn)速。再通過泵與風(fēng)機(jī)中風(fēng)扇轉(zhuǎn)速與產(chǎn)生速度之間的關(guān)系式，如圖2所示。得到風(fēng)扇處瞬時(shí)提供的壓力值。

圖2 風(fēng)扇轉(zhuǎn)速關(guān)系圖Fig.2 Fan Speed Diagram

圖中：H—揚(yáng)程，可根據(jù)物理關(guān)系轉(zhuǎn)化成壓強(qiáng)，m；Q—空氣流量，可計(jì)算出流體速度，m3/s；η—風(fēng)扇轉(zhuǎn)速，與電機(jī)中ω值相同，rad/s。

進(jìn)過分析，上述方程組封閉可以求解。先通過計(jì)算一維流動(dòng)模型得到達(dá)到該風(fēng)速在風(fēng)扇處所需的壓強(qiáng)，再通過風(fēng)扇壓強(qiáng)與轉(zhuǎn)速的關(guān)系得到電機(jī)轉(zhuǎn)速，通過電機(jī)模型用該轉(zhuǎn)速最終計(jì)算得到該時(shí)刻輸入的功率，實(shí)現(xiàn)通過控制輸入功率的方式得到大氣邊界層所需的非定常風(fēng)。

4 風(fēng)洞系統(tǒng)模型求解

4.1 一維流動(dòng)模型求解

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

式中：f—對應(yīng)沒有時(shí)間偏導(dǎo)數(shù)的項(xiàng)；t—時(shí)間步長；下角標(biāo)j—該控制體所在位置；j-1—相鄰的左側(cè)控制體；j+1—相鄰的右側(cè)控制體；上標(biāo)n—當(dāng)前時(shí)刻預(yù)測值作為中間量；n+1—矯正值即下一時(shí)刻的物理量。

通過上述計(jì)算模型可將原偏微分方程組經(jīng)過數(shù)值方法直觀地從當(dāng)前時(shí)刻物理量計(jì)算出下一時(shí)刻的物理量。這樣，可應(yīng)用計(jì)算機(jī)編程實(shí)現(xiàn)迭代計(jì)算每一時(shí)刻的物理量并畫出相應(yīng)曲線。

4.2 電機(jī)風(fēng)扇模型求解

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

代入可得：

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

圖3 計(jì)算模型流程圖Fig.3 Flow Chart of Calculation Model

5 結(jié)果與分析

本次試驗(yàn)在北京三十五中學(xué)風(fēng)洞完成，其風(fēng)洞主要參數(shù)：實(shí)驗(yàn)段高度1m，寬度1.2m，水力直徑1.140m，當(dāng)?shù)貕毫闃?biāo)準(zhǔn)大氣壓力，溫度為25℃，電機(jī)額定功率為75kW，額定轉(zhuǎn)速為750r/min，效率為85.1%，動(dòng)力段直徑為2m，面積為3.1m2，風(fēng)扇葉片展長0.445m，4個(gè)拐角摩擦阻力系數(shù)0.18，蜂窩器阻力系數(shù)0.336，阻尼網(wǎng)阻力系數(shù)1.642。本次試驗(yàn)現(xiàn)場，如圖4所示。

圖4 北京35中學(xué)風(fēng)洞實(shí)驗(yàn)室Fig.4 Wind Tunnel Laboratory of Beijing 35 Middle School

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

通過上述測速方法對試驗(yàn)的風(fēng)洞進(jìn)行控制，設(shè)定穩(wěn)定風(fēng)速為20m/s，當(dāng)風(fēng)速通過負(fù)反饋調(diào)節(jié)PID控制達(dá)到穩(wěn)定后，停止PID控制功能，輸入在原電壓基礎(chǔ)上增加周期為25s振幅為0.3當(dāng)前測得的輸入電壓的正弦電壓，結(jié)果，如圖5所示。

圖5試驗(yàn)風(fēng)速變化曲線Fig.5 Test of Wind Speed Change Curve

圖5 結(jié)果表明電機(jī)經(jīng)PID控制在48s左右達(dá)到穩(wěn)定，穩(wěn)定后實(shí)驗(yàn)段風(fēng)速以正弦波形式隨時(shí)間改變。

應(yīng)用Matlab編程計(jì)算，其中PID控制部分不是研究重點(diǎn)，用線性代替，當(dāng)速度達(dá)到設(shè)定值時(shí)輸入與試驗(yàn)相同的信號，得到結(jié)果，如圖6所示。

圖6 模擬風(fēng)速變化曲線Fig.6 Simulation of Wind Speed Change Curve

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

圖7 試驗(yàn)結(jié)果頻域圖Fig.7 Test Results in Frequency Domain

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

6 結(jié)論

以上給出的通過計(jì)算機(jī)計(jì)算方法較好的模擬出低速風(fēng)洞非定常風(fēng)，模擬結(jié)果與試驗(yàn)對照表明：

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

（2）采用MacCormack二步差分的改進(jìn)方法適用于一維流動(dòng)偏微分方程組的求解，與電機(jī)風(fēng)扇模型耦合求解計(jì)算出滿足給定陣風(fēng)譜的風(fēng)速曲線，與實(shí)驗(yàn)結(jié)果基本吻合，驗(yàn)證了計(jì)算的準(zhǔn)確性。

（3）通過改變輸入功率產(chǎn)生的非定常風(fēng)與該功率在頻譜上的分布基本相同，可通過標(biāo)定的方法改變功率在風(fēng)洞中實(shí)現(xiàn)非定常風(fēng)。