辛俐 胡興軍 張靖龍 王靖宇 高炳釗

(吉林大學(xué) 汽車仿真與控制國家重點實驗室，吉林 長春 130025)

汽車在雨中行駛時，積聚到側(cè)窗玻璃上的雨水形成水膜或小水流，折射和阻擋光線，嚴重損害駕駛的安全性和舒適性[1- 2]。傳統(tǒng)的汽車側(cè)窗水污染問題在實車試驗階段才能發(fā)現(xiàn)，加大了新車開發(fā)的時間和投資成本。因此，在汽車設(shè)計初期對側(cè)窗水污染問題進行研究并提出有效的解決措施十分重要。

側(cè)窗區(qū)域的水污染主要來自：A柱溢流；后視鏡表面水膜隨尾渦脫落撞擊到側(cè)窗玻璃；來流中夾帶的雨滴[3]。汽車A柱的造型和幾何特征對A柱溢流影響較大，因此，通過對A柱造型不斷優(yōu)化，可以有效抑制A柱溢流，提高側(cè)窗視野的清晰度[4]。Foucart等[5]利用仿真模擬了汽車在雨中行駛時車身表面水膜流動的路徑，并利用風(fēng)洞試驗驗證了數(shù)值計算方法的準確性。Gaylard等[6]在A柱上安裝楔形導(dǎo)流裝置，改變氣流和水膜的流動狀態(tài)，減少了A柱邊緣處的雨水積聚。Dasarathan等[7]優(yōu)化了側(cè)窗區(qū)域上方水槽的形式，有效的減少了側(cè)窗污染的面積。Jilesen等[8]通過數(shù)值方法預(yù)測A柱溢流發(fā)生的位置，進一步表明A柱截面形狀的不同會引起氣流剪切力的不同，從而影響側(cè)窗污染模式。國內(nèi)對汽車水污染研究起步較晚，僅辛俐等[9]利用格子玻爾茲曼和拉格朗日方法對側(cè)窗水污染進行了模擬，并對比分析了4種不同A柱截面形狀對側(cè)窗清晰度的影響。從以上研究可知，當(dāng)前對側(cè)窗水污染的A柱優(yōu)化多采用試湊法：首先設(shè)計出不同方案的A柱，然后通過風(fēng)洞試驗或數(shù)值模擬方法對各方案側(cè)窗水污染進行比較，選出最佳A柱形狀。該方法對設(shè)計者的工程經(jīng)驗要求較高，且具有一定的盲目性。通過數(shù)學(xué)方法在一定范圍內(nèi)對設(shè)計變量進行自動尋優(yōu)設(shè)計已經(jīng)廣泛應(yīng)用于汽車空氣動力學(xué)減阻和降噪方面，并取得了較為成功的經(jīng)驗?；诖?，文中嘗試將該數(shù)學(xué)方法結(jié)合子域計算方法應(yīng)用于汽車側(cè)窗水污染處理。

文中以DrivAer為研究對象，選取A柱截面H、L和W3個參數(shù)為設(shè)計變量，采用正交實驗設(shè)計方法建立樣本空間，運用網(wǎng)格變形技術(shù)獲得各樣本點的模型并通過子域仿真方法計算獲得樣本點的側(cè)窗水污染面積，選取側(cè)窗水污染面積為優(yōu)化目標(biāo)，基于Kriging近似模型和多島遺傳算法建立優(yōu)化流程。最終，側(cè)窗區(qū)域水污染面積有了顯著的減小。

1 數(shù)值計算理論

1.1 汽車外流場的計算

汽車雨水模擬首先要完成汽車外流場的計算和分析，在此基礎(chǔ)上進行水滴和水膜運動過程的研究[10]。文中采用PowerFLOW軟件進行數(shù)值模擬。該軟件外流場的計算基于格子玻爾茲曼(LBM)方法。該方法是Boltzmann方法的特殊的離散形式，建立在分子運動論和統(tǒng)計力學(xué)基礎(chǔ)上的計算流體力學(xué)方法[11]。相比于傳統(tǒng)的計算流體力學(xué)(CFD)方法，該方法將流體看成許多微小的粒子，從微觀角度獲取流體分子的運動信息，在時間和空間上完全離散，通過流體粒子的遷移和碰撞來描述流動現(xiàn)象[12]。

