SIMPLE算法求解方腔內(nèi)粘性不可壓流動學習資料

1、精品文檔SIMPLE算法求解方腔內(nèi)粘性不可壓流動目錄一、問題描述 2二、離散格式 3交錯網(wǎng)格 3方程離散 4三、SIMPLE算法基本思想 7邊界條件處理 8虛擬網(wǎng)格處理 9方程求解 11輸出變量處理 12SIMPLE算法流程圖15四、程序中主要變量的意義 15五、計算結(jié)果與討論 17函數(shù)最大值 17變量等值線圖 18主要結(jié)論 22六、源程序22一、問題描述假設0 Wx, y W1的方腔內(nèi)充滿粘性不可壓縮流體，左、右、下壁固定，上壁以u = 16x2(1 -x2 )運動，試求Re =100, 200, 400時的定常解，方腔如圖 1所示。II = -1圖1方腔內(nèi)流動示意圖精品文檔:、離散格式本算

2、例采用求解不可壓縮流動的經(jīng)典算法，即SIMPLE算法，求解方腔內(nèi)粘性不可壓縮流體運動的定常解。SIMPLE算法的全稱為 Semi-Implicit Method for Pressure-Linked Equations , 即求解壓力關(guān)聯(lián)方程的半隱式算法。采用SIMPLE算法時，為了避免中心差分格式將“棋盤”型參量分布誤認為是均勻分布， 需要用交錯網(wǎng)格對計算域進行離散。交錯網(wǎng)格交錯網(wǎng)格如圖2所示，壓力、密度等物理量存儲在控制體（i,j ）的中心，這個控制體稱 為主控制體。速度分量 （u,v）分別存儲在主控制體的 （i +1/2, j ）和（i, j +1/2）位置處,標記為（i, j ）位置

