賀巖松， 曹禮鵬， 張志飛，李 云，陳 釗

(1.重慶大學(xué)汽車工程學(xué)院，重慶 400030； 2.東風(fēng)柳州汽車有限公司，柳州 545005)

前言

汽車高速行駛時，氣動阻力嚴(yán)重影響汽車燃油經(jīng)濟性[1]。 在汽車開發(fā)階段，大幅優(yōu)化汽車外形可減小氣動阻力，而量產(chǎn)車型的氣動阻力優(yōu)化則受到原造型風(fēng)格的限制，優(yōu)化空間有限，面臨著巨大的挑戰(zhàn)。

通過優(yōu)化技術(shù)建立車輛氣動優(yōu)化的數(shù)學(xué)模型可以達(dá)到通過小幅結(jié)構(gòu)改動而減阻的目的。 在方法上，Jameson 等[2]結(jié)合控制理論首次提出了連續(xù)伴隨方法，實現(xiàn)了優(yōu)化過程中梯度的高效求解。 Nielson等[3]在此基礎(chǔ)上推導(dǎo)出了離散伴隨法，使得梯度的求解更加快速精確。 在應(yīng)用方面，劉峰博等[4]使用離散伴隨法實現(xiàn)了對機翼與翼身的定升減阻優(yōu)化設(shè)計。 張亮等[5]使用離散伴隨法對列車頭的氣動阻力與升力進(jìn)行了優(yōu)化。 Han 等[6]和Miretti 等[7]均使用離散伴隨法經(jīng)過多輪迭代實現(xiàn)對整車進(jìn)行氣動減阻。

離散伴隨法的基礎(chǔ)是梯度優(yōu)化法，可得到目標(biāo)函數(shù)的靈敏度，但無法確定最優(yōu)步長且容易陷入局部最優(yōu)。 使用啟發(fā)式算法，可快速得到目標(biāo)函數(shù)的全局最優(yōu)解。 汪怡平等[8]和張英朝等[9]均使用啟發(fā)式算法對代理模型進(jìn)行全局尋優(yōu)，以實現(xiàn)對汽車外形的減阻。 楊易等[10]使用啟發(fā)式算法對汽車的后擾流板進(jìn)行了氣動減阻優(yōu)化。 雖然使用該類優(yōu)化法可快速得到代理模型的全局最優(yōu)解，但建立代理模型的變量選取須依靠經(jīng)驗，存在盲目性。

結(jié)合離散伴隨法與啟發(fā)式算法的優(yōu)點，以某量產(chǎn)SUV 為研究對象，通過離散伴隨法分析汽車表面關(guān)于風(fēng)阻系數(shù)Cd的靈敏度，進(jìn)而確定優(yōu)化部位為前保險杠、后視鏡、尾翼、后保險杠等區(qū)域。 為每一個優(yōu)化部位設(shè)置變形控制體，并以控制體節(jié)點的變形位移作為設(shè)計變量，結(jié)合哈默斯雷試驗設(shè)計法構(gòu)建設(shè)計變量的樣本點空間。 采用自由變形(FFD)技術(shù)得到各樣本點模型并計算出與之對應(yīng)的風(fēng)阻系數(shù)值，再利用Kriging 插值法建立代理模型。 選擇啟發(fā)式算法中的多島遺傳算法對代理模型進(jìn)行全局尋優(yōu)，最終獲得較顯著的減阻效果。 該方法解決了單獨使用離散伴隨法不能快速得到全局最優(yōu)解，以及單獨使用啟發(fā)式算法尋優(yōu)代理模型時優(yōu)化變量選取的盲目性的問題，為以后汽車的優(yōu)化改型提供參考。

1 湍流模型與優(yōu)化方法

1.1 Realizable k-ε 湍流模型

在進(jìn)行穩(wěn)態(tài)汽車外流場模擬時，使用的控制方程為雷諾平均方程，此方程是根據(jù)N-S 方程通過統(tǒng)計平均方法所導(dǎo)出的湍流統(tǒng)計方程，須使用湍流模型來使其封閉[11]。 使用 Realizablek-ε湍流模型來封閉雷諾平均方程，其輸運方程表達(dá)式[12]如下:

