飛行器繞流流場脈動壓力環境預示方法探討

（中國工程物理研究院總體工程研究所，四川 綿陽 621999）

導彈武器裝備在貯存、運輸以及飛行的全壽命剖面中，將受到各種環境的影響。其中，在再入飛行過程中，飛行器機動飛行、大氣擾流、壁面凸起或凹槽結構、外形非光滑轉折等都會在飛行器繞流流場中產生誘導擾流，導致流場參數（壓力、密度、速度、溫度等）脈動，在飛行器表面形成復雜的脈動載荷環境。脈動壓力環境產生的原因可分為附著湍流邊界層、分離區和激波/邊界層的相互干擾區域等3 種[1]。脈動壓力作用于飛行器物面，誘導結構振動，可能導致飛行器結構破壞。同時，脈動壓力可能誘發強烈的噪聲環境，影響機載儀器設備的正常工作，對武器裝備的環境適應性、可靠性、安全性等性能帶來嚴峻考驗。

飛行器繞流流場脈動載荷是誘導飛行器結構振動的主要激勵源，是再入體結構設計中僅次于熱環境的重要動力學環境因素[2]，是飛行器結構響應分析和動力學環境研究的重要依據。開展脈動壓力環境預示研究，對于飛行器結構設計、結構動力學響應分析和環境適應性研究具有重要意義。文中對脈動載荷的形成機理、流動特征、測量方法、數值模擬方法進行了歸納總結，簡要分析了脈動載荷定量評估的難點，并就飛行器脈動壓力環境預示方法進行探討。

1 脈動壓力環境特點

到目前為止，一般認為包括脈動運動在內的湍流瞬時運動服從Navier-Stokes 方程（N-S 方程），而N-S 方程是具有非線性、時間相關特征的三維偏微分方程。

飛行器在大氣中飛行時，由于空氣黏性作用和飛行器表面黏性流體的無滑移邊界條件，飛行器表面附近區域流動形成很強的速度梯度，在飛行器結構表面產生足夠的渦量，形成較強的黏性應力。渦量在黏性作用下，擴散到流體內部，并向下游輸運和耗散[3]。大量研究表明，飛行器近壁區域流動在運動過程中，旋渦不斷經歷形成、輸運、破碎、合并和耗散等變化過程，渦的拉伸、破裂以及渦之間的相互干擾運動使得流場中速度、壓力、溫度和密度等流動參數不斷變化，即產生脈動現象。流動呈現出非線性、隨機、擬序結構等特征，如圖1 和圖2 所示。近壁面湍流擬序結構研究表明，湍流邊界層的“猝發”（burst）、“下掃”現象直接影響壁面脈動壓力環境[4]。飛行器近壁面區域流場呈現湍流狀態，該區域流動參數隨機變化誘導形成了脈動壓力環境。

圖1 湍流渦流動結構Fig.1 Structure of turbulence vortex

圖2 湍流邊界層(來自Van Dyke,1982)Fig.2 Turbulence boundary layer(From Van Dyke,1982)

脈動壓力環境的特征主要有如下幾點。

1）空間尺度范圍大。湍流脈動的空間尺度包含了最大湍渦與最小湍渦尺度之間的所有空間尺度，其中最大湍渦尺度與流動特征密切相關，大渦主要受慣性影響，尺寸與流場大小相當，是引起低頻、高幅脈動壓力的主要原因。最小湍渦尺度由黏性決定，也就是取決于分子作用力，即取決于由Kolmogorov 定義的黏性耗散尺度（Kolmogorov 尺度）η∝(υ3/ε)1/4[4]，其中υ和ε分別為運動黏性系數和湍能耗散率。該尺寸可能只是流場尺度的1/1000 量級，是引起高頻、低幅脈動壓力的主要原因。由于湍流中包含了所有尺度的渦，脈動壓力環境具有寬帶特征。

2）時間尺度范圍大。湍流脈動具有連續的頻譜，脈動時間尺度涵蓋最小湍渦時間尺度至最大湍渦時間尺度之間的所有時間尺度，最大湍渦的時間尺度為L/u′，最小湍渦的時間尺度為η/u′[4]，其中L和u′分別為最大渦尺度和大渦的特征速度，因此脈動時間尺度范圍極大。

