內(nèi)燃機進排氣變網(wǎng)格一維非定常流動計算方法

張琨， 陳昊天， 鄧康耀， 胡志龍， 錢躍華， 劉博

(1.上海交通大學(xué) 動力機械及工程教育部重點實驗室， 上海 200240； 2.上海船舶設(shè)備研究所， 上海 200031； 3.中國船舶集團有限公司第七一一研究所， 上海 201108； 4.中船動力研究院有限公司， 上海 200120)

近年來，隨著內(nèi)燃機強化程度的提高，進排氣的流速更高和能量更大，特別是排氣壓力波變化梯度加大，造成壓力波在進排氣管中的傳播過程更加復(fù)雜[1]。為了更加準(zhǔn)確地預(yù)測進排氣系統(tǒng)的能量損失和轉(zhuǎn)化過程，進排氣系統(tǒng)仿真需有更高精度的計算模型[2]。但是，高精度的計算格式通常需要犧牲計算效率，因此內(nèi)燃機性能仿真在不降低計算效率的基礎(chǔ)上提高進排氣流動模型精度至關(guān)重要。

Benson等[3-4]開展了內(nèi)燃機工作過程仿真，引入特征線方法為進排氣管道一維非定常流動計算方法，該方法較簡便，經(jīng)改進后能達到較高的精度，但流量誤差較大，不能計算帶有EGR多種混合成分的流動計算，因此在隨后的應(yīng)用中逐漸被淘汰[5]。隨著對稱或迎風(fēng)激波捕捉數(shù)值方法用于求解雙曲守恒方程組的興起，其逐漸取代特征線法成為一維氣體動力學(xué)計算的主流格式[6]，并逐步發(fā)展出了Lax-Friedrichs格式、Lax-Wendroff格式、MacCormack格式、ROE格式、TVD格式、ENO格式[7]等多種高精度計算格式。上述迭代格式在內(nèi)燃機進排氣流動仿真過程中，通常網(wǎng)格都是均勻網(wǎng)格，提高計算精度需要增加網(wǎng)格密度或者采用高精度的求解格式，這將導(dǎo)致計算效率降低。同時，實際內(nèi)燃機進排氣過程中，不同區(qū)域不同時刻對于網(wǎng)格分辨率的需求是不一樣的，例如在排氣壓力波峰面經(jīng)過的區(qū)域，需要采用較高分辨率捕捉壓力從而提高仿真精度，而在其他區(qū)域可以采用較粗的網(wǎng)格。因此，傳統(tǒng)迭代均勻網(wǎng)格的使用難以實現(xiàn)計算精度與計算效率的平衡。

為了克服傳統(tǒng)內(nèi)燃機進排氣一維非定常流動不能實時調(diào)整網(wǎng)格密度的缺點，本文在均勻網(wǎng)格TVD格式基礎(chǔ)上，考慮空間位置對通量限制器中的Sweby因子的影響，推導(dǎo)建立適用于內(nèi)燃機進排氣流動計算的非均勻網(wǎng)格TVD迭代格式。而后通過捕捉管道中參數(shù)梯度變化較大的接觸面位置，對該區(qū)域?qū)嵤┚W(wǎng)格加密，從而實現(xiàn)網(wǎng)格自適應(yīng)加密和合并，建立可變網(wǎng)格一維非定常流動數(shù)值仿真計算方法。最后通過建立激波管和內(nèi)燃機整機算例對比變網(wǎng)格與傳統(tǒng)計算方法的計算精度和計算效率。

1 非均勻網(wǎng)格一維非定常流動模型

1.1 基本控制方程及網(wǎng)格劃分

考慮截面變化、有摩擦、傳熱的一維非定常流動控制方程為[1]：

Wt+F(W)x=S(x,W)

(1)

其中：

一維非定常流動方程組的Jacobi矩陣為[7]：

(2)

