胡園園，謝 江，張 武

(1.上海大學(xué) 計(jì)算機(jī)工程與科學(xué)學(xué)院，上海200444; 2.上海大學(xué) 高性能計(jì)算中心，上海 200444)

(*通信作者電子郵箱jiangx@shu.edu.cn)

二維不可壓縮Navier-Stokes方程的并行譜有限元法求解

胡園園1，謝 江1*，張 武2

(1.上海大學(xué) 計(jì)算機(jī)工程與科學(xué)學(xué)院，上海200444; 2.上海大學(xué) 高性能計(jì)算中心，上海 200444)

(*通信作者電子郵箱jiangx@shu.edu.cn)

針對(duì)不可壓縮Navier-Stokes (N-S)方程求解過(guò)程中的有限元法存在計(jì)算網(wǎng)格量大、收斂速度慢的缺點(diǎn)，提出了基于面積坐標(biāo)的三角網(wǎng)格剖分譜有限元法(TSFEM)并進(jìn)一步給出了利用OpenMP對(duì)其并行化的方法。該算法結(jié)合譜方法和有限元法思想，選取具有無(wú)限光滑特性的指數(shù)函數(shù)取代傳統(tǒng)有限元法中的多項(xiàng)式函數(shù)作為基函數(shù)，能夠有效減少計(jì)算網(wǎng)格數(shù)量，提高算法的精度和收斂速度；利用面積坐標(biāo)便于三角形單元計(jì)算的特點(diǎn)，選取三角單元作為計(jì)算單元，增強(qiáng)了適用性；在頂蓋方腔驅(qū)動(dòng)流問(wèn)題上對(duì)該算法進(jìn)行驗(yàn)證。實(shí)驗(yàn)結(jié)果表明，TSFEM較傳統(tǒng)有限元法(FEM)無(wú)論是收斂速度還是計(jì)算效率都有了顯著提高。

不可壓縮N-S方程；OpenMP；方腔驅(qū)動(dòng)流；高精度；無(wú)窮收斂性

0 引言

Navier-Stokes (N-S)方程是流體力學(xué)中最重要的方程之一。數(shù)學(xué)家和物理學(xué)家深信，無(wú)論是微風(fēng)還是湍流，都可以通過(guò)求解N-S方程來(lái)進(jìn)行解釋和預(yù)言，因此研究N-S方程具有廣泛的應(yīng)用價(jià)值。對(duì)N-S方程的研究距今已有200多年的歷史，其弱解又稱(chēng)為L(zhǎng)eray-Hopf弱解。關(guān)于N-S方程強(qiáng)解的局部適定性、存在性與光滑性被列為21世紀(jì)7個(gè)價(jià)值100萬(wàn)美元的數(shù)學(xué)難題之一。數(shù)學(xué)家斷言，如果沒(méi)有新的分析工具和數(shù)學(xué)思想，這個(gè)難題將很難得到解決。但是,到目前為止,證明弱解的唯一性和正則性,即強(qiáng)解的整體存在性,仍是一個(gè)極具挑戰(zhàn)性的問(wèn)題。只有極少數(shù)非常簡(jiǎn)單的流動(dòng)問(wèn)題才能求得其精確解，大多數(shù)還是要用離散的方法求得數(shù)值解。在利用數(shù)值方法對(duì)其進(jìn)行求解時(shí)，計(jì)算格式的時(shí)間步長(zhǎng)受制于CFL(Courant-Friedrichs-Lewy)數(shù)，CFL=λ·Δτ/Δη，λ是與當(dāng)?shù)芈曀俪烧鹊奶卣骱瘮?shù)。當(dāng)流體近似不可壓縮時(shí)，λ趨向于無(wú)窮大，而CFL數(shù)是一個(gè)定值，因此最大時(shí)間步長(zhǎng)Δτ趨向于0。如果用可壓縮N-S方程求解近似不可壓縮流場(chǎng)，計(jì)算效率將會(huì)變得極低，因此只能通過(guò)不可壓縮N-S方程來(lái)計(jì)算近似不可壓流場(chǎng)。不可壓縮N-S方程的主要優(yōu)點(diǎn)是求解未知數(shù)的數(shù)量較少，將連續(xù)性方程和動(dòng)量方程聯(lián)立，求解得到速度和壓力后，代入能量方程就可以得到溫度，計(jì)算效率較高。

從原則上講，不可壓縮N-S方程作為描述粘性流體流動(dòng)的微分方程能夠求解流體運(yùn)動(dòng)的一切問(wèn)題。但是不可壓縮N-S方程是非線(xiàn)性偏微分方程，在實(shí)踐中，解決較為復(fù)雜的問(wèn)題，往往要考慮方程的所有項(xiàng)，因此求解非常困難。

