一種可應用于內燃機瞬態仿真的動網格模型

孫華文， 楊麗紅， 明平劍， 張文平

(1. 國家超級計算天津中心， 天津 300457；2. 國家知識產權局專利局專利審查協作天津中心， 天津 300304；3. 哈爾濱工程大學動力與能源工程學院， 黑龍江 哈爾濱 150001)

?

一種可應用于內燃機瞬態仿真的動網格模型

孫華文1， 楊麗紅2， 明平劍3， 張文平3

(1. 國家超級計算天津中心， 天津 300457；2. 國家知識產權局專利局專利審查協作天津中心， 天津 300304；3. 哈爾濱工程大學動力與能源工程學院， 黑龍江 哈爾濱 150001)

提出了一種基于非結構網格的動態層網格實現算法，結合滑移網格算法構建了基于分塊滑移動態層的非結構化內燃機動網格模型，并基于TBD620柴油機建立了計算模型，所有算法都基于課題組自主研發的通用輸運方程求解軟件實現。流場計算采用適用于可壓縮流場的有限體積法及SIMPLE算法。通過數值算例對所開發的滑移網格模型和動態層網格模型進行了驗證，最后對內燃機缸內瞬態流場進行了仿真。計算結果表明，所發展的非結構化動網格模型可應用于內燃機瞬態流場的仿真。

內燃機； 滑移網格； 動態層算法； 數值模擬

數值模擬作為一種有效的分析設計方法已廣泛應用于內燃機的結構設計及性能分析中，然而現行的內燃機三維仿真分析軟件大多來自國外，國內自主研發的相對較少。內燃機缸內瞬態流場仿真的關鍵技術主要包括運動網格處理、噴霧模型、燃燒及化學反應機理以及缸內流場數值求解方法等，其中處理好網格運動是一切物理模型應用的前提[1]。

對此國內外學者進行了針對性研究。A.Velghe等[2]采用全局的網格重構插值方法對內燃機的瞬態流場進行了仿真；D.Abouri等[3]研究了基于多面體網格的運動網格模型并應用于內燃機瞬態流場的計算；國內蔣炎坤、羅馬吉等[4]研究了帶氣口發動機的瞬態流動計算，并利用SNAPPER技術解決運動邊界問題；劉金武等[5]基于網格的拆解組合給出了內燃機缸內復雜空間三維動態網格的生成方法；秦文瑾等[6]基于KIVA3V進行了4 氣門直噴式汽油機缸內湍流場多周期循環變動的大渦模擬研究，其中網格的運動主要依賴于KIVA3V中的自帶模型?？傮w來說國外的內燃機動網格研究相對成熟，國內學者大多基于國外成熟算法或開源代碼進行研究及改進，基于自主研發軟件開發內燃機仿真平臺研究內燃機動網格的較少。

本研究基于自主研發的通用輸運方程求解軟件[7-8]研究運動網格模型的實現算法，提出一種基于非結構網格的動態層網格實現算法，并針對內燃機運動網格的需求結合滑移網格算法提出了分塊滑移動態層內燃機運動網格模型的解決方案，開發相應模塊，為內燃機瞬態流場的仿真奠定基礎。

1 數學模型

1.1 有限體積法離散

內燃機缸內流場的各基本控制方程最后都可以寫成通用的輸運方程形式：

(1)

方程由4個不同類型的項組成，依次稱為瞬態項、對流項、擴散項和源項。其中源項又可以分為面積源項和體積源項。式中：φ為通用形式的場變量；Γφ為與φ相對應的擴散系數；SφA和SφV為與φ相對應的面積源項和體積源項。方程離散的目的是將方程轉化成代數方程組的形式：

(2)

式中：nf為當前研究單元的面個數；j代表第j個面；φc0和φcj分別為當前所研究單元和第j個相鄰單元在方程中的變量值；ap和aj分別為方程組中當前研究單元和第j個相鄰單元在方程中的系數；bp為方程的源項。

方程的離散采用基于非結構同位網格單元中心的有限體積法，網格的形狀可以是任意形狀或者多種形狀的混合網格。瞬態項采用中心差分的隱式格式求解，對流項采用一階迎風格式，擴散項采用法向交叉方法。根據各項系數的貢獻，可得方程組的系數分別為

ap=apT+apC+apD=

(3)

aj=ajT+ajC+ajD=

(4)

bp=bpT+bpC+bpD=

(5)

1.2 數值求解