Jacobi矩陣J(W)由特征值組成的對角矩陣為Λ，右特征向量矩陣為P，左特征向量矩陣為Q。

圖1 有限體積法非均勻網(wǎng)格劃分示意Fig.1 Schematic diagram of finite volume method for heterogeneous meshing

1.2 源項處理

對于帶有源項的一維非定常流動控制方程組式(1)，在從n時刻推進至n+1時刻時，可將其分裂為2步求解[8]：

Wt+F(W)x=0

(3)

(4)

為了得到二階精度的解，采用二階Runge-Kutta方法：

K1=ΔtS(tn,Wn)

(5)

K2=ΔtS(tn+Δt,Wn+K1)

(6)

(7)

1.3 非均勻網(wǎng)格TVD格式

對于不含源項的一維非定常流動控制方程組(3)，因傳統(tǒng)的TVD格式是針對標(biāo)量方程的，因此首先需要對控制方程組運用特征理論解耦。

定義新的特征變量:

U=QW

(8)

將方程組(3)左側(cè)乘以Q，可以得到:

(QW)t+QJP(QW)x=0

(9)

結(jié)合式(8)得:

Ut+ΛUx=0

(10)

定義新的變量:

Fu=QF

(11)

根據(jù)Harten推導(dǎo)的證明，對于迎風(fēng)迭代格式:

(12)

(13)

式中：h(φ(rj+1/2))為對角矩陣；r為Sweby因子，其含義為網(wǎng)格邊界的上游梯度和下游梯度的比值。式(12)和(13)為TVD格式。在滿足穩(wěn)定性準(zhǔn)則前提下，如果通量限制器φ(r)滿足:

(14)

對于對角矩陣的第k行元素，均勻網(wǎng)格的Sweby因子定義[9]為:

(15)

圖2為均勻網(wǎng)格與非均勻網(wǎng)格的通量求解示意圖，網(wǎng)格間距離對參數(shù)的梯度變化率計算有較大影響，非均勻網(wǎng)格的網(wǎng)格間距離不同，因此非均勻網(wǎng)格和均勻網(wǎng)格的梯度變化率求解方法是不一樣的。

圖2 不同網(wǎng)格劃分通量求解示意Fig.2 Diagram of flux solution for different mesh division

對于非均勻網(wǎng)格，需考慮空間位置對Sweby因子的影響，因此，非均勻網(wǎng)格的Sweby因子定義為:

(16)

按照上述定義，原TVD格式定義中的式(13)也需要相應(yīng)修正為:

(17)

當(dāng)采用均勻網(wǎng)格時，非均勻網(wǎng)格計算式(16)和式(17)退化為均勻網(wǎng)格計算式(15)和(13)。

將式(12)和(17)左乘矩陣P可以獲得原始變量的迭代方程，即：

(18)

(19)

式(18)和(19)為考慮管道面積變化的一維非定常流動TVD格式，公式中單元網(wǎng)格內(nèi)的參數(shù)用平均值表示，下標(biāo)為j+1/2或j-1/2的交界面參數(shù)根據(jù)單元網(wǎng)格中的參數(shù)通過Roe平均獲得。

2 管道變網(wǎng)格數(shù)值仿真方法

一維非定常流動數(shù)值仿真過程中為了保證計算穩(wěn)定性，迭代步長需要滿足CFL準(zhǔn)則：

(20)