文中采用的是D3Q19離散速度模型，該離散模型的離散速度ξα的計算如下

ξα=

(1)

平衡態(tài)分布函數(shù)f(eq)可表示為

(2)

(3)

外流場計算過程中，根據(jù)雷諾數(shù)的大小采用不同的數(shù)值方法。在自由流區(qū)域利用LBM方法求解，在近壁面區(qū)域選用ABLM邊界層模型和DNS方法求解。利用超大渦模擬(VLES)的方法對湍流區(qū)域進行求解。對于可分辨渦結(jié)構(gòu)的大尺度區(qū)域直接求解，而小尺度區(qū)域則采用RNGk-ε方程建模得到。該超大渦模擬方法，有效提高了計算效率和計算準確度，被廣泛應(yīng)用于汽車空氣動力學(xué)和水污染計算[13]。

1.2 離散相求解

對汽車外流場求解穩(wěn)定后，得到初始流場，然后采用拉格朗日離散相模型模擬水滴運動。該方法將雨滴定義為球形顆粒，通過求解雨滴的運動方程來計算其運動的路線[14- 15]。單個雨滴顆粒的運動方程為

(4)

(5)

式中，v為來流速度，vp為顆粒速度，R為雨滴顆粒的半徑，mp為粒子的質(zhì)量，g為重力加速度，μ為流體動力粘度，ρ為流體密度，CD為雨滴阻力系數(shù)，Re為相對雷諾數(shù)，

(6)

(7)

此外，分別采用O’Rourke顆碰撞模型和泰勒類比(TAB)液滴破碎模型計算顆粒運動過程中的碰撞和破碎。同時，采用歐拉壁面液膜模型對水膜進行模擬。

2 數(shù)值方案設(shè)置

2.1 整車模型及建立計算域

文中的研究對象為DrivAer模型，如圖1所示。該模型是依據(jù)實車建立的具有整車細節(jié)的空氣動力學(xué)模型。比通用的空氣動力學(xué)簡化模型更適用于研究復(fù)雜的汽車空氣動力學(xué)現(xiàn)象[16]。采用前處理軟件ANSA中的STL格式劃分面網(wǎng)格。將網(wǎng)格與幾何平面之間的弦偏差設(shè)置為0.012 5 mm，使面網(wǎng)格與幾何表面高度貼合，盡可能的保留幾何特征，并將最大網(wǎng)格尺寸設(shè)置為15 mm。

圖1 DrivAer 幾何模型

圖2所示為仿真模擬時所構(gòu)建的數(shù)值模擬風(fēng)洞模型，其入口設(shè)置為速度入口，速度大小為27.28 m/s。出口設(shè)置為壓力出口，壓力大小為0。汽車模型表面設(shè)置為固定無滑移壁面，數(shù)值風(fēng)洞的左右兩側(cè)和頂面設(shè)置為滑移壁面。為消除地面效應(yīng)的影響，車輛前端與后端地面設(shè)置為滑移壁面，車輛周圍地面設(shè)置為非滑移壁面。

圖2 計算域及邊界條件

為控制網(wǎng)格總體數(shù)量，且能捕捉流場細節(jié)，保證仿真的精確度和可靠性，共設(shè)置10個VR區(qū)(加密區(qū))，對計算敏感區(qū)域進行加密，其中最小加密網(wǎng)格區(qū)尺寸為1.25 mm，隨著加密區(qū)的向外擴展，其尺寸也在成倍增大，最終體網(wǎng)格數(shù)量為8 397.36×104，如圖3所示。其中：圖3(a)為車體縱向中心截面的網(wǎng)格分布，圖中顯示了加密區(qū)10到加密區(qū)4。圖3(b)為Z=0.7 m截面體網(wǎng)格分布及其細節(jié)放大圖，針對后視鏡、A柱以及后視鏡關(guān)鍵區(qū)域進行了局部網(wǎng)格加密。

圖3 計算域網(wǎng)格顯示

2.2 子域計算方法

不同于汽車外流場的計算，汽車水污染的計算不僅包含多相流之間的耦合，而且更需要細化關(guān)鍵區(qū)域的網(wǎng)格，捕捉水滴和液膜的運動，這都會耗費大量的計算資源。為實現(xiàn)計算準確性與計算效率之間的平衡，文中采用車身水污染子域計算方法，計算流程如圖4所示。

圖4 子域計算流程

