|
|
马上注册，结交更多好友，享用更多功能，让你轻松玩转社区。
您需要 登录 才可以下载或查看，没有账号？注册
x
哈工大清洁能源技术论坛
计算流体力学——PHOENICS专版 [返回]
PHOENICS程序应用-理论基础部分
标记论坛所有内容为已读
>> 计算流体力学——PHOENICS专版欢迎您的到来 <<
◆您是本帖第 24 个阅读者◆
* 贴子主题： PHOENICS程序应用-理论基础部分
zhangjie
--------------------------------------------------------------------------------
目录
1. 前言
2. PHOENICS概述
3. 理论部分
3.1.1 控制方程
3.1.2 数值方程
3.1.3 计算方法
4. 结论
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分别由下面模化的湍动能K方程和湍动能耗散率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)
Jy=rvf-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杂志及报告。
|
|
IP卡
狗仔卡
发表于 2003-5-30 15:57:39
提升卡
置顶卡
沉默卡
喧嚣卡
变色卡
显身卡
发表于 2003-5-30 17:42:03
