鄧 韜，楊志剛，2，賈 青

（1.同濟(jì)大學(xué)上海地面交通工具風(fēng)洞中心，上海201804；2.北京民用飛機(jī)技術(shù)研究中心，北京102211）

氣流流經(jīng)汽車，在汽車尾部形成復(fù)雜的非定常流動結(jié)構(gòu)。汽車尾跡的非定常特性對汽車氣動特性有重要影響［1］，因此有必要對復(fù)雜的流動結(jié)構(gòu)進(jìn)行研究。對于復(fù)雜的流動問題，本征正交分解（proper orthogonal decomposition，POD）可從能量角度對流場進(jìn)行分解，并根據(jù)各階模態(tài)和其模態(tài)系數(shù)的信息識別出大尺度的擬序結(jié)構(gòu)［2］。提取汽車尾跡中包含的大尺度擬序結(jié)構(gòu)，分析其特性，可以以此為理論依據(jù)對流場有針對性地施加控制，改善汽車氣動特性。

近年來，國內(nèi)外研究人員將POD方法運(yùn)用到類車體非定常尾跡的分析中。?sth等［3］基于大渦模擬（large eddy simulation，LES）的計(jì)算結(jié)果，采用POD研究了0°Ahmed body展向兩側(cè)渦脫落的機(jī)制。文獻(xiàn)［4-6］通過實(shí)驗(yàn)手段，用粒子圖像測速法對方背類車體尾跡速度場進(jìn)行測量，通過對速度快照進(jìn)行POD降階分析，揭示方背類車體尾跡雙穩(wěn)態(tài)現(xiàn)象的演化規(guī)律。國內(nèi)王國俊等［7］在LES數(shù)值模擬結(jié)果基礎(chǔ)上，用POD對25°Ahmed body斜背分離流進(jìn)行了分解，分析了各階模態(tài)的能量占比及其物理意義。

回顧過去的研究發(fā)現(xiàn)，基于POD的汽車非定常尾跡研究的對象大都集中于Ahmed body［8］和Windsor body［9］等類車體，與真實(shí)汽車相比，其幾何外形過于簡單，研究成果不能完全轉(zhuǎn)換到汽車工業(yè)中。為此，本文以更接近真實(shí)汽車的標(biāo)模DrivAer快背車為研究對象，采用改進(jìn)的延遲分離渦模擬（improved delayed detached-eddy simulation，IDDES）進(jìn)行數(shù)值仿真，與實(shí)驗(yàn)對比進(jìn)行驗(yàn)證；使用POD研究車體尾跡回流區(qū)的非定常特性，進(jìn)行模態(tài)分析，基于此進(jìn)行模態(tài)重構(gòu)的相位平均，提取出大尺度的擬序結(jié)構(gòu)并分析其特點(diǎn)，并對原流場進(jìn)行重構(gòu)。

1 研究對象

Ahmed body和MIRA等類車體幾何外形相對簡單，難以復(fù)現(xiàn)真實(shí)汽車的復(fù)雜流動結(jié)構(gòu)，樣車雖然能夠真實(shí)地展現(xiàn)流動特征，但是將其作為空氣動力學(xué)研究對象時(shí)，又受到較短的使用壽命和有限的獲取渠道的制約。為了縮短過度簡化的類車體和真實(shí)汽車之間的差距，德國慕尼黑工業(yè)大學(xué)空氣動力學(xué)及流體力學(xué)研究所基于BMW 3系和Audi A4這兩款中型車的幾何模型提出一種全新的汽車模型——DrivAer［10］。

DrivAer模型具有3種尾部造型風(fēng)格（兩廂掀背、階背、快背）和2種底盤結(jié)構(gòu)（復(fù)雜底盤、光滑底盤），如圖1所示。本文以縮比為1：2.5的底盤光滑的快背DrivAer模型為研究對象，模型長、寬、高分別為L=1.845m，W=0.812m和H=0.567m。

圖1 DrivAer及1：2.5模型尺寸Fig.1 DrivAer model and dimensions of 1:2.5 model

2 數(shù)值模擬

2.1 計(jì)算域及邊界條件

如圖2所示，計(jì)算域長、寬、高分別為11倍車長、11倍車寬、8倍車高。入口距離車頭4倍車長，入口速度U∞=40m·s-1，對應(yīng)雷諾數(shù)Re=4.87×106（基于車長L），湍流強(qiáng)度為1%；出口為零壓力出口，距離車尾6倍車長。DrivAer模型表面和地面為無滑移壁面條件，其余邊界為對稱邊界。邊界條件的設(shè)置與文獻(xiàn)［10］中無地面模擬的實(shí)驗(yàn)相對應(yīng)。