從式(20)可以看出，管道流動的迭代步長是由網(wǎng)格尺度Δxj、網(wǎng)格速度|uj|和音速aj決定的。分母|uj|+aj是根據(jù)流動狀態(tài)迭代計算出的，其數(shù)值大小不可控，而分子網(wǎng)格尺度可以調(diào)整，通過調(diào)整分子網(wǎng)格尺度Δxj的大小，可以實現(xiàn)迭代步長的調(diào)整。一般情況下，管道數(shù)值仿真的網(wǎng)格尺度是定值，設(shè)置較大的網(wǎng)格尺度Δxj可以增加迭代步長，減少計算所需時間。但是，當(dāng)沿管長方向管道參數(shù)變化較劇烈時，仿真結(jié)果不能準(zhǔn)確反映各網(wǎng)格點的參數(shù)變化情況，仿真結(jié)果存在失真現(xiàn)象。如果設(shè)置較小的網(wǎng)格尺度，則會增加數(shù)據(jù)存儲和計算時間，大大降低計算效率。通過接觸斷面區(qū)域采用局部加密網(wǎng)格，可以有效提高仿真精度和效率。但在內(nèi)燃機實際運行過程中，由于接觸斷面隨壓力波傳播沿管道運動，需要實時捕捉接觸斷面的位置，并對該區(qū)域?qū)嵤┚植考用?。因此，為了滿足計算精度和計算效率的平衡，在非均勻網(wǎng)格TVD格式基礎(chǔ)上，設(shè)計管道變網(wǎng)格計算方法。網(wǎng)格劃分初始為較粗的網(wǎng)格，計算過程中根據(jù)各網(wǎng)格間物理量波動情況，自適應(yīng)實現(xiàn)局部網(wǎng)格尺度的加密或合并。

2.1 基本原理

為了表征管道內(nèi)參數(shù)變化，定義光滑度指數(shù)為：

(21)

式中：φ為密度、速度、濃度等物理量，一維非定常流動基本方程組中，每個方程中均包含密度參數(shù)，可以選取密度梯度為光滑度指標(biāo)參數(shù)；Δxj-1,j表示網(wǎng)格j-1與網(wǎng)格j中點之間的距離，下標(biāo)ref表示該參數(shù)為參考值。光滑度指數(shù)σ大表示該區(qū)域參數(shù)變化劇烈，需要提高網(wǎng)格分辨率，即進行網(wǎng)格加密；當(dāng)光滑度指數(shù)σ較小時，參數(shù)變化平緩，已經(jīng)加密的網(wǎng)格可以進行合并操作，減少數(shù)據(jù)存儲空間及計算時間，提高計算效率。

1)當(dāng)σj≥ε時，當(dāng)前網(wǎng)格需進行加密操作，將當(dāng)前網(wǎng)格沿中點一分為二，變?yōu)?個網(wǎng)格，由于計算采用的是有限體積法，在粗網(wǎng)格中添加細網(wǎng)格后，能夠保證網(wǎng)格交界面上的流量守恒；

2)當(dāng)σj<ε時，當(dāng)前網(wǎng)格需進行合并操作，將當(dāng)前加密網(wǎng)格與其鄰近網(wǎng)格合并為一個網(wǎng)格。

有時，對于粗網(wǎng)格進行一次加密操作后仍不能滿足計算精度的要求，需要進行多級嵌套加密，加密或合并的閾值ε需要根據(jù)加密級數(shù)做相應(yīng)變化，即對于第k級加密：

εk=ε/2k

(22)

綜上所述，變網(wǎng)格加密計算基本步驟為：

1)參數(shù)輸入，給出網(wǎng)格尺度和初值、閾值；

2)初始化，劃分網(wǎng)格并給各網(wǎng)格參數(shù)賦初值；

3)根據(jù)當(dāng)前時刻各網(wǎng)格參數(shù)及閾值，確定該網(wǎng)格是否需要加密或合并，按照前面制定的網(wǎng)格加密和合并實施方案，計算各網(wǎng)格參數(shù)；

系統(tǒng)庫文件提供用戶程序，可直接使用API接口實現(xiàn)使用注冊、控制密鑰、加解密功能，內(nèi)部實現(xiàn)與內(nèi)核驅(qū)動之間的數(shù)據(jù)交互，支持設(shè)備復(fù)用功能[17]。

4)采用非均勻網(wǎng)格TVD求解格式，迭代下一時刻各網(wǎng)格點參數(shù)；

