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:
I and J, the radiation fluxes (in W/m^2) in the positive and negative y (or r) directions, respectively;
K and L, the radiation fluxes (in W/m^2) in the positive and negative z directions, respectively; and
M and N, the radiation fluxes (in W/m^2) in the positive and negative x directions, respectively.
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/y**j}*d(I*y**j)/dy = -(a+s)*I + a*E + s*SUM + j*J/y ...(4.1) {1/y**j}*d(J*y**j)/dy = (a+s)*J - a*E - s*SUM + j*J/y ...(4.2) dK/dz = -(a+s)*K + a*E + s*SUM ...(4.3) dL/dz = (a+s)*L - a*E - s*SUM ...(4.4) dM/dx = -(a+s)*M + a*E + s*SUM ...(4.5) dN/dx = (a+s)*N - a*E - s*SUM ...(4.6)
where j=0 in Cartesian coordinates, but,
in polar coordinates, j=1 and
d/dx={1/y}*d/dtheta.
The term SUM in equations (4.1) to (4.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 ...(4.7)
whereas in two-dimensional systems just four fluxes are present, and in one-dimensional problems only two fluxes.
In equations (4.1) to (4.6):
E = S*T**4 ...(4.8)
where the Stefan-Boltzmann constant S=5.6678E-8 W/m^2/K^4.
Each of the differential flux equations (4.1) to (4.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 ...(4.9)
With these definitions, the six fluxes appearing in equations (4.1) to (4.6) can be eliminated so as to produce the following three second-order ordinary differential equations:
{1/y**j}*d[{y**j/(a+s+j/y)}*dRy/dy]/dy = -a*(Ry-E) - s*[2*Ry-Rx-Rz]/3 ...(4.10) d[{1/(a+s)}*dRz/dz]/dz =-a*(Rz-E) - s*[2*Rz-Ry-Rx]/3 ...(4.11) d[{1/(a+s)}*dRx/dx]/dx =-a*(Rx-E) - s*[2*Rx-Ry-Rz]/3 ...(4.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:
Srad = - d(Qri)/dxi ...(4.12a)
where: Srad is the volumetric energy source in W/m^3; and Qri ( in W/m^2 ), is the i-direction component of 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 (4.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 ...(4.15)
Qry = K - L = - {2/(a+s+j/y)}*dRy/dy ...(4.13)
Qrz = M - N = - {2/(a+s)}*dRz/dz ...(4.14)
Thus, from equation (4.12a) the contribution of the radiation fluxes to the energy equation source term is:
Srad = 2*a*[sum(Ri) - n*E] ...(4.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], Gosman and Lockwood [1973a, 1973b] and Khalil [1982]), and also in fire simulations (see for example Hoffmann and Markatos [1988]).
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:
(a) radiation is transmitted in coordinate directions only;
(b) no interlinkages, apart from scattering, arise between the radiation fluxes in the respective coordinate directions;
(c) the scattering term presumes that the scattering is isotropic with angle, which is probably reasonable only if the total contribution of scattering is not large; and
(d) the model is not readily extended to co-ordinate systems which are neither cartesian nor cylindrical-polar.
It is interesting to consider the diffusion approximation to the composite-flux equations (4.10) to (4.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 (4.13) to (4.15) is:
Qri = - [8*S*T**3/(a+s)]*dT/dxi ...(4.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-relecting 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 = eg*S*Tin**4 ...(4.18)
where eg is the inlet gas emissivity and Tin is the inlet gas temperature.
From equations (4.11) and (4.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 ...(4.19)
which from equation (4.9) for Rz and equation (4.18) reduces to:
Sr = (Kin - Rz) ...(4.20)
Equation (4.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 = ew*Ew + (1-ew)*K ...(4.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 = [ew/(2-ew)]*(Ew - Rz) ...(4.22)
In the foregoing, ew is the wall emissivity and Ew the black-body emissive power at the wall temperature Tw, i.e. Ew=S*Tw**4.
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:
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, ew (the emissivity of the wall), S, Tw, Tin and eg (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 (4.8), may be printed in the RESULT file or viewed via PHOTON and AUTOPLOT.
wbs