一種求解汽車外流場問題的可擴展數(shù)值算法

閆爭爭陳榮亮趙宇波蔡小川，2

1（中國科學院深圳先進技術(shù)研究院工程與科學計算研究室 深圳 518055）

2（美國科羅拉多大學博爾德分校計算機系 博爾德 80309）

一種求解汽車外流場問題的可擴展數(shù)值算法

閆爭爭1陳榮亮1趙宇波1蔡小川1，2

1（中國科學院深圳先進技術(shù)研究院工程與科學計算研究室 深圳 518055）

2（美國科羅拉多大學博爾德分校計算機系 博爾德 80309）

受外型復雜、雷諾數(shù)高等因素影響，汽車外流場流動的數(shù)值計算規(guī)模巨大且難以精確求解。發(fā)展高效并行算法以利用超級計算平臺資源來數(shù)值求解外流問題成為該領域的研究熱點。文章提出一種全隱格式的可擴展并行 Newton-Krylov-Schwarz 算法對某真實汽車的外流場流動問題進行計算。通過與風洞試驗以及主流計算流體力學軟件的計算結(jié)果對比驗證了算法的正確性。并行數(shù)值計算結(jié)果顯示，文章的算法在數(shù)千處理器規(guī)模下仍具有很好的并行可擴展性。

三維非定常不可壓縮流動；Navier-Stokes 方程；并行計算；區(qū)域分解算法中

1 引 言

在汽車工業(yè)中，車體外型的空氣動力學分析涉及車輛設計、測試和制造等諸多環(huán)節(jié)，直接決定了車輛的能耗比、安全性乃至美觀性，同時也對乘客的舒適感以及道路周邊噪聲環(huán)境產(chǎn)生重要影響。獲取車體外流場流體的流動數(shù)據(jù)是進行車型氣動分析的前提，工程師只有在速度、壓力等物理參量的基礎上才能通過一系列分析手段得到車型的阻力系數(shù)、氣動力矩等評價車型的直觀數(shù)據(jù)，因此外流場流動問題對汽車空氣動力學研究意義重大。

最初，風洞試驗是獲取外流場數(shù)據(jù)的主要技術(shù)手段，通過全尺寸或比例模型在風洞中的測試結(jié)果來調(diào)整實際車型的幾何尺寸與構(gòu)形。根據(jù)Buchheim 的經(jīng)驗，開發(fā)一款新的車型需要 1000小時的風洞試驗［1］，顯然這種方法需要非常長的設計周期，嚴重影響車型定型與市場推廣。近二十年來，利用計算機來求解流體流動的控制偏微分方程組，并通過得到的流場來研究流動現(xiàn)象的計算流體力學（Computational Fluid Dynamics，CFD）技術(shù)廣泛應用于汽車設計領域，現(xiàn)已成為不可或缺的技術(shù)手段，其應用程度已經(jīng)成為評價工程設計先進水平的主要標志［2］。與風洞試驗相比，CFD 技術(shù)有助于工程師實現(xiàn)數(shù)字化、靈活化以及低成本的理想設計模式。

在流體力學中，雷諾數(shù) Re 作為一個無量綱量用來表征流體運動中慣性力與粘性力之比，其值越大表明流體運動中對流越占優(yōu)，控制方程離散后得到的方程組非線性越強，對求解造成巨大困難。而汽車外流場流動為典型的高雷諾數(shù)流動，一般大于 106，采用未添加任何近似流動模型進行直接數(shù)值求解方法需要極高的時間、空間分辨率。工程應用中普遍采用近似模擬方法將流體運動加以近似或簡化，使之可在高性能工作站（一般為十數(shù)核處理器）平臺上求解。依據(jù)假設狀況，近似模型方法可歸納為如下幾類：

（1）基于統(tǒng)計理論的雷諾平均 Navier-Stokes方程（Reynolds Averaged Navier-Stokes，RANS）方法［3，4］。該方法是當前工程問題中最常用以及應用最廣泛的流動模型，其原理是將流動視為平均流動與脈動流動的疊加，流場中的各參數(shù)也表示為兩者的疊加值。由于引入脈動量乘積等物理量導致方程組的不封閉，因此需添加湍流模型使方程組可解。通過流體力學理論或者實驗數(shù)據(jù)可得到多種不同的模型，如 k-ε 雙方程模型、 k-ω雙方程模型、Reynolds 應力模型等，需根據(jù)實際情況選擇合適的模型進行計算。RANS 方法通過近似使得計算量大為降低，但缺點是模型選擇依靠經(jīng)驗、計算仿真度較低。

