Simultaneous estimation of strength and position of a heat source in a participating medium using DE algorithm

Journal of Quantitative Spectroscopy & Radiative Transfer 127 (2013) 130–139

Contents lists available at SciVerse ScienceDirect

Journal of Quantitative Spectroscopy & Radiative Transfer journal homepage: www.elsevier.com/locate/jqsrt

Simultaneous estimation of strength and position of a heat source in a participating medium using DE algorithm Ajit K. Parwani, Prabal Talukdar n, P.M.V. Subbarao Department of Mechanical Engineering, Indian Institute of Technology Delhi, New Delhi 110016, India

a r t i c l e in f o

abstract

Article history: Received 8 January 2013 Received in revised form 12 April 2013 Accepted 11 May 2013 Available online 18 May 2013

An inverse heat transfer problem is discussed to estimate simultaneously the unknown position and timewise varying strength of a heat source by utilizing differential evolution approach. A two dimensional enclosure with isothermal and black boundaries containing non-scattering, absorbing and emitting gray medium is considered. Both radiation and conduction heat transfer are included. No prior information is used for the functional form of timewise varying strength of heat source. The finite volume method is used to solve the radiative transfer equation and the energy equation. In this work, instead of measured data, some temperature data required in the solution of the inverse problem are taken from the solution of the direct problem. The effect of measurement errors on the accuracy of estimation is examined by introducing errors in the temperature data of the direct problem. The prediction of source strength and its position by the differential evolution (DE) algorithm is found to be quite reasonable. & 2013 Elsevier Ltd. All rights reserved.

Keywords: Inverse radiation Participating medium Differential evolution algorithm Simultaneous estimation Finite volume method

1. Introduction Inverse techniques are applied when the direct measurements of a desired quantity become difficult. The inverse technique can be used to estimate the boundary or inlet conditions, thermal properties or source term distributions in a medium. The inverse problems are mathematically classified as ill-posed in a general sense, because their solutions may become unstable because of the errors inherent in the measurements used in the analysis. To overcome the instability of the inverse problem different methods have been developed [1]. Inverse problems for the estimation of source terms in a heat transfer medium especially in conduction or/and radiation have been a subject of major interests of different research groups, because of its importance in many practical applications, such as predicting the energy generation in a computer chip or in a microwave heating process, n

Corresponding author. Tel.: +91 11 26596337. E-mail addresses: [email protected], [email protected] (P. Talukdar). 0022-4073/$ - see front matter & 2013 Elsevier Ltd. All rights reserved. http://dx.doi.org/10.1016/j.jqsrt.2013.05.006

predicting and controlling chemical reactions and phasechange etc. A considerable amount of work has been reported for determining the radiative properties [2–13] of a medium or the boundary condition [14,15] from various types of measurements. The estimation of the heat source with known location was also investigated by a number of researchers in conduction or diffusive problems [16–22] and in radiation problems [23–31]. The conjugate gradient method (CGM) with an adjoint equation is however mostly employed for the estimation of source term in radiation problems. The CGM is also called an iterative regularization method, which means the regularization procedure is performed during the iterative processes and thus the determination of optimal regularization conditions is not needed. The CGM derives from the perturbation principles and transforms the inverse problem to the solution of three problems, namely, the direct, sensitivity and the adjoint problems [1,27,28,32,33]. Linhua et al. [25] estimated source term in one-dimensional semitransparent plane-parallel media with opaque and specularly reflecting boundaries utilizing CGM. They employed the discrete ordinates method to solve the radiative transfer

A.K. Parwani et al. / Journal of Quantitative Spectroscopy & Radiative Transfer 127 (2013) 130–139

Nomenclature cp RMS F1, F2 I Ib J k M N Nt qR ! r s s^ So t tf T Tc Tm Tw ! u i;G ! v i;G x, y

specific heat capacity of fluid (J/kg K) RMS error mutation scale factors intensity (W/m2 sr) blackbody intensity (W/m2 sr) objective function, defined by Eq. (6) thermal conductivity (W/m K) number of measured data conduction radiation parameter number of time steps radiative heat flux (W/m2) position vector distance (m) unit vector in given direction heat source function (W/m) time (s) final time (s) temperature (K) temperature obtained with estimated flux (K) measured temperature by sensors (K) wall temperature (K) trial vector mutation vector coordinate directions

equation and its adjoint equation. Li [26] employed the CGM to estimate the unknown source term in a twodimensional rectangular medium with transparent boundaries. A few works in three dimensional geometries can be seen [27–30] for the estimation of the strength of heat source utilizing the CGM. Park and Yoo [27] and Park and Lee [28] estimated the strength of a transient heat source which mimics flame in a furnace from temperature measurements in the domain. Due to the null gradient at the final time step, difficulty arises to estimate the final condition [27,28]. Generally a combination of modified and regular CGM is used to alleviate this difficulty reasonably [27,28]. Apart from simple three dimensional geometries, Liu and his co-workers [29,30] solved inverse problem of source estimation for complex geometries. However there are a few drawbacks with the CGM as they can converge to local minima or maxima and there is a requirement of an initial guess value. CGM shows difficulty in the estimation of sharp peaks or discontinuities [1,32,33] on the functional variation of the source strength. In contrast search-based methods or stochastic methods like genetic algorithm (GA), particle swarm optimisation (PSO) and differential evolution (DE) have outstanding characteristics. These methods have advantages over gradient method for (i) non-requirement of initial guess and (ii) absence of complex formulations of sensitivity, adjoint and gradient equations that are inherently present in the gradient based methods. Moreover these search based methods converge to a global value. In general, GA uses the values of the parameters (to be estimated) from the given range. Variables are

