劉曉剛，明平劍，張格健，張文平

(哈爾濱工程大學(xué) 動力與能源工程學(xué)院，黑龍江 哈爾濱 150001)

?

柴油機(jī)缸內(nèi)三維動網(wǎng)格和流場數(shù)值模擬程序開發(fā)

劉曉剛，明平劍，張格健，張文平

(哈爾濱工程大學(xué) 動力與能源工程學(xué)院，黑龍江 哈爾濱 150001)

現(xiàn)較為成熟的CFD仿真軟件均為國外所有，國內(nèi)缺乏相應(yīng)的知識產(chǎn)權(quán)。為此，本文自主研發(fā)了一套仿真程序，對內(nèi)燃機(jī)缸內(nèi)三維工作過程進(jìn)行數(shù)值模擬。本文建立了動態(tài)層網(wǎng)格模型和匹配的任意邊界模型，采用有限體積方法離散流場控制方程，通過SIMPLEC算法求解NS方程，基于FORTRAN語言開發(fā)了一套內(nèi)燃機(jī)缸內(nèi)流場仿真軟件。對內(nèi)燃機(jī)的缸內(nèi)空氣純壓縮和膨脹過程做了仿真分析，計算結(jié)果表明，缸壓曲線、缸內(nèi)溫度曲線與應(yīng)用商業(yè)軟件PUMPLINX、FLUENT計算結(jié)果基本吻合,證明動網(wǎng)格模型和流場數(shù)值模擬方法正確，本文的程序有進(jìn)一步市場應(yīng)用的前景。

內(nèi)燃機(jī)；運動網(wǎng)格；CFD；數(shù)值模擬；程序開發(fā)；軟件仿真

內(nèi)燃機(jī)缸內(nèi)流場數(shù)值模擬主要是借助CFD手段預(yù)測缸內(nèi)的流動形式，從而在設(shè)計階段可預(yù)先改進(jìn)進(jìn)氣形式或燃燒室形狀，形成更好的缸內(nèi)氣流狀態(tài)以提高燃燒效率和降低污染物排放，在現(xiàn)代內(nèi)燃機(jī)設(shè)計和預(yù)先研究中發(fā)揮了越來越重要的作用[1]。而目前國內(nèi)的內(nèi)燃機(jī)設(shè)計上主要采用市面上較為流行的商用軟件，如KIVA系列軟件[2]、AVL Fire軟件、STAR-CD軟件、ANSYS Fluent軟件和新興的CONVERGE軟件等，購買和升級費用昂貴，沒有獨立的知識產(chǎn)權(quán)，一直受制于人，阻礙了中國內(nèi)燃機(jī)行業(yè)的發(fā)展。

哈爾濱工程大學(xué)自主研發(fā)的基于非結(jié)構(gòu)化網(wǎng)格的通用輸運方程求解軟件[3-5](general transport equation analysis,GTEA)采用有限體積方法離散流場控制方程，通過SIMPLEC算法求解NS方程和分離式變量求解方法，目前可以較為精確的求解二維和三維、層流和湍流的流場流動問題。本文在該軟件的基礎(chǔ)上，根據(jù)內(nèi)燃機(jī)的氣缸往復(fù)直線運動形式建立了動態(tài)層網(wǎng)格模型和匹配的任意邊界模型，并實現(xiàn)了其與流場求解器的耦合計算，形成了一套內(nèi)燃機(jī)缸內(nèi)流場專業(yè)化仿真軟件。

1 數(shù)學(xué)模型

數(shù)值計算中數(shù)學(xué)模型可分為計算域的幾何信息模型和數(shù)學(xué)控制方程模型。

1.1 運動網(wǎng)格模型