（2）基于流動分解技術(shù)的大渦模擬（Large Eddy Simulation，LES）方法［5，6］。大渦數(shù)值模擬來源于對流場中不同尺度的渦的識別與歸類。其原理是將流動分解為大尺度和小尺度兩種運動：大尺度的漩渦具有各向異性特征，是流動的主要表現(xiàn)；小尺度漩渦主要體現(xiàn)能量的耗散情況。在計算過程只使用直接模擬方法計算大尺度渦的運動，小尺度渦則通過近似模型來模擬。這種方法與 RANS 方法相比雖然增加了計算量但卻提高了計算仿真度。

（3）脫體渦模擬（Detached Eddy Simulation，DES）方法［7］。脫體渦模擬結(jié)合了雷諾時均方法與大渦模擬的優(yōu)點，簡單而言即是在附面層附近使用 RANS 方法，在其他區(qū)域使用 LES 方法，其計算量與計算仿真度介于雷諾平均 Navier-Stokes方程與大渦模擬兩者之間。

仿真度與計算能力始終是個不可調(diào)和的矛盾，要想達到可信的仿真度就必須細化網(wǎng)格以精確描述流體的流動特征。CFD 技術(shù)的核心是求解流體控制方程經(jīng)有限元等方法離散后所形成的稀疏方程組。實際計算中，該稀疏方程組的規(guī)模十分巨大且伴隨著強烈的非線性，求解異常困難。近年來，大型并行計算機的出現(xiàn)為計算能力提供了有力的保障，研究適用于求解大規(guī)模非線性稀疏方程組的高可擴展并行算法，并使之與計算機硬件資源匹配，已經(jīng)成為計算機輔助工程領域非常重要的研究方向［8］。

為獲取流場的高分辨率數(shù)值結(jié)果，本文采用未添加流動近似模型的全 Navier-Stokes 方程方法。該方法龐大的計算量對計算機的硬件與算法都有極高的要求。本文基于國家超級計算深圳中心的曙光星云超級計算平臺，通過對大規(guī)模非線性系統(tǒng)進行高效求解器和預處理技術(shù)的研究，提出一種可擴展并行 Newton-Krylov-Schwarz（NKS）算法。非線性方程采用非精確Newton 方法進行求解，在每個 Newton 步，Jacobian 系統(tǒng)通過基于區(qū)域分解的限制加性Schwarz 預條件子處理后，使用以廣義最小殘量（Generalized Minimal RESidual，GMRES）方法為代表的 Krylov 子空間迭代法作為線性求解器進行求解。為驗證算法，本文首先對障礙物繞流標準問題進行研究，所獲得的流場數(shù)值結(jié)果與風洞試驗數(shù)據(jù)及商業(yè)軟件計算結(jié)果進行對比。隨后，作為應用，我們對某真實汽車的外流場問題進行計算。數(shù)值結(jié)果顯示，本文提出的算法計算所得流場數(shù)值相較于采用近似流動模型的 RANS方法所得結(jié)果更為精確。算法的并行測試結(jié)果顯示，本文的算法在擴展至數(shù)千個處理器平臺時仍具有非常良好的可擴展并行性能。

本文的結(jié)構(gòu)如下：第一節(jié)簡要介紹汽車外流場的數(shù)值模擬技術(shù)的現(xiàn)狀與趨勢；第二節(jié)對流體控制方程及本文采取的有限元離散格式進行介紹；NKS 算法描述以及相關(guān)數(shù)值試驗分列于第三、四節(jié)；最后在第五節(jié)做出總結(jié)與展望。

2 流體控制方程及其有限元離散

2.1 流體控制方程

數(shù)值模擬技術(shù)通過計算流體控制方程以獲取流體運動區(qū)域的離散解，汽車周圍流體流動速度一般低于 0.3 Ma，可視為不可壓縮 Newton 流，因此采用三維非定常不可壓 Navier-Stokes 方程描述汽車外流場流動：

以及邊界條件：

2.2 控制方程離散

為了獲取流體流動數(shù)值解，數(shù)值模擬技術(shù)通過將控制偏微分方程按照某種離散格式離散，然后轉(zhuǎn)化為各個計算節(jié)點上的代數(shù)方程組，最后求解得到各節(jié)點上的值來獲取整體流場的近似解，因此離散化與代數(shù)化是數(shù)值計算的前提。對于本文的瞬態(tài)問題，偏微分方程的離散在空間和時間上都有體現(xiàn)，選取合適的時空離散格式對后續(xù)求解算法有重要影響。

