高興龍，王超，符澄,*，孫運強，寇杰

1.空氣動力學(xué)國家重點實驗室，綿陽 621000

2.中國空氣動力研究與發(fā)展中心 設(shè)備設(shè)計與測試技術(shù)研究所，綿陽 621000

0 引 言

在航空航天領(lǐng)域，各類先進軍/民用航空航天飛行器的發(fā)展必然給國防安全和社會經(jīng)濟帶來巨大變化，同時也帶來了大量新的空氣動力學(xué)問題[1-2]，如高速氣動力精確預(yù)測及一體化設(shè)計、湍流減阻、邊界層轉(zhuǎn)捩、地面效應(yīng)等基礎(chǔ)和前沿問題，這些問題的研究都強烈依賴風(fēng)洞試驗，且對風(fēng)洞流場動態(tài)特性、低擾動特性及特種模擬能力要求很高，傳統(tǒng)風(fēng)洞設(shè)備越來越難以滿足日益特殊的空氣動力試驗需求[3-4]。在先進軌道交通領(lǐng)域，磁浮真空管道超高速列車等先進軌道交通技術(shù)被納入我國“十三五”重點研發(fā)計劃。真空管道列車面臨著復(fù)雜的激波邊界層干擾、超聲速氣動減阻、管道激波反射、活塞效應(yīng)等一系列復(fù)雜的空氣動力學(xué)問題，必須在風(fēng)洞中加以解決[5]。目前，傳統(tǒng)跨超聲速風(fēng)洞存在試驗段尺寸（長度和截面積）和氣流相對運動帶來的模擬真實性不足等問題，難以在受限空間內(nèi)開展航空航天飛行器和超高速列車空氣動力特性及力/熱/結(jié)構(gòu)/控制耦合問題研究[6-9]。

“飛行風(fēng)洞”的概念最早由美國研究者于20 世紀(jì)90 年代提出，具備“體動風(fēng)靜”的特殊運行方式和性能優(yōu)勢。NASA 的蘭利研究中心在較早前就開始關(guān)注這一全新概念的地面試驗設(shè)備，稱之為“高升力飛行風(fēng)洞（high-lift flight tunnel）”[10]，并開展了3 期關(guān)鍵技術(shù)研究和論證規(guī)劃工作。由于該設(shè)備的關(guān)鍵技術(shù)攻關(guān)難度較高，至今仍未取得較明顯的進展?！按鸥★w行風(fēng)洞”是利用真空管道列車概念結(jié)合動模型試驗技術(shù)提出的一種新概念風(fēng)洞設(shè)備，其原理是在一段封閉的直線長管道內(nèi)安裝磁浮模型驅(qū)動機構(gòu)，利用電磁懸浮、牽引和導(dǎo)向技術(shù)驅(qū)動模型高速運動，模擬各類飛行器及高速列車運動的物理過程，構(gòu)建接近真實飛行環(huán)境和運動特點的“體動風(fēng)靜”試驗狀態(tài)，可以滿足航空航天飛行器、超高速列車等寬馬赫數(shù)范圍、寬雷諾數(shù)范圍、低噪聲、低湍流度、高真空度（高空）、特殊氣體介質(zhì)、受限空間條件下的空氣動力學(xué)及其交叉學(xué)科地面試驗需求[3,5]。