5)判斷計算是否完成，如果完成則退出計算，如果沒完成則進入第3)步繼續(xù)計算。

2.2 數(shù)據(jù)結(jié)構(gòu)設(shè)計

從上面的敘述不難看出，本文采用的網(wǎng)格加密生成方法為二分法，針對這種類型的數(shù)據(jù)生成方法，選用二叉樹的樹形存儲結(jié)構(gòu)用于計算過程的數(shù)據(jù)生成、存儲計算和數(shù)據(jù)刪除。如圖3所示，對于父網(wǎng)格j，實施加密后，在第二級分裂為2個子網(wǎng)格，網(wǎng)格編號為2j-1和2j。為了保證數(shù)據(jù)結(jié)構(gòu)的完整性和易操作性，在進行合并操作時，只能將同一父網(wǎng)格下的2個子網(wǎng)格合并，即將第2層的2j-1和2j節(jié)點合并為第1層的節(jié)點j。

2.3 網(wǎng)格數(shù)據(jù)轉(zhuǎn)換

圖3 變網(wǎng)格二叉樹網(wǎng)格劃分及數(shù)據(jù)存儲結(jié)構(gòu)示意Fig.3 Diagram of variable mesh binary tree mesh division and data storage structure

網(wǎng)格合并過程中，父網(wǎng)格的參數(shù)根據(jù)2個子網(wǎng)格的平均值求得，即:

(23)

與網(wǎng)格合并過程類似，網(wǎng)格加密過程中，子網(wǎng)格的參數(shù)根據(jù)父網(wǎng)格及其鄰近網(wǎng)格的參數(shù)插值求得，加密過程子網(wǎng)格插值公式為:

(24)

(25)

根據(jù)網(wǎng)格加密插值式(24)和(25)，反推網(wǎng)格合并過程，可以得到與網(wǎng)格合并相同的式(23)，保證了網(wǎng)格加密和合并過程的一致性。

3 變網(wǎng)格數(shù)值仿真應(yīng)用

3.1 激波管算例

一維非定常流動激波管問題具有精確解，可以仿真激波管中的氣體流動驗證非均勻網(wǎng)格TVD格式的準(zhǔn)確性。建立如圖4所示激波管算例，管道長1 m，中間由膜片分割，膜片兩側(cè)網(wǎng)格賦初值如圖所示，管道內(nèi)為理想氣體，絕熱系數(shù)為1.4，氣體初始流速均為0 m/s。

圖4 管道激波管算例示意Fig.4 Diagram of calculation case of pipeline shock wave tube

3.1.1 均勻網(wǎng)格不同網(wǎng)格密度對比

上述算例中，對比設(shè)置不同網(wǎng)格數(shù)(50、100、200)時，管道內(nèi)激波發(fā)生后t=0.5 ms時壓力、密度、速度和溫度分布情況，如圖5所示。從圖5中可以看出，在不同網(wǎng)格數(shù)情況下都能擬合出沿管道方向物理參數(shù)的變化，參數(shù)變化平緩的區(qū)段，其數(shù)值仿真結(jié)果較好，激波面仿真擬合效果較差。網(wǎng)格數(shù)越多，管道計算結(jié)果分辨率越高，模擬結(jié)果與精確解的擬合度越高。隨著網(wǎng)格密度的增加，仿真激波發(fā)生后t=0.5 ms所需的CPU時間分別為0.025、0.036、0.064 s，因此，隨著網(wǎng)格數(shù)增加，CPU計算所需時間增加，并且呈線性增長趨勢。

圖5 不同網(wǎng)格數(shù)管道激波管算例t=0.5 ms時參數(shù)分布Fig.5 The parameter distribution of shock wave tube with different mesh numbers under TVD format (t=0.5 ms)

3.1.2 非均勻網(wǎng)格與均勻網(wǎng)格對比

通過分析上述激波發(fā)生后t=0.5 ms時的算例參數(shù)分布，可以看出，該時刻管道中流動狀態(tài)發(fā)生較大變化的區(qū)域為以下3個區(qū)域：