現(xiàn)如今，求解不可壓縮N-S方程早已不單單是數(shù)學(xué)家關(guān)心的問(wèn)題，其在工程實(shí)際中的應(yīng)用日益廣泛而深入，例如近十年里，對(duì)于飛機(jī)抖振問(wèn)題已經(jīng)逐漸由原來(lái)粘、位流分開(kāi)計(jì)算的方法的研究轉(zhuǎn)向基于不可壓縮N-S方程的方法的研究[1]。故而在求解過(guò)程中，不但計(jì)算外形越來(lái)越復(fù)雜，結(jié)構(gòu)網(wǎng)格類(lèi)型已經(jīng)由單一的C型網(wǎng)格、O型網(wǎng)格、H型網(wǎng)格發(fā)展到塊結(jié)構(gòu)網(wǎng)格、嵌套網(wǎng)格等形式；網(wǎng)格劃分的精細(xì)程度也隨著人們對(duì)計(jì)算結(jié)果可靠性要求的增高而越來(lái)越細(xì)密。但是，單純地從細(xì)分網(wǎng)格出發(fā)去追求計(jì)算精度顯然是不切實(shí)際的，一味地細(xì)分網(wǎng)格只會(huì)因?yàn)橛?jì)算量過(guò)大導(dǎo)致最后計(jì)算無(wú)法進(jìn)行。因此，如何從計(jì)算方法本身出發(fā)，提高計(jì)算方法的效率和準(zhǔn)確性就成了越來(lái)越多的學(xué)者關(guān)心的問(wèn)題[2-4]。

有限元法是根據(jù)變分原理求解數(shù)學(xué)物理問(wèn)題的一種數(shù)值方法。自20世紀(jì)50年代提出以來(lái)，隨著矩陣?yán)碚?、?shù)值分析方法、特別是計(jì)算機(jī)技術(shù)的發(fā)展，有限元法無(wú)論是在基礎(chǔ)理論還是在實(shí)現(xiàn)技術(shù)上都取得了巨大的進(jìn)步[5]。它從最初的固體力學(xué)領(lǐng)域拓展到傳熱學(xué)、電磁學(xué)、流體力學(xué)以及聲學(xué)等其他物理場(chǎng)，從簡(jiǎn)單的靜力分析發(fā)展到動(dòng)態(tài)、非線(xiàn)性、多場(chǎng)耦合等復(fù)雜計(jì)算問(wèn)題[5]。由于有限元法可以對(duì)區(qū)域進(jìn)行任意的劃分，更適合用于復(fù)雜的計(jì)算域，故而在求解微分方程的實(shí)際計(jì)算中得到越來(lái)越多的應(yīng)用[6]。

傳統(tǒng)的有限元方法的通用性最好，但是有限元法在分析波的傳播時(shí)需要使單元大小與波的波長(zhǎng)相當(dāng)，且時(shí)間分辨率也非常小，使得計(jì)算效率較低。因此許多學(xué)者利用混合有限元法來(lái)求解微分方程。由于不可壓縮N-S方程具有速度和壓力兩個(gè)變量的性質(zhì)，利用混合有限元對(duì)其進(jìn)行求解時(shí)，離散方程經(jīng)常不滿(mǎn)足inf-sup穩(wěn)定化條件[7]，并且會(huì)產(chǎn)生較大的數(shù)值震蕩。在此基礎(chǔ)上，李劍等[8]和何銀年等[9]提出了基于局部Gauss積分的穩(wěn)定化方法。該方法不需要穩(wěn)定化參數(shù)，不需要考慮剖分網(wǎng)格邊界的條件，也不需要考慮不同網(wǎng)格單元力的變化；只需要計(jì)算兩個(gè)Gauss積分值之差。這個(gè)方法避開(kāi)了inf-sup條件，穩(wěn)定了一些低次有限元。許進(jìn)超[10]提出了用二層或多層網(wǎng)格的方法處理偏微分方程(PartialDifferentialEquation,PDE)問(wèn)題。該方法是在粗網(wǎng)格上處理一些離散的非線(xiàn)性問(wèn)題，得到一個(gè)初始值，然后在細(xì)網(wǎng)格上求一個(gè)離散線(xiàn)性方程組。Girault等[11]提出了一系列求解N-S方程的二層和多層方法。采用兩層穩(wěn)定有限元方法求解N-S方程,可以保持穩(wěn)定化方法在細(xì)網(wǎng)格上的精度，但是多層網(wǎng)格勢(shì)必會(huì)增加計(jì)算量，且在粗細(xì)網(wǎng)格的交界處處理起來(lái)十分困難，如若處理不當(dāng)，很容易造成二次人為誤差。除了有限元法以外，譜方法因其精度高的特點(diǎn)在求解不可壓縮N-S方程上也有較為廣泛的應(yīng)用。