將圖5所示的側(cè)窗水污染關(guān)鍵區(qū)域從計算域中提取出來作為子域。子域的進口離前擋風(fēng)玻璃約 2 m，側(cè)面包含整個側(cè)窗后視鏡尾渦，距離側(cè)窗玻璃約0.8 m，出口包含整個后視鏡的尾渦到達B柱位置。重新細化子域中網(wǎng)格，將其邊界條件劃分為3個區(qū)域來確保計算準確性，每個區(qū)域的信息都是來自完整計算該區(qū)域時流場的平均信息的映射。子域的出口設(shè)置為壓力出口，流動區(qū)域為整個區(qū)域的相同位置的壓力和速度平均。子域的入口為速度入口，并且在該入口處設(shè)置粒子發(fā)射器。本次模擬大雨工況下的側(cè)窗水污染情況，粒子發(fā)射器雨量設(shè)置為17 mm/h，液滴顆粒粘度為0.001 Pa·s，表面張力為0.072 8 N/m2，顆粒密度為1 000 kg/m3。在粒子運動過程中受重力影響，g=9.8 m/s2。液滴顆粒直徑分布服從Gaussian分布，最大顆粒為3 mm，最小顆粒為1 mm。液滴的水平初速度設(shè)置為27.78 m/s，與汽車行駛速度相同。當(dāng)車身表面薄膜厚度超過0.3 mm時，水膜二次破碎成液滴，重新進入流場中。

圖5 子域計算域

2.3 試驗驗證

2.3.1 汽車外流場仿真結(jié)果及驗證

汽車外流場求解的準確性對雨滴以及水膜運動有著重要的影響。準確的外流場計算是進行汽車水污染計算的第1步。圖6為阻力系數(shù)的發(fā)展曲線。從圖中可知，初始計算階段，阻力系數(shù)的波動較大，流場不穩(wěn)定，但隨著計算步數(shù)的增加(4×104步)，計算趨于收斂，阻力系數(shù)逐漸穩(wěn)定在固定值附近。

圖6 阻力系數(shù)曲線

為驗證汽車外流場數(shù)值仿真數(shù)據(jù)的準確性，將該車的1/4模型在吉林大學(xué)風(fēng)洞實驗室進行了風(fēng)洞試驗，如圖7所示?？諝庾枇ο到y(tǒng)通過六分力平衡系統(tǒng)進行測量。表1為不同風(fēng)速下的試驗結(jié)果與仿真計算結(jié)果的對比，計算結(jié)果和風(fēng)洞實驗結(jié)果的相對誤差在5%之內(nèi)，可見外流場的模擬準確可靠。

(a)側(cè)視圖 (b)正視圖

表1 風(fēng)洞試驗數(shù)據(jù)與仿真數(shù)據(jù)對比

2.3.2 子域計算結(jié)果驗證

文中所用整車側(cè)窗水污染仿真計算方法已于本課題組發(fā)表文獻中進行了試驗驗證[9]?；诖?，現(xiàn)只需驗證子域計算方法與整車側(cè)窗水污染計算時的差異性。將子域計算所得后視鏡尾渦的水相分布和側(cè)窗區(qū)域水相分布與整車計算結(jié)果進行了對比，如圖8和圖9所示?？梢?，子域計算結(jié)果與整車計算結(jié)果水相的分布具有良好的一致性。由于文中采用瞬態(tài)仿真計算方法，存在微小的差異是正常的仿真現(xiàn)象。

圖8 后視鏡尾渦水相分布圖

(a)子域計算

(b)整車計算

對整個計算域進行水管理計算時，需花費 12 288 h(單個CPU)，利用子域方法進行計算時，計算時間僅僅需2 304 h(單個CPU)?？梢?，子域計算時間是整個計算時間的18.75%，子域計算方法可以在保證精度的同時大幅度提高計算效率。

3 側(cè)窗水污染優(yōu)化

3.1 設(shè)計變量選取及試驗設(shè)計