3）脈動量值小。與時均流場參數相比，飛行器繞流流場脈動量值非常小，約比時均參數小4～5 個數量級。實驗研究很容易被背景噪聲所湮沒，同時測試儀器誘導的擾動也可能影響真實的脈動壓力環境；而數值模擬必須保持低色散、低耗散的高精度計算格式，降低數值計算對脈動壓力環境的影響。

4）三維性。湍流不穩定性源自于N-S 方程中非線性慣性項和黏性項之間的有旋、時間相關和全三維的相互作用。由于渦拉伸作用不可能存在于二維空間，相互作用中的有旋性、時間相關和三維性又是彼此關聯的，二維近似處理不能合理解釋湍流現象，因此湍流脈動呈現為強三維性。

5）隨機性和輸運傳遞特性。湍流是一種極不規則的流動現象，其不規則性不僅表現在速度、壓強等流動物理量在時空中的不規則分布，還表現在不可重復性，脈動載荷具有很強的隨機性。由于湍流中湍渦能量的傳遞和輸運特性，而湍流脈動又被認為是流體波動，飛行器壁面脈動壓力環境易受流場區域內、各邊界處的擾動和限制條件的影響。

2 研究方法

理論分析、實驗研究和數值模擬是研究脈動壓力問題的三大手段。在研究早期，由于計算機能力的限制，各國的研究學者多以理論分析和實驗研究為主。隨著數值計算方法和計算機平臺運算能力的不斷發展，數值模擬分析在脈動壓力環境預示中發揮了重要的作用。在流體力學和空氣動力學發展過程中，學者們將3 種方法有機結合，互為補充，加深了對脈動壓力問題的認識和理解。

2.1 理論分析

理論分析具有適用范圍廣、影響因素清晰等優點，能夠用于指導后續的實驗研究和驗證性的數值模擬方法。對于特定的湍流脈動來說，同樣可參照N-S方程進行定量分析。參照雷諾平均方法的思路，將可壓縮瞬時流動的質量、動量方程中各瞬時物理量分解為時均量與相應脈動量之和，推導出脈動壓力滿足的偏微分方程，如式（1）所示。

針對特定的湍流邊界層流動，可對上述脈動壓力方程等式右邊各項進行量階分析，消除某些對脈動壓力影響較小的項，通過積分即可獲得脈動壓力的量階估算式，結合一些實驗數據，即可擬合建立脈動壓力數值計算模型。

文獻[5-6]以二維湍流邊界層為例，通過分析式（1）右邊各項的量階，消除第二項、第五項小量綱項，聯合求解邊界層外緣處的平均流場和湍流脈動動能K的控制方程，最終獲得了該湍流邊界層內脈動壓力的數值計算模型，如式（2）所示。

式（2）中的常數C可由實驗數據擬合給出。國外就錐形、圓柱形等氣動外形的脈動壓力環境開展實驗測量，對脈動壓力方程式（1）進行簡化，獲得了適用于多種氣動外形、流動馬赫數等參數的脈動壓力計算模型，有效地推動了脈動壓力環境預示及動力學環境適應性研究。

復雜湍流流動中，脈動生成、輸運、黏性和壓縮性都會影響脈動壓力環境，定量分析脈動壓力各項量階大小十分困難，甚至變得不可能。此外，某一脈動壓力數值計算模型只適用于某些特定湍流現象，飛行器氣動外形及飛行工況一旦發生變化，計算模型就不再適用。

2.2 實驗研究

實驗研究所得到的測量數據可信度相對較高，在脈動壓力環境研究中發揮了重要作用。20 世紀60 年代以來，國外科學家就飛行器脈動壓力問題進行了大量的實驗研究，在實驗測試方法、技術等方面進行了相應探討[7-11]，如圖3 所示。結合理論推導和實驗數據，獲得了適用于錐形、圓柱形等結構外形的脈動壓力預示模型[12-13]，擬合了一系列針對特定外形的經驗公式[14-15]，推動了飛行器動力學環境激勵預示技術的發展。