圖2 計(jì)算域Fig.2 Computational domain

2.2 網(wǎng)格策略

使用Hypermesh軟件劃分面網(wǎng)格，體網(wǎng)格通過Star-CCM+軟件生成。采用六面體非結(jié)構(gòu)化Trimmer網(wǎng)格對計(jì)算域進(jìn)行離散，并在車身壁面和地面生成邊界層網(wǎng)格，保證車身表面第一層網(wǎng)格法向高度對應(yīng)的量綱一化值y+小于1。

為了更好地模擬汽車周圍的流場并提高計(jì)算效率，對后視鏡、車輪和車體尾跡處的網(wǎng)格進(jìn)行加密；同時(shí)，為研究網(wǎng)格數(shù)量對數(shù)值計(jì)算結(jié)果的影響，通過改變車體尾跡加密區(qū)的網(wǎng)格尺寸，形成了2套網(wǎng)格方案，方案1和方案2體網(wǎng)格總數(shù)分別為2 245萬和4 416萬，圖3為計(jì)算域縱向?qū)ΨQ面網(wǎng)格。

圖3 流場縱向?qū)ΨQ面網(wǎng)格Fig.3 Mesh of flow field on longitudinal symmetric plane

2.3 數(shù)值方法

使用商業(yè)軟件Star-CCM+進(jìn)行數(shù)值計(jì)算，數(shù)值方法采用IDDES（基于SST模型）［11］，該方法是一種混合的RANS-LES模型，同時(shí)具備延遲分離渦模擬（delayed detached-eddy simulation，DDES）和壁面大渦模擬 （wall-modelled largeeddy simulation，WMLES）優(yōu)點(diǎn)。其中WMLES模型減小雷諾數(shù)影響，實(shí)現(xiàn)了高雷諾數(shù)下壁面邊界層的大渦模擬［12-13］。

在IDDES數(shù)值方法中，定義了一種新的亞格子尺度，如式（1）所示，它不同于只涉及網(wǎng)格間距的大渦模擬（large eddy simulation，LES）和分離渦模（detached-eddy simulation，DES）的亞格子尺度，它還包括了壁面距離的影響。

其中dw為壁面距離；hwn為壁面法向網(wǎng)格步長；Cw為與亞格子模型無關(guān)的經(jīng)驗(yàn)常數(shù)；hmax為當(dāng)?shù)刈畲缶W(wǎng)格間距，定義如下：

其中hi（i=x，y，z）為當(dāng)?shù)?個(gè)方向的網(wǎng)格尺寸。IDDES完整的方程可見文獻(xiàn)［13］。

空間離散采用基于有限體積法的混合格式，LES區(qū)域采用中心差分有界格式，RANS區(qū)域采用二階迎風(fēng)格式，并結(jié)合二階隱式方程。采用SIMPLE算法求解壓力-速度耦合，使用定常RANS結(jié)果初始化流場。

方案1和方案2的量綱一化時(shí)間步長ΔtU∞/L（其中Δt為時(shí)間步長）分別為4.3×10-3和2.2×10-3，保證計(jì)算域內(nèi)大部分區(qū)域的庫朗數(shù)C小于1（C=UabsΔt/Δδ，其中Uabs為穿過網(wǎng)格的速度大小，Δδ表示沿速度方向的網(wǎng)格尺寸）。

2.4 仿真結(jié)果驗(yàn)證

用量綱一化的時(shí)間T=L/U∞表征氣流流過一倍車長的時(shí)間，設(shè)定總計(jì)算時(shí)長為25T，取后20T的計(jì)算結(jié)果進(jìn)行平均，并與文獻(xiàn)［14］的實(shí)驗(yàn)結(jié)果進(jìn)行對比，實(shí)驗(yàn)與仿真的邊界條件一致。氣動阻力系數(shù)CD、氣動升力系數(shù)CL定義如下：

其中Fx為阻力，F(xiàn)z為升力，ρ為空氣密度，Ax為迎風(fēng)面積。壓力系數(shù)Cp定義如下：

其中p代表表面壓力，p0為參考壓力。