A區(qū)域：0.1～0.3 m，該區(qū)域速度狀態(tài)由0至250 m/s連續(xù)變化，為左行波到達區(qū)域；

B區(qū)域：0.5～0.7 m，該區(qū)域密度溫度發(fā)生急劇變化，為發(fā)生接觸間斷的區(qū)域；

C區(qū)域：0.7～0.8 m，該區(qū)域速度狀態(tài)由250～0 m/s變化，為右行波到達區(qū)域。

存在物理參數(shù)急劇變化的區(qū)域是仿真計算產(chǎn)生誤差較大的區(qū)域，可以通過對局部區(qū)域網(wǎng)格加密以達到提高計算精度的目的。因此，針對激波管算例，將上述區(qū)域網(wǎng)格加密，管道非均勻網(wǎng)格的網(wǎng)格尺度分布示意如圖6所示。在A區(qū)域(0.1～0.3 m)，網(wǎng)格尺度ΔxA為0.01 m，與網(wǎng)格數(shù)為100的均勻網(wǎng)格網(wǎng)格密度相同；在B區(qū)域和C區(qū)域(0.5～0.8 m)，網(wǎng)格尺度ΔxB,C為0.005 m，與網(wǎng)格數(shù)為200的均勻網(wǎng)格網(wǎng)格密度相同；其他區(qū)域網(wǎng)格尺度為0.02 m，與網(wǎng)格數(shù)為50的均勻網(wǎng)格網(wǎng)格密度相同。采用非均勻網(wǎng)格激波管算例進行劃分的網(wǎng)格總數(shù)為110個。

激波管算例中，計算上述非均勻網(wǎng)格管道內(nèi)工質(zhì)激波發(fā)生后t=0.5 ms時壓力、密度、速度和溫度分布情況，并將結(jié)果與均勻網(wǎng)格(網(wǎng)格數(shù)100和200)方案進行對比，壓力和密度結(jié)果如圖7所示。在左行波區(qū)域A，由于非均勻網(wǎng)格110方案的網(wǎng)格密度與均勻網(wǎng)格100方案相同，所以2個方案在該區(qū)域各狀態(tài)參數(shù)分辨率相同，模擬結(jié)果仿真精度相當(dāng)，但仿真精度都略小于均勻網(wǎng)格200的算例；在接觸間斷區(qū)域B和右行波區(qū)域C，由于非均勻網(wǎng)格110方案的網(wǎng)格密度與均勻網(wǎng)格200方案相同，所以兩方案在該區(qū)域各狀態(tài)參數(shù)分辨率相同，模擬結(jié)果仿真精度相當(dāng)。結(jié)合表1所示的非均勻網(wǎng)格和均勻網(wǎng)格激波管算例誤差可以看出，均勻網(wǎng)格100的仿真精度最低，最大誤差為4.13%；采用與均勻網(wǎng)格100網(wǎng)格數(shù)量相當(dāng)?shù)姆蔷鶆蚓W(wǎng)格方案，不同參數(shù)的預(yù)測精度都有一定提高，非均勻網(wǎng)格的誤差均在3.1%以內(nèi)。綜上所述，采用非均勻網(wǎng)格方案后，在網(wǎng)格數(shù)增加不多的前提下，整體仿真精度較網(wǎng)格數(shù)同量極的均勻網(wǎng)格方案有大幅提高，局部高密度區(qū)域可以達到與2倍均勻網(wǎng)格密度相當(dāng)?shù)木取?/p>

圖6 管道非均勻網(wǎng)格劃分方案示意Fig.6 Diagram of non-uniform meshing scheme for pipeline

