駱佳玲 李偉忠 王文全

摘 要：本文利用投影法和浸入邊界法使用的固定笛卡爾網(wǎng)格技術(shù)，三維泊松方程（poisson equation）采用有限七點差分格式離散，結(jié)合迭代法和直接法將三維的數(shù)值求解問題簡化為二維問題，然后利用二維的直接數(shù)值解法進行快速的迭代求解。采用VC++編寫數(shù)值計算代碼，通過三維數(shù)值算例驗證了三維壓力Poisson方程求解數(shù)值方法的有效性和求解精度，并以三維槽道流場為基準數(shù)值算例，驗證利用投影法的三維流場數(shù)值計算結(jié)果的可靠性和適用性。

關(guān)鍵詞：投影法; 三維Poisson方程; 三維流場

中圖分類號：O35

文獻標志碼：A

投影法是數(shù)值求解N-S（navier stokes）方程的一種十分有效的方法，最初是由CHORIN[1]在1968年提出的，在之后的發(fā)展過程中不斷得到改進，得到一系列的改進后的投影格式，使其求解精度得到了極大的提高，例如，常懷見[2]構(gòu)造了一種基于非結(jié)構(gòu)網(wǎng)格的求解非穩(wěn)態(tài)流動的高精度數(shù)值算法，時間和空間的收斂精度達到了二階以上。胡銳鋒等[3]提出了具有四階空間精度的不可壓縮流動的投影格式。浸入邊界法中采用固定的笛卡爾網(wǎng)格，且在求解中不需要更新網(wǎng)格，在變形大的流固耦合問題計算中具有很大的優(yōu)勢，所以投影法和浸入邊界法的結(jié)合受到了很多研究人員的關(guān)注。王文全等[4]在傳統(tǒng)的投影法基礎(chǔ)上發(fā)展了一種投影浸入邊界法，提高了計算的精度。TIAN等[5]利用有限元法和投影浸入邊界法三維數(shù)值模擬了大變形的流固耦合問題，研究了其在生物流體力學中的具體應(yīng)用。李宇航[6]在投影法與直接力浸入邊界法結(jié)合的基礎(chǔ)上修正體積力，研究了仿飛蛇滑翔平板的空氣動力學機理。

在N-S方程求解中，壓力泊松方程（Poisson equation）的數(shù)值求解是難點之一，常用的數(shù)值求解方法有迭代法和直接法。二維Poisson方程具有五點差分格式[7]，可利用直接法進行求解，三維Poisson方程的差分格式在有些文獻[8]中進行了一般性的介紹，文獻[9]中也給出了具體的三維Poisson方程的七點差分格式， 但是其精度及其實用性沒有得到驗證。二維Poisson方程可利用直接法進行快速穩(wěn)定的求解[10]，但是三維問題得到的線性方程組相對于二維情況變得很大，這對求解方程組的數(shù)值方法和計算機的存儲都具有很大的挑戰(zhàn)。

為此，本文運用七點差分格式并結(jié)合迭代法原理，將三維的數(shù)值求解簡化為二維的數(shù)值求解，結(jié)合方程組求解的迭代求解法和直接求解法，建立投影法的三維數(shù)學模型和數(shù)值求解策略，對三維問題進行快速穩(wěn)定的求解。采用VC++編寫投影法的數(shù)值計算代碼，通過求解三維數(shù)值算例驗證了三維壓力Poisson方程數(shù)值方法的求解精度和準確性，并以三維槽道流場為基準數(shù)值算例，驗證了三維壓力Poisson方程數(shù)值計算結(jié)果的可靠性和適用性。

1 數(shù)值計算方法

11 控制方程

流體為黏性不可壓縮流體時，張量形式的控制方程表示如下[11]，

ui=0（1）

uit+uj·ui=-p+1ReΔui（2）

式中，ui，uj為速度矢量;p為壓強;Re為流動雷諾數(shù)。

12 時間推進和對流擴散項的離散方法

利用投影法對方程（1），（2）進行分步求解，基本步驟如下：

