微分求積模擬二維流體中流函數約束的施加方法研究

王 通， 何 濤,3， 曹曙陽

(1.上海師范大學 建筑工程學院，上海 201418；2. 同濟大學 土木工程防災國家重點實驗室，上海 200092；3. 伯明翰大學 工學院，伯明翰 B15 2TT，英國)

微分求積模擬二維流體中流函數約束的施加方法研究

王 通1， 何 濤1,3， 曹曙陽2

(1.上海師范大學 建筑工程學院，上海 201418；2. 同濟大學 土木工程防災國家重點實驗室，上海 200092；3. 伯明翰大學 工學院，伯明翰 B15 2TT，英國)

采用微分求積法數值求解流函數-渦度方程來模擬二維流體時會遇到流函數的超約束問題，即雖然流函數方程為二階偏微分方程，但在每個固體邊界上都存在兩個約束條件：一個Dirichlet條件和一個Neumann條件。以二維驅動方腔流動為例，對該問題進行深入分析，進而提出一種新的超約束處理方法，即在邊界渦度的計算中考慮Neumann條件，而僅將Dirichlet條件施加于流函數方程。數值結果顯示該方法可行，且計算效率較高。同時給出前人提出的單層法和雙層法進行比較。試算表明單層法對于網格數的奇偶性很敏感，不適于處理該問題。與雙層法對比后發現：該方法計算精度較高，且由于回避了超約束問題而更加方便于使用。

微分求積；流函數-渦度方程；方腔流；邊界條件；超約束

計算機硬件和軟件的迅速發展以及數值算法研究的不斷深入使得數值模擬成為流體力學領域繼理論分析和實驗研究之后的第三大研究手段，并越來越受到重視。相對于計算機的發展，人們往往對數值算法抱以更高的期望，這是因為算法的改進和發展會對數值模擬產生更深刻的影響。微分求積法(Differential Quadrature Method，DQM)自BELLMAN等[1]提出以來就備受關注，該方法不依賴于泛函和變分原理，是一種數學原理簡單、計算精度高、計算量少的高階數值算法，已被廣泛應用于流體力學、結構力學、熱傳導、生物科學等諸多領域[2-7]，并有望發展成為與有限差分、有限元等傳統低階算法同等的求解微分方程的強有力方法。本文研究內容是將DQM用以模擬二維驅動方腔流體，并著重討論流函數邊界約束的處理問題。

對于二維驅動方腔流體這種簡單流動前人已做過大量研究，并取得了豐富成果：GHIA等[8]采用CSI-MG方法求解流函數-渦度方程得到一系列雷諾數(最高達10 000)下的方腔流的解，最大網格數257×257；STRIZ等[9]采用DQM離散純流函數形式的流體控制方程，然后通過Newton-Raphson方法求解所得到的非線性方程組，得到較低雷諾數下的方腔流動；SHU等[10]采用DQM結合SIMPLE算法求解原參數形式(速度和壓力)的流體控制方程來模擬方腔流動，利用很少的網格就得到較理想的結果；BRUNEAU等[11]采用有限差分離散求解原參數形式的流體控制方程，得到了最高值為10 000等多個雷諾數下的方腔流數值解，最大網格數達2 048×2 048；MARCHI等[12]全面回顧了之前針對二維方腔流的研究，并采用有限體積法在1 024×1 024的網格上算得了被公認為迄今為止最準確的方腔流數值解。由上述研究可以看出，若想采用低階算法算得較準確的結果就需要劃分高密度的網格并花費大量的計算資源，這也是發展高精度算法的根本原因。

相對于原參數形式的流體控制方程，流函數-渦度方程在模擬二維流場中具有很大的優勢，它不僅求解變量少、自動滿足連續性條件，而且能將壓力計算和速度計算解耦，從而極大地簡化了計算。正因為如此，雖然它不能被拓展到三維流場，但在求解二維流場時頗受青睞。然而，當采用DQM數值求解流函數-渦度方程來模擬方腔流時卻遇到了流函數超約束的問題，即流函數的邊界條件數量多于流函數方程適定解所需要的約束數量[13]。SHU等[14]在采用類似方法模擬方腔內的自然對流時也遇到了相同的問題，他們給出了兩種處理方法，即單層法和雙層法(詳述見第3節)，并對兩種方法做了對比研究。