磁浮飛行風(fēng)洞采用高速磁浮技術(shù)驅(qū)動試驗?zāi)Ｐ驮诘冉孛?、直線、密閉管道內(nèi)做加速、勻速和減速運動，可以模擬高加速/高減速（高過載）過程及速度突變時的氣動現(xiàn)象。試驗過程中，風(fēng)洞內(nèi)模型高速運動所產(chǎn)生的氣動力和聲波的傳播十分復(fù)雜，同時，受擾動的流體與結(jié)構(gòu)之間還存在相互作用，屬于氣動結(jié)構(gòu)交叉學(xué)科之間的耦合問題[11]。本文提到的氣動結(jié)構(gòu)耦合仿真，與傳統(tǒng)流固耦合存在一定區(qū)別。傳統(tǒng)流固耦合概念為固體力學(xué)研究領(lǐng)域所提出，更側(cè)重于結(jié)構(gòu)的動態(tài)響應(yīng)問題，而本文則更側(cè)重于管道內(nèi)模型高加速/高減速和勻速運行過程的非定常氣動特性求解，并不存在結(jié)構(gòu)變形問題；但同時也必須考慮動模型的高加速/高減速情況與流場之間的耦合作用，傳統(tǒng)CFD 仿真對于高加速/高減速的流場動態(tài)響應(yīng)求解還存在一定局限。此外，本文研究的方法更多是為后期更深入地研究高過載條件下磁浮平臺的結(jié)構(gòu)振動和運動部件的功能及安全性評估等問題提供理論基礎(chǔ)。傳統(tǒng)CFD 仿真方法在此類動模型高速運動耦合問題的處理上還存在一定局限性，須結(jié)合動網(wǎng)格等技術(shù)，利用新型氣動結(jié)構(gòu)耦合仿真方法對動模型試驗過程進行仿真，從氣動結(jié)構(gòu)耦合的角度研究模型在密閉管道內(nèi)高速運動過程的氣動特性問題，為磁浮飛行風(fēng)洞流動控制系統(tǒng)的設(shè)計提供技術(shù)支撐。

時空守恒元和解元方法（Space Time Conservation Element and Solution Element Method,CE/SE）是近年來興起的一種全新的高分辨率守恒型方程數(shù)值積分方法，由Chang[12]于1995 年提出。該方法將時間和空間統(tǒng)一起來同等對待，利用守恒型積分方程，通過定義守恒元（Conservation Element,CE）和解元（Solution Element,SE）使局部和整體都嚴(yán)格滿足守恒律[13]。經(jīng)過開發(fā)，該方法已成為一種新型的可求解二維和三維非穩(wěn)態(tài)歐拉方程的可壓縮體求解方法，并被應(yīng)用于沖擊波、氣動噪聲、磁流體力學(xué)等典型多物理場耦合問題[14-19]。

與傳統(tǒng)“體靜風(fēng)動”風(fēng)洞試驗技術(shù)相比，“體動風(fēng)靜”風(fēng)洞試驗技術(shù)主要基于動模型試驗技術(shù)原理，依靠動力驅(qū)動模型在滑軌上運動，從而開展相關(guān)的空氣動力學(xué)測試[3]。國內(nèi)外在高速動模型試驗技術(shù)與計算模擬方面開展了大量的研究。日本Kobayasi 物理研究所建立了小型高速列車壓力波測量裝置，德國宇航中心的哥廷根隧道模擬試驗平臺可以開展靜態(tài)和瞬態(tài)列車空氣動力學(xué)研究，韓國建筑技術(shù)研究院建立了電動機帶傳動彈射的小型列車動模型試驗平臺。國內(nèi)開展動模型試驗技術(shù)研究的主要有中南大學(xué)軌道交通安全重點實驗室、中國科學(xué)院先進軌道交通力學(xué)研究中心和西南交通大學(xué)試驗中心等。在氣動結(jié)構(gòu)計算方面，中南大學(xué)王樂卿和高廣軍等[3]通過仿真研究了高速列車車身風(fēng)阻制動板氣動外形設(shè)計問題；西南交通大學(xué)李一凡等[4]利用CFD方法研究了導(dǎo)流裝置對高速磁浮列車氣動特性的影響。

本文基于CE/SE 方法和氣動結(jié)構(gòu)分區(qū)耦合技術(shù)，根據(jù)磁浮飛行風(fēng)洞的實際試驗需求，合理建立仿真工況，從氣動結(jié)構(gòu)耦合的角度對動模型在磁浮飛行風(fēng)洞內(nèi)的加速、勻速和減速全過程進行氣動特性仿真評估，分析動模型在密閉管道內(nèi)高速運動的流固耦合問題，獲得受擾動流場各參數(shù)的時間歷程曲線，為磁浮飛行風(fēng)洞流動控制系統(tǒng)和消波措施等的方案設(shè)計提供支撐。

1 問題描述與方法

對磁浮飛行風(fēng)洞的磁浮平臺和軌道模型進行適當(dāng)簡化，計算運行過程中平臺及模型所受到的氣動力。簡化后的管道、軌道橫截面和平臺的模型如圖1 所示，其中，z 為模型運動方向，x 為管道橫截面的水平方向，y 為管道橫截面垂直方向；管道截面可以沿軸向直接拉伸形成三維管道模型，管道的軸向總長度為磁浮飛行風(fēng)洞的長度。

