1
Title:
Different Methods Obtained by PHOENICS Simulation to Improve the Performance of
Pusher

Type Steel Slab Reheating Furnace
Author :
Yong Tang Jarmo Laine Timo Fabritus Jouko Härkki
Affiliation:
Process metall
urgy laboratory, PL4300, Oulu University, 90014 Finland
Fax: +358 8 553 2339
Tel: +358 8 553 2423
Email:
yong.tang@oulu.fi
Website:
http://cc
.oulu.fi/~pometwww/
Date:
6
th
,June,2002
Computer and operating system used:
Intel Pentium 1.7 GHz, 128 RAM, Windows 98 second edition
Abstract:
Phoenics3.1 was used to simulate the gas temperature and flow velocity
distribution in the
pusher

type steel slab

reheating furnace Extended Simple Chemically

Reacting System
(ESCRS)
, k

ε model and composite

flux radiation model were included in the numerical
computing. Author coded some special boundary conditions in the Phoenics GROUND file.
Performance of the reheating furnace in a steel industry company was analysed according to
the simulation results. In order to obtain a better heating performance of slab and reduce the
possibility of scale accumulation near the lower burner in the heating zone, a block wall was
added in front of the burners during the simulation.
Contents li
st:
Title and Abstract Page 1
Target Page 2
Model Description
Page 2
PHOENICS
Settings and Iteration Process Page 6
Results Page 7
Verificat
ion of the Model Page 9
Conclusions
Page 9
Reference and Nomenclature
Page 10
Appendix of Q1 and GROUND files
Page 12
2
Target
The p
usher

type reheating furnace is installed to heat steel plates before they enter the mill.
The uniform slab t
emperature profile is the most concern of the furnace operators.
Consequently, it is necessary to understand the temperature distribution and flow type in the
furnace. This CFD simulation include two parts: modelling the gas temperature and velocity
in the
furnace; investigate the gas flow modification while a block wall is built in front of the
lower burners in the heating zone.
Model Description
a. The Structure of the Pusher

Type Furnace
The furnace modelled is shown in Fig. 1 and its highest
roof is about 6 meters inside.
There are eight burners in the soaking zone. In the heating zone, six burners are arranged on the upper
side while seven burners are set up on the lower part as shown in the picture. Flows from upper part
burners are injected
11
o
downwards. Coke oven gas is used as the fuel.
Since it is symmetric, half of
the furnace was simulated. Three

dimensional modeling and Cartesian coordinate is
employed in this furnace CFD simulation.
Fig.1 Schematic
of the pusher ty
pe slab

reheating furnace
The calculation domain is 28.07 metres in z direction, 4.25 metres in x direction and 5.66
metres in y direction. A rectangle grid is chosen, and so the irregular geometry parts of
furnace, such as transition between heating
zone and convection zone, are simplified by
dividing them into several rectangular cuboids. The total number of cells is about
31
39
68=82212,as shown in Fig.2. In the calculation, some cells are blocked according to
the structure of the furnace. Supportin
g pillars are considered in this simulation if compared
with previous work done by Maki and Yang
[1

3]
.
8.5 m
X
Y
Z
3
Fig.2 Grids of furnace modelled
b. Mathematical models applied in this simulation
According to the models provided by the PHOENICS3.1 and the re
al situation of the
pusher

type furnace, standard
K

model is chosen for the turbulent flow,
Extended Simple
Chemically

Reacting System (ESCRS) is selected to simulate the combustion and the
radiation model is composite

