黃筱云 ，李紹武

(1. 天津大學水利工程仿真與安全國家重點實驗室，天津 300072；2. 長沙理工大學水利工程學院，長沙 410114；3. 長沙理工大學水沙科學與水災害防治湖南省重點實驗室，長沙 410114)

采用歐拉方法模擬不可壓縮黏性流動時，計算精度在很大程度上取決于網(wǎng)格的尺寸．然而，對網(wǎng)格進行全局細化會增加計算機內存用量和計算時間．與同分辨率的均勻網(wǎng)格數(shù)值模型相比，自適應網(wǎng)格數(shù)值模型能夠在保證計算結果精度的前提下有效減少計算網(wǎng)格數(shù)量．

與有限單元法和有限體積法采用的自適應三角網(wǎng)格不同，有限差分法采用的是一種按尺寸大小分級的正方形網(wǎng)格．網(wǎng)格等級結構有2 種：一種是自適應網(wǎng)格細分(adaptive mesh refinement，AMR)結構[1-3]；另一種是樹型網(wǎng)格結構[4]，二維情況稱為四叉樹結構．AMR 結構是將計算區(qū)域劃分為若干大小不一的矩形子區(qū)域，每個子區(qū)域再劃分為大小均一的正方形網(wǎng)格．有限差分計算從低到高逐層進行，每層計算結果作為下一層計算的邊界條件；反過來，自高到低，每層計算結果又用于對下一層計算結果的修正．四叉樹結構則對每層計算區(qū)域形狀不做任何要求，而通過采用特定的有限差分格式直接在葉網(wǎng)格上進行計算．

2003 年，Popinet[5]首先在平衡的樹型網(wǎng)格結構上對不可壓縮理想流體進行了模擬，其中，壓力和速度都定義在網(wǎng)格中心，Euler 方程采用兩步投影法求解，壓力泊松方程則通過多重網(wǎng)格法來計算．此后，Losasso 等[6]提出一種交錯的平衡四叉樹網(wǎng)格結構，即將壓力定義在網(wǎng)格中心，速度分量定義在網(wǎng)格邊上，對不可壓縮理想流體進行了模擬，其中，對流項通過半拉格朗日方法求解，泊松方程則采用近似中心差分進行離散，以達到簡化計算的目的．2006 年，Losasso 等[7]又進一步提出了基于支撐算子方法[8](support operators method，SOM)的離散格式，將泊松方程在樹型網(wǎng)格下解的精度從一階提高到二階，且得到的系數(shù)矩陣保持對稱．到目前為止，應用自適應四叉樹網(wǎng)格結構建立的多數(shù)流體運動模型主要是基于理想流體運動方程，忽略了黏性項，因此不具普適性．Min 等[9]曾將壓力和速度都定義在網(wǎng)格節(jié)點上，提出了一種四叉樹下偏導數(shù)的離散新格式以求解黏性項和壓力泊松方程，但泊松方程離散后的系數(shù)矩陣為非對稱，方程求解難度顯著增加．

本文在Losasso 的基礎上提出一種平衡四叉樹網(wǎng)格下黏性項的新離散方法，并建立了一種新的不可壓縮黏性流體運動數(shù)值模型．該模型可通過定義在網(wǎng)格節(jié)點上渦度值的大小對網(wǎng)格尺度進行自動調整，既保證復雜形狀流場區(qū)域具有足夠網(wǎng)格分辨率，又不顯著增加網(wǎng)格數(shù)量．此外，新模型采用指針和鏈表來描述網(wǎng)格層級關系，采用具有二階精度的無條件穩(wěn)定MacCormack 方法[10]來計算對流項以保證計算結果具有較高精度．最后，通過單渦、方腔流和圓柱繞流等典型算例來驗證新模型的有效性．

1 控制方程

新模型的控制方程為不可壓縮黏性流體的N-S 方程，即

式中：ui為速度張量；ρ為流體密度；p為壓力；τij為黏性應力張量；gi為體積力張量．對于牛頓流體，黏性應力張量可表示為

式中μ為動力黏性系數(shù)．

不可壓縮黏性流體的N-S 方程將通過兩步投影法進行求解[11]．第1 步，計算中間速度，其計算式為

式中Δt 為時間步長．第2 步，將中間速度場投影到散度為0 的平面上以獲得最終速度，即

在式(5)的兩端取散度，并運用連續(xù)性條件式(6)，可得

該方程稱為壓力泊松方程．

