夏軼棟，伍貽兆，呂宏強(qiáng)，宋江勇

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

0 引言

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

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

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

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

守恒型Euler方程

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

式中:Ω為計(jì)算域;?Ω為Ω的邊界。

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

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

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

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

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

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

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

在物面和遠(yuǎn)場分別采用無穿透和無反射邊界條件。對于網(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ū)邊界上的單元，只需通過計(jì)算進(jìn)程間的數(shù)據(jù)傳遞得到。本文的數(shù)據(jù)傳遞方法基于MPI。

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

2.2 并行策略設(shè)計(jì)與實(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)格式計(jì)算數(shù)值通量。

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

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

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

式中:C∈［0，1］。△u通過求解以下線性方程

獲得。△u的形式如下:

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

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

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

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

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

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

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

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

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

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

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

圖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為串行程序計(jì)算時(shí)間，Tp為并行程序計(jì)算時(shí)間，n為處理器計(jì)算進(jìn)程數(shù)，則并行計(jì)算的加速比為S=Ts/Tp，并行效率為E=S/n×100%。

表1給出了采用不同規(guī)模的網(wǎng)格，對上述NACA0012翼型繞流問題進(jìn)行數(shù)值模擬的并行性能分析。第一行數(shù)據(jù)1、2、4、8表示劃分局域網(wǎng)格的個(gè)數(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ù)一定時(shí)，隨著計(jì)算規(guī)模的增大，數(shù)據(jù)通信時(shí)間TM在整體計(jì)算時(shí)間中所占比例降低，加速比和并行效率都有所提高。

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

4 結(jié)論

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

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

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

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

［11］賀立新，張來平，張涵信，等.間斷Galerkin有限元和有限體積混合計(jì)算方法研究［J］.力學(xué)學(xué)報(bào)，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.