3、，再分別以此為中心，劃分速度分量u、v的控制體。采用空間均勻網(wǎng)格,等間距離散整個求解域，如圖 3所示。圖2交錯網(wǎng)格示意圖圖3求解域離散示意圖圖3中陰影部分代表方腔內(nèi)的流動區(qū)域，陰影區(qū)域的邊界代表方腔的上、下、左、右壁面，陰影區(qū)域外面的網(wǎng)格節(jié)點是為邊界處理需要而設定的虛擬網(wǎng)格節(jié)點，后面介紹邊界處理方法時詳細論述。方程離散無量綱化的守恒型不可壓縮 N -S方程為，U = 02U=0-U 一 . UU P-t其積分形式為y UdS = 0u一一 1 ，一dVn U udSpnxdS - 一 ： n、udS = 0V 4ss r x Re s：:v1.dV - n UvdS ( pnydS - 一 ；

4、 n、vdS = 0V ： t SSRe S圖4主控制體圖5速度u控制體圖6 速度v控制體采用有限體積法離散 N -S方程，連續(xù)性方程在主控制體上離散M 1 M 1M 1 M 1Ui,j -Uid,j y Vi,j -M,j/：x=0X方向動量方程在速度u控制體上離散，時間采用前差第(UiMj*Y 汁(Fj-F 沏+(Gj-G電 3x + (PiM；-PiMj* 的=0 LtY方向動量方程在速度v控制體上離散，時間采用前差第錯書 W)+(F啕 3x + (Gj-G電的+ (pM-pMj、x=0 t其中，數(shù)值通量F1=u2-工,，Guv-4Re 二xRe 二 yF 2 =v2工 g2 =uv？R

5、e xRe：y通量fC )gC)分別定義在主控制體的中心和角點，如圖所示，并按照如下方法離散11 M M M 1 M 11 M 1 M 1FUi 1,jUi,j U 1,j Ui,j - 5 u -Ui,j4Re x. 11m MM1 M1 1 M 1 M 1GV1,j - vi,j U1,j . Ui,jui1,j -ui,j4Re y通量F2,G2分別定義在主控制體的中心和角點，如圖所示，并按照如下方式離散F 2 =-VM1 j ,j4M M 1 M 1Vi,jVi 1,jVi,j一G2 4uM1,j4M M 1 M 1 ui,jVi 1, j Vi,j1IRe y1Re xM 1 M 1

6、Vi 1,j -Vi,jM 1 M 1Vi U Vi,j通量F QG加F（2） G?用某些項凍結(jié)于 M時間層，使離散化之后的方程對uM書，vM斗是線性的。將離散化之后的fC）g”口 F（2）,G（2%弋入離首W的x方向和y方向的動量方程,整理之后得離散后的動量方程如下u M 1 a u i,J ui,Jv M 1 a V ai,j vi, jT . u- a p,quVaP,qVM 1 p,qM 1 p,q.uM 1bi,j - Pi 1,j.vM 1bi,j . Pi -M 1 - pi,jM 1 一 pi,j其中bujxyM ui,juai 1,juai,j 1二y gu，x 4 viMM

7、i 1,j , ui,jM M i,j - Vi d,j-u1.M, ,Ma=y 4 ui 1,j1 M-y 二 ui,j.4M 1ui-1 j -i1,jRexu, ai,jA 1,1 / M , M 皿 1-Ax - (Vi j/+ Vi+jrM ,j ,jRe”)+ 2 +/xJ(VM+VM _VM _VM )+ 2)ReixJ* i Vi+，j Vi，j,Vf Re”一bivjM vi,jv1MMai,j1 X - vi,j- vi,j 1IL4v . I 1 M Mai 1,j = y - ui,j .1 , ui,j,4.二Rey 1,v ai,jvai,jv A 1,1 / M

8、M ,2州崎,j”奇以上是SIMPLE算法中離散化的動量方程.1 M M-x - j-Vi,j,41 M=r y - Ui,j 1,4ReAy_M i，jMUi,jRe x工Re x SIMPLE算法基本思想SIMPLE算法是一種解決壓力-速度耦合問題的“半隱式”算法。首先給定 M時刻猜測的速度場UM ,VM ,用于計算離散動量方程中的系數(shù)和常數(shù)項。給定 M+1時刻猜測的壓力場 *估計值p ,迭代求解離散動量方程，得到M+1時刻速度場的估計值 U ,V ,速度場的估計*值u ,v滿足如下離散方程。U *U * U*a u a u b + ：p4- p y = 0ai,jui,j- ap,qup

9、,qi,jpi 1,jpi,jyv *v *v*a v a v b t p p x = 0ai,jvi,j- ap,qvp,q i,j pi,j 1pi,jx *.*因而需要對速度場u , v和壓力場p*般地，速度場u ,v不滿足離散的連續(xù)性方程, 進行修正。M+1時刻的修正值和估計值有如下關(guān)系M ：；1u = uM -1v u vM ：1p其中，u ； v和p分別速度和壓力的修正量，修正量亦滿足離散的動量方程au u au ubu0ai,jui,j-ap,qup,qi,jpi 1,jpi,jyav v x av v - bvx - 0ai, j vi,j一ap,qvp,qi,jpi,j 1p

10、i,jx編號為（i,j ）的速度修正量u：v不僅與壓力修正量p有關(guān)，還與鄰近點的速度修正量有關(guān)。SIMPLE算法的重要假定：速度的改變只與壓力的改變有關(guān)，忽略鄰近點對速度修 正的影響。因而得到如下速度修正量a=T r 1,j - Rjai,j : x , .vi,j = - pi,j 1 - pi,jai,j修正后的速度分量M 1*uUui,jui,jy .M 1*VVvi,jVi,jpi 1,j - pi,jai,jvT pi,j 1 - pi,jai,j將修正后的速度分量代入離散后的連續(xù)性方程，得到壓力修正方程pppai,j pi,j =ap,q pp,qbi ,j其中2.2.2.2p y

11、 p y p x p xai 1，j = -, ai,j =_u- , ai,j 1 =, ai,j .1 ;F-ai,jai,jai,jai, j J12222p y y x x ai,j r ai,jai,jai,jai,jp*bi,j - - yUi,j -Ui-xVi,j -Vi,jj采用迭代法求解壓力修正方程，得到壓力修正量p,代入修正公式得到M+1時刻的速一 M 1 M HM 1M 1 M ; 1M 1度場U , V和壓力場p 。將M+1時刻的速度場u , V和壓力場p作為新的猜測的速度場和猜測的壓力場估計值，采用上述方法計算下一個時刻的速度場和壓力場，直到滿足收斂條件。收斂判據(jù)M

12、ax（bipi ） 名1, j名為很小的正實數(shù)，視計算的精度要求而定。本算例中取 名=1e-8。若b/=0,則=0,此時U41 =U ,從而來自于離散動量的U*,j滿足離散的連續(xù)性方程。因此MaxMpj）8可以作為收斂判據(jù)。邊界條件處理首先對計算區(qū)域離散，并流動邊界之外擴充一個虛擬網(wǎng)格，將真實流動的離散域包圍, 如圖7所示。7虛擬網(wǎng)格擴充示意圖壓力p的存儲位置由圖中符號,代表，速度分量u的存儲位置由圖中符號代表,速度分量v的存儲位置由圖中符號 代表。陰影區(qū)域表示真實的流動區(qū)域，陰影區(qū)域外的網(wǎng)格為虛擬網(wǎng)格。速度分量u在x方向虛擬網(wǎng)格上的存儲位置分別記為uj , und1.j ，在y方向虛擬網(wǎng)格上

13、的存儲位置分別記為Ui,0 , Ui,n書。速度分量v在x方向虛擬網(wǎng)格上的存儲位置分別記為V0,j, Vn書.j ，在y方向虛擬網(wǎng)格上的存儲位置分別記為Vi,，Ui.n+。壓力和壓力修正量在x方向虛擬網(wǎng)格上的存儲位置分別記為po,j, pn丸j, p0,j, p；+,j ,在y方向虛擬網(wǎng)格上的存儲位置分另1J記為pi,o, pi,n+, p：0, p：,n書。虛擬網(wǎng)格處理靜止壁面的虛擬網(wǎng)格按圖 8所示方式處理，可以滿足粘性邊界條件，以下邊界為例說 明。圖8 靜止壁面虛擬網(wǎng)格處理示意圖邊界外的虛擬網(wǎng)格存儲位置上的速度分量，由邊界內(nèi)的真實網(wǎng)格對應存儲位置上的速 度分量沿邊界對稱得到。這樣可以保證虛

14、擬網(wǎng)格與真實網(wǎng)格在對應存儲位置上的速度分量大 小相等，方向相反，二者在邊界上的線性平均值為零，滿足粘性邊界條件。上壁面有速度，為運動壁面，速度分量u #0, v = 0。速度分量v采用前述對稱方法處理。速度分量u的處理方法略微不同。上壁面速度Uup,oundary不隨時間變化，虛擬網(wǎng)格上速度分量ui,n + 由 u up -boundary , ui,n 工線性插值得到，可以保證上壁面滿足邊界條件，即ui,n 書2 父 Uup boundary 一 ui,n 工。存儲位置位于恰好位于靜止壁面上的速度分量恒為零，不必迭代計算，例如u0,j = 0。邊界外虛擬網(wǎng)格中心的壓力和壓力修正量，認為分別等

15、于邊界附近對應真實網(wǎng)格中心的壓力和壓力修正量，即假定區(qū)域邊界附近沿邊界法向的壓力梯度為零。因此，只需求解真實網(wǎng)格上各個存儲位置上變量的值，虛擬網(wǎng)格上各個存儲位置上變量的值可以通過前述處理虛擬網(wǎng)格的方式得到。至此，包括虛擬網(wǎng)格在內(nèi)的所有存儲位置上都被賦予了確定的值。方程求解X方向動量方程 * 迭代法求解x萬向動量萬程，若要得到M+1時刻Uij的值，需要用到 M時刻M M M M MM M M MUi, j，Ui書,j，Uij，Ui ,j書，Ui,j，vi, j ,vi書,j , vi ,j，vi,書j九1位置的速度分重值和 M+1時刻 * * . 一 . . *Pi書,pi兩個位置的壓力值。內(nèi)點

16、上的值，可以直接迭代得到。 求M+1時刻速度分量u1,1 , 需要用到M時刻u1MIUuMcuNuMo，v1M1,v2MVMoNMo九個位置的速度分量值和M+1時 刻P2,1,P1,1兩個位置的壓力值，包括虛擬網(wǎng)格在內(nèi)的離散域滿足計算需求。同理可得， 在其他邊界位置，離散域滿足求解X方向動量方程的需要。Y方向動量方程迭代法求解y萬向動量萬程，若要得到M+1時刻vi i的值，需要用到 M時刻i , jvM vM . vM . vM * vM， UM uM . UMUM書九個位置的速度分量值和M+1時刻1，1書，12，1書，1 j,j，J|H1,j，Ji,j，j,d1ji八 f g n* * .

17、一 . . *Pi書，Pi兩個位置的壓力值。內(nèi)點上的值，可以直接迭代得到。求M+1時刻速度分量v1,1,需要用到M時刻v1,1, v2,1, v0,1，v1,2，v1,0 , U1,1 ,U2,1,U1,0,U2,0九個位置的速度分量值和M + 1時*刻Pl,2, P1,1兩個位置的壓力值，包括虛擬網(wǎng)格在內(nèi)的離散域滿足計算需求。同理可得，在 其他邊界位置，離散域滿足求解y方向動量方程的需要。壓力修正方程迭代法求解壓力修正方程，若要M+1時刻得到P。的值，需要用到 M時刻Pi1,j, Pi1 , P； j * P：, j 4四個位置的壓力修正值和M+1時刻*ui,j , ui 書j, ui,j,

18、 ui,j 書, ui,j，ui,j 書, ui,j, ui-2,j ，*vi,j , vi,j 書, vi,j4, vi書,j , vi 4,j , vi書,j，vi,j, vi,j-2十六個位置的速度分量估計值。內(nèi)點上的值，可以直接迭代得到。求M+1時刻壓力修正值P；,需要用到M時刻P2,1, P；,1, P；2, P；o四個位置的壓力修正值和M+1時刻*U1,1 , U2,1 , U0,1, U1,2, U1,0, U0,2, U0,0, U_1,1*V1,1 , V1,2, V1,0, V2,1, V0,1, V2,0, V0,0, V1, _1十六個位置的速度分量估計值，包括虛擬網(wǎng)格

19、在內(nèi)的離散域滿足計算需求。同理可得, 在其他邊界位置，離散域滿足求解壓力修正方程的要求。將M+1時刻變量的估計值和修正值代入如下公式*U*VM ：；1UM -1V即可得到M+1時刻變量U M中，V M的值。輸由變量處理渦量定義；:u渦量是一個矢量，在平面流動中只有一個方向，垂直于流動平面。流函數(shù)滿足u =：x：:y流函數(shù)等值線為流線c*dx ：:ydy - - vdx udy在交錯網(wǎng)格里，對渦量和流函數(shù)在一個網(wǎng)格內(nèi)做中心差分，渦量和流函數(shù)的存儲位置為真實網(wǎng)格的節(jié)點，如圖 9所示。圖9渦量和流函數(shù)存儲位置示意圖渦量差分格式vi 1, j - vi, j ui 1, j - ui, jco =-x

20、yy方向流函數(shù)差分格式= yui, j yui, j i, j -ix方向流函數(shù)差分格式二- xvi, jxvi, j i -1, j任取其中一個式子計算流函數(shù)即可，或者將兩個式子的計算結(jié)果取平均值。本算例采用的是取平均值的辦法。邊界上的流函數(shù)值為零。需要輸出用于后處理的變量主要有位置坐標x、v,速度分量u、v,壓力p,流函數(shù)e和渦量O。由于SIMPLE算法采用交錯網(wǎng)格，不同性質(zhì)的變量存儲在網(wǎng)格的不同位置，在變量輸出時，有必要將所有變量統(tǒng)一到真實網(wǎng)格的網(wǎng)格節(jié)點上，便于后處理，如圖10所示。圖10變量統(tǒng)一后存儲位置示意圖,，一 一一1網(wǎng)格節(jié)點上的速度分量 u?為節(jié)點上下兩速度分量的平均值，即u?

21、=(ui i + ui i書)2 i, ji, j 1 1網(wǎng)格節(jié)點上的速度分量 。為節(jié)點左右兩速度分量的平均值，即？=1(vi,j +VT,j )網(wǎng)格節(jié)點上的壓力 ？為節(jié)點周圍四個主控制體中心壓力的平均值，即C 1?pi,jPi,j I . pi 1,j - Pi 1,j -14邊界處的壓力平均采用相同的公式，需要用到虛擬網(wǎng)格中心的壓力，角點依然。SIMPLE算法流程圖開始統(tǒng)一變量存儲位置，后處理結(jié)束圖11 SIMPLE算法流程圖四、程序中主要變量的意義主程序中主要變量的意義變量意義counter_max預設最大循環(huán)次數(shù)grid_number、nX方向和y方向的離散網(wǎng)格數(shù)目t、h時間步長和空間

22、步長Re雷諾數(shù)1d收斂判據(jù)，壓力修正方程中bi*絕對值的最大值i , j數(shù)組u、vX方向和y方向速度分量數(shù)組 P、p_wave壓力和壓力修正量數(shù)組 flow_function、 vortex_function流函數(shù)和渦量函數(shù)數(shù)組 u_scalar速度絕對值數(shù)組 d_array、g、u_0計算中需要用到的一些輔助數(shù)組子程序1,即求解x方向動量方程程序中主要變量的意義變量意義nX方向和y方向的離散網(wǎng)格數(shù)目t、h時間步長和空間步長Re雷諾數(shù)a、b、c、d、e迭代求解x方向動量方程所需的輔助變量 n數(shù)組u、v、pX方向、y方向速度分里和壓力場子程序2,即求解y方向動量方程程序中主要變量的意義變量意義n

23、X方向和y方向的離散網(wǎng)格數(shù)目t、h時間步長和空間步長Re雷諾數(shù)a、b、c、d、e迭代求解x方向動量方程所需的輔助變量 數(shù)組u、v、pX方向、y方向速度分里和壓力場子程序3,即求解壓力修正方程程序中主要變量的意義變量意義nX方向和y方向的離散網(wǎng)格數(shù)目t、h時間步長和空間步長Re雷諾數(shù)a、b、c、e、f壓力修正方程中各壓力修正項系數(shù)1d壓力修止方程中k項a_u1、 a_u2輔助變量，與u方向速度修正式有關(guān)a_v1、 a_v2輔助變量，與v方向速度修正式有關(guān)error求解壓力修正方程的迭代收斂判據(jù)數(shù)組u、vX方向、y方向速度分量數(shù)組 p、p_wave壓力和壓力修正量五、計算結(jié)果與討論本算例分別對網(wǎng)格

24、數(shù)目為 20 X 20,40X 40兩種情況，數(shù)值求解了 Re=100,200,400時的方 腔流動的定常解，計算結(jié)果列舉如下。函數(shù)最大值離散網(wǎng)格為20X 20數(shù)值計算得到的流函數(shù)最大值Re100200400流函數(shù)最大值0.29269532920.27242897280.2427700376對應x坐標0.400.450.45對應y坐標0.650.600.60離散網(wǎng)格為20X 20數(shù)值計算得到的上壁面渦函數(shù)最大值Re100200400上壁向渦函數(shù)最大值91.7220275878104.8169276505117.6126516741對應x坐標0.800.750.75對應y坐標1.001.001.

25、00離散網(wǎng)格為20X 20數(shù)值計算得到的方腔中線渦函數(shù)最大值Re100200400中線上渦函數(shù)最大值25.400771350538.822847263853.8979317133對應x坐標0.500.500.50對應y坐標0.951.001.00離散網(wǎng)格為40X 40數(shù)值計算得到的流函數(shù)最大值Re100200400流函數(shù)最大值10.31935425740.31347035100.2972193667對應x坐標0.4250.4500.475對應y坐標0.6250.6000.575離散網(wǎng)格為40X 40數(shù)值計算得到的上壁面渦函數(shù)最大值Re100200400上壁向渦函數(shù)最大值119.06996768

26、98149.6299870265181.4581848045對應x坐標0.8250.8000.800對應y坐標1.001.001.00離散網(wǎng)格為40X 40數(shù)值計算得到的方腔中線渦函數(shù)最大值Re100200400中線上渦函數(shù)最大值32.418443541041.919028264845.2529732845對應x坐標0.500.500.50對應y坐標0.9750.9750.975從表中數(shù)據(jù)可以看到，隨著雷諾數(shù)增大，流函數(shù)最大值減小，最大流函數(shù)的位置也更趨 于方腔的中心。這是因為雷諾數(shù)增大，粘性效應減弱，上壁面運動帶動放腔內(nèi)流體運動的能力降低。上壁面渦量最大值位于上壁面的速度最大處，這是因為在該

27、處速度分量 u沿y方向的的梯度最大，速度分量v沿x方向的梯度為零（壁面粘性邊界條件），因而剪切效應最強。 隨著雷諾數(shù)增大，上壁面渦量最大值增大。中線上渦量最大值出現(xiàn)在上壁面附近，但不位于上壁面。離散網(wǎng)格為 20X 20,雷諾數(shù)為100時，數(shù)值計算結(jié)果表明，中線上最大渦量出現(xiàn) 在y=0.95處；雷諾數(shù)為200和400時，數(shù)值計算結(jié)果表明，中線上最大渦量出現(xiàn)在y=1.00處。離散網(wǎng)格為40X40,對于雷諾數(shù)為100、200和400三種情況，數(shù)值計算結(jié)果表明，中 線上最大渦量出現(xiàn)在 y=0.975處。這說明增大離散網(wǎng)格數(shù)目，即加密網(wǎng)格，可以增加數(shù)值計 算對渦量最大值的捕捉精度。變量等值線圖速度大小等

28、值線圖12 Re=100時速度大小等值線/ 6 -2 o O 0 d精品文檔圖15 Re=100時流函數(shù)等值線精品文檔圖13 Re=200時速度大小等值線圖14 Re=400時速度大小等值線流函數(shù)等值線0.14 0.12 0.1 0 08 MR AM 期旗0 QQH 。0061 Q川口迎 -0 AC&4 -Q.0-125。也。Q 口 -DCI.FLOW-F .03 .023 02S-0 K022Q2 一。博 口他 0 14 0.12I打 0.09 0.00口嶷I 0.02 0.0077J 口 EU1工必包 .00125精品文檔圖18 Re=100時渦量等值線精品文檔J3.00S4 0.0135

29、Fl的X2Q第 M22Z 拘幡1412 QO.口 000.O.。0 I O.Ofi 0.06 。國 0.02 0.0077 0.OT51 QRes200FLOWFfl QMS 力川儂 -0.01250.22 力;D.18 (116 0.14 0.12 fl.1 d.dQ CIjOB 亞3 flfla 卿77 0.00&1圖16 Re=200時流函數(shù)等值線ReMOOFLW-F值 07 0等 016 D14 0 12 01 D0 QiD6QD4 D.QEDDOTT 000&1 0-OC046-0G125FLCW FR-4fiC61GMJ5 .2仲仲141213.302即期曲的012 QO.QQ 0

30、口0,口0.0,?？?。通力力圖17 Re=400時流函數(shù)等值線渦量等值線精品文檔精品文檔4030201000冷|山泊口町口 。心力。 1 1 1 1 n -K- ? 6 5 4 3 J 1 D 1 .c - - - -HMHHM-1圖19 Re=200時渦量等值線圖20 Re=400時渦量等值線從以上各個變量的等值線圖中可明顯看出，相同雷諾數(shù)不同網(wǎng)格密度時， 各個變量相應的等值線形狀基本相同， 但40 X 40網(wǎng)格的等值線明顯比 20 X 20網(wǎng)格的等值線光滑。這說明 加密網(wǎng)格可以提高數(shù)值計算的精度。從流函數(shù)等值線中可以看到，Re=100和200時，方腔左下角處流函數(shù)均勻，Re=400時,方腔

31、左下角處流函數(shù)不均勻，出現(xiàn)了與方腔中心流函數(shù)等值線不一致的流函數(shù)等值線。這說明在增大雷諾數(shù)，流體粘性減弱的情況下，方腔下壁面角點處可能出現(xiàn)回流。因此本算例對Re=1600的情形在40X40網(wǎng)格上做了數(shù)值計算，結(jié)果如圖所示。圖中流函數(shù)等值線分布表 明，在方腔左下角附近，確實出現(xiàn)了回流，但是回流區(qū)的流函數(shù)值很小，說明回流強度很小。40X40FLOW-F0.0077O.OQ61 0-0.0046-0.0084-0.0125o.o.o.o.so.so.o.o.圖21 Re=1600時流函數(shù)等值線主要結(jié)論加密離散網(wǎng)格，可以提高數(shù)值計算的精度。隨著雷諾數(shù)的增加， 流體粘性效應減弱，上壁面驅(qū)動流體的能力減弱

32、，放腔內(nèi)流函數(shù)整體減小，下壁面角點附近出現(xiàn)回流，上壁面渦量增大。壁面上的渦量可以為正， 也可以為負，渦量的最大值出現(xiàn)在速度分量u沿y方向梯度與速度分量v沿x方向梯度只差最大的地方，因而不一定出現(xiàn)在壁面上。六、源程序program mainimplicit none integer , parameter : counter_max=50000 !設置最大的時間推進步數(shù)為五萬步real(8) , parameter : standard=0.00000001!小數(shù)點后取八位integer : i,j,i_max,j_max,n,counter,grid_number!grid_numer 為網(wǎng)格

33、數(shù)目real(8) , allocatableu(:,:),v(:,:),p(:,:),p wave(:,:),d array(:,:),g(:),u 0(:,:),flow function(:,:),flow function_u(:,:),flow_function_v(:,:),vortex_function(:,:),u_scalar(:,:)real(8) : h,t,Re !h 為空間步長，t 為時間步長，Rej Reynolds number real(8) : d,p_wave_maxwrite (*,*) 輸入離散網(wǎng)格數(shù)目： read (*,*)grid_numberwri

34、te (*,*) 輸入 Reynolds number ： read (*,*)Recounter=1 ! 迭代步數(shù)計數(shù)器h=1.0/ real (grid_number)n=grid_numbert=0.001 allocate (u(-1:n+1,0:n+1),v(0:n+1,-1:n+1),p(0:n+1,0:n+1),p_wave(0:n+1,0:n+1),d_array(n,n),g(0:n),u_0(-1:n+1,0:n+1),flow_function(0:n,0:n),vortex_function(0:n,0:n),u_scalar(0:n,0:n),flow_functio

35、n_u(0:n,0:n),flow_function_v(0:n,0:n)!聲明二維數(shù)組的長度，p_wave代表修正壓力u=0! 賦初值v=0do i=1,n-1u(i,n+1)=-16.0*(real (i)*h)*2*(1.0-(real (i)*h)*2)*2-u(i,n)end dop=0p_wave=0d=1.0do while(d=standard)if (countercounter_max) exit! 如果時間推進步數(shù)大于五萬步，跳出循環(huán)u0=ucallsub1_x_momentum(u,v,p,h,t,n,Re)! 求解離散的 x 方向動量方程callsub2_y_mome

36、ntum(u_0,v,p,h,t,n,Re)! 求解離散的y 方向動量方程callsub3_pressure(u,v,p,p_wave,n,t,Re)! 求解壓力修正方程do i=1,ndo j=1,nd_array(i,j)= 10.0*abs(u(i,j)-u(i-1,j)+v(i,j)-v(i,j-1)end doend dod=maxval(d_array)!d 作為收斂判據(jù)do i=1,ndo j=1,nd_array(i,j)=abs(p_wave(i,j) ! 只考慮實際的壓力波動，即邊界內(nèi)的壓力波動end doend do p_wave_max=maxval (d_array)

37、counter=counter+1write (*,*)d=,d,p_wave_max=,p_wave_max,counter-1 end dowrite (*,*)counter-1flow_function_u = 0! 求流函數(shù)flow_function_v = 0do i=1,n-1do j=1,nflow_function_u(i,j)= h*u(i,j)+flow_function_u(i,j-1)end doend dodo j=1,n-1do i=1,nflow_function_v(i,j)= -1.0*h*v(i,j)+flow_function_v(i-1,j)end d

38、oend dodo i=0,ndo j=0,nflow_function(i,j)= 0.5*(flow_function_u(i,j)+flow_function_v(i,j) end doend dodo i=1,n-1! 求渦量函數(shù)do j=1,n-1vortex_function(i,j)= (v(i+1,j)-v(i,j)/h - (u(i,j+1)-u(i,j)/h end doend do! 求邊界上的渦量函數(shù)do j=1,n-1!left wallvortex_function(0,j)= (v(1,j)-v(0,j)/h - (u(0,j+1)-u(0,j)/h end do

39、do j=1,n-1vortex_function(n,j)= (v(n+1,j)-v(n,j)/h - (u(n,j+1)-u(n,j)/h end dodo i=1,n-1vortex_function(i,n)= (v(i+1,n)-v(i,n)/h - (u(i,n+1)-u(i,n)/h end dodo i=1,n-1vortex_function(i,0)= (v(i+1,0)-v(i,0)/h - (u(i,1)-u(i,0)/h end do!right wall!up wall!down wall! 計算角點渦量函數(shù)vortex_function(0,0)= 0.5*(vo

40、rtex_function(0,1) + vortex_function(1,0)vortex_function(0,n)= 0.5*(vortex_function(0,n-1) + vortex_function(1,n)vortex_function(n,0)= 0.5*(vortex_function(n-1,0) + vortex_function(n,1)vortex_function(n,n)= 0.5*(vortex_function(n,n-1) + vortex_function(n-1,n)do j=0,n!將u速度節(jié)點與網(wǎng)格節(jié)點重合do i=0,ng(i)=(u(i,j

41、+1)+u(i,j)*0.5end dodo i=0,nu(i,j)=g(i)end doend dodo i=0,n!將v速度節(jié)點與網(wǎng)格節(jié)點重合do j=0,ng(j)=(v(i+1,j)+v(i,j)*0.5end dodo j=0,nv(i,j)=g(j)end doend dodo j=0,n! 求節(jié)點速度絕對值do i=0,nu_scalar(i,j)= sqrt (u(i,j)*u(i,j) + v(i,j)*v(i,j)end doend dodo i=0,n! 將壓力節(jié)點與網(wǎng)格節(jié)點重合do j=0,np_wave(i,j)=0.25*(p(i,j)+p(i+1,j)+p(i,j

42、+1)+p(i+1,j+1)end do end dop=p_waveopen(unit=10,file=data.dat)write (10,(TITLE = Grid=,I3,Re=,f5.1,)n,Re! 輸出雙引號時需要打出兩個雙引號，否則出錯write (10,*)VARIABLES = X Y U V U-SCALAR P FLOW-F VORTEX-F ! 這種情況不需要打兩個雙引號就可以輸出雙引號write (10,(/)write (10,(ZONE N=,I6, E=,I6, ZONETYPE=FEQuadrilateral,DATAPACKING=POINT)(n+1)*

43、(n+1),n*ndo j=0,ndo i=0,nreal (i)*h, real (j)*h,u(iwrite (10,(f15.10,f15.10,f15.10,f15.10,f15.10,f15.10,f15.10,f15.10),j),v(i,j),u_scalar(i,j),p(i,j),flow_function(i,j),vortex_function(i,j)end doend dodo j=0,n-1! 對網(wǎng)格節(jié)點進行編號，便于tecplot 軟件后處理i= j*(n+1)counter=1do while (counter=n)write (10,(I5,I5,I5,I5)

44、i+counter, i+counter+1, i+counter+1+n+1, i+counter+n+1counter=counter+1end doend doopen(unit=11,file=function_max.txt)d=0.0do j=0,n! 流函數(shù)最大值do i=0,nif (d=flow_function(i,j) thend=flow_function(i,j)i_max=ij_max=jend ifend doend doreal (i_max)*h, real (j_max)*hwrite (11,*)flow_functioin_max: write (11,(f15.10,I5,I5,f15.10,f15.10)d,i_max,j_max,d=0.0上壁面渦函數(shù)最大值do i=0,nif (d=vortex_function(i,n)thend=vortex_function(i,n)i_max=iend ifend dowrite (11,*)up_wall_vortex_functioin_max:write (11,(f15.10,I5,f15.10)d,i_max,r