動態(tài)層模型[5]是利用增加或刪減網(wǎng)格層來適應(yīng)邊界的運動，該方法對網(wǎng)格的劃分及邊界的運動形式要求較高，邊界運動主要為往復(fù)式直線運動，網(wǎng)格需要在動邊界運動方向上分層。程序?qū)崿F(xiàn)時首先要建立結(jié)構(gòu)體數(shù)組來保存每層網(wǎng)格的單元和節(jié)點信息，以便進(jìn)行網(wǎng)格層間拓?fù)湫畔⒌牟檎?，之后再根?jù)計算域是膨脹還是壓縮來做網(wǎng)格層的加減，如下圖1所示，只對最下層的網(wǎng)格做加減層操作，以最大程度保持網(wǎng)格原有質(zhì)量。通過網(wǎng)格節(jié)點坐標(biāo)替換及網(wǎng)格標(biāo)記阻斷相結(jié)合的方法來滿足增減網(wǎng)格層時網(wǎng)格拓?fù)浣Y(jié)構(gòu)的變化，減層時1、2、3、4單元被阻斷，也就是不參與計算，F(xiàn)、G、H、I、J節(jié)點替換A、B、C、D、E節(jié)點坐標(biāo)；加層時增加1、2、3、4單元，增加A、B、C、D、E節(jié)點，并替換F、G、H、I、J節(jié)點坐標(biāo)。

圖1 增減網(wǎng)格層時層動邊界網(wǎng)格拓?fù)浣Y(jié)構(gòu)變化Fig.1 Topological structure change of the moving boundary mesh when the grid layer is changed

匹配的任意邊界網(wǎng)格模型[6]，即商用軟件中的r任意邊界模型，是指由于相對運動等原因在同一位置有兩個邊界存在，再通過其上的邊界面找到對應(yīng)關(guān)系以傳遞流場中的信息，以按內(nèi)部單元面處理。由于動態(tài)層網(wǎng)格在三維情況下只適用于規(guī)則的六面體或者三棱柱網(wǎng)格單元，且本文只對初始網(wǎng)格進(jìn)行編號，所以要將整個計算域劃分為氣缸區(qū)域和燃燒室區(qū)域，在兩個區(qū)域結(jié)合處采用任意邊界面的處理方法。本文為處理方便在網(wǎng)格劃分時將這兩個邊界上的網(wǎng)格一一對應(yīng)(在Gambit軟件中選中Link mesh face選項)，這樣兩個邊界面上的所有網(wǎng)格幾何信息全部相同，單元面間的相交判斷就可以通過邊界單元面的中心坐標(biāo)是否相同來實現(xiàn)，進(jìn)而確定邊界單元面間的對應(yīng)關(guān)系和所屬單元。

1.2 流場控制方程

流體的運動要滿足基本守恒定律，其控制方程組是各種CFD軟件求解的基本數(shù)學(xué)模型，通用輸運方程的積分形式可寫為

(1)

式中:ρ是流體密度，U和Ub分別是流體和網(wǎng)格面運動速度，圖1中紅色運動邊界處的邊界單元面需要計算其Ub，本文采用上一時間段平均移動速度計算，n是單元面的外法線方向單位向量，V和A分別為網(wǎng)格單元體積和單元面的面積，φ為通用形式控制方程的求解變量，Γφ為對應(yīng)的擴(kuò)散系數(shù)，Sφ為相應(yīng)的源項。各方程的具體形式如表1所示。

表1 擴(kuò)散系數(shù)和源項

表中，μeff為有效粘性系數(shù)，求解層流問題時μeff=μl，求解湍流時μeff=μl+μt，H為流體總焓，λ為導(dǎo)熱系數(shù)，cp為定壓比熱，Prt為湍流普朗特數(shù)，τ為流體所受切應(yīng)力，k-ε系列湍流方程同樣可以寫為上述形式，本文算例沒有采用湍流方程，這里不再贅述。

2 數(shù)值方法

2.1 控制方程離散

時間項的離散計算類似于一階常微分方程，本文采用歐拉隱式格式對其進(jìn)行計算，近似可得離散形式為

(2)