譜方法是加權(quán)余量法的一種，源于經(jīng)典的Rize-Galerkin方法，它將未知量用正交函數(shù)代替，在整個(gè)計(jì)算區(qū)域內(nèi)把微分方程離散成代數(shù)方程組。這些正交函數(shù)往往是一些特殊二階常微分方程的特征函數(shù)，在數(shù)學(xué)上稱(chēng)之為譜。譜方法最大的魅力在于它具有“無(wú)窮階”收斂性，即如果原方程的解無(wú)窮光滑，那么用適當(dāng)?shù)淖V方法所求得的近似解以N-1的任意次速度收斂于精確解，這里的N為所選取的基函數(shù)的個(gè)數(shù)[12]。已有的譜方法一般適用于有界直角區(qū)域問(wèn)題和周期問(wèn)題，然而，在不可壓縮N-S方程的許多實(shí)際應(yīng)用領(lǐng)域中的問(wèn)題往往都是復(fù)雜邊界或是無(wú)界區(qū)域問(wèn)題和非周期問(wèn)題。因此，如何將有限元法與譜方法相結(jié)合，得到更為通用的高效求解方法是本文較為關(guān)注的問(wèn)題。

譜有限元法是譜方法與有限元方法的結(jié)合，它既有譜方法的高精度和收斂快的特點(diǎn)，也有有限元方法可以模擬任何復(fù)雜介質(zhì)模型的特點(diǎn)。譜有限元法中最廣為人知的是直接動(dòng)態(tài)剛度矩陣法、確切成員方法和譜元法[13]。譜有限元法和有限元法的主要區(qū)別是，有限元法是基于位移多項(xiàng)式的假設(shè)，譜有限元法是通過(guò)其插值函數(shù)來(lái)求解波動(dòng)方程的精確值。與傳統(tǒng)有限元法相比，譜有限元法的先進(jìn)性體現(xiàn)在可以通過(guò)單一的分析給出時(shí)域和頻域兩方面的結(jié)果[14]。

利用數(shù)值方法求解不可壓縮N-S方程，隨著網(wǎng)格劃分?jǐn)?shù)目越多，計(jì)算規(guī)模會(huì)越來(lái)越大，當(dāng)計(jì)算規(guī)模達(dá)到足夠大時(shí)，計(jì)算將會(huì)由于機(jī)器硬件資源的限制而無(wú)法進(jìn)行。近年來(lái)，隨著多核處理器的發(fā)展，并行計(jì)算在一定程度上得以普及，設(shè)計(jì)合適的并行算法，既可以提高計(jì)算速度，又可以緩解內(nèi)存的不足。譜有限元法和有限元法一樣，計(jì)算單元之間具有一定的獨(dú)立性，因此可以采用并行計(jì)算，將大規(guī)模重復(fù)性的單元分析分配在不同的節(jié)點(diǎn)上同時(shí)進(jìn)行，將會(huì)大大縮短計(jì)算的時(shí)間，提高計(jì)算的效率。

本文的主要工作是通過(guò)選取具有無(wú)限光滑特性指數(shù)函數(shù)作為基函數(shù)，構(gòu)造了基于面積坐標(biāo)的三角網(wǎng)格剖分譜有限元法(Triangular mesh Spectral Finite Element Method based on area coordinate， TSFEM)，并用其求解二維定常不可壓縮N-S方程，對(duì)于方程中的非線(xiàn)性項(xiàng)采用迎風(fēng)緊致格式計(jì)算，能夠有效地抑制混淆誤差[15]。本文模擬了頂蓋驅(qū)動(dòng)的二維方腔流動(dòng)，分析了實(shí)驗(yàn)所得到的結(jié)果，并與標(biāo)準(zhǔn)有限元法得到的結(jié)果相比較，得出并行TSFEM在求解二維定常不可壓縮N-S方程時(shí)具有精度高、適用性好和計(jì)算效率高的特點(diǎn)。

1 不可壓縮N-S方程的譜解法

據(jù)連續(xù)介質(zhì)力學(xué)知識(shí)，描述不可壓縮粘性流的控制微分方程，包括動(dòng)量平衡方程、連續(xù)方程、本構(gòu)方程以及能量方程。