表1為不同網(wǎng)格方案下氣動力系數(shù)CD和CL的仿真結(jié)果，與實(shí)驗(yàn)結(jié)果對比發(fā)現(xiàn)，方案2的CD仿真與實(shí)驗(yàn)結(jié)果相差—0.003，誤差為1.15%，優(yōu)于方案1的0.004（1.53%），主要體現(xiàn)在車身背部壓力的區(qū)別上，如圖4所示，方案2車身背部高壓區(qū)面積大于方案1。對于CL仿真與實(shí)驗(yàn)值的對比，由圖5車身上下壓力系數(shù)分布可看出，與方案1相比，方案2上表面尾部壓力系數(shù)較高，下表面前軸區(qū)域壓力系數(shù)較小，總的效果是CL較小，更接近實(shí)驗(yàn)值。

從圖6可見，仿真值和實(shí)驗(yàn)值在流向上的變化趨勢一致，量值上相差較小。對比2種方案的仿真結(jié)果，車身中截面上側(cè)差異很?。卉嚿碇薪孛嫦聜?cè)X/L=－0.8附近以及車身Z=0.32H平面X/L=－0.9～－0.7處，方案2仿真結(jié)果與實(shí)驗(yàn)值更接近。

綜合對比仿真與實(shí)驗(yàn)結(jié)果，加密車體尾跡區(qū)的體網(wǎng)格可提高計(jì)算精度，因此下文的POD分析基于方案2的非定常仿真結(jié)果。

表1 不同網(wǎng)格方案氣動阻力系數(shù)對比Tab.1 Drag coefficient for different mesh strategics

圖4 車身背部壓力系數(shù)分布Fig.4 Mean of pressure coefficient distribution for the back of DrivAer Model

圖5 車身壓力系數(shù)分布Fig.5 Mean of pressure coefficient distribution of DrivAer Model from vertical direction

3 POD分析

3.1 POD原理

POD于1967年由Lumley［15］首先應(yīng)用到湍流相干結(jié)構(gòu)的研究中，基本思想是把時(shí)間和空間上連續(xù)的物理量構(gòu)成的場分解成只與時(shí)間和空間相關(guān)的函數(shù)展開序列，且它們在均方意義上最優(yōu)，這樣，只需要較少的項(xiàng)數(shù)就可以比較準(zhǔn)確地描述該物理過程，節(jié)省計(jì)算成本。

POD實(shí)際應(yīng)用中，往往出現(xiàn)空間監(jiān)測點(diǎn)數(shù)量遠(yuǎn)多于時(shí)刻數(shù)的問題，為此，Sirovich［2］于1987年提出了snapshots POD方法，每一時(shí)刻的瞬時(shí)流場稱為一張快照，用于定格流場變化過程中的瞬時(shí)狀態(tài)。給定 速 度 場q（ξ，t）（ξ=ξ1，ξ2，…，ξn；t=t1，t2，…，tm），其中n為空間監(jiān)測點(diǎn)數(shù)，m為時(shí)刻數(shù)。先計(jì)算平均速度場v（ξ），再用每張快照減去平均場得到對應(yīng)速度脈動分量u（t），并以列向量表示：

根據(jù)POD的基本思想，目的是找出最優(yōu)的基向量來表示u（t），該組基向量是完備且正交的。實(shí)際上，這個(gè)問題可以轉(zhuǎn)化為求特征向量φj和特征值λj的問題：

其中，j=1，2，…，m，C為u（t）的協(xié)方差矩陣。

式中，U為m張快照對應(yīng)的u（t）組成的矩陣：

則POD模態(tài)ψj可由特征向量φj映射得到：

由上述方法求得的特征值λj之和代表了系統(tǒng)的總能量，則每一階模態(tài)的能量占比為

因此，可以選取較低的維數(shù)r，使得

則原速度場可近似重構(gòu)為

式中，每階模態(tài)的模態(tài)系數(shù)ai（t）可由以下內(nèi)積求得：

圖6 車身中截面和Z=0.32 H平面壓力系數(shù)Fig.6 Pressure coefficient on longitudinal symmetric plane and Z=0.32 H plane of DrivAer Model

以方案2的非定常仿真結(jié)果為基礎(chǔ)，運(yùn)用snapshots POD方法對DrivAer快背車的非定常尾跡進(jìn)行分析。選取計(jì)算域中截面Y=0上汽車的尾跡回流區(qū)為本文POD分析的區(qū)域，如圖7中矩形方框所示，范圍為X/h=－0.28～ 1.52，Z/h=0～1.81，其中h=0.288m，為車尾高度。POD快照的采樣頻率為10 kHz，空間監(jiān)測點(diǎn)數(shù)為22 771個(gè)。

圖7 中截面Y=0瞬態(tài)流線圖及POD分析區(qū)域Fig.7 Streamlines on longitudinal symmetric plane and the region for POD analysis.

