基于本征正交分解的缸內流場仿真校核方法探究

顧鵬，葛鵬輝，許敏

(上海交通大學汽車工程研究院，上海 200240)

活塞點燃直噴式(Spark Ignition Direct Injection)汽油機相比進氣道噴射汽油機在燃油經濟性和顆粒物排放方面具有很大的優勢，并得到了廣泛的研究和應用[1-2]。研究表明，缸內流場對發動機燃油霧化、火焰傳播路徑和速率具有重要的影響[3]，而缸內流場結構又受到諸多因素如進氣道形狀、進氣歧管傾角和缸蓋結構等的影響[4]。如果能夠開發出較為準確的發動機三維模型，這將極大地節約發動機開發成本。因此，發動機的三維模型校核也顯得尤為重要。

與傳統的RANS仿真不同的是，大渦模擬和試驗的結果都存在著強烈的循環差異，對發動機的三維模型校核帶來了極大的障礙。而本征正交分解(Proper orthogonal decomposition)通過基函數(mode)為大渦模擬和試驗結果提供了一個比較的基礎[5]。Lumley 等[6]最早將本征正交分解的方法引入到湍流場的分析中，將湍流場中能量較大的大尺度渦團結構分離提取出來。由于發動機缸內流場同樣具有高度湍流特性，因此該方法逐漸被用于發動機的流場分析中。Sirovich[7]提出的“snapshot”的本征正交分解方法，在與原先的直接法等效的情況下，大大提升了計算效率，也廣泛應用于流場分析。在發動機流場試驗和仿真分析方面，很多研究者都基于不同的試驗和仿真工況做了很多有益的工作。秦文瑾、許敏等[8]證明了本征正交分解可以作為分析發動機缸內流場模擬結果準確性的有效工具。Fogleman等[9]通過相位固定的本征正交分解的方法，對不同曲軸轉角的缸內流場變化進行了分析。陳豪等[10]指出目前對實際發動機流場結果本征正交分解的研究相對匱乏，并闡述了試驗與仿真缸內流場比較的應用思路。此外，一些研究者對試驗數據的完善，對缸內流場仿真校核的工作也具有重要意義。Druault等[11]通過本征正交分解對PIV的流場結果進行時間跨度上的插值，Borée 等[12]通過本征正交分解的方法闡述了時間和空間位置上的關聯。

本研究在前人工作的基礎上，著眼于基于本征正交分解的缸內流場仿真校核方法，從能量和模態的角度，全面合理地評價試驗和仿真缸內流場結果的差異性，借此提出并建立缸內流場仿真校核的流程策略。

1 數值模型描述

1.1 發動機仿真方法

關于缸內流場循環的分析，大量基礎性的研究是基于雷諾分解，建立對湍流脈動量的分析。雷諾方程的核心價值在于，將對湍流的研究轉化為了建立可以求解的、與雷諾應力這個變量相關的輸運方程。1887年提出的Boussinesq假設，通過與分子運動動量傳遞的類比，將雷諾應力轉化為渦黏性系數，而渦黏性系數進一步地被分解為對湍流渦團的長度和速度尺度這兩個變量的求解。在后續的缸內流場研究中區分零方程模型、單方程模型、多方程模型等，一定意義上就是在于求解這兩個變量的微分方程的個數。Boussinesq假設之后，混合長度理論卡門相似理論等將渦團的長度和速度尺度以經驗或者代數的方法表達，在早期的缸內流場分析中起到了很大的推動作用，但是它在時間、空間變化上的表征仍存在許多不足。因此，從1942年Kolmogorov和1945年Prandtl的研究開始，人們逐漸開始建立κ方程、渦黏度方程(Spalart-Allmaras)等單方程模型。現在主流的CFD軟件中，頻繁使用κ-ε,κ-ω等雙方程模型，對于這些模型，研究者們還在不斷地根據具體的使用情況進行修正，例如應用較廣的壓縮性修正、強旋流修正等。此外，在強調各向異性的缸內流場分析中，對雷諾應力直接進行分析的RSM模型也應用極廣。

近幾十年來，研究者們逐漸認識到缸內流場可以分解為一系列的有序流場結構，現在商業應用極廣的LES(大渦模擬)實際上就是基于空間尺度的劃分，將對應于不同尺度的渦團劃分為可解尺度和亞網格尺度，兩者分別使用直接求解的方法和假設模型的方法進行分析。大渦模擬的核心前提是大尺度渦團與仿真時設置的初始和邊界條件相關性較大，而小尺度渦團受它們的影響較小，大尺度渦團的各向異性和小尺度渦團的各向同性促成了大渦模擬將它們分開處理的合理性。