2 四叉樹網(wǎng)格結構與變量定義空間布置

等級結構網(wǎng)格如圖1 所示，其中四叉樹網(wǎng)格屬于非結構網(wǎng)格，網(wǎng)格形狀為正方形．每個正方形網(wǎng)格對應四叉樹的一個結點，可以有上一層父網(wǎng)格和下一層子網(wǎng)格(如圖2 所示，NE、NW、SW 和SE 分別表示根網(wǎng)格的4 個象限)，尺寸最大的網(wǎng)格所在層為0層，對應的節(jié)點為根結點，沒有子網(wǎng)格的為末端的葉網(wǎng)格．每個0 層網(wǎng)格所包含的全部網(wǎng)格的集合稱為一個四叉樹結構，若干個四叉樹形成的網(wǎng)格結構群稱為四叉樹森林．每個四叉樹網(wǎng)格有4 個角點，相鄰角點的連線稱為網(wǎng)格的側邊，而側邊上相鄰網(wǎng)格節(jié)點的連線稱為網(wǎng)格的邊(圖3)．只包含一條邊的側邊稱為最小公共邊．包含相同公共邊的2 個網(wǎng)格互為鄰居．

為敘述方便，引入指針來說明建立四叉樹網(wǎng)格結構的基本思路．每個網(wǎng)格有4 個指針指向其4 個子網(wǎng)格，1 個指針指向其父網(wǎng)格．同時，每個網(wǎng)格還有分別指向其4 個角點和4 個側邊的指針．為了方便尋找鄰居，每個網(wǎng)格節(jié)點有4 個分別指向其4 個象限處葉網(wǎng)格的指針，很明顯，T 型節(jié)點會有2 個指針指向同一個葉網(wǎng)格．每個最小公共邊有2 個指針指向相鄰的葉網(wǎng)格．如圖4 所示，網(wǎng)格Ⅱ、Ⅲ、Ⅳ和Ⅴ有共同的父網(wǎng)格Ⅰ，網(wǎng)格Ⅳ有4 個子網(wǎng)格Ⅵ、Ⅶ、Ⅷ和Ⅸ．網(wǎng)格Ⅰ和子網(wǎng)格Ⅳ有共同角點C．若網(wǎng)格Ⅱ、Ⅲ、Ⅴ和Ⅶ是葉網(wǎng)格，網(wǎng)格節(jié)點E 的4 個指針將分別指向這4 個網(wǎng)格．若網(wǎng)格Ⅸ和Ⅴ是葉網(wǎng)格，網(wǎng)格Ⅸ的右側邊1 是最小公共邊，則最小公共邊1 的2 個指針分別指向Ⅸ和Ⅴ．在遍歷四叉樹網(wǎng)結構時，采用鏈表連接同層的葉網(wǎng)格是十分方便的．

圖1 等級結構網(wǎng)格Fig.1 Hierarchical structure of grids

圖2 四叉樹網(wǎng)格結構示意Fig.2 Quadtree and quadtree subdivision

圖3 側邊、邊和角點Fig.3 Side，edge and vertex

為了確保唯一性，定義某一方位上層數(shù)不大于n的所有相鄰網(wǎng)格中尺寸最小者為n 層網(wǎng)格在該方位上的鄰居．例如，圖4 中網(wǎng)格Ⅴ左側與網(wǎng)格Ⅳ、Ⅸ和Ⅶ相接，但只有Ⅳ被認為是它的鄰居．相鄰葉網(wǎng)格尺寸之比不超過2 時為平衡四叉樹結構．圖5 給出平衡前、后兩種網(wǎng)格的對比，其中實線代表平衡前的網(wǎng)格，虛線代表通過平衡增加的網(wǎng)格．對矩形計算區(qū)域，可先將其劃分成大小相同的若干個0 層正方形網(wǎng)格，再對每個網(wǎng)格進行逐層細化，直到滿足分辨率要求．

圖4 四叉樹的等級結構Fig.4 Hierarchical structure of a quadtree

圖5 平衡的四叉樹網(wǎng)格結構Fig.5 A balanced quadtree

為離散方便，速度和壓力采用交錯定義，即壓力定義在葉網(wǎng)格中心，速度分量定義在最小公共邊上(圖6(a))．網(wǎng)格側邊上的速度可以通過所包含的最小公共邊上的速度遞歸得到，如圖6 中根網(wǎng)格右側邊上的速度u0可從葉網(wǎng)格右側邊上速度遞歸得到，即有u0＝0.5,u2＋0.25,u5＋0.125(u11＋u12)．

