求解非定常不可壓Navier-Stokes方程的一種高精度并行算法

倪詩浩，田振夫

(復旦大學 航空航天系，上海 200433)

?

求解非定常不可壓Navier-Stokes方程的一種高精度并行算法

倪詩浩，田振夫

(復旦大學 航空航天系，上海 200433)

摘要：采用一階投影法，建立了一種基于MPI求解非定常不可壓N-S方程的高精度并行算法.該算法在空間上可達到4階精度，其中，對流項中的1階導數(shù)和粘性項中的2階導數(shù)分別采用WENO格式和4階對稱型寬格式進行離散，而Poisson方程則采用4階精度的緊致格式進行迭代求解.通過對2維Taylor渦列和雙周期雙剪切邊界層流動問題及3維回轉(zhuǎn)體繞流問題的數(shù)值計算，驗證了算法的可靠性及其并行效率.

關鍵詞：并行算法； 標準消息傳遞界面； WENO格式； 投影法； 高精度

隨著計算機的高速發(fā)展，高性能計算日益成為人們用來研究各種數(shù)學物理問題的重要手段.隨著計算流體力學(Computational Fluid Dynamics， CFD)的不斷發(fā)展，為了更精細地研究復雜多尺度的湍流問題以及各類具有復雜外形的流動問題，計算量的需求越來越大，網(wǎng)格規(guī)模甚至達到了千萬級至上億級，單個CPU計算能力已難以滿足計算需要，同時也為了追求更加高效和強大的計算能力，開展并行計算也就成了必然.

本文基于分布存儲式MPI并行方法建立數(shù)值算法.MPI的全稱是Message Passing Interface，即標準消息傳遞界面，是一種基于消息傳遞的分布存儲式并行計算模式.分布存儲式并行計算機中，每一個CPU都有自己獨立的地址空間，但相互之間訪問不能直接進行，各個處理器之間必須通過顯式的消息傳遞來交換信息.對于MPI并行算法而言，通信是至關重要的，而復雜的消息傳遞語句的編寫是進行MPI并行程序設計的難點所在，但是其優(yōu)勢也是十分顯著的，它有很好的擴展性，因此被廣泛應用.

近年來，許多學者都針對粘性不可壓Navier-Stokes(N-S)方程設計了一些新算法.鞠立力和張林波等[1]針對重點構(gòu)造壓力Poisson方程快速解法設計了一種并行算法，并研究了其并行效率和可擴展性；郭曉虎等[2]以方腔自然對流為例，介紹了一種穩(wěn)定性、計算精度和并行性能都較佳的求解不可壓N-S方程的新方法；李寧和羅紀生[3]主要針對三對角矩陣的求解設計了一種MPI并行算法，并研究了平板邊界層問題，驗證了算法的有效性；王連生等[4]在Karniadakis[5]提出的不可壓縮算法的基礎上，提出了一種基于MPI的主從式并行算法，并對槽道湍流進行了數(shù)值模擬，取得了良好的效果.Gupta和Kalita[6]提出了一種基于流函數(shù)-速度形式下求解2維不可壓N-S方程的新方法，該方法解決了原先計算渦量時的一些難題，比如固壁邊界上的渦量處理，在計算方腔流動等問題時也取得了非常好的效果；Pandit等[7]介紹了一種在非規(guī)則網(wǎng)格下求解2維不可壓N-S方程的高階緊致方法，該方法能夠高效求解第一、第二類邊界條件下N-S方程的定常解，同時對自由流動和固壁邊界剪切流動也具有良好的捕捉能力；Zhang和Jackson[8]介紹了一種計算不可壓N-S方程的新方法，空間上引入了WENO(Weighted Essentially Non-Oscillatory)格式，時間方向采用了低耗散低色散Runge-Kutta方法，該方法在計算間斷初值問題和復雜強剪切流動問題時顯示出了良好的性能.

本文基于分布式并行運算的思想，采用有限差分算法，建立了一套基于原始變量投影法的黏性不可壓N-S方程的高精度數(shù)值方法.通過對2維Taylor渦列問題、雙周期雙剪切邊界層流動問題和3維回轉(zhuǎn)體繞流流動問題進行數(shù)值計算，對并行計算的效率進行研究，并考察算法的可靠性.

