本發(fā)明屬于離散信號計(jì)算領(lǐng)域，尤其涉及一種基于四階緊致差分格式的不可壓縮流動模擬方法。

背景技術(shù)：

在完全交錯網(wǎng)格中基于四階緊致差分格式的不可壓縮流動模擬方法，一方面通過采用四階緊致差分格式提高數(shù)值方法的精度，從而提高數(shù)值模擬的準(zhǔn)確性；另一方面，在完全交錯網(wǎng)格中的投影方法能夠精確保證求解結(jié)果滿足質(zhì)量守恒方程，并且不會出現(xiàn)非交錯網(wǎng)格中的壓力振蕩現(xiàn)象。然而，基于四階緊致格式的不可壓縮流動模擬方法中，壓力泊松方程的求解是很大的難點(diǎn)，這是本發(fā)明要解決的關(guān)鍵問題。此方法屬于通用計(jì)算流體力學(xué)方法(cfd)，能夠應(yīng)用于各種不可壓縮流動問題的仿真模擬，比如管道流動、槽道流動、機(jī)翼繞流、建筑物繞流等。

綜上所述，現(xiàn)有技術(shù)存在的技術(shù)問題是：傳統(tǒng)不可壓縮流動求解(包括各種商用cfd軟件)一般采用的二階格式方法，對于計(jì)算氣動力來講精度是足夠的；但如今工業(yè)界要求對由湍流引起的氣動噪聲、壓力脈動等現(xiàn)象有更高的預(yù)測精度，而二階格式本身精度較低、并有較大數(shù)值耗散，在應(yīng)用中有很大局限性；本發(fā)明采用四階緊致差分格式，能夠顯著提高計(jì)算精度和降低誤差，非常適用于湍流的數(shù)值模擬以及在氣動噪聲、壓力脈動預(yù)測等領(lǐng)域的應(yīng)用。

技術(shù)實(shí)現(xiàn)要素：

為解決現(xiàn)有技術(shù)存在的問題，本發(fā)明的目的在于提供一種基于四階緊致差分格式的不可壓縮流動模擬方法。

本發(fā)明是這樣實(shí)現(xiàn)的，一種完全交錯網(wǎng)格中基于四階緊致差分格式的不可壓縮流動模擬方法，該完全交錯網(wǎng)格中基于四階緊致差分格式的不可壓縮流動模擬方法包括：

在時(shí)間方向，通過微分算子的近似分裂實(shí)現(xiàn)半隱式crank-nicolson格式離散粘性項(xiàng)，對流項(xiàng)離散則采用顯式adams-bashforth格式；

對于精確投影步所產(chǎn)生的壓力poisson方程，將壓力poisson方程變換至fourier空間并得到四階緊致差分算子的修正波數(shù)，通過代數(shù)運(yùn)算獲得fourier空間中的壓力解后，再經(jīng)過fourier反變換得到物理空間中精確滿足速度散度為零的壓力場。

進(jìn)一步，adams-bashforth格式離散對流項(xiàng)和渦粘系數(shù)項(xiàng)，得到的離散化方程為

其中

進(jìn)一步，通過微分算子的近似分裂實(shí)現(xiàn)半隱式crank-nicolson格式離散粘性項(xiàng)的方法,包括：

通過block-lu分解將與δp解耦后得到的方程為

或

進(jìn)一步，待求解的x方向動量方程為

進(jìn)一步，待求解的y方向動量方程為

進(jìn)一步，待求解的z方向動量方程為

進(jìn)一步，采用四階緊致格式離散壓力梯度項(xiàng)，則有

即

其中

是n3×n3矩陣，矩陣；通過求解上述線性方程組得到壓力梯度項(xiàng)的四階離散結(jié)果。

進(jìn)一步，所述壓力poisson方程為：

dgδp＝q

其中，右端項(xiàng)