空間離散的具體實現(xiàn)形式分為兩方面：在數(shù)學上是對控制方程的空間項采用一定數(shù)值格式使之代數(shù)化；在物理模型上是根據(jù)數(shù)學離散格式的要求進行計算區(qū)域的網(wǎng)格剖分?？紤]到高階的有限元雖然在計算精度上有所提高，但其在離散化后所形成的矩陣中非零元素過多，從而影響算法整體的并行可擴展性。結(jié)合并行計算的特點，我們采用低階的 p1-p1有限元格式對空間進行離散。計算網(wǎng)格主要分為結(jié)構(gòu)化網(wǎng)格與非結(jié)構(gòu)化網(wǎng)格兩種類型。結(jié)構(gòu)化網(wǎng)格對復雜幾何外形適應性較差，因此難以用于車輛外流場等復雜的實際工程問題。而非結(jié)構(gòu)網(wǎng)格中節(jié)點與單元可控性強，可根據(jù)問題對網(wǎng)格密度進行優(yōu)化處理，提高網(wǎng)格的質(zhì)量，目前已成為諸多工程問題普遍采用的網(wǎng)格類型。

有限維試函數(shù)及全函數(shù)空間可以表示為：

多數(shù)空間離散時通常按照物理量排列，此類按照各相關(guān)變量（壓力、速度）的排序方式會在結(jié)構(gòu)上影響最終求解的 Jacobian 矩陣，呈現(xiàn)出鞍點結(jié)構(gòu)，不利于并行計算。本文中，空間離散按照單元順序排列方法構(gòu)造。此方法作為預條件子，在后續(xù)的區(qū)域分解過程中將原有整體的鞍點結(jié)構(gòu)分散在局部的子矩陣中，減少子區(qū)域間共享的變量數(shù)量，從而減少通信量，保證了并行效率。

時間項方面，本文采用了全隱格式進行離散。由于全隱格式減輕甚至消除了 CFL 條件對時間步長的限制，在進行大規(guī)模數(shù)值計算時具有一定優(yōu)勢。具體地，針對上述不可壓縮 Navier-Stokes 方程，選取隱式向后 Euler 有限差分格式對時間項離散，即將控制方程（1）離散為如下非線性代數(shù)方程：

其中，Δt為時間步長；S(xn)為在第 n 時間步經(jīng)空間離散后形成的半離散系統(tǒng)。至此，問題化解為在每個時間步求解非線性方程組：

3 并行 Newton-Krylov-Schwarz 算法

在每個時間步內(nèi)，經(jīng)有限元離散后形成的系統(tǒng)（3）規(guī)模十分巨大，且因?qū)α髡純?yōu)而具有很強的非線性。求解此類大規(guī)模的病態(tài)問題時，一般均需借助大型并行計算機及相應的可擴展性算法。實際上，多數(shù)基于偏微分方程的工程應用與科學計算問題，如污染擴散模擬、數(shù)值天氣預報等，最終都會形成稀疏的線性或非線性方程組。由于在現(xiàn)實生活中這些問題多以非線性的方式存在且往往規(guī)模巨大，因此大規(guī)模稀疏非線性方程組的求解方法成為數(shù)值計算領域的研究熱點。非線性方程組的求解通常包含線性化、線性求解器構(gòu)造以及預處理等多個步驟，本文據(jù)此進行相關(guān)非線性方程組求解算法的研究。

目前非線性方程組線性化方法中最常用且有效的是 Newton 類方法。Newton 類方法具有結(jié)構(gòu)簡單、自校正等諸多優(yōu)點，已廣泛應用于工程與科學計算各個領域。Newton 類方法包括精確 Newton方法、擬 Newton 方法以及非精確 Newton 方法三種主要類型。精確 Newton 方法需要對 Jacobian 矩陣進行精確求解，因此該方法收斂速度最快，但當問題規(guī)模較大時，整體的計算量非常巨大，計算成本高且計算效果差。Newton 方法一般只要進行若干次迭代即可達到收斂精度，為了保留此優(yōu)點，雖然可以構(gòu)造更高收斂階的方法，但許多數(shù)值算法還是在 Newton 方法基礎上改進 Jacobian 矩陣處理方式以提高求解的性能。其中，擬 Newton 方法通過將 Jacobian 矩陣近似為易于求解的其他矩陣，從而達到減少整體計算量的目的。某些情況下，擬 Newton 方法中近似矩陣的求解成本依然較高。而非精確 Newton 方法雖然需要求解原始的 Jacobian 矩陣，但是不需要精確求解，而是通過調(diào)節(jié)求解精度來達到計算量與收斂速度的平衡。