1）速度預算步。不考慮壓強對時間離散求解得到預估速度u*i，

u*i=uni-Δt（C+V）（3）

式中，C表示對流項;V表示黏性擴散項。在求解方程（3）時，對C和V的數(shù)值離散化，參照文獻[10]中二維的數(shù)值離散化，將其直接推廣至三維問題。

2）速度校正步。考慮壓強項和瞬態(tài)時間項離散求解得校正速度u′i，

u′i=u*i-Δtpn+1xi（4）

式中，pn+1表示下一時間步的壓力值，可由壓力泊松方程求出，

2pn+1x2i=xiu*iΔt=fi，j，k（5）

求解（5）式時使用Neumann邊界條件，

pn+1xi=-u′i-u*iΔt=Ci（6）

13 三維泊松方程的數(shù)值求解方法

在整個數(shù)值求解計算過程中，求解三維壓力Poisson方程是關(guān)鍵也是難點。我們考慮一長方體區(qū)域，在x、y、z方向上分別劃分為M、N、P等份，三個方向上的網(wǎng)格間距分別為Δx、Δy、Δz。利用標準的七點差分格式空間離散三維Poisson方程，即：

pi+1，j，k-2pi，j，k+pi-1，j，kΔx2+pi，j+1，k-2pi，j，k+pi，j-1，kΔy2+

pi，j，k+1-2pi，j，k+pi，j，k-1Δz2=fi，j，k（7）

將包含z方向變化的數(shù)值項視為已知數(shù)，移至等式右邊整理得，

pi-1，j，kΔx2+pi，j-1，kΔy2-21Δx2+1Δy2+1Δz2pi，j，k+pi，j+1，kΔy2+pi+1，j，kΔx2=fi，j，k-pi，j，k+1+pi，j，k-1（Δz）2=Ci，j，k（8）

第一時間步內(nèi)分別在各個xy平面上進行直接迭代求解，后續(xù)時間步內(nèi)不用迭代法，右端項的壓力項采用上一時間步的值用直接法直接求解，邊界內(nèi)部采用中心差分格式離散，詳細的邊界處理步驟見文獻[10，12-13]。

離散化后，通過整理可以得到與二維同階的大型稀疏方程組，表示如下：

AP=D （9）

式中，

P=[P0，j，k，P1，j，k，…，Pi，j，k，…，PM-1，j，k，PM，j，k]T

Pi，j，k=[pi，0，k，pi，1，k，…，pi，j，k，…，pi，N-1，k，pi，N，k]T，i=0，1，…，M

D=[D0，j，k，D1，j，k，…，Di，j，k，…，DM-1，j，k，DM，j，k]T

D0，j，k=C0，0，k+2aΔxC1+2bΔyC2C0，1，k+2aΔxC1C0，j，k+2aΔxC1C0，N-1，k+2aΔxC1C0，N，k+2aΔxC1-2bΔyC2

Di，j，k=Ci，0，k+2bΔyC2Ci，1，kCi，j，kCi，N-1，kCi，N，k-2bΔyC2，i=1，2，…，M-1

DM，j，k=CM，0，k-2aΔxC1+2bΔyC2CM，0，k-2aΔxC1CM，j，k-2aΔxC1CM，N-1，k-2aΔxC1CM，0，k-2aΔxC1-2bΔyC2

A=B 2C

C Β C

C B C

C B C

2C B（M+1）×（M+1）

B=d 2b

b d b

b d b

b d b

2b d（N+1）×（N+1）

C=a a

a

a

a（N+1）×（N+1）

a=1Δx2，b=1Δy2，c=1Δz2，d=-21Δx2+1Δy2+1Δz2

在每一次對大型稀疏線性方程組的求解時，k為常數(shù)，是二維求解問題，k值在0，P上變化。流場在時間上推進過程中，只需要在t=0～Δt的第一步求解時進行迭代求解，可以減小因多次迭代求解產(chǎn)生的誤差累積，通過解+1次線性方程組，三維問題就可得到求解。

2 數(shù)值算例