本發(fā)明提供的基于四階緊致差分格式的不可壓縮流動模擬方法，與現(xiàn)有四階精度的求解方法相比，本發(fā)明通過采用fft變換求解壓力泊松方程從而大大提高了計(jì)算效率、而采用緊致差分格式使得邊界條件的處理更為簡單。對二維taylor-oseen渦算例的計(jì)算結(jié)果表明，本發(fā)明所提出的方法具有四階空間精度，即隨著網(wǎng)格尺寸的減小，離散誤差以網(wǎng)格尺寸的四次方減小。對reτ＝180的槽道湍流dns模擬結(jié)果與譜方法dns結(jié)果吻合很好，最大誤差不超過5％。

附圖說明

圖1是本發(fā)明實(shí)施例提供的完全交錯網(wǎng)格中基于四階緊致差分格式的不可壓縮流動模擬方法流程圖。

具體實(shí)施方式

為了使本發(fā)明的目的、技術(shù)方案及優(yōu)點(diǎn)更加清楚明白，以下結(jié)合實(shí)施例，對本發(fā)明進(jìn)行進(jìn)一步詳細(xì)說明。應(yīng)當(dāng)理解，此處所描述的具體實(shí)施例僅僅用以解釋本發(fā)明，并不用于限定本發(fā)明。

下面結(jié)合附圖對本發(fā)明的應(yīng)用原理作詳細(xì)描述。

如圖1所示，一種完全交錯網(wǎng)格中基于四階緊致差分格式的不可壓縮流動模擬方法，包括：

s101：在時(shí)間方向，通過微分算子的近似分裂實(shí)現(xiàn)半隱式crank-nicolson格式離散粘性項(xiàng)，對流項(xiàng)離散則采用顯式adams-bashforth格式；

s102：對于精確投影步所產(chǎn)生的壓力poisson方程，將壓力poisson方程變換至fourier空間并得到四階緊致差分算子的修正波數(shù)，通過代數(shù)運(yùn)算獲得fourier空間中的壓力解后，再經(jīng)過fourier反變換得到物理空間中精確滿足速度散度為零的壓力場。

進(jìn)一步，adams-bashforth格式離散對流項(xiàng)和渦粘系數(shù)項(xiàng)，得到的離散化方程為

其中

進(jìn)一步，通過微分算子的近似分裂實(shí)現(xiàn)半隱式crank-nicolson格式離散粘性項(xiàng)的方法,包括：

通過block-lu分解將與δp解耦后得到的方程為

或

進(jìn)一步，待求解的x方向動量方程為

進(jìn)一步，待求解的y方向動量方程為

進(jìn)一步，待求解的z方向動量方程為

進(jìn)一步，采用四階緊致格式離散壓力梯度項(xiàng)，則有

即

其中

是n3×n3矩陣，矩陣；通過求解上述線性方程組得到壓力梯度項(xiàng)的四階離散結(jié)果。

進(jìn)一步，所述壓力poisson方程為：

dgδp＝q

其中，右端項(xiàng)

下面結(jié)合實(shí)施例對本發(fā)明的應(yīng)用原理作進(jìn)一步描述。

1、控制方程為三維非定常不可壓縮navier-stokes方程

以及連續(xù)方程

在三維直角坐標(biāo)系中，navier-stokes方程和連續(xù)方程的完整形式為：

2、本發(fā)明變量定義

空間離散采用交錯網(wǎng)格。三個方向的速度以及壓力定義為u(n1,n2-1,n3)、υ(n1,n2,n3)、w(n1,n2-1,n3)、p(n1-1,n2-1,n3-1)。x、y和z三個方向的網(wǎng)格間距分別為δx、δy和δz，其中x和z方向采用均勻網(wǎng)格，而y方向采用非均勻網(wǎng)格。

定義hj為垂向相鄰網(wǎng)格單元中心的間距，即

定義與流向和展向速度垂向二階導(dǎo)數(shù)離散相關(guān)的系數(shù)

以及(akselvoll&moin，1995)

與垂向速度沿垂向二階導(dǎo)數(shù)離散相關(guān)的系數(shù)為

定義空間離散算子h、l、g、d。其中，h為非線性對流項(xiàng)的離散算子：

l是粘性項(xiàng)的離散拉普拉斯算子：

g是梯度算子，

d是散度算子：

3、本發(fā)明半隱式離散(四階緊致差分)