式中:k為湍流動能；ε為湍流耗散率；為平均速度；μ為動力黏度；ρ為流體密度；YM為壓力修正項；υ為運動黏度；Pk和Pb為湍動能為生成項；Sk和S∈為源項為湍流黏度為平均應(yīng)變率張量[13]。 其它常數(shù)值為:C1∈＝1.44；σ∈＝1.2；σk ＝1.0；C2＝1.9；Cμ ＝0.09；當(dāng)Pb≥ 0 時，C3∈＝1，當(dāng)Pb＜0 時，C3∈＝0。

1.2 氣動優(yōu)化的離散伴隨法原理

在汽車氣動優(yōu)化過程中，目標(biāo)函數(shù)風(fēng)阻力D可寫為網(wǎng)格節(jié)點上流場變量U與設(shè)計變量x的二元函數(shù)，即

同時，U與x滿足控制方程在求解過程中殘差為 0，即

上述優(yōu)化問題可表述為:以控制方程殘差值為約束，尋找風(fēng)阻力D的極小值所對應(yīng)的x值的約束優(yōu)化問題，即

引入拉格朗日算子λ，將上述約束優(yōu)化問題變?yōu)闊o約束優(yōu)化問題，即

將D關(guān)于x求偏導(dǎo)并整理可得

拉格朗日算子λ可為任意值，滿足

則式(7) 變?yōu)?/p>

式(8) 即為伴隨方程，式(9) 為靈敏度關(guān)系，λ為伴隨變量。

式(9) 表明目標(biāo)函數(shù)D對設(shè)計變量x的靈敏度關(guān)系不再依賴于而只需要先求解一次流動控制方程式(4)，然后再求解一次伴隨方程式(9) 即可得到D關(guān)于x的靈敏度關(guān)系，且計算量與設(shè)計變量x的數(shù)目無關(guān)[14]。

1.3 哈默斯雷試驗設(shè)計原理

哈默斯雷試驗設(shè)計屬于準(zhǔn)蒙特卡洛法的一種，該方法利用基于哈默斯雷點的準(zhǔn)隨機數(shù)生成器來進(jìn)行均勻采樣。 哈默斯雷點為設(shè)計空間內(nèi)一定數(shù)量下的一組均勻采樣點。 對n維設(shè)計變量，運用哈默斯雷試驗設(shè)計從設(shè)計空間中抽取N個樣本點，第p個樣本點可以表示如下[15]:

式中Ri為大于等于 2 的素數(shù)(2、3、5、7…)，且R1＜R2＜…＜Rn-1。 函數(shù)關(guān)系φR的計算方式如下:

式中ai為[0，R-1]上的整數(shù)。 非負(fù)整數(shù)p可以展開為關(guān)于基R的多項式:

其中m ＝int[logRN]

2 仿真分析與試驗驗證

2.1 整車模型與計算域建立

根據(jù)實車尺寸建立整車的幾何模型，打開進(jìn)氣格柵，保留汽車底盤原有幾何特征，忽略了雨刮、輪胎花紋、線束和螺釘?shù)燃?xì)小幾何特征，如圖1(a)所示。 再建立出計算流體動力學(xué)(CFD)仿真所需的計算域，計算域的尺寸為長11L，寬9W，高5H，其中L、W、H分別為汽車的長、寬、高，如圖1(b)所示。

圖1 整車幾何模型與仿真計算域

機艙內(nèi)的冷卻系統(tǒng)——冷凝器、中冷器、散熱器，使用多孔介質(zhì)模型來進(jìn)行模擬。 冷卻系統(tǒng)試驗數(shù)據(jù)如表1 所示，表中v為氣流速度， Δp為壓力損失。

表1 冷卻系統(tǒng)試驗數(shù)據(jù)

由表1 中的試驗數(shù)據(jù)，根據(jù)公式?p ＝L(Piv2＋Pυv) 擬合出冷卻系統(tǒng)的多孔慣性系數(shù)Pi與多孔黏性系數(shù)Pυ，如表2 所示，式中L為特征長度。

