The flux radiation model incorporated into PHOENICS is the composite-flux model of Spalding [1980], which originates from the works of Hamaker [1947a,1947b] and Schuster [1905].
In this model, the discretisation of angle is such that the effects of radiation are accounted for by reference to the positive and negative radiation fluxes in each of the coordinate directions, as follows:
Each of these fluxes may be interpreted as the total diffusely-distributed radiation crossing a normal plane.
For one-dimensional problems only 2 fluxes are employed, one outward and one inward.
Two-dimensional problems are handled by 4 fluxes, and three-dimensional problems by 6 fluxes.
The differential equations describing the variations of the six radiation fluxes are:
{1/yj}*d(I*yj)/dy = -(a+s)*I + a*E + s*SUM + j*J/y Eqn 5.1
{1/yj}*d(J*yj)/dy = (a+s)*J - a*E - s*SUM + j*J/y Eqn 5.2
dK/dz = -(a+s)*K + a*E + s*SUM Eqn 5.3
dL/dz = (a+s)*L - a*E - s*SUM Eqn 5.4
dM/dx = -(a+s)*M + a*E + s*SUM Eqn 5.5
dN/dx = (a+s)*N - a*E - s*SUM Eqn 5.6
where j=0 in Cartesian coordinates, but in polar coordinates, j=1 and d/dx={1/y}*d/dθ.
The term SUM in equations (5.1) to (5.6) is the mean of all the fluxes in the system. For example, in three-dimensional systems:
SUM = (I + J + K + L + M + N)/6 . Eqn 5.7
whereas in two-dimensional systems just four fluxes are present, and in one-dimensional problems only two fluxes.
In equations (5.1) to (5.6):
E = σT4 Eqn 5.8
where the Stefan-Boltzmann constant σ=5.6678 10-8 Wm-2K-4.
Each of the differential flux equations (5.1) to (5.6) expresses the diminution of a flux with distance in consequence of absorption and scattering, and its augmentation by emission and scattering from other directions.
The composite radiation fluxes are defined as:
Rx = (I+J)/2 ; Ry = (K+L)/2 ; and Rz = (M+N)/2 . Eqn 5.9
With these definitions, the six fluxes appearing in equations (5.1) to (5.6) can be eliminated so as to produce the following three second-order ordinary differential equations:
{1/yj}*d[{yj/(a+s+j/y)}*dRy/dy]/dy = -a*(Ry-E) - s*[2*Ry-Rx-Rz]/3 . Eqn 5.10 d[{1/(a+s)}*dRz/dz]/dz =-a*(Rz-E) - s*[2*Rz-Ry-Rx]/3 . Eqn 5.11 d[{1/(a+s)}*dRx/dx]/dx =-a*(Rx-E) - s*[2*Rx-Ry-Rz]/3 . Eqn 5.12
The form of the scattering term in the foregoing equations is valid for three-dimensional systems only.
For one-dimensional problems the scattering term is zero; and in two-dimensional systems the term reduces to s*[Ri-Rj]/2, where Ri is the composite flux in the solved-for direction, and Rj is the flux in the other coordinate direction.
The contribution of radiative heat transfer to the energy equation H1 or TEM1 is a source term involving the divergence of the radiative heat-flux vector, Qri, as follows:
Sr = - ∇Qr Eqn 5.12a
where Sr is the volumetric energy source in W/m3; and Qr ( in W/m2 ), is the radiative heat flux vector.
However, it is also more easily understood and computed as the negative of the sum of all the E-related sources of Rx, Ry and Rz, as in equation (5.16) below.
Once the composite radiation fluxes are determined, the net radiative heat fluxes in the three coordinate directions, can be evaluated from:
Qrx = I - J = - {2/(a+s)}*dRx/dx . Eqn 5.13
Qry = K - L = - {2/(a+s+j/y)}*dRy/dy . Eqn 5.14
Qrz = M - N = - {2/(a+s)}*dRz/dz . Eqn 5.15
Thus, from equation (5.12a) the contribution of the radiation fluxes to the energy equation source term is:
Sr = 2a[Σ(Ri) - n*E] . Eqn 5.16
where the summation includes all the composite fluxes present in the system, and n is the number of these fluxes.
The composite flux model has been widely applied in combustion chambers and furnaces (see for example Patankar and Spalding [1974]a>, Gosman and Lockwood [1973a, 1973b] and Khalil [1982]), and also in fire simulations (see for example Hoffmann and Markatos [1988]a>).
While the composite-flux model has the advantage of simplicity and computational economy and has resulted in many satisfactory predictions, it does suffer from a number of limitations, including:
It is interesting to consider the diffusion approximation to the composite-flux equations (5.10) to (5.12). This is valid when the absorption coefficient a is large compared with the inverse flow dimension, so that local radiative equilibrium exists (E=R=Ri). The result from equations (5.13) to (5.15) is:
Qr = - (8/3)σT3(a+s)-1]∇T . Eqn 5.17
If s=0, the flux absorption coefficient is related to the Rosseland mean absorption coefficient aR (see Siegel and Howell [1992]) via a=1.5*aR.
At these boundaries, the radiative heat flux is zero by definition; this is the default boundary condition in PHOENICS for the radiation equations. The user is reminded that a perfectly-reflecting boundary is one which has a surface emissivity of zero.
At non-reflecting boundaries, such as openings or free boundaries, the outgoing radiation leaves the calculation domain without reflection, and only the incoming radiation needs to be prescribed. Thus, for example, the incoming flux at a low inlet boundary may be defined as:
Kin = εgσTin4 . Eqn 5.18
where εg is the inlet gas emissivity and Tin is the inlet gas temperature.
From equations (5.11) and (5.14) the required inlet-boundary source term per unit low area for the composite flux Rz is:
Sr = -{1/(a+s)}dRz/dz |in = Qz/2 . Eqn 5.19
which from equation (5.9) for Rz and equation (5.18) reduces to:
Sr = (Kin - Rz) . Eqn 5.20
Equation (5.19) also defines the required boundary source term for a wall, so that if the net radiative wall heat flux Qz is known, then Sr can simply be set equal to one half of the prescribed wall heat flux. The convention is that positive values of Qz denote that radiant energy is entering the system.
If the wall temperature is prescribed rather than Qz, then the required expression for Qz may be derived by considering the radiation flux leaving the wall:
L = εwEw + (1-εw)K . Eqn 5.21
where the first contribution is the radiation emitted and the second is the radiation reflected. Then, via the definition of Qz and Rz, the following relation is obtained:
Sr = Qz/2 = [εw/(2-εw)](Ew - Rz) . Eqn 5.22
In the foregoing, εw is the wall emissivity and Ew the black-body emissive power at the wall temperature Tw, i.e. Ew=σTw4.
The composite-flux radiation model is activated by inserting the following PIL command in the Q1 file:
RADIAT(FLUX,ABSORB,SCAT,Energy)
where:
The RADIAT command activates solution of Rx (=RADX), Ry (=RADY) and Rz (=RADZ) as required, but for a more detailed description the user is referred to the PHOENICS Encyclopaedia entry 'RADIAT'.
The user is advised to employ absolute temperature in degrees K, although the implementation does support the use of the PIL variable TEMP0, which may be used to define a reference temperature, i.e. TEMP0=273K.
The boundary conditions at walls and non-reflecting boundaries must be introduced explicitly by the use of PATCH and COVAL statements.
Thus, for example, if the wall
temperature is known at a wall boundary, the following settings are appropriate:
PATCH(WALLR,EAST,NX,NX,1,NY,1,NZ,1,1) COVAL(WALLR,RADX,EMIW/(2.0-EMIW),GSIGMA*TWAL**4)
whereas, if the net radiative heat influx is known, the COVAL setting should be:
COVAL(WALLRA,RADX,FIXFLU,0.5*QRAD)
For a non-reflecting boundary, such as a flow inlet at a west boundary, the following settings are required:
PATCH(RADXIN,WEST,1,1,1,NY,1,NZ,1,1) COVAL(RADXIN,RADX,1.0,GSIGMA*EMIGIN*TIN**4)
where, in all of the foregoing examples, EMIW, GSIGMA, TWAL, TIN and EMIGIN are user-defined PIL variables representing, respectively, εw (the emissivity of the wall), σ, Tw, Tin and ε g (the emissivity of the gas at inlet).
Finally, when STORE(EMPO) appears in the Q1 file, the emissive power of the gas E, as defined by equation (5.8), may be printed in the RESULT file or viewed via PHOTON and AUTOPLOT.
A.D.Gosman and F.C.Lockwood, 'Incorporation of a flux model for radiation into a finite-difference procedure for furnace calculations', 14th Symp. Combustion, p661, (1973a).
A.D.Gosman and F.C.Lockwood, 'Prediction of the influence of turbulent fluctuations on flow and heat transfer in furnaces', Mech. Eng. Dept. Report HTS/73/52, Imperial College, London, (1973b).
H.C.Hamaker,'Radiation and heat conduction in light-scattering material: Part I', Philips Res. Rep., Vol.2, p55, (1947a).
H.C.Hamaker,'Radiation and heat conduction in light-scattering material: Part II', Philips Res. Rep., Vol.2, p103, (1947b).
N.Hoffmann and N.C.Markatos, 'Thermal radiation on fires in enclosures', Appl. Math. Modelling, Vol.12, p129, (1988).
E.E.Khalil, 'Modelling of furnaces and combustors', Abacus Press, (1982).
E.E.Khalil, D.B.Spalding and J.H.Whitelaw, 'The calculation of local flow properties in two-dimensional furnaces', Int. J.Heat Mass Transfer, Vol.18, p775, (1975).
S.V. Patankar and D.B.Spalding, 'Simultaneous Prediction of Flow Patterns and Radiation for Three-Dimensional Flames', presented at the 1973 Seminar of the International Center for Heat and Mass Transfer, Yugoslavia.
S.Rosseland, 'Theoretical astrophysics', Oxford Univ. Press, Clarendon, London and New York, (1936).
R.Siegel and J.R.Howell, 'Thermal Radiation Heat Transfer', 3rd Edition, Hemisphere Publishing Corporation, (1992).
A.Schuster, 'Radiation through a foggy atmosphere', J.Astrophys., Vol. 21, p1, (1905).
D.B.Spalding, 'Idealisations of radiation', In Mathematical Modelling of Fluid-Mechanics, Heat-Transfer and Chemical-Reaction Processes, Lecture 9, HTS/80/1, Imperial College, Mech. Engng., Dept., London, (1980)