圖1 磁浮飛行風(fēng)洞磁浮平臺及管道截面組合模型圖Fig.1 Combined model between maglev platform and tunnel section of maglev flight tunnel

1.1 流體求解方法

首先給出一維的波傳播方程：

式中：a ＞0，為對流項常數(shù)；u 為波傳播速度；t 為時間。該方程令單元空間離散，在不同空間位置和給定時間產(chǎn)生離散點。CE/SE 方法通過構(gòu)建二維歐拉空間域 E2將時間作為附加空間坐標(biāo)考慮，在 E2內(nèi)的時間和空間上應(yīng)用高斯散度定理，從而將微分形式的式（1）改寫為積分形式[13]：

式中：S(V)為 E2內(nèi)任意時間-空間域 V 的邊界；=(au,u)；=dσ，dσ和分別為S(V)表面單元外法向的面積和單位向量。在 E2內(nèi)建立一個單元體積，該體積內(nèi)的單元在時間和空間上守恒且被統(tǒng)一處理，這是構(gòu)建守恒元CE 的關(guān)鍵。之后建立SE，保證SE 單元內(nèi)的變量足夠小，從而可以在 E2內(nèi)坐標(biāo)上某點中心的時空區(qū)域內(nèi)部邊界附近對流體變量進行近似泰勒級數(shù)展開：

式中：u*為SE 單元內(nèi)部流體變量；uq為時空區(qū)域守恒元q 點的流體變量；x 為空間坐標(biāo)，xq為q 點的空間坐標(biāo)，tq為q 點的時間坐標(biāo)。通過流體的對流擴散方程將時間和空間導(dǎo)數(shù)相關(guān)，從而僅存在uq(x,t)和它的空間導(dǎo)數(shù)這兩個未知量。

為了令系統(tǒng)封閉，定義2 個CE 的方程，之后獲取積分形式的對流方程，如公式（2）所示，則CE/SE格式可以保證SE 在時間和空間上統(tǒng)一守恒。為保證流體守恒，沿CE 所形成的直線域進行時間和空間積分：

1.2 氣動結(jié)構(gòu)耦合平臺

對于氣動結(jié)構(gòu)耦合問題，須將流體求解器和結(jié)構(gòu)求解器進行一致求解，尋找信息傳遞變量，建立耦合平臺。這里可以將CE/SE 方法和結(jié)構(gòu)動力學(xué)方法進行耦合，以實現(xiàn)高速運動過程的氣動結(jié)構(gòu)耦合計算仿真。上述2 個方法可以獨立進行計算，因此可以分別對流場和結(jié)構(gòu)進行網(wǎng)格劃分和時間步長設(shè)置，耦合系統(tǒng)自動追蹤2 個域中的最小時間步長進行計算。氣動結(jié)構(gòu)耦合須將拉格朗日結(jié)構(gòu)嵌入到流體域中，并在每個時間步內(nèi)將節(jié)點位移和速度信息作為向流體求解器傳遞的變量。在初始求解時間步內(nèi)，CE/SE 方法首先用初始速度和壓力求解可壓縮流場，獲得界面壓力；之后將壓力作為邊界條件施加到結(jié)構(gòu)域上，對結(jié)構(gòu)域進行求解后，得到節(jié)點的位移和速度；再將邊界節(jié)點上的位移和速度返回到流體場，從而使流體和結(jié)構(gòu)進入收斂迭代的過程；當(dāng)滿足收斂條件后，再進入下一個時間步的計算。整個系統(tǒng)的流固耦合迭代過程如圖2 和3 所示。

圖2 流固耦合平臺Fig.2 Fluid structure interaction platform

圖3 基于CE/SE 算法的流固耦合求解流程圖[20]Fig.3 Solving flowchart of FSI by CE/SE algorithm[20]

1.3 仿真模型及參數(shù)設(shè)置

結(jié)構(gòu)有限元模型主要由磁浮飛行風(fēng)洞磁浮平臺和試驗?zāi)Ｐ徒M成。流體模型如圖4 所示。流體網(wǎng)格域的長度為1 000 m。根據(jù)CE/SE 方法的特點，流場采用結(jié)構(gòu)化網(wǎng)格，最小網(wǎng)格尺寸為0.01 m，結(jié)構(gòu)和流體網(wǎng)格共計約1 200 萬單元。結(jié)構(gòu)域和流體域網(wǎng)格分別如圖5 和6 所示。