求解大規(guī)模稀疏線性方程組分為直接法與迭代法兩類。直接法按照元素或行處理矩陣，以矩陣分解或消去為手段，具有高精度、高穩(wěn)定性等優(yōu)點，在求解規(guī)模較小的問題時能達到快速求解的目的。然而，當計算規(guī)模增大時，如浮點運算次數(shù)、所需存儲量等計算成本開銷大幅增加。隨著三維有限元網(wǎng)格高分辨率的要求，自由度以千萬為量級的問題愈加平常。因此在計算超大規(guī)模系統(tǒng)時，直接方法已不現(xiàn)實。而迭代法能夠以較小的計算機存儲及開銷代價求解大規(guī)模的問題，且相較而言其在并行求解中比直接法具有更易控制的優(yōu)點，因此迭代法得到廣泛的應用，已經(jīng)成為求解大規(guī)模非線性方程組的不二選擇。目前 Krylov 子空間迭代方法具有占用內(nèi)存小、求解速度快等優(yōu)點而被廣泛作為線性解法器。

無論是在線性還是在非線性迭代法中，預條件子都起到加速收斂的作用，在數(shù)值算法構(gòu)造中占據(jù)重要地位。其作用是以較低的額外代價換取迭代次數(shù)的大幅降低和求解時間的有效減少，預條件子的構(gòu)造往往依賴于所研究的物理問題［10］，對預條件子的改進是提高整個數(shù)值算法性能的關(guān)鍵。Dryja 等［11］的研究表明，區(qū)域分解方法既可以作為迭代法也更適合于作為其他高效迭代求解器的預條件子。同時區(qū)域分解預條件子由于其“分而治之”的特點，尤其適合于當前流行的分布式并行計算機環(huán)境。

基于上述方法，本文提出了求解汽車外流場問題的 NKS 算法。使用非精確 Newton 算法作為非線性求解器，在每個非精確 Newton 步里使用由基于區(qū)域分解方法的 Schwarz 預處理算子來加速Krylov 算法求解 Jacobian 方程。其詳細步驟如下：

（1）使用前一個時間步的解作為初始值

（2）對 k=0，1，… 執(zhí)行以下步驟，直至收斂：

① 構(gòu)造 Jacobian 矩陣；

在本研究中，我們使用該算法求解 2.2 節(jié)中控制方程離散后形成的非線性系統(tǒng)（3）。

本章介紹的算法基于阿貢國家實驗室開發(fā)的開源可移植并行軟件包 PETSc（Parallel Extensible Toolkits for Scientific Computing）［13］實現(xiàn)。PETS基于基礎線性代數(shù)子程序庫 BLAS、線性代數(shù)包LAPACK 以及消息傳遞接口庫 MPI 等構(gòu)成線性求解器、非線性求解器以及時間步進積分器三個主要的組件，為用戶提供了包括求解大規(guī)模稀疏線性方程組的多種 Krylov 子空間方法在內(nèi)的豐富的算法以及預條件子，在高性能計算平臺上具有強大的偏微分方程數(shù)值求解能力。目前 PETSc已經(jīng)成功應用于優(yōu)化問題、生物醫(yī)學、計算流體動力學等多類工程與科學計算問題。

本文的數(shù)值算法的應用算例在國家超級計算深圳中心的曙光星云超級計算機上展開。該計算平臺的科學計算分區(qū)擁有 640 個計算節(jié)點，每個節(jié)點配置雙路六核心 Intel 5650@2.76 GHz 處理器與 24 GB 內(nèi)存，共計 7680 個 CPU 處理核心，網(wǎng)絡互聯(lián)方式為 InfiniBand4 萬兆網(wǎng)絡，操作系統(tǒng)為 Linux。

4 數(shù)值算例

在本節(jié)，我們采用 NKS 算法所獲取的結(jié)果與風洞試驗以及計算流體力學程序包 CFX 的計算結(jié)果進行對比來驗證算法的正確性。