在時(shí)間方向采用adams-bashforth格式離散對流項(xiàng)、crank-nicolson格式離散粘性項(xiàng)，adams-bashforth格式離散渦粘系數(shù)項(xiàng)，得到的離散化方程為

其中

與全隱式離散類似，通過block-lu分解將與δp解耦后得到的方程為

或

進(jìn)一步展開可得

p^n+1/2＝p^n-1/2+δp,

并且

則滿足的方程為

最終得到

進(jìn)一步，將l分解為l＝l1+l2+l3，分別代表對x、y和z的導(dǎo)數(shù)離散項(xiàng)，上述方程可得到如下具有二階時(shí)間精度的近似

本發(fā)明待求解的x方向動量方程為

分解之后為：

其中，

(1)左端項(xiàng)離散結(jié)果

首先離散方程

其中的二階中心差分離散結(jié)果為：

因此

j＝1：

j＝n2-1：

其次離散方程

采用二階導(dǎo)數(shù)的四階緊致格式

即

其中，

和是n1×n1矩陣。在方程兩端乘以可得到

而是循環(huán)三對角矩陣，可直接求解從而得到u^***。

最后離散方程

同樣采用四階緊致格式離散

即

其中，

和是n3×n3矩陣。在方程兩端乘以可得到

也是循環(huán)三對角矩陣，可直接求解得到u^*。

(2)右端項(xiàng)離散結(jié)果

首先對

進(jìn)行離散。

首先采用二階中心差分格式離散通過線性插值可得到

則的離散結(jié)果為：

對于首先插值得到離散點(diǎn)上的y方向速度v，即

再采用中心差分得到

對于通過四階緊致插值格式得到

即

并且

是n1×n1矩陣，是n1×n1矩陣。通過求解上面的線性方程組得到

然后采用四階緊致格式離散

即

其中

f1^x和p1^x是n1×n1矩陣。求解上述線性方程組，得到的值。

對于采用四階緊致差分離散

采用四階緊致格式離散

即

其中

f1^x和p1^x是n1×n1矩陣。求解上述線性方程組，得到的值。

最后離散首先通過四階緊致插值格式得到

即

并且

是n3×n3矩陣，是n1×n1矩陣，是n3×n3矩陣，是n1×n1矩陣。通過求解上面的線性方程組得到

而后采用四階緊致格式離散

即

其中

f1^z是n3×n3矩陣，p1^z是n3×n3矩陣。求解上述線性方程組，得到的值。

對于首先利用四階緊致格式插值得到離散點(diǎn)處的w，即

以及

離散由渦粘系數(shù)產(chǎn)生的導(dǎo)數(shù)項(xiàng)

首先離散

得出離散點(diǎn)a和b上的采用四階緊致格式

即

其中

f1^x和p1^x是n1×n1矩陣。求解上述線性方程組，得到

同樣的網(wǎng)格離散離散點(diǎn)上的和νt已知，采用四階緊致格式離散的一階導(dǎo)數(shù)：

即

其中

f1^x和p1^x是n1×n1矩陣。求解上述線性方程組，得到

離散

首先插值得到a和b點(diǎn)的νt，得到

j＝1:

j＝n2m:

同樣的網(wǎng)格，離散

j＝n2m

j＝1

插值出網(wǎng)格點(diǎn)a、b上的νt和的值，再進(jìn)一步離散可得

離散

首先對νt進(jìn)行兩次四階插值，得出a,b兩點(diǎn)的

即

并且

是n1×n1矩陣，是n3×n3矩陣，是n1×n1矩陣，是n3×n3陣，通過求解上面的線性方程組得到

同樣的四階格式離散a、b點(diǎn)

即

其中

f1^z和p1^z是n3×n3矩陣。求解上述線性方程組，得到

最終的四階差分離散公式為：

即

其中

f1^z和p1^z是n3×n3矩陣。求解上述線性方程組，得

離散

首先插值a、b點(diǎn)的νt和

即

和是n1×n1矩陣。求解上述線性方程組，得到在a、b點(diǎn)的值。

渦粘系數(shù)插值公式為

則最終離散結(jié)果為：

最后離散

插值出a、b點(diǎn)上的νt和四階緊致格式離散

四階緊致格式插值νt

則最終離散結(jié)果為