3.2 POD結(jié)果分析

3.2.1 收斂性判斷

為了使POD的結(jié)果收斂，同時(shí)排除快照數(shù)量對結(jié)果的影響，分別基于5 000、6 000和7 000張快照進(jìn)行POD分析，前20階模態(tài)能量占比如圖8所示。由圖可知，6 000張與7 000張的結(jié)果相差不大，對POD分析結(jié)果的影響可忽略不計(jì)。取6 000張快照進(jìn)行POD分析，滿足收斂性要求。

3.2.2 POD模態(tài)分析

POD前100階的能量占比如圖9所示，低階模態(tài)代表流場中能量較高的大尺度流動結(jié)構(gòu)，在流場發(fā)展中占支配地位，為整個(gè)流動的載體。前4階模態(tài)能量占比分別為6.78%、5.61%、2.98%和2.62%，與其他模態(tài)有較好的區(qū)分度，其中前2階能量占比較為突出；5階以上模態(tài)能量過度比較平滑；20階以上模態(tài)能量均小于1%，且逐漸趨近于零。圖10為渦量在Y方向的分量Ωy的前4階模態(tài)及其模態(tài)系數(shù)隨時(shí)間的變化情況。

圖8 不同數(shù)量快照下前20階模態(tài)能量占比Fig.8 Energy contribution of POD modes 1～20 with different number of snapshots

圖9 前100階模態(tài)能量占比Fig.9 Energy contribution of POD modes 1～100

模態(tài)表征渦量的空間分布特點(diǎn)。從圖10可看出，Ωy模態(tài)1和模態(tài)2的渦量分布特點(diǎn)相似，渦量均主要分布在車體尾跡下剪切層區(qū)域，且渦量集中區(qū)沿著剪切層方向正負(fù)交替分布，而上剪切層區(qū)域渦量較少；模態(tài)3和模態(tài)4的渦量在尾跡上下剪切層的分布比較均勻。

圖10 渦量Ωy前4階模態(tài)及其模態(tài)系數(shù)Fig.10 The first four POD modes of vorticity Ωyand mode coefficients

模態(tài)系數(shù)表征模態(tài)的幅值隨時(shí)間的變化規(guī)律。如圖10所示，模態(tài)的階數(shù)越低，模態(tài)系數(shù)的振幅越大，對應(yīng)的能量占比也越大。模態(tài)系數(shù)隨時(shí)間的變化均圍繞著零點(diǎn)上下波動，第1階模態(tài)和第2階模態(tài)的周期性較強(qiáng)，而3階和4階模態(tài)系數(shù)的周期性較弱。對前4階模態(tài)系數(shù)進(jìn)行功率譜密度（power spectral density，PSD）分析，其中頻率f用車尾高度h和來流速度U∞進(jìn)行量綱一化，以斯特勞哈爾數(shù)（Strouhal number）St=f h/U∞表示，如圖11所示。

圖11 前4階POD模態(tài)系數(shù)功率譜密度分析Fig.11 Frequency spectrum analysis of the first four POD mode coefficients

圖12 模態(tài)系數(shù)a1與a2相位分析Fig.12 Phase analysis of mode coefficients a1and a2

從圖11可以發(fā)現(xiàn)，模態(tài)階數(shù)越低，PSD峰值越高，更直觀地表征了與能量的關(guān)系；模態(tài)1和模態(tài)2的模態(tài)系數(shù)能量均集中在低頻區(qū)域，PSD峰值均出現(xiàn)在St=0.216（f=30Hz），表明模態(tài)1和模態(tài)2對應(yīng)同一種低頻的擬序流動結(jié)構(gòu)，主頻為0.216。進(jìn)一步對前2階模態(tài)系數(shù)作相位分析，如圖12所示，相位分析表明，St=0.216下，二者相位差為1.43rad（0.455 π）。

3.2.3 前2階模態(tài)重構(gòu)及相位平均

為更好地解前2階模態(tài)對應(yīng)的擬序流動結(jié)構(gòu)，研究其在一個(gè)周期內(nèi)的演化規(guī)律，對前2階模態(tài)進(jìn)行重構(gòu)并進(jìn)行相位平均。相位平均的方法參照文獻(xiàn)［16］，用式（7）得到的特征值λ1和λ2分別對模態(tài)系數(shù)a1和a2進(jìn)行量綱一化，過濾其能量信息。