與大多數(shù)計算流體力學程序包采用有限體積方法不同，CFX 采用了基于有限單元的有限體積方法，在保證單元的守恒特性的基礎上結(jié)合了有限元方法的精確特性。CFX 提供了隱格式的求解技術(shù)，改進了“壓力項假設-求解-修正壓力項”的傳統(tǒng)迭代算法，增加了計算的穩(wěn)定性與收斂速度。為使實驗結(jié)果具有可比性且有對比價值，NKS 與 CFX 進行數(shù)值計算時均采用同樣的幾何模型、流體材料屬性、計算區(qū)域、有限元網(wǎng)格以及初始和邊界條件。

4.1 障礙物繞流標準問題

為了支持和開展基于高性能計算機的流體數(shù)值模擬研究，德國研究基金通過優(yōu)先支持項目資助了多種典型流動的數(shù)值模擬研究工作。通過開展大量的風洞試驗獲取相關(guān)流場信息，為驗證數(shù)值模擬算法提供了豐富有效的對比數(shù)據(jù)，障礙物繞流問題便是其中之一。車輛的外部流場屬于典型的繞流問題，我們采用 NKS 算法對其中兩類繞流標準問題進行計算，所獲取的結(jié)果與風洞試驗以及 CFX 的計算結(jié)果進行對比來驗證算法的正確性。

4.1.1 障礙物標準繞流問題模型

本小節(jié)的數(shù)值算例為圓柱體與長方體外部繞流［14］，如圖 1 所示。

圖1 障礙物標準繞流問題模型及計算區(qū)域尺寸Fig. 1 Dimensions of the obstacles and the computational domain

本例中 Newton-Krylov-Schwarz 算法各參數(shù)設置如下：線性與非線性相對誤差分別設置為1.0×10-4和 1.0×10-6；CFX 中，采用二階向后Euler 格式對時間項進行離散。殘差收斂條件設置為均方根誤差設為 1.0×10-6，在 16 核處理器條件下進行計算。

4.1.2 流場數(shù)值結(jié)果

在外流場分析中，氣動力及其系數(shù)如阻力系數(shù) Cd、升力系數(shù) Cl等是衡量氣動性能的最重要的參數(shù)，氣動力的計算公式如下：

阻力 Fd：

升力 ：

其中，S 為障礙物表面積；nx、ny分別為 S 在x、y 方向的法向分量；vt為速度沿著 S 的切向向量方向的分量。

與阻力及升力相對應的阻力系數(shù) Cd與升力系數(shù) Cl定義如下：

可見標準繞流問題的雷諾數(shù)較低。

為了對比其他流場信息，除了氣動系數(shù)外，我們還在障礙物迎風面及背風面設置兩個觀測點 a（xa，ya，za）=（0.45，0.20，0.205）與 b（xb，yb，zb） =（0.55，0.20，0.205），使用兩點間壓力差= P（xa，ya，za，t）-P（xb，yb，zb，t）來觀測與對比流場的數(shù)值結(jié)果。對于兩個障礙物算例，本文對比了其穩(wěn)態(tài)與瞬態(tài)兩種情況，列于表 1。

表1 障礙物繞流測試算例Table 1 Test cases of flow around the obstacles

NKS、CFX 以及風洞（Wind Tunnel，WT）試驗對比結(jié)果列于表 2。從對比結(jié)果可以看出，本文提出的算法與商業(yè)軟件計算結(jié)果接近，同樣都在風洞試驗所得數(shù)據(jù)誤差范圍之內(nèi)，證明了數(shù)值方案以及算法的正確性。

表2 針對障礙物擾流問題 NKS 算法、CFX 以及風洞試驗結(jié)果對比Tables 2 Comparison of the NKS algorithm， CFX and wind tunnel for the benchmark problems

4.1.3 并行可擴展性結(jié)果

本小節(jié)給出基于 NKS 算法計算圓柱形障礙物標準繞流問題的并行可擴展性結(jié)果。兩套不同規(guī)模的網(wǎng)格，分別約含 5.1×106個四面體單元（DOF =3.6×106）和 1.3×107個四面體單元（DOF= 9.1×106），被用于可擴展性的測試。取固定時間步長，NKS 算法中線性與非線性相對誤差分別設置為 1.0×10-6和 1.0×10-12。為減小計算量，取前 20 時間步數(shù)值結(jié)果的平均值。

表 3 列出了不同處理器條件下計算的可擴展并行計算數(shù)值結(jié)果。其中，np為處理器個數(shù)；Newton、GMRES 和 Time 分別表示為每個時間步內(nèi)非線性迭代次數(shù)、線性迭代次數(shù)以和計算時間（單位為秒）。

