張琨， 陳昊天， 鄧康耀， 胡志龍， 錢(qián)躍華， 劉博

(1.上海交通大學(xué) 動(dòng)力機(jī)械及工程教育部重點(diǎn)實(shí)驗(yàn)室， 上海 200240； 2.上海船舶設(shè)備研究所， 上海 200031； 3.中國(guó)船舶集團(tuán)有限公司第七一一研究所， 上海 201108； 4.中船動(dòng)力研究院有限公司， 上海 200120)

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

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

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

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

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

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

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

(1)

其中：

一維非定常流動(dòng)方程組的Jacobi矩陣為[7]：

(2)

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

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

1.2 源項(xiàng)處理

對(duì)于帶有源項(xiàng)的一維非定常流動(dòng)控制方程組式(1)，在從n時(shí)刻推進(jìn)至n+1時(shí)刻時(shí)，可將其分裂為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格式

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

定義新的特征變量:

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)的證明，對(duì)于迎風(fēng)迭代格式:

(12)

(13)

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

(14)

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

(15)

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

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

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

(16)

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

(17)

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

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

(18)

(19)

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

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

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

(20)

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

2.1 基本原理

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

(21)

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

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

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

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

εk=ε/2k

(22)

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

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

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

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

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

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

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

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

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

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

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

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

(23)

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

(24)

(25)

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

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

3.1 激波管算例

一維非定常流動(dòng)激波管問(wèn)題具有精確解，可以仿真激波管中的氣體流動(dòng)驗(yàn)證非均勻網(wǎng)格TVD格式的準(zhǔn)確性。建立如圖4所示激波管算例，管道長(zhǎng)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)格密度對(duì)比

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

圖5 不同網(wǎng)格數(shù)管道激波管算例t=0.5 ms時(shí)參數(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)格對(duì)比

通過(guò)分析上述激波發(fā)生后t=0.5 ms時(shí)的算例參數(shù)分布，可以看出，該時(shí)刻管道中流動(dòng)狀態(tài)發(fā)生較大變化的區(qū)域?yàn)橐韵?個(gè)區(qū)域：

A區(qū)域：0.1～0.3 m，該區(qū)域速度狀態(tài)由0至250 m/s連續(xù)變化，為左行波到達(dá)區(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變化，為右行波到達(dá)區(qū)域。

存在物理參數(shù)急劇變化的區(qū)域是仿真計(jì)算產(chǎn)生誤差較大的區(qū)域，可以通過(guò)對(duì)局部區(qū)域網(wǎng)格加密以達(dá)到提高計(jì)算精度的目的。因此，針對(duì)激波管算例，將上述區(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)格激波管算例進(jìn)行劃分的網(wǎng)格總數(shù)為110個(gè)。

激波管算例中，計(jì)算上述非均勻網(wǎng)格管道內(nèi)工質(zhì)激波發(fā)生后t=0.5 ms時(shí)壓力、密度、速度和溫度分布情況，并將結(jié)果與均勻網(wǎng)格(網(wǎng)格數(shù)100和200)方案進(jìn)行對(duì)比，壓力和密度結(jié)果如圖7所示。在左行波區(qū)域A，由于非均勻網(wǎng)格110方案的網(wǎng)格密度與均勻網(wǎng)格100方案相同，所以2個(gè)方案在該區(qū)域各狀態(tài)參數(shù)分辨率相同，模擬結(jié)果仿真精度相當(dāng)，但仿真精度都略小于均勻網(wǎng)格200的算例；在接觸間斷區(qū)域B和右行波區(qū)域C，由于非均勻網(wǎng)格110方案的網(wǎng)格密度與均勻網(wǎng)格200方案相同，所以?xún)煞桨冈谠搮^(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ù)測(cè)精度都有一定提高，非均勻網(wǎng)格的誤差均在3.1%以?xún)?nèi)。綜上所述，采用非均勻網(wǎng)格方案后，在網(wǎng)格數(shù)增加不多的前提下，整體仿真精度較網(wǎng)格數(shù)同量極的均勻網(wǎng)格方案有大幅提高，局部高密度區(qū)域可以達(dá)到與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過(guò)程消耗CPU時(shí)間分別為0.036、0.064、0.040 s。非均勻網(wǎng)格110與均勻網(wǎng)格100方案的網(wǎng)格數(shù)相當(dāng)，仿真耗時(shí)較為接近。均勻網(wǎng)格200方案由于網(wǎng)格數(shù)較多，計(jì)算量較大，仿真耗時(shí)是非均勻網(wǎng)格110方案的1.6倍。

圖7 非均勻和均勻網(wǎng)格激波管算例t=0.5 ms時(shí)參數(shù)分布對(duì)比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時(shí)間基礎(chǔ)上，提高仿真精度1.03個(gè)百分點(diǎn)；與加密均勻網(wǎng)格方案相比，能夠在滿足仿真精度相當(dāng)?shù)幕A(chǔ)上，降低計(jì)算耗時(shí)37.5%。

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

搭建相同的管道激波管算例，變網(wǎng)格格式初始網(wǎng)格數(shù)設(shè)置為50，閾值設(shè)置為5，加密最高級(jí)數(shù)設(shè)置為3，即最多能加密至200個(gè)網(wǎng)格。對(duì)比均勻網(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)格求解格式對(duì)于壓力、密度的數(shù)值仿真精度能夠達(dá)到與均勻網(wǎng)格200相當(dāng)?shù)乃?，各參?shù)誤差均在2.2%以?xún)?nèi)。

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

圖9 變網(wǎng)格求解格式管道激波管算例t=0.5 ms時(shí)各加密級(jí)網(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計(jì)算激波管算例所需的CPU時(shí)間分別為0.025、0.064、0.039 s。變網(wǎng)格求解格式所消耗的CPU時(shí)間比均勻網(wǎng)格數(shù)50略高。與均勻網(wǎng)格數(shù)200的算例相比，變網(wǎng)格求解格式所消耗的CPU時(shí)間僅是均勻網(wǎng)格數(shù)200的60.94%。采用變網(wǎng)格求解格式能夠在保持?jǐn)?shù)值仿真精度不降低的前提下，提高計(jì)算效率。

3.2 內(nèi)燃機(jī)排氣管內(nèi)的流動(dòng)

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

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

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

3.2.1 排氣壓力波對(duì)比

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

圖10 不同轉(zhuǎn)速第3缸排氣壓力波仿真與試驗(yàn)結(jié)果對(duì)比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ù)測(cè)精度低于高轉(zhuǎn)速工況。1 500 r/min轉(zhuǎn)速工況計(jì)算誤差最大，均勻網(wǎng)格和變網(wǎng)格格式分別為10.62%和5.47%，采用變網(wǎng)格方案排氣壓力波預(yù)測(cè)精度最大可提高5.15個(gè)百分點(diǎn)。

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

3.2.2 整機(jī)性能對(duì)比

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

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

3.2.3 計(jì)算效率對(duì)比

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

表6 不同求解格式整機(jī)性能仿真消耗CPU時(shí)間

4 結(jié)論

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

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

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