每個(gè)時(shí)刻的重構(gòu)用a*1、a*2表示，做出其散點(diǎn)圖如圖13所示。從圖13a可見，a*1、a*2隨時(shí)刻數(shù)m呈螺旋式變化，圖13b顯示其在二維平面（a*1，a*2）上的投影分布在半徑為2圓心為（0，0）的圓內(nèi)，則a*1和a*2可根據(jù)相位角θ進(jìn)行儲存，保存其相位信息，如圖14所示；a*1和a*2在相位上均呈周期性分布，可用正弦函數(shù)對進(jìn)行曲線擬合，其中相位角θ定義為：

其中j=1，2，…，6 000，表示時(shí)刻數(shù)，每個(gè)時(shí)刻對應(yīng)一個(gè)θ值。

圖13 量綱一化模態(tài)系數(shù)a*1、a*2散點(diǎn)圖Fig.13 Scatter plot of normalized mode coefficients a*1anda*2

將區(qū)間［0，2π］以長度π/4進(jìn)行8等分，可得到8個(gè)相位區(qū)間：［0，π/4］，［π/4，π/2］，…，［7π/4，2π］。根據(jù)θ值將m個(gè)重構(gòu)結(jié)果儲存到8個(gè)相位區(qū)間里，并對每個(gè)區(qū)間取平均值，得到渦量Ωy重構(gòu)的8個(gè)相位平均，如圖15所示。圖中可看出，在一個(gè)周期內(nèi)，上剪切層渦量強(qiáng)度較弱，變化較小，而下剪切層渦量強(qiáng)度較大，且變化較為劇烈，意味著較高的能量耗散，說明前2階模態(tài)對應(yīng)的低頻擬序流動結(jié)構(gòu)由車體尾跡下剪切層周期性運(yùn)動主導(dǎo)。

圖14 a*1、a*2相位信息Fig.14 Phase information ofa*1anda*2

3.2.4 基于POD結(jié)果的流場重構(gòu)

根據(jù)平均流場，POD模態(tài)及對應(yīng)的模態(tài)系數(shù)重構(gòu)流場，圖16為T=11.64時(shí)刻原始流場與不同階數(shù)POD模態(tài)重構(gòu)的結(jié)果的渦量云圖對比。從圖中可以看出，隨著使用的POD模態(tài)階數(shù)的不斷提高，對原始流場的細(xì)節(jié)捕抓得更清楚，前273階模態(tài)（能量占比為93%）與原流場差別很小，基本可以抓住流場的主要特征，且增加階數(shù)到前323階（能量占比為95%）重構(gòu)結(jié)果差別很小。由此可見，前273階模態(tài)就可以很好地重構(gòu)流場，階數(shù)由6 000階減少到了273階，極大地降低了自由度。

4 結(jié)論

圖15 前2階模態(tài)重構(gòu)的渦量Ωy的相位平均Fig.15 Phase-averaged of reconstructed vorticity Ωy field based on the first two modes

（1）采用IDDES對DrivAer快背車進(jìn)行數(shù)值仿真，對氣動力系數(shù)與實(shí)驗(yàn)的對比發(fā)現(xiàn)，加密車體尾跡區(qū)的體網(wǎng)格可提高計(jì)算精度，車身中截面和車身Z=0.32H平面壓力系數(shù)與實(shí)驗(yàn)結(jié)果變化趨勢相同，量值相差較小，表明本文數(shù)值方法是正確可行。

（2）對車體尾跡中截面回流區(qū)進(jìn)行POD分析，前2階模態(tài)能量占比分別為6.78%和5.61%，合計(jì)占總能量的12.39%，模態(tài)系數(shù)的頻譜分析表明，前2階模態(tài)對應(yīng)同一種低頻的擬序流動結(jié)構(gòu)，主頻為0.216（30Hz）；相位分析表明前2階模態(tài)相位差為0.455π。通過重構(gòu)渦量場并進(jìn)行相位平均，可以直觀地觀察到前2階模態(tài)對應(yīng)的擬序流動結(jié)構(gòu)的演化過程，發(fā)現(xiàn)該流動結(jié)構(gòu)由車體尾跡下剪切層周期性運(yùn)動主導(dǎo)。后續(xù)研究可以以該擬序結(jié)構(gòu)的特征為參考，對流場有針對性地施加控制，探究改善汽車氣動特性的方法。

圖16 原始流場與不同階數(shù)POD重構(gòu)結(jié)果對比Fig.16 Comparison of the original flow field and reconstructed vorticity Ωyfield based on the different modes

（3）對于本文研究的汽車尾跡回流區(qū)，前273階POD模態(tài)（能量占比93%）可以很好地重構(gòu)原流場，極大地降低了自由度。