三維頂板斜向驅(qū)動(dòng)方腔流有限元并行計(jì)算

王吉飛，萬(wàn)德成

(1. 上海交通大學(xué) 船舶海洋與建筑工程學(xué)院 海洋工程國(guó)家重點(diǎn)實(shí)驗(yàn)室，上海 200240； 2. 高新船舶與深海開(kāi)發(fā)裝備協(xié)同創(chuàng)新中心，上海 200240)

三維頂板斜向驅(qū)動(dòng)方腔流有限元并行計(jì)算

王吉飛1,2，萬(wàn)德成1,2

(1. 上海交通大學(xué) 船舶海洋與建筑工程學(xué)院 海洋工程國(guó)家重點(diǎn)實(shí)驗(yàn)室，上海 200240； 2. 高新船舶與深海開(kāi)發(fā)裝備協(xié)同創(chuàng)新中心，上海 200240)

基于PETScFEM開(kāi)源代碼，采用分步有限元算法和區(qū)域分解法，并行計(jì)算了不同雷諾數(shù)下的三維頂板斜向驅(qū)動(dòng)方腔流問(wèn)題。計(jì)算結(jié)果表明，當(dāng)頂板沿其對(duì)角線方向運(yùn)動(dòng)時(shí)，流體與下游側(cè)壁發(fā)生斜向碰撞后在下游對(duì)角處匯聚并形成射流，該射流在與底面、上游側(cè)壁碰撞后形成較為復(fù)雜的渦流結(jié)構(gòu)。雷諾數(shù)大小對(duì)三維頂板斜向驅(qū)動(dòng)方腔流的渦流場(chǎng)結(jié)構(gòu)形態(tài)具有重要影響。并行性能分析表明區(qū)域分解法能有效地提高三維粘性方腔流的計(jì)算速度。

三維方腔流；斜向角；有限元；區(qū)域分解；并行計(jì)算

頂板驅(qū)動(dòng)方腔流可模擬由某一壁面運(yùn)動(dòng)或外部流場(chǎng)運(yùn)動(dòng)所引起的方腔內(nèi)部環(huán)流。這類(lèi)流動(dòng)在航空、環(huán)境、工業(yè)、交通運(yùn)輸、海洋工程等領(lǐng)域廣泛存在，研究其流動(dòng)特性具有重要的工程應(yīng)用價(jià)值。Aidun 等[1]指出頂板驅(qū)動(dòng)方腔流在涂鍍和熔紡等工業(yè)裝置中起到重要的作用。在液壓管路運(yùn)輸中，適當(dāng)部位設(shè)置一個(gè)較大的方腔可減少或吸收液壓泵的壓力和流量脈動(dòng)對(duì)系統(tǒng)的影響，防止或減少液壓沖擊現(xiàn)象的發(fā)生。在航空領(lǐng)域，前緣縫翼、后緣縫翼、起落架箱和釘狀洞穴等所形成的方腔流動(dòng)是機(jī)身噪音的重要來(lái)源。在潛艇噪聲中，弦外孔和排水孔等所產(chǎn)生的內(nèi)部環(huán)流也是其水動(dòng)力噪聲的來(lái)源之一。研究方腔流對(duì)流動(dòng)機(jī)理的探討也有非常重要的意義，因?yàn)槠浜?jiǎn)單的設(shè)置能展示幾乎所有重要的流體力學(xué)現(xiàn)象，如角渦、展向渦、分叉、轉(zhuǎn)捩、湍流等。方腔流模型的簡(jiǎn)單性使得其結(jié)果的分析和比較更容易，同時(shí)極易擴(kuò)展到全雷諾數(shù)范圍。

隨著實(shí)驗(yàn)技術(shù)和計(jì)算機(jī)性能的不斷提高，人們對(duì)三維頂板驅(qū)動(dòng)方腔流進(jìn)行了越來(lái)越多的研究，其中對(duì)頂板沿平行于側(cè)壁方向運(yùn)動(dòng)引起的方腔流研究得最為廣泛。Prasad和Koseff[2]對(duì)頂板驅(qū)動(dòng)方腔流進(jìn)行了開(kāi)拓性的實(shí)驗(yàn)研究。Jiang等[3]，Wong等[4]以及其他學(xué)者對(duì)中、低雷諾數(shù)下的方腔流進(jìn)行了數(shù)值模擬，Leriche等[5]，Hachem等[6]，Zang等[7]對(duì)高雷諾數(shù)的方腔流進(jìn)行了數(shù)值模擬等。此外，Shankar等[8]對(duì)方腔流的研究進(jìn)展做了很好的總結(jié)。

