A probabilistic study of the influence of parameter uncertainty on thermal radiation heat transfer

International Journal of Heat and Mass Transfer 53 (2010) 1461–1472

Contents lists available at ScienceDirect

International Journal of Heat and Mass Transfer journal homepage: www.elsevier.com/locate/ijhmt

A probabilistic study of the influence of parameter uncertainty on thermal radiation heat transfer M.M.R. Williams Applied Modelling and Computational Group, Department of Earth Science and Engineering, Imperial College of Science, Technology and Medicine, Prince Consort Road, London SW7 2BP, UK

a r t i c l e

i n f o

Article history: Received 21 October 2009 Accepted 3 December 2009 Available online 6 January 2010 Keywords: Data uncertainty Probability Polynomial chaos

a b s t r a c t The influence of uncertainty in the boundary conditions on the solution of the radiative transfer equation is studied using probabilistic methods. We consider the problem of a heated slab, the surfaces and volume of which are at fixed temperatures. The basic problem stemmed from the need to study a hypothetical accident scenario in a fast breeder nuclear reactor in the cover gas region which may contain a sodium aerosol. The emissivities of the surface of such a region will be uncertain and estimates must be made. We demonstrate how the uncertainties in these data parameters are transmitted to the solution of the associated transfer equation and define sensitivity coefficients for the configuration factors. We also consider a number of different material surfaces, the emissivities of which are only known within certain limits, and obtain sensitivity coefficients for them as well as mean values, variances and, in some cases, the probability distribution of the configuration factors. A variational method is used to solve the basic equation, but for more general problems where such methods are inapplicable we describe a technique based on polynomial chaos theory. Ó 2009 Elsevier Ltd. All rights reserved.

1. Introduction An important source of error in calculations employing the radiative transfer (RT) equation arises from uncertain data values. These errors originate from lack of knowledge about the absorption and scattering coefficients and also of the emissivities present in the boundary conditions. It is vital, therefore, to know how changes in these parameters affect the outcome of the calculations. The problem of uncertainties in the absorption and scattering coefficients has already been examined [1] and so, in this paper, we wish to concentrate on the problem of boundary conditions where, for practical reasons, it is not possible to measure directly the emissivity of an industrial surface and only estimates within a certain range are available. In particular we are concerned with a problem that arose from a study of a hypothetical accident scenario in a fast breeder nuclear reactor where a breach in the pressure vessel may lead to the cover gas region being filled with a sodium aerosol which can absorb and emit radiation. This, together with the cover gas surfaces, leads to a radiant heat transfer problem which may be studied by the introduction of configuration factors. In this respect, it is important to know how the uncertainties in the data affect the heat transfer; namely, what are the uncertainties in the heat flux in relation to uncertainties in the emissivities? When the uncertainties are relatively small, these matters can be been handled by per-

E-mail address: [email protected] 0017-9310/$ - see front matter Ó 2009 Elsevier Ltd. All rights reserved. doi:10.1016/j.ijheatmasstransfer.2009.12.004

turbation theory [2,3]. However, there is another effective method of dealing with such a problem which does not have the limitations of perturbation theory and which may be used for arbitrarily large parameter uncertainties; namely, polynomial chaos expansions and in certain cases direct statistical averaging. In order to proceed, therefore, we deliberately introduce uncertainty into the boundary conditions of the radiative transport equation. This is done by assuming that the parameters are functions of a random variable which, in turn, allows us to prescribe a probability distribution for the values of the parameters. In the present case, a uniform distribution will be employed instead of the often used Gaussian. This is done for several reasons but the main one is that Gaussians can give rise to physically impossible data values, e.g. negative values, if the variance of the distribution is large. We avoid that pitfall by the present method which can bracket the uncertainty to within an upper and lower bound with an appropriate distribution in the interval. The method of attack will be to consider a simple geometry, i.e. plane parallel surfaces bounding a heated medium. The associated RT equation will be formulated and cast into a form for solution by the variational method [4]. However, because many problems can arise which are not conveniently cast into variational form, we also discuss a more general technique based upon an expansion of the unknown radiation intensity in terms of a series of polynomial chaos functions, more details of which will be discussed below. Essentially, the polynomial chaos expansion (PCE) [5–10] will be employed to transform the stochastic RT equation into a set of coupled deterministic transport equations for the

1462

M.M.R. Williams / International Journal of Heat and Mass Transfer 53 (2010) 1461–1472

Nomenclature I(r, X, m; n) radiation intensity as a function of position r, direction X, frequency m and random variables n . I0(s; n) scalar intensity at optical path s s = (ja + js)a optical path a physical width of system j(m0 ;n) = ja(m;n) + js(m;n) absorption plus scattering coefficient 0 f(m ? m, XX ;n) phase function for differential scattering w = js/(ja + js) albedo T physical temperature

coefficients in the polynomial chaos expansion (PCE). From these expansion coefficients, the mean and variance and higher statistical moments may be obtained. In particular, we shall calculate the mean and variance of the configuration factors and also form probability distributions for them by a simulation technique. Cases for a transparent medium will be given and also for a non-scattering medium for which exact results are available. However, by means of the variational method discussed above, we also obtain the magnitude of the influence of scattering on the configuration factors and the statistical uncertainty introduced by this. It is important to note that a number of other fields of technology have benefited from the use of PCE to study the effects of data uncertainty. For example, the uncertainty of modes of vibration in structural mechanics due to errors in elastic constants [11], the flow of groundwater in the geosphere due to rock fractures [12] and uncertainties in fluid flow due to both geometric uncertainty and basic data [13]. 2. The radiative transfer equation

r e q En(s) #

s Cn(. . .) Pi (n)

Stefan–Boltzmann constant emissivity heat flux exponential integral function of order n rT4/p optical path related to the polynomial chaos functions Legendre polynomials

Ið0; l; nÞ ¼ e1 ðnÞ

rT 41 þ 2ð1  e1 ðnÞÞ p

Z

1

dl0 l0 Ið0; l0 ; nÞ;

l>0

0

ð4Þ and

Iðs0 ; l; nÞ ¼ e2 ðnÞ

rT 42 þ 2ð1  e2 ðnÞÞ p

Z

1

dl0 l0 Iðs0 ; l0 ; nÞ;

l

0

<0

ð5Þ