采用非結構化網格的SIMPLE算法進行流場求解，采用Rhie-Chow[9]插值方法處理壓力速度失耦問題，通過動量方程預測流場，再根據壓力修正方程結果對壓力場、速度場和單元面流量進行修正。適用于任意多面體網格及移動網格，可求解可壓及不可壓流動。基于SIMPLE算法求解包含移動邊界流場的基本流程見圖 1。

2 運動網格算法

2.1 滑移網格算法

由于網格區域間的滑移將會產生滑移面上網格間對應關系的不匹配(見圖 2)，要想知道網格間的對應關系，就要進行網格單元面之間的相交判斷，本研究采用AABB(Axis Aligned Bounding Box)算法[10]進行計算。下面以二維情況為例說明AABB算法的基本原理。

二維時網格的單元面為線段，圖3示出兩相交單元，其中線段AB及線段CD需要進行相交判斷，并假設線段AB為主要邊界上單元面。AABB方法主要是采用建立包圍盒來包圍所需判斷單元面，二維時包圍盒為四邊形，圖4中的四邊形1及四邊形2分別為線段AB及線段CD的包圍盒，由于包圍盒的四邊分別與坐標軸平行，因而只需存儲四邊形中兩個節點坐標便可以確定此包圍盒。為方便后續比較，采用存儲XY坐標均最小的點及XY坐標均最大的點的坐標來確定包圍盒，計算得到擴展包圍盒BOX后，通過比較BOX與BOX_CD的最大最小點坐標來判斷BOX是否包含BOX_CD，如果包含則表明單元面AB與單元面CD相交。由此便完成了AABB算法的二維單元面相交判斷，該方法同樣也可以應用于三維情況。

通過對兩個不匹配邊界的單元面進行循環便可得到不匹配邊界上網格的對應關系，形成新的內部單元面，具體流程見圖5。

2.2 動態層網格模型

動態層模型是利用增加或刪減網格層來適應邊界的運動，該方法對網格的劃分及邊界的運動形式要求較高，邊界運動主要為往復式直線運動，網格需要在動邊界運動方向上分層，因而常用于結構化網格，也稱為層動網格。對于非結構網格存在網格間拓撲不規則的問題，主要通過網格層間的拓撲結構進行網格層的查找、增加及刪減。

自主軟件讀入網格文件后，非結構網格間的基本拓撲關系便為已知信息，包括單元面兩側單元序號、構成單元的單元面序號等，且節點在單元中均有一定的排序，存儲方式可參見文獻[11]。本研究網格動態層的分層存儲便是基于已知的非結構網格間的拓撲結構進行的，重新查找存儲形成一個結構體數組來記錄每層網格的信息。如圖6網格含有一個層動邊界，DE,EF,GH,HI代表著單元面的序號，所形成的結構體數組見表1。

獲得層動網格信息后，在邊界運動過程中，通過網格節點坐標替換及網格標記阻斷相結合的方法來滿足增減網格層時網格拓撲結構的變化。以圖 7中網格2為例，當網格2由于層動邊界運動被壓扁需要刪除時，只需將節點E及節點F的坐標分別賦給節點C和節點D，節點E及節點F的坐標恢復初始位置，之后將網格2標記為無效單元便完成網格2的刪除，這樣遍歷層動邊界上的網格便可達到網格層的刪除。采用網格層存儲算法得到的網格層信息中已將需要相互替換節點的單元面查找出并將節點排序一一對應，因此在網格刪除及添加時直接替換即可。這種方法不需要網格及網格面的重新排序，整個過程網格結構不發生改變，且可以處理包含多個層動邊界的情況。增加網格層時同上。所采用增加或刪減網格層的判據為給定一個(0,1)之間的系數a，當層動邊界網格高h(1+a)·h0時進行網格層的增加。

2.3 三維內燃機運動網格模型

以TBD620柴油機為例，帶有進排氣道的缸內流場模型見圖8，其中左側氣閥為進行閥，右側氣閥為排氣閥，且均處于全開狀態，其中心截面見圖9。

圖10中區域劃分依據分塊滑移動態層法的需要進行，其中區域1，2，3，4，7，8，11為層動區域，因而這些區域的網格劃分需要在層動邊界的移動方向上分層。區域1，2，3，4，7，8的下邊界均為層動邊界，上邊界為內部邊界，即與相連區域公用此邊界；區域11的下邊界及上邊界中的氣閥下底面邊界為層動邊界。滑移邊界主要有區域11與區域1，2，7，8的連接邊界，區域7與區域5、區域8與區域6的連接邊界，區域5，6與區域10的連接邊界，區域3的左邊界，區域4的右邊界。區域3與區域5、區域4與區域6及區域11與區域12的連接邊界為靜止的不匹配邊界面。進氣閥上表面區域5及區域6為跟隨氣閥運動區域，這兩區域的劃分是為了適應進氣閥表面的形狀。