但在現(xiàn)實(shí)情況中，頂板運(yùn)動(dòng)方向或外部流場(chǎng)方向并不一定平行于方腔側(cè)壁方向，而是與其有一定偏角。該偏角的存在使得方腔內(nèi)部渦流場(chǎng)結(jié)構(gòu)發(fā)生重大變化。研究頂板斜向驅(qū)動(dòng)方腔流具有工程實(shí)用價(jià)值，能更好地防止或減少管路輸運(yùn)中的液壓沖擊現(xiàn)象，更好地降低飛機(jī)機(jī)身或潛艇噪聲等，同時(shí)對(duì)流動(dòng)機(jī)理的探討也有重要的意義。國(guó)內(nèi)外學(xué)者對(duì)頂板斜向驅(qū)動(dòng)方腔流的研究較少，主要集中在飛機(jī)機(jī)身降噪方面。Disimile等、Czech等分別對(duì)飛機(jī)著陸裝置(可簡(jiǎn)化為頂板驅(qū)動(dòng)方腔流)的噪聲進(jìn)行了實(shí)驗(yàn)研究，結(jié)果表明當(dāng)外部流方向與方腔側(cè)壁成一定偏角時(shí)，不管流動(dòng)是否振蕩，方腔中均存在較強(qiáng)的聲學(xué)壓力。Pocitsky[9]基于商業(yè)軟件Fluent研究了頂板沿其對(duì)角線方向運(yùn)動(dòng)所引起的方腔流，其雷諾數(shù)包括100，400，700，2 000和40 000，結(jié)果顯示頂板斜向驅(qū)動(dòng)方腔流具有更復(fù)雜的渦流結(jié)構(gòu)。為了便于與前人的數(shù)值計(jì)算結(jié)果進(jìn)行比較，本文著重研究頂板沿其對(duì)角線方向運(yùn)動(dòng)所引起的方腔流。

根據(jù)控制方程的離散方式不同，計(jì)算流體力學(xué)數(shù)值離散方法大體上可分為三個(gè)分支[10]：有限差分法、有限體積法和有限元法。相比于其他方法，有限元方法具有如下特點(diǎn)和優(yōu)勢(shì)[11,12]：1)有限元法具有理論完整可靠，形式單純、規(guī)范，程序標(biāo)準(zhǔn)化，精度和收斂性得到保證等優(yōu)點(diǎn)；2)相比于有限差分法，有限元法對(duì)于求解區(qū)域的單元剖分沒(méi)有特別的限制，因此特別適合處理具有復(fù)雜邊界流場(chǎng)的區(qū)域；3)相比于有限體積法，有限元法可以通過(guò)提高單元插值多項(xiàng)式的次數(shù)來(lái)提高解的精度，且對(duì)單元正則性要求不高?；诖?，本文采用有限元法對(duì)不可壓粘性流場(chǎng)進(jìn)行數(shù)值模擬。

在對(duì)不可壓粘性流體進(jìn)行有限元分析時(shí)，速度和壓力插值空間函數(shù)需滿足inf-sup相容性條件，如不滿足，則可能產(chǎn)生虛假的壓力振蕩。為解決壓力穩(wěn)定性問(wèn)題，學(xué)者們相繼提出了不同的方法，如速度和壓力的不同階插值，罰函數(shù)法，附加壓力穩(wěn)定項(xiàng)等方法。由Chorin等和Temam等提出的基于POISSON投影的分步算法則可采用無(wú)需滿足相容性條件的速度壓力插值空間函數(shù)。Guermond等[13]進(jìn)一步驗(yàn)證，當(dāng)速度和壓力采用同階插值函數(shù)時(shí)，只有在時(shí)間步長(zhǎng)不太小時(shí)才能得到穩(wěn)定的壓力場(chǎng)。Codina等[14]得到了類(lèi)似的結(jié)論，并進(jìn)一步提出了具有壓力穩(wěn)定性的分步算法。由于有限元分步算法具有計(jì)算效率高，精度得到保證，編程較為簡(jiǎn)單等優(yōu)點(diǎn)[15]，文中采用該方法對(duì)三維方腔流問(wèn)題進(jìn)行數(shù)值模擬。