(1)

方程(1)中，v是非負(fù)常數(shù)，稱(chēng)為粘性常數(shù)，且有

(2)

設(shè)所求的近似解為：

(3)

k,l=0,1,2,…

(4)

對(duì)于Fourier譜方法，要求：

(5)

其中：L是包含u和u關(guān)于空間變量倒數(shù)的算子，Lu=-v▽2u+(u·▽)u+▽p+f。將式(3)代入式(5)得到：

k∈[-N,N-1]

由式(4)得：

根據(jù)式(1)得：

2 TFSEM

譜有限元法的基本思想是：首先利用有限元的思想將求解區(qū)域劃分為若干單元，在每個(gè)單元內(nèi)將微分方程進(jìn)行變分，得到微分方程的積分弱形式；然后將特定參量和未知量在每個(gè)單元內(nèi)用高階正交多項(xiàng)式(如Chebyshev、Legendre和Fourier多項(xiàng)式等)進(jìn)行譜近似，得到每個(gè)單元的線(xiàn)性方程組，進(jìn)而得到單元?jiǎng)偠染仃?；利用有限元的技術(shù)進(jìn)行單元的疊加，將單元?jiǎng)偠染仃嚱M合成總體剛度矩陣；最后得到總求解域的線(xiàn)性代數(shù)方程組，加入邊界條件，從而求得整個(gè)區(qū)域的未知變量。有限元法只能通過(guò)增加單元來(lái)提高精度，譜有限元法不僅可以通過(guò)增加單元，還能通過(guò)提高每個(gè)單元的插值多項(xiàng)式階數(shù)來(lái)提高求解精度。此外，采用面積坐標(biāo)構(gòu)建三角形單元能夠有效地克服譜方法對(duì)復(fù)雜邊界問(wèn)題實(shí)用性不強(qiáng)的弱點(diǎn)。這兩方面的改進(jìn)使得譜有限元法相比其他數(shù)值方法具有更大的優(yōu)勢(shì)。

2.1 三角形單元譜有限元法

平面上有限元的計(jì)算單元通常有兩種：三角形單元和四邊形單元。一般來(lái)說(shuō)，有限元對(duì)區(qū)域的剖分要求很少，僅僅要求是相容的、正則的或者擬一致的這樣一些輕微條件。相對(duì)于三角形單元，四邊形單元分割起來(lái)相對(duì)簡(jiǎn)單，計(jì)算精度較高，比較適合較為規(guī)則的區(qū)域；但是三角形單元適應(yīng)性更好，求解具有復(fù)雜幾何構(gòu)形或具有復(fù)雜邊界的流動(dòng)問(wèn)題時(shí)一般不會(huì)出現(xiàn)壞單元，若將區(qū)域分割成四邊形單元很容易出現(xiàn)壞單元，導(dǎo)致計(jì)算不收斂。因此對(duì)于復(fù)雜區(qū)域問(wèn)題多采用三角形單元對(duì)其進(jìn)行劃分。

要在三角形單元上使用譜方法，必須建立起三角形區(qū)域上的正交多項(xiàng)式。

對(duì)于三角形單元，若用直角坐標(biāo)定義形函數(shù)，計(jì)算剛度矩陣將十分復(fù)雜；而改用面積坐標(biāo)后，公式可大為簡(jiǎn)化且積分運(yùn)算非常簡(jiǎn)單。故在此本文選擇用面積坐標(biāo)來(lái)表示三角形單元[17]。

設(shè)P為三角形單元中任一點(diǎn)，將P與三角形的三個(gè)頂點(diǎn)連接起來(lái)形成三個(gè)子三角形，如圖1所示。

P點(diǎn)的坐標(biāo)p(li,lj,lm)可由三個(gè)比值來(lái)唯一確定:

li=Δi/Δ

lj=Δj/Δ

lm=Δm/Δ

其中，Δ是△ijm的面積，Δi、Δj、Δm分別是△pjm、△pim、△pij的面積，且有

Δi+Δj+Δm=Δ

假設(shè)Ω是R2中的有限區(qū)域，f∈(L2(Ω))2，定長(zhǎng)的齊次Dirichlet邊值問(wèn)題如下。

求u=(u1,u2)T和p滿(mǎn)足：

(6)

其中：v是正常數(shù)，稱(chēng)為粘性常數(shù)；u稱(chēng)為流體的速度；p是流體的壓力；f是流體所受外力。

圖1 面積坐標(biāo)圖

2.2 混合變分

對(duì)方程(6)進(jìn)行混合變分，定義L2(Ω)的一個(gè)封閉子空間為：