圖4 有限元流體模型示意圖Fig.4 Layout of finite element flow model

圖5 磁浮平臺和試驗?zāi)Ｐ徒Y(jié)構(gòu)域網(wǎng)格圖Fig.5 Structure mesh of maglev platform and test model

圖6 管道橫截面流體域網(wǎng)格圖Fig.6 Flow mesh of tunnel section

仿真模型采用的長度單位為m，時間單位為s，質(zhì)量單位為kg。平臺和結(jié)構(gòu)模型假設(shè)為不變形，材料采用剛體材料。

2 仿真結(jié)果及評估分析

2.1 仿真工況

磁浮飛行風(fēng)洞管道總長度為 1000 m，分為加速段、勻速段和減速段，根據(jù)試驗條件建立2 個仿真工況：工況1 的最大速度為170 m/s，工況2 的最大速度為340 m/s，加速段和減速段為勻加速和勻減速。工況參數(shù)如表1 所示，所有的仿真工況均考慮重力因素。仿真模型模擬的試驗阻塞比為20%，即模型運動過程中阻力面積占管道橫截面可用面積的20%。

表1 仿真工況參數(shù)Table 1 List of simulation conditions

2.2 氣動結(jié)果分析評估

工況1、2 模型運行速度隨時間變化曲線如圖7、8 所示。取風(fēng)洞管道內(nèi)部縱向位置為0、450、500、600 和1 000 m 的5 個橫截面，輸出密閉管道內(nèi)橫截面所受壓力隨時間變化的曲線，工況1 各橫截面的結(jié)果如圖9 所示。從450、500 和600 m 橫截面的情況可以看出：這3 個橫截面均經(jīng)歷了2 個壓力波峰值，第一個壓力波峰出現(xiàn)的時刻略早于模型到達所選橫截面的時刻，說明模型頭部壓縮波的運動速度明顯快于模型本身；第一個壓力波達到峰值后馬上急劇下降，之后緩慢增大，這應(yīng)該是模型頭部產(chǎn)生的壓縮波運動經(jīng)過所選橫截面后，尾部產(chǎn)生的膨脹波又出現(xiàn)在所選橫截面導(dǎo)致的；第二個壓力波峰出現(xiàn)時，模型處于減速階段，可能是壓縮波遇到管道末端壁面產(chǎn)生的反射波重新運動到所選橫截面位置導(dǎo)致的。1 000 m 處管道末端橫截面只產(chǎn)生了一個波峰，在模型減速運動階段出現(xiàn)。

圖7 工況1 模型運行速度隨時間變化曲線Fig.7 Model moving velocity vs.time for condition 1

圖8 工況2 模型運行速度隨時間變化曲線Fig.8 Model moving velocity vs.time for condition 2

圖9 工況1 管道內(nèi)各橫截面壓力變化曲線Fig.9 Pressure vs.time curves of tunnel sections for condition 1

工況2 各橫截面的結(jié)果如圖10 所示。從圖10可以看出：450、500 和600 m 橫截面處的情況與工況1 類似，均出現(xiàn)了了模型前端的壓縮波和反射波導(dǎo)致的2 次峰值。但工況2 管道末端1 000 m 橫截面處出現(xiàn)了多個峰值，應(yīng)該是壓縮波在末端壁面多次反射的結(jié)果。

圖10 工況2 管道內(nèi)各橫截面的壓力變化曲線Fig.10 Pressure vs.time curves of tunnel sections for condition 2

將2 個工況壓力隨時間變化的曲線與模型運動的速度變化曲線進行對比，可以看出：模型勻速運行過程中，試驗段橫截面上的壓力波均出現(xiàn)了波動峰值，不利于氣動試驗的開展，須采取一定的消波措施削弱壓力波。

圖11 為工況1 模型運動過程中所受氣動阻力隨時間變化的曲線。從圖中可以看出：模型所受氣動阻力在加速段與勻速段相差不大，在減速段達到最大，最大阻力約為初始阻力（t=0 時）的1.1 倍。圖12 為工況2 模型運動過程中所受氣動阻力隨時間的變化曲線。從圖中可以看出：在加速段接近勻速段時，模型所受的氣動阻力開始逐漸減小，在減速段，氣動阻力達到最大，氣動阻力峰值約為加速段所受氣動阻力均值的2 倍。