131

Greek symbols β δn ε κ ω ϕ ρ s sb θ

extinction coefficient (m−1) dirac delta function convergence criteria absorption coefficient (m−1) random variable (K) azimuthal angle (rad) density of fluid (kg/m3) standard deviation Stefan–Boltzmann constant (W/m2 K4) polar angle (rad)

Subscripts G i j m

current iteration index for population index for time step measurement point

Superscripts n

representation of position

coded in binary strings. The population is operated upon by three main operators namely reproduction, crossover and mutation—to create a new population of points. On the other hand PSO does not use crossover and mutation [34–37]. PSO has been successfully applied to solve a lot of practical applications such as function optimization, artificial neural network training, pattern recognition, fuzzy control and some other fields in recent years [35]. Although DE uses a similar population based computing strategy like GA, but contrary to GA, in this method, a real parameter representation is used and an individual is formed by a vector array of all the variables in the problem. There have been works seen in the literature where these search-based methods have been used for finding different properties [8–13,34–36] of the medium as well as source strength [22,31,38]. Li and Yang [8], Verma and Balaji [9] and Das et al. [10] estimated various radiation properties using GA. For estimation of transient source in a conduction problem, Liu [22] used modified GA. Recently Lobato et al. used DE approach for inverse problems of estimation of radiative properties [11–13] in a participating medium and estimation of source strength and position [38] in a conduction problem. Parwani et al. [31] also used hybrid DE approach for estimation of source strength in a conduction–radiation problem. The estimation of strength and position of internal heat source exhibits a practical strong appeal in areas such as chemical and mechanical engineering. Some of the available works in the literature determine either the location or the strength of a heat source [39]. There are also a few literature where simultaneous estimation of location and

132

A.K. Parwani et al. / Journal of Quantitative Spectroscopy & Radiative Transfer 127 (2013) 130–139

strength [38,40–46] of the source was determined. These studies are limited to conduction heat transfer only. Karami and Hematiyan [39] considered 2D heat conduction case to estimate the location or the strength of multiple point heat sources. Only one of these two variables was estimated, the other being known utilizing the boundary element method (BEM). Lobato et al. [38] applied DE and Neto and Ozisik [40] applied CGM in one dimensional heat conduction problem while Khachfe and Jarny [41] applied the finite element method (FEM) associated to CGM in 2D heat conduction problem to estimate the location and strength of point heat sources simultaneously. Lefèvre and Niliot [42–46] also estimated simultaneously the location and strength of heat sources in heat conduction problems by applying BEM. Until now, no work has been seen for simultaneous estimation of location and strength of heat sources in a conduction–radiation problem using search based methods including DE considered in this work. In the present work such an effort has been made for simultaneous estimation of the unknown position as well as the strength of a transient heat source in a 2D participating medium with no prior information on the functional form of the strength of the heat source. The analysis is also extended for the estimation of location and strength of two heat sources. The DE is a structural algorithm or a stochastic method proposed by Storn and Price [47] for optimization problems. Besides its good convergence properties and suitability for parallelization, DE's main assets are its conceptual simplicity and ease of use. This approach is an improved version of Goldberg GA for faster optimization. Unlike CGM there is no such difficulty with DE for estimating the strength of a source at the final time step or estimating function with discontinuities or peaks. This is because in these methods, a set of solutions called population is used to converge to the optimal solution. These methods aim at finding the optima from a population of points in parallel rather than from a single point. In addition, these methods use only objective function's information instead of derivatives or other auxiliary information of the problems. The applicability and performance of DE approach for simultaneous estimation of strength and position of heat source in inverse radiation problems are tested in the present work. The DE approach is proved to be accurate for simultaneous estimation of strength and position of heat source.

2. Problem statement and formulation 2.1. Direct problem A 2D enclosure of dimension 1 m  1 m as shown in Fig. 1 is considered. Heat transfer within the enclosure is contributed by conduction as well as radiation. All boundaries of the medium are considered as black with known temperatures. The optical and thermophysical properties, such as the absorption coefficient (κ ¼0.5 m−1), extinction coefficient (β¼0.5 m−1), thermal conductivity (k¼0.05 W/ m K) and heat capacity (ρcp ¼1.2 kJ/m3 K) of the medium, are assumed to be uniform everywhere. The governing

Fig. 1. Problem description.

equations for the temperature field are described in the following subsection. 2.1.1. Energy equation The transient energy equation considering conduction, radiation and heat generation can be expressed as ρcp

∂T ¼ k∇2 T−∇⋅qR þ SoðtÞδn ðx−xn Þδn ðy−yn Þ ∂t

ð1Þ

where ∇⋅qR is the divergence of radiative heat flux, So(t) denotes the strength of heat source and δn is the dirac delta function which determines the position of heat source (xn, yn). The relevant boundary and initial conditions are At all boundaries; T ¼ T w ¼ 300 K

ð2aÞ

At t ¼ 0; T ¼ 300 K in the entire medium