在進行側(cè)窗水污染優(yōu)化之前，分析側(cè)窗水污染的主要來源及敏感區(qū)域是十分必要的。對于DrivAer模型，A柱溢流是側(cè)窗水污染的主要來源[9]。基于此，選取A柱截面的3個參數(shù)作為設(shè)計變量，如圖10所示。根據(jù)轎車A柱的大小以及A柱變形對人機工程、美學(xué)、視野盲區(qū)等因素的影響，給定設(shè)計變量的H、L、W的取值范圍分別為[0，15]，[0，18]，[5，-10]，單位為mm，對于A柱截面參數(shù)的移動，水平方向規(guī)定向左為正，向右為負。豎直方向規(guī)定向上為正，向下為負。

圖10 A柱截面示意圖

試驗設(shè)計的方法較多，一般常用的有正交試驗設(shè)計、拉丁超立方試驗設(shè)計和優(yōu)化的拉丁超立方試驗設(shè)計等[17- 18]。每種試驗設(shè)計方法的理論根據(jù)存在差異，導(dǎo)致設(shè)計變量的樣本點分布、目標(biāo)函數(shù)響應(yīng)的方式不同。因此，需要根據(jù)計算工況的情況選擇合適的試驗設(shè)計方法?；谒廴居嬎愕膹?fù)雜性，應(yīng)在設(shè)計因子均勻分布于設(shè)計空間的基礎(chǔ)上，盡可能的減小計算量。將A柱截面的3個參數(shù)定義為3個因子，并根據(jù)參數(shù)的取值范圍，對每個因子取4個水平。根據(jù)所確定的因子數(shù)和水平數(shù)，選用正交試驗設(shè)計方法可以得到合理的L16(43)試驗設(shè)計表。根據(jù)正交表中的16組樣本點取值，采用ANSA軟件的網(wǎng)格變形技術(shù)對A柱截面變形得到計算網(wǎng)格模型。然后，分別對這16個網(wǎng)格模型進行CFD仿真計算，得到側(cè)窗水污染面積，如表 2所示。其中側(cè)窗區(qū)域的污染面積Sw可以通過下式獲得，

該類型街道白天時段交通壓力持續(xù)較大，未存在明顯的高峰和非高峰，06:00—21:00交通運行狀況一直處于中度偏嚴重擁堵的狀態(tài)，交通需求持續(xù)較大，僅在中午12點左右交通狀況微弱改善，總體上看晚高峰的擁堵情況較早高峰更差. 具有這一特征的街道多為全天交通發(fā)生和吸引強度均較大的區(qū)域，如東華門街道內(nèi)有天安門廣場、故宮博物院等歷史遺址，還包括很多政府機構(gòu)、醫(yī)院、學(xué)校以及商業(yè)街，白天會持續(xù)發(fā)生和吸引大量的交通流，全天交通需求量居高不下，全天擁堵時間較長.

(8)

式中，Rw為側(cè)窗總像素，Rs為側(cè)窗污染像素數(shù)，S為側(cè)窗總面積。

表2 正交實驗表及計算結(jié)果

3.2 近似代理模型建立及可信度分析

近似代理模型可以根據(jù)設(shè)計變量和結(jié)果之間的映射關(guān)系，擬合為經(jīng)驗公式，能快速實現(xiàn)優(yōu)化設(shè)計[19]。相比于神經(jīng)網(wǎng)絡(luò)近似模型和響應(yīng)面模型，該模型可以覆蓋更多的樣本點的取值。通常由全局代理模型和局部代理模型組成[20]：

y(x)=f(x)+z(x)

(9)

式中，y(x)為待求解的響應(yīng)函數(shù)，f(x)為通過樣本點擬合的近似函數(shù)，z(x)為模擬局部偏差的近似函數(shù)，可表示為均值為0和方差為σ2的隨機分布誤差。z(x)的協(xié)方差矩陣可表示為

cov[z(xi)，z(xj)]=σ2R[R(xi，xj)]

(10)

式中，R為相關(guān)矩陣，R(xi，xj)是隨機樣本點中任意兩個樣本點的相關(guān)函數(shù)。該相關(guān)函數(shù)的種類很多，此次計算采用常用的高斯函數(shù)作為相關(guān)函數(shù)，可表示為

(11)

式中：m表示設(shè)計變量的個數(shù)；xik和xjk表示樣本點xi和xj第k個設(shè)計元素；θk表示待擬合相關(guān)函數(shù)的未知相關(guān)系數(shù)，其極大似然估計可以用來擬合Kriging模型，表示為

(12)