圖6 交錯四叉樹網(wǎng)格上變量的定義Fig.6 Definition of variables of a staggered quadtree mesh

3 四叉樹網(wǎng)格下N-S方程數(shù)值求解

3.1 對流項

新模型采用具有二階精度無條件穩(wěn)定的MacCormack 方法離散對流項，其基本原理如下所述.

第2 步，根據(jù)預測結果采用半拉格朗日方法進行反向運算得到原始時刻的校正值為

第3步，估計預測值的誤差e，e=進而獲得的最終結果=整個計算過程如圖7 所示．在上述過程中，誤差矯正可能導致數(shù)值不穩(wěn)定，可采用限制器來修正最終結果，即將最終結果限制在t?值范圍之內．

圖7 MacCormack方法的計算步驟示意Fig.7 Illustration of the MacCormack method

按照上述原理，在新模型中，對流項的計算過程如下：

(1) 計算t 時刻網(wǎng)格節(jié)點上的速度以及最小公共邊上切向速度；

(2) 按照式(8)計算t＋Δt 時刻最小公共邊上法向速度的預測值；

(3) 根據(jù)預測值，計算t＋Δt 時刻網(wǎng)格節(jié)點上的速度以及最小公共邊上的切向速度；

(4) 按照式(9)反向計算出t 時刻最小公共邊上法向速度的校正值；

(5) 計算誤差，獲得t＋Δt 時刻最小公共邊上法向速度的最終值．

網(wǎng)格節(jié)點上的速度可通過其周圍最大葉網(wǎng)格所在層的速度插值得到，T 型節(jié)點上的速度則通過所在側邊端點上的速度線性插值得到．最小公共邊的切向速度通過其端點上速度的中心插值得到．

3.2 黏性項的離散

在均勻網(wǎng)格中，黏性項一般可采用中心差分進行離散．然而，在樹型網(wǎng)格中，由于相鄰葉網(wǎng)格尺寸不一定相同，這種方法不便直接運用．本文提出了一種新方法，其計算細節(jié)如下．

(1) 葉網(wǎng)格中心的速度梯度通過對該網(wǎng)格側邊上的速度中心差分得到，如葉網(wǎng)格1 中心的正應力? u ?x 和 ? v ? y分別為和，葉網(wǎng)格2中心的正應力ux? ?和v y? ?分別為和速度則等于字母上標t、b、l 和r 分別表示網(wǎng)格上、下、左、右側邊，數(shù)字下標表示網(wǎng)格(圖8(a))．

(2) 網(wǎng)格節(jié)點上的速度梯度通過與該節(jié)點相鄰的最大葉網(wǎng)格所在層上速度的中心差分得到，而T型節(jié)點上的速度梯度通過所在側邊端點上的剪應力中心插值得到，如網(wǎng)格節(jié)點A 上速度梯度 u y? ? 和、? v ? x 可通過和節(jié)點B 上速度梯度可通過節(jié)點A 和C 上相應值的中心插值得到(圖8(a))．

(3) 根據(jù)式(3)，分別計算出葉網(wǎng)格中心和節(jié)點上的正應力和剪應力．

我很小很小的時候，在澡盆里洗澡，洗完了，澡盆被端走，地上有一個圓圓的水印，我就指著水印說：“太陽！太陽！”據(jù)說我當時這樣說的時候，是十分激動的。夏天，我赤著腳在地上走，腳上有水，地上就有腳印，我又指著腳印說：“小船！小船！”看來我小時候是有些想象力的，而我現(xiàn)在想象力要比那時糟得多。

(4) 若相鄰葉網(wǎng)格大小一致，則最小公共邊上的正壓力梯度通過對葉網(wǎng)格正壓力的中心差分獲得；否則，令定義在同一方位上與相同葉網(wǎng)格相鄰的兩個最小公共邊上正應力梯度相等，都等于上一級網(wǎng)格側邊上的正應力梯度，該值可通過其相鄰的大葉網(wǎng)格和小葉網(wǎng)格的父網(wǎng)格上正應力的中心差分獲得，而父網(wǎng)格上正應力通過其子網(wǎng)格上正應力的中心插值得到，如網(wǎng)格1 左側邊上的正壓力梯度為

網(wǎng)格10 和11 下側邊的正應力梯度與網(wǎng)格2 上側邊的正應力梯度相等，等于對網(wǎng)格5 和2 上正應力的中心差分(圖8(a))，即