表3 障礙物繞流問題中 NKS 算法的并行性能Table 3 Parallel performance of the NKS algorithm for the benchmark problem

為了更加清晰地展示可擴展性性能，本文將并行計算數(shù)值結(jié)果以加速比的形式列圖 2 中。加速比通常用來衡量并行系統(tǒng)或程序并行化的性能和效果。在本文中加速比是指固定計算規(guī)模，以最低核數(shù)所需計算時間基準，增加并行處理器個數(shù)后計算時間與基準時間之比的倒數(shù)。 理想加速比是指線性的加速情況。

圖2 障礙物繞流模擬中每時間步的平均計算時間及加速比Fig. 2 The average compute time per time step and the speedup for benchmark problem

并行計算數(shù)值結(jié)果顯示非精確 Newton 方法使非線性迭代次數(shù)非常少，非線性迭代次數(shù)幾乎與處理器個數(shù)無關(guān)。線性迭代次數(shù)隨核數(shù)的增加而略微增加。計算時間與加速比等數(shù)值結(jié)果顯示，針對低雷諾數(shù)問題，在 DOF=3.6×106的算例中，當核數(shù)小于 1024 時，并行效率在 65%以上；當核數(shù)達到 1024 時，計算效率會有所下降。這并不是算法使然，而是由計算規(guī)模不足造成的。當進行分區(qū)計算時，如果子區(qū)域的規(guī)模較小，那么在整體求解時間中，子區(qū)域間通信時間所占的比例將會提高，導致計算效率降低。整體而言，本文提出的并行 NKS 算法具有非常好的可擴展性，當處理器個數(shù)達到 2048 時并行效率仍在 50% 左右。

4.2 汽車外流場問題

4.2.1 汽車幾何模型與計算區(qū)域

在分析車型空氣動力學性能前，首先需要使用數(shù)字化建模技術(shù)對原始設計車型進行幾何建模。本文采用三維 CAD 軟件依據(jù)某真實車型全尺寸數(shù)據(jù)進行模型構(gòu)造。如圖 3 所示，整車長 L =4.6 m，車寬 W=2.0 m，車高 H=1.4 m，軸距B=1.6 m。

圖3 汽車三維 CAD 模型及計算區(qū)域示意圖Fig. 3 The CAD model of car and the computational domain

圖4 汽車外流計算網(wǎng)格示意圖Fig. 4 A view of the inside of the meshed domain

流體取 25℃ 標準大氣壓下的空氣，密度 ρ= 1.185 kg/m3，動力粘度系數(shù) μ=1.831×10-5kg/ms。設進口速度 Vin=30 t，時間步長=0.01 s，計算1 s 內(nèi)流場的變化情況。算法中的線性和非線性求解器的相對收斂誤差分別取為 10-4和 10-6。分別取 =30 m/s 為特征速度，車高 H=1.41 m為特征長度，計算雷諾數(shù)

4.2.2 流場數(shù)值結(jié)果

以阻力系數(shù) Cd為代表的空氣動力學系數(shù)是影響車型設計的最主要參數(shù)，設計人員往往據(jù)此來評判車型的能耗比以及氣動安全性能。本文中首先對 0～1 s 時間內(nèi) NKS 算法與 CFX 計算得到阻力系數(shù) Cd進行比較，瞬態(tài)數(shù)值結(jié)果列于圖 5。阻力系數(shù)對比圖顯示，兩種算法計算得到的阻力系數(shù)基本一致。

圖5 汽車阻力系數(shù)數(shù)值結(jié)果對比Fig. 5 The comparison of drag coefficient

隨后，我們?nèi)?t=1 s 時刻兩種算法的流場結(jié)果進行比較。首先我們比較了車身表面靜壓（單位為帕），如圖 6 所示。

圖6 汽車車體前身壓力分布對比云圖（左：CFX，右：NKS）Fig. 6 The comparison of pressure contours（left：CFX， right： NKS）

靜壓的大小與分布與汽車所受的氣動力直接相關(guān)，NKS 與 CFX 的計算結(jié)果均顯示車身表面靜壓呈對稱分布。靜壓高值分布于車體前身、前擋風玻璃、后視鏡以及車輪的迎風面。車前身其余部分壓力稍高于大氣壓，車體頂部、車窗與車體邊緣處靜壓在較小。