文獻[15]給出了附體湍流邊界層流動的脈動壓力預示公式，如式（3）所示。式（3）根據飛行器物面壓力系數Cp(x)=[Pe(x)－P∞]/q∞預示均方根脈動壓力系數，

圖3 風洞實驗研究[10]Fig.3 Wind-tunnel test[10]:a) exploded view of experimental model of pressure-fluctuation cone;b) magnified cone head

式中：下標∞、e、W 分別為無窮遠、附面層邊緣和表面參數；x為沿壁面子午線距駐點的距離；γ為比熱比；Pr 為普朗特常數。文獻[15]同時也給出了功率譜密度與交叉功率譜密度的預示方法，并對某一球頭雙錐帶控制翼的機動再入飛行器擾流流場的脈動壓力環境進行預示，預示結果與實驗數據吻合較好。

由于脈動載荷的小量值、輸運傳遞、三維等特性，如果利用高動態壓力傳感器等類似的儀器進行實驗測量，測量儀器的布置將對原有湍流流場引入新的擾動，如何降低測量儀器對湍流流動的擾動是脈動壓力環境測量方法、測試技術中需要重點關注的問題。如果利用PIV 等流動顯示技術對脈動壓力環境進行測量，示蹤粒子的跟隨性、激光的頻率都會成為限制脈動壓力環境測量的因素。模型尺寸、流場擾動和測量精度等因素往往會成為實驗開展的限制因素，同時實驗研究還會遇到人力、物力、財力消耗大以及實驗周期長的問題。

2.3 計算流體動力學（Computational Fluid Dynamics，CFD）

隨著計算流體力學技術和計算機平臺的飛速發展，越來越多的湍流問題可通過數值計算進行求解，利用數值求解再現并解釋湍流流動現象。湍流的模擬是脈動壓力環境數值計算的關鍵，因此脈動壓力的數值求解方法與當前湍流數值模擬的方法相對應，主要包括：直接數值模擬（Direct Numerical Simulation，DNS）、湍流模式、大渦模擬（Large-Eddy Simulation，LES）、計算氣動聲學（Computational Aeroacoustics，CAA）以及基于上述方法建立的混合數值模擬方法。

2.3.1 直接數值模擬方法

直接數值模擬無需湍流模型化，通過直接求解三維非定常N-S 方程來模擬湍流瞬時流場。因此該方法的數值模擬可以精確獲取每一瞬時流場上的全部信息，提供很多在實驗上目前還無法測量的流場參量。

由于DNS 方法直接對N-S 方程進行數值求解，為了精確求解和刻畫湍流流動現象，離散網格必須能夠捕捉流場中所有尺度的流動結構，計算區域網格的尺寸在大到足以包含最大尺度渦的同時，還要小到足以刻畫最小尺度渦的特征，大小尺度的比值隨著Re的提高而迅速增大，整個計算域上的網格數量至少為Re9/4。另外，數值計算格式必須具有較高的精度，減少數值色散和耗散度，數值模擬過程中的時間長度要大于大渦的時間尺度，且計算的時間步長要小于小渦的時間尺度，總的計算量正比于Re3。在工程應用中的高雷諾數湍流流動中，湍流渦的時間、空間尺度范圍較大。為了模擬最小尺度渦的流動，需要劃分很精細的近壁面網格，消耗大量的計算時間和計算機容量，而這遠遠超過現有的計算機能力。

由于DNS 需要耗費龐大的計算資源和巨額的計算時間，在工程應用所關注的頻域范圍內（如隨機振動所關心的頻段為10～2000 Hz，噪聲關心的頻段為50～10 000 Hz），實際流動雷諾數為千萬量級，其計算工作量更加驚人。目前難以直接應用于求解飛行器真實尺度的工程問題，但這不妨礙DNS 成為湍流機理研究的有力工具。隨著計算機能力和求解方法的不斷改進，DNS 目前不僅能求解壓縮拐角、膨脹角等簡單幾何的激波邊界層干擾問題[17-18]，也能求解來流馬赫數為6，以頭半徑定義的來流雷諾數為10 000 的三維錐體轉捩問題[19]（如圖4 所示），求解問題的雷諾數也逐漸擴展到106量級[20]?？梢灶A見，在未來的脈動壓力環境研究中，DNS 將發揮越來越重要的作用。