圖11 工況1 模型運動過程氣動阻力隨時間變化曲線Fig.11 Aerodynamic drag vs.time curve during model moving for condition 1

圖12 工況2 模型運動過程氣動阻力隨時間變化曲線Fig.12 Aerodynamic drag vs.time curve during model moving for condition 2

2.3 管道內(nèi)波系傳播結(jié)果

圖13 和14 為模型在直線密閉管道內(nèi)加速、勻速和減速運動過程中周圍非定常流場演化的三維速度和Schlieren 數(shù)云圖。從圖中可以看出：波陣面的前移速度明顯快于模型的運動速度，在模型未到達管道末端時，波陣面已先于模型到達管道末端壁面并產(chǎn)生反射。

圖13 模型運動不同時刻速度云圖Fig.13 Velocity contours for model moving at different times

圖14 模型運動不同時刻Schlieren 數(shù)分布云圖Fig.14 Schlieren contours for model moving at different times

2.4 管道末端壓力分布

為開展消波措施設(shè)計，通過仿真計算獲得管道末端的壓力波變化情況。在管道模型的末端截面不同位置選取若干點（圖15），并輸出這些點的流場參數(shù)時間歷程曲線，如圖16 所示。

圖15 管道末端截面選取單元位置示意圖Fig.15 Layout for choosing elements’ positions at end wall

圖16 工況1 管道末端截面選取單元的壓力變化曲線Fig.16 Pressure vs.time of choosing elements for condition 1

表2 為本文氣動結(jié)構(gòu)耦合仿真方法（FSI）與傳統(tǒng)CFD 方法得到的風(fēng)洞沿程壓力和氣動阻力對比。從表2 可以看出：2 種方法的結(jié)果差異基本在5%以內(nèi)[21]，說明了本文方法的可行性；而整體上CFD 仿真結(jié)果小于FSI 仿真結(jié)果，可能是CFD 方法未考慮氣動與結(jié)構(gòu)耦合效應(yīng)的影響。

表2 FSI 與CFD 仿真結(jié)果[21]對比Table 2 Comparison results of FSI and CFD[21] simulation

3 管道消波措施分析

為確定多孔介質(zhì)參數(shù)，通過數(shù)值仿真進行分析。仿真模型如圖17 所示，給定二維模型，劃分為高壓區(qū)、低壓區(qū)和多孔介質(zhì)區(qū)，長度取1 m；高、低壓區(qū)壓差設(shè)置為5 000 Pa，通過高、低壓區(qū)壓差形成壓縮波，壓縮波強度為2 500 Pa，向右傳播；在低壓區(qū)中間平面處設(shè)置壓力監(jiān)測點，監(jiān)測面平均壓力隨時間的變化規(guī)律。

圖17 多孔介質(zhì)參數(shù)數(shù)值計算模型Fig.17 Numerical calculation model of parameters for porous medium

根據(jù)一維等熵理論，壓縮波在傳播過程中遇到壓力間斷面會產(chǎn)生反射壓縮波。本文計算了不同黏性系數(shù)（即不同流阻率）多孔介質(zhì)條件下反射壓縮波的強度。為了便于對比，首先計算了多孔介質(zhì)黏性系數(shù)為0（即無多孔介質(zhì)）時的壓縮波反射工況，結(jié)果如圖18 所示，壓縮波與右端壁面接觸后，向左方向反射時，無強度和能量損失。

圖18 無多孔介質(zhì)時壓縮波的傳播特性Fig.18 Compression wave propagation characteristics without porous medium

常規(guī)吸聲棉的黏性系數(shù)為5.6 × 108m-2，給定多孔介質(zhì)黏性系數(shù)，設(shè)置孔隙率為1，計算結(jié)果如圖19所示：壓縮波在與多孔介質(zhì)接觸后，形成強度約1 800 Pa 的反射壓縮波，說明該多孔材料黏性系數(shù)過大，造成了壓縮波在接觸面的強反射；后續(xù)壓力緩慢上升到4 890 Pa，說明多孔介質(zhì)對壓縮波有弛豫作用，且對壓縮波有110 Pa 對應(yīng)能量的耗散。