大渦模擬的思路如下：首先消除小渦的脈動影響，得出局部的平均場，即為大尺度場，而亞網格尺度場則為原始流場與該平均場的差值。大尺度場根據網格，由N-S方程直接求解，而對于亞網格尺度場，常用的模型有Smagorinsky模型、WALE模型、亞網格單方程模型以及動態結構模型。其中，對于Smagorinsky模型和亞網格單方程模型，在后續LES的研究和計算中模型系數進行了修正，以滿足更多的情況，進而產生了動態Smagorinsky模型[13]和亞網格單方程模型。

不管是直接數值模擬(DNS)，還是雷諾平均(RANS)或者大渦模擬(LES)，都因為其假定或者切入角度的問題，各自存在計算精度、計算能力、結果呈現等方面的問題，在實際的情況中研究者或者商業工程師應根據具體情況適當選用。而大渦模擬對于大尺度渦團脈動的保留，相對合理的計算能力要求以及對循環變動的分辨，使得它在未來的發動機仿真中將有極大的發展空間。本研究采用的仿真方法即為大渦模擬法。

本研究中大渦模擬的發動機是G4VDI GM 的四沖程直噴發動機，轉速1 300 r/min下，渦流控制閥完全開啟，形成一個低渦流比的工況，使得缸內流場湍流強度更強，這對于仿真模型校核更有意義。與試驗相匹配的發動機結構見圖1。基本網格大小為4 mm，在下止點處設置了109萬個網格。針對仿真具體位置和情況，對網格進行優化。在閥附近最小網格為0.25 mm，圖2示出網格優化的一些細節位置和內容。

圖1 發動機結構

圖2 網格優化設置

仿真所有的數據來源與試驗嚴格對應。試驗采用二維粒子圖像測速方法來獲得同樣工況下的流場數據，在1臺排量0.55 L、四沖程、4氣門單缸光學發動機上進行(見圖3),發動機具體運行參數見表1。

圖3 光學發動機臺架

缸徑/mm86進氣溫度/℃25活塞行程/mm94.6轉速/r·min-11 300連桿長度/mm160進氣門開啟時刻(ATDC)/(°)-366壓縮比11∶1進氣門關閉時刻(ATDC)/(°)-114進氣渦流比0.55排氣門開啟時刻(ATDC)/(°)131進氣壓力/kPa40排氣門關閉時刻(ATDC)/(°)372

1.2 本征正交分解方法

POD(本征正交分解)對于湍流的研究不同于傳統的統計方法，不僅提供了周期平均流場和脈動流場宏觀的分析，而且提供了更多內燃機循環變動的細節。POD在流場循環變動的分析中，將多循環同一時刻的流場提取為一系列基函數，而基函數對應的流場結構可以重新組合為原始循環的流場。鑒于是基于特征值的提取，POD的前幾個基函數已經包含了大多數的流場動能，即通過分析前幾個POD 模態，可以得到對于缸內流場循環變動的很多細節信息。

POD在很多領域都有著關鍵的應用，也存在許多的推導方法。本研究將結合PCA(主成分分析)的思路，理解POD計算過程在缸內流場分析中的真實含義，分析POD在缸內流場分析中的原理。

圖4左側為發動機缸內流場的速度矢量圖，右側為類比建立的一個n×n的樣本，共有n2個節點。每一個節點包含對應矢量ui,j的信息。假定uk是第k個循環的速度場，其包含了該樣本所有節點的矢量信息。

圖4 速度場矢量示意

(1)

定義總體u為循環速度場uk的集合，行數為N，N即為樣本總數，每個樣本的維度為n2。

u=u1,u2,…uNT。

(2)

PCA的核心思想是去除冗余的維度和降噪，這兩點與缸內流場分析中POD的思路不謀而合。對于總體u的分析，首先計算u1,u2,…uN的期望E(u),然后計算離散的協方差矩陣C。

(3)

為了分析方便，在計算協方差矩陣之前，一般會進行數據中心化處理，即令E(u)=0，則有：

(4)

以三維的協方差矩陣為例可以發現，協方差矩陣的對角線實際上體現了各自維度的方差，非對角線元素體現了維度間的相互影響。PCA的思想是去除冗余的維度和降噪，實際上一方面盡可能使保留的維度之間相關性小，另一方面使保留的維度方差較大。從信號與數學意義上分析，根據最大方差理論，方差的大小體現的是包含能量的大小。

(5)