表2 冷卻系統(tǒng)擬合數(shù)據(jù)

2.2 體網(wǎng)格劃分

在實際CFD 仿真計算中，六面體網(wǎng)格相對于多面體與四面體網(wǎng)格來說，具有網(wǎng)格質(zhì)量好、收斂快且精度高的特點[16]。 所以本文中選擇六面體網(wǎng)格填充計算域，并使用棱柱層網(wǎng)格來模擬近壁面邊界層的流動。

壁面第一層網(wǎng)格的厚度對計算結(jié)果的精度和殘差收斂性有著較大的影響。 在CFD 仿真中，常用壁面Y＋值來檢驗第一層網(wǎng)格厚度是否合理。Y＋是第一層網(wǎng)格節(jié)點的無量綱化壁面距離。 由于汽車外流場是一個多尺度的復(fù)雜流動過程，既有高雷諾數(shù)的流動也有低雷諾數(shù)的流動，于是在軟件Starccm＋中選擇“全Y＋”的壁面處理方式來自適應(yīng)地模擬汽車外表近壁面的流動情況。 試驗測量和數(shù)值模擬結(jié)果表明:對于低雷諾數(shù)流動情況，需要Y＋＜5 時才能精確模擬近壁面流動情況[11]；對于高雷諾數(shù)流動情況，為了減少近壁面網(wǎng)格層數(shù)，須使用壁面函數(shù)來模擬近壁面流動，此時推薦Y＋位于30 ～100 之間，且靠近30 一端。 考慮到汽車外流場的復(fù)雜性，文獻(xiàn)[17]指出Y＋位于0 ～300 之間即可得到滿足工程精度需要的模擬結(jié)果。 本文模型表面Y＋分布與網(wǎng)格模型如圖2 和圖3 所示。

圖2 模型表面Y＋分布云圖

圖3 計算域網(wǎng)格模型

2.3 邊界條件與求解設(shè)置

依據(jù)風(fēng)洞試驗條件進(jìn)行CFD 仿真邊界條件的設(shè)置。 計算域入口設(shè)置為速度入口，風(fēng)速為33.33 m/s，出口為壓力出口，參考壓力為0。 計算域具體邊界設(shè)置如表3 所示。 之后選取 Realizablek-ε湍流模型對整車外流場進(jìn)行計算模擬，求解參數(shù)設(shè)置如表4 所示。

表3 邊界條件參數(shù)設(shè)置

表4 求解參數(shù)設(shè)置

2.4 仿真結(jié)果與試驗驗證

根據(jù)以上的設(shè)置，對整車模型進(jìn)行外流場穩(wěn)態(tài)計算，得到了風(fēng)阻系數(shù)值為0.340 1 以及汽車對稱截面上的20 個壓力測點的壓力值變化趨勢。

通過風(fēng)洞試驗對CFD 模型進(jìn)行驗證，風(fēng)洞試驗條件如下:風(fēng)洞為3/4 開口回流風(fēng)洞，噴口面積為27 m2， 試驗環(huán)境溫度為 26 ～27 ℃，相對濕度 35%，大氣壓力103.1 ～103.3 kPa，入口湍流強度小于0.3%。 試驗得到了風(fēng)阻系數(shù)值為0.342 9 以及與仿真對應(yīng)壓力測點的壓力值變化趨勢。

表5 風(fēng)阻系數(shù)仿真值與試驗值對比

風(fēng)阻系數(shù)仿真值與試驗值結(jié)果如表5 所示，表中顯示其仿真值與試驗值的誤差為0.82%，滿足工程精度要求的5%以內(nèi)。 測點壓力的仿真值與試驗值的變化趨勢如圖4 所示，圖中顯示其仿真值與試驗值的變化趨勢一致。 因此，認(rèn)為模型精度可靠，且仿真結(jié)果有效。

圖4 測點壓力值變化趨勢對比圖

2.5 氣動特性分析

