基于車身繞流的低雷諾數湍流模型改進研究1)

強光林 楊 易 陳 陣 谷正氣 張 勇

(湖南大學汽車車身先進設計制造國家重點實驗室，長沙 410082)

引言

汽車車身繞流是三維效應強、非定常、多尺度的復雜流動問題.目前由于缺少完備的實驗手段能夠準確測量汽車周邊流場，揭示其機理，同時尚無足夠的計算資源可對車身繞流進行直接數值求解，因此至今人們還沒能全面、深入地掌握車身繞流規律[1]，并發展與之相適應的湍流模型.研究表明，不合理的湍流模型是數值求解誤差的最主要來源，無法準確預測車身繞流轉捩是常用湍流模型數值求解誤差的主要原因[2].

本文通過梳理汽車繞流的流動機理及物理本質，分析汽車周邊流動轉捩的類型，將汽車各個復雜局部流動等效為近似簡易模型.在此基礎上，以LRNk-ε 模型[3]為原型，通過引入流線曲率因子及響應閾值改進了該模型對湍流耗散率的過強依賴性及全應力發展預測不足等問題，提升該模型對轉捩的捕捉能力，最終獲得一種改進的低雷諾數湍流模型.同時，引入約束大渦模擬思想，結合改進的低雷諾數湍流模型，構造適用于汽車繞流瞬態計算的轉捩LRN CLES 模型.

1 汽車繞流局部近似模型

氣流流經汽車車身，在汽車前端分離，形成主分離渦，主分離渦卷起一個包裹車頭的橢球形漩渦，在頂面及側面發生再附著現象.由于地面效應，汽車底部氣流形成了特有的復雜邊界條件，流動特性沿地面垂直方向變化劇烈.解決汽車外流場這種三維效應強、非定常、會流分離區多和渦結構復雜的流動問題，需深入剖析其局部流動[4].將汽車繞流模塊化為典型局部流動，一方面有助于揭示汽車周邊流動機理；另一方面使湍流模型的改進更具針對性.

本文將汽車繞流大致分為如下幾類：汽車前緣底部類似于突縮流；后緣底部類似于半限制射流；車身幾何外形突變的位置，如后視鏡與A 柱間隙等類似于突縮、突擴流；前風擋與發動機艙拐角處類似于前臺階流動；車頂后緣幾何斷面突然擴張處，類似于后臺階流動；車頂部及發動機艙蓋板處類似于平板流動.根據上述分類，分別建立平板流、前臺階流、射流、后臺階流及突縮流模型，構建合適邊界條件，將車身繞流等效為5 類局部近似模型.

采用大渦模擬方法[5]對平板流、前臺階流、射流、后臺階流及突縮流的繞流機理進行數值研究,如圖1，常用可實現的k-ε 湍流模型[6]計算結果被引入作為對比.對比分析發現：由于層流湍流轉捩導致表面摩擦系數發生突變，LES 和可實現的k-ε 湍流模型在突變處有顯著不同，因此可推測湍流模型對轉捩的預測不足是造成計算偏差的主要原因.

圖1 近似簡易模型計算結果Fig.1 Computational results for approximate simple models

圖1 近似簡易模型計算結果(續)Fig.1 Computational results for approximate simple models(continued)

除可實現的k-ε 模型外，引入類車身平板流及類車尾后臺階流評估汽車繞流數值求解中常用湍流模型對轉捩的模擬能力，標準k-ε 模型[7]、RNGk-ε 模型[8]、k-ω SST 模型[9]、RSM 雷諾應力模型[10]及LRNk-ε 模型[11]計算結果被引入作為對比，如圖2 所示.Cf的走勢圖更加驗證了湍流模型對轉捩預測的不足是造成仿真偏差的主要原因.標準k-ε 模型及RNGk-ε 模型對層流湍流轉捩的擬合偏差較大，這是由于以上兩種模型引入了半經驗式的壁面函數[12].雖然k-wSST 模型避免了壁面函數的經驗式求解，但對邊界層轉捩依然無法準確預測.RSM 雷諾應力模型則耗費較大計算資源，且較難收斂.常用湍流模型中，對Cf預測最為準確的是LRNk-ε 模型，其低雷諾數邊界層模擬機制與轉捩較為相似，但對轉捩的預測仍然存在一定偏差，特別是類車尾后臺階流動.

