夏軼棟，伍貽兆，呂宏強，宋江勇

(1.南京航空航天大學(xué)航空宇航學(xué)院，江蘇南京 210016;2.中國飛行試驗研究院，陜西西安 710089)

0 引言

間斷有限元法［1－3］在 1973 年由 Lesaint和 Raviart［4－5］首先提出，后經(jīng) Cockburn 和 Shu［6－8］等逐步發(fā)展。在當(dāng)今計算流體力學(xué)領(lǐng)域，間斷有限元法由于集合了傳統(tǒng)有限元法和有限體積法的特點而成為目前科學(xué)計算領(lǐng)域的研究熱點之一［9－11］。在單元內(nèi)部，間斷有限元法跟傳統(tǒng)有限元法一樣用高階多項式來提高精度。而在單元的邊界處，間斷有限元法又可以方便地采用有限體積法中的思路來基于Riemann問題構(gòu)造數(shù)值通量［5，12］，以此來實現(xiàn)逆風(fēng)格式。

高階間斷有限元法應(yīng)用上的瓶頸之一為計算效率問題［13］。由于在相同的網(wǎng)格上相對于有限體積法包含更多未知變量，這種方法需要龐大的存儲空間，且計算時間較長。在串行計算機上，這種局限對于大型問題和高精度求解尤為突出。因此為實現(xiàn)高階間斷有限元法的高效性，引入并行計算勢在必行。

本文根據(jù)間斷有限元法的數(shù)據(jù)結(jié)構(gòu)，基于METIS網(wǎng)格分區(qū)技術(shù)，設(shè)計了并行計算策略，實現(xiàn)了高階間斷有限元法并行計算程序。在非結(jié)構(gòu)網(wǎng)格上對二維Euler方程的亞聲速情況進行了數(shù)值模擬實驗，并對其加速比和并行效率進行分析。在程序設(shè)計中，采用動態(tài)內(nèi)存分配合理使用有限存儲空間。經(jīng)與串行計算結(jié)果對比，本文的高階間斷有限元法并行計算程序得到了較好的加速比和并行效率。這使得采用高階間斷有限元法計算更為復(fù)雜的問題成為可能。

1 間斷有限元法數(shù)值離散

守恒型Euler方程

式(1)兩邊乘測試函數(shù)V，并在計算域內(nèi)積分，運用分部積分后，弱解形式為

式中:Ω為計算域;?Ω為Ω的邊界。

將計算域劃分為網(wǎng)格{e}，其中e表示網(wǎng)格單元。在每個單元內(nèi):

式中φi(X)為基函數(shù)。將式(3、4)代入式(2)，得:

式中?e為單元e的邊界。式(5)必須對于任何單元和任何Vh都滿足。而在每個單元內(nèi)部，Vh是N個基函數(shù)的線性組合，因此式(5)可寫為:

式(6)稱為p階間斷有限元離散［9］，p為式(3)和式(4)中基函數(shù)的最大階數(shù)。

2 分區(qū)并行計算方法

2.1 網(wǎng)格分區(qū)和邊界處理

高階間斷有限元法的數(shù)據(jù)結(jié)構(gòu)特點是，每個網(wǎng)格單元流場信息的計算只涉及其相鄰單元。為實現(xiàn)其并行計算，使用軟件METIS對網(wǎng)格進行分區(qū)，將全局網(wǎng)格劃分為若干局域網(wǎng)格，局域網(wǎng)格之間的數(shù)據(jù)僅與相鄰局域網(wǎng)格分區(qū)邊界上的單元有關(guān)。METIS劃分的局域網(wǎng)格具有相近的單元數(shù)，可以保證多進程并行計算的負(fù)載平衡。

在物面和遠場分別采用無穿透和無反射邊界條件。對于網(wǎng)格分區(qū)所產(chǎn)生的分區(qū)邊界，通過引入虛擬網(wǎng)格單元的方法解決(如圖1所示)。1號局域網(wǎng)格的虛擬單元(用●表示)與2號局域網(wǎng)格分區(qū)邊界上的單元(用■表示)相對應(yīng);2號局域網(wǎng)格的虛擬單元(用●表示)與1號局域網(wǎng)格分區(qū)邊界上的單元(用■表示)相對應(yīng)。因此虛擬單元的所有流場數(shù)據(jù)均來自其對應(yīng)相鄰局域網(wǎng)格分區(qū)邊界上的單元，只需通過計算進程間的數(shù)據(jù)傳遞得到。本文的數(shù)據(jù)傳遞方法基于MPI。

圖1 分區(qū)邊界示意圖Fig.1 Illustration of partition boundary

2.2 并行策略設(shè)計與實現(xiàn)

直接采用有限體積方法在單元邊界處定義數(shù)值通量［9］，即式(6)中的通量函數(shù)F(Uh)·n用數(shù)值通量函數(shù)H(U－，U+，n)來代替，其中 U－和 U+分別表示該單元變量和其相鄰單元的變量，單元Ωe的流場信息的求解只與相鄰單元相關(guān)。本文采用Local Lax-Friedrichs(LLF)格式計算數(shù)值通量。