flux. Flow in the furnace is consid
ered as a single phase, steady
state, and Newtonian fluid. The following is the simple mathematic descriptions of models
mentioned above.
Continuity equation:
0
)
(
u
(1)
Momentum conservation equation:
f
u
p
g
u
u
2
(2)
Energy conservation equation:
h
S
kT
H
u
)
(
(3)
Chemical species conservat
ion equation:
l
l
l
R
C
C
u
)
(
(4)
K

equation:
)
(
)
(
)
(
k
i
k
e
i
i
i
P
x
K
x
K
u
x
(5)
)
(
)
(
)
(
2
1
C
P
C
K
x
x
u
x
k
i
e
i
i
i
¶
¶
¶
¶
¶
¶
(6)
The combustion model used is ES
CRS in PHOENICS and the reaction is assumed as
irreversible, i.e. the rate of the reverse reaction is presumed to be very low
[4]
. The diffusion
coefficients of fuel, oxidant and product are presumed equal to each other, and to the
diffusivity of heat. Con
sequently, their Prandtl/Schmidt numbers are equal, which is a good
approximation for turbulent flow.
In this reheating furnace, fuel is coke oven gas with component of H
2
, CO and CH
4
. (Table
1). The oxidizer is air (21%O
2
and 79% N
2
). The chemical rea
ctions defined as
2 steps 3
reactions finite

rate,
are expressed below:
2
2
4
4H
2CO
O
2CH
(7)
2
2
CO
2
O
2CO
(8)
O
H
2
O
2H
2
2
2
(9)
According to the Eddy

Break

Up (EBU) model, this kind of hydroc
arbon combustion
reaction is assumed rather fast in the furnace and the combustion stage is controlled by the
4
mixing rate of reactant in a turbulent scale. So the reaction rate of the fuel species can be
deduced from:
k
s
M
M
Min
C
R
O
Fu
EBU
)
/
,
(
2
(10)
The Min (…) argument implies that R is dependent on the species in the shortest supply.
Significantly, the reaction rate and further enthalpy source are linked with the momentum and
continuity equations mentioned above. The mas
s fractions of fuel and oxidizer will be
obtained from those equations automatically when the fuel input data is added to the model in
PHOENICS.
Table 1. The composition of gas fuel used in the steel

reheating furnace
CO %
H
2
%
CH
4
%
N
2
%
7.4
54
28
1
0.6
The thermodynamic data of the gas, such as the specific heat and heats of formation, are
deduced by the CHEMKIN interface in PHOENICS, while the gas composition and
temperature are solved.
The composite

flux model of Spalding is chosen
[5]
fo
r radiation calculation, which has
been widely used in combustion chambers and furnaces. The difficult part of radiation
calculation is how to treat the radiation angle. There are different ways for such radiation
angle discretisation. One method is to div
ide the radiation direction from arbitrary angles into
six directions in Cartesian coordinates, as used in composite

flux model (six

flux). The
following assumptions are given for the radiation model: radiation is transmitted in coordinate
directions only;
no interlinkages, apart from scattering, arise between the radiation fluxes in
the respective coordinate directions; the scattering term presumes that the scattering is
isotropic with angle, which is probably reasonable only if the total contribution of s
cattering
is not large. Thus, the differential equations describing the variations of the six radiation
fluxes in Cartesian coordinates are:
SUM
K
I
K
I
K
K
dy
dI
S
b
a
S
a
)
(
(11)
SUM
K
I
K
J
K
K
dy
dJ
S
b
a
S
a
)
(
(12)
SUM
K
I
K
K
K
K
dz
dK
S
b
a
S
a
)
(
(13)
SUM
K
I
K
L
K
K
dz
dL
S
b
a
S
a
)
(
(14)
SUM
K
I
K
M
K
K
dx
dM
S
b
a
S
a
)
(
(15)
SUM
K
I
K
N
K
K
dx
dN
S
b
a
S
a
)
(
(16)
The term
SUM
in
equation
(11)

(16) is the average of all the fluxes in the modelled
system. For three

dimensions:
6
/
)
(
N
M
L
K
J
I
SUM
(17)
5
The composite radiation fluxes are defined as:
R
y
= (I+J)/2 ; R
z
= (K+L)/2 ; and R
x
= (M+N)/2 (18)
The contribution of
radiation
to the energy equation source terms is given by:
b
a
b
z
y
x
i
i
a
radiation
I
N
M
L
K
J
I
K
I
R
K
S
6
3
2
,
,
,
(19)
The absorption coefficient K
a
=0.1, the scattering coefficient K
s
=0.1, and the emissivity
coefficient of the walls which enters through the boundary conditions, is
w
=0.85.
C. Boundary conditions
In the simulation, the slabs are
in contact with each other and form a whole block region
along the longitudinal direction, as shown in Fig.1.
Since the main target is to understand the gas flow and temperature distribution in the
reheating furnace, the following boundary conditions s
hould be known before the calculation:
surface temperature of slabs, temperature of walls in the furnace, the mass flux (or velocity),
temperature and composition of gas flow from the burners. The outlet boundary conditions for
the furnace also need to be
specified.
The surface temperature of slabs was measured in the furnace modelled. Because the
furnace is considered steady, the temperature measured at different time can be referred to as
the surface temperature of the slabs at different sections alon
g the furnace.
The inlet boundary conditions such as mass flux, velocity, chemical composition, enthalpy
and
turbulent kinetic energy k, dissipation rate
can be obtained by certain calculations from
the specific inlet gas composition, temperature and
mass inflow, which are coded in Phoenics
Q1 file. Since the circle inlet is assumed as square, the equivalent length of each square side
can be described as follows:
2
)
2
(
D
L
square
(20)
Here D is the
diameter
of the circular inlet; L
squa
re
is the equivalent length of each square
inlet side.
The kinetic energy K and dissipation rate
are obtained from the equations given below:
2
2
3
U
T
K
Int
(21)
2
07
.
0
09
.
0
2
3
75
.
0
square
L
K
(22)
Where T
Int
is turbulence i
ntensity, its value is 2% in this furnace simulation. U is the gas
velocity from the inlet.
The composition of gas fuel has been given in Table 1, and Table 2 provides the burner
conditions for the inlet flow.
Table 2 Conditions of the burners inlet
flow
6
Position of
burners
Number
Of
burners
Burner’s
Diameter
[mm]
Coke
gas
[Nm
3
/h]
Air
[Nm
3
/h]
Flow Velocity
[m/s]
The
temperature of
air and coke
gas inflow is
350 and 50
Celsius degree
separately.
U
z
U
y
Soaking zone
8
380
900
3800
4.3
0
Uppe
r heating
zone
6
380
2200
11200
10.3
2.0
Lower heating
zone
7
300
2100
10800
8.6
0
There are two outlets (exhaust pipes) in the reheating furnace. In the modelling process the
external pressure is assumed to be atmosphere pressure (1.01
10
5
Pa), o
r the relative pressure
is 0 Pa. The discharge door area is about 1.856m
2
.The estimation relative pressure in the
furnace near the discharge door is 0.3 Bar (0.2999
10
5
Pa). The discharge door is assumed be
open in the model.
Several temperature measur
ement points are projected in the furnace, which are just near
the inside walls. The measured temperature is used as the reference inner surface temperature
of wall without further consideration of heat transfer through the wall. The walls in the
furnace w
ere divided into several parts. Each part was given an estimated boundary
temperature according to the measured value.
Since the heat capacity of the gas varies with temperature and composition, the gas
enthalpy near the wall is obtained in the Ground
file of Phoenics. This method is simply
described as follows.
The gas component in the furnace includes N
2
, H
2
, CH4, H2O and CO
2
. Their heat
capacity, for example Cp
N2
(T), Cp
H2
(T), Cp
CH4
(T), Cp
H2O
(T)
and Cp
CO2
(T)
,
is the function of
temperature T. Aft
er the composition data of gas near the wall, such as Y
N2
,Y
H2
, Y
CH4
,Y
H2O
and Y
CO2
are obtained in the iteration process, the heat capacity of gas at temperature T is:
2
2
2
2
4
4
2
2
2
2
CO
CO
O
H
O
H
CH
CH
H
H
N
N
gas
Y
Cp
Y
Cp
Y
Cp
Y
Cp
Y
Cp
Cp
(23)
The enthalpy of gas attached to the boundary wall (ga
s near the wall), H
gas,wb
can be
obtained from:
form
refer
gas
bw
gas
H
H
T
Cp
H
,
(24)
PHOENICS
Settings and Iteration Process
Phoenics Version 3.1 of MS

DOS was used in this simulation.
a. SATELLITE
Pleas
e see the appendix1 of Q1 files that satellite read. The satellite module operates in the
batch model and stop at the first STOP line. There are no other special requirements for
SATELLITE.
b. GROUND
7
The thermal boundary conditions are determined in t
he calculation and coded in the
GROUND file as shown in the appendix 2. After the GROUND file is compiled and re

link,
private executables (earexe.exe) is created. Type “ run77 earexe” to start private EARTH.
c. Iteration
More than 2000 sweeps was ap
plied in the iteration for this simulation. Variables, such as
P1,KE, EP and H1 use the linear relaxation and , velocity components, chemical species CO,
H2 and CH4, use the false step time relaxation. At the beginning of the iteration, the
relaxation coef
ficients are big and then set small values while the iteration proceeding. The
number
of F

array locations available is about 15 MB.
It took more than 12 hours for this
calculation in the IBM compatible personal computer with Intel Processor PIV 1.7 GHZ. T
he
convergence is good in the calculation. Errors of all the variables concerned except chemical
composition of gas drop three orders during the iteration. (from 10
4

6
to 10
1

3
). However, the
calculation errors of chemical composition of gas only fall 2 or
ders. When all the values of
the monitored spot are not changed significantly, convergence is archived. Several special
monitor spots were chosen in the calculation, all the results showed a satisfied convergence.
Some calculation grids were modified to
study its influence on the simulation results. No
significant result changes occur in this kind of modification.
Results
a.
Flow Pattern and Gas Temperature Distribution
Fig.3 shows the gas flow pattern and gas temperature distribution along the l
ongitudinal
direction, which crosses the burners i
n
the heating zone. Flow pattern obtained in this
numerical model is like the type of the physical model carried out by Matsunaga
[6]
. Reverse
flow near the roof close to the burners, exists either in the he
ating zone or soaking zone.
There is also reverse flow in the corner between floor and lower burners wall. However,
simulation results indicate that there is a strong reverse flow under the slab and near the lower
burners, which did not appear in the water
model. This difference owes to the floor step in the
current furnace, while bottom floor of heating zone is flat in the previous physical model
according to the report.
Fig.3 (a) shows that flames from the lower burners slant upwards when they bounce
the
step of floor in heating zone. Flames of upper burners inject downward because of the burners
installation angle. This arrangement of flames is good for the efficient heating and would not
overheat the surface of slab. The gas temperature in heating zo
ne is more than 1260
o
C except
the regions near the convection zone.
b. The flow modification when a block wall is installed in the heating zone
This CFD simulation
results
indicate situation could be improved if proper block wall is
installed in fro
nt of the lower burners in the furnace.
Fig.4 (a) is the gas flow distribution in the pusher

type reheating furnace. Fig.4 (b) shows
the flow vector of plane x=8 with one block wall in the furnace.
Flow patterns modeled
suggest
that
reverse flow streng
th is weakened when one block wall is installed near the
burners in the heating zone. Weak reverse flow is a good sign in reducing scale accumulating
near the lower burners.
8
(a)
Gas temperature distribution
(b)
Gas flow pattern in the furnace
Fig.3
—
Gas
flow pattern and temperature along longitudinal furnace
(Cross burners in the heating zone)
(a) Without block wall
(b) With a block wall
Fig.4
Gas velocity distribution near the burners in the heating zone
9
Verification
Industry verification was carried out in this furnace since there are several measuring points
installed in the furnace.
Gas temperature close to the flame cannot be meas
ured because the flame temperature
injected from the burners is very high and work conditions are serious. More attentions are
focused on the mixture gas temperature above the plane of slab, gas temperature near the roof
and sidewall. Fig.5 is the comparis
on of modeled results with measured gas temperature at
different positions in the furnace. The head of curve through grids x=17 and y=29 crosses area
near roof in soaking zone. The middle part of curve of grids through x=17 and y=39 pass the
area near the
roof in the heating zone. Tail of curve with grids x=17 and y=30 is near the roof
of the convection zone. The results show that calculated gas temperatures in the soaking zone
and fore part of the heating zone are very close to the measured values. However
, in the
transition areas from the heating zone to the convection zone, calculated values are higher
than measured gas temperature. Overestimate radiation in these transition areas maybe the
reason for this kind of deviation.
Although the calculation
error for chemical components in the furnace is a little high,
comparison between the oxygen distribution in the furnace and modelled oxygen values is
implemented. It seems that calculated oxygen near the furnace wall is not so far away from
the measured d
ata, as shown in Fig.6.
Conclusions
Computational fluid dynamics is applied to predict gas flow pattern and temperature
distribution in the pusher

type reheating furnace. Since until now available measurements are
limited for the details of the gas
temperature distribution and flow pattern, CFD simulation
plays an important role in the investigation. In the present, momentum, combustion and
radiation models are combined together. Results of the model were verified with the plant
data. Situations rela
ted to heating improvement and reducing reverse flow under the slab near
the lower burners is modeled.
Fig.5
—
Comparison of calculated gas temperature with measured results at different positions in the furnace
10
Fig.6
—
Comparison of modeled O
2
distribution with measured values in the furnace
Reference
1. Anne Mäki:
Aihioiden Kuumennusuunin Numeerinen Virtausmallinnus.
Lisensiaatintyö,Oulu
University,2001
2.Yongxiang Yang, Ari Jokilaakso:
Modeling non

isothermal fl
ows and air leakage in a steel
reheating furnace.
Helsinki University of Technology ,Espoo, Finland, 1997,pp 1

7
3. Anne Mäki, Pekka J.Österman and Matti J.Luomala:
Numerical Study of the Pusher type reheating
furnace.
Scand.J. Metall,2002,Vol.31, pp81

87
4
..Combustion
. Phoenics Encyclopaedia, online electric manual of Phoenics software, 1997.
5.
Ra
diation.
Phoenics Encyclopaedia, online electric manual of Phoenics software, 1997.
6.
S.Matsunaga and B.Hiraoka: Transactions ISIJ, 1972,Vol.12, pp.72

78
Nom
enclature
Cp
specific heat at constant pressure
l
C
c
hemical component fraction of species
l
D
diameter of circle fuel burner in the furnace
f
momentum source term in momentum governing equation
EBU
C
empirical constant for eddy

break

up model
H
enthalpy of fluid
refer
H
reference enthalpy of mixture gas at temperature of 298 K
form
H
chemical formation energy of mixture gas with composition derived from iteration
b
I
black

body emissive power at the absolute temperature of mixture gas in the furnace
11
J
I
,
radiation fluxes in
the positive and negative y direction, respectively
k
turbulent kinetic energy of fluid flow
a
K
absorption coefficient of mixture gas in furnace
S
K
scattering coefficient of mixture gas i
n furnace
L
K
,
radiation fluxes in the positive and negative z direction, respectively
fu
m
mass concentration of fuel
2
O
m
mass fraction of oxygen in the furnace
N
M
,
radi
ation fluxes in the positive and negative
z
direction, respectively
P
pressure of fluid in furnace
k
P
volumetric production rate of
k
by shear forces
l
R
source term in chemical concentration conservation equation for species
l
h
S
source term in energy conservation equation
s
stoichiometric requirement for chemical reaction
diffusion coefficient of species
e
kinematics viscosity
,
k
,
1
C
,
2
C
empirical constant in
k
equation
Appendix 1:
Q1 file sample
TALK=F;RUN( 1, 1);VDU=VGAMOUSE
INTEGER(NPATCH)
NPATCH=350
Group 1. Run Title
TEXT(Furnace, consider poles)
****************************************************************************
Model used: ESCRS(EBU)
for combustion, six composite

flux model for radiation,
Thermal data obtained by CHENKIM, K

e model for flow
*****************************************************************************
**********************************************************
******************
Define reference pressure (PRESS, a), reference temperature (RTEMP,K),
gas mole constant(RGAS,J/(mol*K)
****************************************************************************
REAL(RPRESS,RTEMP,RGAS)
RPRESS=1.01325*1.0E5;
RTEMP=298;RGAS=8.3143
****Define some pressure in the furnace OPRESS,Pa)**************************
REAL(OPRESS)
OPRESS=1.0*1.0e5
********Define the mole mass (g/mol)*****************************************
REAL(MO2,MH2,MN2,CH4,MCO,MCH4)
MO2=32.
00;MH2=2.02;MN2=28.01;MCH4=16.04;MCO=28.01
12
**********Temperature of gas and oxygen********************************
REAL(TFU,TOX)
TFU=273+50;TOX=273+350
********************************************************************
Real variable for
coefficients of heat capacity calculation of gas
********************************************************************
REAL(CH4A,CH4B,CH4C,CH4D)
……………………………….(Omitted)
******Methane,CH4, T>298K ***************
CH4A=12.447;CH4B=76.689;CH4C=1.448;
CH4D=

18.004
**********Carbon monoxide, CO********
COA=25.694;COB=8.293;COC=1.109;COD=

1.477
******Hydrogen, H2, T<400K*************
H2A=16.920;H2B=61.459;H2C=0.590;H2D=

79.559
****** Hydrogen, H2, T>400K*************
H2A2=28.280;H2B2=0.41
8;H2C2=0.820;H2D2=1.469
*******Nitrogen N2, T<400K**************
N2A=29.192;N2B=

1.121;N2C=0.000;N2D=3.092
***Nitrogen N2,T>400K************
N2A2=22.552;N2B2=13.209;N2C2=23.130;N2D2=

3.389
*****Oxygen**************
O2A=31.321;O2B=3.895;O2C=

3.105;O2D=

0.335
*******Calculation of heat capacity of CH4 injected from the inlet*****
CCH4=CH4A+CH4B*10**

3*TFU+CH4C*10**5*TFU**

2+CH4D*10**

6*TFU**2
CPCH4=CCH4/MCH4*1e3
*******Calculation of heat capacity of CO injected from the inlet*****
C
CO=COA+COB*10**

3*TFU+COC*10**5*TFU**

2+COD*10**

6*TFU**2
CPCO=CCO/MCO*1e3
*******Calculation of heat capacity of H2 injected from the inlet*****
C1H2=H2A+H2B*10**

3*TFU+H2C*10**5*TFU**

2+H2D*10**

6*TFU**2
CP1H2=C1H2/MH2*1e3
*******Calculation of
heat capacity of N2 injected from the inlet,T<400k*****
C1N2=N2A+N2B*10**

3*TFU+N2C*10**5*TFU**

2+N2D*10**

6*TFU**2
CP1N2=C1N2/MN2*1e3
*******Calculation of heat capacity of N2 injected from the inlet,T>400k*****
C2N2=N2A2+N2B2*10**

3*TOX+N2C2*10**5*
TOX**

2+N2D2*10**

6*TOX**2
CP2N2=C2N2/MN2*1e3
*******Calculation of heat capacity of O2 injected from the inlet*****
CO2=O2A+O2B*10**

3*TOX+O2C*10**5*TOX**

2+O2D*10**

6*TOX**2
CPO2=CO2/MO2*1e3
***************************************************
**********************
Calculate the boundary conditions from the inlets of the up burners in the
Heating zone
*************************************************************************
****************************************************
**********************
Define the gas volume flow rate from the inlets (m^3/h),fuel and oxygen
**************************************************************************
REAL(QFUELy,QOXy)
****************
State the variables of volume rate
of each component in the gas injected
****************
REAL(INCOy,INH2y,INCH4y,INN2y,INO2y)
***************************************************************************
13
State the variables of mass proportion of each component in the gas injec
ted
****************************************************************************
REAL(mINCOy,mINH2y,mINCHy,mINN2y,mINO2y,XMy)
REAL(WGIy,ABRy,CHECKy)
************
Define the whole flow rate, composition of gas and air in all 6 burners
(m3
/s) 350 C
************
fuel 2200 (m3/h) CO=7.4, H2=54.0, CH4=28.0, N2=10.6 t

%
air 11200 (m3/h) O2=79, N2=21
******
REAL(xFCOy,xFH2y,xFCH4y,xFN2y,xOO2y,xON2y)
xFCOy=.074;xFH2y=.54;xFCH4y=.28;xFN2y=.106;xOO2y=.21;xON2y=.79
QFUELy=2200;QOXy=
11200
*****Calculate the composition of mixed flow(gas and air)*******
INCOy=QFUELy*(xFCOy)/(QFUELy+QOXy);INH2y=QFUELy*xFH2y/(QFUELy+QOXy)
INCH4y=QFUELy*xFCH4y/(QFUELy+QOXy);INO2y=QOXy*xOO2y/(QFUELy+QOXy)
INN2y=1.

(INCOy+INH2y+INCH4y+INO2y)
*******Calculate the mass portion in the mixed flow***************
XMy=INCOy*MCO+INH2y*MH2+INCH4y*MCH4+INN2y*MN2+INO2y*MO2
*******Mass portion*************
mINCOy=INCOy*MCO/XMy
mINH2y=INH2y*MH2/XMy
mINCHy=INCH4y*MCH4/XMy
mINN2y=INN2y*MN2/XMy
mI
NO2y=INO2y*MO2/Xmy
************************************
Calcualte enthalpy
************************************
****Mass portion in the fuel********
REAL(FUCOy,FUCH4y,FUH2y,FUN2y,FUMy)
FUMy=xFCOy*MCO+xFH2y*MH2+xFCH4y*MCH4+xFN2y*MN2
FUCOy=xFCOy
*MCO/FUMy
FUH2y=xFH2y*MH2/FUMy
FUCH4y=xFCH4y*MCH4/FUMy
FUN2y=xFN2y*MN2/FUMy
*******Mass portion in the air***********
REAL(OXO2y,OXN2y,OXMy)
OXMy=xOO2y*MO2+xON2y*MN2
OXO2y=xOO2y*MO2/OXMy
OXN2y=xON2y*MN2/OXMy
REAL(HLy,HLFUy,HLOXy,FUy,OXy,HLTOTy)
*****Fuel mass portion in the mixed flow (gas and air)*******
FUy=QFUELy*FUMy/(QFUELy*FUMy+QOXy*OXMy)
*******Oxygen mass portion in the mixed flow*****************
OXy=QOXy*OXMy/(QFUELy*FUMy+QOXy*OXMy)
*******Enthalpy of fuel gas in the mixed flo
w********************
HLFUy=(FUCOy*CPCO+FUCH4y*CPCH4+FUH2y*CP1H2+FUN2y*CP1N2)*TFU*FUy
*******Enthalpy of air in the mixed flow****************
HLOXy=(OXO2y*CPO2+OXN2y*CP1N2)*OXy*TOX
HLTOTy=HLOXy+HLFUy
*****Reference enthalpy at 298 K ***********
***************
REAL(HREFy,CPy)
CPy=mINCOy*CPCO+mINH2y*CP1H2+mINCHy*CPCH4+mINO2y*CPO2+mINN2y*CP1N2
HREFy=(CPy*298)
14
*****Reference chemical formation enthalpy***************
REAL(STDHy,STDSy,HTOTy,FSTDSy,OSTDSy)
STDHy=(mINCOy/MCO*110.541+mINCHy/MCH4*
74.873)*1000*1000
******Enthalpy from the inlet************
HTOTy=HLTOTy

HREFy

STDHy
******Calculate the dynamics variables of the inlet**********
****Define the diameter of inlet (m) and pii
REAL(DIAMy,PII)
DIAMy=.38;PII=3.141
*****Calculat
e the equivalent diameter of a inlet (m)**************
ABRy=((DIAMy/2)**2*PII)**0.5
*****Calculate the flow volume rate ********
REAL(FWGIy,OWGIy)
FWGIy=((QFUELy)/3600/(ABRy**2*6))*(TFU*RPRESS)/(RTEMP*OPRESS)
OWGIy=((QOXy)/3600/(ABRy**2*6))*(TO
X*RPRESS)/(RTEMP*OPRESS)
WGIy=OWGIy+FWGIy
*******Calculate the mole mass and mass portion of fuel and air****
REAL(YMSUMy,OMSUMy,FMSUMy)
YMSUMy=mINCOy/MCO+mINH2y/MH2+mINCHy/MCH4+mINN2y/MN2+mINO2y/MO2
FMSUMy=FUCOy/MCO+FUH2y/MH2+FUCH4y/MCH4+FUN2
y/MN2
OMSUMy=OXN2y/MN2+OXO2y/MO2
****Calculate the mixed gas sensity (kg/m^3)*************
REAL(RHOINy,FRHOy,ORHOy)
FRHOy=OPRESS/(RGAS*TFU*FMSUMy)/1000
ORHOy=OPRESS/(RGAS*TOX*OMSUMy)/1000
RHOINy=(QFUELy*FRHOy+QOXy*ORHOy)/(QFUELy+QOXy)
****Calc
ulate the mass flow rate (kg/(m^2*s))**************
REAL(PRINy)
PRINy=WGIy*RHOINy
****Check Enthalpy*******
real(H1N2OXy,H1N2FUy)
H1N2FUy=(QFUELy*.106/(QFUELy*.106+QOXy*.21))*mINN2y*CP1N2*TFU
H1N2OXy=(QOXy*.21/(QFUELy*.106+QOXy*.21))*mINN2y
*CP1N2*TOX
REAL(H1FUy)
H1FUy=mINCOy*CPCO*TFU+mINH2y*CP1H2*TFU+mINCHy*CPCH4*TFU+H1N2FUy
REAL(H1OXy)
H1OXy=mINO2y*CPO2*TOX+H1N2OXy
REAL(H1INy)
H1INy=H1OXy+H1FUy
*********Calculate turbulent energy KEy and dissipate rate EPy********
REAL(KEy,EPy,TI
NT)
*******Define Turbulent intensity TINT****************
TINT=0.02
KEy=3/2*(TINT*WGIy)**2
EPy=0.09**0.75*KEy**(3/2)/(0.07*(DIAMy/2))
*********************************************************************************
Calculate the boundary c
onditions from the inlets of the lower burners in the
Heating zone
***********************************************************************************
……………………………….(Omitted)
***********************************************************
**********************
Calculate the boundary conditions from the inlets of the burners in the
Soaking zone
***********************************************************************************
15
…………………………………..(Omitted)
******
R
adiation defining
***************
REAL(CPGAS)
CPGAS=1005
****************************************
Define Stefan

Bolzman constant, absorption coefficient, scatter coefficient
and emissivity coefficient
**************************************
***
REAL(GSIGMA, ABSORB, SCAT,EMISS,EG)
GSIGMA=5.6697E

08; ABSORB=0.10; SCAT=0.10; EMISS=0.85;EG=0.20
************************************************************************
Groups 3, 4, 5 Grid Information
* Overall number of cells, RSET(M,N
X,NY,NZ,tolerance)
RSET(M,34,39,59)
* Set overall domain extent:
* xulast yvlast zwlast
name
XSI= 4.250000E+00; YSI= 5.660000E+00; ZSI= 2.807000E+01
RSET(D,CHAM )
* Set objects: x0 y0 z0
* dx
dy dz
name
XPO= 0.000000E+00; YPO= 4.340000E+00; ZPO= 0.000000E+00
XSI= 3.750000E+00; YSI= 1.320000E+00; ZSI= 3.730000E+00
RSET(B,BLK1 )
……………………………….(Omitted)
Group 7. Variables: STOREd,SOLVEd,NAMEd
ONEPHS = T
* Non

default variable names
* Solved variables list
SOLVE(P1 ,U1 ,V1 ,W1 ,H1)
* Stored variables list
STORE(PRPS)
STORE(CP1)
STORE(VIST,DEN1,TMP1,SPH1,YSUM,IMB1)
* Additional solver options
SOLUTN(P1 ,Y,Y,Y,N,N,Y)
SOLUTN(H1 ,Y,Y,Y,
P,P,P)
*Turbulent model
TURMOD(KEMODL)
*Radiation model
RADIAT(ABSORB,SCAT,CP1)
Group 8. Terms & Devices
TERMS (H1 ,N,Y,Y,Y,Y,Y)
Group 9. Properties
ENUL=4.2E

5
PRNDTL(H1)=7.154E

01
*** START OF EXTENDED SCRS MODEL SETTINGS
PRESS0=OPRE
SS
INTEGER(NSPEC,NELEM);NSPEC=7;NELEM=4
INTEGER(NCSTEP,NCREAC);NCSTEP=2;NCREAC=3
SCRS(SYSTEM,NCSTEP,NCREAC,NELEM,FRATE*)
SCRS(SPECIES,CH4,O2,H2,CO,H2O,CO2,N2)
STORE(S1RS,S2RS,S3RS,MMWT)
** Define fuel & oxidizer composition & temperatures
SCRS(FUI
N,mINCHy,mINO2y,mINH2y,mINCOy,0.0,0.0,mINN2y,TFU)
SCRS(OXIN,mINCHt,mINO2t,mINH2t,mINCOt,0.0,0.0,mINN2t,TOX)
SCRS(PROP,CHEMKIN,SCRS)
MESG(2 step 3 reactions finite

rate EBU model
16
MESG(2CH4 + O2 > 2CO+4H2
MESG(2CO + O2 > 2CO2
MESG(2H2 + O2 > 2H2O
*** END OF EXTENDED SCRS MODEL SETTINGS
Group 11.Initialise Var/Porosity Fields
FIINIT(W1 ) = 3
FIINIT(KE ) = 1.000E

02 ;FIINIT(EP ) = 1.079E

02
FIINIT(H1 ) = 2.084E+06 ;FIINIT(RADX) = 1.000E+04
FIINIT(RADY) = 1.000E+04 ;FIINIT(RADZ) = 1
.000E+04
FIINIT(F ) = 1.000
FIINIT(H2 ) = 1.0e

5
FIINIT(CO2) = 1.0e

5
FIINIT(CH4 ) = 1.0e

5
FIINIT(H1 ) = 1500*1005
CONPOR(BLK1 , 0.00,VOLUME,#1,#17,#15,#18,#1,#2)
CONPOR(BLK4 , 0.00,VOLUME,#1,#17,#11,#18,#4,#4)
……………………………….(Omit
ted)
CONPOR(POL2 , 0.00,VOLUME,

#8,

#9,

#2,

#7,

#10,

#10)
……………………………….(Omitted)
INIADD = F
Group 13. Boundary & Special Sources
PATCH (RADISO ,VOLUME,1,31,1,39,1,68,#1,#1)
COVAL (RADISO ,H1 , GRND1 , GRND1 )
COVAL (RADISO ,RADX,
GRND1 , GRND1 )
COVAL (RADISO ,RADY, GRND1 , GRND1 )
COVAL (RADISO ,RADZ, GRND1 , GRND1 )
PATCH (SCRSPRSO,PHASEM,1,31,1,39,1,68,#1,#1)
COVAL (SCRSPRSO,H2 , FIXFLU , GRND2 )
COVAL (SCRSPRSO,CH4 , GRND2 , GRND2 )
COVAL
(SCRSPRSO,CO , FIXFLU , GRND2 )
PATCH (SCRSSRSO,PHASEM,1,31,1,39,1,68,#1,#1)
COVAL (SCRSSRSO,H2 , GRND5 , GRND5 )
COVAL (SCRSSRSO,CO , GRND5 , GRND5 )
******Inlet (Burner)**************
INLET(SCRSFy1,LOW,#14,#15,#17,#17,#5,#5,#1
,#NREGT)
VALUE(SCRSFy1,P1,PRINy)
VALUE(SCRSFy1,V1,WGIy*

.190808)
VALUE(SCRSFy1,W1,WGIy*.981627)
VALUE(SCRSFy1,KE,KEy)
VALUE(SCRSFy1,EP,EPy)
VALUE(SCRSFy1,F,1.)
VALUE(SCRSFy1,CO,mINCOy)
VALUE(SCRSFy1,H2,mINH2y)
VALUE(SCRSFy1,CH4,mINCHy)
VALUE(SCRS
Fy1,H1,HTOTy)
……………………………….(Omitted)
*******Outlet**************
OUTLET(VUOTO1 ,LOW ,1,30,#7,#7,1,1,#1,#1)
COVAL (VUOTO1 ,P1 , FIXVAL , 0.400E+00)
COVAL (VUOTO1 ,KE , 0.000E+00, SAME )
COVAL (VUOTO1 ,EP , 0.000E+00, SAME
)
COVAL (VUOTO1 ,H1 , 0.000E+00, SAME )
COVAL (VUOTO1 ,F , 0.000E+00, SAME )
COVAL (VUOTO1 ,O2 , FIXVAL, 0.233 )
COVAL (VUOTO1 ,N2 , FIXVAL, 0.766 )
17
PATCH (OUTLET ,EAST ,#18,#18,#1,#13,#37,#37,#1,#1)
COVAL (OUTLET ,P1 , FIXVAL
, 0.000E+00)
COVAL (OUTLET ,KE , 0.000E+00, SAME )
COVAL (OUTLET ,EP , 0.000E+00, SAME )
COVAL (OUTLET ,H1 , 0.000E+00, SAME )
COVAL (OUTLET ,F , ONLYMS ,SAME)
**************Boundary wall*************
store(hgas,hre,hst)
RG(1)=1573
PATCH (TOPW1 ,NWALL ,#1,#17,#18,#18,#6,#11,#1,#1)
COVAL (TOPW1 ,U1 , GRND3 , 0.000E+00)
COVAL (TOPW1 ,W1 , GRND3 , 0.000E+00)
COVAL (TOPW1 ,KE , GRND3 , GRND3 )
COVAL (TOPW1 ,EP , GRND3 , GRND3 )
COVAL
(TOPW1 ,H1 , GRND3 , GRND)
……………………………….(Omitted)
RG(121)=1473
PATCH (POL1

LW,LWALL ,#8,#9,#2,#7,#8,#8,#1,#1)
COVAL (POL1

LW,U1 , GRND3 , 0.000E+00)
COVAL (POL1

LW,V1 , GRND3 , 0.000E+00)
COVAL (POL1

LW,KE , GRND3 , GRND3 )
CO
VAL (POL1

LW,EP , GRND3 , GRND3 )
COVAL (POL1

LW,H1 , GRND3 , GRND )
……………………………….(Omitted)
******Radiation boundary conditions***********
PATCH (TOPWA1 ,NORTH ,#1,#17,#18,#18,#6,#11,#1,#1)
COVAL (TOPWA1 ,RADY, 6.667E

01, 5.66978E

08
*(1573**4) )
……………………………….(Omitted)
PATCH (SLABWWA1,WEST ,#17,#17,#8,#8,5,15,#1,#1)
COVAL (SLABWWA1,RADX, 6.667E

01, 5.66978E

08*(1473**4) )
……………………………….(Omitted)
PATCH (POL12EW ,EAST ,#9,#9,#7,#7,#34,#34,#1,#1)
COVAL (POL12EW ,RADX, 6.667E

01, 5.
66978E

08*(873**4) )
Group 15. Terminate Sweeps
LSWEEP = 2000
SELREF = T
RESFAC = 1.000E

02
Group 17. Relaxation
RELAX(P1 ,LINRLX, 8.00000E

01)
RELAX(U1 ,FALSDT, 5.0000E

01)
RELAX(V1 ,FALSDT, 5.000E

01)
RELAX(W1 ,FALSDT, 5.0000E

01)
RELAX(KE ,LINRLX, 5.00000E

01)
RELAX(EP ,LINRLX, 5.00000E

01)
RELAX(H1 ,LINRLX, 8.000000E

01)
RELAX(CO,FALSDT,2.00E

1)
RELAX(H2,FALSDT,2.00E

1)
RELAX(CH4,FALSDT,2.00E

1)
KELIN = 3
**********************************************************
**
************************************************************
Group 19. EARTH Calls To GROUND Station
GENK = T
ASAP = T
RADIA = 1.000E

01 ;RADIB = 1.000E

01
************************************************************
Group 20. Pre
liminary Printout
18
ECHO = T
************************************************************
Group 21. Print

out of Variables
OUTPUT(CH4 ,Y,N,Y,Y,N,N)
OUTPUT(O2 ,Y,N,Y,Y,N,N)
OUTPUT(H2O ,Y,N,Y,Y,N,N)
OUTPUT(CO2 ,Y,N,Y,Y,N,N)
OUTPUT(N2 ,Y,N,Y,Y,N,
N)
OUTPUT(KE ,N,N,N,Y,Y,Y)
OUTPUT(EP ,N,N,N,Y,Y,Y)
OUTPUT(RADX,N,N,N,Y,Y,Y)
OUTPUT(RADY,N,N,N,Y,Y,Y)
OUTPUT(RADZ,N,N,N,Y,Y,Y)
************************************************************
************************************************************
Group 24. Dumps For Restarts
************************************************************
MENSAV(S,RELX,DEF,1.0119E

01,4.0870E+01,1)
MENSAV(S,PHSPROP,DEF,200,0,1.1890E+00,1.0000E

05)
MENSAV(S,FLPRP,DEF,K

E,CONSTANT)
RESTRT(ALL)
STOP
Appendix 2:
Ground file sample
CXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX USER SECTION STARTS:
C
PARAMETER (NLG=100, NIG=200, NRG=200, NCG=100)
C
COMMON/LGRND/LG(NLG)/IGRND/IG(NIG)/RGRND/RG(NRG)/CGRND/CG(NCG)
LOGICAL LG
CHARACTER*4 CG
C
C 2
User dimensions own arrays here, for example:
C DIMENSION GUH(10,10),GUC(10,10),GUX(10,10),GUZ(10)
C
PARAMETER (JDIM=150, IDIM=150)
DIMENSION HLPN2(JDIM,IDIM),HLPO2(JDIM,IDIM),HLPH2(JDIM,IDIM),
+ HLPCH4(JDIM,IDIM),HLPCO(JDIM,
IDIM),HLPCO2(JDIM,IDIM),
+ HLPH2O(JDIM,IDIM),HLPTMP1(JDIM,IDIM),HLP1(JDIM,IDIM),
+ HLP2(JDIM,IDIM),HLP3(JDIM,IDIM)
REAL MO2,MH2,MN2,MCO,MCH4,CH4A,CH4B,CH4C,CH4D,COA,COB,COC,
+ COD,H2A,H2B,H2C,H2D,H2A2,H2B2,H2C2,H2
D2,N2A,N2B,N2C,N2D,
+ N2A2,N2B2,N2C2,N2D2,O2A,O2B,O2C,O2D,CCH4,CCO,CH2,CH2O,
+ CN2,CCO2,CO2,CPCH4,CPCO,CPH2,CPN2,CPO2,CPH2O,CPCO2,HCH4,
+ HH2O,HCO2,H2OA,H2OB,H2OC,H2OD,H2OA2,H2OB2,H2OC2,
+ H2OD2,CO2A,CO2B,CO2C,CO2D,MH2O,
MCO2,TGAS,HWAL,HST,HREF,
+ HGAS
C******Mole mass kg/mol
MO2=32.00/1000
……………………………….(Omitted)
c******The coefficients in the formula of gas heat capacity
c Methane,CH4, T>298
CH4A=12.447
……………………………….(Omitted)
c****
**Standard free enthalpy J/mol
HCH4=74.873*1000
……………………………….(Omitted)
c*****Calculate reference enthalpy at 298 K
c heat capacity of O2 at 298 K
CO2=(O2A+O2B*10.**(

3.)*298.+O2C*10.**5.*298.**(

2.)+
19
+ O2D*10.**(

6.)*298.**2
.)/MO2
c heat capacity of N2 at 298 K
……………………………….(Omitted)
C******Temperature of wall
TW1=RG(1)
……………………………….(Omitted)
C

GROUP 13. Boundary conditions and special sources
C Index for Coefficien
t

CO
C Index for Value

VAL
13 CONTINUE
GO TO (130,131,132,133,134,135,136,137,138,139,1310,
11311,1312,1313,1314,1315,1316,1317,1318,1319,1320,1321),ISC
130 CONTINUE
C

SECTI
ON 12

value = GRND
c gas enthalpy closed to the wall
C Mass portion of O2
CALL GETYX(LBNAME('O2'),HLPO2,JDIM,IDIM)
C Mass portion of N2
CALL GETYX(LBNAME('N2'),HLPN2,JDIM,IDIM)
C Mass portion of H2
CALL GETYX(LBNAME('H2'),HLPH2,JDIM,IDIM)
C Mass portion of CH4
CALL GETYX(LBNAME('CH4'),HLPCH4,JDIM,IDIM)
C Mass portion of CO
CALL GETYX(LBNAME('CO'),HLPCO,JDIM,IDIM)
C Mass portion of CO2
CALL GETYX(LBNAM
E('CO2'),HLPCO2,JDIM,IDIM)
C Mass portion of H2O
CALL GETYX(LBNAME('H2O'),HLPH2O,JDIM,IDIM)
C Temperature of gas
CALL GETYX(LBNAME('TMP1'),HLPTMP1,JDIM,IDIM)
c*****Wall 1, Topwall*****
IF((NPATCH(1:5).EQ.'TOPW1').AND.(IND
VAR.EQ.LBNAME('H1')))THEN
C wall temperature
TW=TW1
L0FVAL=L0F(VAL)
DO 3109 II=1,NX
DO 3109 JJ=1,NY
IICELL=JJ+NY*(II

1)
C mass portion of O2
YO2=HLPO2(JJ,II)
C mass portion of N2
YN2=HLPN2(JJ,II)
C mass portion of H2
YH2=HLPH2(JJ,II)
C mass portion of CH4
YCH4=HLPCH4(JJ,II)
C mass portion of CO
YCO=HLPCO(JJ,II)
C mass portion of CO2
YCO2=HLPCO2(JJ,II)
C
mass portion of H2O
YH2O=HLPH2O(JJ,II)
C Gas temperature
TGAS=HLPTMP1(JJ,II)
c** Calculate heat capacity of different gas component
c Heat capacity of O2 at the wall temperature
CPO2=(O2A+O2B*10
.**(

3.)*TW+O2C*10.**5.*TW**(

2.)
20
+ +O2D*10.**(

6.)*TW**2)/MO2
c Heat capacity of N2 at the wall temperature
IF(TW.LE.400.0)THEN
CPN2=(N2A+N2B*10.**(

3.)*TW+N2C*10.**5.*TW**(

2.)
+ +N2D*10.**(

6.)*TW
**2)/MN2
ELSE
CPN2=(N2A2+N2B2*10.**(

3.)*TW+N2C2*10.**5.*
+ TW**(

2.)+N2D2*10.**(

6.)*TW**2)/MN2
ENDIF
c Heat capacity of H2 at the wall temperature
IF(TW.LE.400.0)THEN
CPH2=(H2A+H2
B*10.**(

3.)*TW+H2C*10.**5.*TW**(

2.)
+ +H2D*10.**(

6.)*TW**2)/MH2
ELSE
CPH2=(H2A2+H2B2*10.**(

3.)*TW+H2C2*10.**5.*TW**(

2.)
+ +H2D2*10.**(

6.)*TW**2)/MH2
ENDIF
c Heat capacity of CH4 at
the wall temperature
CPCH4=(CH4A+CH4B*10.**(

3.)*TW+CH4C*10.**5.*TW**(

2.)
+ +CH4D*10.**(

6.)*TW**2)/MCH4
c Heat capacity of CO at the wall temperature
CPCO=(COA+COB*10.**(

3.)*TW+COC*10.**5.*TW**(

2.)
+
+COD*10.**(

6.)*TW**2)/MCO
c Heat capacity of H2O at the wall temperature
IF(TW.LE.600.0)THEN
CPH2O=(H2OA+H2OB*10.**(

3.)*TW+H2OC*10.**5.*TW**(

2.)
+ +H2OD*10.**(

6.)*TW**2)/MH2O
ELSE
CPH2O=(H2OA2+H2OB2*10.**(

3.)*TW+H2OC2*10.**5.*TW**(

2.)
+ +H2OD2*10.**(

6.)*TW**2)/MH2O
ENDIF
c Heat capacity of CO2 at the wall temperature
CPCO2=(CO2A+CO2B*10.**(

3.)*TW+CO2C*10.**5.*TW**(

2.)
+
+CO2D*10.**(

6.)*TW**2)/MCO2
C**
c Gas enthalpy at wall temperature
HGAS=(yO2*CPO2+yH2*CPH2+yN2*CPN2+yCH4*CPCH4+yCO*CPCO+yH2O
+ *CPH2O+yCO2*CPCO2)*TW
c Reference enthalpy of gas at temperature of 298K
HREF=(
yO2*CO2+yH2*CH2+yN2*CN2+yCH4*CCH4+yCO*CCO+yH2O*CH2O+
+ yCO2*CCO2)*298.
c Chemical formation enthalpy
HST=(yCH4*HCH4/MCH4+yCO*HCO/MCO+yH2O*HH2O/MH2O+yCO2*
+ HCO2/MCO2)
c Enthalpy of the boundary wall
HWAL=HGAS

HREF

HST
F(L0FVAL+IICELL)=HWAL
HLP1(JJ,II)=HGAS/TW
HLP2(JJ,II)=HST
HLP3(JJ,II)=HREF/298.
3109 CONTINUE
ENDIF
…………………………………(Other boundary walls)
c
CALL SETYX(LBNAME('HGAS'),HLP1,JDI
M,IDIM)
CALL SETYX(LBNAME('HST'),HLP2,JDIM,IDIM)
CALL SETYX(LBNAME('HRE'),HLP3,JDIM,IDIM)
RETURN
1312 CONTINUE
Comments 0
Log in to post a comment