胡園園，謝 江，張 武

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

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

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

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

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

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

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

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

0 引言

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

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

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

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

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

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

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

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

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

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

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

(1)

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

(2)

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

(3)

k,l=0,1,2,…

(4)

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

2.1 三角形單元譜有限元法

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

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

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

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

P點的坐標(biāo)p(li,lj,lm)可由三個比值來唯一確定:

li=Δi/Δ

lj=Δj/Δ

lm=Δm/Δ

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

Δi+Δj+Δm=Δ

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

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

(6)

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

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

2.2 混合變分

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

同時，定義另外一個空間

記

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

求u∈V,p∈M，使得

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

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

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

根據(jù)變分問題

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

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

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

3 數(shù)值算例

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

3.1 串行數(shù)值實驗

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

Re=UupL/v

其中Uup為方腔上蓋速度,取方腔的邊長為特征速度L。

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

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

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

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

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

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

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

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

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

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

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

3.2 并行實驗設(shè)計

與有限元法計算剛度矩陣一樣，TSFEM先分別計算各個單元剛度矩陣，然后將單元剛度矩陣集成到整體剛度矩陣。不同單元剛度矩陣計算之間沒有直接聯(lián)系，非常適合針對不同單元剛度矩陣設(shè)計并行算法。采用OpenMP將上述算法并行化，對驅(qū)動方腔內(nèi)剪切流動問題進(jìn)行了計算。OpenMP是一個為在共享存儲的多處理機上編寫并行程序而設(shè)計的應(yīng)用編程接口，由一個小型的編譯器命令集組成，包括一套編譯制導(dǎo)語句和一個用來支持它的函數(shù)庫。

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

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

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

圖5 整體剛度矩陣并行計算流程

3.3 并行實驗性能分析

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

4 結(jié)語

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

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

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

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

4)運用并行計算可以在一定程度上縮短計算時間，提高計算效率。

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

References)

[1] 樊楓,徐國華,史勇杰.基于N-S方程的剪刀式尾槳前飛狀態(tài)氣動力計算研究[J].空氣動力學(xué)學(xué)報,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].數(shù)值計算與計算機應(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)定化有限元方法求解定常不可壓縮流動問題[J].工程數(shù)學(xué)學(xué)報,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)定化有限元方法[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].中國科學(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ū)動方腔流動問題的直接數(shù)值模擬[J].中國科學(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邊值問題[D].長春：吉林大學(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] 陳雪江,秦國良,徐忠.譜元法和高階時間分裂法求解方腔頂蓋驅(qū)動流[J].計算力學(xué)學(xué)報,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。 基金項目:國家自然科學(xué)基金重大研究計劃項目(91330116)。

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

1001-9081(2017)01-0042-06

10.11772/j.issn.1001-9081.2017.01.0042

TP301.6

A