計(jì)算能力不足是計(jì)算流體力學(xué)發(fā)展所面臨的一大挑戰(zhàn)。為解決大規(guī)模計(jì)算問(wèn)題，并行化是大幅度提高計(jì)算效率的最有效手段。常用的并行化方法有線性代數(shù)方程組的并行計(jì)算和區(qū)域分解法兩種。并行計(jì)算工具箱PETSc(Portable, Extensible Toolkit for Scientific Computation)提供了并行計(jì)算線性代數(shù)方程組的庫(kù)函數(shù)。軟件包 PETScFEM 基于并行計(jì)算工具箱 PETSc，將非重疊區(qū)域分解法應(yīng)用于計(jì)算流體力學(xué)的有限元分析中，實(shí)現(xiàn)了不可壓粘性流體力學(xué)問(wèn)題的并行計(jì)算。

本文基于PETScFEM開(kāi)源代碼[16-17]，采用分步有限元算法對(duì)不可壓粘性流體進(jìn)行求解，并結(jié)合區(qū)域分解法并行計(jì)算了不同雷諾數(shù)下的頂板斜向驅(qū)動(dòng)方腔流問(wèn)題。數(shù)值模擬時(shí)首先對(duì)計(jì)算區(qū)域進(jìn)行有限元離散，不可壓 Navier-Stokes 方程中的速度和壓力采用同階的線性形函數(shù)和權(quán)函數(shù)，時(shí)間離散采用后退歐拉法；采用基于 Poisson 投影的分步算法對(duì)速度和壓力進(jìn)行解耦，并形成三個(gè)子方程；計(jì)算網(wǎng)格通過(guò)非重疊的區(qū)域分解法進(jìn)行剖分，并分配到不同的CPU中并行計(jì)算，MPI(Message Passing Interface)庫(kù)函數(shù)用于各子區(qū)域間并行計(jì)算數(shù)據(jù)通信，離散后的線性代數(shù)方程組采用 GMRES(Generalized Minimal RESidual)結(jié)合 Jacobi 預(yù)處理方法進(jìn)行計(jì)算。在之前的工作中[18]，已采用該方法對(duì)頂板正向驅(qū)動(dòng)方腔流問(wèn)題進(jìn)行了數(shù)值模擬，計(jì)算結(jié)果表明該方法計(jì)算效率高，精度較好。在此基礎(chǔ)上，文中的主要工作是對(duì)雷諾數(shù)為100、400、1 000和2 000的頂板斜向驅(qū)動(dòng)方腔流進(jìn)行并行數(shù)值模擬。

本文結(jié)構(gòu)安排如下：首先描述不可壓粘性流體控制方程及分步有限元算法，其次介紹非重疊區(qū)域分解并行算法，然后對(duì)不同雷諾數(shù)下的頂板斜向驅(qū)動(dòng)方腔流進(jìn)行數(shù)值模擬，給出速度矢量圖、壓力等值面、流線、渦量等值線以及中心垂線上的速度剖面等流場(chǎng)結(jié)構(gòu)信息，并進(jìn)行并行性能分析，最后給出結(jié)論。

1 數(shù)值算法

1.1 問(wèn)題描述及控制方程

圖1 三維頂板斜向驅(qū)動(dòng)方腔流計(jì)算模型Fig. 1 Computational model of 3D lid-driven cubic cavity flow at yaw

粘性流體流動(dòng)可用 Navier-Stokes 方程進(jìn)行求解。對(duì)于不可壓縮流體流動(dòng)，其控制方程包含動(dòng)量守恒方程和質(zhì)量守恒方程：

?tu+u·u+

1.2 時(shí)間和空間離散

為了方便描述有限元空間離散式，引入如下記號(hào)：

a(u,v):=ν(u,v),b(q,v):=(q,v),c(u,v,w):=(u·v,w)

(?tu,v)+c(u,u,v)+a(u,v)-b(p,v)=0 ?v∈V,

b(q,u)=0 ?q∈Q

采用梯形法則對(duì)方程進(jìn)行整體時(shí)間離散，即同時(shí)求解速度和壓力。簡(jiǎn)單起見(jiàn)，將時(shí)間 [0,T] 劃分為N個(gè)均勻的時(shí)間步，時(shí)間步長(zhǎng)為δt。令f為關(guān)于時(shí)間的通用函數(shù)，fn為函數(shù)f在tn=nδt時(shí)刻的值，fn+θ:=θf(wàn)n+1+(1-θ)fn，δtfn:=(fn+1-fn)/δt，其中θ∈[0,1]，則弱解方程的時(shí)間離散式為：

當(dāng)θ=1/2時(shí)，該方程對(duì)應(yīng)于二階的 Crank-Nicolson 時(shí)間離散，當(dāng)θ=1時(shí)，該方程對(duì)應(yīng)于向后 Euler 時(shí)間離散。