均勻網(wǎng)格100、均勻網(wǎng)格200、非均勻網(wǎng)格110等3種方案在仿真激波發(fā)生0.5 ms過程消耗CPU時間分別為0.036、0.064、0.040 s。非均勻網(wǎng)格110與均勻網(wǎng)格100方案的網(wǎng)格數(shù)相當(dāng)，仿真耗時較為接近。均勻網(wǎng)格200方案由于網(wǎng)格數(shù)較多，計算量較大，仿真耗時是非均勻網(wǎng)格110方案的1.6倍。

圖7 非均勻和均勻網(wǎng)格激波管算例t=0.5 ms時參數(shù)分布對比Fig.7 Comparison of parameter distribution between the non-uniform mesh and the uniform mesh solution format shock wave tube cases (t=0.5 ms)

表1 非均勻和均勻網(wǎng)格激波管算例誤差

綜上所述，采用非均勻網(wǎng)格，與同等均勻網(wǎng)格密度方案相比，能夠?qū)崿F(xiàn)在消耗相同CPU時間基礎(chǔ)上，提高仿真精度1.03個百分點；與加密均勻網(wǎng)格方案相比，能夠在滿足仿真精度相當(dāng)?shù)幕A(chǔ)上，降低計算耗時37.5%。

3.1.3 均勻網(wǎng)格與變網(wǎng)格對比

搭建相同的管道激波管算例，變網(wǎng)格格式初始網(wǎng)格數(shù)設(shè)置為50，閾值設(shè)置為5，加密最高級數(shù)設(shè)置為3，即最多能加密至200個網(wǎng)格。對比均勻網(wǎng)格(網(wǎng)格數(shù)為50和200)和變網(wǎng)格(初始網(wǎng)格數(shù)50)等不同數(shù)值仿真方法的仿真精度和仿真效率。

變網(wǎng)格格式與均勻網(wǎng)格格式數(shù)值仿真結(jié)果如圖8所示。從圖中可以看出，與均勻網(wǎng)格數(shù)為50求解格式相比，變網(wǎng)格格式仿真精度有明顯提高，尤其是在參數(shù)數(shù)值變化較為劇烈的激波接觸面，變網(wǎng)格格式更加逼近精確解。結(jié)合表2變網(wǎng)格與均勻網(wǎng)格激波管算例仿真誤差可以看出，變網(wǎng)格求解格式對于壓力、密度的數(shù)值仿真精度能夠達到與均勻網(wǎng)格200相當(dāng)?shù)乃剑鲄?shù)誤差均在2.2%以內(nèi)。

圖8 變網(wǎng)格與均勻網(wǎng)格求解格式管道激波管算例t=0.5 ms時參數(shù)分布對比Fig.8 Comparison of parameter distribution between variable grid and uniform grid solution format shock wave tube cases (t=0.5 ms)

表2 不同網(wǎng)格求解格式管道激波管算例誤差

圖9所示為使用變網(wǎng)格格式求解管道激波管算例，在t=0.5 ms時各加密級網(wǎng)格分布圖。沿管道長度方向生成了3類網(wǎng)格密度級別不同的區(qū)域，其中第1級網(wǎng)格的網(wǎng)格密度與初始網(wǎng)格相同，網(wǎng)格離散長度為0.02 m；第2級網(wǎng)格經(jīng)歷了一次網(wǎng)格加密過程，網(wǎng)格離散長度為0.01 m；第3級網(wǎng)格的網(wǎng)格密度最密集，相應(yīng)的網(wǎng)格離散長度也最小，其離散長度為0.005 m。結(jié)合圖8(b)可以發(fā)現(xiàn)，在密度變化梯度較大的3個區(qū)域，均生成了網(wǎng)格最密集的第3級網(wǎng)格，而在密度參數(shù)變化平緩的區(qū)域，仍然保留為初始的第1級網(wǎng)格，第1級和第3級網(wǎng)格間有少量第2級網(wǎng)格過渡。在初始網(wǎng)格數(shù)為50的基礎(chǔ)上，計算過程中實施了3級網(wǎng)格加密，仿真結(jié)束時，沿管道長度方向有第1級網(wǎng)格17個，第2級網(wǎng)格22個，第3級網(wǎng)格88個，3級網(wǎng)格總計有127個，與仿真精度相當(dāng)?shù)木鶆蚓W(wǎng)格數(shù)為200的算例相比網(wǎng)格數(shù)減少了73個，相應(yīng)的數(shù)據(jù)存儲量也降低了36.5%。