where e1 (n) and e2 (n) are the emissivities of surfaces 1 and 2, respectively, and are regarded as uncertain by virtue of the random variables n. We note that each surface is at a different temperature (T1, T2) and that the volume in between is also at a different and spatially varying temperature T(s). Clearly the optical thickness of the system is s0 = (ja + js)a, where a is the actual physical thickness. We also define the heat flux across the system as

qðs; nÞ ¼

Z

1

dllIðs; l; nÞ

ð6Þ

1

The radiative transfer (RT) equation for the radiation intensity I(r, X, m; n), may be written as [14]

X  rIðr; X; m; nÞ þ jðm; nÞIðr; X; m; nÞ Z Z ¼ dm0 dX0 js ðm0 ; nÞ  f ðm0 ! m; X  X0 ; nÞIðr; X0 ; m0 ; nÞ þ Sðr; X; m; nÞ

ð1Þ

In Eq. (1) the extinction coefficient j(m; n) = ja(m; n) + js(m; n), where ja and js are the absorption and scattering coefficients at 0 0 frequency m, respectively. f(m ? m, XX ; n) is the energy and angle exchange phase function. S(r, X, m; n) is a source term which itself may contain uncertainties and is the Planck function and/or an independent source of photons. The n(n1, n2,. . . nN) are independent random variables which will describe the data uncertainty. In many cases we may consider the grey approximation to Eq. (1) in which it is averaged over frequency and appropriately averaged scattering and absorption coefficients used [14]. In this paper we are concerned only with problems involving parallel plates and hence Eq. (1) then becomes for isotropic scattering,

@Iðs; l; nÞ rT 4 ðsÞ w l þ Iðs; l; nÞ ¼ ð1  wÞ þ I0 ðs; nÞ @s 2 p

ð2Þ

where l is the cosine of the angle between the photon trajectory and the s axis, with s = (ja + js)x. The albedo w is defined as w = js/(js + ja). The scalar intensity is defined as

I0 ðs; nÞ ¼

Z

and we are particularly interested in q(0; n) and q(s0; n). At this point we refer the reader to our earlier work in which Eqs. (2), (4), and (5) are converted to integral form [15,16] and expressions for the configuration factors obtained. We summarise the results below. Defining # = rT4/2p, Ai ¼ ei rT 4i =p and Bi = 2(1  ei), we have for the scalar intensity

I0 ðs; nÞ ¼ gðs; nÞ þ

w 2

Z s0

ds0 Kðs; s0 ; nÞI0 ðs0 ; nÞ

ð7Þ

0

where

gðs; nÞ ¼ 4ð1  wÞ# þ 2e1 P½#1  ð1  wÞ#K 2 ðsÞ þ 2e2 P½#2  ð1  wÞ#K 1 ðsÞ

ð8Þ

Kðs; s0 ; nÞ ¼ E1 ðjs  s0 jÞ þ 2P½ð1  e1 ÞE2 ðsÞE2 ðs0 Þ þ ð1  e2 ÞE2 ðs0  sÞE2 ðs0  s0 Þ þ 4ð1  e1 Þ  ð1  e2 ÞPE3 ðs0 Þ½E2 ðsÞE2 ðs0  s0 Þ þ E2 ðs0 ÞE2 ðs0  sÞ ð9Þ with

P ¼ ½1  4ð1  e1 Þð1  e2 ÞE23 ðs0 Þ1

ð10Þ

K 1 ðsÞ ¼ E2 ðs0  sÞ þ B1 E3 ðs0 ÞE2 ðsÞ

ð11Þ

K 2 ðsÞ ¼ E2 ðsÞ þ B2 E3 ðs0 ÞE2 ðs0  sÞ

ð12Þ

0

We note that K(s, s ; n) is symmetric. The values for the heat fluxes at the surfaces are

1

dlIðs; l; nÞ

ð3Þ

qð0Þ ¼ e1 ð#1  Jð0ÞÞ

ð13Þ

1

The associated boundary conditions for each surface at s = 0 and s = s0 are

and

qðs0 Þ ¼ e2 ðJðs0 Þ  #2 Þ

ð14Þ

M.M.R. Williams / International Journal of Heat and Mass Transfer 53 (2010) 1461–1472

where

Jð0Þ ¼

A2 E3 ð 0 Þ þ B2 A1 E23 ð 0 Þ 1  B1 B2 E23 ð 0 Þ

s 1 þ s 1  B1 B2 E23 ðs0 Þ  ½E2 ðs0 Þ þ B2 E3 ðs0 ÞE2 ðs0  s0 Þ

Jðs0 Þ ¼

s

A1 E3 ðs0 Þ þ B1 A2 E23 ðs0 Þ 1  B1 B2 E23 ðs0 Þ

þ

Z s0

which is a standard result, but we remind the reader that the emissivities are now random variables due to their inherent uncertainty. For this simple case, DS1 = DS2 = 0 and

ds0 S0 ðs0 Þ

0

ð15Þ

1

Z s0

1  B1 B2 E23 ðs0 Þ

0

ds0 S0 ðs0 Þ

S1S2 ¼

1 1 þ e1 e2  1 1

ð16Þ

rT 4 w þ I ðsÞ p 2 0

ð17Þ

with

The functions En(z) are the exponential integral functions. Also, for notational simplicity, we have omitted the random variables n in the above equations. The stochastic equations derived above can be solved by Monte Carlo simulation and the surface heat fluxes obtained, however this is a time-consuming process and there are practical reasons for recasting the equations into a more conventional form as we will see below.

If the medium between the plates is non-scattering then we may by definition set w = 0. In that case we have an exact solution for the scalar intensity because the integral term in Eq. (7) is absent and we may write

I0 ðsÞ ¼ 4# þ 2e1 P½#1  #K 2 ðsÞ þ 2e2 P½#2  #K 1 ðsÞ Thus from Eqs. (13)–(17) the heat fluxes are

qð0Þ ¼ e1 P½1  4ð1  e2 ÞE23 ðs0 Þ#1  2e1 e2 PE3 ðs0 Þ#2  e1 P½1  2E3 ðs0 Þ½1 þ 2ð1  e2 ÞE3 ðs0 Þ#

where

ð26Þ

and

þ e2 P½1  2E3 ðs0 Þ½1 þ 2ð1  e1 ÞE3 ðs0 Þ# A configuration factor representation of the above problem has some value and the background theory may be found in Siegel and Howell [17]. In the present case there are three distinct configuration factors. The first, which we designate by S1S2, is the fraction of radiant heat leaving one surface that is incident on the other. In this case we can set #2 = # = 0 and #1 = 1 in Eq. (14) to get