圖19 多孔介質(zhì)黏性系數(shù)為5.6 × 108 m-2 時壓縮波的傳播特性Fig.19 Compression wave propagation characteristics at viscous factor 5.6 × 108 m-2

降低多孔介質(zhì)黏性系數(shù)至1 × 107m-2時，計算結(jié)果如圖20 所示。從圖中可以看出：壓縮波與多孔介質(zhì)接觸后產(chǎn)生的反射波壓力緩慢增大了300 Pa，壓縮波與右端壁面接觸后的反射壓縮波通過多孔介質(zhì)的強度約為1 400 Pa，且后續(xù)約有300 Pa 的壓力緩升。此時壁面反射壓縮波強度過高，且最后壓力升高到5 000 Pa，說明此時多孔材料黏性較低，對能量耗散也較小。

圖20 黏性系數(shù)1 × 107 m-2 時壓縮波傳播特性Fig.20 Compression wave propagation characteristics at viscous factor 1 × 107 m-2

圖21 給出了無多孔介質(zhì)、有高黏性多孔介質(zhì)、有低黏性多孔介質(zhì)時監(jiān)測點的時域壓力變化曲線?？梢钥闯觯憾嗫捉橘|(zhì)對壓縮波具有顯著的弛豫作用。由于壓縮波強度與壓力梯度成線性關(guān)系，所以多孔介質(zhì)可以顯著降低壓縮波強度。多孔介質(zhì)黏性越高，壓縮波與其接觸后的反射強度越高，能量耗散也越高。因此，要達到較好的降壓效果，應(yīng)選擇黏性系數(shù)合理的多孔介質(zhì)，使得壓縮波盡可能地通過多孔介質(zhì)，同時又有不錯的弛豫作用。

圖21 無多孔介質(zhì)、有高/低黏性多孔介質(zhì)時監(jiān)測點的時域壓力變化曲線Fig.21 Time-domain pressure change curve of monitoring point without porous medium and with high/low viscosity porous medium

經(jīng)過計算分析，確定了一個較優(yōu)的黏性系數(shù)：2 × 107m-2。如圖22 所示，此時入射壓縮波和反射壓縮波強度分配合理，且都有不錯的弛豫效果。

圖22 多孔介質(zhì)黏性系數(shù)為2 × 107 m-2 時的監(jiān)測點壓力變化曲線Fig.22 Variation curve of monitoring point ’s pressure at viscous factor 2 × 107 m-2

綜上所述，尾端消波方案考慮填充黏性系數(shù)為2 × 107m-2、孔隙率≥0.95 的多孔介質(zhì)。有關(guān)此類問題的研究較少，沒有可以借鑒的方向，為達到最優(yōu)設(shè)計目標(biāo)，可通過仿真和實驗方法開展關(guān)鍵技術(shù)研究解決以下問題：

1） 此種多孔介質(zhì)非常規(guī)棉類多孔介質(zhì)，須專門設(shè)計制造。

2） 考慮到壓縮波的強度，對多孔介質(zhì)的抗沖擊性、抗疲勞性也有相應(yīng)的要求。

3） 目前確定的黏性系數(shù)僅針對強度為 2 500 Pa的壓縮波較優(yōu)，對于其他強度的壓縮波是否合適尚有待驗證，還可以繼續(xù)優(yōu)化。

4） 尾端消波方案還可以嘗試通過添加空氣層、使用阻抗?jié)u變的多孔介質(zhì)來優(yōu)化。

5） 多孔介質(zhì)的固定結(jié)構(gòu)須專門設(shè)計，應(yīng)既能固定多孔介質(zhì)，又不造成明顯的壓縮波反射。

4 結(jié) 論

本文基于CE/SE 方法對動模型管道內(nèi)模型高速運動過程進行了氣動結(jié)構(gòu)耦合仿真，建立的每個工況均計算收斂，且可以獲得模型運動過程中完整的流場結(jié)果。通過仿真結(jié)果評估了管道內(nèi)壓力波的傳播規(guī)律和模型運動過程所受氣動力的變化規(guī)律。開展了多孔介質(zhì)參數(shù)的仿真設(shè)計，計算結(jié)果可為管道消波的設(shè)計方案比選提供參考。