本文將基于二維驅動方腔流動深入分析流函數的超約束問題，進而提出一種新的處理方法，然后由數值計算驗證其可行性，最后對比本文方法與已有的單層和雙層法，并通過數值結果對比來證明本文方法的高效性。

1 微分求積法

微分求積法本質上是一種特殊的加權殘值法，等價于混合配點法[15-16]，基本思想是將函數在給定網點處的導數值近似用該導數自變量方向上全部網點處函數值的加權和表示。考慮一個光滑的二維函數f(x,y)，它在計算域內連續可微，根據微分求積原理有

(1a)

(1b)

根據插值基函數類型的不同，DQM又可以分成PDQ(基于多項式插值的DQM)和FDQ(基于傅里葉級數展開的DQM)兩類。PDQ一般適用于非周期性問題的求解；而FDQ則適用于周期和非周期性的問題，但對于周期性問題的效果更好。根據基函數選取的不同，權系數的計算也不同。本文研究僅采用PDQ，并采用拉格朗日插值函數作為基函數，由此，一階導數權系數的計算公式為

當i≠j時

(2a)

(2b)

當i=j時

(3a)

(3b)

對于二階甚至更高階導數的權系數，可采用下面的矩陣公式進行計算：

(4a)

(4b)

式中：[A(r)]為函數f(x,y)對x的r階偏導數所對應的權系數矩陣；[B(s)]為函數f(x,y)對y的s階偏導數所對應的權系數矩陣。

網格點的分布形式對計算結果的精度影響很大，一般情況下不等分網格比等分網格的計算精度高，所以本文將采用不等分網格進行計算。

2 問題描述及模型離散

二維驅動方腔流動的坐標系及速度邊界條件如圖1所示。流函數-渦度形式的流體控制方程為

(5a)

(5b)

式中：Ψ、ω、u、v、Re分別為流函數、渦度、x向速度、y向速度、雷諾數，且

(6a)

(6b)

圖1 二維驅動方腔流動圖示

將圖1所示的速度邊界條件代入式(6)就得到流函數的約束條件：

(7a)

(7b)

(7c)

由式(7)可以看出，流函數Ψ在每個邊界上都存在兩個邊界條件：一個Dirichlet條件和一個Neumann條件。而流函數方程(5a)卻僅僅是一個二階偏微分方程(只需每個邊界上有一個Dirichlet條件就能滿足適定解的要求)，這就是所謂的流函數的超約束問題。

(8a)

(8b)

根據微分求積原理將式(5)在空間上離散后得到

(9a)

(9b)

采用同樣的方法也可以將式(6)～(8)進行空間離散，這里不再贅述。

3 流函數約束施加方法

SHU等[14]采用微分求積法模擬二維方腔內的自然對流時也遇到類似的超約束問題，他們提出了單層法和雙層法等兩種處理方法：所謂單層法就是將兩個約束條件全部施加在處于邊界上的那層網格點上，具體做法是在計算流函數二階導數的權系數時嵌入Neumann條件；而雙層法就是將兩個邊界約束分別施加于不同的兩層網格點上，具體做法是將Dirichlet條件施加于邊界點上，而將Neumann條件施加于緊鄰邊界的那層網格點上。單層法和雙層法都是把流函數的兩個約束條件全部施加于流函數方程，雖然流函數方程僅僅是二階微分方程。然而從微分方程的理論角度講，這樣做可能是不合理和不準確的。我們在試算中發現，單層法對于網格數的奇偶性很敏感，不適于處理本文問題(詳細分析見第5節)，這在一定層度上也印證了上述判斷。