采用四階緊致格式離散壓力梯度項(xiàng)，即

即

其中

f1^m是n1×n1矩陣，p1^mn1×n1矩陣。求解上述線性方程組得到壓力梯度項(xiàng)的四階離散結(jié)果。

4、壓力poisson方程

壓力poisson方程為

dgδp＝q

其中，右端項(xiàng)

其四階緊致差分離散結(jié)果為

p′和q沿x和z方向的離散傅里葉變換結(jié)果為

以及逆變換

首先采用四階緊致格式離散x方向的壓力梯度，即

將用離散傅里葉展開，得到

則在傅里葉空間中，對x方向壓力梯度的四階緊致離散結(jié)果可表示為

再采用四階緊致格式離散x方向壓力梯度的散度，即

則在傅里葉空間中，x方向壓力梯度的散度的四階緊致離散結(jié)果可表示為

其中，

類似的，z方向壓力梯度的散度的四階緊致離散結(jié)果可表示為

其中，

將上述各項(xiàng)傅里葉變換展開式代入壓力poisson方程，得到

其中，對的二階中心差分離散結(jié)果為

j＝1(neumann條件)：

j＝n2-1(neumann條件)：

則傅里葉空間內(nèi)的壓力possion方程可以離散為

其線性方程組的系數(shù)矩陣為三對角矩陣，采用thomas追趕法求解。求得后，便可通過離散傅里葉逆變換得到p′i,j,k。

5、本發(fā)明提供的三對角線性方程組求解方法

(1)三對角方程

對于系數(shù)矩陣為三對角矩陣的線性方程組，比如

采用thomas追趕法(tdma算法)求解，本質(zhì)上仍然是gauss消去法。算法由forwardelimination和backwardsubstitution兩步組成。

在forwardelimination步，通過gauss消元將系數(shù)矩陣變?yōu)橐粋€上三角矩陣，具體步驟為：

在backwardsubstitution步，求解得到

(2)循環(huán)三對角方程

在周期邊界條件下，會出現(xiàn)系數(shù)矩陣為循環(huán)三對角矩陣的線性方程組，即

這里給出兩種循環(huán)三對角方程組的求解方法。

a.解法一

首先將循環(huán)三對角矩陣做lu分解，得到

通過等式左右兩端對比可知

β1＝b1,γ1＝c1/β1,

βi＝bi-aiγi-1,γi＝ci/βi,i＝2,...,n-1

p1＝cn,q1＝a1/β1

pi＝-pi-1γi-1,qi＝-aiqi-1/βi，i＝2,...,n-2

pn-1＝an-pn-2γn-2,

pn＝bn-(p1q1+p2q2+...+pn-1qn-1).

forwardsubstitution:

lz＝d

z1＝d1/β1,

zi＝(di-aizi-1)/βi,i＝2,...,n-1

backwardsubstitution:

ux＝z,

xn＝zn,

xn-1＝zn-1-qn-1xn,

xi＝zi-γixi+1-qixn,i＝1,...,n-2.

b.解法二

根據(jù)sherman-morrison公式，令

則有

a＝a′+uυ^t.

采用thomas追趕法求解下列三對角方程組

a′y＝d,a′q＝u,

則原方程組的解為

6、保證流量守恒的處理方法

假設(shè)壓力沿流向的變化可以表示為

p＝p′+x·px,

px是流向平均壓力梯度。p′滿足與p相同的possion方程(邊界條件不同)，即

以及

δp＝δp′+x·δpx,

gp＝gp′+px.

流量守恒條件可以表示為

∫vρδu·dv＝0,

表明不同時(shí)刻的體積流量相同。

利用

δu＝δu^*-δtgδp,

代入到流量守恒條件中，可得

從而

得到后，n+1時(shí)刻的速度更新為

壓力更新為

保持展向流量守恒的方法與流向相同。

以上所述僅為本發(fā)明的較佳實(shí)施例而已，并不用以限制本發(fā)明，凡在本發(fā)明的精神和原則之內(nèi)所作的任何修改、等同替換和改進(jìn)等，均應(yīng)包含在本發(fā)明的保護(hù)范圍之內(nèi)。