ðsÞ I0 ð

ð25Þ

qðs0 Þ ¼ e2 P½1  4ð1  e1 ÞE23 ðs0 Þ#2 þ 2e1 e2 PE3 ðs0 Þ#1

3. Configuration factors

w S1S2 ¼ 2e1 e2 PE3 ðs0 Þ þ e2 P 2

ð24Þ

3.2. A non-scattering medium

 ½E2 ðs0  s0 Þ þ B1 E3 ðs0 ÞE2 ðs0 Þ

S0 ðsÞ ¼ ð1  wÞ

1463

Z s0 0

ds

0 ðsÞ I0 ð 0 ÞK 1 ð 0 Þ

s

s

ð18Þ

sÞ is the solution of the integral equation

Z w s0 0 ðsÞ ðsÞ I0 ðsÞ ¼ 2e1 PK 2 ðsÞ þ ds Kðs; s0 ÞI0 ðs0 Þ 2 0

S1S2 ¼ 2e1 e2 PE3 ðs0 Þ

ð28Þ

DS1 ¼ e1 P½1 þ 2ð1  e2 ÞE3 ðs0 Þ½1  2E3 ðs0 Þ

ð29Þ

DS2 ¼ e2 P½1 þ 2ð1  e1 ÞE3 ðs0 Þ½1  2E3 ðs0 Þ

ð30Þ

and that consequently we may write the heat fluxes as

ð20Þ

ð31Þ

and

qðs0 Þ ¼ DS2ð#  #2 Þ þ S1S2ð#1  #2 Þ

The second configuration factor is the fraction of radiant heat emanating from the volume heat source defined by T, which is incident on a surface at T1 or T2. In this case we have DS1 and DS2, respectively. These are defined by setting #1 = #2 = 0 and # = 1 and noting that DS1 = q(0) and DS2 = q(s0). In that case we have

DS1 ¼ e1 ð1  wÞP½1 þ 2ð1  e2 ÞE3 ðs0 Þ½1  2E3 ðs0 Þ Z w s0 0 ðdÞ 0 þ e1 P ds I0 ðs ÞK 2 ðs0 Þ 2 0

It is clear from Eqs(18), (20), and (21) that the configuration factors are

qð0Þ ¼ DS1ð#1  #Þ þ S1S2ð#1  #2 Þ ð19Þ

ð27Þ

ð32Þ

In non-scattering media, the advantage of configuration factors is that one may write heat exchange between surfaces very easily in terms of them [17]. It should be noted that for the special cases derived above, we have explicit expressions for the configuration factors. Thus if the random form of the emissivities is known, a straightforward quadrature will give the mean and variance and indeed any higher statistical moment. Details of this averaging procedure will be discussed in detail below.

and

DS2 ¼ e2 ð1  wÞP½1 þ 2ð1  e1 ÞE3 ðs0 Þ½1  2E3 ðs0 Þ Z w s0 0 ðdÞ 0 þ e2 P ds I0 ðs ÞK 1 ðs0 Þ 2 0

3.3. Scattering medium: variational method

ð21Þ

ðdÞ

where now I0 ðsÞ satisfies the integral equation ðdÞ

I0 ðsÞ ¼ 4ð1  wÞ  2Pð1  wÞ½e1 K 2 ðsÞ þ e2 K 1 ðsÞ Z w s0 0 ðdÞ þ ds Kðs; s0 ÞI0 ðs0 Þ 2 0

ð22Þ

It is useful to consider some special limiting cases of the general solution. 3.1. Transparent medium In the limit of zero optical path, i.e. where s0 ? 0, we find

#1  #2 qð0Þ ¼ qðs0 Þ ¼ 1 1 e1 þ e2  1

For the general case of a scattering medium, it is necessary to solve the integral equations and then insert the solutions into the expressions for S1S2, DS1 and DS2. We shall consider this matter below, but first we note that in our earlier work on this subject [15], we employed a variational method to obtain explicit expressions for the configuration factors. This involved forming a stationary functional for these quantities and using a trial function to obtain lower bounds for them [4]. This procedure is very accurate and the detailed results may be found in ref [15], but the general expressions are given below. Thus, once again, we are able to perform our statistical averages directly. From the variational method we find for the configuration factors

S1S2 ¼ 2e1 e2 PE3 ðs0 Þ þ we1 e2 P2 F st1 ðe1 ; e2 Þ ð23Þ

ð33Þ

where the stationary functional Fst1 with the appropriate trial functions may be found in [15]. For the other case,

1464

M.M.R. Williams / International Journal of Heat and Mass Transfer 53 (2010) 1461–1472

DS1 ¼ e1 P½1 þ 2ð1  e2 ÞE3 ðs0 Þ½1  2E3 ðs0 Þ 2

 e1 P wF st2 ðe1 ; e2 Þ

ð34Þ

with the functional Fst2 also available from [15]. An analogous result arises for DS2. It should be noted that the parameters in S1S2, DS1 and DS2 are all functions of e1(n1) and e2(n2) where n1 and n2 are independent random variables, the properties of which will be discussed below. In order to obtain the mean and variance of the configuration factors CF, it is necessary to calculate the averages

hCF n i ¼

Z

dnCF ðnÞn

R

where the mean is hCFi and the variance hCF 2 i  hCFi2 . Values of these quantities will be obtained below. 4. Polynomial chaos expansions

ðN þ ‘Þ! N!‘!

ð37Þ

Thus when N = 2 as we have seen, each U‘ corresponds to a combination (‘1, ‘2). Clearly as N increases, the number of terms in the expansion increases rapidly. This will place some practical limitations on the efficiency of the numerical work. Let us now apply this technique to the calculation of S1S2 using Eqs. (18) and (19). Before doing so, however, we write the emissivities as follows:

1 2

1 2

ei ðni Þ ¼ ðei0 þ ei1 Þ þ ðei1  ei0 Þni

A method for solving stochastic equations in general was introduced by Wiener [7] who called it homogeneous chaos. In one sense this involves a Fourier expansion with the expansion polynomials being an orthogonal set of functions of random variables rather than deterministic ones. Wiener’s ideas have been extended over the years and it is found that the solution of the integral Eq. (7) can be written

hei i ¼

1 ðei0 þ ei1 Þ 2

ð39Þ

1 ðei1  ei0 Þ2 12