令Vh為V的近似有限元空間，Qh為Q的近似有限元空間，Vh空間中的函數(shù)要求為分段連續(xù)多項(xiàng)式，Qh空間中函數(shù)的連續(xù)性不要求，則方程(3)～(4)的有限元空間離散式為：

在基于壓力 Poisson 方程的分步算法中，當(dāng)時(shí)間步長(zhǎng)較大時(shí)，速度和壓力插值函數(shù)不需要滿足 LBB 條件，因此在本文的計(jì)算中速度和壓力插值函數(shù)均采用一階線性多項(xiàng)式。

方程(5)～(6)的矩陣形式可表示為：

式中：U和P分別表示結(jié)點(diǎn)速度和壓力數(shù)組，M為質(zhì)量矩陣，K為包含擴(kuò)散項(xiàng)和對(duì)流項(xiàng)的矩陣，G為梯度矩陣，D為散度矩陣。

1.3 分步算法

方程(7)～(8)等價(jià)于如下的方程組：

1.4 區(qū)域分解法

圖2 非重疊區(qū)域分解法計(jì)算模型(2個(gè)子區(qū)域)Fig. 2 Computational model of non-overlapping domain decomposition

令A(yù)ii,i=1,2,…,n為各子區(qū)域Ωi內(nèi)部結(jié)點(diǎn)所形成的矩陣，ALL=diag[A11,A22,…,Ann]為分塊對(duì)角矩陣，AII為各子區(qū)域間交界面結(jié)點(diǎn)所形成的矩陣。相應(yīng)地，令x=(xL,xI)T，y=(yL,yI)T，則代數(shù)方程組可分裂為：

式中：ALI和AIL為子區(qū)域內(nèi)部結(jié)點(diǎn)和交界面結(jié)點(diǎn)的連接關(guān)系矩陣。該方程組的求解可以分為兩步完成，首先求解方程(17)獲得交界面結(jié)點(diǎn)的變量值xI，各子區(qū)域間的數(shù)據(jù)通信通過(guò) MPI 庫(kù)函數(shù)來(lái)完成，然后再并行求解方程(18)，得到各子區(qū)域內(nèi)部結(jié)點(diǎn)的變量值xL：

2 數(shù)值模擬結(jié)果

首先模擬雷諾數(shù)為1 000時(shí)的頂板斜向驅(qū)動(dòng)方腔流，進(jìn)行了網(wǎng)格收斂性驗(yàn)證，并給出速度矢量圖、壓力等值面、流線、渦量等值線等流場(chǎng)信息。其次模擬雷諾數(shù)為100，400，2 000等不同頂板斜向驅(qū)動(dòng)方腔流，比較不同雷諾數(shù)對(duì)頂板斜向驅(qū)動(dòng)方腔流的影響。

2.1 雷諾數(shù)為1 000時(shí)的頂板斜向驅(qū)動(dòng)方腔流

首先模擬雷諾數(shù)為1 000時(shí)的頂板斜向驅(qū)動(dòng)方腔流。計(jì)算采用三套不同精度的均勻六面體網(wǎng)格進(jìn)行網(wǎng)格收斂性驗(yàn)證，粗、中、細(xì)網(wǎng)格的單元數(shù)分別為 36×36×36、48×48×48、64×64×64。計(jì)算時(shí)間步長(zhǎng)取為 0.05，迭代1 000步均能得到收斂的結(jié)果。計(jì)算結(jié)果表明，中等精度網(wǎng)格已能得到較好的結(jié)果，而計(jì)算時(shí)間卻比細(xì)網(wǎng)格少很多。圖3所示結(jié)果為中等精度網(wǎng)格計(jì)算結(jié)果。

雷諾數(shù)為1 000時(shí)頂板斜向驅(qū)動(dòng)方腔流的流線分布如圖3 所示。從正視圖可以看出，流體由頂板驅(qū)動(dòng)，以45°角撞擊x=1的下游側(cè)壁和y=1的下游側(cè)壁，撞擊后的流體在由兩側(cè)壁構(gòu)成的角落處匯聚，并形成一股強(qiáng)大的射流射向底面。從后視圖可以看出，射流撞向底面后在底面呈90°扇形散開(kāi)。流體在遇到x=0的上游側(cè)壁和y=0的上游側(cè)壁后轉(zhuǎn)而向上流動(dòng)，最終匯入頂板附近流體，完成一次循環(huán)。從側(cè)視圖可以看出，由于頂板驅(qū)動(dòng)流體的匯聚和射流的形成，使得主渦變得傾斜而狹長(zhǎng)，且有流體從下游側(cè)壁中心附近匯入主渦中。