現在我們換一個角度來分析流函數的超約束問題：首先，把式(5a)代入式(5b)，消去渦度變量后得到關于流函數的四階偏微分方程，即純流函數形式的流體控制方程，而流函數在每個邊界上恰好有兩個約束條件，滿足適定解的要求，所以不存在超約束的問題；其次，流函數-渦度方程與純流函數方程是等價的，流函數的約束條件也相同，僅僅由于形式不同，流函數-渦度方程出現了超約束問題；最后，由于流函數方程與渦度方程是耦合的，若其中一個方程受到某種約束，另一方程也應受到相同約束。基于上述分析，我們給出一種新的處理方法：將Neumann條件施加于渦度方程，即在采用式(8)的微分求積離散式來計算邊界渦度時考慮Neumann條件，而僅將Dirichlet條件施加于流函數方程，從而回避了流函數的超約束問題。

4 方法驗證

本節將通過對二維驅動方腔流體的數值求解來檢驗新方法的可行性。求解方法與SHU等[14]所用方法類似：采用LU分解來求解微分求積離散所得到的線性代數方程組，對于渦度方程在時間域采用四階Runge-Kutta法進行計算，網格采用不等分形式的Chebyshev-Gauss-Lobatto節點。流場變量初始值都設置為0，迭代計算當達到下式標準時結束，

(10)

式中：Resij表示渦度方程在網格點(xi,yj)處的殘差。計算平臺是ThinkPad SL400，雙核主頻分別為2.2 GHz，內存2 GB RAM。

圖2和圖3所示為Re=100、400、1 000，對應網格分別為19×19、25×25、37×37時，由本文方法計算所得x向速度u沿豎向中軸線以及y向速度v沿水平中軸線的分布，同時還給出已有結果進行對比。比較發現，本文方法的計算結果與已有結果吻合的很好。為進一步驗證該方法的精度和效率，表1和表2分別給出Re=1 000時由本文方法在不同網格密度下計算所得u沿豎向中軸線以及v沿水平中軸線的離散數值，并給出MARCHI等的經典結果進行對比。由表中數據可以看出，當網格密度較小時，計算結果偏差較大，但隨著網格密度的加大，計算結果迅速靠近參照值，而當網格密度增大到一定程度后，再增大網格密度對計算結果的改進不大。以上結果充分證明了該方法的可行性。

圖2 u沿豎向中軸線的分布

圖3 v沿水平中軸線的分布

表1 本文方法算得速度u在豎向中軸線上的數值(Re=1 000)

表2 本文方法算得速度v在水平中軸線上的數值(Re=1 000)

5 對比分析

由第3節可知，單層法和雙層法都是將流函數的所有約束條件全部施加于流函數方程，操作較為復雜；而本文方法將Dirichlet條件和Neumann條件分別施加于流函數方程和渦度方程，從根本上回避了超約束問題，而且實施方便。另外，單層和雙層法將Neumann條件強加于流函數方程的做法從理論上講是不合理的，可能會造成計算發散或數值奇異。研究中我們發現，雙層法的計算效果較好，而單層法對于x向和y向網格點數的奇偶性很敏感，對此我們嘗試了以下四種組合：N和M都是奇數；N為奇數，M為偶數；N為偶數，M為奇數；N和M都是偶數。試算發現：前兩種組合下，計算發散；后兩種組合下，只有當時間步長很小，并且對方腔頂部兩角點附近的速度u做一些光滑處理后才能得到收斂解[17]，但數值精度較差。這或許是強加所有約束條件于流函數方程的不合理性造成的，關于這一點本文不做進一步的深入研究。由此可見，單層法不適于處理本文問題，所以在下面的對比研究中，我們僅考慮雙層法。

由上述分析可知本文方法較雙層法操作更加方便，不僅如此，我們將從計算精度上進一步證明該方法的優越性。表3和表4所示為Re=1 000時分別由本文方法和雙層法在不同網格密度下計算所得u沿豎向中軸線以及v沿水平中軸線的離散數值，并給出MARCHI等的經典結果進行比較。可以看出，隨著網格密度的加大，兩種方法的計算精度都迅速提高，但在相同網格密度下，本文方法的精度明顯高于雙層法。

表3 本文方法與雙層法算得u在豎向中軸線上的對比(Re=1 000)

表4 本文方法與雙層法算得v在水平中軸線上的對比(Re=1 000)