1控制方程的離散和數(shù)值算法

1.1控制方程

本文采用有限差分方法來求解原始變量下的粘性不可壓N-S方程.控制方程為

·V=0，

方程組的整體求解采用了投影法，投影法是根據(jù)Helmholtz-Hodge分解理論，引入一個中間速度場u*，對速度和壓力解耦后分步計算.

首先，第一步是預算步：

上式中，V表示速度矢量，上一步可計算得到中間速度場V*；

第二步計算壓力場：

結(jié)合連續(xù)性方程，不難得到

第三步就是利用求得的壓力校正獲得下一時間步的速度Vn+1：

Vn+1=V*-Δt·pn+1.

投影法是一種廣泛應用的計算方法，在不可壓縮流動的計算中具有十分重要的意義.近些年來不少學者也提出了基于投影法的高精度算法，如Tian[9]等提出過一種高精度緊致差分投影算法，王坪和田振夫[10]提出過一種2階投影方法.

1.2數(shù)值算法

在處理方程中的對流項時，我們采用了數(shù)值穩(wěn)定性較好的WENO格式[11]，WENO格式在計算間斷解和分辨強剪切問題中具有良好的數(shù)值穩(wěn)定性.因此，Zhang等[8]把WENO格式也引入到不可壓流動的計算中.為實現(xiàn)空間高精度，本文采用了數(shù)值穩(wěn)定性好、分辨率高的5階WENO格式(下文簡寫為WENO5)作為計算不可壓對流項的計算格式.

對于粘性項的處理，我們采用了4階寬格式：

第二步的Poisson方程求解，本文采用了文獻[12]提出的一種4階9點緊致格式，其離散格式如下：

其中：i=1，2，3，4分別為與中心點相鄰的4個點，i=5，6，7，8點分別是9點中的角點.

計算中，并行化的思想是將整體區(qū)域劃分為多塊，由各個處理器的內(nèi)存存儲，由其CPU進行計算，相互之間的通信方式則采用MPI的標準.

如圖1所示，整塊計算域被劃分成了12塊分區(qū)域，每塊區(qū)域設定為由一個處理器進行計算和存儲.但是，值得一提的是，由于我們采用的是分布存儲式的計算結(jié)構(gòu)，信息交換是這種并行化計算的關鍵部分，而通信是需要時間的，因此，并不是分區(qū)數(shù)量越多越好，需要平衡通信與計算的時間.如果分區(qū)過多，通信時間甚至超過了計算時間，這樣的并行化是沒有意義的.同時，也可以看到，并行化的效率和并行化的規(guī)模并不是成正比的，因為通信時間會隨并行規(guī)模的變大而增加.

在本文的格式中，我們需要傳遞的是分區(qū)邊界相鄰的數(shù)據(jù)，而這些數(shù)據(jù)在與該分區(qū)相鄰的其他分區(qū)里.而計算格式中，5階WENO格式需要兩側(cè)相鄰的3個點，2階導數(shù)的4階寬格式需要相鄰的2個點，Poisson方程的緊致格式迭代只需要相鄰的點.根據(jù)計算需求來看，我們把任一分區(qū)看作需要獲得傳值的中心分區(qū)，則需要給中心分區(qū)傳輸相鄰3層物理量的值，這樣，每一個計算域內(nèi)部分區(qū)的邊界就可以計算了.

2數(shù)值算例

本節(jié)中，將通過3個數(shù)值算例來考察我們建立的高精度并行算法的數(shù)值特性.

2.1Taylor渦列流動問題

Taylor渦列問題是一個具有解析解的經(jīng)典算例，可以用來檢驗算法的有效性.該問題以原始變量形式的粘性不可壓N-S方程作為控制方程，初始條件取為：

u(x，y，0)=-cos(Nx)sin(Ny)，v(x，y，0)=sin(Nx)cos(Ny)，

計算區(qū)域為0≤x≤2π，0≤y≤2π.該問題的邊界條件取為周期性邊界.其解析解為：

其中：N為整型參數(shù)；Re為流動雷諾數(shù).

在參數(shù)N=2，雷諾數(shù)分別為Re=100和Re=1000的條件下，時間步長Δt=0.002，計算時間T=2，Taylor渦列的計算結(jié)果如表1所示.計算中采用的WENO5的理論精度應該在3到5階之間，而在表1中實際計算結(jié)果顯示，當前算法的數(shù)值精度可以達到4階精度以上.