圖3 雷諾數(shù)為1 000時(shí)流線分布的不同視角(左：正視圖，中：后視圖，右：側(cè)視圖)Fig. 3 Streamlines at Re=1 000: front view (left); back view (middle); oblique view (right)

從圖5可看到方腔流不同剖面的速度矢量投影。在CP剖面中，由于流體的堆積和射流的形成，可以看到主渦在平面的左下方，而右下角可看到次級(jí)渦的存在。在CP2剖面中，渦心位于剖面中央，平均速度較CP平面小。在CP3剖面中，渦心位于平面右上角，且大部分區(qū)域的速度矢量向上。從PP剖面和MP剖面的速度矢量投影中，可以看到二次渦的存在。本文結(jié)果與Povitsky的計(jì)算結(jié)果[9]基本一致。

圖6 所示為雷諾數(shù)為1 000時(shí)頂板斜向驅(qū)動(dòng)方腔流速度大小等值面。從圖6(a)中可以看到，由頂板帶動(dòng)的流體在由x=1的下游側(cè)壁和y=1的下游側(cè)壁形成的垂直角落里匯聚形成一股速度較大的射流。圖6(b)則顯示了由頂板帶動(dòng)的流體在撞擊下游側(cè)壁后匯聚的過(guò)程。圖6(c)顯示了射流在撞擊底面之后朝兩個(gè)方向散開(kāi)的低速射流。圖6(d)、6(e)和6(f)分別顯示了頂板驅(qū)動(dòng)流體在到達(dá)底面之后呈扇形面散開(kāi)，在遇到x=0的上游側(cè)壁和y=0的上游側(cè)壁后轉(zhuǎn)而向上流動(dòng)，由于流體運(yùn)動(dòng)空間的增大和壁面粘性的阻滯作用，這一過(guò)程中流體的速度逐漸降低。

圖4 不同剖面的速度矢量和流線投影及其法向渦量分量(Re=1 000)Fig. 4 Perspective 3D solution summary at Re=1 000

圖5 不同剖面的速度矢量投影分布(Re=1 000)Fig. 5 2D planar projections of velocity vector at Re=1 000 on different planes

圖7 所示為雷諾數(shù)為1 000時(shí)頂板斜向驅(qū)動(dòng)方腔流的壓力等值面。圖8所示為雷諾數(shù)為1 000時(shí)不同精度網(wǎng)格計(jì)算結(jié)果(沿 CP 剖面中垂線的水平速度分量)比較。從圖8中可以看到，從粗網(wǎng)格到細(xì)網(wǎng)格，計(jì)算結(jié)果逐漸收斂，中等精度網(wǎng)格與細(xì)網(wǎng)格計(jì)算結(jié)果基本保持一致。由于中等精度網(wǎng)格已能得到較好的計(jì)算精度，而網(wǎng)格單元數(shù)卻比細(xì)網(wǎng)格少很多，可以減少大量的計(jì)算時(shí)間，因此在本文其它的工況計(jì)算中均采用中等精度網(wǎng)格進(jìn)行計(jì)算。從速度剖面看，水平速度分量在底面附近存在一個(gè)峰值，該峰值由底面射流產(chǎn)生。在距底面0.1～0.3的區(qū)域速度絕對(duì)值均較大，由流體的匯聚擠壓引起。

圖6 頂板斜向驅(qū)動(dòng)方腔流速度大小等值面(Re=1 000)Fig. 6 Iso-surfaces for different velocity magnitudes at Re=1 000

圖7 壓力等值面(Re=1 000)Fig. 7 Pressure iso-surfaces at Re=1 000

圖8 沿CP 剖面中垂線的水平速度分量(Re=1 000)Fig. 8 Horizontal component velocity profile at Re=1 000

2.2 不同雷諾數(shù)對(duì)頂板斜向驅(qū)動(dòng)方腔流的影響

本文還模擬了雷諾數(shù)為100、400和2 000時(shí)的頂板斜向驅(qū)動(dòng)方腔流，并考察不同雷諾數(shù)對(duì)頂板斜向驅(qū)動(dòng)方腔流的影響。本節(jié)計(jì)算采用48×48×48的均勻六面體網(wǎng)格，單元特征尺度為0.020 83。計(jì)算時(shí)間步長(zhǎng)取為0.05，迭代1 000步均能得到收斂的結(jié)果。