其中

(5) 最小公共邊上的剪應力梯度通過其端點上剪應力的中心差分獲得，在同一方位上與相同葉網(wǎng)格相鄰的兩個最小公共邊上剪應力梯度也相等．

為了保證計算的穩(wěn)定性，時間步長取

式中：(Δx)min為最小網(wǎng)格尺寸；v 為運動黏性系數(shù)．

圖8 四叉樹中不同層的網(wǎng)格Fig.8 Meshes in different levels of quadtree

3.3 壓力泊松方程的離散

將式(7)寫成

對于相鄰葉單元大小一致的情況，定義在側邊上的壓力梯度可直接采用中心差分格式進行離散．以圖9 為例，網(wǎng)格1 上側邊的壓力梯度為

但是，對于相鄰葉單元大小不一致的情況，不能直接使用中心差分格式．這里采用一種基于SOM 的新方法．按此方法，定義在葉網(wǎng)格1 和2 之間的網(wǎng)格側邊上的壓力梯度為

式中Δ為大小葉網(wǎng)格尺寸的平均值．葉網(wǎng)格1 右側邊上的壓力梯度則通過對其包含的最小公共邊上的壓力梯度的加權平均得到，即

式中下標s 和L 表示小網(wǎng)格和大網(wǎng)格．為了保證計算結果具有二階精度，小葉網(wǎng)格2 左側邊上的壓力梯度采用與大葉網(wǎng)格1 右側邊上一致的壓力梯度，即

通過上述處理，壓力泊松方程的系數(shù)矩陣的對稱性也能得到保證，可以采用ICCG 算法求解．

圖9 離散泊松方程時變量的定義Fig.9 Definition of variables of discretized Poisson equation

4 網(wǎng)格的適時調整

為了保證計算精度和減少計算量，在計算過程中需要對網(wǎng)格尺度進行適時調整，方法是根據(jù)定義在網(wǎng)格節(jié)點上的渦度τ值的大小來進行判斷．渦度

在計算區(qū)域Ω內，若父網(wǎng)格C 中所有網(wǎng)格節(jié)點上的τ值都小于閾值τ0，則將其子網(wǎng)格進行合并；若葉網(wǎng)格C 上存在大于閾值τ0的渦度，則該網(wǎng)格需要細分．合并和細分的過程將分步進行，即按四叉樹的結構先從上往下對滿足條件的父網(wǎng)格進行合并，再從下往上對滿足條件的葉網(wǎng)格進行細化．調整完畢后，再將生成的四叉樹網(wǎng)格進行平衡化處理．網(wǎng)格細分時，網(wǎng)格側邊上新節(jié)點的速度通過側邊端點上速度中心插值得到，網(wǎng)格中心處新節(jié)點的速度則通過網(wǎng)格四周節(jié)點上速度線性插值得到，新生成的最小公共邊上的法向速度則通過其端點上相應速度線性插值得到．

5 算 例

5.1 泊松方程

本算例將討論新模型中泊松方程的求解精度．考慮泊松方程

其真解為

取尺寸為1×1 的計算區(qū)域，初始網(wǎng)格層數(shù)是4，初始最大網(wǎng)格尺寸是0.125×0.125，初始網(wǎng)格劃分如圖10 所示．分別用不同分辨率的網(wǎng)格進行計算，顯然網(wǎng)格分辨率越高計算結果越精確．根據(jù)計算結果可以獲得數(shù)值計算的誤差和精度(表1)，從中可以看出，格式的精度超過二階，其中，分辨率n2～m2表示最粗網(wǎng)格的分辨率是n2，最細網(wǎng)格的分辨率是m2．

圖10 四叉樹網(wǎng)格結構(分辨率82～642)Fig.10 Sketch of quadtree mesh(in 82—642)

表1 泊松方程的計算精度Tab.1 Computational accuracy of the Poisson equation

5.2 單渦流

本算例通過模擬區(qū)域?＝[0，π]×[0，π]內雷諾數(shù)為1 的單渦流動來分析新模型在四叉樹網(wǎng)格下的數(shù)值精度．單渦流的真解為

體積力的表達式為

初始網(wǎng)格劃分如圖11 所示，網(wǎng)格層數(shù)為4 層，最大網(wǎng)格為0.785,4×0.785,4．計算時長為π．表2 給出了最終時刻速度數(shù)值誤差(L2誤差為二階范數(shù)誤差，L∞為無窮范數(shù)誤差)和精度結果，可以看出新模型給出的速度具有一階以上精度．圖12 給出了單渦流流線的數(shù)值計算結果．