圖9 變網(wǎng)格求解格式管道激波管算例t=0.5 ms時各加密級網(wǎng)格分布Fig.9 Mesh distribution of each encryption layer in the variable-scale solution format pipeline shock wave tube cases

均勻網(wǎng)格50、均勻網(wǎng)格200、變網(wǎng)格50計算激波管算例所需的CPU時間分別為0.025、0.064、0.039 s。變網(wǎng)格求解格式所消耗的CPU時間比均勻網(wǎng)格數(shù)50略高。與均勻網(wǎng)格數(shù)200的算例相比，變網(wǎng)格求解格式所消耗的CPU時間僅是均勻網(wǎng)格數(shù)200的60.94%。采用變網(wǎng)格求解格式能夠在保持?jǐn)?shù)值仿真精度不降低的前提下，提高計算效率。

3.2 內(nèi)燃機排氣管內(nèi)的流動

為實現(xiàn)內(nèi)燃機整機性能仿真，除了進排氣管道流動仿真模型外，還需其他組件和邊界模型。參照相關(guān)文獻，其他組件和邊界模型選用已經(jīng)廣泛使用的模型，如表3所示。為便于實現(xiàn)與進排氣流動模型的耦合，邊界模型均基于特征線方法建立。

表3 其他組件模型Table 3 Other component models

搭建HC4132柴油機整機試驗臺架及測試系統(tǒng)，相關(guān)細節(jié)見文獻[11]，選取HC4132柴油機外特性1 500、1 700、1 900和2 100 r/min等4個工況點，其中2 100 r/min為柴油機標(biāo)定轉(zhuǎn)速，1 500 r/min為最大扭矩點轉(zhuǎn)速，測量不同轉(zhuǎn)速工況排氣壓力波以及整機性能指標(biāo)。搭建HC4132柴油機整機仿真模型，對比管道采用均勻網(wǎng)格和變網(wǎng)格不同轉(zhuǎn)速下排氣壓力波、整機性能預(yù)測精度以及計算效率。均勻網(wǎng)格算例網(wǎng)格數(shù)設(shè)置為10；管道變網(wǎng)格算例初始網(wǎng)格數(shù)設(shè)置為5，閾值設(shè)置為5，加密最高級數(shù)設(shè)置為3。

3.2.1 排氣壓力波對比

HC4132柴油機排氣系統(tǒng)采用的是2缸連接一根排氣管，即第1缸與第4缸共用一根排氣管，第2缸與第3缸共用一根排氣管，試驗過程中在第3缸排氣管處安裝傳感器用于測量排氣壓力波變化。對比管道采用均勻網(wǎng)格和變網(wǎng)格不同轉(zhuǎn)速下排氣壓力波仿真結(jié)果與試驗結(jié)果如圖10所示。第3缸排氣管測得的壓力波，其壓力波最大峰值出現(xiàn)在第3缸排氣階段；由于第3缸與第2缸鄰近，且兩缸連接在同一根排氣管上，因此第2缸排氣過程對第3缸排氣管處壓力波影響也較大。不同轉(zhuǎn)速工況，采用變網(wǎng)格計算方案，第3缸排氣壓力波波峰波谷處預(yù)測精度均有明顯提高。

圖10 不同轉(zhuǎn)速第3缸排氣壓力波仿真與試驗結(jié)果對比Fig.10 Comparison of calculated and measured exhaust pressure wave for the third cylinder under different speed