綜上所述，對于汽車繞流各局部流動，轉捩是決定流場發展的關鍵因素.因此，針對汽車繞流的數值模擬，提高湍流模型適應性的關鍵為增強其對轉捩的預測能力.

2 構建改進的LRN k-ε 穩態湍流模型

圖2 常用湍流模型計算結果對比Fig.2 Comparison of the computational results of commonly used turbulence models

由于低雷諾數湍流模型將N-S 方程解析至壁面[13]，避免了經驗式壁面函數的引入，有助于揭示近壁面區流動機理，并對轉捩有巨大預測潛力，本文引入低雷諾數k-ε 湍流模型[11]，并對其進行轉捩修正.全應力限制方法[14]使k-ω 模型獲得了良好模擬轉捩的能力，將此方法應用于LRNk-ε 湍流模型.同時，引入流線曲率因子提高模型對分離流轉捩的模擬能力，構建改進的LRNk-ε 湍流模型.

2.1 基于全應力限制方法的湍流耗散率推導

Jones 和Launder[11]提出的基于時均N-S 方程的低雷諾數湍流模型如下式所示

式中，μ 為分子黏性系數，ρ 為流體密度，μt是湍流黏度，常數η0=4.377,β=0.012，Cμ=0.09，f1=1，cε2=1.92，時均應變率Si j通過cε1引入以處理高應變率及彎曲流.ReT是湍流雷諾數.

擴散系數中的分子擴散系數及湍流擴散系數反映了分子黏性的影響；附加項2μ(?k1/2/?xi)2被引入以反映黏性底層中湍動能耗散的各向異性；附加項2(μμt/ρ)()2被引入以使湍動能k的計算結果與某些實驗測定值更相符.

由于常用湍流模型對流動分離預測的不足表現為：過晚模擬出分離及過早模擬出再附著，相關實驗及DNS[15]結果也證實了這種不足[14].針對以上問題，引入全應力限制方法，消除湍流模型對湍流耗散率的過強依賴性，避免了由原始湍流模型對湍動能耗散率的過低估計帶來的負面影響，具體改進過程如下.

由湍動能k的表達式

綜合考慮全應力及雷諾應力的影響，可知

該式可限制額外湍流應力的生成.

由Boussinesq 假設可知

進一步推導可得

綜上所述，湍動能耗散率應被約束如下

ε0由原始湍動能耗散率輸運方程得出，k由原始湍動能耗散方程得出，對于大多數流動，φ 取2.5 可保證湍流模型對轉捩的良好預測能力.這項改進未引入額外輸運方程，可由原始湍流模型直接計算獲得，對湍流模型在湍流充分發展區的計算結果影響較小.

2.2 考慮流線彎曲效應的湍流黏度構造

全應力限制方法可使LRNk-ε 模型獲得良好預測轉捩的能力，然而對于分離流轉捩，還需發展更為有效的方法.相關實驗及DNS計算結果表明[14]：在流動分離區域，即使僅存在極小分離泡[16]，擾動也會被迅速放大，湍動能則會迅速增長.隨著湍動能的較快增長，即使當地雷諾數較低，流動也會迅速發生轉捩.因此，要提高湍流模型對分離流及分離流轉捩的模擬能力，須構造對流動分離較為敏感的感應量，流動分離通常伴隨著高彎曲率，因此，彎曲率被引入作為感應量.

流線s的彎曲率Cs為

式中，q為速度矢量，q=(u,v,w)；q為速度標量q=，矢量nq被定義為nq=(u/q,v/q,w/q)；Cs代表三維流線曲率的標量值，也是渦半徑的倒數值，可用于表征渦的長度尺度.Cs在3 個方向上的彎曲應變率分量分別為

構造一個關聯Cs的無量綱量Sl，將Cs引入湍流黏度，Sl可表示為

Sl的物理意義是湍流長度(時間)尺度與漩渦長度(時間)尺度的比值，也可作為湍流與漩渦之間的生成項與湍流耗散項的比值.

通過參數α，湍流黏度被表征為關聯Sl的方程

經計算對比發現，χ 應與Rl在同一數量級，取為0.1 較合適，C0取為10 較為合適.

以上改進中，流線曲率被無量綱化為Sl，Sl則被用于構造α.流動雷諾數較低時，ReT較小，當層流邊界層分離，Sl迅速增大，α 將促發轉捩，使得模型具備模擬分離流轉捩的能力.一旦層流轉捩為湍流，ReT劇增，Sl的影響就微乎其微了.流線曲率的引入結合全應力限制方法可得出能準確預測轉捩的改進的LRNk-ε 湍流模型，這個模型不僅能準確預測分離流轉捩[17]，還可準確模擬自然轉捩及旁路轉捩[18].