ð40Þ

and

r2ei ¼

4.2. The integral equation

ðsÞ

I0 ðs; n1 ; n2 Þ ¼

N X

U‘ ðn1 ; n2 Þ ¼ P‘1 ðn1 ÞP‘2 ðn2 Þ

ðsÞ

N X

/ðsÞ n ðsÞUn ðn1 ; n2 Þ ¼ 2e1 ðn1 ÞPðn1 ; n2 ÞK 2 ðs; n1 ; n2 Þ

ð35Þ

ð36Þ

ð41Þ

where we have re-instated the random variable dependence in I0 . Inserting (41) into (19), we find

Z w s0 0 ds Kðs; s0 ; n1 ; n2 Þ 2 0 N X 0  /ðsÞ n ðs ÞUn ðn1 ; n2 Þ

þ

n¼0

where the polynomial chaos functions Un(n) are defined in [5] and the polynomial expansion coefficients (PCE), un(s), are deterministic. Thus the randomness is moved to the polynomials Un(n). Although Wiener, in his original paper, only dealt with Gaussian random variables, it has been shown that the same ideas apply to a large number of orthogonal polynomials, e.g. those of Laguerre and Jacobi. These generalised polynomials have been shown to satisfy orthogonality relations analogous to their deterministic counterpart, with an appropriate weight function [18] corresponding to the probability distribution function of the random variables. In general, as many orders of nn will be present in the expansion (35) as there are random variables in the problem. As an explicit example, let us assume that the problem has only two independent random variables. In that case, the polynomial chaos expansion uses the functions U‘(n1, n2). At this point it is necessary to specify the probability distribution of the random variables, which in turn depends on the nature of the problem. For the case of random emissivities, as discussed in this paper, we will assume that the uncertainty can be described by upper and lower limits, with a uniform distribution in the interval. Thus we assume that n1 and n2 are uniformly distributed in the range (1, 1), and the appropriate polynomials are those of Legendre with a weight function p(n) = 1/2 [19,20] such that

/ðsÞ n ðsÞUn ðn1 ; n2 Þ

n¼0

n¼0

/n ðsÞUn ðnÞ

ð38Þ

Let us now expand the radiation intensity as in Eq. (35) for the S1S2 equation, viz,

4.1. Polynomial chaos

I0 ðs; nÞ ¼

P¼

where ei0 6 ei 6 ei1 . Thus the mean value and variance of the ei’s are

As we have seen above, if it is possible to obtain explicit results for the configuration factors then direct statistical averaging may be carried out easily. However, in many cases, the integral equation and the quantity of interest may not be conveniently expressed in terms of a variational principle and so a numerical solution of the integral equation is required. To reduce the stochastic integral equation to a form convenient for numerical evaluation it is recommended that the concept of polynomial chaos be introduced, which we briefly outline in Section 4.1.

1 X

where ‘ = ‘1 + ‘2 and P‘(n) are the Legendre polynomials. For a given value of ‘ we take all combinations of ‘1 and ‘2 that satisfy ‘ = ‘1 + ‘2. It may be shown that if N is the total number of independent variables ni, then the number of terms in the PCE is

ð42Þ

n¼0

Now multiplying Eq. (42) by Um(n1, n2) and integrating over n1 and n2, we find

hU2m ðn1 ; n2 Þi/ðsÞ m ðsÞ ¼ h2e1 ðn1 ÞPðn1 ; n2 ÞK 2 ðs; n1 ; n2 ÞUm ðn1 ; n2 Þi N Z s0 wX 0 ds0 /ðsÞ þ n ðs Þ 2 n¼0 0  hUm ðn1 ; n2 ÞKðs; s0 ; n1 ; n2 ÞUn ðn1 ; n2 Þi

ð43Þ

where m = 1, 2,. . .N and the angular bracket is defined by the double integral hUm ðn1 ; n2 ÞBðn1 ; n2 ÞUn ðn1 ; n2 Þi 

1 4 

Z

1

dn1

1

Z

1

dn2 Um ðn1 ; n2 ÞBðn1 ; n2 ÞUn ðn1 ; n2 Þ

ð44Þ

1

We also define N2m ¼ hU2m i. It is clear that on re-arrangement, Eq. (43) reduces to the symmetric form

wm ðsÞ ¼ fm ðsÞ þ

N Z s0 wX ds0 wn ðs0 ÞK m;n ðs; s0 Þ 2 n¼0 0

ð45Þ

1465

M.M.R. Williams / International Journal of Heat and Mass Transfer 53 (2010) 1461–1472

with wm ¼ N m /ðsÞ m and

1 h2e1 ðn1 ÞPðn1 ; n2 ÞK 2 ðs; n1 ; n2 ÞUm ðn1 ; n2 Þi Nm 1 hUm ðn1 ; n2 ÞKðs; s0 ; n1 ; n2 ÞUn ðn1 ; n2 Þi K m;n ðs; s0 Þ ¼ Nn Nm

fm ðsÞ ¼

l ð46Þ ð47Þ

Thus we have reduced the stochastic integral Eq. (19) to a set of coupled deterministic integral equations for the polynomial chaos expansion amplitudes. With these equations, one goes back to the definition of S1S2, Eq. (18), to get

S1S2ðn1 ; n2 Þ ¼ 2e1 ðn1 Þe2 ðn2 ÞPðn1 ; n2 ÞE3 ðs0 Þ þ e2 ðn2 ÞPðn1 ; n2 Þ Z s0 N wX 0 0  Un ðn1 ; n2 Þ ds0 /ðsÞ n ðs ÞK 1 ðs ; n1 ; n2 Þ 2 n¼0 0

/m ð0; lÞ ¼ a1;m

N X rT 41 þ2 b1;m;n p n¼0

Z

1

dl0 l0 /n ð0; l0 Þ;

l

0

>0

ð55Þ

and

ð48Þ

hS1S2i ¼ 2he1 ðn1 Þe2 ðn2 ÞPðn1 ; n2 ÞiE3 ðs0 Þ N Z s0 wX 0 ds0 /ðsÞ þ n ðs Þ 2 n¼0 0

ð54Þ

with the boundary conditions (4) and (5) becoming

/m ðs0 ; lÞ ¼ a2;m

To find the mean and variance of S1S2, we then statistically average as follows

N X rT 42 þ2 b2;m;n p n¼0

Z

1

dl0 l0 /n ðs0 ; l0 Þ;