汽車頭部與尾部的壓力分布如圖5 所示，在汽車頭部存在著較大的正壓區(qū)分布，而在汽車尾部則全是負(fù)壓區(qū)，且在后保險杠兩側(cè)存在著兩處異常的低壓區(qū)。 頭部與尾部的兩股力同時作用阻礙了汽車的前行。

圖5 汽車表面壓力分布云圖

汽車對稱截面上的湍動能分布如圖6 所示，在汽車尾流區(qū)存在著較大的湍動能區(qū)域，該區(qū)域消耗汽車行駛能量，阻礙汽車前行。

圖6 湍動能云圖

3 氣動減阻優(yōu)化

3.1 靈敏度識別

由1.2 節(jié)中的理論可知，通過離散伴隨法可識別出汽車表面關(guān)于風(fēng)阻系數(shù)的靈敏度。 首先對汽車外流場進(jìn)行一次耦合流初始計算，之后在初始計算的基礎(chǔ)之上進(jìn)一步求解伴隨方程，即可得到汽車表面關(guān)于風(fēng)阻系數(shù)的靈敏度關(guān)系。

靈敏度識別結(jié)果如圖7 所示，圈選區(qū)域為綜合考慮之后選取的優(yōu)化部位，除了尾翼外，其它部位均為對稱選取，在圖中只標(biāo)出了一側(cè)。 圖中區(qū)域A 代表沿著表面的正法線方向移動可以減小風(fēng)阻系數(shù)，區(qū)域B 則反之。

圖7 風(fēng)阻系數(shù)靈敏度云圖

3.2 設(shè)計變量選取

使用軟件Hypermesh 的網(wǎng)格變形功能，在選取的優(yōu)化部位上建立變形控制體，并通過tcl 腳本驅(qū)動變形控制體的控制節(jié)點進(jìn)行自動變形。 由于在區(qū)域4 處添加變形控制體會與區(qū)域5 處形成干涉，于是僅考慮在此處加裝側(cè)向擾流板，而不對其進(jìn)行變形。 變形體設(shè)置的具體情況如圖8所示。

圖8 變形控制體設(shè)置(圖中7表示垂直平面向外，?反之)

圖8 中黑色箭頭所指的方向為變形過程中控制節(jié)點移動的方向，X1～X6為本次優(yōu)化所選取的6 個設(shè)計變量，設(shè)計變量的取值范圍如表6所示。

表6 變量取值范圍 mm

3.3 試驗設(shè)計

常用的試驗設(shè)計方法有全因子設(shè)計、正交試驗設(shè)計、拉丁超立方試驗設(shè)計和哈默斯雷試驗設(shè)計等。對于相同的試驗水平和因素，前兩種試驗設(shè)計方法相較于后兩者需要更多的試驗次數(shù)，而哈默斯雷試驗設(shè)計在處理高維非線性問題時，比拉丁超立方試驗設(shè)計采樣點的空間分布更加均勻[18]。

根據(jù)設(shè)計變量的取值范圍，采用哈默斯雷試驗設(shè)計方法抽取30 組樣本點，并據(jù)此編制tcl 腳本，用于驅(qū)動變形體控制節(jié)點按照樣本點取值進(jìn)行移動變形。 分別對變形得到的30 個樣本模型進(jìn)行仿真計算，得到與之對應(yīng)的風(fēng)阻系數(shù)Cd值，結(jié)果如表7所示。

表7 試驗設(shè)計樣本點與結(jié)果

3.4 代理模型建立與全局尋優(yōu)

根據(jù)表7 中數(shù)據(jù)值，采用Kriging 插值法來建立代理模型。 代理模型的精確性與泛化性可采用交叉驗證集的決定性系數(shù)R2來進(jìn)行評價，R2值越接近1表明模型擬合的精度越高，但此時模型很有可能出現(xiàn)過擬合的情況，導(dǎo)致泛化性變差。 因此，認(rèn)為代理模型的R2＞0.8 即滿足精度要求，且R2能在0.9～1之間最好。

