3.2 數值計算方法

為了數值求解控制偏微分方程（Partial Differential Equations，PDEs），引入了偏微分方程的近似方程。這些近似將偏導數轉化為有限差分式，將偏微分方程改寫為代數方程。隨后，就可以利用計算機在目標區域內的離散點處求解這些近似后的代數方程。有數量眾多的數值方法可供選擇，它們在計算精度、效率和穩定性方面各具特點。在本書中，我們不打算涵蓋關于CFD數值方法的大量內容，而只是介紹常用的內燃機模擬程序KIVA和CONVERGE中數值處理方法的主要特點。對CFD數值方法系統化的描述可以在許多CFD書籍中找到，例如張德良的書[6]。

在KIVA程序中[5，7-9]，氣相求解過程是基于一種被稱為任意拉格朗日-歐拉（Arbitrary Lagrangian-Eulerian，ALE）的有限容積法[10，11]。空間差分是在具有任意六面體的有限差分網格上形成的，該網格將計算區域細分為若干個六面體小單元。使用任意形狀網格對內燃機模擬來說具有顯著優點，因為任意網格可以順應曲線邊界，并且可以隨著燃燒室幾何形狀的移動而變化。在這種方法中，除速度外，所有氣體流動物理量都存儲在單元中。氣體速度存儲在以單元頂點為中心的動量單元中。因此，速度位于頂點，而其它量位于單元中心。

KIVA程序中的時間差分格式在很大程度上是隱式的。隱式差分用于所有擴散項和與壓力波傳播相關的項。耦合隱式方程組的求解采用壓力連接方程的SIMPLE半隱式方法（Semi-Implicit Method for Pressure-Linked Equations，SIM-PLE）[12]。瞬時解以一系列被稱為周期或時間步長的有限時間增量函數進行。在每個時間步長中，因變量的值是根據上一個時間步長中的值來計算的。在ALE方法中，每個時間步長被分為兩相：拉格朗日相和再分區相。在拉格朗日相中，頂點隨著流體速度的變化而移動，并且沒有穿過單元邊界的對流輸運。在再分區階段，流場被凍結，頂點被移動到用戶指定的新位置，流場被重新映射或再分區到新的計算網格上。這種再映射是通過跨計算單元邊界的變量輸運來完成的，這些計算單元被認為是相對于流場移動的。準二階迎風QSOU（Quasi-Second-Order Upwind）方法[5]用于計算再分區相的對流輸運。

KIVA中時間步長是根據準確性而非穩定性標準計算得出的。此特點與某些商用CFD軟件不同。后者以穩定性為重，以避免客戶因計算崩潰而產生抱怨。KIVA程序中使用了一種自動時間步長控制方法，該方法要求時間步長Δt滿足相關條件。Δt的第一個精度條件是：

式中，f_a的典型值是0.5；Δx是一個單元的平均尺寸。Δt的第二個精度條件是：

式中，；λ是應變率張量的特征值。此條件限制了由于網格運動而導致的單元變形量。

Δt的另外兩個精度條件，是考慮到化學放熱以及與噴霧的質量和能量交換，從而可以精確地耦合氣體流場和源項的需要。它們由以下公式表示：

式中，f_ch和f_sp的典型值為0.1。此外，時間步長Δt_gr用于限制時間步長的增長量。其表達式為

柯朗穩定性條件還用于限制對流輸運的時間步長Δt_c。該條件由下式給出：

式中，V_i是計算單元i的體積；δV_i是該單元計算的通量體積；f_con的典型取值為0.2，以減少時間步長來保證計算精度。

然后，第n+1個計算周期的主時間步長為

式中，時間步長Δt_mx是計算初始輸入的最大時間步長；Δt_mxca是基于計算初始輸入的最大曲軸轉角的最大時間步長。

利用基于蒙特卡羅方法和離散粒子方法的隨機粒子方法求解噴霧動力學方程。該方法已被證實非常有效和準確。在離散粒子方法中，連續液滴的概率分布函數f近似為離散分布函數f′：

