二维椭圆型流动传热通用程序.docx
《二维椭圆型流动传热通用程序.docx》由会员分享，可在线阅读，更多相关《二维椭圆型流动传热通用程序.docx（20页珍藏版）》请在冰豆网上搜索。
二维椭圆型流动传热通用程序
二维椭圆型流动传热通用程序
变量表及算例说明
（本材料仅供教学参考）
西安交通大学
CFD&NHT/EHT研究中心
陶文铨教授
2002/10/15
一、FORTRAN变量表3
二、关于程序的主要说明6
（1）二维椭圆型流动和传热问题通用计算机算法方面的特点6
（2）各程序的主要功能7
（3）三种坐标系统8
（4）网格系统与节点命名方法9
十一个例题的已知条件与求解内容12
例题1直角坐标中二维稳态无内热源的导热12
例题2空心圆柱内的稳态热传导12
例题3正方形管道内的充分发展对流换热13
例题4内壁上有直肋的环形通道内充分发展对流换热13
例题5给定流场条件下温度场的计算14
例题6二维突扩通道中的流动与换热14
例题7方形通道内的复杂充分流动15
例题8旋转圆盘上的冲击流动15
例题9轴对称燃烧内的瞬间燃烧过程16
例题10带中心射流的通道内的紊流换热16
例题11有质量源的流动问题17
FORTRAN变量表
ListofFORTRANVariables
ACOF
QuantitycalculatedbysubroutineDIFLOWtogivethecombinedconvectionanddiffusioneffect.
AIM（I,J）
ThecoefficientOn.
AIP（I,J）
ThecoefficientaE.
AJM（I,J）
ThecoefficientOs.
AJP（I,J）
ThecoefficientOn.
AP（I,J）t
Thecoefficientap;alsoSpinGAMSOR.
APT
Theunsteadyterm"/匚t.
AREA
Localvariable,usuallytheareaofaC.V.face.
ARHO
a
Localvariable,（area）（）.
ARX（J）
TheareaofthemainC.V.facenormaltothexdirection.
ARXJ（J）
ThepartofARX（J）thatoverlapsontheC.V.forV（I,J）.
ARXJP（J）
ThepartofARX（J）thatoverlapsontheC.V.forV（I,J+1）.
BL
BLC
Coefficientsusedintheblockcorrection.
BLM
BLP
CON（I,J）
Theconstanttermbinthediscriminationequation;alsostandsforSCinGAMSOR.
DENOM
Temporarystorage.
DIFF
DiffusionconductanceD.
DT
Thetimestept.
DU（I,J）
deInfluencingU（I,J）.
DV（I,J）
dnInfluencingV（I,J）.
F（I,J,NF）
Various.
FL
TemporarystorageleadingtoFLOW.
FLM
TemporarystorageleadingtoFLOW.
FLOW
MassflowratethroughaC.V.face.
FLP
TemporarystorageleadingtoFLOW.
FV（J）
Interpolationfactorswhichgivethemassflow
FVP（J）
VrAtamaingridpoint,I,JasFV（J）V「（I,J）+FVP（I,J）
：
?
vr（i,j+1）
FX（I）
Interpolationfactors,whichgivetheinterface.
FXM（I）
DensityRHOM（atthelocationofU（I,J））asFX（I）■-RHO（I,J）+FXM（I）
RHO（I-1,J）.
FY（J）
Interpolationfactors,whichgivetheinterface.
FYM（J）
DensityRHOM（atthelocationofV（I,J））asFY（J）RHO（I,J）+FYM（J）
RHO（I,J-1）.
GAM（I,J）
Thediffusioncoefficient-.
I
Indexdenotingthepositioninx.
IBEG
TemporaryvaluesusedinPRINT.
IEND
IFST
ThefirstvalueofIforwhichtheprint-outisarranged;usedinPRINT.
II
Temporaryindex.
IPREF
ThevalueofIforthegridpoint,whichisusedasareferenceforpressure.
1ST
Thefirstinternal-pointvalueofI.
ISTF
IST-1;usedinSOLVE.
ITER
Acounterforiterations.
IT1
TemporaryvaluesusedinSOLVE.
IT2
J
Indexdenotingthepositioniny.
JFL
TemporaryindexusedinPRINT.
JFST
SimilartoILST.
JJ
Temporaryindex
JLST
SimilartoILST
JPREF
SimilartoIPREF.
JST
Thefirstinternal-pointvalueofJ.
JSTF
JST-1;usedinSOLVE.
JT1
TemporaryvaluesusedinSOLVE.
JT2
LAST
Themaximumnumberofiterationsallowedbytheuser.
LBLK（NF）
When.TRUE.TheblockcorrectionforF（I,J,NF）isused.
LPRINT（NF）
When.TRUE.,F（I,J,NF）isprinted.
LSOLVE（NF）
When.TRUE.,WesolveforF（I,J,NF）.
LSTOP
When.TRUE.Computationstops.
L1
ThevalueofIforthelastgridlocationinthexdirection.
L2
（L1-1）.
L3
（L1-2）.
MODE
Indexforthecoordinatesystem;=1forXy,=2for「X,=3for「二.
M1
ThevalueofIforthelastgridlocationintheydirection.
M2
（M1-1）.
M3
（M1-2）.
N
TemporarystorageforNF.
NF
6
Indexdenotingaparticular.
NFMAX
ThelargestvalueofNFforwhichstorageisassigned.
NGAM
NFMAX+3;GAM（I,J）canbeconsideredasF（I,J,NGAM）.
NP
NFMAX+1;P（I,J）canbeconsideredasF（I,J,NP）.
NRHO
NFMAX+2;RHO（I,J）canbeconsideredasF（I,J,NRHO）.
NTIMES（NF）
ThenumberofrepetitionsofthesweepsinSOLVEforthevariableF（I,J,NF）.
P（I,J）
Thepressurep.
PC（I,J）
1
Thepressurecorrectionp.
PREF
Thepressureatthereferencepoint.
PT（I）orPT（J）
QT（I）orQT（J）
TransformedcoefficientsintheTDMA.
R（J）
TheradiusrforamaingridpointI,J.
REL
1.0-RELAX（NF）.
RELAX（NF）
RelaxationfactorforF（I,J,NF）.
RHO（I,J）
p
Thedensity.
RHOCON
Thevalueofforaconstant-densityproblem.
RMN（J）
ThevalueofradiusrforthelocationtowhichV（I,J）refers.
SMAX
I
Thelargestabsolutevalueofthe“masssource”PusequaAiohne
1
SSUM
1
Thealgebraicsumofallthe“masssourcespiietfisition.
SX（J）
ScalefactorforthexdirectionatthemaingridlocationsY（J）.
SXMN（J）
ScalefactorforthexdirectionatinterfacelocationsYV（J）.
TEMP
Temporarystorage.
TIME
Timetforunsteadyproblems.
TITLE（NF）
AlphabetictitleforF（I,J,NF）.
U（I,J）
Thex-directionvelocityu.
V（I,J）
They-directionvelocityv.
VOL
VolumeoftheC.V.
X（I）
Thevaluesofxatgridpoints.
XCV（I）
Thex-directionwidthsofmainC.V..
XCVI（I）
ThepartofXCV（I）thatoverlapsontheC.V.forU（I,J）.
XCVIP（I）
ThepartofXCV（I）thatoverlapsontheC.V.forU（I+1,J）
XCVS（I）
Thex-directionwidthofthestaggeredC.V.forU（I,J）.
XDIF（I）
ThedifferenceX（I）-X（I-1）.
XL
Thex-directionlengthofthecalculationdomain.
XU（I）
ThelocationsoftheC.V.faces;i.e.thelocationofU（I,J）.
丫（J）
Thevaluesofyatgridpoints.
YCV（J）
They-directionwidthsofmainC.V.
YCVR（J）
Thearear-yforamainC.V.
YCVRS（J）
Thearear:
yfortheC.V.forv（I,J）.
YCVS（J）
They-directionwidthofthestaggeredC.V.forV（I,J）.
YDIF（J）
ThedifferenceY（J）-Y（J-1）.
YL
They-directionlengthofthecalculationdomain.
YV（J）
ThelocationsoftheC.V.faces;i.e.thelocationofV（I,J）.
二关于程序的主要说明
（1）二维椭圆型流动和传热问题通用计算机程序算法方面的特点
1、采用原始变量法，即以速度U、V及压力P作为直接求解的变量
2、守恒型的差分格式，离散方程系对守恒型的控制方程通过对控制容积作积分而得出
的，无论网格疏密程度如何，均满足在计算区域内守恒的条件；
3、采用区域离散化方法B，即先定控制体界面、再定节点位置
4、采用交叉网格，速度U、V与其他变量分别存储于三套网格系统中；
5、不同的项在空间离散化过程中去不同的型线假设，源项采用局部线性化方法；扩散一
—对流项采用乘方格式（但很容易转化为中心差分、迎风差分或混合格式）；街面上
的扩散系数采用调和平均法，而密度与流速则用线性插值；
6、不稳态问题采用全隐格式，以保证在任何时间步长下均可获得具有物理意义的解；
7、边界条件采用附加源项法处理；
8、耦合的流速与压力采用SIMPLE算法来求解；
9、迭代式的求解方法，对非线性问题，整个求解过程具有迭代性质；对于代数方程也采
用迭代法求解；
10、采用交替方向先迭代法求解代数方程并补以块修正技术以促进收敛。
（2）各方程块的主要功能
1、主程序MAIN:
规定整个计算过程的流程，决定是否停止计算。
2、子程序SETUP1:
设置与网格系统有关的，在计算过程中保持不变的几何参数及设置
U、V、P'P、p、CON（Sc）、AP（Sp）的初值。
3、子程序SETUP2:
1.建立每一轮迭代中变量U、V及P'和其他变量离散方程的系数as、aw、a”、as、ap及b项。
所有变量公用一套数组逐个求解，逐个确定各自目标的系数，
以节省内存。
在建立系数过程中调用DIFLOW子程序；2调用SOLVE子程序求解代
数方程；3迭代次数指标及时间步长增值；
4、子程序DIFLOW：
计算乘方格式中的系数A（|p|）;
5、子程序SOLVE:
采用交替方向线迭代法并辅以块修正技术求解代数方程；
6、子程序GRID：
设置为建立网格系统所必须的量，包括
（1）设置x,y方向求解区域的宽度XL及YL
（2）设置x,y方向上各自的节点数L1及M1
（3）设置x,y方向上控制体积的界面位置，XU（I）,I=2,L1;YV（J）,J=2,M1。
如为均分网格可确定子程序UGRID,如为非均分网格需由用户在GRID中把XU（I）,YV（J）一一设置好。
（4）规定坐标系，MODE=1，2，3分别相应于直角、圆柱轴对称及极坐标
（5）当MODE不等于1时，规定径向起始点半径R（I）
7、子程序START:
设置初值，包括
（1）、对不稳态问题给出初始条件
（2）、对稳态问题给出迭代求解的假定值
（3）、已知的边界值也可在此块中，在设置初值或假定值时一并送入，但对边界条件随时间而异的不稳态不能在此块中赋值，因在整个计算过程中此会仅执行一次。
8、子程序DENSE:
规定密度场，对于密度为常数的问题，此块中可不设任何语句，但必须保留ENTRY及RETUIRN语句
9、子程序BOUND:
规定边界条件，包括
（1）设置各变量的边界条件
（2）对各特殊问题所需的量，如Nu,fRe等等可在此块中计算，但也可在OUTPUT中设置
10、子程序OUTPUT：
打印输出，包括
（1）每做一轮迭代，输出一行信息，以观察收敛情形；
（2）调用PRINT子程序，实现二维物理量场的输出；
（3）用户所需输出的其他特殊变量
11、子程序GAMSOR：
设置扩散系数及源项，包括
（1）设置控制方程中规定的源项Sc,Sp
（2）对第二，三类边界条件设置与边界相邻的控制体中的附加源项，Sc,ad及Sp,ad
（3）内节点及边界节点规定扩散系数，扩散采用附加源项法时令边界扩散系数为零。
12、子程序UGRID:
为均分网格设置界面位置
13、子程序PRINT:
输出二维物理场
（3）二种坐标系统直角坐标系，Z方向厚度为1;M0DE=1。
x,l
轴对称圆柱坐标，计算对0=1弧度的区域进行；M0DE=2;Y（J）可以从任何起点算起，但R（J）须从对称轴开始；R
（1）维计算区的最小半径。
两个区域XL,YL相等，但R⑴不同
极坐标系，Z方向厚度为1;M0DE=3;计算区域的0角须小于2n;Y（J）可从任何点起算，但R（J）比从中心点开始；R
（1）为计算区域边界的最小半径。
MODE
X（I）
Y（J）
R（J）
SX（J）
1
x
y
1.0
1.0
2
x
y
r
1.0
3
0
y
r
r
⑷网格系统与节点命名方法
4
V（2,M1）
U（L1,M2）
：
U（2,M2）]
Fi**iT
■1•I•丨・p（L1i,M1）*“
F*
j
L!
iIL!
XCV（I）
（c）x方向的主要几何参数
ARXJP（J）
M1y（M1）yv（M1）
XCVS（3）
M2
U（2,J）
XDIF（3）
（d）速度U（3,J）的控制容
ARX（J）
ARXJ（J）
YDIF（J）ycvs（j）
ycv（J-1）
yv（3）
2y
（2）
1y（i）yv
（2）
V（I,J）|I
11»»—
I,J
]
3y（3）
（e）y方向的主要几何参
边界压力P（1,2）、P（2,1）由内点值境外
推而得；
参考压力点p（1,1）按下式计算：
P（1,1）=P（1,2）+P（2,1）-P（2,2）
变量内点下标起始值
变量
1ST
JST
①,P',P
2
2
U
3
2
V
2
3
三、^一个例题的已知条件与求解内容
例题1直角坐标中二维稳态无内热源的导热
已知：
求解区域如图1所示，四个边界上的温度由下式决定；T=x+y+xy
求：
该区域中的温度分布。
图1
例题2空心圆柱内的稳态热传导
已知：
求解区域如图2所示。
左边界为给定温度，右边界为对流，上边界为绝热，下边界为给定热流。
整个计算区域内都有源项。
求：
该区域中的温度分布。
计算方法：
用附加源项法处理第2,3类边界条件
例题3正方形管道内的充分发展对流换热
已知：
常物性流体在均匀壁温的正方形截面管道内作充分发展的对流换热（如图3）。
求：
截面上的速度分布、温度分布、fRe及Nu。
计算方法：
充分发展对流换热问题的处理特点
例题4内壁上有直肋的环形通道内充分发展对流换热
已知：
常物性流体在带直肋的环形通道内作充分发展层流对流换热，内表面为均匀壁温（周向）,轴向则呈线性变化，外表面绝热。
Ri=1,R2=2,a=15o。
求：
截面上速度分布、温度分布及fRe、Nu。
计算方法：
耦合问题的一种处理方法。
例题5给定流场条件下温度场的计算
已知：
流体流过一个直角的两表面，速度场为u=Ax,V=-AY,a=10,Tin=500,Tw=100
求计算区域中温度的分布。
Tin
Tout
绝热
例题6二维突扩通道中的流动与换热
已知：
由两平行平板组成的突扩通道，尺寸如图所示。
流动为层流。
入口流速均匀，Vin=100,
入口温度均匀，Tin=100,通道壁温均匀，Tw=300。
流体Pr=0.7,分子粘性尸1，密度按下式变化，
密度分布及压力场。
计算方法：
出口边界条件的处理
图6
例题7方形通道内的复杂充分流动
已知：
常物性流体在如下图示方形管道内作充分发展的流动与换热（重力项中的密度采用
Boussines假设）.p3=10：
T仁0,T2=1,PR=0.7,尸1.0dp/dz=const计算中取-3000,进入充分发展时右边进入的热量等于左边导出的热量•
求截面上速度U、V的分布，Z方向的分量W的分布，温度分布及压力场。
计算方法：
如何用一个二维程序来计算三个速度分量。
例题9轴对称燃烧内的瞬间燃烧过程
已知：
燃料与空气如图9所示进入轴对称燃烧室，设：
1、仅有三种组分：
燃料，氧气及燃烧产物
2、氧气与燃料的扩散系数相等
3、反应在瞬间内完成
4、对各个组分Cp均相同
U1=2O,U2=5O,尸1,Pr=0.7燃料焓Hfu=3X1O【Vo=O.2U2求：
燃料室中的速度场、温度场与浓度场。
例题10带中心射流的通道内的紊流换热
已知：
如图10，带一般中心射流的流体进入一平行板通道。
流动为紊流，卩=10「6
T=100
V=10
例题11
已知：
一深为0.2米的湖泊尺寸如图11所示。
湖泊中有一源与汇。
湖面上与空气发生换热。
Tair=100,h=2。
Uin=100,Tin=500,湖的四个边均视为绝热。
源处的入口强度为4。
湖深方向当
作均匀处理（即按X——Y坐标中的二维问题处理）。
求：
湖泊中的速度、温度分布
计算方法：
质量源项方法的应用。
图11
升级会员 