**Contents**

- The main features
- The rationale
- How
*W*_{gap}is calculated - Implementation in PHOENICS
- An early example
- More recent examples
- Concluding remarks regarding IMMERSOL
- References

*T*_{1}, the temperature of the first phase*e.g.*air;*T*_{2}, the temperature of the second phase, if present,*e.g.*a cloud of solid paricles suspended within the air; and*T*_{3}, the 'radiosity temperature' defined below.

If *R *stands for the radiosity, *i.e.* the sum all radiation fluxes traversing the volume, regardless of direction and wavelength, then
*T*_{3} is related to it by the equation:

where σ is the Stefan-Boltzmann constant, ( = 5.670373x10

R= σT_{3}^{4}

T_{3}=(R/σ)^{1/4}Eq (1)

∂(

c_{1}ρ_{1}T_{1})/∂t+ div(c_{1}ρ_{1}v_{1}T_{1}) = div(λ_{1}gradT_{1}) +S_{1,2}+S_{1,3}Eq (2)

and

∂(

c_{2}ρ_{2}T_{2})/∂t+ div(c_{2}ρ_{2}v_{2}T_{2}) = div(λ_{2}gradT_{2}) +S_{2,1}+S_{2,3}Eq (3)

wherein *S*_{1,2} and *S*_{2,1} represent energy transfers per unit volume netween phases 1 and 2, and *v*, the velocity vector.
*S*_{1,2} and *S*_{2,1} volumetric radiative heat absorption and emission. Here the λs represent the thermal conductivities, laminar plus turbulent, of the respective phases. The significance of the other symbols are: *c* specific heat capacity, ρ density, *t* time and
** v** velocity vector.

*T*_{3} obeys a similar equation but with fewer terms. Specifically it has a zero on the left-hand side, because radiation is not convected In either time or space. It is:

0 = div(λ

_{3}gradT_{3}) +S_{3,1}+S_{3,2}(4)

The IMMERSOL presumption is that they are related to the three temperatures *via* the following equations:

S_{1,3}= -S_{3,1}= 4 ε'_{1}σ(T_{3}^{4}-T_{1}^{4}) (5)

and

S_{2,3}= -S_{3,2}= 4 ε'_{2}σ(T_{3}^{4}-T_{2}^{4}) (6)

wherein ε'_{1} and ε'_{2} are the emissivities of the respective phases per unit length. These quantities are supposed numerically equal to the absorptivities, which measure the proportion of the radiation flux which is absorbed per meter of its passage through the medium in question.

λ

_{3}= 4 σT_{3}^{3}/ { 0.75( ε'_{1}+s'_{1}+ ε'_{2}+s'_{2}) + 1/W_{gap}} (7)

At open boundaries of the domain, nett radiation fluxes are ordinarily prescribed.

In that book, and elsewhere, two other accepted ideas are described which IMMERSOL has incoporated, namely those of the **'optically-thick'** and **'optically-thin' **limits. Here 'thick' and 'thin' compare
the size of the gap between the solid walls enclosing the transparent medium with what can be called the 'mean free path' of radiation, *i.e.* the
reciprocal of ε'+*s*'.

What is distinctive about IMMERSOL is that it is valid both for and **between** those limits.

Both extremes arise in practice. Within a large coal-fired furnace, the cloud of burning particles and gaseous combustion products can be regarded as optically thick; for so much solid surface is present per unit volume that rays emanating from the middle of the furnace must be multiply reflected before they escape to the water-cooled walls. The air within the rooms and corridors of inhabited buildings, by contrast, constitutes an optically thin medium; wall-to-wall radiation suffers no impediment.

For optically thick media, there exists a formula which connects the radiative heat flux ** q**, in W/m

= -(4/3)(ε'+qs')^{-1}σ grad(T^{4}) (8)

Here *T* is the local temperature of the transparent medium.

If the equation is expressed in terms of an effective thermal conductivity,
λ_{eff}, involving grad *T* rather than
grad(*T*^{4}),
the expression for λ_{ref} becomes:

λ

_{eff}= (16/3) (ε'+s')^{-1}σT^{3}(9)

At the other extreme, when the medium is so thin as not to participate at all in the radiative heat transfer between two solid sufaces, at temperatures *T*
_{hot} and *T*_{cold}, say, the heat flux *q* is well known to obey the formula:

={1 + (1-εq_{hot})/ε_{hot}+ (1-ε_{cold})/ε_{cold}}^{-1}σ(T_{hot}^{4}-T_{cold}^{4})

={1/ε_{hot}+ 1/ε_{cold}- 1}^{-1}σ(T_{hot}^{4}-T_{cold}^{4}) (10)

Equation (8) is of the flux-proportional-to-gradient kind which PHOENICS is well-equipped to solve. Equation (10) is of the less amenable action-at-a distance kind. The question arises how can the latter be made more like the former?

= 4σqT^{3}(T_{hot}-T_{cold}) (11)

where *T* stands for either temperature

Since the temperature gradient equals (*T*_{hot} -
*T*_{cold})/*W*_{gap}, the effective conductivity
which corresponds to equation (11) is simply:

λ

_{eff}= 4W_{gap}σT^{3}(12)

where *W*_{gap} stands for the distance between the solid surfaces. So the conductivity increases with inter-wall distance; as it must do if the heat flux is to be independent of that distance.

It is interesting to compare the value of this conductivity with the thermal conductivities of common materials such as:

atmospheric air | water at 0 degC | steel |

0.0258 | 0.569 | 43.0 |

In the same units, and with a wall gap equal to one metre, the values of λ_{eff} at various temperatures in degrees Celsius are:

T_{3} |
20 | 100 | 500 | 1000 | 1500 | 2000 | |

λ_{eff}
| 5.706 | 11.77 | 104.8 | 467.9 | 1264.1 | 2663.6 |

λ

_{eff}^{-1}= (3/16)(ε'+s')/(σT^{3}) (13)

and, for the thinnest-possible totally-empty medium, eq. (12) yields:

λ

_{eff}^{-1}= 1/(4W_{gap}σT^{3}) (14)

It is therefore not unreasonable to suppose that, for intermediate
conditions, the two multipliers of 4σ *T*^{3} should
be **added** so as to create a more-generally-valid single resistivity formula thus:

λ

_{eff}^{-1}= { (3/4) (ε'+s') + 1/W_{gap}}/(4σT^{3}) (15)

This is the origin of equation (7), introduced in section 4.1(d) above.

As is usual in PHOENICS, the conductivities pertaining to the cell are stored at each grid node; Therefore the radiative heat flux q_{r,B-C} crossing the boundary between cell B and cell C will be deduced from the formula:

q

_{r,B-C}= (T_{3,B}-T_{3,C})/{(x_{I}-x_{B})/λ_{3,B}+ (x_{C}-x_{I})/λ_{3,C})} Eqn (15a)

where *x* is the horizontal co-ordinate and the λ's are given by equation (15) above.

However, the calculation of the radiant flux at the S interface
requires more careful study because the surface emissivity can cause
a discontinuity of
*T*_{3} gradient there, as is illustrated in the
following figure, which shows the postulated profiles of both *T*_{3} and *T*_{1} because of their inescapable interaction.

Here it is postulated that *T*_{3} and
*T*_{1} are equal to each other within the solid but,
whereas the latter has a finite gradient everywhere, the latter may
have an infinite one at the interface between the phases.

The fluxes of energy in question are as follows:

- conduction from A to S, namely:
q

_{A-S}= (*T*_{1,A}-*T*_{1,S})/R_{A-S}Eqn 15(b)

R_{A-S}= (*x*_{S}-*x*_{A})/λ_{1,A}Eqn 15(c) - conduction and convection from S to B, namely:
q

_{c,S-B}= (*T*_{1,S}-*T*_{1,B})/R_{c,S-B}Eqn 15(d)

R_{c,S-B}= (*x*_{B}-*x*_{S})/λ_{1,B}Eqn 15(e) - radiation from S to B, namely:
q

_{r,S-B}= (*T*_{3,S}-*T*_{3,B})/R_{r,S-B}Eqn 15(f)

R_{r,S-B}= (*x*_{B}-*x*_{S})/λ_{3,B}+(1-ε_{S})/[ε_{S}4σT_{3}^{3}] Eqn 15(g)wherein the term involving ε

_{S}is inserted so as to conform with equation (10) above.

q

_{A-S}= q_{c,S-B}+ q_{r,S-B}Eqn 15(h)

and with *T*_{1S}=*T*_{3S}, the necessary formula can be
deduced from eqn 15(h) as follows:

T_{1,S}=T_{3,S}= (T_{1,A}/R_{A-S}+T_{1,B}/R_{c,S-B}+T_{3,B}/R_{r,S-B})/(1/R_{A-S}+1/R_{c,S-B}+1/R_{r,S-B}) Eqn(16)

where the resistances R_{A-S}, R_{c,S-B} and R_{r,S-B}
were defined by equations 15(c), 15(e) and 15(g) earlier. Inspection of the R_{r,S-B} term in equation 15(g) shows that
(*x*_{B}-*x*_{S})/λ_{3,B} is the
conventional resistance and that (1-ε_{S})/(ε_{S}4σT_{3}^{3})
is the extra resistance caused by the non-unity emissivity.

From equations 15, 15(f) and 15(g) above, the radiative heat flux at the wall can now be written as:

q

_{r,S-B}= Γ_{S-B}(T_{3,S}-T_{3,B}) Eqn (17)

where

Γ

_{S-B}= 4σT_{3,m}^{3}/[(x_{B}-x_{S})(0.75(ε^{'}+s^{'})+1/W_{gap}) + (1-ε_{w})/ε_{w}] Eqn (18)

and T_{3,m} is a mean temperature between S and B, which is evaluated from:

T

_{3,m}^{3}=(T_{S}^{2}+T_{B}^{2})(T_{S}+T_{B})/4 Eqn (19)

This expression follows from linearising the wall heat flux, as follows:

q

_{r,S-B}= σβ(T_{S}^{4}-T_{B}^{4}) = σβ(T_{S}^{2}+T_{B}^{2})(T_{S}^{2}-T_{B}^{2}) q_{r,S-B}= 4σβT_{3,m}^{3}(T_{S}-T_{B}) Eqn (20)

where by comparison with equation (18) above it is clear β is defined by:

β = [(x

_{B}-x_{S})(0.75(ε^{'}+s^{'})+1/W_{gap}) + (1-ε_{w})/ε_{w}]^{-1}Eqn (21)

Few actions are needed in order to activate IMMERSOL. They comprise:

SOLVE(T3) to ensure that T3 is solved

DISWAL to activate the WDIS and WGAP calculation

STORE(WGAP) to enable this quantity to be accessed

together with such information as provides material properties and initial and boundary conditions.

Examples may be found in the IMMERSOL-related cases, in the advanced-radiation-option library.

From PHOENICS-3.4 onwards, however, the setting of radiation-related properties has been made more similar to the setting of other properties, such as density and specific heat.

Thus:

- there now exist two new PIL variables with meaningful names, namely:
- EMISS,
- the emission (and absorption) coefficient of a transparent medium, with the dimensions of 1/length; and
- SCATT,
- the scattering coefficient of the transparent medium, also with the dimensions of 1/length

- If these quantities are set to non-negative values, those values
are used for all regions which are not occupied by solids, unless
further instructions are given.
- Negative values are treated as zero, being without either physical
or formula-choosing significance[
**unlike**the case of RHO1, etc]. - If the commands:

STORE(EMIS) and/or STORE(SCAT)

are placed in the Q1, for print-out or other purposes, the fields of values are printed in the result file, and stored in the phi file, in the usual manner.Moreover

**different**values can then be placed in different parts of the domain by the conventional use of:

FIINIT( ),

PATCH(name,INIVAL, ),

INIT(name, EMIS, ) and

INIT(name, SCAT, ) - If InForm-style property formulae are the provided in the Q1,
it is the values which the formulae provide which are used in the
radiation calculations.
- PLANT-generated or user-generated Fortran coding can also be used to provide non-uniform emissivity and scattering coefficients.

The **differences** from the treatment of other properties are:

- no entries have been provided for emissivity and scattering
coefficients in the PROPS file, because it is usually small contaminants
ot the pure materials which affect the radiation, for example water
droplets (mist) in air, or soot particles in combustion gases; and
- no GX-style formulae, accessed via GRNDx settings, are provided, because the availability of InForm.

Examples will be found in the input-library cases:

r201.htm,
r202.htm,
r203.htm and
r204.htm.
However, it should be noted that, when EMIS is set for a cell containing a
non-transparent solid, it becomes dimensionless; and it then represents the
emissivity/absorptivity of the **surface** of the solid.

For example, it may be desired to represent an aperture in the bounding surface, though which radiation can enter and leave; or a domain boundary may be required to simulate a solid wall of prescribed temperature, for which the emissivity must also be defined.

To meet this requirement, PHOENICS has been caused to react
appropriately to PATCH and COVAL statements, for which the patch-name
begins with the characters 'IM' **and** the dependent variable in
question is T3.

How such patches are used is explained in the section of the lecture on IMMERSOL which may be viewed by clicking here,

**(a) PLATE**

For a specified wall temperature and emissivity, the radiative heat flux at the wall is computed from:

q

_{w}= σβ(T_{w}^{4}-T_{3,p}^{4}) = 4σβT_{3,m}^{3}(T_{w}-T_{3,p}) Eqn (22)

where the wall resistance term β is given by:

β = [δ(0.75(ε

^{'}+s^{'})+1/W_{gap}) + (1-ε_{w})/ε_{w}]^{-1}Eqn (23)

and T_{3,m} is a mean temperature between wall w and the near-wall node p, which is evaluated from:

T

_{3,m}^{3}=(T_{w}^{2}+T_{3,p}^{2})(T_{w}+T_{3,p})/4 Eqn (24)

**(b) Participating BLOCKAGE**

Likewise, at the solid-fluid interface of a participating BLOCKAGE, the radiative heat flux at the wall for a specified emissivity is computed from:

q

_{w}= σβ(T_{s}^{4}-T_{3,g}^{4}) = 4σβT_{3,m}^{3}(T_{s}-T_{3,g}) Eqn (25)

where the wall resistance term β is given by:

β = [δ(0.75(ε

^{'}+s^{'})+1/W_{gap}) + (1-ε_{w})/ε_{w}]^{-1}Eqn (26)

and T_{3,m} is a mean temperature between solid node and the fluid node g, which is evaluated from:

T

_{3,m}^{3}=(T_{s}^{2}+T_{3,g}^{2})(T_{s}+T_{3,g})/4 Eqn (27)

**(c) Non-participating BLOCKAGE**

In PHOENICS 2018 and earlier, no action is taken at the adiabatic surface so that both the convective and radiative heat flux are zero.

In PHOENICS 2019 and onwards, both the radiative and convective heat flux are evaluated at the surface, and surface temperature is computed from the requirement that the total heat transfer rate is zero.

**(d) THIN PLATE**

The boundary conditions applied on either side of a thin plate are explained in detail here

- The flow is steady and laminar.
- The Reynolds number is 10
- The pressure is atmospheric.
- The dimensions are small and the air is unpolluted, so that it does not participate at all in the radiative process.
- The emissivities of the surfaces are:-
- bottom surface... 0.9
- top surface.... 0.9
- first solid...... 0.1
- second solid... 0.2

The practically interesting question is: What temperatures will be attained in the solid?

First some images showing the geometry and some results:

The grid and velocity vectors

The distributions of WGAP and of PRPS

The true (TEM1) and radiation (T3) temperatures

Since it is not easily seen from the above contour plot that TEM1 and T3 are very close within the solid, the following plot is provided of the difference between them. This is easily calculated by placing the following In-Form line in the Q1 file.

(stored tdif is tem1-t3)

True temperature minus radiation temperature.

The upward (y) and rightward (x) radiative heat fluxes

IMMERSOL has allowed this problem to be set up and solved with great ease, the advantage of which can be best appreciated by those who have attempted to solve such problems by the methods advocated in many text-books, involving use of 'view-factor' calculations.

- Comparisons with text-book solutions (A.F.Mills):
- Other examples
- Use of IMMERSOL in connection with the calculation of stresses in solids
- Application to an electronics-cooling problem.

- IMMERSOL was conceived and introduced into PHOENICS during
1996.
- Wherever tested, it has performed well, in respect of:
- accuracy, where exact solutions are known;
- plausibility, where they are not;
- economy, in all circumstances.

- Since it is the ONLY model which can be economically used in complex circumstances, including multi-phase, chemically reactive and body-fitted coordinates, it seems likely to become the most widely used of all the radiation models in PHOENICS.

- Spalding, D.B: Trends, Tricks and Try-ons in CFD/CHT. In Advances in Heat Transfer, Vol.45, pp1-78, Editors: E.M.Sparrow, Y.I.Chou, J.P.Abraham, J.M.Gorman, Academic Press, Elsevier. (2013)