因此，對該矩陣C進行對角化處理并求特征值。對應于C的最大特征值的特征向量，就是POD的主要模態。也正因如此，POD的主要模態對應于大的方差，包含了大多數的能量，并且POD的降維過程也去除了不必要的維度。本研究從PCA的思路推導出這一點，也闡釋了POD的處理過程對于缸內流場分析的實際含義。

2 缸內流場校核方法論

2.1 基于能量的缸內流場校核

將試驗的50個循環(B1～B50)和仿真的50個循環在不破壞次序的情況下組成一個新的流場矩陣，對其進行本征正交分解。因此，試驗和仿真得到共同的模態(100個)。雖然在物理意義上無法區分試驗和仿真的模態，但正是由于試驗和仿真有相同的模態，才使得模態前面的能量系數足以代表某個模態下試驗和仿真的能量大小。

(6)

(7)

COV(變異系數)在考慮平均值的基礎上體現了系數的離散程度，發動機研究中常用它來分析多循環間的變動性，比如平均有效指示壓力變異系數。

(8)

這種基于能量的校核方法在一些仿真分析中已經開始廣泛使用，并對于試驗和仿真的分析工作提供了很多有力的支持。陳豪等[10]利用本征正交分解的方法分析了缸內進氣和壓縮過程中的循環變動以及氣流的發展變化，并提及了試驗與仿真缸內流場比較的應用想法。K. Liu等[14]利用本征正交分解的方法著重分析了大渦模擬的結果，并借此討論了流場的循環差異。但該方法只考慮了主導模態的能量校核，并沒有全面地對仿真模型進行校核，因此需要增加對仿真和試驗各自的主導模態的校核。

2.2 基于模態相關性系數的缸內流場校核

在發動機缸內流場試驗和仿真結果的對比中，面對的難題是缺乏一個行之有效、單一的方法來進行衡量。若是從能量的角度來分析和比較主要模態前面的系數，缺陷在于試驗和仿真結果的差異是模態系數的差異和模態本身的共同作用。當某些循環中存在能量較大的流場結構時，該結構會嚴重影響模態，進而使得模態系數偏向試驗或者仿真，不能更好地校核仿真數據。如果能夠增加對模態的分析和比較，將大大提高模型校核的準確性。因此，對模態之間的相關系數進行分析，能夠極好地解決上述問題。具體思路是將試驗和仿真結果的差異轉化為兩個系數之間的比較，而這兩個系數是在考慮了所有主要模態的基礎上得到的。這種基于模態特征方法的比較流程見圖5。

一般而言，選擇前5個模態來進行對比，因為從能量上講，前5個模態包含了流場的大多數能量和渦團結構。在實際應用中，可以根據工況和計算結果進行調整。圖6示出了前1～5個模態能量系數的和，可以發現，前5個模態包含的能量非常接近于1，均在0.998以上。因此，本方法選取前5個模態來表征缸內流場是合理的。將試驗和仿真的所有循環進行混合POD分析，并將其結果作為衡量仿真和試驗結果的基準。與此同時，將試驗的50個循環和仿真的50個循環分別進行POD分析，計算它們與混合POD主要模態之間的相關性系數后進行比較，從而得出試驗和仿真結果的差異。

圖5 基于模態相關性系數的校核流程

圖6 本征正交分解主要模態包含的能量之和

對于模態之間相關性分析，使用以下幾種相關系數：

SI(Structure index):

(9)

MI(Magnitude index):

(10)

KER(Kinetic energy ratio):

(11)

MI*(Magnitude index*):

(12)

式中：u，v分別為試驗和仿真的主要模態的速度場矩陣。

(13)

至此，通過流程圖中系數1和系數2的對比，可以得出基于模態相關性系數的試驗仿真結果差異。

3 缸內流場校核的結果分析

本研究數據來源由SIDI發動機的大渦模擬(LES)和粒子圖像測速試驗結果兩部分組成。發動機轉速為1 300 r/min，模擬的是部分負荷工況，進氣壓力為40 kPa。仿真基于CONVERGE 2.2 動態結構湍流模型進行計算，在去除第一個循環的基礎上，共計算了50個循環的仿真結果。試驗和仿真流場為圖7所示的中心B-B平面，本研究以曲軸轉角為270°為例進行分析。

圖7 試驗和仿真流場平面示意

3.1 基于能量校核的結果分析

對流場試驗和仿真數據進行能量的校核，即為根據模態前的能量系數進行分析。對試驗和仿真數據進行有序的混合，前50個為試驗工況，后50個為仿真工況。混合后進行POD分析，得出第1模態試驗和仿真的對應系數。圖8示出第1模態兩組系數的波動情況對比。