圖9所示為雷諾數(shù)為100、400和2 000時(shí)頂板斜向驅(qū)動(dòng)方腔流的流線分布。從圖中可以看出，隨著雷諾數(shù)的增加，流線由簡(jiǎn)單變?yōu)閺?fù)雜，特別是主渦流線，當(dāng)雷諾數(shù)為2 000時(shí)變得雜亂。

圖10所示為不同雷諾數(shù)下CP剖面內(nèi)速度矢量投影。從圖中可以看出，隨著雷諾數(shù)的增加，主渦渦心位置不斷變化，當(dāng)雷諾數(shù)為100時(shí)，主渦位于中部偏上的位置，當(dāng)雷諾數(shù)為400時(shí)，主渦中心在右上角，當(dāng)雷諾數(shù)為1 000時(shí)，主渦變?yōu)閮A斜狹長(zhǎng)的卵形，渦心在平面的左下方，當(dāng)雷諾數(shù)達(dá)到2 000時(shí)，渦心位置更靠近底面。此外，當(dāng)雷諾數(shù)為400時(shí)，可以看到較為明顯的下游二次渦。當(dāng)雷諾數(shù)達(dá)到2 000時(shí)，可以明顯地觀察到上游二次渦的形成。

圖9 不同雷諾數(shù)(從上至下：Re=100，400，2 000)流線分布的不同視角(從左至右：正視圖，后視圖，側(cè)視圖)Fig. 9 Streamlines for different Reynolds numbers

圖10 不同雷諾數(shù)下CP剖面內(nèi)速度矢量投影Fig. 10 2D planar projections of velocity vector on CP plane for different Reynolds numbers

圖11所示為不同雷諾數(shù)下 MP 剖面內(nèi)速度矢量投影(本文結(jié)果)。從圖中可以看出，隨著雷諾數(shù)的增加，二次流變得越來(lái)越復(fù)雜，當(dāng)雷諾數(shù)為1 000時(shí)，可以觀察到3對(duì)二次渦，而當(dāng)雷諾數(shù)增加到2 000時(shí)，可以觀察到6對(duì)二次渦。圖12所示為Povitsky的計(jì)算結(jié)果[9]。將本文結(jié)果與參考結(jié)果進(jìn)行比較，可以看出對(duì)應(yīng)的漩渦個(gè)數(shù)及漩渦形態(tài)均保持一致。

圖11 不同雷諾數(shù)下MP剖面內(nèi)速度矢量投影(本文結(jié)果)Fig. 11 2D planar projections of velocity vector on MP plane for different Reynolds numbers (our results)

圖12 不同雷諾數(shù)下MP剖面內(nèi)速度矢量投影(參考結(jié)果)Fig. 12 2D planar projections of velocity vector on MP plane for different Reynolds numbers (reference results)

圖13所示為不同雷諾數(shù)下沿CP剖面中垂線的水平速度分量。從圖中可以看到，隨著雷諾數(shù)的增加，邊界層的厚度逐漸減小。當(dāng)雷諾數(shù)為100時(shí)，粘性作用較強(qiáng)，壁面的阻滯效應(yīng)較大，速度剖面變化較為平緩。當(dāng)雷諾數(shù)為400時(shí)，底面附近已形成射流，因此出現(xiàn)了較大的峰值。當(dāng)雷諾數(shù)為1 000時(shí)，底面附近的射流峰值有所減小，但流體的堆積擠壓現(xiàn)象較為明顯，導(dǎo)致距底面0.1～0.3附近的速度絕對(duì)值較大。當(dāng)雷諾數(shù)達(dá)到2 000時(shí)，流體堆積擠壓現(xiàn)象更為嚴(yán)重，使得流體在距底面0.2附近出現(xiàn)峰值。

圖13 不同雷諾數(shù)下沿CP剖面中垂線的水平速度分量Fig. 13 Horizontal velocity component along the vertical midline of the CP plane at different Reynolds numbers

不同雷諾數(shù)下沿 CP 剖面中垂線的水平速度分量數(shù)值如表1所示。

2.3 并行性能

本文采用網(wǎng)格非重疊區(qū)域分解法結(jié)合消息傳遞的模式進(jìn)行并行計(jì)算，非重疊網(wǎng)格區(qū)域劃分的方式如圖14所示。以雷諾數(shù)為1 000的頂板斜向驅(qū)動(dòng)方腔流計(jì)算為例，取時(shí)間步長(zhǎng)為0.05，迭代1 000步，不同進(jìn)程數(shù)(從1到8)的計(jì)算時(shí)間及加速比如表2所示，相對(duì)應(yīng)的曲線如圖15所示。從圖中可以看出，隨著進(jìn)程數(shù)的增加，計(jì)算時(shí)間不斷減少，并行加速比不斷增大，但由于進(jìn)程間數(shù)據(jù)傳遞量的增加，加速比的增長(zhǎng)速度比線性增長(zhǎng)速度為小。