圖4 采用DNS 方法求解三維錐體的轉捩問題[18]Fig.4 The transition solution of three-dimensional blunt cone using DNS [18]

2.3.2 湍流模式理論方法

由于N-S 方程的非線性，通過解析方法來精確刻畫三維全部湍流信息極其困難。在工程應用中，通常更關注湍流所引起平均流場的變化，即時均參量，因而工程數值仿真通常將非定常N-S 方程作平均處理。在三維N-S 方程計算模型中，雷諾平均法是目前使用最廣的湍流數值模擬方法。該方法將湍流流動中的任何變參量都分解為平均值和脈動值等2 部分。雷諾平均法的核心不是直接求解瞬時的N-S 方程，而是求解雷諾平均N-S 方程，即RANS 方程。

由于RANS 方程中包含了脈動量乘積的時均值未知數，方程個數少于未知數個數，方程無法封閉。為了使方程封閉，需依據湍流理論知識、實驗數據或直接數值模擬結果，根據經驗數據、物理類比對雷諾應力作出各種假設，建立基于經驗和半經驗的本構關系，即通過建立湍流模型封閉RANS 方程。目前常用的湍流模型有雷諾應力模型和渦粘模型2 類。

無論采用何種平均方法和湍流模型，脈動運動的時空變化細節在平均的流場結果中都會被一概抹平，數值仿真結果無法獲得全部的脈動運動信息，難以直接對脈動壓力環境進行預示。工程上多采用數值模擬與經驗公式相結合的方法[21]，通過RANS 方法求解時均流場，以流場信息作為輸入，利用工程經驗公式獲得脈動壓力環境特性的統計參數，如圖5 所示。

2.3.3 大渦模擬方法

由于DNS 方法對流場全尺度的渦均進行模擬刻畫，需要耗費大量的計算資源和時間?？紤]到大渦與平均流場之間的相互作用比較強烈，大渦的形態和強度與流場特征密切相關，各向異性較強；而小渦主要是大渦之間非線性作用下的產物，受平均流場或邊界特性的影響相對較小，可近似認為具有各向同性特征。因此LES 將比網格尺度大的湍流運動通過DNS求得，而小尺度渦對大尺度渦運動的影響則通過一定模型，通常是亞格子模型（多為Smagorinsky 模型，簡稱SGS 模型）來刻畫。Ferzinger 對LES 給出了一個較為精確的定義：LES is same as any simulation of a turbulent flow in which the large-scale motions are explicitly resolved while the small-scale motions are represented approximately by a model (in engineering nomenclature) or parameterization (in the geosciences)。LES 是介于DNS 與RANS 之間的一種湍流數值模擬方法，與DNS 相比，LES 的計算量顯著降低。

圖5 基于半經驗公式和時均流場的預示方法[21]Fig.5 The combination of time-averaged flow field and semi-empirical formula[21]:a) structure of time-averaged flow field;b) axial distribution of RMS fluctuating pressure coefficient;c) distribution of power spectral density

從飛行器動力學的環境適應性角度來說，關注的脈動壓力主要集中在10～10 000 Hz 這一頻段范圍內。湍流流動中脈動信息的頻率特性主要與大渦和小渦的時間尺度有關：其最高頻率由最小湍渦的時間尺度η/u′≈(ν/ε)1/2決定，當來流速度約5000 m/s、特征長度約1 m、大氣密度約0.5 kg/m3時，其最小湍渦對應的頻率約為1012Hz 量級，遠大于工程應用所關注的最高頻率。因此脈動壓力環境可利用大渦模擬進行數值求解。

LES 需選擇合適的計算格式，格式選擇不當可能會引起嚴重的數值耗散，甚至淹沒亞格子（sub-grid）的作用。目前LES 在理論研究中多采用耗散性較小的高階無耗散中心格式或緊致格式。此外，LES 的SGS 模型只對小尺度渦（各向同性的相似渦）進行?；幚?，該尺度以上的湍渦信息仍需通過數值求解瞬時N-S 方程獲得。DNS、LES 和RANS 等3 種數值求解方法的對比如圖6 所示。