表4為不同轉(zhuǎn)速均勻網(wǎng)格和變網(wǎng)格第3缸排氣壓力波循環(huán)平均誤差，均勻網(wǎng)格和變網(wǎng)格循環(huán)平均誤差低轉(zhuǎn)速工況預(yù)測精度低于高轉(zhuǎn)速工況。1 500 r/min轉(zhuǎn)速工況計算誤差最大，均勻網(wǎng)格和變網(wǎng)格格式分別為10.62%和5.47%，采用變網(wǎng)格方案排氣壓力波預(yù)測精度最大可提高5.15個百分點。

表4 不同轉(zhuǎn)速不同網(wǎng)格格式排氣壓力波循環(huán)平均誤差

3.2.2 整機性能對比

HC4132柴油機主要性能參數(shù)不同求解格式仿真結(jié)果與穩(wěn)態(tài)試驗結(jié)果誤差如表5所示。從表中可以看出，不同求解格式仿真結(jié)果中功率、渦輪進口溫度等參數(shù)預(yù)測較準(zhǔn)確，各轉(zhuǎn)速工況誤差均在2.82%以內(nèi)，燃油消耗、中冷前溫度、渦輪進口壓力等參數(shù)預(yù)測精度雖然偏低，但最大誤差也在4.68%以內(nèi)。不同管道求解格式整機性能預(yù)測精度相當(dāng)，均勻網(wǎng)格最大誤差為4.68%，變網(wǎng)格最大誤差為4.14%，采用變網(wǎng)格整機性能預(yù)測精度較均勻網(wǎng)格提高0.54個百分點。

表5 HC4132不同求解格式整機性能仿真誤差

3.2.3 計算效率對比

程序運行平臺處理器主頻為4.20 GHz，分別進行4個工況計算，設(shè)置仿真20個循環(huán)，管道采用不同求解格式進行計算消耗CPU時間如表6所示，表中差值為變網(wǎng)格比均勻網(wǎng)格計算增加的耗時。變網(wǎng)格方案雖然初始網(wǎng)格數(shù)較少，但計算過程中需要生成和合并網(wǎng)格造成計算耗時增加，不同轉(zhuǎn)速工況采用變網(wǎng)格方案均略高于均勻網(wǎng)格方案。仿真計算20個循環(huán)，變網(wǎng)格方案平均增加耗時2.08 s。

表6 不同求解格式整機性能仿真消耗CPU時間

4 結(jié)論

1)激波管算例中，均勻網(wǎng)格TVD格式網(wǎng)格數(shù)越多，管道計算分辨率越高，模擬結(jié)果與精確解擬合度越高，但隨著網(wǎng)格數(shù)增加，CPU消耗時間呈線性增長；采用非均勻網(wǎng)格，與同等網(wǎng)格密度方案相比，能夠在消耗同等CPU時間基礎(chǔ)上，提高仿真精度1.03個百分點，與加密網(wǎng)格方案相比，在滿足仿真精度相當(dāng)?shù)幕A(chǔ)上，可以降低計算耗時37.5%。

2)初始網(wǎng)格數(shù)為50的變網(wǎng)格計算方法與網(wǎng)格數(shù)為200的均勻網(wǎng)格計算方法相比，兩者數(shù)值仿真精度相當(dāng)，誤差均在2.2%以內(nèi)，但前者數(shù)據(jù)存儲量比后者降低了36.5%，所消耗的CPU時間僅是后者的60.94%。因此，采用變網(wǎng)格數(shù)值仿真方法后，可以在保證數(shù)值仿真精度不降低的前提下，降低數(shù)據(jù)存儲量，減少CPU消耗時間，提高計算效率。

3)整機性能仿真中，采用變網(wǎng)格計算20個循環(huán)平均增加CPU耗時2.08 s，排氣壓力波預(yù)測精度較均勻網(wǎng)格方案最大可提高5.15個百分點，整機性能預(yù)測精度較均勻網(wǎng)格最大可提高0.54個百分點。