l

0

<0

ð56Þ

where

ai;m ¼

hei Um i 2 mi

hU

and bi;m;n ¼

hUn ð1  ei ÞUm i hU2m i

ð57Þ

The corresponding heat flux is given by

0

 hUn ðn1 ; n2 Þe2 ðn2 ÞPðn1 ; n2 ÞK 1 ðs ; n1 ; n2 Þi

ð49Þ

and

qðs; nÞ ¼

N Z X n¼0

1

dll/n ðs; lÞUn ðnÞ

ð58Þ

1

There are many methods available to solve the set of Eqs. (54)–(56) including well-established computer codes for radiation and neutron transport [22].

hS1S22 i ¼ 4he1 ðn1 Þ2 e2 ðn2 Þ2 Pðn1 ;n2 Þ2 iE23 ðs0 Þ N Z s0 X 0 ds0 /ðsÞ þ 2wE3 ðs0 Þ n ðs ÞhUn ðn1 ; n2 Þe1 ðn1 Þ n¼0

@/n ðs; lÞ rT 4 w þ /n ðs; lÞ ¼ ð1  wÞ d þ / ðsÞ @s p n;0 2 0n

5. Numerical results and discussion

0

2 2

 e2 ðn2 Þ P ðn1 ;n2 ÞK 1 ðs0 ; n1 ;n2 Þi N X N Z s0 w2 X 0 ds0 /ðsÞ þ n ðs Þ 4 n¼0 m¼0 0 * + Z s0 e2 ðn2 Þ2 P2 ðn1 ; n2 ÞK 1 ðs0 ;n1 ; n2 Þ: 00 ð50Þ ds00 /ðsÞ ð s Þ  m K 1 ðs00 ;n1 ; n2 ÞUn ðn1 ;n2 ÞUm ðn1 ;n2 Þ 0 The known values of /ðsÞ n ðsÞ are then inserted into (49) and (50). Although this may appear to be a numerically tedious process, the polynomial chaos expansion has the powerful virtue of avoiding any of the closure approximations which are usually encountered when solving stochastic equations [21]. An analogous procedure can be carried out to obtain DS1 and DS2 from Eqs. (20)–(22). In this paper it is not necessary to evaluate these integral equations but, to illustrate the rate of convergence of the series expansion (41), we shall expand the intensity for the w = 0 (non-scattering) case in terms of polynomial chaos functions and this will be described in Section 4.3. Finally, we note from the properties of the polynomial chaos functions that the mean intensity is given by ðsÞ

ðsÞ

hI0 ðs; nÞi ¼ /0 ðsÞ

ð51Þ

and its variance by

r2IðsÞ ðsÞ ¼

N X 2 2 ½/ðsÞ n ðsÞ hUn i

ð52Þ

n¼1

4.3. The integro-differential equation Instead of using the integral form for the scalar intensity, it is sometimes more convenient to apply PCE to the integro-differential form of the transport equation as given by Eq. (2) and its associated boundary conditions, Eqs. (4) and (5). Thus if we write

I0 ðs; l; n1 ; n2 Þ ¼

N X

/n ðs; lÞUn ðn1 ; n2 Þ

ð53Þ

n¼0

insert it into Eq. (2), multiply by Um(n) and statistically average, we find that the PCE terms are given by

This section will be concerned with numerical results for three different physical situations: (a) the transparent medium, where there is no material between the surfaces, (b) a purely absorbing medium, i.e. non-scattering, for which exact results are available, (c) an assessment of the influence of scattering on the configuration factors. Finally, we assess the convergence rate of the polynomial chaos expansion for the non-scattering case. In all instances, for the numerical work, we use the IMSL Library of subroutines [23]. The data for the emissivities used in our calculations is taken form the CRC Handbook [24] and for later use we give below in Table 1 the upper and lower values of the emissivities cited in [24]. 5.1. The transparent medium To obtain a clear understanding of the influence of uncertainty, we first consider the very simple case as given by Eq. (24). Table 2 gives values for the configuration factor S1S2 and its statistical variation. We stress here the importance of the relative standard deviation as this is a direct measure of the effect of uncertainty. For example, Relstd multiplied by 100 is the percentage deviation from

Table 1 Upper and lower bounds on the emissivities of various materials [24]. Material

elower

eupper

Al Be Ce Cr Cu Fe Mg Ni Th Sn Zr Porcelain

0.22 0.07 0.58 0.60 0.60 0.63 0.10 0.85 0.2 0.32 0.18 0.25

0.40 0.37 0.80 0.80 0.80 0.98 0.43 0.96 0.57 0.60 0.43 0.50

1466

M.M.R. Williams / International Journal of Heat and Mass Transfer 53 (2010) 1461–1472

Table 2 S1S2 and its statistical properties. Material

hS1S2i

MP

R

Lower

Upper

Relstd

S

Al Be Ce Cr Cu Fe Mg Ni Th Sn Zr Porcelain

0.1809 0.1141 0.5251 0.5372 0.5372 0.6709 0.1431 0.8263 0.2299 0.2945 0.1751 0.2268

0.167 0.065 0.508 0.523 0.523 0.623 0.091 0.821 0.175 0.265 0.147 0.201

0.986 0.923 0.997 0.998 0.998 0.996 0.937 0.999 0.964 0.986 0.973 0.983

0.124 0.0363 0.408 0.428 0.428 0.461 0.0530 0.739 0.111 0.190 0.0989 0.142

0.250 0.227 0.667 0.667 0.667 0.962 0.274 0.923 0.399 0.429 0.274 0.333

0.1445 0.3669 0.1004 0.0905 0.0905 0.1519 0.3346 0.0454 0.2633 0.1669 0.2093 0.1742

0.8620 0.9321 1.091 1.097 1.097 1.210 0.9308 1.294 0.9489 0.9497 0.8846 0.9052

hS1S2i = statistical average of S1S2. R = hS1S2i/deterministic* value of S1S2. MP = most probable value of S1S2. Lower = lower bound of S1S2. Upper = upper bound of S1S2. Relstd = relative standard deviation of S1S2 = rðS1S2Þ=hS1S2i. S = sensitivity = Relstd(S1S2)/Relstd(e). *The deterministic value is obtained by using the average values of the emissivities, i.e. (elower + eupper)/2, directly in S1S2, DS1 or DS2.