圖8 試驗和仿真數據的第1模態系數對比

由圖8可見，仿真和試驗的模態系數均為正值，說明其對應的模態結構具有一個相同的流動方向。仿真數據的第1模態系數平均值為290.69，標準差為21.76，對應的變異系數(COV)為7.49%；試驗數據的第1模態系數平均值為253.53，標準差為29.92，對應的變異系數(COV)為11.8%。仿真的主導模態的能量偏高，但其對應的標準差和變異系數偏小，說明仿真的循環變動要小于試驗的循環變動。

從系數平均值的對比上看，仿真值相對于試驗值存在著12%左右的差異。這一方面與發動機本身性能有關，試驗數據的波動性體現了發動機本身的循環間差異；另一方面，從試驗和仿真對比的角度來講，仿真結果需要通過調整參數等方法來確保仿真結果的模態系數和試驗結果更為接近，進而增強仿真結果的說服力。

圖9示出模態2，3的試驗和仿真系數對比。從物理意義上講，除了包含大多數能量的第1模態，其他模態在一定程度上體現了較小的渦團結構特征。為了更直觀地說明，圖10示出了試驗和仿真結果中某個循環的實際流場圖。不難發現，除了左側流場，尤其是左上側體現缸內流場進氣流態之外，其他的區域存在著較小的渦團結構，而這些小渦團結構又存在著極大的差別。秦文瑾等[15]指出，模態2顯示的流場結構整體性明顯下降．而局部小渦團則明顯增多，模態3中該現象更加明顯。說明模態越高，捕獲局部隨機脈動小渦團的能力越強。圖10中的圓圈部分的次級渦團結構是模態2，3的主要特征點。從模態系數趨勢和區別上來看，模態2的能量系數在量級上是很接近的，但小渦團模態的方向卻是相反的，這一點與很多因素有關，如對第1模態結果的補償，對圓圈部分次級渦團結構的捕獲等。模態3的系數對于仿真和試驗均具有隨機性，體現了高階小渦團的無序性。

圖9 試驗和仿真數據的第2，3模態系數對比

圖10 試驗和仿真單循環的速度場矢量圖

3.2 基于模態相關性校核的結果分析

從模態間相關系數的角度來分析，分別計算了試驗50個循環，仿真50個循環，混合100個循環這3種情況的POD結果。取試驗和仿真50個循環的平均流場，繪制流場速度場(見圖11a和圖11b),展示試驗和仿真流場的形態差異。為了進一步量化這種差異，依照基于模態相關性校核的思路和方法，對試驗、仿真、混合的POD主要模態進行分析。

以混合POD的100個模態為基準，選取前5個模態，分別與仿真和試驗的前5個模態進行相關性分析，進而比較其相關性系數。圖11c,11d,11e分別為試驗、仿真、混合三者的第1模態差異，用以形象化理解該方法結果分析的含義。混合流場的第1模態在流場的左上和左下有明顯的大速度區域，而試驗的第1模態側重體現左下的大速度區域，仿真的第1模態側重體現左上的大速度區域，圖11c,11d,11e的形態也驗證了圖8中仿真的第1模態系數平均值要大于試驗的第1模態系數平均值。這種方法對于前5個模態的相關性分析，實質上是在討論相比等權重的混合POD，試驗和仿真單獨的POD流場模態所占的權重。理想情況下，這兩個相關性系數應該都接近1或者其他臨界數(視不同相關性系數定義)，而實際情況下，兩個系數與1的差異以及兩個系數間的差異體現了試驗和仿真流場的差異。

圖11 試驗和仿真平均流場以及POD第1模態流場矢量

進一步對試驗、仿真、混合的主要模態分別進行相關性計算，得出前面在模態校核流程中提到的系數1和系數2(見表2)。SI,MI的系數差異體現了試驗和仿真的差異,由于本征正交分解過程中，模態只保留了結構信息，不包含能量信息，因此，KER的值為1,MI*與1的差值在10-6左右。通過這種全局分析的方法能夠得到流場主要模態的差異，進而展現了流場結構的差異。這種方法實質上是在等權重的混合POD下，試驗和仿真單獨的POD流場所占的權重。系數1，2的理論值是1，并且接近1的程度和兩個系數之間的差異體現了試驗和仿真結果的差異。因此，可以基于某一個相關性系數(如SI)建立例如

(14)