式中，Δt為計算時間步長、(ρVPφ)o為上一時刻值，F(xiàn)(φ,t)為通用輸運方程中除時間項外的所有其他項當(dāng)前時刻值之和，所以這種格式稱作為隱式格式。

通過對對流項的離散，其可以近似等于各單元面上的對流通量之和：

(3)

式中：nf為當(dāng)前單元所含有的單元面?zhèn)€數(shù)，fj為通過該單元第j個單元面的流量，當(dāng)流速與單元面外方向n相同時為正，否則為負(fù)。單元面上變量值φj的取值取決于對流項采用的格式，常采用的一階迎風(fēng)格式：

(4)

則離散后的對流項最終可寫為

(5)

擴(kuò)散項的離散要稍微復(fù)雜一些，將變量梯度分為法向n和交叉向τ兩部分，略去中間計算過程，最終可得：

(6)

式(6)中右端第一項為法向擴(kuò)散部分，第二項為交叉向擴(kuò)散部分，近似離散為

(7)

(8)

(9)

通過有限體積法離散流體通用控制方程，最終可得代數(shù)方程組目標(biāo)形式為

(10)

式中：φP和φPj分別為當(dāng)前單元和第j個相鄰單元的變量值，ap和aj分別為代數(shù)方程中當(dāng)前單元和相鄰單元變量系數(shù)，bp為代數(shù)方程源項，本文通過代數(shù)多重網(wǎng)格和共軛梯度法求解代數(shù)方程組。

2.2 時間項與對流項特殊處理

為方便程序計算，要對時間項和對流項進(jìn)行特殊處理。采用一階迎風(fēng)格式離散對流項，可進(jìn)一步寫成如下形式：

(11)

可知，當(dāng)fj>0時，fj=fout；當(dāng)fj<0時，fj=fin。同樣經(jīng)過上述處理，忽略源項，離散后的連續(xù)方程可寫為

(12)

將其代入式(11)，對流項最終可得：

(13)

式中:fin只能等于零或-fj，所以可讓fin=-max(0, fj)。

離散后的時間項可進(jìn)一步整理得：

(14)

最后將時間項和對流項合并在一起，可得：

(15)

根據(jù)上述各項對代數(shù)方程組的貢獻(xiàn)，經(jīng)整理各系數(shù)可表示為

這樣就會對每個方程形成代數(shù)方程組，為使方程求解穩(wěn)定，本程序采用求解小量形式進(jìn)行求解，式(10)可變?yōu)?/p>

(16)

式中：φ'=φ-φo。GTEA采用分離式變量求解方法，即逐一求解各個方程組，本文采用的共軛梯度法和代數(shù)多重網(wǎng)格法來求解代數(shù)方程組，進(jìn)而求得各變量。

2.3 非結(jié)構(gòu)化網(wǎng)格SIMPLEC算法

動量預(yù)測是SIMPLE系列算法[7]的第一步，其主要思想是根據(jù)初始給定的壓力p*和ρ*，通過求解動量方程得出預(yù)測步的U*。將動量方程中的壓力梯度項提出進(jìn)行單獨處理，可得離散后動量方程的代數(shù)方程形式為

(17)

則可通過上式求得預(yù)測步的各單元中心的速度。

壓力修正的主要思想是利用動量預(yù)測得出的單元中心的速度，再用求解動量方程相似的形式來計算網(wǎng)格面上的速度，即動量加權(quán)或Rhie-Chow插值的方法；并采取壓力、速度、密度修正的方法，用壓力修正量表示速度和密度的修正量，最后修正單元面的流量，將其代入連續(xù)方程。SIMPLEC算法得到的壓力修正量和速度修正量之間的關(guān)系為

(18)

省略中間計算,最后得出壓力修正方程的形式如下

(19)

動量修正的主要目的是通過壓力修正方程求出的壓力修正量來修正速度和壓力。雖然速度修正量與壓力修正量之間的關(guān)系是通過動量方程求出，但通過壓力修正方程之后，修正后的速度是嚴(yán)格滿足連續(xù)方程的，修正后的壓力也同樣滿足動量方程。程序中修正的公式為