同時(shí)，定義另外一個(gè)空間

記

由此可以得到方程(6)的弱形式：

求u∈V,p∈M，使得

用譜方法的思想進(jìn)行方程的離散，選取Chebyshev多項(xiàng)式的極值點(diǎn)作為插值點(diǎn):

2.3 選取插值基函數(shù)

選取Φk(x)=eikx,x∈[-1,1],k∈[0,N]作為插值基函數(shù)，則u(u1,u2)的插值表達(dá)式為：

根據(jù)變分問(wèn)題

a(ui,vi)+a1(w;ui,vi)+b(v,p)=(f,vi)

就可以分別得到各項(xiàng)的表達(dá)式，具體求解過(guò)程見(jiàn)文獻(xiàn)[18]。最后得到原微分方程的離散方程組：

譜有限元法與有限元法一樣要進(jìn)行問(wèn)題的劃分，將一個(gè)規(guī)模較大的求解問(wèn)題劃分成多個(gè)規(guī)模較小的單元求解問(wèn)題，并且單元內(nèi)的計(jì)算與其他單元之間具有一定的獨(dú)立性，因此譜有限元法同樣適用于做并行計(jì)算，這將大大縮短計(jì)算所需要的時(shí)間，進(jìn)一步提高計(jì)算的效率。

3 數(shù)值算例

方腔頂蓋驅(qū)動(dòng)流(簡(jiǎn)稱(chēng)方腔流)是一個(gè)經(jīng)典的非定常解問(wèn)題。幾十年來(lái)，很多數(shù)值方法的研究者都以驅(qū)動(dòng)方腔流動(dòng)作為模型[18-19]。一方面因其在工業(yè)生產(chǎn)中有著廣泛的應(yīng)用；另一方面，求解方腔頂蓋驅(qū)動(dòng)流問(wèn)題可以驗(yàn)證數(shù)值方法正確性，從而評(píng)價(jià)一種數(shù)值方法的優(yōu)劣。

3.1 串行數(shù)值實(shí)驗(yàn)

設(shè)計(jì)算區(qū)域?yàn)閇0,1]×[0,1]，選取Dirichlet邊界條件，即在邊界上函數(shù)值為零。流動(dòng)雷諾數(shù)為：

Re=UupL/v

其中Uup為方腔上蓋速度,取方腔的邊長(zhǎng)為特征速度L。

根據(jù)前人的數(shù)值模擬結(jié)果可以得到簡(jiǎn)單的結(jié)論，驅(qū)動(dòng)方腔流動(dòng)在Re<5 000的時(shí)候是定常解，在5 00010 000可以得到湍流解。

本文將利用在雷諾數(shù)等于100,400,1 000,3 200,5 000,7 200,10 000時(shí)候的標(biāo)準(zhǔn)解來(lái)評(píng)估本文算法的分辨率和數(shù)值精度(如表1所示)。

上邊界速度u=1,v=0，壓力取Dirichlet邊界條件▽p=0；其他三條邊界采用無(wú)滑移邊界條件u=0,v=0，壓力取Dirichlet邊界條件▽p=0。

計(jì)算結(jié)束判據(jù)為：

或者itrTimes<1 000，這里的itrTimes指的是迭代次數(shù)。

表1 不同雷諾數(shù)下實(shí)驗(yàn)的時(shí)間步長(zhǎng)與網(wǎng)格數(shù)

為了避免四個(gè)頂點(diǎn)的壓力值出現(xiàn)矛盾邊界值,在x=0,y=1點(diǎn)上令p=0,其他三個(gè)頂點(diǎn)上則使用線(xiàn)性外插得到的值作為邊界值。

本文計(jì)算了和Garoosi等[20]和Ghia等[21]相同雷諾數(shù)的驅(qū)動(dòng)方腔流動(dòng)。圖2和圖3分別給出了雷諾數(shù)為100，400，1 000，3 200，5 000，7 200和10 000情況下的流線(xiàn)圖和收斂圖。從圖2可以看出本文計(jì)算結(jié)果與參考文獻(xiàn)[20-21]給出的計(jì)算結(jié)果相吻合。

圖2 不同雷諾數(shù)情況下的流線(xiàn)

圖3 不同雷諾數(shù)情況下的收斂