每個粒子p由N_p個具有相同位置x_p、速度、尺寸r_p、溫度和振蕩參數y_p和的液滴組成。粒子和液滴的軌跡重合，粒子在其所在的計算單元中與氣體交換質量、動量和能量。

差分方程可以用來近似控制方程。然后使用上述數值方法計算差分方程。差分方程的計算過程需要一系列計算步驟。計算程序的細節在這里不過多展開，可以在文獻[5]中找到。

為了求解有限差分方程，必須給出計算域內的一組網格點及其邊界。這種網格系統的創建被稱為網格生成。由于內燃機復雜的邊界幾何形狀（如曲面的活塞形狀、傾斜的氣門），網格生成一直是內燃機數值模擬的瓶頸。直到20世紀90年代中期，才報道了針對直噴柴油機[13]和直噴汽油機[14]的具有真實氣門運動的進氣流動的開創性模擬結果。

在KIVA3[7]中，使用塊狀結構網格，通過間接尋址定義連通性。由于不再需要維持大量的失去作用的計算網格，使其可以高效地對復雜幾何結構進行模擬。KIVA3V[8]增加了可對氣缸蓋上具有任意數量的垂直或傾斜氣門進行模擬的功能。它將氣門視為實體，通過所謂的移動技術使其在網格中移動。

雖然網格生成技術的進步使得數值模擬一個完整的內燃機進氣和排氣過程成為可能，但是網格生成仍然是耗時和煩瑣的。福特研究實驗室開發了一種先進的方法來快速生成帶移動氣門的發動機網格，將某新型發動機的網格生成時間從3～4周縮短到了2～3天[15]。

Senecal及其同事[16]在內燃機網格生成方面進行了突破性創新，他們創建了商業CFD軟件CONVERGE并不斷發布更新版本[17]。由于CONVERGE已在全球范圍內獲得廣泛的應用，我們在此簡要介紹其網格生成和數值模擬方法。

CONVERGE使用了一種改進的笛卡兒單元切割方法，該方法不再需要計算網格隨目標幾何體形狀改變而改變，同時仍然可以正確地表示幾何體。這種方法有兩個顯著的優點。首先，選擇網格的類型是為了計算效率，而不是幾何形狀。因而可以使用簡單的正交網格，這大大簡化了程序的數值計算量。第二，網格生成時間和復雜度大大降低，因為復雜的幾何體只需要映射到底層的正交網格上。用戶只需提供一個包含幾何體表面的文件，該表面由三角形曲面構成封閉表面。這個文件很容易從大多數CAD軟件包中以STL格式輸出。

CONVERGE運行時在內部生成實際網格，并使用提供的三角表面去切割與曲面相交的單元。任何單元都可以被任意數量的表面三角形切割，從而使網格能夠精確地表示邊界。使用簡單正交網格允許選擇多種網格處理選項，以提高模擬的精度和速度。網格處理選項包括細化整體網格分辨率、在模擬過程中從粗網格到細網格的網格尺寸映射、由用戶在關鍵流動區域加密網格，以及基于預定義亞網格尺度的自適應加密網格。CONVERGE在網格生成的這些特點使其在內燃機數值模擬中具有獨特的優勢。

在CONVERGE中，采用有限容積法求解守恒方程。使用Rhie和Chow算法[18]將所有計算變量設置在計算單元的中心，并維持變量的同位設置，消除不需要的網格排列。動量方程通過Issa提出的PISO拆分運算式的隱式壓力算法（Pressure Implicit with Splitting of Operators）[19]求解。在CONVERGE中實現的PI-SO算法首先經過一個預估過程求解動量方程。完成預估過程后，導出并求解壓力方程，接下來對動量方程進行修正。這個修正和解析過程可根據需要迭代計算多次，直到達到精度要求。在完成第一次預估和修正后，依次求解其他輸運方程。

CONVERGE還具有多重網格壓力求解器和時間步長控制算法。它可以在共享和分布式內存計算機上并行運行。更詳細的數值方法參見CONVERGE手冊[17]。