(20)

(21)

而最后修正的密度是通過狀態(tài)方程和已修正后的速度計算得出。

圖2所示為GTEA軟件的計算流程圖。

圖2 程序計算流程圖Fig.2 Program calculation flow chart

3 算例驗證

為驗證本文內(nèi)燃機(jī)缸內(nèi)動網(wǎng)格模型及其流場計算的正確與否，先后設(shè)計了以下的動態(tài)層網(wǎng)格模型流場算例和匹配的任意邊界網(wǎng)格模型流場算例，并將計算結(jié)果與商業(yè)軟件的計算結(jié)果進(jìn)行了對比。

3.1 動態(tài)層網(wǎng)格模型驗證

1)測試算例

任意選取了一個氣缸模型，活塞的半徑為0.017m，沖程為0.05m，曲軸轉(zhuǎn)速為1 000r/min，回轉(zhuǎn)半徑為0.024 75m，連桿長度為0.1m。計算的壓力初值為100 000Pa，溫度初值為300K。采用六面體網(wǎng)格單元，共劃分了550個網(wǎng)格，初始網(wǎng)格模型如圖3所示，計算時間為180°BTDC到180°ATDC，計算得到的平均壓力和溫度曲線如圖4所示。

由圖4可知，計算所得的缸內(nèi)平均壓力和平均溫度曲線與商業(yè)軟件PUMPLINX計算所得結(jié)果幾乎完全吻合。

圖3 測試算例網(wǎng)格模型Fig.3 Grid model of test case

圖4 氣缸平均壓力、溫度曲線對比Fig.4 Cylinder mean pressure and mean temperature curve comparison

圖5 TBD620內(nèi)燃機(jī)氣缸模型Fig.5 Internal combustion engine cylinder model of TBD620

2)TBD620氣缸算例

TBD620內(nèi)燃機(jī)活塞的半徑為0.085m，沖程為0.2m；曲軸轉(zhuǎn)速為1 800r/min，回轉(zhuǎn)半徑為0.097 5m，連桿長度為0.35m。計算的壓力初值為30 000Pa，溫度初值為380K。采用六面體網(wǎng)格單元，共劃分了25 752個網(wǎng)格。初始缸內(nèi)計算域網(wǎng)格模型如圖5所示，計算得到的平均壓力和溫度曲線如圖6所示。

由圖6可知，計算所得的缸內(nèi)壓力和溫度曲線與商業(yè)軟件FLUENT的計算所得結(jié)果幾乎完全吻合。

3.2 燃燒室流場仿真驗證

為了模擬完整的缸內(nèi)過程流場仿真，還要加上燃燒室部分，建立的TBD620完整的缸內(nèi)燃燒室模型如圖7所示，氣缸部分為六面體網(wǎng)格，燃燒室部分為金字塔和四面體混合單元，結(jié)合面處為一一對應(yīng)的四邊形面單元，網(wǎng)格數(shù)為10 786。

如圖8所示為缸內(nèi)壓縮和膨脹過程所得到的壓力和溫度曲線與FLUENT計算結(jié)果對比圖，可看出結(jié)果還是基本完全重合，說明了任意邊界網(wǎng)格模型與流場耦合正確。

圖6 氣缸平均壓力、溫度曲線對比Fig.6 Cylinder mean pressure and mean temperature curve comparison of TBD620

圖7 TBD620內(nèi)燃機(jī)氣缸Fig.7 Complete model of engine of TBD620

如圖9所示為10°BTDC和10°ATDC時的氣缸縱截面矢量圖和全域速度矢量圖的仰視圖，可以發(fā)現(xiàn)活塞在上行到上止點之前，四周的氣體流向燃燒室，會在燃燒室的兩個凹坑內(nèi)形成非常明顯滾流漩渦，同理當(dāng)活塞下行到上止點后，燃燒室內(nèi)的氣體由流出到活塞定于氣缸蓋的空間，會在燃燒室凹坑內(nèi)形成兩個逆向的滾流漩渦，可見這種氣體流動形式可以使燃料與空氣混合完全，符合柴油機(jī)工作時的實際情況有利于燃料燃燒。