圖11 四叉樹網(wǎng)格結構(分辨率42～322)Fig.11 Sketch of quadtree mesh(in42—322)

圖12 單渦流流線數(shù)值計算結果Fig.12 Snapshot of the streamline of single vortex flow

表2 單渦流動速度的計算精度Tab.2 Computational accuracy of velocities in single vortex flow

5.3 方腔流

本算例對Ghia 的二維方腔流算例[12]進行模擬．計算區(qū)域大小為1×1(圖13)．計算區(qū)域上邊界為剪切邊界，剪切速度大小恒為1，其他邊則為無滑移邊界．雷諾數(shù)取為1,000，網(wǎng)格調整閾值τ0同樣設為0.04．初始時刻，整個計算區(qū)域被均勻地分成64×64個正方形網(wǎng)格．預設網(wǎng)格總層數(shù)為3．計算時長為50,s.

計算出的流線和網(wǎng)格變化如圖13 所示，可以看出，由于頂部邊界的剪切作用，附近的網(wǎng)格始終保持最細狀態(tài)；在整個流動發(fā)展過程中，部分區(qū)域網(wǎng)格隨渦度增加自動加密，又隨渦度減小而恢復初始分辨率．

圖13 不同時刻流線和網(wǎng)格分布(分辨率為642～2562)Fig.13 Snapshots of the streamline and meshes at different times(in 642—2562)

圖14給出了穩(wěn)定狀態(tài)時不同分辨率下獲得的中軸線上速度分布與Ghia 的計算結果(分辨率為1282)的對比情況，可以看出，網(wǎng)格分辨率為642的計算結果與Ghia 的計算結果偏差較大，而分辨率為2562的結果與Ghia 的計算結果最接近，自適應格分辨率為642～2562的結果與分辨率為2562的結果相差無幾．表3 給出了各類網(wǎng)格下網(wǎng)格數(shù)量、泊松方程求解矩陣階數(shù)和計算時間的比較，可見自適應網(wǎng)格下網(wǎng)格數(shù)量和消耗的計算時間都顯著降低．雖然分辨率為642～2562的自適應網(wǎng)格下系數(shù)矩陣的最大階數(shù)與分辨率為1282的均勻網(wǎng)格下系數(shù)矩陣的階數(shù)接近，但由于自適應下系數(shù)矩陣階數(shù)不斷變化，且相同階數(shù)下矩陣中非零元素的個數(shù)多于均勻網(wǎng)格下的個數(shù)，迭代步數(shù)會有所增加，計算時間也相應增加．此外，網(wǎng)格的自動調整也會消耗一定的計算時間．

圖14 中軸線上速度分布的計算結果與Ghia計算結果的比較Fig.14 Comparisons of velocities along central axes between numerical results and solution of Ghia

表3 自適應網(wǎng)格與均勻網(wǎng)格下方腔流計算所需網(wǎng)格與時間比較Tab.3 Comparisons of the number of meshes and the computational time in adaptive mesh and uniform mesh

5.4 圓柱繞流

計算區(qū)域如圖15 所示，無量綱長度為25，無量綱寬度為20．圓柱位于(5，10)處，其直徑D＝1，Ls＝10，Le＝5，Lr＝20．計算區(qū)域左邊界為入流邊界，入流速度大小為1．右邊界為出流邊界，上下兩個邊界為自由滑移邊界，圓柱表面為無滑移邊界．四周壓力邊界條件設為零梯度邊界條件．

圖15 計算區(qū)域示意Fig.15 Sketch of computational domain

本算例雷諾數(shù)取為200，網(wǎng)格調整閾值τ0仍為0.04．初始時刻整個計算區(qū)域分成100×80 個正方形網(wǎng)格．固壁邊界按階梯邊界近似處理．為了較好地刻畫圓柱，初始時刻對圓柱周圍一定范圍內的網(wǎng)格進行了加密處理，最大網(wǎng)格尺寸為0.025，預設最大層數(shù)為4．表4 給出了圓柱幾何特征值比與加密次數(shù)的關系．從圖16 可以看出，網(wǎng)格能夠很好地在渦旋處進行細分，最小網(wǎng)格變化與渦旋運動保持一致，在遠離渦街的區(qū)域，網(wǎng)格保持不變，這樣可使網(wǎng)格數(shù)量大大減少．圖17 給出了t＝62.5,s 時刻流線分布情況．