ð2bÞ

The divergence of the radiative heat flux required in the energy equation (Eq. (1)) is given by Z ! ð3Þ ∇⋅qR ¼ κð4πI b − Ið r ; s^ ÞdΩÞ 4π

where I is the intensity, Ib is the blackbody intensity and Ω is the solid angle. To calculate the divergence of radiative heat flux (Eq. (3)), the solution of radiative transfer equation (RTE) is required. The RTE for a gray, absorbing and emitting medium (neglecting scattering) in the direction s^ can be written as ! ! ! ! ∇⋅ðs^ IÞ ¼ −βð r ÞIð r ; s^ Þ þ κð r ÞI b ð r Þ ð4Þ ! In Eqs. (3) and (4), r and s^ are the position vector and the unit vector describing the direction of intensity, respectively and β is the extinction coefficient which is equal to the absorption coefficient κ since scattering is not considered in this work. The discretization of the solution domain is the same for both the energy equation (Eq. (1)) as well as for the RTE (Eq. (4)) as both the equations are discretized using the finite volume method (FVM). The step scheme is used for spatial discretization of the RTE. The readers can refer to the work of Talukdar et al. [48,49] for the details of the discretized equations.

A.K. Parwani et al. / Journal of Quantitative Spectroscopy & Radiative Transfer 127 (2013) 130–139

The boundaries are isothermal (Eq. (2)) and black and can be calculated as ! Ið r ; s^ Þ ¼ I b ¼ sb T 4w =π ð5Þ The direct problem can solved to obtain the temperature field when the strength So(t) and the position (xn,yn) of the heat source are known. The numbers of control volumes considered in the spatial discretization of both the energy equation and RTE are 40  40 and the number of control angles required in angular discretization of the RTE is considered to be 10  8(θ  φ). The size of control volume is Δx¼Δy¼ 2.5 cm for the geometry illustrated in Fig. 1. The total experiment duration tf is taken as 120 s with 20 number of time steps. The computational procedure for calculation of temperature field is as follows: (1) Assume a temperature field. (2) Calculate Ib using the given temperature field. (3) Solve the RTE (Eq. (4)) using the FVM to obtain the radiation intensity I. (4) Calculate the divergence of the radiative heat flux from Eq. (3). (5) Solve Eq. (1) with the divergence of the radiative heat flux as a source term to obtain the temperature field. (6) Return to step 2 with the updated temperature field. (7) Iterate until a convergence in the temperature field is obtained. 2.2. Inverse problem In the inverse problem, the strength of the source So(t) and its position (xn, yn) are unknown and are to be estimated. The additional information which is required for the solution of this inverse problem is some temperature data in the solution domain. In a practical problem, these temperature data are found by measuring the temperature of the enclosure at different locations using some sensors. In this work, instead of measured data, some temperature data are taken from the solution of the direct problem described in Section 2.1. The solution of the present inverse problem is to be sought in such a way that the objective function J(So,(xn, yn)) is minimized. The objective function J is expressed as Z tf M JðSo; ðxn ; yn ÞÞ ¼ ∑ ½T c ðxm ; ym ; tÞ−T m ðxm ; ym ; tÞ2 dt ð6Þ m¼1

0

where M is the number of the mounted sensors and Tc(xm, ym,t) are the temperatures at the measurement locations with estimated source strength So(t) and position (xn, yn). Tm(xm,ym,t) are the measured temperatures taken at discrete locations at fixed intervals of time. For the sake of analysis, we assume the measurement to be continuous and some function of time. 3. Differential evolution DE is a parallel direct search method which utilizes Np D-dimensional parameter vectors [47] as a population for each generation G, i.e. for each iteration of the optimization. Here dimension D corresponds to the number of time

133

steps Nt of source strength and the x and y dimensional values for the source position. The initial population (total of Np) is chosen randomly between the lower and the upper bounds for the source strength for Nt number of time steps as h i ! S oi;G ¼ Soi;1;G ; Soi;2;G ⋯ Soi;Nt;G ; i ¼ 1; 2; 3⋯Np ð7aÞ Similarly the initial population for source position (independent of number of time steps Nt) is generated as !n x i;G ¼ ½xni;G ; yni;G ;

i ¼ 1; 2; 3⋯Np

ð7bÞ

The population should try to cover the entire parameter space uniformly. The value of the fitness function is calculated by substituting the variables in the objective ! !n function (Eq. (6)). S obest;G and x best;G represent the best vectors or the individuals which give minimum objective function values in Np population of source strength and position respectively in generation G. The crucial idea behind DE is a scheme for generating trial parameter vectors. DE algorithm consists of three genetic operators: mutation, crossover and selection. After obtaining the best vector from the initial population, the population is operated by mutation operators (Eq. (8)) and crossover operators (Eq. (9)) to generate new individuals and selection operator (Eq. (10)) determines the suitable individual which gives minimum value of objective function, and in this way population gets the better individuals. The mutation operator used in this work for generating new individuals for the strength and position of the source is [50] ! ! ! ! v i;G ¼ S oi;G þ F 1 ⋅ð S obest;G − S oi;G Þ ! ! þF 2 ⋅ð S or1;G − S or2;G Þ

ð8aÞ

!n !n !n !n v i;G ¼ x i;G þ F 1 ⋅ð x best;G − x i;G Þ !n !n þF 2 ⋅ð x r1;G − x r2;G Þ