由圖2中(a)～(b)可以看出，利用TFSEM求解方腔流問(wèn)題，在低密度網(wǎng)格的劃分的條件下，仍然可以獲得和高密度劃分下傳統(tǒng)的有限元方法較為一致的計(jì)算結(jié)果。這表明譜有限元法解的精度很高；由圖2(c)～(f)和圖3(c)～(f)可以看出在低密度網(wǎng)格劃分條件下，TFSEM能夠快速收斂。這表明TFSEM在計(jì)算二維不可壓縮N-S方程時(shí)具有譜方法的無(wú)窮收斂性。

3.2 并行實(shí)驗(yàn)設(shè)計(jì)

與有限元法計(jì)算剛度矩陣一樣，TSFEM先分別計(jì)算各個(gè)單元?jiǎng)偠染仃嚕缓髮卧獎(jiǎng)偠染仃嚰傻秸w剛度矩陣。不同單元?jiǎng)偠染仃囉?jì)算之間沒(méi)有直接聯(lián)系，非常適合針對(duì)不同單元?jiǎng)偠染仃囋O(shè)計(jì)并行算法。采用OpenMP將上述算法并行化，對(duì)驅(qū)動(dòng)方腔內(nèi)剪切流動(dòng)問(wèn)題進(jìn)行了計(jì)算。OpenMP是一個(gè)為在共享存儲(chǔ)的多處理機(jī)上編寫(xiě)并行程序而設(shè)計(jì)的應(yīng)用編程接口，由一個(gè)小型的編譯器命令集組成，包括一套編譯制導(dǎo)語(yǔ)句和一個(gè)用來(lái)支持它的函數(shù)庫(kù)。

OpenMP采用fork-join并行執(zhí)行模式，如圖4所示：OpenMP程序首先由主線(xiàn)程執(zhí)行。需要并行計(jì)算時(shí)，主線(xiàn)程派生出子線(xiàn)程來(lái)執(zhí)行并行任務(wù)。并行代碼執(zhí)行結(jié)束，子程序退出或者阻塞，不再工作，控制流程回到單獨(dú)的主線(xiàn)程中。定義在成對(duì)的fork和join之間的區(qū)域，稱(chēng)為并行域。

圖4 fork-join并行執(zhí)行模式

圖5給出了整體剛度矩陣并行計(jì)算流程，采用如圖5所示的并行計(jì)算流程。實(shí)驗(yàn)環(huán)境：CPU為4核；每個(gè)核分4個(gè)線(xiàn)程，共計(jì)16個(gè)線(xiàn)程，將計(jì)算任務(wù)平均分配給每個(gè)線(xiàn)程。

圖5 整體剛度矩陣并行計(jì)算流程

3.3 并行實(shí)驗(yàn)性能分析

圖6給出了不同了雷諾數(shù)下的加速比，通過(guò)觀(guān)察圖6可以清楚地看出，隨著雷諾數(shù)的增加，加速比在總體上呈現(xiàn)上升的趨勢(shì)，但是當(dāng)雷諾數(shù)Re>7 200時(shí)，加速比達(dá)到峰值并趨于穩(wěn)定。如何進(jìn)一步提高并行算法效率是下一步研究的問(wèn)題。

4 結(jié)語(yǔ)

譜有限元法即具有有限元法適應(yīng)性強(qiáng)的特點(diǎn)又具有譜方法的精度和無(wú)窮收斂的特性。本文構(gòu)造了二維不可壓縮N-S方程的并行譜有限元解法，其具體的優(yōu)點(diǎn)體現(xiàn)在如下幾點(diǎn)：

1)插值精度較高，對(duì)于具有波動(dòng)性的方程，只需要很少的點(diǎn)就可也模擬出原函數(shù)。

2)基函數(shù)是無(wú)限光滑的指數(shù)函數(shù)，具有無(wú)限可微性，適用于大雷諾數(shù)下湍流的計(jì)算。

3)基函數(shù)形式簡(jiǎn)單，與待插函數(shù)形式無(wú)關(guān)。

4)運(yùn)用并行計(jì)算可以在一定程度上縮短計(jì)算時(shí)間，提高計(jì)算效率。

本文使用譜基函數(shù)，推導(dǎo)出求解N-S方程的Galerkin積分表達(dá)式并通過(guò)并行的方式計(jì)算方腔流來(lái)進(jìn)行驗(yàn)算，得到了較好的結(jié)果，證實(shí)了此格式的適用性。但是當(dāng)實(shí)驗(yàn)網(wǎng)格數(shù)進(jìn)一步增多，進(jìn)程與進(jìn)程之間的通信量將會(huì)增大，計(jì)算時(shí)間反而會(huì)增加，使得計(jì)算效率低下。因此，如何減少并行計(jì)算之間的通信，使得計(jì)算能夠在更大規(guī)模的問(wèn)題上進(jìn)行，從而提高該算法在并行計(jì)算上的通用性將會(huì)成為接下來(lái)研究的重點(diǎn)。