這樣的基準數，通過極其大量的試驗，去挖掘得出普適意義上的臨界值。其中，Indexexp，Indexsim為試驗和仿真的某一種相關性系數，用以構建能明確化仿真模型適用范圍的基準數。在基準數的某一個范圍內，認為試驗和仿真結果接近，可接受；在另外一個基準數的范圍，則認定試驗和仿真結果匹配精確度高。這樣的工作是可行的并且對于試驗和仿真差異的界定具有極大的意義。

表2 模態相關性系數(系數1和系數2)

4 結束語

結合PCA(主成分分析)信號分析的角度，深入推導并剖析了本征正交分解在處理缸內流場過程中的物理意義，闡述了發動機缸內流場的主要模態實際上是原始流場經過“去除冗余的維度”和“降噪”的結果，對應于協方差矩陣對角線較大的方差值，即能量。提出并分析了能量和模態相關性兩個方面的校核方法，其中創新性地提出了從模態相關性角度的計算方法和思路，對于過去僅從能量方面分析主要模態的系數，是一種合理的補充，有利于大渦模擬結果的全面校核。

基于高速粒子圖像測速和大渦模擬數值計算的結果，對發動機缸內流場進行了模態和能量角度的校核分析。一方面，分析了發動機缸內流場大渦模擬和試驗結果的主要模態的系數波動，模態1的系數差異體現了包含大部分能量的流場結構上的試驗和仿真的差異，模態2，3的系數差異體現了局部小渦團的特征差異。另一方面，結合缸內流場速度矢量圖、主要模態的相關性系數進行了定性和定量的分析。相關性系數體現的是等權重的混合POD下，試驗和仿真單獨的POD流場模態所占的權重。兩個相關性系數的值的大小和差值，可以有效地體現仿真校核結果。

在模態相關性系數的校核方面，提出了進一步深化的方向，即試驗，仿真模態的相關性系數1，2的大小本身以及差值，體現了試驗仿真結果的客觀差異，通過構建準則數，并基于大量實踐，可以得出判斷仿真結果優劣的準則數范圍。這對于目前仿真結果可信度難以量化的情況，是一種可以考慮的實踐思路。

參考文獻：

[1] 蔣堅,高希彥.汽油缸內直噴式技術的研究與應用[J].內燃機工程,2003(5):39-44，58.

[2] 李相超,張玉銀,許敏,等.直噴汽油機缸內噴霧濕壁問題研究[J].內燃機工程,2012,33(5):17-23.

[3] 田少雄,孔令遜,許敏,等.缸內渦流比對冷起動燃燒火焰的影響探究[J].內燃機工程,2017,38(3):1-8.

[4] 王天友,劉大明,張學恩,等.可變氣門升程下汽油機缸內氣體流動特性的研究[J].內燃機學報,2008,26(5):420-428.

[5] Sick V,Chen H,Abraham P S,et al.Proper-orthogonal decomposition analysis for engine research[C]//9th Congress,Gasoline Direct Injection Engines.Essen：[s.n.],2012:1-12.

[6] Lumley J L.The structure of inhomogeneous turbulent flows[J].Atmospheric turbulence and radio wave propagation,1967:166-178.

[7] Sirovich L.Turbulence and the dynamics of coherent structures.I.Coherent structures[J].Quarterly of applied mathematics,1987,45(3):561-571.

[8] 秦文瑾,許敏,孔令遜.直噴式汽油機缸內渦量場的本征正交分解[J].內燃機學報,2016,34(3):246-252.

[9] Fogleman M A.Low-dimensional models of internal combustion engine flows using the proper orthogonal decomposition[C]//Dissertation Abstracts International.[S.l.]：John L.Lumley,2005.

[10] 陳豪.直噴汽油機缸內過程穩定性機理的可視化研究[D].上海：上海交通大學,2014.

[11] Druault P,Guibert P,Alizon F.Use of proper orthogonal decomposition for time interpolation from PIV data[J].Experiments in Fluids,2005,39(6):1009-1023.

[12] Borée J.Extended proper orthogonal decomposition: a tool to analyse correlated events in turbulent flows[J].Experiments in Fluids,2003,35(2):188-192.

[13] Germano M,Piomelli U,Moin P,et al.A dynamic subgrid-scale eddy viscosity model[J].Physics of Fluids A:Fluid Dynamics,1991,3(7):1760-1765.

[14] Liu K,Haworth D C.Development and assessment of POD for analysis of turbulent flow in piston engines[C].SAE Paper 2011-01-0830.

[15] 秦文瑾.內燃機缸內湍流特性及其循環變動的大渦模擬與本征正交分解研究[D].大連：大連理工大學,2014.