由以上可知，盡管LES 的計算量比DNS 要小幾個數量級（NLES≈(0.4/Re1/4)NDNS）[16]，但也比RANS大出幾個數量級。由于再入飛行器實際流動的雷諾數通常為千萬量級，將LES 應用于工程脈動壓力環境計算，仍將耗費巨大資源和時間。工程中應用LES，必須有足夠的減小計算量的手段。為了彌補RANS方法捕捉分離流動的不足，同時提高LES 方法求解邊界層流動的效率，RANS/LES 混合方法應運而生，并得到迅速發展。該類方法頗具代表性的是由Spalart[22]提出的脫體渦模擬法（Detached Eddy Simulation，DES），DES 結合RANS 和LES 的優點，對S-A 湍流模型進行長度尺度修正，近壁面通過RANS 求解附面層流動，遠離物面的主流和分離區采用LES 方法進行求解，在降低計算量的同時，計算精度也令人滿意。國內外學者采用該方法對鈍體外形的亞聲速、跨聲速以及超聲速流場[23-26]進行了數值模擬，背風面計算結果與風洞實驗數據吻合較好，計算結果如圖7 所示。

圖6 3 種求解方法對比[16]Fig.6 Comparison of three solution[16]

圖7 通過DES 求解得到的瞬時和時均流場的馬赫數云圖[24]Fig.7 Plots of Mach number for an in stantaneous and time-averaged flow field obtained by DES solution[24]:a) instantaneous flow field;b) time-averaged flow field

2.4 計算氣動聲學方法

飛行器繞流流場脈動是流動噪聲的表現形式，即湍流流動同時形成湍流噪聲[27]。湍流噪聲隨湍流流動的產生而產生，并隨湍流流動的消失而消失，因而脈動壓力環境也可從氣動聲學角度進行研究。氣動聲學的研究在20 世紀50 年代由Lighthill 從流體力學基本方程得到了對流場的外圍聲場的波動方程[28-31]，建立聲比擬理論，并逐漸形成氣動聲學這一新學科。

目前，CAA 領域的計算方法主要包括兩類[31]：一類是直接模擬，即在一定計算域內運用經過特殊處理的CFD 方法直接計算近場的流場和聲場，包括上述介紹的DNS、LES 及RANS，然后用Kirchhoff 積分外推到遠場；另一類是基于Lighthill 氣動聲學理論的聲比擬方法，此法將流場與聲場分開計算，近場流場由LES、RANS 計算獲得，遠場由各種理論模型得到。目前主要的方法有Lighthill 聲比擬、渦聲理論、Hardin 和Pope 的分裂方法、有聲源項的線化歐拉方法、聲擾動方程等。采用CAA 方法求解得到的增升裝置瞬態攝動壓強流場如圖8 所示。

CAA 主要關注湍流噪聲（脈動壓力）的非定常流動機理、聲源確定、聲與流體的相互作用、聲波的傳遞等研究。利用CAA 數值方法求解飛行器近場聲場，獲得飛行器動力學環境（振動）相應的環境激勵源——脈動壓力環境，即CAA 可應用于飛行器繞流流場脈動壓力環境研究。如通過RANS 方法計算近場流場聲源，利用Lighthill 方程計算其脈動密度，進而獲得脈動壓力環境，聲學方程見式（4）。