式(6)最終的離散形式可寫為:

式中:u為向量，包含式(4)中的所有未知變量;M為質(zhì)量矩陣;R(u)為殘值。

對式(7)的求解，本文采用的是牛頓-塊高斯賽德爾法(Newton-Block GS)，即:

式中:C∈［0，1］?！鱱通過求解以下線性方程

獲得?！鱱的形式如下:

上述線性方程式(9)的數(shù)據(jù)結(jié)構(gòu)非常適合并行:△u中向量按單元序號依次排列，單元Ωi的△ui求解只與相鄰單元有關(guān)，?R/?u 是稀疏塊矩陣［14－16］。因此在塊矩陣?R/?u的第i行中，只有第i列塊矩陣和與第i號單元相鄰的單元所對應(yīng)的塊矩陣為非零矩陣。

針對式(9)的數(shù)據(jù)結(jié)構(gòu)，設(shè)計了并行的 Newton-Block GS迭代策略。在每次分區(qū)并行Newton迭代中，首先分配動態(tài)空間，將上一次Block GS迭代后各局域網(wǎng)格分區(qū)邊界上的單元△u新值傳遞給相鄰局域網(wǎng)格的虛擬單元，作為相鄰局域網(wǎng)格分區(qū)邊界單元的本次Block GS迭代的舊值，最后釋放動態(tài)空間。各進程按上述規(guī)則同步傳遞數(shù)據(jù)，并通過迭代直至收斂，進入下一次Newton迭代。

另外，在求解 H(U－，U+，n)前，將本次 Newton迭代時局域網(wǎng)格分區(qū)邊界上U－新值傳遞給其相鄰局域網(wǎng)格的虛擬單元，作為其相鄰局域網(wǎng)格分區(qū)邊界上單元的U+新值。各計算進程在每次Newton迭代后需將所有局域網(wǎng)格上求得的殘值R(u)經(jīng)傳遞后求和，同時判斷是否收斂。圖2為多進程并行計算的牛頓法迭代流程。

3 并行計算數(shù)值結(jié)果與性能分析

3.1 NACA0012翼型繞流問題數(shù)值模擬

并行計算環(huán)境應(yīng)用MPICH2，建立在lenovo R510計算機工作站(雙四核Xeon 5310處理器，內(nèi)存8Gb)，程序采用C++語言編寫。

使用單元數(shù)為800的二維非結(jié)構(gòu)網(wǎng)格。來流馬赫數(shù)Ma和迎角α分別取為0.63和2°。利用上述介紹的并行算法對NACA0012翼型繞流問題進行數(shù)值模擬。首先用串行程序計算，然后將全局網(wǎng)格劃分為2個、4個和8個局域網(wǎng)格分別進行并行計算，并給出數(shù)值結(jié)果。

圖2 牛頓法迭代流程圖Fig.2 Flowchart of Newtown iteration

圖3給出了NACA0012翼型非結(jié)構(gòu)網(wǎng)格的8分區(qū)示意圖。圖中，粗線為局域分區(qū)邊界。局域網(wǎng)格的單元個數(shù)分別為 102，100，101，98，100，99，98，102。

圖3 NACA0012翼型8分區(qū)網(wǎng)格Fig.3 8 partition grids of NACA0012 airfoil

將串行計算結(jié)構(gòu)和8分區(qū)并行計算結(jié)果進行比較。圖4給出了p=4的等馬赫線圖。上圖為串行計算所得，下圖為8分區(qū)并行計算所得。通過比較，流場等馬赫線圖高度吻合。

圖5為翼型上下表面壓力系數(shù)分布。上為p=4的并行計算結(jié)果與串行計算結(jié)果比較，曲線完全重合。下為并行計算p=0～4時得到的翼型上下表面壓力系數(shù)分布，可見高階情況下的精度提高效果非常明顯。

圖4 亞聲速下NACA0012翼型繞流等馬赫線圖(Ma=0.63，α =2°)Fig.4 Mach isolines around NACA0012 airfoil(Ma=0.63，α =2°)

圖5 NACA0012翼型壓力系數(shù)分布(Ma=0.63，α=2°)Fig.5 Distribution around NACA0012 airfoil(Ma=0.63，α =2°)

圖6對并行計算和串行計算的數(shù)值模擬收斂過程進行比較，其中橫坐標(biāo)為Newton迭代步數(shù)，縱坐標(biāo)為殘值R。結(jié)果表明:兩者殘值曲線重合。p=0時以給定的初始解為遠場來流，經(jīng)16步Newton迭代后殘值降到10－10以下，然后以此收斂解作為p=1時的初始解，經(jīng)12步迭代后殘值降到10－10以下。最后當(dāng)p=4時經(jīng)9步迭代后達到收斂。上述比較驗證了高階間斷有限元法并行計算的正確性，同時表明，Newton-Block GS法具有快速收斂和迭代效率高的優(yōu)點。

圖6 NACA0012翼型數(shù)值模擬的收斂過程(Ma=0.63，α =2°)Fig.6 Convergence history of NACA0012 airfoil numerical simulation(Ma=0.63，α =2°)