21 三維泊松方程的數(shù)值求解驗證

為了驗證13節(jié)提出的三維Poisson方程數(shù)值求解方法的有效性和準確性，考慮Ω={（x，y，z）|0≤x≤1，0≤y≤1，0≤z≤1}內(nèi)的數(shù)值邊值問題，函數(shù)的表達式及其邊界條件表示如下：

Δu=-3π2sin（xπ）cos（yπ+π/3）cos（zπ+π/4）

u|x=0=0

u|x=1=0

u/y|y=0=-πsin（xπ）sin（π/3）cos（zπ+π/4）

u/y|y=1=πsin（xπ）sin（π/3）cos（zπ+π/4）

u/z|z=0=-πsin（xπ）cos（yπ+π/3）sin（π/4）

u/z|z=1=πsin（xπ）cos（yπ+π/3）sin（π/4） （10）

函數(shù)u（x，y，z）的解析表達式如下，

u（x，y，z）=sin（xπ）cos（yπ+π/3）cos（zπ+π/4）（11）

如圖1所示為x=05，z=0面上不同網(wǎng)格劃分下計算精度的比較，網(wǎng)格數(shù)劃分大于50×50×50數(shù)值結(jié)果與精確解相一致。精確解與數(shù)值解以及它們的相對誤差見表1，可以發(fā)現(xiàn)隨著網(wǎng)格數(shù)量的增加，計算結(jié)果與精確解的誤差減小。

22 三維槽道流動

計算域和邊界條件如圖2所示，左面為速度入口，使用Dirichlet壓力邊界條件[4]，右面為出口，上下左右四個面為固體壁面，使用Neumann壓力邊界條件[4]，網(wǎng)格步長為Δx=Δy=Δz=0025，時間步長為Δt=0001。

圖3—圖6為Re=10，網(wǎng)格劃分50×50×50時流場中x，y，z方向的速度等值線圖及縱向各截面的壓強等值線圖。

由圖3—圖5可以看出，槽道中x方向的流動速度較大，所以x方向的速度分布在橫向各截面上都呈現(xiàn)對稱分布，中部速度為最大值約等于10，槽道壁面上速度為零，滿足邊界無滑移條件;槽道中y，z方向的流動速度較小，從圖4可以發(fā)現(xiàn)槽道下游右側(cè)y方向的速度較大，由圖5中可看出，z方向的速度變化是上游右側(cè)的速度大于下游右側(cè)的速度，隨著深度的變化，上游右側(cè)的速度變?yōu)樾∮谙掠斡覀?cè)的速度;整體上通過觀察比較圖4和圖5可以發(fā)現(xiàn)，離入口截面越遠的x方向下游y，z方向的速度較大，能清晰明顯的看出流場的三維效應(yīng)。圖6為縱向各截面壓強等直線圖，從圖中可看出，上游下部到下游上部壓強值從大到小呈梯度變化，在y方向上越往下游，這種壓強的梯度變化越明顯，壓強值的這種變化規(guī)律主要是由于流體速度的變化的影響，因為從圖4和圖5中可看出在流場中下游右側(cè)的三維效應(yīng)比較明顯，速度變化率大，壓強較小，符合三維槽道流場的真實物理變化規(guī)律。

3 結(jié)論

利用程序設(shè)計語言VC++編寫三維投影法的數(shù)值計算程序，并對具有精確解的三維Poisson方程和三維純流體的槽道進行數(shù)值模擬計算，結(jié)論如下：

1）迭代法比較容易實現(xiàn)，三維數(shù)值問題求解的大型方程組維數(shù)遠超出二維情況，直接數(shù)值求解和計算機存儲變得更加困難，據(jù)此結(jié)合迭代法原理將三維問題化為二維問題，利用二維問題大型稀疏線性方程組的快速求解和數(shù)據(jù)結(jié)構(gòu)優(yōu)化存儲方法，既能減小存儲空間，又能實現(xiàn)快速準確求解三維復雜流場問題，且由迭代法引起的計算誤差也非常小，求解保持了良好的數(shù)值穩(wěn)定性。