為進一步說明，圖 11示出層動區域及滑移邊界的應用示意圖。虛線為靜止的不匹配網格邊界，主要有兩個，一個是氣道區域與氣閥區域的銜接處，一個是燃燒室區域與氣缸區域的銜接處；實線為滑移邊界，主要在氣閥內部區域及外部區域的交界處以及氣閥外部區域與氣缸區域的交界處。圖中區域2為跟隨氣閥運動區域，區域1為適應閥桿運動的層動區域，區域3為適應氣閥上表面運動的層動區域，區域5為適應活塞運動的氣缸層動區域，區域4為適應氣閥底部及活塞運動的層動區域。

假定初始時刻活塞位于上止點，進排氣閥均為關閉狀態，采用上述方法得到的不同曲軸轉角時TBD620柴油機三維網格模型的運動結果見圖12，該款柴油機的基本參數見表2，活塞的運動規律根據轉速、行程、連桿長度及曲軸轉角計算得到。

表2 TBD620柴油機基本參數

可以看到，所采用的分塊滑移動態層法可以很好地適應TBD620柴油機工作過程中氣閥及活塞的運動，整個過程中除運動邊界上的網格需要調整外，其余網格均為靜止狀態。

3 數值應用

3.1 頂蓋驅動方腔流動

頂蓋驅動方腔流動算例是檢驗CFD程序的標準算例，本研究采用該算例驗證所應用不匹配邊界面算法在層流計算時的準確性。采用兩種不同的網格劃分形式，如圖 13所示，其中網格1為均勻網格劃分，不包含不匹配邊界面，網格2采用不均勻網格劃分，包含不匹配邊界面。圖13c為網格2中的實線圍成區域的放大圖，并給出了不匹配邊界處的邊界設置，其中包含兩對不匹配邊界面，分別為邊界1與邊界2，邊界1與邊界3，其中邊界1與邊界2的劃分形式相同，邊界1與邊界3劃分形式不同，存在不對應關系。方腔邊長為1 m，采用的流體介質為水，密度為1 000.0 kg/m3，黏性系數為1.0×10-3kg/(m·s)，上邊界速度為0.000 1 m/s，即計算雷諾數為100，邊界均為無滑移壁面邊界條件。

圖14示出采用兩種網格計算得到的壓力云圖分布對比。可以看到，不匹配邊界面的存在并未對結果產生影響。圖15示出兩種網格X方向及Y方向中心線的速度分布與文獻[12]標準值的對比，兩種網格的計算結果均吻合較好，說明了本研究采用的流場求解方法及不匹配邊界面算法的準確性，也驗證了所開發程序的計算能力。圖16示出采用兩種網格計算得到的壓力殘差對比，可以看到少量的不匹配邊界面并未對殘差的總體收斂趨勢造成影響，但可以看到網格2的收斂還是慢于網格1，說明不匹配邊界面對計算的收斂產生了影響。

3.2 缸內空氣壓縮膨脹過程模擬

對缸內空氣的壓縮膨脹過程進行模擬，取缸內平均壓力值與理論值進行比較。假設壓縮膨脹過程為絕熱過程，則理論壓力可由下式計算得到：

pn=εγp0。

式中：pn為當前壓力；ε為當前容積與初始容積的比值；γ為理想氣體絕熱指數；p0為初始壓力。

圖17、圖18分別示出本算例二維及三維初始網格，缸徑為170mm，高度為200mm，行程為195mm，底面邊界為運動邊界，向頂面作周期為360°的壓縮膨脹運動，在此算例中壓縮比為40，采用動態層網格模型適應底面邊界的運動進行流場計算。圖19示出了計算得到的缸內平均壓力值與理論值的對比，可以看到計算值與理論值吻合較好，從而驗證了動態層模型與現有求解器結合的正確性。

3.3 二維內燃機缸內瞬態流場仿真

采用二維內燃機運動網格模型進行TBD620柴油機瞬態流場的仿真計算。假定初始時刻(0°曲軸轉角)活塞位于存在氣閥重疊角的上止點，配氣定時如圖 20所示。給定進口為速度進口，來流速度為2 m/s，出口為壓力出口，初始時刻缸內空氣密度為1.1 kg/m3，壓力為0.1 MPa，溫度為300 K，計算采用的曲軸轉角步長為0.5°。計算得到不同曲軸轉角下缸內流場的流線分布(見圖21)，可以清楚地看到不同曲軸轉角下缸內流場的分布及燃燒室內氣相漩渦的形成過程，符合真實情況，說明了本研究所開發的非結構內燃機運動網格模型的實用性。