表1 不同雷諾數(shù)下沿 CP 剖面中垂線的水平速度分量Tab. 1 Horizontal velocity component along the vertical midline of the CP plane at different Reynolds numbers

圖14 不同塊數(shù)的區(qū)域分解網(wǎng)格Fig. 14 Domain decomposition meshes for different domains

進(jìn)程數(shù)12345678時(shí)間/s4323020520156601286010540951791288171加速比12.10672.76053.36164.10154.54244.73605.2907

圖15 并行性能分析Fig. 15 Analysis of the parallel performance

3 結(jié) 語(yǔ)

采用分步有限元算法結(jié)合非重疊區(qū)域分解法并行計(jì)算了不同雷諾數(shù)(100、400、1 000 和2 000)下的頂板斜向(斜向角為 45°)驅(qū)動(dòng)方腔流問(wèn)題。從計(jì)算結(jié)果可以看出，頂板斜向驅(qū)動(dòng)方腔流具有較為復(fù)雜的渦流結(jié)構(gòu)。由頂板驅(qū)動(dòng)的流體以45°角撞擊x=1的下游側(cè)壁和y=1的下游側(cè)壁后匯聚成一股速度較大的射流，該射流在底面呈扇形散射，在遇到x=0的上游側(cè)壁和y=0的上游側(cè)壁后轉(zhuǎn)向上流動(dòng)，最終匯入頂板附近流體并進(jìn)入循環(huán)。在這一循環(huán)過(guò)程中，由于主流運(yùn)動(dòng)方向與側(cè)壁呈45°偏角，使得二次流現(xiàn)象較為明顯和復(fù)雜。且隨著雷諾數(shù)的增加，流動(dòng)更為劇烈，整個(gè)渦流場(chǎng)結(jié)構(gòu)也更為復(fù)雜。

計(jì)算表明，采用基于Poisson投影的分步算法后，應(yīng)用傳統(tǒng)的有限元方法即能很好地求解不可壓粘性方腔流問(wèn)題。同時(shí)，結(jié)合網(wǎng)格區(qū)域分解和消息傳遞的并行計(jì)算模式，可以得到很好的并行計(jì)算性能。開(kāi)源軟件PETScFEM為計(jì)算不可壓粘性流體問(wèn)題提供了很好的工具，且為進(jìn)一步開(kāi)發(fā)提供了很好的平臺(tái)。

[1] AIDUN C K, TRIANTAFILLOPOULOS N G, BENSON J D. Global stability of a lid-driven cavity with through flow: flow visualization studies [J]. Physics of Fluids, 1991, 3: 2081-2091.

[2] PRASSAD A K, KOSEFF J R. Reynolds number and end-wall effects on a lid-driven cavity flow [J]. Physics of Fluids, 1989,1(2): 208.

[3] JIANG B N, LIN T L, POVINELLI L A. Large scale computation of incompressible viscous flows by least-squares finite element method [J]. Computer Methods in Applied Mechanics and Engineering, 1994, 114(3-4): 213-231.

[4] WONG K L, BAKER A J. A 3D incompressible Navier-Stokes velocity-vorticity weak form finite element algorithm [J]. International Journal for Numerical Methods in Fluids, 2002, 38(9): 9-23.

[5] LERICHE E. Direct numerical simulation in a lid-driven cubical cavity at high reynolds number by a chebyshev spectral method [J]. Journal of Scientific Computing, 2006: 335-345.

[6] HACHEM E, RIVAUX B, KLOCZKO T, et al. Stabilized finite element method for incompressible flows with high Reynolds number [J]. Journal of Computational Physics, 2010, 229: 643-665.

[7] ZANG Y, STREET RL, KOSEFF JR. A dynamic mixed subgrid-scale model and its application to turbulent recirculating flows [J]. Physics of Fluids, 1993, 5: 3186.

[8] SHANKAR PN, DESHPANDE MD. Fluid mechanics in the driven cavity [J]. Annual Review of Fluid Mechanics, 2000, 32: 93-136.

[9] POVITSKY A. Three-dimensional flow in cavity at yaw [J]. Nonlinear Analysis, 2005, 63(5-7): 1573-1584.