接下來分析兩種方法的計算效率問題。由第3節可知，雙層法是把兩個流函數約束分別施加在兩層網格上，那么最終需要求解的離散變量個數為(N-4)× (M-4)，而本文方法需要求解(N-2)×(M-2)個離散變量，計算量明顯高于雙層法。另外，采用雙層法時所允許的最大時間步也明顯大于本文方法。這兩個因素決定了雙層法得計算效率明顯高于本文方法，如表5所示。但在現有的計算機水平下，本文方法的計算效率是完全可以接受的，并且該方法回避了超約束問題，操作相對更加方便。

表5 本文方法與雙層法計算效率對比(Re=1 000)

對于本文所討論的問題，若要嚴格反映流體邊界條件，應該將流體控制方程升階到四階的純流函數方程進行求解，類似工作已由STRIZ等開展，但本文目的是研究如何采用微分求積法數值求解流函數-渦度函數來模擬二維方腔流動，而為了驗證本文方法的準確性，我們直接將本文數值結果與前人經典結果進行對比，并未與STRIZ等的結果比較。實際上，當采用微分求積法數值求解純流函數方程來模擬二維方腔流動時同樣會遇到在一個邊界點上存在兩個邊界條件的問題，但并非超約束，此時可采用WANG等[18]的嵌入法來處理邊界條件，另外若在時域上也采用本文所選擇的時間推進格式，則上述微分求積數值求解純流函數方程來模擬方腔流的方法等價于本文微分求積數值求解流函數-渦度方程來模擬方腔流的方法。但這不是本文所討論的問題，所以針對這一問題，我們不再深入展開。

6 結論

基于二維驅動方腔流體，本文深入探討了微分求積數值求解流函數-渦度方程中所遇到的流函數的超約束問題，并提出一種新的處理方法，即把Dirichlet條件和Neumann條件分別施加于流函數方程和渦度方程，從而回避了超約束問題。同時將本文方法與已有的單層法和雙層法進行對比。研究發現，單層法對于網格數的奇偶性很敏感，不適于處理本文問題。對比本文方法與雙層法發現：雖然本文方法的計算效率相對較低，但計算精度更高，且由于回避了超約束問題而更加方便于使用。

[1] BELLMAN R E, CASTI J. Differential quadrature and long-term integration[J]. Journal of Mathematical Analysis and Applications, 1971, 34(2): 235-238.

[2] BERT C W, MALIK M. Differential quadrature method in computational mechanics: a review[J]. Applied Mechanics Reviews, 1996, 49(1): 1-28.

[3] SHU C. Differential quadrature and its application in engineering[M]. London: Springer, 2000.

[4] 程昌鈞, 朱正佑. 微分求積方法及其在力學應用中的若干新進展[J]. 上海大學學報(自然科學版), 2009, 15(6): 551-559.

CHENG Changjun, ZHU Zhengyou. Recent advances in differential quadrature method with applications to mechanics[J]. Journal of Shanghai University (Natural Science Edition), 2009, 15(6): 551-559.

[5] 李鴻晶, 王通. 結構地震反應分析的逐步微分積分方法[J]. 力學學報, 2011, 43(2): 430-435.

LI Hongjing, WANG Tong. A time-stepping method of seismic response analysis for structures using differential quadrature rule[J]. Chinese Journal of Theoretical and Applied Mechanics, 2011, 43(2): 430-435.

[6] 武吉梅, 陳媛, 王硯, 等. 基于微分求積法的印刷運動薄膜動力穩定性分析[J]. 振動與沖擊, 2015, 34(20): 57-60.

WU Jimei, CHEN Yuan, WANG Yan, et al. Study of kinetic stability for membrane based on differential quadrature method[J]. Journal of Vibration and Shock, 2015, 34(20): 57-60.

[7] 董宇, 楊翊仁, 魯麗. 基于微分求積法的軸向流作用下二維板復雜響應研究[J]. 振動與沖擊, 2015, 34(6): 46-51.

DONG Yu，YANG Yiren，LU Li. On the complicated responses of a two-dimensional plate in axial liquid flow by differential quadrature method[J]. Journal of Vibration and Shock, 2015, 34(6): 46-51.