the mean and we see from Table 2 that this varies between 9% and 33% according to the material; these are not insignificant values and show the importance of our statistical study in assessing the confidence that may be placed in a result. In Table 3 we have studied the case when each surface is of a different material, e.g. surface 1 is Al and surface 2 is Fe.

Table 3 Configuration factor for aluminium and iron surfaces. 0.2868 0.995 0.320 0.195 0.397 0.1635 0.9753 1.3026

hS1S2i R MP Lower Upper Relstd S1 S2

In Fig. 1 we have the associated probability distribution functions (pdf) for Al, Be, Cu, and Fe. Also we give the pdf for the combined Al–Fe case. The shapes of the pdfs are interesting. We note that, in the case of identical surfaces, the distribution is always quasi-triangular with the peak giving the most probable value of S1S2. On the other hand, if the surfaces have different properties, then the pdf is quasi-quadrilateral. This radical change of shape is due to the superposition effect of two different pdfs. In a more complex geometry where there are more than two surfaces, we would expect the pdf to show even more fine structure. A parameter of some practical interest is the sensitivity S as this gives a direct correlation between data values and S1S2. In most cases in Table 2, S is around unity but there are two materials, viz: Fe and Ni where the sensitivity is 20% larger and so, from this, we may deduce that uncertainties in the data give even larger uncertainties in the configuration factor. Another quantity of interest is the ratio R, which is a measure of the difference between the true statistical average and the deterministic value. We observe that R is always less than unity and this shows that the deterministic approximation is always an overestimate of S1S2, but not by a substantial amount; about 8% at the most. Note though that the actual spread of values of S1S2 can be large, as seen from the ‘upper’ and ‘lower’ values. Thus, in practice, if products were being produced in large numbers, we might expect a large range of uncertainty in an ensemble member; particularly if the surfaces are of differing materials. 5.2. The non-scattering medium Eqs. (28)–(30) denote the configuration factors for a non-scattering, i.e. w = 0, medium. To illustrate the uncertainties that arise for this case we consider only two materials; Al and Fe with an optical path between the surfaces of s0 = 1. Three problems are discussed: (a) both surfaces are Al, (b) both surfaces are Fe and (c) one surface is Al and the other Fe. In Table 4, we give the corresponding statistical values of the configuration factors; namely, their mean and relative standard deviation. In the Al–Fe case, surface 1 is Al and surface 2 is Fe. For hS1S2i, we note that the Al and Fe values for the mean bracket the Al-Fe value. On the other hand, for the hDS1i and hDS2i values we find the interesting result that the Al-Fe value of

20 ε1=0.40 ε0=0.22

TRANSPARENT MEDIUM

Al

15 Be

Cu

ε1=0.37

δ1=0.8

ε0=0.07

ε0=0.6

Al-Fe

10

ε1=0.40, 0.98 ε0=0.22, 0.63

Fe 5

ε1=0.98 ε0=0.63

0 0

0.2

0.4

0.6

0.8

1.0

S1S2 Fig. 1. Probability distributions for configuration factors S1S2 for aluminium, beryllium, copper and iron.

1467

M.M.R. Williams / International Journal of Heat and Mass Transfer 53 (2010) 1461–1472

values of S1S2 for the Al–Fe case lies between those for Al and Fe alone. However, the curious feature discussed above regarding DS1 and DS2 shows up well in the figure. Namely, we see that the pdf of DS1 for Al–Fe is shifted slightly below that for Al alone. Similarly the pdf for DS2 is shifted slightly above that for Fe alone. These two distributions are shown in more detail in Fig. 2b and 2c. In Table 5 we give the corresponding sensitivity factors for the three cases in Table 4. It is clear that S1S2 is most sensitive to uncertainty in the emissivity, there being an amplification of the uncertainty in S1S2 by up to about 70% in one case. For DS1 and DS2, the sensitivities are around unity when both surfaces are the same. However in the Al–Fe case, we note that there is amplification in one case and diminishment in the other. The conclusion may be drawn that the effect of uncertainty on the quantity of interest is not always directly proportional to the driving uncertainty.

Table 4 Configuration factors and their statistical moments for a non-scattering medium. Material

hS1S2i

hDS1i

hDS2i

Rel(S1S2)

Rel(DS1)

Rel(DS2)

Al Fe Al–Fe

0.02156 0.1424 0.05508

0.2851 0.6564 0.2540

0.2851 0.6564 0.7279

0.2363 0.1770 0.2081

0.1663 0.1265 0.1690

0.1663 0.1265 0.1226

NB: RelðXÞ  rðXÞ=hXi.

hDS1i is less than that for Al alone. Similarly the hDS2i value for Al– Fe is larger than that for Fe alone. This feature is also found for deterministic values so it does not arise from any statistical argument. Indeed, using Eqs. (28) and (29), this inequality may be shown to be a general result. That is, if e2 > e1, then DS1(e1, e1) > DS1(e1, e2) and DS2(e2, e2) < DS2(e1, e2). The results above are better viewed in terms of probability distribution functions (pdfs) and Fig. 2a shows P(S1S2), P(DS1) and P(DS2). Here it is clear how the

100 ABSORBING MEDIUM

S1S2 Al 80

w=0 τ0=1.0

60

DS1 Al-Fe DS2 Al-Fe S1S2 Al-Fe DS1 Fe S1S2 Fe DS1 Al S1S2 Al

S1S2 Al-Fe

40

S1S2 Fe

20

DS1 Al-Fe

DS1 Al

DS2 Fe

DS2 Al-Fe

0 0

0.2

0.4

0.6

0.8

DS1 or S1S2 Fig. 2a. Probability distributions for configuration factors S1S2, DS1, DS2, Al, Fe, Al–Fe.

10 ABSORBING MEDIUM 8

DS1 Al-Fe 0.128787

6

DS1 Al

DS1 Al-Fe DS1 Al

4

2

0 0.15

0.20

0.25

0.30

DS1 Fig. 2b. Probability distributions for configuration factors DS1, Al, Al–Fe.

0.35

1468

M.M.R. Williams / International Journal of Heat and Mass Transfer 53 (2010) 1461–1472

5

ABSORBING MEDIUM 4

3

DS2 Fe

DS2 Al-Fe

2

DS2 Al-Fe DS1 Fe

1

0 0.45

0.50

0.55

0.60

0.65

0.70

0.75

0.80

DS2 Fig. 2c. Probability distributions for configuration factors DS2, Fe, Al–Fe.