表1　WENO5格式計算得到的Taylor渦列速度u與解析解的誤差和收斂階

注： 表1中，41×41與81×81表示網(wǎng)格數(shù).

表2　81×81網(wǎng)格下不同計算核心數(shù)下的計算效率

2.2雙周期雙剪切邊界層問題

雙周期雙剪切邊界層問題是一個具有大梯度的非定常流動問題，可以測試數(shù)值算法在非定常問題中的數(shù)值特性.該問題中，在初始時刻，剪切邊界層中受到小擾動，經(jīng)過一段時間的演化，這種小擾動最終將導致形成大的渦結(jié)構(gòu)，問題的控制方程為非定常的粘性不可壓縮N-S方程.速度分量u，v的初始條件為：

v(x，y，0)=σsin(x).

圖2是對流項1階導數(shù)采用不同數(shù)值格式下的雙周期雙剪切邊界層問題在T=10時刻的渦量圖.從圖中可以看到，雖然4階中心格式和6階中心格式精度也非常高，但是其計算結(jié)果存在嚴重的數(shù)值振蕩.而WENO5格式和WENO7格式(即7階WENO格式)理論精度分別在3到5階和5到7階之間，精度與中心格式相差不大，但是WENO格式的分辨率卻要顯著高于中心格式，它能精細地分辨出高Re數(shù)下雙周期雙剪切邊界層的流場結(jié)構(gòu)，因此我們的算法中采用WENO格式是有效可信的.

2.33維回轉(zhuǎn)體繞流問題

本問題考慮的是一個水滴型回轉(zhuǎn)體模型在水下運動的流場計算，回轉(zhuǎn)體的線型采用經(jīng)典的SUBOFF模型，雷諾數(shù)Re=1.2×107，計算域取O型計算域，回轉(zhuǎn)體上游邊界和下游邊界分別離回轉(zhuǎn)體頭部駐點和尾部駐點為3.125倍回轉(zhuǎn)體長度，回轉(zhuǎn)體兩側(cè)沿流向邊界離回轉(zhuǎn)體側(cè)面約為3.817倍回轉(zhuǎn)體長度.無量綱化后，初始速度為1.

總體來說，計算結(jié)果是有效可信的.而并行化計算也克服了原先串行計算因單個CPU計算能力不足而無法計算的問題.在這個算例中，計算網(wǎng)格規(guī)模為282×200×240，這樣的規(guī)模在單個CPU下幾乎無法計算，而采用并行算法，本文計算時使用了1200個CPU，要獲得發(fā)展到穩(wěn)定的流場需要3個星期左右，而很多湍流問題的網(wǎng)格規(guī)模都在千萬以上，甚至上億級，并行計算往往需要幾千核甚至萬核級才能在可接受的時間內(nèi)完成計算.通過這個算例我們可以看到，并行化計算對于計算流體力學的進一步發(fā)展還是有很大的推動作用.

3結(jié)語

本文主要提出了一種基于MPI求解非定常不可壓N-S方程的高精度并行算法，設計并行化程序來求解原始變量形式下的黏性不可壓縮N-S方程.通過采用有限差分方法對2維Taylor渦列和雙周期雙剪切邊界層流動問題以及3維回轉(zhuǎn)體繞流問題進行了數(shù)值計算.通過驗證，可以看到，采用在MPI環(huán)境下分布式并行計算的思想，本文的算法可以達到約65%的并行效率，顯著提高了計算效率，同時，以高雷諾數(shù)下3維回轉(zhuǎn)體繞流問題為例，對于該大規(guī)模復雜問題的模擬，原先的串行模式難以實現(xiàn)，而在并行計算下得到了實現(xiàn).而本文設計的這種顯式并行計算的方法經(jīng)過計算算例的驗證也被證明是可信有效的.當然，如何減少通信時間和合理分區(qū)以獲得更高的并行效率依然是一個值得探討的問題.

參考文獻：

[1]鞠立力，張林波.解非定常不可壓N-S方程的壓力修正投影法的并行算法研究 [J].計算數(shù)學，1998，20(3)： 325-336.