3 構建轉捩LRN CLES 瞬態湍流模型

汽車繞流的瞬態數值求解對計算資源的消耗巨大，特別是近壁面附近[19].混合RANS/LES 方法[20]則可針對不同的流動問題，引入與之相適應的RANS 模式[21].考慮到近壁面處LES 模式的巨大計算資源消耗，本節對RANS/LES 混合模型中的RANS 模式進行轉捩修正.考慮到發展較為成熟的CLES 模型[22]能較好地模擬圓柱繞流轉捩，將其引入以處理汽車繞流瞬態問題.同時，汽車繞流為高應變率及大曲率流動，存在較多的流動分離，轉捩現象更為復雜，將改進的LRNk-ε 模型引入以構造約束大渦模擬RANS 區域的C-SGS 模型[23]，進一步提高對轉捩的捕捉能力，最終可獲得轉捩LRN CLES模型.

3.1 約束亞格子應力模型建立

由于汽車繞流屬于低馬赫數流動，引入不可壓縮流的大渦模擬控制方程

其中，頂標～代表低通過濾，通常以Δ 為過濾尺度；ρ 為流體的密度；為濾波后大尺度渦的速度；為濾波后的壓力；υ 為運動黏性系數；τij為亞格子尺度應力，代表過濾后的小尺度渦對大尺度渦的影響，可表示為

基于大渦模擬求解的關鍵為亞格子模型[24]，其中應用最為廣泛的是由Lilly 及Liu 提出的DSLM 模型[25]，由于DSLM 模型不能較好的處理近壁面流動，此處引入約束大渦模擬思想，對RANS 區域C-SGS 模型的推導如下.

根據RANS 及LES 控制方程在數學上的高度一致性，RANS 區域總雷諾應力Ri j、總亞格子應力的時均值〈τij〉及濾波后的雷諾應力應存在以下關系

符號〈〉 代表相應物理參量的時均值，Ri j可由RANS 模化為，則通過直接求解N-S 方程得出

由此可得

總亞格子應力可分解為時均值〈τij〉及脈動值(τij)′

(τij)′可由Lilly 及Liu 提出的DSLM 模型模化為()′

由DSLM 模型確定，由此可得RANS 區域;CSGS 應力為

C-SGS 模型僅用于求解近壁面區流動或者總雷諾應力可由RANS 模化的流動區域，LES 模型則對其余流動區域直接求解.C-SGS 模型的準確性很大程度上取決于，的準確性又依賴于RANS 模型，針對汽車繞流問題，RANS 模型的準確性則取決于其對轉捩的預測能力.

3.2 基于改進的LRN k-ε 模型的CLES 建立

針對汽車繞流問題，C-SGS 模型主要作用于近壁面區，該區域內可由DSLM 模型直接計算，Ri j則由第二節提出的能準確預測轉捩的改進的LRNk-ε 模型進行計算，計算過程如下

式中，μt為湍流黏性系數；常數Cμ=0.09，Φ=2.5；Si j為時均應變率；ε0由原始湍動能耗散率輸運方程得出，k由原始湍動能輸運方程得出

其中C0=10，Cs代表三維流線曲率的標量值，也是渦半徑的倒數值，可用于表征渦的長度尺度，Cs由無量綱量數Sl引入湍流黏度.

4 改進模型驗證與分析

4.1 風洞實驗設置

實驗采用1/3 縮比的汽車模型如圖3 所示，實驗遵循相似準則.風洞試驗段尺寸為2.5 m×3 m×17 m，可進行最高1:3 的縮比汽車模型實驗，湖南大學HD-2 風洞[26]設計最大風速為58 m/s，可滿足雷諾數相似原則，進入“自模擬”狀態.

4.2 改進的LRN k-ε 穩態湍流模型驗證分析

本節將具備主流三廂轎車造型的某實車模型及其HD-2 風洞實驗數據引入以驗證改進的LRNk-ε 模型及標準化數值求解體系.

圖3 試驗模型Fig.3 Test model

4.2.1 網格類型選擇