在本次代理模型建立的過程中，隨機選擇樣本點的20%，即6 個作為交叉驗證集，剩下的24 個點用于模型的擬合。 擬合結(jié)果顯示，交叉驗證集的決定性系數(shù)R2＝0.92。 因此，認(rèn)為代理模型的精確性與泛化性滿足要求。

為了得到設(shè)計變量的最優(yōu)組合，選擇多島遺傳算法來對代理模型進(jìn)行全局尋優(yōu)。 多島遺傳算法參數(shù)設(shè)置如表8 所示，優(yōu)化的目標(biāo)函數(shù)與設(shè)計空間如下。

目標(biāo)函數(shù):min{Cd}

設(shè)計空間:0≤X1≤10，0≤X2≤10，0≤X3≤18，0≤X4≤15，0≤X5≤100，0≤X6≤100

表8 多島遺傳算法參數(shù)設(shè)置

根據(jù)表8 的參數(shù)設(shè)置，對代理模型進(jìn)行全局尋優(yōu)，得到代理模型預(yù)測的最優(yōu)結(jié)果組合以及最小風(fēng)阻系數(shù)值。 將預(yù)測的最優(yōu)結(jié)果組合帶入CFD 模型中進(jìn)行仿真計算，得到實際的最小風(fēng)阻系數(shù)值，具體數(shù)據(jù)如表9 所示。

表9 優(yōu)化結(jié)果

由表可見，最小風(fēng)阻系數(shù)預(yù)測值與計算值的誤差僅為0.15%，這進(jìn)一步證明了代理模型以及優(yōu)化結(jié)果的可靠性。 風(fēng)阻系數(shù)仿真值在優(yōu)化前為0.340 1， 優(yōu)化后 為 0.328 9， 下 降 了0.011 2， 約 為 11 count(1 count＝0.001)， 減阻率為 3.29%。

3.5 優(yōu)化結(jié)果分析

優(yōu)化前后汽車尾部的壓力分布如圖9 所示。 由圖可見，優(yōu)化后汽車后窗、后背門等處的壓力恢復(fù)較優(yōu)化前有所回升，分布在后保險杠兩側(cè)的異常低壓區(qū)域也已經(jīng)消失。

圖9 優(yōu)化前后尾部壓力云圖對比

優(yōu)化前后汽車對稱截面湍動能分布如圖10 所示，可見優(yōu)化后汽車尾部湍流動能區(qū)域有明顯減小，表明優(yōu)化后汽車尾流對汽車能量的消耗減少，從而風(fēng)阻減小。

圖10 優(yōu)化前后湍動能云圖對比

被優(yōu)化區(qū)域在優(yōu)化前后的輪廓形狀如圖11 所示，虛線框為優(yōu)化前輪廓，實線框為優(yōu)化后輪廓。 由圖可見，被優(yōu)化區(qū)域的實際改變量較小，未對汽車的外型風(fēng)格造成明顯改變。 表明本優(yōu)化策略可以很好地應(yīng)用于量產(chǎn)車的氣動優(yōu)化。

圖11 優(yōu)化前后輪廓對比

4 結(jié)論

以某量產(chǎn)SUV 為研究對象，對其進(jìn)行了氣動分析，之后利用離散伴隨法結(jié)合哈默斯雷試驗設(shè)計、FFD技術(shù)和Kriging 插值法對其進(jìn)行了靈敏度識別與代理模型建立，最后采用多島遺傳算法對代理模型進(jìn)行了全局優(yōu)化。 根據(jù)靈敏度識別結(jié)果，通過添加側(cè)向擾流板并對6 個設(shè)計變量進(jìn)行全局尋優(yōu)，最終使整車風(fēng)阻系數(shù)下降了11 count，減阻率達(dá)3.29%。

所使用的優(yōu)化策略可很好地滿足僅小幅改動汽車外形卻須明顯減阻的實際需求。 同時，還解決了單獨使用離散伴隨法時不能快速得到全局最優(yōu)解，以及單獨使用啟發(fā)式算法尋優(yōu)代理模型時優(yōu)化變量選取的盲目性問題，為快速有效確定優(yōu)化方案提供了理論支持。