為了更加直觀地觀測汽車外流場流體流動情況，本文也給出車體對稱面 z=0 m 上的二維流線分布。圖 7 描述了 NKS 與 CFX 計算的汽車對稱面上的渦流結(jié)構(gòu)。其中，NKS 所計算結(jié)果中，在車體前部地面、發(fā)動機罩與前擋風玻璃交界處，以及尾部地面上有較明顯的渦流，車尾區(qū)域亦有數(shù)個尺度較小的渦流。而 CFX 計算結(jié)果中卻無法清晰顯示，原因是 CFX 采用了基于統(tǒng)計理論的 RANS 流動模型，相同網(wǎng)格尺寸下其計算精度會有影響。兩種算法的結(jié)果比較顯示，未添加近似流動模型的 NKS 算法可以計算得出更詳細的流場信息。

圖7 車體對稱面二維流線圖（上：CFX，下：NKS）Fig. 7 Surface streamline patterns on the symmetry plane for the car （upper： CFX， lower： NKS）

4.2.3 并行可擴展性結(jié)果

本小節(jié)給出基于 NKS 算法求解汽車外流場的并行可擴展性數(shù)值結(jié)果，實驗所取網(wǎng)格自由度約為 1.03×107，Newton 非線性迭代與 Krylov 子空間方法線性迭代收斂誤差分別設置為 10-6和10-12，數(shù)值結(jié)果取前 20 個時間步的平均值。

表4 列出了不同處理器條件下計算的可擴展并行計算數(shù)值結(jié)果。為了更加清晰地展示可擴展性性能，本文將并行數(shù)值結(jié)果以加速比的形式列于圖 8 中。

表4 汽車外流場模擬中 NKS 算法的并行性能Table 4 Parallel performance of the NKS algorithm for the simulation of flow around car

圖8 每時間步平均計算時間與加速比Fig. 8 The average compute time per time step and the speedup for the simulation of flow around car

從并行測試結(jié)果可以看到，由于采用了前一時間步的結(jié)果作為初始值，每個時間步的非線性迭代次數(shù)非常少，非線性迭代次數(shù)幾乎與處理器個數(shù)無關(guān)。隨著計算核數(shù)的增加，線性迭代次數(shù)略微增加。計算時間等并行計算數(shù)值結(jié)果顯示，本文提出的算法具有非常好的可擴展性，處理器個數(shù)達到 2048 時并行效率仍在 50% 左右。

5 結(jié) 論

復雜幾何造型、高雷諾數(shù)導致汽車外流場的精確數(shù)值求解因計算規(guī)模巨大、非線性強而極難展開，是目前工程計算領域中非常具有挑戰(zhàn)性的問題。本文提出一種并行的 Newton-Krylov-Schwarz 算法來求解控制方程有限元離散后形成的大規(guī)模稀疏非線性方程組。基于該算法，本文首先對較低雷諾數(shù)的標準繞流問題進行求解，流場結(jié)果與風洞試驗數(shù)據(jù)以及 CFX 計算結(jié)果進行對比，驗證了算法的準確性。隨后作為應用，對真實尺寸的汽車外流場流體流動進行數(shù)值模擬。流場數(shù)值結(jié)果以及并行計算數(shù)值結(jié)果表明，本文提出的算法在數(shù)千核處理器條件下仍具有良好的并行可擴展性能。

在本文提出的并行 Newton-Krylov-Schwarz算法中，Newton 型方法嚴重依賴于初值，選取合適的初值對整體計算的迭代次數(shù)與計算時間等并行性能有重要影響［17］。為加快收斂，多水平網(wǎng)格算法是可選擇的算法之一，其原理是先在較粗網(wǎng)格獲取數(shù)值解，再將其插值到較細網(wǎng)格上，作為細網(wǎng)格計算的初值。該算法能夠有效地降低非精確 Newton 方法的迭代次數(shù)以及各步中的線性迭代次數(shù)［18］。因此，下一步研究計劃是結(jié)合車輛外流場的流動特點，將多水平網(wǎng)格技術(shù)應用于大雷諾數(shù)非定常不可壓縮 Navier-Stokes 方程的求解中，獲取更好的并行性能。此外，算法的很多細節(jié)，如網(wǎng)格并行分區(qū)與管理、通信問題等仍有待于進一步的深入研究。

［1］ 張英朝. 汽車空氣動力學數(shù)值模擬技術(shù) ［M］. 北京: 北京大學出版社， 2011.

［2］ 閻超， 都彥喆. CFD 技術(shù)及其在大飛機研制中的應用 ［J］. 航空制造技術(shù)， 2008， 14: 42-44.