References)

[1] 樊楓,徐國(guó)華,史勇杰.基于N-S方程的剪刀式尾槳前飛狀態(tài)氣動(dòng)力計(jì)算研究[J].空氣動(dòng)力學(xué)學(xué)報(bào),2014,32(4):527-533.(FAN F, XU G H, SHI Y J.Computation research on aerodynamic force of scissors tail-rotor in forward flight based on the N-S equations [J].Acta Aerodynamica Sinica, 2014, 32(4): 527-533.)

[2] REKATSINAS C S, NASTOS C V, THEODOSIOU T C, et al.A time-domain high-order spectral finite element for the simulation of symmetric and anti-symmetric guided waves in laminated composite strips [J].Wave Motion, 2015, 53:1-19.

[3] SAMARATUNGA D, JHA R, GOPALAKRISHNAN S.Wavelet spectral finite element for modeling guided wave propagation and damage detection in stiffened composite panels [J].Structural Health Monitoring, 2016, 15(3): 317-334.

[4] AXELSSON O, FAROUQ S, NEYTCHEVA M.A preconditioner for optimal control problems, constrained by Stokes equation with a time-harmonic control [J].Journal of Computational & Applied Mathematics, 2016, 310: 5-18.

[5] 王勖成.有限單元法[M].北京：清華大學(xué)出版社,2003:5-9.(WANG M C.Finite Element Method [M].Beijing: Tsinghua University Press, 2003: 5-9.)

[6] RANNACHER R.Finite Element Methods for the Incompressible Navier-Stokes Equations [M]// Fundamental Directions in Mathematical Fluid Mechanics.Berlin: Springer, 2000: 191-293.

[7] 楊建宏.定常Navier-Stokes方程的三種兩層穩(wěn)定有限元算法計(jì)算效率分析[J].數(shù)值計(jì)算與計(jì)算機(jī)應(yīng)用,2011,32(2):117-124.(YANG J H.Computation efficiency on three kinds of two-level stabilized finite element methods for stationary Navier-Stokes equations [J].Journal of Numerical Methods and Computer Applications, 2011, 32(2): 117-124.)

[8] 于佳平,史峰,李劍,等.應(yīng)用四邊形單元穩(wěn)定化有限元方法求解定常不可壓縮流動(dòng)問(wèn)題[J].工程數(shù)學(xué)學(xué)報(bào),2012,29(2):309-316.(YU J P, SHI F, LI J, et al.Numerical applications of the new stabilized quadrilateral finite element method for stationary incompressible flows [J].Chinese Journal of Engineering Mathematics, 2012, 29(2): 309-316.)

[9] 黃鵬展,何銀年,馮新龍.解Stokes特征值問(wèn)題的一種兩水平穩(wěn)定化有限元方法[J].應(yīng)用數(shù)學(xué)和力學(xué),2012,33(5):588-597.(HUANG P Z, HE Y N, FENG X L.A two-level stabilized finite element method for the Stokes eigenvalue problem [J].Applied Mathematics and Mechanics, 2012, 33(5): 588-597.)

[10] 許進(jìn)超.數(shù)值格式的穩(wěn)定性,相容性和收斂性[J].中國(guó)科學(xué)(數(shù)學(xué)),2015,45(8):1205-1216.(XU J C.On the stability, consistency and convergence of numerical schemes [J].Science in China(Series A), 2015, 45(8): 1205-1216.)

[11] CHACON REBOLLO T, GIRAULT V, MURAT F, et al.Analysis of a coupled fluid-structure model with applications to hemodynamics [J].SIAM Journal on Numerical Analysis, 2016, 54(2): 994-1019.

[12] BIRGERSSON F, FERGUSON N S, FINNVEDEN S.Application of the spectral finite element method to turbulent boundary layer induced vibration of plates [J].Journal of Sound & Vibration, 2003, 259(4): 873-891.

[13] GOPALAKRISHNAN S, CHAKRABORTY A, MAHAPATRA D R.Spectral Finite Element Method [M].London: Springer, 2008:177-217.

[14] 劉宏,傅德薰,馬延文.迎風(fēng)緊致格式與驅(qū)動(dòng)方腔流動(dòng)問(wèn)題的直接數(shù)值模擬[J].中國(guó)科學(xué)(A輯),1993,23(6):657-665.(LIU H, FU D X, MA Y W.Direct numerical simulation on the problem of upwind compact scheme and driven cavity flow [J].Science in China (Series A), 1993, 23(6): 657-665.)