Table 5 Sensitivities for the configuration factors.

Table 6 Mean values and relative standard deviations for an Al–Fe system.

Material

S1S2

DS1

DS2

S1

S2

S1

S2

S1

S2

Al Fe Al–Fe

1.410 1.410 1.242

1.410 1.410 1.658

0.989 0.993 0.997

0.989 0.993 1.332

0.989 0.993 0.729

0.989 0.993 0.974

w

hS1S2i

hDS1i

hDS2i

std(S1S2)

std(DS1)

std(DS2)

0.0 0.5 0.9 0.95 0.99 0.999

0.0551 0.0873 0.1735 0.1988 0.2277 0.2383

0.2436 0.2060 0.0832 0.0488 0.0127 0.00205

0.6322 0.5778 0.2241 0.1268 0.0271 0.00208

1.146E2 1.672E2 2.703E2 2.938E2 3.198E2 3.304E2

4.072E2 3.366E2 1.262E2 7.021E3 1.421E3 1.336E4

7.725E2 6.105E2 1.661E2 8.715E3 1.813E3 1.609E4

5.3. The scattering medium We now consider the case of a scattering medium, using the variational results given above. To illustrate how the uncertainty affects the configuration factors we assume that s0 = 1 and compare values of w = 0, 0.5, 0.95, 0.99 and 0.999. The latter case

means a very weakly absorbing medium. Table 6 shows the mean values and relative standard deviations for these cases for one surface of Al and the other of Fe. The values in Table 6 show, as expected from the deterministic case, that S1S2 increases with decreasing absorption and DS1 and

40 SCATTERING MEDIUM

τ0=1.0

w=0 30

w=0.999 w=0.99 w=0.9 w=0.5 w=0

w=0.5 20

w=0.9 w=0.99 w=0.999

10

0 0

0.05

0.10

0.15

0.20

0.25

0.30

0.35

S1S2 Fig. 3. Probability distribution function for configuration factor P(S1S2) Al–Fe surfaces.

1469

M.M.R. Williams / International Journal of Heat and Mass Transfer 53 (2010) 1461–1472

w=0.99 300

30 w=0.99

200

20

100

10

w=0.9

SCATTERING MEDIUM

w=0.99 w=0.9 w=0.5 w=0

w=0.5 w=0

0

0 0

0.1

0.2

0.3

0.4

DS1 Fig. 4. Probability distribution function for configuration factor P(DS1): Al–Fe surfaces.

w=0.99 250 25 w=0.99 200

20

150

15

100

10

50

5

0

0

w=0.9

SCATTERING MEDIUM w=0.99 w=0.9 w=0.5 w=0

w=0.5 w=0

0

0.2

0.4

0.6

0.8

DS2 Fig. 5. Probability distribution function for configuration factor P(DS2): Al–Fe surfaces.

DS2 decrease. The interesting conclusion arises from the variation of the standard deviation (std) with absorption. Here we find for S1S2 that the std increases as absorption decreases, but in the case of DS1 and DS2 the std decreases. Physically, this means that the effects of uncertainties in the boundary conditions becomes more influential for S1S2 as the medium becomes more scattering, but less influential for DS1 and DS2. Thus for DS1 and DS2 the diffusive nature of scattering reduces the proportion of direct flights between surfaces and the net effect is due to volume rather than surface uncertainty. The reverse is the case for S1S2. This phenomenon can be seen more clearly in Figs. 3–5. In Fig. 3 we show the pdf P(S1S2) for the range of w values and note that as w increases, the pdf becomes wider and flatter. On the other hand, in Figs. 4 and 5 for P(DS1) and P(DS2) the pdf becomes narrower and more peaked as well as moving to smaller values. This effect

Table 7 Polynomial chaos expansion coefficients for S1S2. n

/nðsÞ ðs Þ

0 1 2 3 4 5 6 7 8 9 10 11 12

0.1574 1.3147E002 0.1284 1.0965E004 1.0341E002 2.8670E004 8.3647E007 6.4818E005 2.8968E004 7.4229E007 6.1769E009 3.2589E007 1.6849E005

1470

M.M.R. Williams / International Journal of Heat and Mass Transfer 53 (2010) 1461–1472

0.9

0.8

Average/deterministic

0.7 w=0 Average/deterministic Deterministic S1S2 intensity Average S1S2 intensity

0.6

0.5

Deterministic S1S2 intensity

0.4

0.3

0.2 Average S1S2 intensity 0.1 0

0.1

0.2

0.3

0.4

0.5

0.6

τ

0.7

0.8

0.9

1.0

Fig. 6. Average and deterministic S1S2 flux across system: Al–Fe surfaces.

2

σI 0.05

0.52

0.04

0.51

w=0 Relative standard deviation 2 Variance σ1

0.03

0.50

0.02

0.49 Variance

0.01

0.48 Relative standard deviation

0

0.47 0

0.1

0.2

0.3

0.4

0.5

τ

0.6

0.7

0.8

0.9

1.0

Fig. 7. Variance and relative standard deviation of S1S2 flux, Al–Fe surfaces.

is due to the fact that as scattering increases, the associated angular distribution of photons becomes more isotropic and hence the uncertainty increases. The converse of zero scattering, i.e. a purely absorbing medium, leads to a highly singular form of photon angular distribution of the form d(l  l*) because there is no mechanism for a change of direction. Thus for highly scattering media the errors caused by uncertainties in emissivities is less important for DS1 and DS2 and more important for S1S2. 5.4. Polynomial chaos for a non-scattering medium As explained in Section 4, we now examine the rate of convergence of the polynomial chaos expansion by using the exact solu-

tion available when w=0. In that case we write the intensity from Eq. (19) as ðsÞ

I 0 ð sÞ ¼

N X

/ðsÞ n ðsÞUn ðn1 ; n2 Þ

ð59Þ

n¼0

whence the polynomial expansion coefficients are

/ðsÞ n ðsÞ ¼

1 hU2n i

h2e1 PK 2 ðsÞUn i

ð60Þ

Now from the DS Eqs. (20) - (22) we find by analogy with (59) and (60) that

1471

M.M.R. Williams / International Journal of Heat and Mass Transfer 53 (2010) 1461–1472

ðdÞ

I0 ðsÞ ¼

N X

/ðdÞ n ðsÞUn ðn1 ; n2 Þ