［3］ Guilmineau E. Computational study of flow around a simplified car body ［J］. Journal of Wind Engineering and Industrial Aerodynamics， 2008，96（6）: 1207-1217.

［4］ Han T. Computational analysis of three-dimensional turbulent flow over a bluff body in ground proximity ［J］. AIAA Journal， 1989， 7: 1213-1219.

［5］ Krajnovi? S， Davidson L. Flow around a simplified car， part 1: large eddy simulation ［J］. Journal of Fluids Engineering， 2005， 127（5）: 907-918.

［6］ Howard RJA， Pourquie M. Large eddy simulation of an Ahmed reference model ［J］. Journal of Turbulence， 2002， 3（5）: 1-18.

［7］ Guilmineau E， Chikhaoui O， Deng GB， et al. Cross wind effects on a simplified car model by a DES approach ［J］. Computers & Fluids， 2013， 78: 29-40.

［8］ 邱劍， 顧兆林. 利用高階分區(qū)并行算法實現(xiàn)直接數(shù)值模擬 ［J］. 計算力學學報， 2008， 25（1）: 20-24.

［9］ Franca LP， Frey SL. Stabilized finite element method: II. The incompressible Navier-Stokes equations ［J］. Computer Methods in Applied Mechanics and Engineering， 1992， 99: 209-233.

［10］ Knoll D， Keyes DE. Jacobian-free Newton-Krylov methods: a survey of approaches and applications ［J］. Journal of Computational Physics， 2004， 193: 357-397.

［11］ Dryja M， Widlund OB. Towards a Unified Theory of Domain Decomposition Algorithms for Elliptic Problems ［M］. New York: New York University，Courant Institute of Mathematical Sciences，Division of Computer Science， 1989.

［12］ Cai XC， Sarkis M. A restricted additive Schwarz preconditioner for general sparse linear systems ［J］. SIAM Journal on Scientific Computing， 1999，21（2）: 792-797.

［13］ Balay S， Buschelman K， Gropp WD， et al. PETSc Users Manual ［Z］. Argonne National Laboratory，2012.

［14］ Sch?fer M， Turek S. Benchmark Computations of Laminar Flow Around a Cylinder ［M］. Vieweg+ Teubner Verlag， Springer， 1996.

［15］ 王夫亮. 側(cè)風作用下的汽車氣動特性研究 ［D］.吉林: 吉林大學， 2009.

［16］ Karypis G， Aggarwal R， Schloegel K， et al. ParMETIS home page ［OL］. ［2014-7-9］. http:// glaros.dtc.umn.edu/gkhome/metis/parmetis/ overview.

［17］ Leader JJ. Numerical Analysis and Scientific Computation ［M］. Boston: Pearson Addison Wesley， 2004.

［18］ Smith BF， Bj?rstad PE， Gropp WD. Domain Decomposition: Parallel Multilevel Methods for Elliptic Partial Differential Equations ［M］. New York: Cambridge University Press， 1996.

A Scalable Numerical Method for Simulating the External Flows around Cars

YAN Zhengzheng1CHEN Rongliang1ZHAO Yubo1CAI Xiaochuan1，2

1（ Research Centre for Science and Engineering Computation， Shenzhen Institutes of Advanced Technology，Chinese Academy of Sciences， Shenzhen 518055， China ）

2（ Department of Computer Science， University of Colorado Boulder， Boulder， CO 80309， USA ）

The high-fidelity numerical simulation of unsteady flows around cars are a very challenging computational problem because of the huge computation caused by high Reynolds number and complex geometry. In this paper， we presented a scalable parallel Newton-Krylov-Schwarz based fully implicit algorithm for the full unsteady incompressible flows around cars and compared our results with results obtained from commercial CFD software. Our algorithm shows very good parallel scalability on a supercomputer with over two thousand processors.

three-dimensional unsteady incompressible flows； Navier-Stokes equations； parallel computing； domain decomposition method

圖分類號 O 246A

2014-04-23

2014-07-09

中國科學院知識創(chuàng)新工程重要方向項目（KJCX2-EW-L01）；廣東省科技計劃國際合作項目（2011B050400037）

閆爭爭（通訊作者），博士研究生，研究方向為可擴展并行計算，計算流體力學，E-mail：zz.yan@siat.ac.cn；陳榮亮，助理研究員，研究方向為流體優(yōu)化算法；趙宇波，正研級高級工程師，研究方向為計算流體動力學數(shù)值算法；蔡小川，研究員，研究方向為偏微分方程數(shù)值算法與區(qū)域分解方法。