[15] 向新民.譜方法的數(shù)值分析[M].北京：科學(xué)出版社,2000:48-90.(XIANG X M.Numerical Analysis of Spectral Method [M].Beijing: Science Press, 2000: 48-90.)

[16] 黃冬冬.用Fourier譜方法求解二維橢圓方程Dirichlet邊值問(wèn)題[D].長(zhǎng)春：吉林大學(xué),2010.(HUANG D D.The Fourier spectral method for solving two-dimensional elliptic partial differential equations with Dirichlet boundary condition [D].Changchun: Jilin University, 2010.)

[17] 王健平.有限譜有限元法求解二維不可壓縮Navier-Stokes方程[J].力學(xué)研究,2014,3(3):33-42.(WANG J P.A finite spectral finite element method for incompressible Navier-Stokes equations [J].International Journal of Mechanics Research, 2014, 3(3): 33-42.)

[18] GRILLET A M, YANG B, KHOMAMI B, et al.Modeling of viscoelastic lid driven cavity flow using finite element simulations [J].Journal of Non-Newtonian Fluid Mechanics, 1999, 88(1/2): 99-131.

[19] 陳雪江,秦國(guó)良,徐忠.譜元法和高階時(shí)間分裂法求解方腔頂蓋驅(qū)動(dòng)流[J].計(jì)算力學(xué)學(xué)報(bào),2002,19(3):281-285.(CHEN X J, QIN G L, XU Z.Spectral element method and higher-order fractional step method for lid driven flow in closed square cavity [J].Chinese Journal of Computational Mechanics, 2002, 19(3): 281-285.)

[20] GAROOSI F, BAGHERI G, RASHIDI M M.Two phase simulation of natural convection and mixed convection of the nanofluid in a square cavity [J].Powder Technology, 2015, 275: 239-256.

[21] OSSWALD G A, GHIA K N, GHIA U.Simulation of dynamic stall phenomenon using unsteady Navier-Stokes equations [J].Computer Physics Communications, 1991, 65(1/2/3): 209-218.

This work is supported by the Major Research Plan of the National Natural Science Foundation of China (91330116).

HU Yuanyuan, born in 1991.Her research interests include bioinformatics, high performance computing.

XIE Jiang, born in 1971, Ph.D., associate professor.Her research interests include bioinformatics, high performance computing.

ZHANG Wu, born in 1957, Ph.D., professor.His research interests include high performance computing, bioinformatics, computational fluid mechanics.

Solution of two dimensional incompressible Navier-Stokes equation by parallel spectral finite element method

HU Yuanyuan1, XIE Jiang1*, ZHANG Wu2

(1.SchoolofComputerEngineeringandScience,ShanghaiUniversity,Shanghai200444,China;2.HighPerformanceComputingCenter,ShanghaiUniversity,Shanghai200444,China)

Due to a large number of computational grids and slow convergence existed in the numerical simulation of Navier-Stokes (N-S) equation, Triangular mesh Spectral Finite Element Method based on area coordinate (TSFEM) was proposed.And further, TSFEM was paralleled with OpenMP.Spectral method was combined with finite element method, and the exponential function with infinite smoothness was selected as the basis function to replace the polynomial function in the traditional finite element method, which can efficiently reduce the amount of computational grids as well as improve the convergence and accuracy of the proposed algorithm.Because area coordinates can facilitate the calculation of triangular units, which were selected as the computing units to enhance the applicability of the algorithm.The lid-driven cavity flow was used to verify the TSFEM.The experimental results show that, compared with the traditional Finite Element Method (FEM), the TSFEM greatly improves the convergence rate and the calculation efficiency.

incompressible Navier-Stokes (N-S) equation; OpenMP; lid-driven cavity flow; high-precision; infinite convergence

2016-08-16;

2016-08-26。 基金項(xiàng)目:國(guó)家自然科學(xué)基金重大研究計(jì)劃項(xiàng)目(91330116)。

胡園園(1991—)，女，江蘇淮安人，CCF會(huì)員，主要研究方向：生物信息學(xué)、高性能計(jì)算； 謝江(1971—)，女，湖北恩施人，副教授，博士，CCF會(huì)員，主要研究方向：生物信息學(xué)、高性能計(jì)算； 張武(1957—)男，江西武寧人，教授，博士，CCF會(huì)員，主要研究方向：高性能計(jì)算、生物信息學(xué)、計(jì)算流體學(xué)。

1001-9081(2017)01-0042-06

10.11772/j.issn.1001-9081.2017.01.0042

TP301.6

A