圖8 氣缸平均壓力、溫度曲線對比Fig.8 Cylinder mean pressure and mean temperature curve comparison of TBD620

圖9 流場速度矢量圖Fig.9 Flow velocity vector diagram

4 結(jié)論

本文在動態(tài)層網(wǎng)格模型和任意邊界網(wǎng)格模型的基礎(chǔ)上，采用有限體積法離散流場控制方程，SIMPLEC算法求解NS方程，開發(fā)了一套三維缸內(nèi)動網(wǎng)格流場數(shù)值模擬程序，即自主程序GTEA。并利用氣缸算例和燃燒室算例進(jìn)行驗證，將GTEA仿真模擬所得的缸內(nèi)平均壓力值和缸內(nèi)平均溫度值與商業(yè)軟件PUMPLINX和FLUENT計算結(jié)果對比，發(fā)現(xiàn)本文所得的缸內(nèi)曲線與商業(yè)軟件所得結(jié)果幾乎完全吻合，并且給出了TBD620柴油機(jī)的缸內(nèi)流場速度矢量圖，就目前來講，本自主程序可為相關(guān)仿真研究提供一定的指導(dǎo)和參考意義。

[1]舒歌群, 馬維忍, 許世杰, 等. 噴霧夾角對柴油機(jī)性能影響的數(shù)值模擬[J]. 工程熱物理學(xué)報, 2008, 29(7): 1239-1242. SHU Gequn, MA Weiren, XU Shijie, et al. Simulation of the effect of spray angle on diesel performance[J]. Journal of engineering thermophysics, 2008, 29(7): 1239-1242.

[2]AMSDEN A A, O'ROURKE P J, BUTLER T D. KIVA-II: a computer program for chemically reactive flows with sprays[R]. Technical Report LA-11560-MS. Los Alamos, New Mexico: Los Alamos National Laboratory, 1989.

[3]明平劍. 基于非結(jié)構(gòu)化網(wǎng)格氣液兩相流數(shù)值方法及并行計算研究與軟件開發(fā)[D]. 哈爾濱: 哈爾濱工程大學(xué), 2008. MING Pingjian. Development of numerical modeling for gas-liquid two-phase flows based on unstructured grids and parallel computing[D]. Harbin: Harbin Engineering University, 2008.

[4]雷國東. 非結(jié)構(gòu)網(wǎng)格FVM在復(fù)雜幾何結(jié)構(gòu)的湍流反應(yīng)流計算中的應(yīng)用研究[D]. 哈爾濱: 哈爾濱工程大學(xué), 2008. LEI Guodong. The application and research of the unstructured grid FVM in the turbulent reaction flows simulation with complex geometries[D]. Harbin: Harbin Engineering University, 2008.

[5]張志榮, 冉景煜, 張力, 等. 內(nèi)燃機(jī)缸內(nèi)氣體CFD瞬態(tài)分析中動態(tài)網(wǎng)格劃分技術(shù)[J]. 重慶大學(xué)學(xué)報: 自然科學(xué)版, 2005, 28(11): 97-100, 121. ZHANG Zhirong, RAN Jingyu, ZHANG Li, et al. Techniques of dynamic mesh in the CFD transient analysis of an engine[J]. Journal of Chongqing university: natural science edition, 2005, 28(11): 97-100, 121.

[6]孫華文. 基于非結(jié)構(gòu)動網(wǎng)格的內(nèi)燃機(jī)瞬態(tài)流場仿真[D]. 哈爾濱: 哈爾濱工程大學(xué), 2014. SUN Huawen. Numerical simulation of the transient flow in diesel engine based on unstructured dynamic meshes[D]. Harbin: Harbin Engineering University, 2014.