圖8 CAA 方法求解的瞬態攝動壓強流場[32]Fig.8 Instantaneous perturbations pressure flow field using CAA[32

式中：ρ′=ρ-ρ0，。

CAA 采用的數值方法與CFD 方法存在緊密聯系，特別是它們在湍流脈動信息數值處理方面，同DNS 和LES 類似，CAA 在數值耗散、色散、高精度離散格式、多尺度等方面具有很高的要求。此外，CAA邊界條件（需要無反射邊界條件）的要求十分苛刻。利用CAA 數值求解飛行器繞流流場脈動壓力環境仍需耗費較大計算資源和計算時間，工程實現存在較大難度。

3 脈動壓力工程研究與應用可行途徑

脈動壓力環境是飛行器飛行任務剖面主要的動力學環境激勵源，為飛行器振動噪聲環境適應性研究提供數據輸入，直接影響飛行器動力學結構響應。由于再入飛行器飛行馬赫數高，繞流流場雷諾數高達108，如果利用DNS 數值求解飛行器繞流流場脈動壓力環境，其計算量約比雷諾數105的流動大9 個量級，將極大耗費計算資源和計算時間，短期內難以實現工程研究與應用，但應用前景誘人。

為了實現計算效率和計算精度的平衡，LES 方法和以DES 方法為代表的混合方法不斷發展，該方法已應用于求解雷諾數約為107的湍流流動。然而，由于LES 和DES 方法分別利用亞格子模型和RANS 方法求解湍流邊界層，它更適用于模擬遠場流動分離、尾跡渦等脫體渦分離流動。對于飛行器再入過程中湍流邊界層占主導的脈動壓力環境問題而言，難以實現準確的脈動壓力環境預示。

對于CAA 而言，多采用基于聲比擬的方法，將流場和聲場分別求解，通過流場結果構造聲源，再利用聲傳播方程獲得遠離物面流場的脈動壓力環境。也就是說CAA 方法在不額外增加較大計算量的前提下，能夠準確地獲得遠場的脈動壓力，但對于物面處的脈動壓力環境預示多取決于求解流場采用的方法精度，對于預示物面的脈動壓力提升有限。

根據以上分析可知，由于飛行器繞流流場雷諾數高，脈動壓力環境具有多尺度、量值相對較小等特征，湍流脈動數值模擬方法和技術尚處于湍流流動機理分析和研究階段，短期內無法解決飛行器繞流流場脈動壓力環境工程問題。分析現有的飛行器脈動壓力環境預示技術，可能存在如下幾條途徑。

1）風洞實驗測量預示方法。國內外就飛行器脈動壓力環境的風洞實驗測量技術和方法進行了大量研究，成功運用于脈動壓力環境預示。由于風洞實驗狀態有別于飛行器飛行狀態（雷諾數、密度等），風洞實驗預示的脈動壓力環境與飛行器真實飛行脈動壓力環境存在一定差異，實驗測量數據處理和修正尚需開展深入研究。

2）基于半經驗公式和時均流場的預示方法。針對典型飛行器的氣動外形和飛行工況，國內外已經研究、推導和構建了豐富的基于時均流場的脈動壓力半經驗公式預示方法，極大地推進了飛行器脈動壓力環境預示研究進展。由于這類方法在很大程度上依賴于半經驗公式對氣動外形和飛行參數的適應性，需針對特定飛行器和飛行工況的半經驗公式進行適應性研究和改進，提高預示精度和準確性。

3）基于DNS 的數值模擬方法。由于飛行器再入過程中速度高，雷諾數大，物體表面的渦尺度較小，只有精確捕捉小尺度的渦才能實現脈動壓力環境的準確預示。從數值模擬方法原理來看，目前只有DNS方法才能實現物面小尺度渦的準確刻畫。但由于其龐大的計算量，目前無法解決高雷諾數復雜三維外形的脈動壓力問題。隨著計算能力的提升，數值方法的改進，目前已經可以求解三維錐體轉捩問題，雷諾數也逐漸擴展到106量級，DNS 方法未來在脈動壓力環境問題研究中有著光明的前景。

4 結論

1）飛行器表面壓力脈動源于流動黏性、無滑移邊界條件和剪切流等流動因素，其本質屬于湍流流動問題，在宏觀上表現為氣動噪聲（聲波），一般認為其物理流動仍服從N-S 方程。

2）風洞實驗是目前獲取飛行器表面脈動壓力環境的有力手段，但風洞實驗工況與真實飛行環境仍有一定差異，需進行測試數據的處理與修正。

3）基于半經驗公式和時均流場的預示方法能夠快速評估獲得飛行器表面的脈動壓力環境，在工程問題中得到了廣泛應用，然而半經驗公式與飛行器氣動外形與飛行工況密切相關，難以獲得普適的關系式。

4）直接數值模擬方法能夠精確刻畫飛行器表面流場的渦系結構，在湍流機理研究方面能夠發揮重要作用，由于該方法龐大的計算量，目前難以解決高雷諾數問題，但隨著計算機能力的提升和數值方法的改進，直接數值模擬方法在未來的脈動壓力環境預示研究中將發揮重要作用。