[2]郭曉虎，田振夫，張林波.求解不可壓Navier-Stokes方程的一種新方法 [J].計算物理，2004，21(6)： 484-494.

[3]李寧，羅紀生.基于MPI的不可壓縮N-S方程并行計算方法的研究 [J].計算機工程與應用，2007，43(9)： 8-10.

[4]王連生，肖紅林，郭明明.基于MPI的不可壓縮N-S方程并行計算方法 [J].計算機應用與軟件，2011，28(8)： 279-283.

[5]KARNIADAKIS G E， ISRAELI M， ORSZAG S A. A high-order splitting method for the incompressible Navier-Stokes equations [J].JournalofComputationalPhysics， 1991，97(2)： 414-443.

[6]GUPTA M M， KALITA J C. A new paradigm for solving Navier-Stokes equations： Stream function-velocity formulation [J].JournalofComputationalPhysics， 2005，207(1)： 52-68.

[7]PANDIT S K， KALITA J C， DAHAL D C. A fourth-order accurate compact scheme for the solution of steady Navier-Stokes equations on non-uniform grids [J].Computers&Fluids， 2008，37(2)： 121-134.

[8]ZHANG J， JACKSON T L. A high-order incompressible flow solver with WENO [J].JournalofComputationalPhysics， 2009，228(7)： 2426-2442.

[9]TIAN Z F， LIANG X， YU P X. A high order compact finite difference algorithm for solving the incompressible Navier-Stokes equations [J].InternationalJournalforNumericalMethodsinEngineering， 2011，88(6)： 511-532.

[10]王坪，田振夫.基于投影法求解不可壓縮流的緊致格式 [J].工程數(shù)學學報，2010，27(2)： 225-232.

[11]JIANG G， SHU C W. Efficient implementation of weighted ENO schemes [J].JournalofComputationalPhysics， 1996，126(1)： 202-228.

[12]TIAN Z F， GE Y B. A fourth-order compact finite difference scheme for the steady stream function-vorticity formulation of the Navier-Stokes/Boussinesq quations [J].InternationalJournalforNumericalMethodsinFluids， 2003，41(5)： 495-518.

[13]李樹民，朱國林，王開春.基于分區(qū)的隱式求解二維不可壓N-S方程的并行實現(xiàn) [J].空氣動力學學報，2002，20(z1)： 39-44.

[14]JIMENEZ J M， HULTMARK M， SMITS A J. The intermediate wake of a body of revolution at high Reynolds numbers [J].JournalofFluidMechanics， 2010，659： 516-539.

[15]HUANG T， LIU H L， GROVES N，etal. Measurements of flows over an axisymmetric body with various appendages in a wind tunnel： The DARPA SUBOFF experimental program [C]. Proceedings of 19thSymposium on Naval Hydrodynamics， Seoul， Korea， 1992.

文章編號：0427-7104(2016)03-0315-06

收稿日期：2014-12-16

基金項目：國家自然科學基金(91330112，11372075)

作者簡介：倪詩浩(1991—)，男，碩士研究生；田振夫(1963—)，男，教授，通訊聯(lián)系人，E-mail： zftian@fudan.edu.cn.

中圖分類號：O 35

文獻標志碼：A

A High-order Parallel Numerical Method for Solving Unsteady Incompressible Navier-Stokes Equations

NI Shihao， TIAN Zhenfu

(DepartmentofAeronauticsandAstronautics，F(xiàn)udanUniversity，Shanghai200433，China)

Abstract：On the basis of the projection method， a high-order parallel algorithm for solving the unsteady incompressible Navier-Stokes(N-S) equations using MPI standard is presented. This algorithm can reach fourth-order accuracy in space. WENO(Weighted Essentially Non-Oscillatory) scheme is used to discretize the first derivative of the non-linear convective terms and the fourth-order symmetrical scheme is used in the discretization of the second derivatives of viscous terms respectively. Pressure Poisson equation is approximated using a fourth-order compact difference scheme. For assessing the effectiveness and parallel efficiency of the algorithm， a two-dimensional Taylor vortex column problem， a doubly-periodic double shear layer problem and a three-dimensional flow around a solid of revolution are numerically solved.

Keywords：parallel algorithm； MPI； WENO scheme； projection method； high-order