引入課題組前期某最簡實車模型風洞實驗數據[27]，基于改進的LRNk-ε 模型，分別利用非結構網格及結構網格進行仿真分析，如圖4 所示.由表1結果可知，即在不考慮模型網格劃分難易的情況下，非結構也能較快收斂，且得到較為準確的計算結果.因此，后續實車仿真模型采用非結構網格劃分.

圖4 最簡實車模型的非結構(左)與結構網格(右)Fig.4 Unstructured(left)and structural grid(right)of the simplest automobile model

表1 結構與非結構網格計算結果Table 1 Computational results based on structural and unstructured grids

4.2.2 模型簡化及仿真參數設置

本實驗采用的實車模型較為復雜，省略前格柵、輪轂、門把手、燈組、尾氣管、B柱及車底細節，保留其余特征.簡化后的計算模型尺寸如圖5 所示，長1588 mm、寬643 mm、高505 mm，最小離地間隙50 mm.

綜合計算效率及阻塞效率，選取計算域高為5 倍模型高度，寬為6 倍的模型寬度，長為12 倍模型長度(車前端距入口3 倍車長，后端距出口8 倍邊界長度)[28].入口邊界設置為速度入口，來流速度設置為30 m/s，湍流度為0.028 0，出口邊界設置為壓力出口，相對壓力為零.計算域頂部及兩側壁面設置為對稱邊界，計算域底部設置為移動壁面，汽車模型壁面設置為固壁無滑移壁面，選取改進的LRNk-ε 和可實現的k-ε 模型進行仿真計算.具體的仿真計算域及邊界條件設置如圖5 所示.控制方程離散格式設置為二階迎風格式，壓力與速度耦合設置為SIMPLE 隱式求解法.

圖5 模型試驗尺寸和仿真計算域及邊界條件Fig.5 Model size,simulation computational domain and boundary conditions

4.2.3 網格無關性驗證

仿真采用非結構網格劃分，分別設置3 種網格劃分策略，驗證網格無關性.流動急劇變化部位被加密，邊界層取10 層，增長率為1.1，經數次調整第一層網格高度，大部分模型表面y+值均分布于合理范圍內，此時第一層網格高度為0.012 mm.計算在8 個8 核CPU、64GB 內存的DELL Power 服務器上進行計算.

表2 不同網格阻力系數和計算時間Table 2 Drag coefficient and computational time based on different grid numbers

綜合考慮表2 中計算結果的誤差和收斂時間，選擇第2 種網格劃分策略，網格劃分情況如圖6所示.

圖6 該實車模型計算域網格Fig.6 Model computing grids

4.2.4 計算結果分析

收斂耗時25 h.保證其他求解條件不變，基于可實現的k-ε 模型的計算結果及HD-2 風洞試驗數據被引入作為對比.計算及實驗所得該實車模型的氣動阻力系數Cd值及氣動升力系數Cl值如表3 及圖7 所示.

表3 氣動阻力系數及氣動升力系數(β=0°)Table 3 Aerodynamic drag coefficient and aerodynamic lift coefficient(β=0°)

由圖表可知，本文提出的改進的LRNk-ε 模型可較好的計算該模型的氣動阻力及氣動升力，與風洞試驗結果誤差較小，特別是對汽車氣動升力的計算，較常用的可實現的k-ε 模型有較大提升.

圖7 實車模型氣動阻力與升力系數Fig.7 Aerodynamic drag and lift coefficient of the model

圖8～圖11 分別表示改進的 LRNk-ε 模型計算所得表面靜壓分布、車頭及車尾處流場分布，該結果與風洞試驗較為吻合.由該實車模型表面壓力分布圖8 及車頭流線圖9 可知，氣流最先到達汽車頭部，受到阻滯作用形成駐點，產生較大正壓區，并貼附汽車表面形成邊界層.之后，氣流在行駛至發動機艙蓋板的過程中，受正壓梯度影響，與前臺階流較為相似.氣流流至汽車尾部的氣流邊界層則類似于后臺階流.本文汽車尾部過渡曲率較小，逆壓梯度作用較小，邊界層行駛至汽車尾端才脫離物面，并形成自由剪切流，在尾部形成漩渦，之后匯入尾流，實驗與仿真在汽車尾部均出現兩個漩渦，實驗對比結果較好.圖10 和圖11 中可以看出實驗與仿真之間仍然有一些偏差，一方面是由于實驗條件的限制，實驗與仿真保證了相似準則，但無法確保完全一致，并且實驗會受到不可控的一些外部因素的干擾；另一方面是由于本文提出的改進的模型仍然基于一部分經驗公式，因此還需經過多個算例進行驗證，發現其不足，并進行修正，提高其普適性.

