顧鵬,葛鵬輝,許敏
(上海交通大學汽車工程研究院,上海 200240)
活塞點燃直噴式(Spark Ignition Direct Injection)汽油機相比進氣道噴射汽油機在燃油經濟性和顆粒物排放方面具有很大的優(yōu)勢,并得到了廣泛的研究和應用[1-2]。研究表明,缸內流場對發(fā)動機燃油霧化、火焰?zhèn)鞑ヂ窂胶退俾示哂兄匾挠绊慬3],而缸內流場結構又受到諸多因素如進氣道形狀、進氣歧管傾角和缸蓋結構等的影響[4]。如果能夠開發(fā)出較為準確的發(fā)動機三維模型,這將極大地節(jié)約發(fā)動機開發(fā)成本。因此,發(fā)動機的三維模型校核也顯得尤為重要。
與傳統(tǒng)的RANS仿真不同的是,大渦模擬和試驗的結果都存在著強烈的循環(huán)差異,對發(fā)動機的三維模型校核帶來了極大的障礙。而本征正交分解(Proper orthogonal decomposition)通過基函數(shù)(mode)為大渦模擬和試驗結果提供了一個比較的基礎[5]。Lumley 等[6]最早將本征正交分解的方法引入到湍流場的分析中,將湍流場中能量較大的大尺度渦團結構分離提取出來。由于發(fā)動機缸內流場同樣具有高度湍流特性,因此該方法逐漸被用于發(fā)動機的流場分析中。Sirovich[7]提出的“snapshot”的本征正交分解方法,在與原先的直接法等效的情況下,大大提升了計算效率,也廣泛應用于流場分析。在發(fā)動機流場試驗和仿真分析方面,很多研究者都基于不同的試驗和仿真工況做了很多有益的工作。秦文瑾、許敏等[8]證明了本征正交分解可以作為分析發(fā)動機缸內流場模擬結果準確性的有效工具。Fogleman等[9]通過相位固定的本征正交分解的方法,對不同曲軸轉角的缸內流場變化進行了分析。陳豪等[10]指出目前對實際發(fā)動機流場結果本征正交分解的研究相對匱乏,并闡述了試驗與仿真缸內流場比較的應用思路。此外,一些研究者對試驗數(shù)據的完善,對缸內流場仿真校核的工作也具有重要意義。Druault等[11]通過本征正交分解對PIV的流場結果進行時間跨度上的插值,Borée 等[12]通過本征正交分解的方法闡述了時間和空間位置上的關聯(lián)。
本研究在前人工作的基礎上,著眼于基于本征正交分解的缸內流場仿真校核方法,從能量和模態(tài)的角度,全面合理地評價試驗和仿真缸內流場結果的差異性,借此提出并建立缸內流場仿真校核的流程策略。
關于缸內流場循環(huán)的分析,大量基礎性的研究是基于雷諾分解,建立對湍流脈動量的分析。雷諾方程的核心價值在于,將對湍流的研究轉化為了建立可以求解的、與雷諾應力這個變量相關的輸運方程。1887年提出的Boussinesq假設,通過與分子運動動量傳遞的類比,將雷諾應力轉化為渦黏性系數(shù),而渦黏性系數(shù)進一步地被分解為對湍流渦團的長度和速度尺度這兩個變量的求解。在后續(xù)的缸內流場研究中區(qū)分零方程模型、單方程模型、多方程模型等,一定意義上就是在于求解這兩個變量的微分方程的個數(shù)。Boussinesq假設之后,混合長度理論卡門相似理論等將渦團的長度和速度尺度以經驗或者代數(shù)的方法表達,在早期的缸內流場分析中起到了很大的推動作用,但是它在時間、空間變化上的表征仍存在許多不足。因此,從1942年Kolmogorov和1945年Prandtl的研究開始,人們逐漸開始建立κ方程、渦黏度方程(Spalart-Allmaras)等單方程模型?,F(xiàn)在主流的CFD軟件中,頻繁使用κ-ε,κ-ω等雙方程模型,對于這些模型,研究者們還在不斷地根據具體的使用情況進行修正,例如應用較廣的壓縮性修正、強旋流修正等。此外,在強調各向異性的缸內流場分析中,對雷諾應力直接進行分析的RSM模型也應用極廣。
近幾十年來,研究者們逐漸認識到缸內流場可以分解為一系列的有序流場結構,現(xiàn)在商業(yè)應用極廣的LES(大渦模擬)實際上就是基于空間尺度的劃分,將對應于不同尺度的渦團劃分為可解尺度和亞網格尺度,兩者分別使用直接求解的方法和假設模型的方法進行分析。大渦模擬的核心前提是大尺度渦團與仿真時設置的初始和邊界條件相關性較大,而小尺度渦團受它們的影響較小,大尺度渦團的各向異性和小尺度渦團的各向同性促成了大渦模擬將它們分開處理的合理性。
大渦模擬的思路如下:首先消除小渦的脈動影響,得出局部的平均場,即為大尺度場,而亞網格尺度場則為原始流場與該平均場的差值。大尺度場根據網格,由N-S方程直接求解,而對于亞網格尺度場,常用的模型有Smagorinsky模型、WALE模型、亞網格單方程模型以及動態(tài)結構模型。其中,對于Smagorinsky模型和亞網格單方程模型,在后續(xù)LES的研究和計算中模型系數(shù)進行了修正,以滿足更多的情況,進而產生了動態(tài)Smagorinsky模型[13]和亞網格單方程模型。
不管是直接數(shù)值模擬(DNS),還是雷諾平均(RANS)或者大渦模擬(LES),都因為其假定或者切入角度的問題,各自存在計算精度、計算能力、結果呈現(xiàn)等方面的問題,在實際的情況中研究者或者商業(yè)工程師應根據具體情況適當選用。而大渦模擬對于大尺度渦團脈動的保留,相對合理的計算能力要求以及對循環(huán)變動的分辨,使得它在未來的發(fā)動機仿真中將有極大的發(fā)展空間。本研究采用的仿真方法即為大渦模擬法。
本研究中大渦模擬的發(fā)動機是G4VDI GM 的四沖程直噴發(fā)動機,轉速1 300 r/min下,渦流控制閥完全開啟,形成一個低渦流比的工況,使得缸內流場湍流強度更強,這對于仿真模型校核更有意義。與試驗相匹配的發(fā)動機結構見圖1。基本網格大小為4 mm,在下止點處設置了109萬個網格。針對仿真具體位置和情況,對網格進行優(yōu)化。在閥附近最小網格為0.25 mm,圖2示出網格優(yōu)化的一些細節(jié)位置和內容。
圖1 發(fā)動機結構
圖2 網格優(yōu)化設置
仿真所有的數(shù)據來源與試驗嚴格對應。試驗采用二維粒子圖像測速方法來獲得同樣工況下的流場數(shù)據,在1臺排量0.55 L、四沖程、4氣門單缸光學發(fā)動機上進行(見圖3),發(fā)動機具體運行參數(shù)見表1。
圖3 光學發(fā)動機臺架
缸徑/mm86進氣溫度/℃25活塞行程/mm94.6轉速/r·min-11 300連桿長度/mm160進氣門開啟時刻(ATDC)/(°)-366壓縮比11∶1進氣門關閉時刻(ATDC)/(°)-114進氣渦流比0.55排氣門開啟時刻(ATDC)/(°)131進氣壓力/kPa40排氣門關閉時刻(ATDC)/(°)372
POD(本征正交分解)對于湍流的研究不同于傳統(tǒng)的統(tǒng)計方法,不僅提供了周期平均流場和脈動流場宏觀的分析,而且提供了更多內燃機循環(huán)變動的細節(jié)。POD在流場循環(huán)變動的分析中,將多循環(huán)同一時刻的流場提取為一系列基函數(shù),而基函數(shù)對應的流場結構可以重新組合為原始循環(huán)的流場。鑒于是基于特征值的提取,POD的前幾個基函數(shù)已經包含了大多數(shù)的流場動能,即通過分析前幾個POD 模態(tài),可以得到對于缸內流場循環(huán)變動的很多細節(jié)信息。
POD在很多領域都有著關鍵的應用,也存在許多的推導方法。本研究將結合PCA(主成分分析)的思路,理解POD計算過程在缸內流場分析中的真實含義,分析POD在缸內流場分析中的原理。
圖4左側為發(fā)動機缸內流場的速度矢量圖,右側為類比建立的一個n×n的樣本,共有n2個節(jié)點。每一個節(jié)點包含對應矢量ui,j的信息。假定uk是第k個循環(huán)的速度場,其包含了該樣本所有節(jié)點的矢量信息。
圖4 速度場矢量示意
(1)
定義總體u為循環(huán)速度場uk的集合,行數(shù)為N,N即為樣本總數(shù),每個樣本的維度為n2。
u=u1,u2,…uNT。
(2)
PCA的核心思想是去除冗余的維度和降噪,這兩點與缸內流場分析中POD的思路不謀而合。對于總體u的分析,首先計算u1,u2,…uN的期望E(u),然后計算離散的協(xié)方差矩陣C。
(3)
為了分析方便,在計算協(xié)方差矩陣之前,一般會進行數(shù)據中心化處理,即令E(u)=0,則有:
(4)
以三維的協(xié)方差矩陣為例可以發(fā)現(xiàn),協(xié)方差矩陣的對角線實際上體現(xiàn)了各自維度的方差,非對角線元素體現(xiàn)了維度間的相互影響。PCA的思想是去除冗余的維度和降噪,實際上一方面盡可能使保留的維度之間相關性小,另一方面使保留的維度方差較大。從信號與數(shù)學意義上分析,根據最大方差理論,方差的大小體現(xiàn)的是包含能量的大小。
(5)
因此,對該矩陣C進行對角化處理并求特征值。對應于C的最大特征值的特征向量,就是POD的主要模態(tài)。也正因如此,POD的主要模態(tài)對應于大的方差,包含了大多數(shù)的能量,并且POD的降維過程也去除了不必要的維度。本研究從PCA的思路推導出這一點,也闡釋了POD的處理過程對于缸內流場分析的實際含義。
將試驗的50個循環(huán)(B1~B50)和仿真的50個循環(huán)在不破壞次序的情況下組成一個新的流場矩陣,對其進行本征正交分解。因此,試驗和仿真得到共同的模態(tài)(100個)。雖然在物理意義上無法區(qū)分試驗和仿真的模態(tài),但正是由于試驗和仿真有相同的模態(tài),才使得模態(tài)前面的能量系數(shù)足以代表某個模態(tài)下試驗和仿真的能量大小。
(6)
(7)
COV(變異系數(shù))在考慮平均值的基礎上體現(xiàn)了系數(shù)的離散程度,發(fā)動機研究中常用它來分析多循環(huán)間的變動性,比如平均有效指示壓力變異系數(shù)。
(8)
這種基于能量的校核方法在一些仿真分析中已經開始廣泛使用,并對于試驗和仿真的分析工作提供了很多有力的支持。陳豪等[10]利用本征正交分解的方法分析了缸內進氣和壓縮過程中的循環(huán)變動以及氣流的發(fā)展變化,并提及了試驗與仿真缸內流場比較的應用想法。K. Liu等[14]利用本征正交分解的方法著重分析了大渦模擬的結果,并借此討論了流場的循環(huán)差異。但該方法只考慮了主導模態(tài)的能量校核,并沒有全面地對仿真模型進行校核,因此需要增加對仿真和試驗各自的主導模態(tài)的校核。
在發(fā)動機缸內流場試驗和仿真結果的對比中,面對的難題是缺乏一個行之有效、單一的方法來進行衡量。若是從能量的角度來分析和比較主要模態(tài)前面的系數(shù),缺陷在于試驗和仿真結果的差異是模態(tài)系數(shù)的差異和模態(tài)本身的共同作用。當某些循環(huán)中存在能量較大的流場結構時,該結構會嚴重影響模態(tài),進而使得模態(tài)系數(shù)偏向試驗或者仿真,不能更好地校核仿真數(shù)據。如果能夠增加對模態(tài)的分析和比較,將大大提高模型校核的準確性。因此,對模態(tài)之間的相關系數(shù)進行分析,能夠極好地解決上述問題。具體思路是將試驗和仿真結果的差異轉化為兩個系數(shù)之間的比較,而這兩個系數(shù)是在考慮了所有主要模態(tài)的基礎上得到的。這種基于模態(tài)特征方法的比較流程見圖5。
一般而言,選擇前5個模態(tài)來進行對比,因為從能量上講,前5個模態(tài)包含了流場的大多數(shù)能量和渦團結構。在實際應用中,可以根據工況和計算結果進行調整。圖6示出了前1~5個模態(tài)能量系數(shù)的和,可以發(fā)現(xiàn),前5個模態(tài)包含的能量非常接近于1,均在0.998以上。因此,本方法選取前5個模態(tài)來表征缸內流場是合理的。將試驗和仿真的所有循環(huán)進行混合POD分析,并將其結果作為衡量仿真和試驗結果的基準。與此同時,將試驗的50個循環(huán)和仿真的50個循環(huán)分別進行POD分析,計算它們與混合POD主要模態(tài)之間的相關性系數(shù)后進行比較,從而得出試驗和仿真結果的差異。
圖5 基于模態(tài)相關性系數(shù)的校核流程
圖6 本征正交分解主要模態(tài)包含的能量之和
對于模態(tài)之間相關性分析,使用以下幾種相關系數(shù):
SI(Structure index):
(9)
MI(Magnitude index):
(10)
KER(Kinetic energy ratio):
(11)
MI*(Magnitude index*):
(12)
式中:u,v分別為試驗和仿真的主要模態(tài)的速度場矩陣。
(13)
至此,通過流程圖中系數(shù)1和系數(shù)2的對比,可以得出基于模態(tài)相關性系數(shù)的試驗仿真結果差異。
本研究數(shù)據來源由SIDI發(fā)動機的大渦模擬(LES)和粒子圖像測速試驗結果兩部分組成。發(fā)動機轉速為1 300 r/min,模擬的是部分負荷工況,進氣壓力為40 kPa。仿真基于CONVERGE 2.2 動態(tài)結構湍流模型進行計算,在去除第一個循環(huán)的基礎上,共計算了50個循環(huán)的仿真結果。試驗和仿真流場為圖7所示的中心B-B平面,本研究以曲軸轉角為270°為例進行分析。
圖7 試驗和仿真流場平面示意
對流場試驗和仿真數(shù)據進行能量的校核,即為根據模態(tài)前的能量系數(shù)進行分析。對試驗和仿真數(shù)據進行有序的混合,前50個為試驗工況,后50個為仿真工況?;旌虾筮M行POD分析,得出第1模態(tài)試驗和仿真的對應系數(shù)。圖8示出第1模態(tài)兩組系數(shù)的波動情況對比。
圖8 試驗和仿真數(shù)據的第1模態(tài)系數(shù)對比
由圖8可見,仿真和試驗的模態(tài)系數(shù)均為正值,說明其對應的模態(tài)結構具有一個相同的流動方向。仿真數(shù)據的第1模態(tài)系數(shù)平均值為290.69,標準差為21.76,對應的變異系數(shù)(COV)為7.49%;試驗數(shù)據的第1模態(tài)系數(shù)平均值為253.53,標準差為29.92,對應的變異系數(shù)(COV)為11.8%。仿真的主導模態(tài)的能量偏高,但其對應的標準差和變異系數(shù)偏小,說明仿真的循環(huán)變動要小于試驗的循環(huán)變動。
從系數(shù)平均值的對比上看,仿真值相對于試驗值存在著12%左右的差異。這一方面與發(fā)動機本身性能有關,試驗數(shù)據的波動性體現(xiàn)了發(fā)動機本身的循環(huán)間差異;另一方面,從試驗和仿真對比的角度來講,仿真結果需要通過調整參數(shù)等方法來確保仿真結果的模態(tài)系數(shù)和試驗結果更為接近,進而增強仿真結果的說服力。
圖9示出模態(tài)2,3的試驗和仿真系數(shù)對比。從物理意義上講,除了包含大多數(shù)能量的第1模態(tài),其他模態(tài)在一定程度上體現(xiàn)了較小的渦團結構特征。為了更直觀地說明,圖10示出了試驗和仿真結果中某個循環(huán)的實際流場圖。不難發(fā)現(xiàn),除了左側流場,尤其是左上側體現(xiàn)缸內流場進氣流態(tài)之外,其他的區(qū)域存在著較小的渦團結構,而這些小渦團結構又存在著極大的差別。秦文瑾等[15]指出,模態(tài)2顯示的流場結構整體性明顯下降.而局部小渦團則明顯增多,模態(tài)3中該現(xiàn)象更加明顯。說明模態(tài)越高,捕獲局部隨機脈動小渦團的能力越強。圖10中的圓圈部分的次級渦團結構是模態(tài)2,3的主要特征點。從模態(tài)系數(shù)趨勢和區(qū)別上來看,模態(tài)2的能量系數(shù)在量級上是很接近的,但小渦團模態(tài)的方向卻是相反的,這一點與很多因素有關,如對第1模態(tài)結果的補償,對圓圈部分次級渦團結構的捕獲等。模態(tài)3的系數(shù)對于仿真和試驗均具有隨機性,體現(xiàn)了高階小渦團的無序性。
圖9 試驗和仿真數(shù)據的第2,3模態(tài)系數(shù)對比
圖10 試驗和仿真單循環(huán)的速度場矢量圖
從模態(tài)間相關系數(shù)的角度來分析,分別計算了試驗50個循環(huán),仿真50個循環(huán),混合100個循環(huán)這3種情況的POD結果。取試驗和仿真50個循環(huán)的平均流場,繪制流場速度場(見圖11a和圖11b),展示試驗和仿真流場的形態(tài)差異。為了進一步量化這種差異,依照基于模態(tài)相關性校核的思路和方法,對試驗、仿真、混合的POD主要模態(tài)進行分析。
以混合POD的100個模態(tài)為基準,選取前5個模態(tài),分別與仿真和試驗的前5個模態(tài)進行相關性分析,進而比較其相關性系數(shù)。圖11c,11d,11e分別為試驗、仿真、混合三者的第1模態(tài)差異,用以形象化理解該方法結果分析的含義?;旌狭鲌龅牡?模態(tài)在流場的左上和左下有明顯的大速度區(qū)域,而試驗的第1模態(tài)側重體現(xiàn)左下的大速度區(qū)域,仿真的第1模態(tài)側重體現(xiàn)左上的大速度區(qū)域,圖11c,11d,11e的形態(tài)也驗證了圖8中仿真的第1模態(tài)系數(shù)平均值要大于試驗的第1模態(tài)系數(shù)平均值。這種方法對于前5個模態(tài)的相關性分析,實質上是在討論相比等權重的混合POD,試驗和仿真單獨的POD流場模態(tài)所占的權重。理想情況下,這兩個相關性系數(shù)應該都接近1或者其他臨界數(shù)(視不同相關性系數(shù)定義),而實際情況下,兩個系數(shù)與1的差異以及兩個系數(shù)間的差異體現(xiàn)了試驗和仿真流場的差異。
圖11 試驗和仿真平均流場以及POD第1模態(tài)流場矢量
進一步對試驗、仿真、混合的主要模態(tài)分別進行相關性計算,得出前面在模態(tài)校核流程中提到的系數(shù)1和系數(shù)2(見表2)。SI,MI的系數(shù)差異體現(xiàn)了試驗和仿真的差異,由于本征正交分解過程中,模態(tài)只保留了結構信息,不包含能量信息,因此,KER的值為1,MI*與1的差值在10-6左右。通過這種全局分析的方法能夠得到流場主要模態(tài)的差異,進而展現(xiàn)了流場結構的差異。這種方法實質上是在等權重的混合POD下,試驗和仿真單獨的POD流場所占的權重。系數(shù)1,2的理論值是1,并且接近1的程度和兩個系數(shù)之間的差異體現(xiàn)了試驗和仿真結果的差異。因此,可以基于某一個相關性系數(shù)(如SI)建立例如
(14)
這樣的基準數(shù),通過極其大量的試驗,去挖掘得出普適意義上的臨界值。其中,Indexexp,Indexsim為試驗和仿真的某一種相關性系數(shù),用以構建能明確化仿真模型適用范圍的基準數(shù)。在基準數(shù)的某一個范圍內,認為試驗和仿真結果接近,可接受;在另外一個基準數(shù)的范圍,則認定試驗和仿真結果匹配精確度高。這樣的工作是可行的并且對于試驗和仿真差異的界定具有極大的意義。
表2 模態(tài)相關性系數(shù)(系數(shù)1和系數(shù)2)
結合PCA(主成分分析)信號分析的角度,深入推導并剖析了本征正交分解在處理缸內流場過程中的物理意義,闡述了發(fā)動機缸內流場的主要模態(tài)實際上是原始流場經過“去除冗余的維度”和“降噪”的結果,對應于協(xié)方差矩陣對角線較大的方差值,即能量。提出并分析了能量和模態(tài)相關性兩個方面的校核方法,其中創(chuàng)新性地提出了從模態(tài)相關性角度的計算方法和思路,對于過去僅從能量方面分析主要模態(tài)的系數(shù),是一種合理的補充,有利于大渦模擬結果的全面校核。
基于高速粒子圖像測速和大渦模擬數(shù)值計算的結果,對發(fā)動機缸內流場進行了模態(tài)和能量角度的校核分析。一方面,分析了發(fā)動機缸內流場大渦模擬和試驗結果的主要模態(tài)的系數(shù)波動,模態(tài)1的系數(shù)差異體現(xiàn)了包含大部分能量的流場結構上的試驗和仿真的差異,模態(tài)2,3的系數(shù)差異體現(xiàn)了局部小渦團的特征差異。另一方面,結合缸內流場速度矢量圖、主要模態(tài)的相關性系數(shù)進行了定性和定量的分析。相關性系數(shù)體現(xiàn)的是等權重的混合POD下,試驗和仿真單獨的POD流場模態(tài)所占的權重。兩個相關性系數(shù)的值的大小和差值,可以有效地體現(xiàn)仿真校核結果。
在模態(tài)相關性系數(shù)的校核方面,提出了進一步深化的方向,即試驗,仿真模態(tài)的相關性系數(shù)1,2的大小本身以及差值,體現(xiàn)了試驗仿真結果的客觀差異,通過構建準則數(shù),并基于大量實踐,可以得出判斷仿真結果優(yōu)劣的準則數(shù)范圍。這對于目前仿真結果可信度難以量化的情況,是一種可以考慮的實踐思路。
參考文獻:
[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] 陳豪.直噴汽油機缸內過程穩(wěn)定性機理的可視化研究[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] 秦文瑾.內燃機缸內湍流特性及其循環(huán)變動的大渦模擬與本征正交分解研究[D].大連:大連理工大學,2014.