[7]PATANKAR S V, SPALDING D B. A calculation procedure for heat, mass and momentum transfer in three-dimensional parabolic flows[J]. International journal of heat and mass transfer, 1972, 15(10): 1787-1806.

[8]劉永豐. 基于非結(jié)構(gòu)網(wǎng)格的缸內(nèi)兩相反應(yīng)流數(shù)值模擬方法研究及軟件開發(fā)[D]. 哈爾濱: 哈爾濱工程大學(xué), 2013. LIU Yongfeng. Development of numerical simulation for the process of two phases reacting flow in cylinder based on unstructured grids[D]. Harbin: Harbin Engineering University, 2013.

[9]RHIE C M, CHOW W L. Numerical study of the turbulent flow past an airfoil with trailing edge separation[J]. AIAA journal, 1983, 21(11): 1525-1532.

[10]王福軍. 計算流體動力學(xué)分析-CFD軟件原理與應(yīng)用[M]. 北京: 清華大學(xué)出版社, 2004: 1-2, 25-28. WANG Fujun. Computational fluid dynamics analysis-theory and application of CFD software[M]. Beijing: Tsinghua University Press, 2004: 1-2, 25-28.

[11]PERAIRE J, VAHDATI M, MORGAN K, et al. Adaptive remeshing for compressible flow computations[J]. Journal of computational physics, 1987, 72(2): 449-466.

[12]BATINA J T. Unsteady Euler airfoil solutions using unstructured dynamic meshes[J]. AIAA journal, 1990, 28(8): 1381-1388.

Development of numerical simulation program for 3D moving grid and flow field in the cylinder of a diesel engine

LIU Xiaogang, MING Pingjian, ZHANG Gejian, ZHANG Wenping

(College of Power and Energy Engineering, Harbin Engineering University, Harbin 150001, China)

As mature CFD simulation software is currently developed overseas and domestic companies lack intellectual property rights, to ameliorate this problem, Harbin Engineering University developed a set of simulation programs that, after many years of accumulation, have proven to be remarkably successful. To continue conducting a numerical simulation of the three-dimensional working process in the cylinder of an internal combustion engine, a dynamic layer model and a matching model with an arbitrary boundary were firstly established. Moreover, the finite volume method was used to discretize the flow field control equations, and the Navier-Stokes equation was solved using the SIMPLEC algorithm. Based on the FORTRAN language, a set of simulation software of the internal flow field inside the cylinder of a diesel engine was then developed. The air compression and expansion process inside the cylinder were subsequently simulated and analyzed. Moreover, results showed that the cylinder pressure curve and cylinder temperature curve agree well with results of commercial software PUMPLINX and FLUENT, testifying that the moving grid model and the flow field simulation method are correct and that the numerical simulation program has potential as a market application.

I.C. Engine; moving grid; CFD; numerical simulation; program development; software simulation

2015-09-08.

日期：2016-10-12.

國家自然科學(xué)基金項目(51206031,51479038);國防預(yù)研基金項目(034315003).

劉曉剛(1990-)，男，碩士； 明平劍(1980-)，男，副教授，博士生導(dǎo)師.

明平劍, E-mail: pingjianming@hrbeu.edu.cn.

10.11990/jheu.201509023

TK402

A

1006-7043(2016) 11-1485-07

劉曉剛，明平劍，張格健，等. 柴油機(jī)缸內(nèi)三維動網(wǎng)格和流場數(shù)值模擬程序開發(fā)[J]. 哈爾濱工程大學(xué)學(xué)報, 2016, 37(11): 1485-1491. LIU Xiaogang, MING Pingjian, ZHANG Gejian, et al. Development of numerical simulation program for 3D moving grid and flow field in the cylinder of a diesel engine[J]. Journal of Harbin Engineering University, 2016, 37(11): 1485-1491.

網(wǎng)絡(luò)出版地址：http://www.cnki.net/kcms/detail/23.1390.u.20161012.0927.004.html