基于離散伴隨法與代理模型的整車(chē)氣動(dòng)阻力優(yōu)化*

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

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

前言

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

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

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

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

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

1.1 Realizable k-ε 湍流模型

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

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

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

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

同時(shí)，U與x滿足控制方程在求解過(guò)程中殘差為 0，即

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

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

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

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

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

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

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

1.3 哈默斯雷試驗(yàn)設(shè)計(jì)原理

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

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

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

其中m ＝int[logRN]

2 仿真分析與試驗(yàn)驗(yàn)證

2.1 整車(chē)模型與計(jì)算域建立

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

圖1 整車(chē)幾何模型與仿真計(jì)算域

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

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

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

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

2.2 體網(wǎng)格劃分

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

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

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

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

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

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

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

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

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

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

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

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

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

圖4 測(cè)點(diǎn)壓力值變化趨勢(shì)對(duì)比圖

2.5 氣動(dòng)特性分析

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

圖5 汽車(chē)表面壓力分布云圖

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

圖6 湍動(dòng)能云圖

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

3.1 靈敏度識(shí)別

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

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

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

3.2 設(shè)計(jì)變量選取

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

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

圖8 中黑色箭頭所指的方向?yàn)樽冃芜^(guò)程中控制節(jié)點(diǎn)移動(dòng)的方向，X1～X6為本次優(yōu)化所選取的6 個(gè)設(shè)計(jì)變量，設(shè)計(jì)變量的取值范圍如表6所示。

表6 變量取值范圍 mm

3.3 試驗(yàn)設(shè)計(jì)

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

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

表7 試驗(yàn)設(shè)計(jì)樣本點(diǎn)與結(jié)果

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

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

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

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

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

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

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

由表可見(jiàn)，最小風(fēng)阻系數(shù)預(yù)測(cè)值與計(jì)算值的誤差僅為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)化前后汽車(chē)尾部的壓力分布如圖9 所示。 由圖可見(jiàn)，優(yōu)化后汽車(chē)后窗、后背門(mén)等處的壓力恢復(fù)較優(yōu)化前有所回升，分布在后保險(xiǎn)杠兩側(cè)的異常低壓區(qū)域也已經(jīng)消失。

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

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

圖10 優(yōu)化前后湍動(dòng)能云圖對(duì)比

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

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

4 結(jié)論

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

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