已構(gòu)建的近似模型的精度是否符合要求可以通過Kriging擬合值與CFD計算值之間的擬合程度來表示。圖11所示為體現(xiàn)Kriging模型所構(gòu)建代理模型精度的散點圖，其中縱坐標(biāo)為仿真計算得到的數(shù)值，橫坐標(biāo)是近似代理模型獲得的預(yù)測值，中間黑色斜線表示真實值與預(yù)測值相等。因此，散點越靠近黑色實線則表示預(yù)測值和CFD計算值的誤差越小，近似代理模型的擬合程度越高。從圖中可以看出，散點基本都落在黑色實線的附近。而且，通過計算，代理模型的決定系數(shù)R2=0.98，遠大于0.9，可見基于Kriging模型構(gòu)建的側(cè)窗水污染近似代理模型具有很高的預(yù)測精度，可以用于后續(xù)優(yōu)化計算。

圖11 近似模型預(yù)測值與CFD計算真實值的比較

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

表3給出了原始模型和優(yōu)化模型之間的設(shè)計參數(shù)和優(yōu)化目標(biāo)的比較。經(jīng)過一些系列優(yōu)化后，H增加了7.2 mm，L增加了1.2 mm，W增加了11.2 mm，與原始模型相比，側(cè)窗水污染面積下降了66.2%。

表3 原始模型與優(yōu)化模型結(jié)果對比

根據(jù)多島遺傳算法得到的最優(yōu)解，采用網(wǎng)格變形技術(shù)對原始模型進行變形，并通過CFD仿真計算，求得Sw=0.046 0 m2，相對誤差為1.7%，說明該全局優(yōu)化方法具有較高的可信度。

圖12為原始模型和優(yōu)化后模型的側(cè)窗水膜厚度分布對比。圖13為原始模型與優(yōu)化后模型的側(cè)窗區(qū)域表面剪切力分布對比。由圖12可以看出，優(yōu)化后模型A柱表面以及側(cè)窗關(guān)鍵區(qū)域的水膜明顯減小。擋風(fēng)玻璃的雨水在高速氣流剪切力和刮雨器的作用下向A柱與擋風(fēng)玻璃交界處流動。積聚在該處的水膜一部分在重力的作用下向下流動，由發(fā)動機艙處的排水系統(tǒng)排出。一部分隨氣流剪切力向側(cè)窗移動。如圖13(a)所示，原始模型A柱表面氣流剪切力較高且方向與水膜運動方向一致，水膜在氣流剪切力的驅(qū)動下，向側(cè)窗區(qū)域流動。側(cè)窗玻璃上方存在明顯的低剪切力區(qū)域，這就意味著，水膜一旦到達該處后，更易受重力的驅(qū)動向下運動，影響駕駛員的關(guān)鍵視野區(qū)。而優(yōu)化后的模型，前側(cè)A柱表面剪切力方向與水膜流動方向相反，迫使水膜積聚在該處，阻礙水膜向側(cè)窗運動，如圖13(b)所示。優(yōu)化后模型A柱頂部區(qū)域氣流剪切力方向向上，進一步阻礙了水膜在重力作用的向下流動。因此優(yōu)化后模型的側(cè)窗水污染明顯減少。

(a)原始模型

(b)優(yōu)化模型

(a)原始模型

(b)優(yōu)化模型

4 結(jié)論

(1)文中以DrivAer A柱截面尺寸H、L和W3個參數(shù)作為設(shè)計變量，選取側(cè)窗水污染面積為優(yōu)化目標(biāo)，采用正交試驗方法進行樣本點抽樣并得到對應(yīng)的CFD仿真結(jié)果，運用Kriging近似代理模型和多島遺傳算法，實現(xiàn)了側(cè)窗水污染優(yōu)化設(shè)計。

(2)與原始模型相比，優(yōu)化后模型的側(cè)窗水污染減小了66.2%。相比于傳統(tǒng)的試湊法，該方法具有顯著的成效，為汽車開發(fā)早期階段進行側(cè)窗水污染研究提供了參考。

(3)采用子域數(shù)值計算方法對側(cè)窗區(qū)域的水污染進行模擬研究，保證準確性的同時大幅度的降低了計算量，保證了為汽車水污染優(yōu)化中大量樣本點的計算提供了可行的方案。

(4)優(yōu)化后的模型，A柱表面的剪切力方向與水膜運動的方向相反。這種現(xiàn)象抑制了水膜溢流到側(cè)窗。