3.2 并行算法性能分析

衡量并行算法性能的主要指標(biāo)是并行加速比和并行效率。設(shè)Ts為串行程序計算時間，Tp為并行程序計算時間，n為處理器計算進程數(shù)，則并行計算的加速比為S=Ts/Tp，并行效率為E=S/n×100%。

表1給出了采用不同規(guī)模的網(wǎng)格，對上述NACA0012翼型繞流問題進行數(shù)值模擬的并行性能分析。第一行數(shù)據(jù)1、2、4、8表示劃分局域網(wǎng)格的個數(shù)即采用處理器的核數(shù)。

表1 NACA0012翼型數(shù)值模擬的并行性能(Ma=0.63，α =2°)Table 1 Convergence timings，speedup and parallel efficiency for NACA0012 airfoil numerical simulation(Ma=0.63，α =2°)

通過數(shù)據(jù)分析可以得出，當(dāng)處理器核數(shù)一定時，隨著計算規(guī)模的增大，數(shù)據(jù)通信時間TM在整體計算時間中所占比例降低，加速比和并行效率都有所提高。

另外，當(dāng)求解問題的計算規(guī)模一定時，隨著處理器數(shù)核數(shù)的增加，盡管加速比有所提高，但各處理器間的總數(shù)據(jù)通信時間所占總計算時間比例增加，因此并行效率逐步降低。

4 結(jié)論

本文實現(xiàn)了非結(jié)構(gòu)網(wǎng)格上高階間斷有限元法的并行計算。實驗表明，該方法的并行求解具有收斂快、精度高以及大幅節(jié)省計算時間的特點。這使得采用高階間斷有限元法處理更復(fù)雜問題時，在保證高精度解的同時，能夠有效縮短計算時間。

另外，本文設(shè)計的高階間斷有限元法的并行計算策略可直接推廣到Navier-Stokes方程和三維情況中去。這方面的研究工作目前正在進行中。

［1］LIN S Y，CHIN Y S.Discontinuous Galerkin finite element method for Euler and Navier-Stokes equations［J］.AIAA Journal，1993:2016-2023.

［2］BASSI F，RRBAY S.High-order accurate discontinuous finite element solution of the 2D Euler equations［J］.Journal of Computational Physics，1997，138(2):251-285.

［3］RASETARINERA P，HUSSAINI M Y.An efficient implicit discontinuous spectral Galerkin method［J］.Journal of Computational Physics，2001:718-738.

［4］LESAINT P.Sur la resolution des systemes hyperboliques du premier order par des methods delements finis［D］.Unversite Paris 6，1975.

［5］LESAINT P，RAVIART P A.On a finite element method to solve the neutron transport equation［M］.//C.de Boor.Mathematical Aspects of Finite Elements in Partial Differential Equations.Academic Press，New York，1974.

［6］COCKBUTRN B，HOU S，SHU C W.TVD Runge-Kutta local projection discontinuous Galerkin finite element method for conservation lawsⅡ:general framework［J］.Mathematics of Computation，1989:52-411.

［7］COCKBUTRN B，HOU S，SHU C W.TVD Runge-Kutta local projection discontinuous Galerkin finite element method for conservation lawsⅢ:one dimensional systems［J］.Journal of Computational Physics，1989:84-90.

［8］COCKBUTRN B，HOU S，SHU C W.TVD Runge-Kutta local projection discontinuous Galerkin finite element method for conservation laws IV:the multidimensional case［J］.Mathematics of Computation，1990:45-581.

［9］呂宏強，伍貽兆，周春華，等.稀疏非結(jié)構(gòu)網(wǎng)格上的亞聲速流高精度數(shù)值模擬［J］.航空學(xué)報，2009，30(2):200-204.

［10］賀立新，張來平，張涵信，等.任意單元間斷Galerkin有限元計算方法研究［J］.空氣動力學(xué)學(xué)報，2007，25(2):157-162.

［11］賀立新，張來平，張涵信，等.間斷Galerkin有限元和有限體積混合計算方法研究［J］.力學(xué)學(xué)報，2007，39(1):15-22.

［12］COCKBURN B，KAMIADAKIS G E，SHU C W.Discontinuous Galerkin methods:theory，computation and applications［M］.Berlin:Springer，1999.

［13］LUO H，BAUM J D，LOHNER R.A p-multigrid discontinuous Galerkin method for the Euler equations on unstructured grids［J］.Journal of Computational Physics，2006，211(2):767-783.

［14］LU H Q.High-order finite element solution of elastohydrodynamic lubrication problems［D］.Leeds，UK:University of Leeds，2006.

［15］HARTMANN R.Adaptive finite element methods for the compressible Euler equations［D］.Heidelberg，Germany:University of Heidelberg，2002.

［16］LU H Q，BERZINS M，GOODYER C E，et al.High order discontinuous Galerkin method for elastohydrodynamic lubrication line contact problems［J］.Communications in Numerical Methods in Engineering，2005，21(11):643-650.