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

王 通， 何 濤,3， 曹曙陽

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

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

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

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

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

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

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

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

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

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

1 微分求積法

微分求積法本質(zhì)上是一種特殊的加權(quán)殘值法，等價于混合配點法[15-16]，基本思想是將函數(shù)在給定網(wǎng)點處的導(dǎo)數(shù)值近似用該導(dǎo)數(shù)自變量方向上全部網(wǎng)點處函數(shù)值的加權(quán)和表示?？紤]一個光滑的二維函數(shù)f(x,y)，它在計算域內(nèi)連續(xù)可微，根據(jù)微分求積原理有

(1a)

(1b)

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

當i≠j時

(2a)

(2b)

當i=j時

(3a)

(3b)

對于二階甚至更高階導(dǎo)數(shù)的權(quán)系數(shù)，可采用下面的矩陣公式進行計算：

(4a)

(4b)

式中：[A(r)]為函數(shù)f(x,y)對x的r階偏導(dǎo)數(shù)所對應(yīng)的權(quán)系數(shù)矩陣；[B(s)]為函數(shù)f(x,y)對y的s階偏導(dǎo)數(shù)所對應(yīng)的權(quán)系數(shù)矩陣。

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

2 問題描述及模型離散

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

(5a)

(5b)

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

(6a)

(6b)

圖1 二維驅(qū)動方腔流動圖示

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

(7a)

(7b)

(7c)

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

(8a)

(8b)

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

(9a)

(9b)

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

3 流函數(shù)約束施加方法

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

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

4 方法驗證

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

(10)

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

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

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

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

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

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

5 對比分析

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

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

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

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

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

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

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

6 結(jié)論

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

[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] 程昌鈞, 朱正佑. 微分求積方法及其在力學(xué)應(yīng)用中的若干新進展[J]. 上海大學(xué)學(xué)報(自然科學(xué)版), 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] 李鴻晶, 王通. 結(jié)構(gòu)地震反應(yīng)分析的逐步微分積分方法[J]. 力學(xué)學(xué)報, 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] 武吉梅, 陳媛, 王硯, 等. 基于微分求積法的印刷運動薄膜動力穩(wěn)定性分析[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] 董宇, 楊翊仁, 魯麗. 基于微分求積法的軸向流作用下二維板復(fù)雜響應(yīng)研究[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]. 力學(xué)季刊, 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]. 計算力學(xué)學(xué)報, 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

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

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

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

TV131；O302

A

10.13465/j.cnki.jvs.2017.08.027