4 結束語

基于自主開發的通用輸運方程求解軟件發展了一種新的基于非結構網格的動態層模型實現算法，并應用分塊滑移動態層思想構建了非結構化內燃機動網格模型。通過頂蓋驅動方腔流動算例和缸內空氣壓縮膨脹過程模擬的算例驗證了所開發滑移網格和層動網格算法下流體流動數值模擬的正確性。最后，以TBD620柴油機模型為例建立了三維內燃機缸內流場運動網格模型，對二維內燃機瞬態流場進行了仿真，計算結果表明本研究提出的網格運動模型可較好地適應內燃機的瞬態流場計算中邊界運動，且可以便捷地加入到現有流動求解器中。

[1] 張志榮, 冉景煜, 張力浦. 內燃機缸內氣體CFD瞬態分析中動態網格劃分技術術[J]. 重慶大學學報,2005,28(11):98-121.

[2] Velghe A, Gillet N, Bohbot J. A high efficiency parallel unstructured solver dedicated to internal combustion engine simulation[J]. Computers & Fluids，2011,45:116-121.

[3] Abouri D, Zellat M, Desoutter G, et al. Advances in combustion modeling in STAR-CD: a new technique for automatic grid and mesh motion generation applied to Diesel combustion and Emissions analysis[C]//19th International Multidimensional Engine User’s Meeting at the SAE Congress. Detroit：SAE，2009.

[4] 蔣炎坤, 李仁旺, 羅馬吉，等. 發動機瞬態數值模擬中氣口網格生成技術研究[J]. 華中科技大學學報,2001,29(5):14-16.

[5] 劉金武, 龔金科, 鐘志華，等. 內燃機缸內復雜空間三維動態網格生成技術[J]. 計算機輔助設計與圖形學學報，2006,18(4):488-492.

[6] 秦文瑾, 解茂昭, 賈明，等. 4氣門直噴式汽油機缸內湍流場多周期循環變動的大渦模擬[J]. Journal of internal combustion engine，2012,30(3):234-240.

[7] 明平劍. 基于非結構化網格氣液兩相流數值方法及并行計算研究與軟件開發[D]. 哈爾濱：哈爾濱工程大學，2008.

[8] 雷國東. 非結構網格FVM在復雜幾何結構的湍流反應流計算中的應用研究[D]. 哈爾濱：哈爾濱工程大學，2008.

[9] Rhie C M, Chow W L. Numerical Study of the Turbulence Flow Past an Airfoil with Trailing Edge Separation[J]. AIAA Journal,1983,21(11):1525-1532.

[10] James M Van Verth, Lars M Bishop. Essential Mat-hematics for Games and Interactive Applications[M]. 2nd Edition.[S.l.]：Elsevier Inc，2008.

[11] 劉永豐. 基于非結構網格的缸內兩相反應流數值模擬方法研究及軟件開發[D]. 哈爾濱：哈爾濱工程大學，2013.

[12] Ghia U, Ghia K N, Shin C T. High Resolutions for Incompressible Flow Using the Navier-Stokes Equations and a Multigrid Method[J]. Journal of Computational Physics，1982,48:387-411.

[編輯： 潘麗麗]

Dynamic Mesh Model Applied to ICE Transient Simulation

SUN Huawen1， YANG Lihong2， MING Pingjian3， ZHANG Wenping3

(1. National Supercomputer Center in Tianjin， Tianjin 300457， China；2. Patent Examination Cooperation Tianjin Center of The Patent Office， Tianjin 300304， China；3. College of Power and Energy Engineering， Harbin Engineering University， Harbin 150001， China)

A new dynamic layer mesh algorithm based on unstructured mesh was introduced and the dynamic mesh model of internal combustion engine(ICE) was built by combining the sliding mesh algorithm based on the sliding dynamic layer. The calculation model of TBD620 diesel engine was further established and all the referred algorithms were realized through the self-developed general transport equation solver. The finite volume method and SIMPLE algorithm for the compressible fluid were utilized in the simulation. Moreover, the sliding mesh model and dynamic layer mesh model were verified through the numerical examples and finally the in-cylinder transient flow field of ICE was simulated. The results show that the introduced dynamic mesh method can realize the transient flow field simulation of ICE.

internal combustion engine(ICE)； sliding mesh； dynamic layer algorithm； numerical simulation

2015-12-14；

2016-03-21

孫華文(1988—)，男，碩士，研究方向為計算流體力學；sun_hua_wen@126.com。

10.3969/j.issn.1001-2222.2016.04.002

TK432

B

1001-2222(2016)04-0007-07