圖8 表面壓力分布Fig.8 Surface pressure distribution

圖8 表面壓力分布(續)Fig.8 Surface pressure distribution(continued)

圖9 車頭流線圖Fig.9 Streamline diagram at the head

圖10 車頭渦量圖Fig.10 Vorticity diagram at the head

圖11 車尾流線圖Fig.11 Streamline diagram at the tail

4.3 轉捩LRN CLES 瞬態湍流模型驗證分析

汽車風振噪聲產生于氣流的脈動壓力，具有較強的瞬態特性，是典型的瞬態汽車繞流問題.汽車周邊脈動壓力的數值求解對湍流模型要求較高，本節引入某款轎車進行的實車道路試驗驗證轉捩LRN CLES 模型對汽車瞬態繞流的求解能力.

4.3.1 模型簡化及仿真參數設置

簡化后的計算模型如圖12 所示，具體尺寸如圖13，長4682 mm，寬1879 mm，高1624 mm，最小離地間隙150 mm.

對車輛進行仿真的外部計算域為長方體，參數設置和4.2.2 節一致.入口邊界設置為速度入口，來流速度設置為30 m/s，湍流度為0.024 3，出口邊界設置為壓力出口，相對壓力為零.計算域頂部及兩側壁面設置為對稱邊界，計算域底部(即地面)設置為移動壁面，汽車模型壁面設置為固壁無滑移壁面，與試驗相對應，將左后側窗設置為Int Wall，模擬該側窗的開啟環境.選取轉捩LRN CLES 模型進行仿真計算.具體的仿真計算域及邊界條件設置如圖13 所示.

對于瞬態求解，控制方程離散格式設置為邊界中心差分格式，壓力與速度耦合設置為PISO 隱式求解法.

圖12 汽車風振噪聲數值計算模型Fig.12 Numerical calculation for buffeting noise of the automobile model

圖13 模型尺寸和計算域及邊界條件Fig.13 Model size,calculation domain and boundary conditions

4.3.2 網格劃分

采用非結構網格劃分后，該實車模型模擬網格總數約為2800 萬，網格設置與4.2.3 節情況相同，網格劃分情況如圖14 所示.

4.3.3 計算結果分析

此瞬態計算時間步長取為1.0×10-4s，每個時間步內迭代計算20 次.由道路實驗結果可知，聲壓級(SPL)第一個階次峰值出現的頻率為17.6 Hz，因此周期約為0.056 8 s，為實現充分瞬態效應，取流動總時長為2.1 s.計算總耗時32 d 3 h.

圖14 該實車模型氣動噪聲計算域網格Fig.14 Aerodynamic noise computing grids

本課題組前期研究中的SAS計算結果[29]被引入作為對比.計算結果如圖15 所示，計算所得前3 個階次聲壓級峰值及其出現頻率如表4 所示.由圖表可知，轉捩LRN CLES 模型及SAS 模型[30]計算所得聲壓級峰值均低于實驗值，這是由于數值計算僅考慮了計算域流場環境，而道路實驗受諸多周邊環境影響，背景噪聲較大.但轉捩LRN CLES 模型與實驗結果更為接近.驗證了轉捩LRN CLES 模型良好的瞬態性能.

圖15 仿真計算和實驗聲壓級Fig.15 SPL of the simulation and test

表4 聲壓級峰值頻率Table 4 Peak frequency of SPL

5 結論

本文首先通過典型局部流動對汽車繞流機理進行深入分析，研究發現轉捩是決定流場發展的關鍵因素.針對汽車繞流穩態及瞬態求解，分別提出了一種能夠準確預測轉捩現象的改進的LRNk-ε 及轉捩LRN CLES 模型，并將其應用于實車模型繞流的數值研究以驗證其準確性.

通過對某實車氣動力及風振噪聲的數值研究，并與HD-2 風洞試驗及實車道路試驗進行了對比分析，結果表明：相同網格條件下，相對于常用的湍流模型，基于本文提出的改進的LRNk-ε 及轉捩LRN CLES 湍流模型的標準化高精度數值求解體系能夠更準確模擬復雜實車的繞流特性，特別是對氣動升力的數值計算，精度高于常用湍流模型.