1. 前言
PHOENICS程序是世界著名的計算流體與計算傳熱學(CFD/NHT)軟件，它是英國皇家學會D.B.SPALDING教授及40多位博士20多年心血的典范之作。PHOENICS已廣泛應用于航空航天、船舶、汽車、暖通空調、環境、能源動力、化工等各個領域。在核電方面，利用PHOENICS不僅可節約大量經費，更為核電的安全可靠運行提供了可靠保證。我院于1995年10月從英國引進PHOENICS程序，經過這幾年的摸索，對其已有了大概的了解，也解決了一些小問題。但對PHOENICS這樣一個通用程序來講，它的功能還遠遠沒有得到充分的發揮，為了對PHOENICS程序進行較為深入地應用開發，讓其發揮應有的作用，本課題從理論和應用兩個角度出發，了解其內在結構，開發其各項功能，使其盡早應用于工程實際當中.
2. PHOENICS概述
PHOENICS是Parbolic,Hyperbolic or Ellicpic Numerical Integration Code Series的縮寫。它可以用來模擬流體流動、傳熱、化學反應及相關現象。程序有前處理、求解器、后處理模塊構成.
PHOENICS程序語言是標準ANSI FORTRAAN77語言，與機器無關，程序總共大約110,000條語句，2000個子程序。
PHOENICS發展歷史:
PHOENICS-81 1981年
PHOENICS-1.4 1987年
PHOENICS-1.5 1989年
PHOENICS-1.6 1991年
PHOENICS-1.6.6 1992年
PHOENICS-2.0 1993年
PHOENICS-2.1 1994年
PHOENICS-2.2 1996年
PHOENICS-3.0 1997年
PHOENICS-3.1 1998年
3. 理論基礎
3.1 控制方程(數學模型)
PHOENICS可以求解的各類問題包括:穩態、瞬態、拋物型、橢圓型、0維、1維、2維、3維方程，可接受前處理部分用戶定義的網格、材料性質、初始條件、邊界條件等。它們的求解方法基本相同，為了說明問題，本文僅以直角坐標下二維不可壓紊流運動為例說明PHOENICS程序所求解的控制方程組及其計算方法。
3.1.1 數學方程
等粘度流體的不可壓平均N-S方程組為:
r(uj )=- + (u( r ) (3-1)
=0 (3-2)
其中r 為雷諾應力項，在方程中它是未知項，它有自己的表達式，稱為湍流模型，對湍流現象的理解不同，就有不同的湍流模型，湍流模型的表達式與平均N-S方程組形成了封閉的方程組，本文采用常用的k-e湍流模型。常用的湍流模型都是建立在渦粘性概念的基礎上，雷諾應力與渦粘性的關系為:
r =mt( (3-3)
代入(3-1)式可得等粘度的不可壓平均N-S方程
r(uj )=- + ((m+mt)( ) (3-4)
=0
其中 mt=rcek2/e (3-5)
ce實驗常數 在大雷諾數流動情況下，k,e分別由下面?；耐膭幽躃方程和湍動能耗散率e方程確定。
K: em/sk )+em(( -e (3-6)
e: em/se )+c1(e/k)em(( -c2e 2/k (3-7)
其中 ce、c1、c2、sk、se 均為常數。 定常的K、e方程為:
K: uj em/sk )+em(( -e (3-6)
e: uj em/se )+c1(e/k)em(( -c2e 2/k (3-7) 按二維形式展開的控制方程組在直角坐標系內可表示為: C: r =0 (3-8)
Mx:r =- (3-9)
My:r =- (3-10)
K: r = (3-11)
e: r = (3-12)
[G]= (3-13)
根據具體問題給出具體的邊界條件和初始條件，即組成了完整的控制方程組及其定解條件。
3.1.2 離散方程 將控制方程(3-8)---(3-13)用統一的輸運方程表示為:
r = +S (3-20)
式中f為個方程中的因變量，如在Mx方程中為u，而在連續方程中f為1等，S表示除式(3-20)左邊的對流項及右邊的擴散項外的所有項之和，稱為輸運方程的源項，擴散項在各方程中也不同，如動量方程中G=mt+m，e方程中G=mt/se等。 各方程中源項分別為:
Mx: S=- G=m+mt
My: S=- G=m+mt
c : S=0 G=0
K : S=mt[G]-re G=mt/sK
e : S= G=mt/se (3-21)
令Jx=ruf-G (3-22)
將(3-22)代入(3-20)，則(3-20)簡化為
(3-23) 至此，我們將控制方程(3-8)至(3-13)轉換成統一的輸運方程(3-23)形式。 PHOENICS 采用控制容積法進行方程離散。所謂控制容積法即在有限小的（由網格確定)控制容積內，對方程(3-22_作一次積分，使方程降階后，再按一定的差分格式離散方程。圖3-1為X-Y平面上的網格與控制容積的關系。
Y N

n
W w P e E
s
S
? X 圖中P表示中心節點，N、S、W、E為該節點周圍最近的四個節點。虛線內為該節點的控制容積，n、s、w、e為控制容積各個面上的中點(對于均勻網格)，但若網格非均勻，上述交點可不為中點。在圖3-1所示的控制容積內對方程3-23進行積分得:
Je-Jw+Jn-Js= DxDy (3-24)
Je、Jw、Jn和Js是整個控制容積面上的積分總流量，即Je個界面e上的 依此類推，S為 (3-25)
類似地，我們在整個控制容積內積分連續方程可得:
Fe-Fw+Fn-Fs=0 (3-26)
式中Fe、Fw、Fn及Fs是通過控制容積面的流量，如果在點e的ru代表整個界面e 上的值，我們就可以導出:
Fe=(ru)eDy
類似地: Fw=(ru)wDy
Fn=(ru)nDy
Fs=(ru)sDy
(3-24)-F*(3-26)得:
(Je-FeF)-(Jw-FwF)+(Jn-FnF)-(Js-FsF)= DxDy (3-28)
為說明問題，本文我們采用混合格式離散方程即:
Je-FeFp=aE(Fp-FE)
Jw-FwFp=aW(Fw-Fp)
Jn-FnFp=aN(Fp-FN)
Js-FsFp=aS(Fs-Fp)
式中:
aE=DeA(|pe||)+[[-Fe,0]]
aW=DwA(|pw||)+[[-Fw,0]]
aN=DnA(|pn||)+[[-Fn,0]]
aS=DsA(|ps||)+[[-Fs,0]] De= Dw=
Dn= Ds=
貝克列數定義為:
Pe=Fe/De Pw=Fw/Dw Pn=Fs/Ds
函數A(|P|)采用混合方安推薦的公式:
A(|P|)=[[0,1-0.5|p|]]
符號[[x,y]]=max(x,y)
將源項S盡可能地線性化為 S=SC+SPFP
至此，我們可以把二維的離散化方程寫成:
apFp=aeFE+awFw+aNFN+asFs+b
b=ScDyDx
式中:
ap=aE+aw+aN+as-spDxDy
至此，方程離散化已完成。
理論上講，到目前為止已經可以進行計算求解了，但目前的離散方法計算很容易產生鋸齒形壓力場，而這又是不合理的，一般解決該問題的方法是才用交錯網格法。所謂交錯網格即:把速度u、v及壓力p分別儲存在三套不同網格上的網格系統，u控制容積與主控制容積之間x方向有半個網格步長的錯位，而v控制容積與主控制容積之間在y方向上有半個步長的錯位。
在交錯網格中一般F變量的離散過程及結果與3.1.2 節所述相同。但對動量方程而言，則帶來一些新的特點:
a.積分用的控制容積不是主控容積而是u、v各自的控制容積。
b. 壓力梯度項從源項中分離出來。例如對ue的控制容積:
?(pp-pe)Dy
這里假設在ue的控制容積的東、西界面上壓力是各自均勻的，分別為pE、pp。于是關于ue的離散方程具有以下形式:
aeue=?anbunb+b+(pp-pe)Ae
類似地，對vn的控制容積作積分可得:
anvn=?anbvnb+b+(pp-pN)An 3.1.3 計算方法 3.1.3.1 SIMPLE算法的計算步驟 采用SIMPLE算法實施關于u、v、p代數方程的分離式求解時，計算步驟如下: (1) 假定一個速度分布，記為u0，v0，以次計算動量離散方程的系數及常數項;
(2) 假定一個壓力場p*;
(3) 依次求解兩個動量方程，得u*、v* ;
(4) 求解壓力修正值方程，得p’ ;
(5) 據p’改進速度值 ;
(6) 利用改進后的速度場求解那些通過源項物性等與速度場耦合的F變量。如果F并不影響流場，則應在速度場收斂后再求解 ;
(7) 利用改進后的速度場重新計算動量離散方程的系數，并用改進后的壓力場作為下一層次迭代計算的初值。重復上述步驟，直到獲得收斂的解。 PHOENICS程序計算方法采用的是SIMPLEST算法，與SIMPLE相比，它主要有以下兩個特點: (1) 對流項采用迎風格式，因為這是一個絕對穩定的格式，且擴散項與對流項的影響系數可以分離開來，不象指數(或乘方)格式那樣綜合在一起，至于由迎風差分所引起的假擴散問題，則采用逐步加密網格、以獲得與網格稀密程度無關的解這種做法加以克服。 (2) 把相鄰點的影響系數表示成對流分量cnb及擴散分量dnb之和，并把對流部分全部歸入源項，于是ue的動量方程為:
aeue=?( dnb+cnb)unb+b+Ae(pP-pE)
=?dnbunb+(b+?cnbu*nb )+Ae(pP-pE)
由此可見，當擴散項略而不計時，動量方程實際上采用Jacobi的點迭代。點迭代的收斂速度是比較慢的，但是由于對流項與壓力之間的耦合關系等原因，正希望利用這一特性以防止迭代發散。這種混合式的計算方法有利于促進強烈非線性問題的迭代過程收斂，SIMPLEST的計算步驟與SIMPLE基本相同。 PHOENICS求解時可采用點迭代、線迭代、面迭代等方法迭代求解。
圖3-2 為PHOENICS采用全場求解方法時的計算步驟: DO ISITEP = 1, LSPTEP
DO ISWEEP = 1,LSWEEP
DO IZ = 1,NZ
Apply previous sweep’s pressure &
velocity corrections DO IC = 1,LITC
Solve scalars in order
KE, EP, H1, C1, C2, .... C35
ENDDO
Solve velocities in order
V1, U1, W1
Construct and store pressure correction sources
and coefficients
ENDDO
Solve and store pressure corrections whole field
ENDDO
ENDDO
圖 3-2
4. 結論 本文僅對PHOENICS程序的控制方程及計算方法進行了簡單的介紹，以給出PHOENICS程序求解問題的大概方法。PHOENICS程序是一個大型通用計算程序，可計算的領域很多，視個人理論基礎的不同，計算的結果和應用的范圍差別很大。而要想完全掌握PHOENICS程序的理論部分，最好在掌握了PHOENICS一般理論基礎上，結合課題逐步對其理論及方法研究掌握。PHOENICS程序附有完整的幫助系統，在使用當中遇到問題時可隨時查閱。這些仍不能滿足要求時可查閱PHOENICS雜志及報告。