駱佳玲 李偉忠 王文全

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

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

中圖分類號(hào)：O35

文獻(xiàn)標(biāo)志碼：A

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

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

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

1 數(shù)值計(jì)算方法

11 控制方程

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

ui=0（1）

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

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

12 時(shí)間推進(jìn)和對(duì)流擴(kuò)散項(xiàng)的離散方法

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

1）速度預(yù)算步。不考慮壓強(qiáng)對(duì)時(shí)間離散求解得到預(yù)估速度u*i，

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

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

2）速度校正步?？紤]壓強(qiáng)項(xiàng)和瞬態(tài)時(shí)間項(xiàng)離散求解得校正速度u′i，

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

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

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

求解（5）式時(shí)使用Neumann邊界條件，

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

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

在整個(gè)數(shù)值求解計(jì)算過程中，求解三維壓力Poisson方程是關(guān)鍵也是難點(diǎn)。我們考慮一長(zhǎng)方體區(qū)域，在x、y、z方向上分別劃分為M、N、P等份，三個(gè)方向上的網(wǎng)格間距分別為Δx、Δy、Δz。利用標(biāo)準(zhǔn)的七點(diǎn)差分格式空間離散三維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ù)值項(xiàng)視為已知數(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）

第一時(shí)間步內(nèi)分別在各個(gè)xy平面上進(jìn)行直接迭代求解，后續(xù)時(shí)間步內(nèi)不用迭代法，右端項(xiàng)的壓力項(xiàng)采用上一時(shí)間步的值用直接法直接求解，邊界內(nèi)部采用中心差分格式離散，詳細(xì)的邊界處理步驟見文獻(xiàn)[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

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

2 數(shù)值算例

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

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

Δ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）的解析表達(dá)式如下，

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

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

22 三維槽道流動(dòng)

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

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

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

3 結(jié)論

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

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

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

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

參考文獻(xiàn)：

[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òng)與換熱的高精度算法研究[D]. 上海： 華東理工大學(xué)， 2015.

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

[4] 王文全， 李偉忠， 閆妍，等. 一種投影浸入邊界法的快速直接求解方法[J]. 計(jì)算力學(xué)學(xué)報(bào)， 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] 李宇航. 仿飛蛇滑翔平板空氣動(dòng)力學(xué)性能的數(shù)值研究[D]. 北京： 中國(guó)科學(xué)院大學(xué)，2018.

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

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

[9] 豆桂芳， 吳振遠(yuǎn)， 杜艷林. 三維泊松方程的七點(diǎn)差分格式[J].工程地球物理學(xué)報(bào)，2009，6（6）：802-805.

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

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

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

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

（責(zé)任編輯：于慧梅）

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