[8] GHIA U, GHIA K N, SHIN C T. High-Resolutions for incompressible flow using the Navier-Stokes equations and a multigrid method[J]. Journal of Computational Physics, 1982, 48(3): 387-411.

[9] STRIZ A G, CHEN W. Application of the differential quadrature method to the driven cavity flow problem[J]. International Journal of Non-Linear Mechanics, 1994, 29(5): 665-670.

[10] SHU C, WANG L, CHEW Y T. Comparative studies of three approaches for GDQ computation of incompressible Navier-Stokes equations in primitive variable form[J]. International Journal of Computational Fluid Dynamics, 2004, 18(5): 401-412.

[11] BRUNEAU C H, SAAD M. The 2D lid-driven cavity problem revisited[J]. Computers & Fluids, 2006, 35(3): 326-348.

[12] MARCHI C H, SUERO R, ARAKI L K. The lid-driven square cavity flow: numerical solution with a 1 024×1 024 grid[J]. Journal of the Brazilian Society of Mechanical Sciences and Engineering, 2009, 31(3): 186-198.

[13] DAVIS M B, CAREY G F. Iterative solution of the stream function-vorticity equations using a multigrid solver with finite elements[J]. Communications in Numerical Methods in Engineering, 1993, 9(7): 587-594.

[14] SHU C, XUE H. Comparison of two approaches for implementing stream function boundary conditions in DQ simulation of natural convection in a square cavity[J]. International Journal of Heat and Fluid Flow, 1998, 19(1): 59-68.

[15] 劉劍, 王鑫偉. 基于微分求積法的逐步積分法與常用時間積分法的比較[J]. 力學季刊, 2008, 29(2): 304-309.

LIU Jian, WANG Xinwei. Comparisons of successive integration method based on differential quadrature with commonly used time integration schemes[J]. Chinese Quarterly of Mechanics, 2008, 29(2): 304-309.

[16] 汪芳宗, 廖小兵, 謝雄. 微分求積法的特性及其改進[J]. 計算力學學報, 2015, 32(6): 765-771.

WANG Fangzong, LIAO Xiaobing, XIE Xiong. Characteristic of the differential quadrature method and its improvement[J]. Chinese Journal of Computational Mechanics, 2015, 32(6): 765-771.

[17] KU H C, TAYLOR T D, HIRSH R S. Pseudospectral methods for solution of the incompressible Navier-Stokes equations[J]. Computers & Fluids, 1987, 15(2): 195-214.

[18] WANG X, LIU F, WANG X, et al. New approaches in application of differential quadrature method to fourth-order differential equations[J]. Communications in Numerical Methods in Engineering, 2005, 21(2): 61-71.

Methods on applying stream-function restraints in differential quadrature modelling of two-dimensional flow

WANG Tong1, HE Tao1,3, CAO Shuyang2

(1. College of Civil Engineering, Shanghai Normal University, Shanghai 201418, China;2. State Key Laboratory for Disaster Reduction in Civil Engineering, Tongji University, Shanghai 200092, China;3. School of Engineering, University of Birmingham, Birmingham B15 2TT, UK)

The 2D lid-driven cavity flow was simulated by applying the differential quadrature method to solve the stream function-vorticity equations. There were two boundary conditions, one Dirichlet and one Neumann, for the stream function equation at each solid boundary though the stream function equation was just second order. Analysis on this over- specified problem was carried out, based on which a new applying method was proposed: the Neumann condition was considered in calculating the vorticity at the boundary while only the Dirichlet condition was applied in the stream function equation. Validity of this method was verified by comparing its numerical results with benchmark data. Two other existing methods, the one-layer approach and the two-layer approach were shown as contrasts. Trial calculations indicate that the one-layer approach is sensitive to the parity of grid numbers and is not suitable for the present problem. Comparisons between the new method and the two-layer approach show that the former is not only more accurate but also more convenient to be used in practice for avoiding the over-specified problem.

differential quadrature method; stream function-vorticity equations; cavity flow; boundary condition; over-specified

國家自然科學基金青年基金項目(51508333)

2016-04-12 修改稿收到日期：2016-08-11

王通 男，博士，講師，1981年2月生

TV131；O302

A

10.13465/j.cnki.jvs.2017.08.027