表4 圓柱幾何特征值比與加密次數(shù)的關系Tab.4 Variation of the geometrical characteristic value ratios of cylinder with the refinement times

圖16 不同時刻網(wǎng)格分布情況Fig.16 Snapshots of cell configuration at different times

圖17 t＝62.5,s 時的流線分布Fig.17 Snapshot of streamlines at t＝62.5,s

圖18給出了計算的拖曳力系數(shù)CD和升力系數(shù)CL隨時間變化的過程．拖曳力系數(shù)CD和升力系數(shù)CL計算式為

式中：FD和FL分別為拖曳力和升力；fρ為流體密度；∞u 為入流邊界的速度．在穩(wěn)定的渦街狀態(tài)，計算出的拖曳力系數(shù)最大值為1.390，最小值為1.233，平均值為1.301，與Wille[13]的試驗結果值1.3 十分接近．升力系數(shù)則在0.623 和－0.611 之間振蕩．

圖18 計算出的拖曳力系數(shù)和升力系數(shù)變化過程Fig.18 Numerical results of drag and lift coefficients of flow over a circle cylinder

6 結 論

(1) 通過數(shù)值試驗，證明所提出的自適應四叉樹網(wǎng)格下的水流模型中的泊松方程解可以達到二階精度，速度精度在一階以上．

(2) 方腔流算例的計算結果表明，自適應網(wǎng)格的網(wǎng)格數(shù)量比同分辨率均勻網(wǎng)格下的網(wǎng)格數(shù)量減少2/3以上，系數(shù)矩陣的階數(shù)要減少3/4，計算時間要減少近一半，而計算的中軸線上流速分布基本接近．

(3) 圓柱繞流算例中所得的拖曳力系數(shù)和升力系數(shù)與經(jīng)典結果基本一致．

(4) 所有算例的計算結果均表明，模型可以很好地根據(jù)設定的網(wǎng)格加密及合并控制指標完成網(wǎng)格的細分和合并．

［1］Berger M J，Oliger J. Adaptive mesh refinement for hyperbolic partial differential equations[J]. Journal of Computational Physics，1984，53(3)：484-512.

［2］Howell L H，Bell J B. An adaptive mesh projection method for viscous incompressible flow[J]. SIAM Journal on Scientific Computing，1997，18(4)：996-1013.

［3］Martin D F，Colella P，Graves D. A cell-centered adaptive projection method for the incompressible Navier-Stokes equations in three dimensions[J]. Journal of Computational Physics，2008，227(3)：1863-1886.

［4］De Berg M，van Kreveld M，Overmars M，et al. Computational Geometry Algorithms and Applications[M].Berlin：Springer，1999.

［5］Popinet S. Gerris：A tree-based adaptive solver for the incompressible Euler equations in complex geometries[J]. Journal of Computational Physics，2003，190(2)：572-600.

［6］Losasso F，Gibou F，F(xiàn)edkiw R. Simulating water and smoke with an octree data structure[C]//Proceedings of SIGGRAPH 2004 Conference， Association for Computing Machinery's Special Interest Group on Computer Graphics and Interactive Techniques ．San Antonio ，New York：ACM Press，2004：457-462.

［7］Losasso F，F(xiàn)edkiw R，Osher S. Spatially adaptive techniques for level set methods and incompressible flow[J]. Computers and Fluids，2006，35(10)：995-1010.

［8］Lipnikov K，Morel J，Shashkov M. Mimetic finite difference methods for diffusion equations on nonorthogonal non-conformal mesh[J]. Journal of Computational Physics，2004，199(2)：587-597.

［9］Min C H，Gibou F. A second order accurate projection method for the incompressible Navier-Stokes equations on non-graded adaptive grids[J]. Journal of Computational Physics，2006，219(2)：912-929.

［10］Selle A，F(xiàn)edkiw R，Kim Byungmoon，et al. An unconditionally stable MacCormack method[J]. SIAM Journal on Scientific Computing，2008，35(2)：350-371.

［11］Chorin A J. Numerical solution of the Navier-Stokes equations[J]. Mathematics of Computation ， 1968 ，22(104)：745-762.

［12］Ghia U，Ghia K N，Shin C T. High-resolutions for incompressible flow using the Navier-Stokes equations and multigrid method[J]. Journal of Computational Physics，1982，48(3)：387-411.

［13］Wille R. Karman vortex streets[J]. Advances in Applied Mechanics，1960，6：273-287.