[10] 任安祿. 不可壓縮粘性流場(chǎng)計(jì)算方法 [M]. 北京: 國(guó)防工業(yè)出版社, 2003.(REN Anlu. Numerical methods in incompressible viscous flow [M]. Beijing: National Defence Industry Press, 2003.(in Chinese))

[11] 章本照，印建安，張宏基. 流體力學(xué)數(shù)值方法 [M]. 北京: 機(jī)械工業(yè)出版社, 2003. (ZHANG Benzhao, YIN Jianan, ZHANG Hongji. Numerical methods in fluid dynamics [M]. Beijing: China Machine Press, 2003.(in Chinese))

[12] 王獻(xiàn)孚，周樹(shù)信，陳澤梁，等. 計(jì)算船舶流體力學(xué) [M]. 上海: 上海交通大學(xué)出版社, 1991. (WANG Xianfu, ZHOU Shuxin, CHEN Zeliang, et al. Computational ship hydrodynamics [M]. Shanghai: Shanghai Jiao Tong University Press, 1991.(in Chinese))

[13] GUERMOND J L, QUARTAPELLE L. On stability and convergence of projection methods based on pressure poisson equation [J]. International Journal for Numerical Methods in Fluids, 1998, 26: 1039-1053.

[14] CONDINA R. Pressure stability in fractional step finite element methods for incompressible flows [J]. Journal of Computational Physics, 2001, 170(1): 12-40.

[15] 劉淼兒，任玉新，張涵信. 數(shù)值求解不可壓縮流動(dòng)的投影方法研究進(jìn)展 [J]. 力學(xué)進(jìn)展, 2006, 36(4): 591-598. (LIU Miaoer, REN Yuxin, ZHANG Hanxin. Review on the projection methods in the numerical solution of the incompressible flow [J]. Advances in Mechanics, 2006, 36(4): 591-598.(in Chinese))

[16] PAZ R R, NIGRO N M, STORTI M A. On the efficiency and quality of numerical solutions in CFD problems using the interface strip preconditioner for domain decomposition methods [J]. International Journal for Numerical Methods in Fluids, 2006, 52(1): 89-118.

[17] SONZOGNI V E, YOMMI A M, NIGRO N M, et al. A parallel finite element program on a Beowulf cluster [J]. Advances In Engineering Software, 2002, 33(7-10): 427-443.

[18] WANG J F, WAN D C. Parallel simulation of 3D lid-driven cubic cavity flows by finite element method [C]// Proceedings of the 21st International Offshore and Polar Engineering Conference. Hawaii:[s.n.], 2011: 644-651.

Parallel simulation of 3D lid-driven cubic cavity flows at yaw by finite element method

WANG Jifei1,2, WAN Decheng1,2

(1. State Key Laboratory of Ocean Engineering, School of Naval Architecture, Ocean and Civil Engineering, Shanghai Jiao Tong University, Shanghai 200240, China; 2. Collaborative Innovation Center for Advanced Ship and Deep-Sea Exploration, Shanghai 200240, China)

The fractional step finite element method and domain decomposition method are applied to simulate in parallel the 3D lid-driven cubic cavity flows at yaw based on the open source codes PETScFEM. When the lid moves along its diagonal, the driven fluid impinges in the spatial angle formed by the downstream side walls, which leads to formation of the jet flow. A system of vortices is caused when the jet flow impinges the bottom wall and upstream side walls. Different Reynolds numbers are investigated, which show the significant influence for this flow. Parallel performance analysis reveals that the domain decomposition method can efficiently speed up the simulation of 3D lid-driven cavity flows.

3D cavity flow, yaw angle, finite element method, domain decomposition, parallel simulation

O351

A

10.16483/j.issn.1005-9865.2015.02.001

1005-9865(2015)02-0001-12

2014-01-06

國(guó)家自然科學(xué)基金資助項(xiàng)目(51379125，11432009,51490675)；上海高校特聘教授(東方學(xué)者)崗位跟蹤計(jì)劃項(xiàng)目(2013022)；國(guó)家重點(diǎn)基礎(chǔ)研究發(fā)展計(jì)劃(973計(jì)劃)項(xiàng)目(2013CB036103)；工信部高技術(shù)船舶科研項(xiàng)目；上海交通大學(xué)高性能計(jì)算中心(HPC)資助項(xiàng)目

王吉飛(1984-)，男，重慶人，博士研究生，從事計(jì)算流體力學(xué)研究。E-mail：wangjifei2000@126.com

萬(wàn)德成。E-mail：dcwan@sjtu.edu.cn