2）對壓力Poisson方程進行數(shù)值模擬，驗證了數(shù)值算法的網(wǎng)格依賴性和迭代誤差的可控性，通過比較計算結(jié)果與解析解，驗證了該數(shù)值算法的有效性和可靠性。

3）以純流體的三維槽道流場為基準數(shù)值算例，得到的計算結(jié)果與真實的物理規(guī)律現(xiàn)象基本相符，驗證了基于投影法的三維流場數(shù)值求解結(jié)果的可行性。

參考文獻：

[1]CHORIN A J. Numerical solution of the Navier-Stokes equations[J]. Mathemaics of Computation，1968，22（104）：745-762.

[2] 常懷見. 基于非結(jié)構(gòu)網(wǎng)格的非穩(wěn)態(tài)流動與換熱的高精度算法研究[D]. 上海： 華東理工大學， 2015.

[3] 胡銳鋒， 王萍， 王麗敏. 完全交錯網(wǎng)格中基于四階緊致差分格式的不可壓縮流動精確投影方法研究[C]// 第九屆全國流體力學學術(shù)會議. 南京： 中國力學學會， 2016.

[4] 王文全， 李偉忠， 閆妍，等. 一種投影浸入邊界法的快速直接求解方法[J]. 計算力學學報， 2014， 31（6）：775-780.

[5] TIAN F B， DAI H， LUO H， et al. Fluid-structure interaction involving large deformations： 3D simulations and applications to biological systems[J]. Journal of Computational Physics， 2014， 258：451-469.

[6] 李宇航. 仿飛蛇滑翔平板空氣動力學性能的數(shù)值研究[D]. 北京： 中國科學院大學，2018.

[7] 陸金甫， 關(guān)治. 偏微分方程數(shù)值解法[M]. 2版. 北京： 清華大學出版社， 2004 .

[8] 梁志輝， 李文杰. 三維球形域泊松方程的差分方程中自然邊界條件的處理[J]. 內(nèi)蒙古民族大學學報，2009， 24（1）：12-14.

[9] 豆桂芳， 吳振遠， 杜艷林. 三維泊松方程的七點差分格式[J].工程地球物理學報，2009，6（6）：802-805.

[10]梁林. 基于投影浸入邊界法的剛體-流體相互作用研究[D]. 昆明： 昆明理工大學，2013.

[11]查德維克， PETER. 連續(xù)介質(zhì)力學：簡明理論和例題[M]. 天津： 天津大學出版社， 1992.

[12]胡建偉， 湯懷民. 微分方程數(shù)值解法[M]. 2版. 北京： 科學出版社，2007.

[13]魯晶津， 吳小平， KLASUS S. 三維泊松方程數(shù)值模擬的多重網(wǎng)格方法[J]. 地球物理學進展，2009，24（1）：154-158.

（責任編輯：于慧梅）

A Three-dimensional Flow Field Numerical Solution

Based on Method of Projection

LUO Jialing1， LI Weizhong2，WANG Wenquan*1，3

（1.Faculty of Civil Engineering and Mechanics， Kunming University of Science and Technology， Kunming 650500，China; 2.Faculty of Metallurgical and Energy Engineering， Kunming University of Science and Technology， ，Kunming 650500，China;3. College of Water Resource & Hydropower， Sichuan University， Chengdu 610065，China）

Abstract：

Projected method and fixed Cartesian grid technology of immersed bounday method were used in the numerical solution. The three-dimensional poisson equation was discretized by limited difference scheme with seven grid points. The iterative method and direct method were used to transform the three-dimensional problem into two-dimensional problem. Then the two-dimensional direct numerical method was used for rapid iterative solution. VC++ language was used to write the numerical solution program.First， a three dimensional numerical example was used to verify the accuracy and precision of the numerical method for solving the three dimensional pressure poisson equation.And then a three dimensions channel flow was solved as a benchmark to verify reliability and applicability of results of 3D flow compute using projected method.

Key words：

projected method; three-dimensional Poisson equation; three-dimensional flow field