iance and higher statistical moments. Using these values the mean and variance can be obtained from Eqs. (51) and (52). Fig. 6 shows the mean intensity across the system and it is compared with the value that would have been obtained if the average values of the emissivities had been inserted directly into the intensity equation (the deterministic approximation). It is clear that significant errors arise if the correct statistical procedure is not used. According to Fig. 6, the statistically averaged mean intensity is some 70–80% of the deterministic value. Fig. 7 shows the variance and the relative standard deviation (relstd). It is clear that the relstd is about 50% which indicates large fluctuations about the mean. The values of the PCE coefficients for the DS1–DS2 equations are similar to those in Table 7 for S1S2. Figs. 8 and 9 show the DS1–DS2 intensities. Firstly we note that the difference between the statistical mean and the deterministic case is less for DS1 and DS2; the

ð61Þ

n¼0

Where

/ðsÞ n ð sÞ ¼

1 hU2n i

h½4  2Pðe1 K 2 ðsÞ þ e2 K 1 ðsÞÞUn i

ð62Þ

To illustrate these results, we use the data for Al–Fe and in Table 7 *  find for s0 = 1.0, the values for /ðsÞ n ðs Þ where s = s0/2. The polynomial chaos expansion coefficients are seen to decrease rapidly although not uniformly. Five terms in the series are usually adequate for acceptable accuracy. The convergence for other values of s* are similar to those in Table 7. We also note that the coefficients can be of either sign and have no particular physical significance except when combined to find the mean, var-

4

Average DS1-DS2 intensity

3

Deterministic DS1-DS2 intensity

w=0 Average DS1-DS2 intensity Deterministic DS1-DS2 intensity average/deterministic 2

average/deterministic 1

0

0.1

0.2

0.3

0.4

0.5

τ

0.6

0.7

0.8

0.9

1.0

Fig. 8. Average and deterministic DS1 intensity across system: Al–Fe surfaces.

0.20

0.15 w=0 Relative standard deviation variance 0.10

Relative standard deviation 0.05 variance

0

0.1

0.2

0.3

0.4

0.5

τ

0.6

0.7

0.8

0.9

Fig. 9. Variance and relative standard deviation of DS1 intensity: Al–Fe surfaces.

1.0

1472

M.M.R. Williams / International Journal of Heat and Mass Transfer 53 (2010) 1461–1472

average being around 8% greater than the deterministic intensity. The variance is also significantly less than the corresponding quantity for S1S2, i.e. 5–15% compared with 50%. The values reported here indicate that when Eq. (45) or Eqs. (54)–(56) are evaluated for the general case, only a relatively few terms will be needed; rarely more than 5. This is the importance of Table 7 as it enables an estimate of the work necessary for the more general case to be obtained in advance. It should also be noted that whilst we have used a uniform pdf for the emissivities, other probability distributions could be used, e.g. log-uniform, normal or log-normal should the situation merit it. References [1] M.M.R. Williams, M. Eaton, A probabilistic study of the influence of parameter uncertainty on solutions of the radiative transfer equation, J. Quant. Spectrosc. Rad. Trans., in press. [2] M.L. Williams, Perturbation theory for nuclear reactor analysis, in: Y. Ronen (Ed.), Handbook of Nuclear Reactor Calculations, vol. III, CRC Press, 1986. [3] D.G. Cacuci, Sensitivity and Uncertainty Analysis, Chapman and Hall, 2003. [4] M.M.R. Williams, Mathematical Methods in Particle Transport Theory, Butterworth, 1971. [5] R.G. Ghanem, P.D. Spanos, Stochastic Finite Elements, Dover Publications, 1991. [6] A.F. Emery, Some thoughts on solving the radiative transfer equation in stochastic media using polynomial chaos and Wick products as applied to radiative equilibrium, J. Quant. Spectrosc. Rad. Trans. 93 (2005) 61. [7] N. Wiener, The homogeneous chaos, Am. J. Math. 60 (1938) 897.

[8] S.-K. Choi, R.V. Grandhi, R.A. Canfield, Structural reliability under non-Gaussian stochastic behaviour, Comput. Struct. 82 (2004) 1113. [9] D. Xiu, G.E. Karniadakis, Modelling uncertainty in steady state diffusion problems via generalized polynomial chaos, Comput. Methods Appl. Mech. Eng. 191 (2002) 4927. [10] S.K. Sachdeva, P.B. Nair, A.J. Keane, Hybridization of stochastic reduced basis methods with polynomial chaos expansions, Prob. Eng. Mech. 21 (2006) 182. [11] R. Li, R.G. Ghanem, Adaptive polynomial chaos expansions applied to statistics of extremes in non-linear random vibration, Prob. Eng. Mech. 13 (1998) 125. [12] R.G. Ghanem, Scales of fluctuation and the propagation of uncertainty in random porous media, Water Res. Res. 34 (1998) 2123. [13] O.P. Le Maitre, O.M. Knio, H.N. Najm, R.G. Ghanem, A stochastic projection method for fluid flow, J. Comput. Phys. 173 (2001) 481. [14] M.M. Modest, Radiative Heat Transfer, Academic Press, 2003. [15] M.M.R. Williams, Radiant heat transfer through an aerosol suspended in a transparent gas, IMA J. Appl. Math. 31 (1983) 37. [16] M.M.R. Williams, Radiant heat transfer through an aerosol suspended in a transparent gas: addendum, IMA J. Appl. Math. 33 (1984) 101. [17] R. Siegel, J.R. Howell, Thermal Radiation Heat Transfer, second ed., McGrawHill Book Co., 1981. [18] D. Xiu, G.E. Karniadakis, The Wiener–Askey polynomial chaos for stochastic differential equations, SIAM J. Sci. Comput. 24 (2002) 619. [19] M.M.R. Williams, Polynomial chaos functions and stochastic differential equations, Ann. Nucl. Energy 33 (2006) 774. [20] M.M.R. Williams, Polynomial chaos functions and neutron diffusion, Nucl. Sci. Eng. 155 (2007) 109. [21] A.T. Bharucha-Reid, Probabilist. Methods Appl. Math., vol. 1, Academic Press, 1968. [22] G. Bell, S. Glasstone, Nuclear Reactor Theory, Van Nostrand, 1970. [23] Mathematical Subroutines IMSL Library, Visual Numerics Inc., 1991. [24] R.C. Weast, M.J. Astle, W.H. Beyer, CRC Handbook of Chemistry and Physics, 66th ed., CRC Press Inc., 1985.