A major limitation of the composite-flux model is that it is not readily extended to BFCs.

In order to remove this deficiency PHOENICS has been equipped with an alternative diffusional radiation model, in which the radiative heat transfer is solved by way of an equation for the "composite radiosity".

The radiosity represents the average of the incoming and outgoing radiation fluxes over all directions of the solid angle.

This radiosity equation was devised intuitively by Spalding [1994], but the resulting equation has the same form as that used for the irradiance in the P-1 spherical-harmonics approximation to the RTE ( see for example Ozisik [1973], Viskanta and Menguc [1987], and Siegel and Howell [1992] ).

The equation to be solved for the composite radiosity R (in W/m^2) is:

d/dxi [{4/3(a+s)}*dR/dxi ] + 4*a*(E - R) = 0 (5.1)

which follows in some sense from simplification of the composite- flux model by ignoring the directional aspects of the radiation.

In effect the RTE has been integrated over all directions of the solid angle, and so scattering source terms make no appearance due to cancellation.

The physical meaning of this equation is that it represents the net rate of loss or gain of radiant energy per unit volume. The term a*E represents the local rate of emission, while the term a*R represents the local rate of absorption of radiation per unit volume.

The net radiative heat fluxes in the three coordinate directions, are given by:

Qri = - [4/{3*(a+s)}]*dR/dxi (5.2)

so that from equations (4.12a) and (5.1) the contribution of the radiation to the energy equation source term is:

Srad = 4*a*[R - E] (5.3)

As with the composite-flux model, the radiosity model is very accurate for optically thick media, i.e. if the optical thickness is greater than say 2.

The model will yield inaccurate results for thinner media, especially near boundaries, and also if the radiation field is anisotropic.

For an optically-thick medium, equation (5.1) yields R=E, so that equation (5.2) reduces to the diffusion approximation

Qri = - [16 * S * T**3 / {3 * (a+s)}] * dT/dxi (5.4)

If s=0, then the absorption coefficient a is exactly equal to the Rosseland mean absorption coefficient aR ( see Siegel and Howell [1992] ).

Finally, it is interesting to compare the radiosity equation (5.1) with that for the irradiance G in the P1-approximation, as presented by, for example, Liu and Swithenbank [1991]:

d/dxi [{1/(3*Ke)}*dG/dxi] + Ka*(4*(n**2)*E - G) = 0 (5.5)

where:

- Ke is an effective extinction coefficient,
- n is the refractive index of the medium, and
- Ka is the absorption coefficient.

- the radiosity is related to the irradiance via R=G/(4*(n**2)),
- Ka=a, and
- the effective extinction coefficient Ke is related to a and s by Ke=(a+s).

In the P-1 model, the radiative heat flux vector is given by

Qri= - [1/{3*Ke}]*dG/dxi (5.6)

so that equations (5.1) and (5.2) reveal that the refractive index has been taken as unity in the radiosity model, which is applicable for all gases ( see Ozisik [1973] ).

Liu and Swithenbank [1991] define Ke = Ka + (1-g)*Ks where Ks is a scattering coefficient and g is a symmetry factor which represents the amount of radiation scattered in the forward direction.

The quantity (1-g)*Ks is an effective scattering coefficient; and when the symmetry factor g=1, 0, -1 it corresponds to complete forward scattering, isotropic scattering, and complete backward scattering, respectively.

The boundary conditions required for the radiosity are very similar to those described earlier for the composite-flux model in Section 4.5. For symmetry planes and perfectly-reflecting boundaries, the treatment is identical.

For wall boundaries with a prescribed temperature, the boundary source term for R per unit area is given by:

Sr = 2*[ew/(2-ew)]*(Ew - R) (5.7)

which is the same as that given by equation (4.22) for the flux model, but follows from the diffusional-slip boundary condition presented by Deissler [1964] for the Rosseland diffusion model.

It follows that, if the net radiative wall heat flux Qr is known, then Sr can simply be set equal to Qr.

It should be mentioned that the wall boundary condition (5.7) is equivalent to the Marshak-type boundary condition used for the irradiance in the P1-approximation ( see for example Ozisik [1973] and Liu and Swithenbank [1991] ).

For non-reflecting boundaries, such as inlets and free boundaries, the following boundary-source term per unit area is proposed:

Sr = (eg*S*Tin**4 - R) (5.8)

which is the same as equation (4.20) for the flux model.

The composite-radiosity radiation model is activated by inserting the following PIL command in the Q1 file:

RADIAT(RADI,ABSORB,SCAT,Energy)

where the arguments have the same meanings as given earlier in Section 4.6 for the composite-flux model. The command activates solution of the radiosity variable SRAD and in common with the flux model:

- the use of absolute temperature is advised, although the use of TEMP0 is allowed;
- the STORE(EMPO) facility is available; and
- the boundary conditions at walls and non- reflecting boundaries must be introduced explicitly by the user.

For a wall boundary, if the wall temperature is known, the following PIL settings are required:

PATCH(WALLR,EAST,NX,NX,1,NY,1,NZ,1,1)

COVAL(WALLR,SRAD,2*EMIW/(2.0-EMIW),GSIGMA*TWAL**4)

whereas if the net radiative heat influx is known:

COVAL(WALLRA,SRAD,FIXFLU,QRAD)

For a non-reflecting boundary, such as a flow inlet, the following settings are required:

PATCH(RADIN,WEST,1,1,1,NY,1,NZ,1,1) COVAL(RADIN,SRAD,1.0,GSIGMA*EMIGIN*TIN**4)

where in all of the foregoing examples, EMIW, GSIGMA, TWAL, TIN and EMIGIN are user-defined PIL variables (see Section 4.6).

wbs