ð8bÞ

where the indexes r1 and r2 represent the random and mutually different integers respectively generated within the range (1, Np) and also different from index i. F1 and F2 are the mutation scale factors usually less than 1. Fig. 2 shows a two dimensional example showing how an individual at the worst position comes to a better position after mutation operation. After mutation, the trial individuals are generated using the following scheme in the crossover procedure: ( if randð0; 1Þ≤CR vi;j;G ui;j;G ¼ Soi;j;G otherwise i ¼ 1; 2; …; Np and j ¼ 1; 2; …; D; ( uni;G ¼

vn i;G xn i;G

if randð0; 1Þ≤CR otherwise

i ¼ 1; 2; …; Np

ð9aÞ ð9bÞ

where rand(0,1) is a random number within range (0, 1) and CR is a crossover parameter taken to be 0.5 in this work and presents the probability of creating parameters for the trial vector from a mutant vector. It is to be noted that an individual already at a better position may or may

134

A.K. Parwani et al. / Journal of Quantitative Spectroscopy & Radiative Transfer 127 (2013) 130–139

not get an improved position by mutation and crossover operations. Therefore to ensure optimization to be global minimum selection procedure is used. In the selection procedure for each target individual of ! !n strength S oi;G and target individual of position x i;G the fitness values of the trial individual of strength and ! !n position u i;G and u i;G are compared, and the individual with the minimum values of objective function (Eq. (6)) is selected for the next generation. The selection is given as 8 ! !n ! !n !n
ð10Þ

3.1. Stopping criteria The stopping criterion is specified as JðSo; ðxn ; yn ÞÞ oε

ð11Þ

where ε is a small predefined constant. In this work ε is taken as 5  10−5. The observed temperature data may contain measurement errors. Therefore, we do not expect the functional (Eq. (6)) to be equal to zero at the final iteration step. We use the discrepancy principle as the stopping criterion. The following expression is obtained for the stopping criteria ε [1]: ε ¼ Ms2 t f

ð12Þ

where s is the standard deviation of measurements. The overall computational procedure is given in the flow chart (Fig. 3). The parameters used for DE algorithm in this work are Np¼ 10, F1 ¼ 0.8 and F2 randomly varying in range (0, 1). The upper and lower limits of source strength for inverse problem are taken as 25,000 W/m and 0 W/m respectively at all time steps. The upper and lower limits of source position taken are 1 m and 0 m respectively at all time steps.

4. Results and discussion

Fig. 2. An example of two-dimensional objective function showing its ! contour lines and the process for generating mutation vector v i;G ! marked “o” for an individual vector S oi;G marked “x1”.

The objective of this work is to estimate simultaneously the strength and position of a transient heat source in a 2D participating medium without any prior information on the functional form of the source strength. The unknown strength of the source is assumed to be any of the different forms of function (sinusoidal, triangular or step) listed below to examine the accuracy of the proposed model. SoðtÞ ¼ 2  104 sin ðπðj−1Þ=19Þ

Fig. 3. Computational procedure.

for

j ¼ 1; 2; …; 20

ð13Þ

A.K. Parwani et al. / Journal of Quantitative Spectroscopy & Radiative Transfer 127 (2013) 130–139

ð14Þ  SoðtÞ ¼

25; 000

for j ¼ 6; 7; …; 14

8000

otherwise

ð15Þ

where j is the index of time step. The effects of the number of measurements on the accuracy of estimations are investigated. Fig. 1 shows schematic of the arrangement of nine sensors (M¼9) with spacing of 0.25 m between each two sensors in the x and y directions. In another sensor arrangement of M ¼15, these are arranged in three columns with five sensors in each column such that spacing between each two sensors in the y direction is 0.125 m and in the x direction is 0.25 m. Similarly for M¼21, seven sensors are arranged in three columns such that spacing between each two sensors in the y direction is 0.1 m. In experiments, there could be errors in the measurement of temperature. Since in the present work, these temperature data are taken from the solution of direct problem, they are error free. In order to be more realistic, errors are added randomly in these temperature data so as to represent them more like an experimental data. Thus, the measured temperature data Tm are represented by adding random errors to the exact temperature Te computed from the solution of the direct problem as T m ¼ T e þ ωs

20 18 16 14 12 exact source strength

10 8

estimated source strength

6 4

X

2 0 0

20

40

60 time (s)

80

100

120

Fig. 4. The number of iterations G and RMS value of exact and estimated heat source function So varying sinusoidal with time with different numbers of noiseless measurements M.

ð16Þ

where ω is a random variable of normal distribution with zero mean and unit standard deviation. The value of ω is calculated by the Fortran subroutine ‘random_number’ and chosen over the range −2.576oωo2.576. Standard deviation, s, can be related to relative measurement error [27,28] as follows: erel ¼ 2:576ðs=T e Þ

for satisfying convergence criterion with M¼9, 15 and 21 are 38,154 s, 37,330 s and 38,426 s respectively. The estimated position of heat source (xn, yn) in the medium with the number of measurements is indicated in Table 1. The relative errors in the x and y coordinates and RMS error for the estimation of position of source are also calculated. It is seen that higher number of measurements do not reduce RMS errors because in the inverse analysis

So (kW/m)

8 7900 for j ¼ 1; …4 or 17; …20 > < 4 for j ¼ 5; …; 10 SoðtÞ ¼ 5  10 ððj−1Þ=19Þ; > : 5  104 ð1−ðj−1Þ=19Þ; otherwise

135

ð17Þ

The value of s is adjusted such that erel is zero, 1% and 5%. The RMS error is calculated to compare the quality of estimation of strength of the source as sffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi D 1 RMS ¼ ð18Þ ∑ ½Soex ðtÞ−Soest ðtÞ2 D i¼1 where Soex and Soest are the exact and estimated strength of source values respectively. As a first test case of the present work, an idealized situation is considered in which there is no measurement error (erel ¼0). Fig. 4 presents the estimated results obtained with different number of measurements for sinusoidal profile of a heat source function So(t). In the same figure, the number of generations and RMS error are also indicated. The estimated profiles are almost coinciding to the actual profile and maximum RMS error is just 0.21 with both M¼ 9 and 15. The minimum RMS error obtained with M¼21 is 0.17. The simulations are performed on a workstation with 2.67 GHz Intel (R) Xeon (R) X5550 processor and 8GB of RAM. The CPU time taken

Table 1 Relative errors (absolute) and RMS values for the estimation of source position. Exact position of source So (x* ¼0.4625, y* ¼0.2125) x*

y*

Relative error— Relative error— RMS x (%) y (%)

Predicted M¼9 M ¼ 15 M ¼ 21

position 0.5074 0.4969 0.4891

of source 0.2617 0.2380 0.2608

So for the case shown in Fig. 4 9.7081 23.1529 7.4378 12 5.7514 22.7294

Predicted M¼9 M ¼ 15 M ¼ 21

position 0.5090 0.5058 0.4884

of source So for the case shown in Fig. 5(a) 0.2512 10.0541 18.2118 0.0302 0.2489 9.3622 17.1294 0.0283 0.2494 5.6000 17.3647 0.0225

Predicted M¼9 M ¼ 15 M ¼ 21

position 0.5099 0.5069 0.5069

of source So for the case shown in Fig. 5(b) 0.2622 10.2486 23.3882 0.0343 0.2405 9.6000 13.1765 0.0262 0.2435 9.6000 14.5882 0.0271

Predicted M¼9 M ¼ 15 M ¼ 21

position 0.5074 0.4891 0.4969

of source 0.2608 0.2380 0.2617

So for the case shown in Fig. 6(a) 9.7081 22.7294 0.0330 5.7514 12.0000 0.0184 7.4378 23.1529 0.0300

Predicted M¼9 M ¼ 15 M ¼ 21

position 0.5074 0.5000 0.5027

of source 0.2555 0.2541 0.2554

So for the case shown in Fig. 6(b) 9.7081 20.2353 0.0311 8.1081 19.5765 0.0280 8.6919 20.1882 0.0294

Predicted N ¼ 0.004 N ¼ 0.01 N ¼ 0.1

position 0.5058 0.5001 0.5091

of source So for the case shown in Fig. 7 0.2489 9.3622 17.1294 0.2550 8.1297 20.0000 0.2526 10.0757 18.8706

0.0333 0.0214 0.0276

0.0283 0.0284 0.0307

136

A.K. Parwani et al. / Journal of Quantitative Spectroscopy & Radiative Transfer 127 (2013) 130–139

with DE algorithm, the locations are randomly chosen and actually they may not be coinciding with the coordinates of the centroid of a control volume used in the FVM. But, in the numerical model it has been assumed that whenever the randomly chosen source position falls inside a control volume, the centroid of that particular control volume is assumed to be the modified position of the point source and the calculations of the direct and inverse problem are done based on that modified source location. However, at the final iteration, the estimated position of source (Fig. 4 and all the subsequent figures) is the actual estimation obtained for the source position (without modifying to centroid of a control volume) determined by the DE algorithm. The next case examines the effect of sharp peak or discontinuity present in the function of source strength. This is seen to be one of the most difficult cases to be estimated with CGM [1,32,33]. Fig. 5(a) and (b) presents the results of the triangular (Eq. (14)) and step functions (Eq. (15)) of heat source with time respectively for different numbers of measurements. It is seen that the quality of the results is nearly similar to the sinusoidal varying function presented in Fig. 4. The maximum RMS error is 0.22 in triangular function and 0.24 in step function with M ¼9. The estimated position of heat source with the number of measurements is shown in Table 1 and it

indicates that even functions with peaks or discontinuities can be well estimated by the DE algorithm. The measurement errors on the accuracy of inverse estimation are examined. Fig. 6(a) shows the exact and estimated profiles of the sinusoidal heat source function So(t) with measurement error erel ¼1% and (b) shows the estimation with measurement error erel ¼ 5%. The result shows a better prediction with M¼15 where the RMS error is close to 14 and 114 for erel ¼1% and 5% respectively. The result shows that with M ¼21, the RMS errors are maximum and are close to 90 and 293 for erel ¼1% and 5% respectively The estimated position of heat source as shown in Table 1 with all the measurements appears to be very close to the exact position in both the cases. The next consideration is the effect of variation of conduction–radiation parameter (N ¼kβ/4sbT3) on the accuracy of estimation. Fig. 7 shows the exact and estimated profiles of the triangular function of heat source with 15 noiseless measurements and for different values of N. In addition to the default value of N ( ¼0.004), two other values of N (¼ 0.01 and 0.1) are considered in this study. When the value of N increases, the contribution of conduction heat transfer becomes much more dominant than radiation heat transfer. It is observed from the results that as the value of N increases the RMS error increases. This is

18

26 24 22 20 18 16 14 12 10 8 6 4 2 0

16

exact source strength

12 10

exact source strength

8

estimated source strength

6 estimated source strength

4 X

2

X 0

So (kW/m)

So (kW/m)

14

0

20

40

60 time (s)

80

100

0

120

20

40

60 time (s)

80

100

120

22

26 24 22 20 18 16 14 12 10 8 6 4 2 0

20 18 16 So (kW/m)

So (kW/m)

20

exact source strength

estimated source strength

20

12

exact source strength

10 8

estimated source strength

6 4

X 0

14

2 40

60 time (s)

80

100

120

Fig. 5. The number of iterations G and RMS value of exact and estimated heat source function So varying with time for (a) triangular function and (b) function having two discontinuities with different numbers of noiseless measurements M.

0 0

20

40

60 time (s)

80

100

120

Fig. 6. The number of iterations G and RMS value of exact and estimated heat source function So varying sinusoidal with error introduced in measurements M (a) erel ¼1% and (b) erel ¼5%.

So (kW/m)

A.K. Parwani et al. / Journal of Quantitative Spectroscopy & Radiative Transfer 127 (2013) 130–139

26 24 22 20 18 16 14 12 10 8 6 4 2 0

137

Table 2 Relative errors (absolute) and RMS values for the estimation of source position. Exact position of source So ( x* ¼ 0.4875, y* ¼ 0.7375) x*

y*

Relative error—x (%)

Relative error—y (%)

Predicted position of source So for the case shown in Fig. 4 M ¼ 9 0.5056 0.7607 3.7128 3.1458 M ¼ 15 0.4972 0.7601 1.9897 3.0644 M ¼ 21 0.4980 0.7587 2.1538 2.8746

exact source strength estimated source strength

RMS

0.0147 0.0123 0.0118

X 0

20

40

60 time (s)

80

100

120

Fig. 7. The number of iterations G and RMS value of exact and estimated heat source function So varying with time in triangular form with different values of conduction radiation parameter N and 15 noiseless measurements.

x

x

x

x

x

x So2

x

x

x

x

x

x

x

x

x

So1

18

20

16

18

14

16 14

12 So (kW/m)

So (kW/m)

20

exact source strength 10 8 estimated source strength

6

12 10 8 6

4

estimated source strengths

4

X

X

2

2

0

0 0

20

40

60 time (s)

80

100

120

Fig. 8. The number of iterations G and RMS value of exact and estimated heat source function So varying sinusoidal with different numbers of noiseless measurements M.

due to the fact that the sensitivity value becomes smaller when the conduction heat transfer becomes more dominant compared to radiative heat transfer. This was also pointed out by Kim and Baek [15]. To check the generality of the DE algorithm for the estimation of strength of heat source at any position in the enclosure, a case is considered in which the heat source is placed at xn ¼ 0.4875 m and yn ¼0.7375 m. Fig. 8 shows the result for this case with noiseless measurements (erel ¼0) and Table 2 shows the estimated position of heat source. It can be concluded from these results that the DE algorithm has the ability to estimate the strength of heat source at any position. DE algorithm can also be applied for simultaneous estimation of two or more heat sources. Fig. 9(a) shows schematic of two heat sources whose strengths are varying

0

20

40

60 time (s)

80

100

120

Fig. 9. (a) Schematic of two heat sources So1 and So2. (b) The exact and estimated sinusoidal heat source strength profiles estimated simultaneously with 15 noiseless measurements.

sinusoidal with time. These strength and positions are estimated with 15 noiseless measurements. The second heat source So2 varies with time in functional form of half the value of the first heat source So1 which varies with time in functional form as given by Eq. (13). Fig. 9(b) shows the results for the estimation of source strengths and Table 3 shows the estimation of source positions for this case. Since double dimensions are required in searching for simultaneous estimation of two heat sources, it takes 149,368 s CPU time and 1512 number of iterations to reach the convergence. 5. Conclusions The DE algorithm is successfully applied for the solution of inverse radiation problem to estimate the position

138

A.K. Parwani et al. / Journal of Quantitative Spectroscopy & Radiative Transfer 127 (2013) 130–139

Table 3 Relative errors (absolute) and RMS values for the simultaneous estimation of position of two sources. n

n

Exact position of sources So1 ( x1 ¼0.4625, y1 ¼ 0.2125) and So2 n n (x2 ¼ 0.7125, y2 ¼ 0.3375) xn

yn

relative error—x (%)

relative error—y (%)

Predicted position of source for the case shown in Fig. 9 So1 0.5025 0.2507 8.65 17.98 So2 0.7602 0.3661 6.69 8.47

RMS

0.0276 0.0278

and strength of an unknown source simultaneously in a 2D participating media. Several test cases involving different functions of heat source, different numbers of thermocouples and effect of measurement errors are considered. The following are the main conclusions:

 The DE algorithm can accurately estimate the strength  

 



of a function of heat source even with peaks or discontinuities. The DE algorithm does not show any limitation in predicting the strength of source in the final time which is not the case with regular CGM. Higher the number of noiseless measurements lower is the RMS error. It is also observed when the error in the measurements (temperature data from direct problem) are added, 15 sets of measurements were found to give better estimation. Higher the thermal conductivity higher is the RMS error. The estimated position of heat source approaches very near to the exact position with all the numbers of measurements considered and even with the consideration of error. DE algorithm can also be applied for simultaneous estimation of two or more heat sources. The overall results show that the inverse solutions obtained by DE remain stable and are found to be quite reasonable.

Acknowledgment This work is partially supported by a fund from the SERC division of Department of Science and Technology, Government of India. The authors greatly appreciate the financial contribution toward this research. References [1] Ozisik MN, Orlande HRB. Inverse heat transfer: fundamentals and applications. New York: Taylor and Francis; 2000. [2] Sanchez R, McCormick NJ. Numerical evaluation of optical singlescattering inverse transport methods. J Quant Spectrosc Radiat Transfer 1982;28:169–84. [3] Dunn WL. Inverse Monte Carlo solutions for radiative transfer in inhomogeneous media. J Quant Spectrosc Radiat Transfer 1983;29: 19–26.

[4] Kamiuto K, Seki J. Study of the P1 approximation in an inverse scattering problem. J Quant Spectrosc Radiat Transfer 1987;37: 455–9. [5] Kamiuto K, Iwanmoto M, Nishimura T, Sato M. Albedos and asymmetry factors of the phase functions for packed-sphere systems. J Quant Spectrosc Radiat Transfer 1991;46:309–16. [6] Kamiuto K. Emerging-flux fitting method for inverse scattering problems. J Quant Spectrosc Radiat Transfer 1993;49:263–7. [7] Silva Neto AJ, Ozisik MN. An inverse problem of simultaneous estimation of radiation phase function, albedo and optical thickness. J Quant Spectrosc Radiat Transfer 1995;53:397–409. [8] Li HY, Yang CY. A genetic algorithm for inverse radiation problems. Int J Heat Mass Transfer 1997;40(7):1545–9. [9] Verma S, Balaji C. Multi-parameter estimation in combined conduction–radiation from a plane parallel participating medium using genetic algorithms. Int J Heat Mass Transfer 2007;50:1706–14. [10] Das R, Mishra SC, Ajith M, Uppaluri R. An inverse analysis of a transient 2-D conduction–radiation problem using the lattice Boltzmann method and the finite volume method coupled with the genetic algorithm. J Quant Spectrosc Radiat Transfer 2008;109: 2060–77. [11] Lobato FS, Steffen Jr. V, Silva Neto AJ. A comparative study of the application of differential evolution and simulated annealing in inverse radiative transfer problems. J Braz Soc Mech Sci Eng 2010;32(5):518–26. [12] Lobato FS, Steffen Jr. V, Silva Neto AJ. Self-Adaptive differential evolution based on the concept of population diversity applied to simultaneous estimation of anisotropic scattering phase function, albedo and optical thickness. Comput Model Eng Sci 2010;69(1): 1–18. [13] Lobato FS, Steffen Jr V, Silva Neto AJ. Estimation of space dependent single scattering albedo in radiative transfer problems using the differential evolution algorithm. In: Inverse problems design and optimization symposium, IPDO 2010. João Pessoa, Brazil; 2010. [14] Barichello LB, Garcia RDM, Siewert CE. On inverse boundarycondition problems in radiative transfer. J Quant Spectrosc Radiat Transfer 1997;57:405–10. [15] Kim KW, Baek SW. Inverse radiation–conduction design problem in a participating concentric cylindrical medium. Int J Heat Mass Transfer 2007;50:2828–37. [16] Silva Neto AJ, Özisik MN. Inverse problem of simultaneously estimating the timewise-varying strengths of two plane heat sources. J Appl Phys 1993;73(5):2132–7. [17] Silva Neto AJ, Özisik MN. Two-dimensional inverse heat conduction problem of estimating the time-vayring strength of a line heat source. J. Appl Phys 1992;71:5357–62. [18] Yang CY. A sequential method to estimate the strength of the heat source based on symbolic computation. Int J Heat Mass Transfer 1998;41(14):2245–52. [19] Yang CY. The determination of two heat sources in an inverse heat conduction problem. Int J Heat Mass Transfer 1999;42:345–56. [20] Niliot CL. The boundary element method for the time varying estimation of point heat sources: application to a two dimensional diffusion system. Numer Heat Transfer, Part B 1998;33:301–21. [21] Huang CH, Cheng SC. A three-dimensional inverse problem of estimating the volumetric heat generation for a composite material. Numer Heat Transfer, Pt A 2001;39:383–403. [22] Liu FB. A modified genetic algorithm for solving the inverse heat transfer problem of estimating plan heat source. Int J Heat Mass Transfer 2008;51:3745–52. [23] Siewert CE. An inverse source problem in radiative transfer. J Quant Spectrosc Radiat Transfer 1993;50:603–9. [24] Siewert CE. A radiative-transfer inverse-source problem for a sphere. J Quant Spectrosc Radiat Transfer 1994;52:157–60. [25] Linhua L, Heping T, Qizheng Y. Inverse radiation problem in onedimensional semitransparent plane-parallel media with opaque and specularly reflecting boundaries. J Quant Spectrosc Radiat Transfer 2000;64:395–407. [26] Li HY. Inverse radiation problem in two-dimensional rectangular media. J Thermophys Heat Transfer 1997;11:556–61. [27] Park HM, Yoo DH. A multidimensional inverse radiation problem of estimating the strength of a heat source in participating media. Int J Heat Mass Transfer 2001;44:2949–56. [28] Park HM, Lee WJ. The solution of inverse radiation problems using an efficient computational technique. J Quant Spectrosc Radiat Transfer 2002;73:41–54. [29] Liu LH, Tan HP. Inverse radiation problem in three-dimensional complicated geometric systems with opaque boundaries. J Quant Spectrosc Radiat Transfer 2001;68:559–73.

A.K. Parwani et al. / Journal of Quantitative Spectroscopy & Radiative Transfer 127 (2013) 130–139

[30] Liu LH, Tan HP, He ZH. Inverse radiation problem of source term in three-dimensional complicated geometric semitransparent media. Int J Therm Sci 2001;40(6):528–38. [31] Parwani AK, Talukdar P, Subbarao PMV. Performance evaluation of hybrid differential evolution approach for estimation of the strength of a heat source in a radiatively participating medium. Int J Heat Mass Transfer 2013;56:552–60. [32] Parwani AK, Talukdar P, Subbarao PMV. Estimation of inlet temperature of a developing fluid flow in parallel plate channel. Int J Therm Sci 2012;57:126–34. [33] Parwani AK, Talukdar P, Subbarao PMV. Estimation of transient boundary flux for a developing flow in a parallel plate channel. Int J Numer Methods Heat Fluid Flow, in press. [34] Becceneri JC, Stephany S, de Campos Velho HF, Silva Neto AJ. Solution of the inverse problem of radiative properties estimation with particle swarm optimization techniques. In: Proceedings of the Inverse Problems Eng Sem (IPES). USA: Iowa State University; 2006. [35] Lee KH, Baek SW, Kim KW. Inverse radiation analysis using repulsive particle swarm optimization algorithm. Int J Heat Mass Transfer 2008;51(11–12):2772–83. [36] Qi H, Wang DL, Wanga SG, Ruan LM. Inverse transient radiation analysis in one-dimensional non-homogeneous participating slabs using particle swarm optimization algorithms. J Quant Spectrosc Radiat Transfer 2011;112:2507–19. [37] Farahmand A, Payana S, Hosseini Sarvari SM. Geometric optimization of radiative enclosures using PSO algorithm. Int J Therm Sci 2012;60:61–9. [38] Lobato FS, Steffen Jr V, Silva Neto AJ. Simultaneous estimation of location and timewise varying strength of a plane heat source using differential evolution. In: Proceedings of the 2nd international conference on engineering optimization. Lisbon, Portugal; 2010. [39] Karami G, Hematiyan MR. A boundary element method of inverse non-linear heat conduction analysis with point and line heat sources. Commun Numer Methods Eng 2000;16:191–203. [40] Silva Neto AJ, Ozisik MN. Simultaneous estimation of location and timewise varying strength of a plane heat source. Numer Heat Transfer 1993;24:467–77.

139

[41] Abou khachfe R, Jarny Y. Determination of heat sources and heat transfer coefficient for two-dimensional heat flow: numerical and experimental study. Int J Heat Mass Transfer 2001;44(7):1309–22. [42] Lefèvre F, Niliot CL. Multiple transient point heat sources identification in heat diffusion: application to experimental 2D problems. Int J Heat Mass Transfer 2002;45:1951–64. [43] Lefèvre F, Niliot CL. The BEM for point heat sources estimation: application to multiple static sources and moving sources. Int J Therm Sci 2002;41:526–45. [44] Lefèvre F, Niliot CL. A boundary element inverse formulation for multiple point heat sources estimation in a diffusive system: application to a 2D experiment. Inverse Prob Eng 2002;10(6): 539–57. [45] Niliot CL, Lefèvre F. A parameter estimation approach to solve the inverse problem of point heat sources identification. Int J Heat Mass Transfer 2004;47:827–41. [46] Niliot CL, Lefèvre F. Multiple transient point heat sources identification in heat diffusion: application to numerical two- and threedimensional problems. Numer Heat Transfer, Pt B 2001;39:277–301. [47] Storn R, Price K. Minimizing the real functions of the ICEC'96 contest by Differential Evolution. In: Proceedings of the IEEE international conference on evolutionary computation (ICEC'96). Nagoya, Japan; 1996. p. 842–4. [48] Talukdar P, Steven M, Issendorff FV, Trimis D. Finite volume method in 3-D curvilinear coordinates with multiblocking procedure for radiative transport problems. Int J Heat Mass Transfer 2005;48: 4657–66. [49] Talukdar P, Issendorff FV, Trimis D, Simonson CJ. Conduction– radiation interaction in 3D irregular enclosures using the finite volume method. Heat Mass Transfer 2008;44:695–704. [50] Brest J, Zumer V, Maucec MS. Self-adaptive differential evolution algorithm in constrained real-parameter optimization. In: Proceedings of the IEEE congress on evolutionary computation (CEC 2006). BC, Canada: Vancouver; 2006. p. 215–22.