Application of inverse analysis to determine the geometric configuration of filament heaters for uniform heating

International Journal of Thermal Sciences 105 (2016) 1e12

Contents lists available at ScienceDirect

International Journal of Thermal Sciences journal homepage: www.elsevier.com/locate/ijts

Application of inverse analysis to determine the geometric configuration of filament heaters for uniform heating rio Brittes b, Francis H.R. França a, * Larissa D. Lemos a, Roge a b

Departament of Mechanical Engineering, Federal University of Rio Grande do Sul, Rua Sarmento Leite, 425, 90050-170 Porto Alegre, RS, Brazil Federal University of Santa Maria, Av. Presidente Vargas, 1958, 96506-302 Cachoeira do Sul, RS, Brazil

a r t i c l e i n f o

a b s t r a c t

Article history: Received 1 November 2015 Received in revised form 23 February 2016 Accepted 26 February 2016 Available online xxx

In typical designs of ovens, it is often desired to determine the position and geometrical configurations of heaters to attain uniform heating of the materials subjected to thermal processes. In this study, uniform temperature and heat flux are prescribed in one of the surfaces of a radiative system, while the shape and position of filament heaters are left to be found. This is equivalent to the inverse problem of the conventional trial-and-error approach. The oven is modeled as a three-dimensional enclosure formed with diffuse-gray walls, and filled with transparent medium. The inverse problem is solved as an optimization problem. The solution is obtained by the generalized extreme optimization (GEO) algorithm, a stochastic global optimization method. The proposed methodology can lead to filament heaters capable of satisfying the two conditions on the design surface with maximum deviation of about 2% and average deviation less than 0.5%. © 2016 Elsevier Masson SAS. All rights reserved.

Keywords: Uniform heating Filament heater Inverse design Stochastic algorithm

1. Introduction Heat treatment is widely used in the industry to restore or change certain properties of materials under high temperature conditions. In general, it is required uniform rise of the temperature of the material to avoid internal thermal stresses, which could reduce the material strengthen or cause cracking. However, the local energy balance on the material requires uniform heat flux on its surface for the uniform rise of the temperature to be attained. It follows from this that both the temperature and heat flux are simultaneously prescribed in the material. This imposes a major difficulty in the standard mathematical formulation of the problem, which is based on the specification of one single boundary condition, temperature or heat flux, on each surface that forms the system. In the design of thermal systems for thermal treatment, two conditions are imposed on the surface of the material under treatment, namely the design surface. The problem consists of determining the conditions in the heaters to satisfy the process. In the conventional approach, the heaters are tentatively designed, and one of the thermal conditions on the design surface is imposed, for instance, the temperature. Then, the heat flux on the design

* Corresponding author. E-mail address: [email protected] (F.H.R. França). http://dx.doi.org/10.1016/j.ijthermalsci.2016.02.015 1290-0729/© 2016 Elsevier Masson SAS. All rights reserved.

surface is determined. If the obtained condition is not equal to the specification, a new design for the heaters is attempted, and the computations are rerun until the prescribed heat flux is obtained within some approximation. Due to the inherent non-linearity of the problem, it is very difficult to make the next design based on the results of the previous attempt. In the inverse formulation, the two conditions on the design surface are the input of the problem, and the conditions on the heaters become the output. While the inverse solution avoids the trial-and-error approach of the direct solution, its mathematical formulation is typically ill-posed, and requires special mathematical methods to recover physically acceptable solutions [1]. In the past twenty years, inverse design of radiative systems has been the subject of several studies that addressed different aspects of the problem. Regarding the conditions that are sought for the heating elements, the studies can be classified in three categories. In the first one, the geometric configurations of the heating elements (that is, their locations and shapes) are specified [2e17]; the problem consists in determining the power in the heating elements to attain the two conditions on the design surface. In the second category, the power input in the heating elements are specified based on the global energy balance in the system, but the geometry of the heating elements is left unconstrained, which is then determined from the inverse analysis [18e22]. Finally, the third approach involves the determination of both the power and the

2

L.D. Lemos et al. / International Journal of Thermal Sciences 105 (2016) 1e12

geometry of the heating elements [23e26]. Another possible classification regards the adopted methodology for the solution of the inverse problem, which could be based on the explicit or the implicit formulation. The explicit formulation involves the direct inversion of the system of equations formed by the application of the radiative energy balance in the system, leading to a system of equations that is typically ill-conditioned. This requires regularization methods to recover meaningful results. Several regularization methods have been adopted, such as the truncated singular value decomposition (TSVD) [2,3,7,10,12,26], conjugate gradient method [2,6e8,11,13,16], Tikhonov method [7,9], LevenbergeMarquardt method [4,17]. In the implicit formulation, the problem is solved by minimizing the deviation between the specified conditions on the design surface and the conditions obtained for a certain design of the heating system. This approach involves the solution of the direct problem in which the conditions in the heating elements are specified, and the heat flux on the design surface is computed and compared with the specification. Most studies in the literature have used stochastic methods to solve this minimization problem, including the micro-genetic algorithm [14,15,18,19,24], particle swarm optimization [20,21], random search [22], and the generalized extremal optimization (GEO) [25,26]. The major limitation of the explicit formulation is perhaps the inherent difficulty of its application to problems where the geometry of the heating elements is unknown, since there seems to be no simple way to establish explicit relations for the unknown geometric parameters of the heating elements. The implicit formulation, on the other hand, requires the solution of the direct problem several times, and therefore the computation demands are usually high. In the hybrid approach in Refs. [26], the positions of the heating elements were sought by optimization (GEO algorithm) while their heat inputs were determined by regularization (TSVD method), considerably reducing the computational time in comparison to the purely optimization problem. Inverse geometry designs can involve either the optimization of the geometric configuration of the radiative enclosure [20e22] or the location and shape of the heating elements to attain the specified conditions on the design surface [18,19,23e26]. In general determining the placement of the heating elements can be a simpler design solution, for it allows the use of the same enclosure configuration for different processes, only changing the placement and shape of the heating elements. However, the solutions in Refs. [19,23e26] involved several punctual heating elements, in some cases operating at different power inputs, which may be not the preferable design in practice. In the present study, it is proposed a methodology to determine the shape and location of heating filaments to allow uniform heating of the design surface. Filament heaters are widely employed in the industry, and can operate with electrical power or as heating tubes with internal recirculation of flue gases. The inverse problem is solved via implicit formulation, in which the optimization problem is solved with the GEO algorithm. The solution is presented for a 3D rectangular enclosure where the heat transfer is dominated by thermal radiation. In the first part of the solution, the filament is assembled with a number of contiguous square elements. In the final part of the solution, it is shown that simply replacing those square elements with a more realistic slender filament would not affect considerably the final results. The accuracy of the proposed methodology is assessed by the comparison of the imposed heat flux on the design surface and the one computed for the different filament configurations. 2. Physical and mathematical formulation Fig. 1(a) shows a three-dimensional rectangular radiative enclosure with length, height and width equal to L, H and W,

respectively. A similar system was considered in Refs. [25,26], but considering a set of independent, square-shaped heaters. In this particular test case, the design surface, in which both the temperature and heat fluxes are specified, is placed on the bottom of the enclosure. The filament heaters are placed on the top surface of the enclosure, but their shape configuration and position are unknown. The remaining surfaces of the enclosure are assumed to be perfectly adiabatic. Moreover, all surfaces are assumed to be diffuse-gray. As seen in the figure, the design surface is placed at some distance from the side walls. Otherwise, the borders of the design surface would be more affected by the side walls than by the heaters, and it could be impossible to find a proper design solution. Finally, it is assumed that the medium in the interior of the enclosure is nonparticipating, and that the heat transfer is solely governed by thermal radiation exchanges between the surfaces. Previous studies [3,6,10,16,17,19] have demonstrated that inclusion of combined heat transfer mode and gas radiation does not impose a major difficulty to the application of the inverse analysis. Thus, the simplifications adopted in the present study provide a simpler scenario to demonstrate the proposed methodology, which can be extended to include other effects. The present design problem can be stated in the following way. Consider that uniform temperature and heat flux are specified on the design surface shown in Fig. 1(a). Determine the shape and position of the filament heaters for these specifications to be attained. Due to the symmetry in the geometry and boundary conditions, only one-quarter of the domain needs to be solved, as seen in Fig. 1(b). One filament is placed in each quarter, which means that a total of four filament heaters would be necessary in this design. However, since each filament would require the same heat input due to symmetry, the same power system could be used to operate the heaters. Moreover, as will be seen in the results section, in some of the configurations the ends of the filaments might meet in the symmetry line, so in those cases they could be treated as a single filament. This makes the present approach of more practical interest than those in Refs. [25,26], in which a relatively large number of heaters, normally 10 heaters in each quarter of the domain, was necessary to achieve the desired precision. For the solution of the radiative exchange in the enclosure, the radiosity method [27] is employed. The surfaces that form the enclosure are subdivided into J surface elements, for which one equation is set according to the specified boundary condition on each surface element j. When the temperature Tj is specified, the radiosity or the outgoing radiative heat flux, qo,j, is given by the following relation: J
X qo; j ¼ εj eb; j þ 1  εj Fjj* qo; j*

(1)

j*¼1

where eb;j ¼ sTj4 is the blackbody emissive power, s ¼ 5.67  108 W/(m2$K4) is the constant of StefaneBoltzmann, εj is the total hemispherical emissivity of surface element j, and Fjj* is the view factor. When the net heat flux qr,j on surface element j is specified, the radiosity is given by:

qo; j ¼ qr; j þ

J X

Fjj* qo; j*

(2)

j*¼1

Equations (1) and (2) are written for the elements that are placed in the one-quarter part of the domain in Fig. 1(b). However, the elements located in the other three quarters of the enclosure also exchange radiation with the element j. To take this into account, the view factor Fjj* in the above equations is actually a

L.D. Lemos et al. / International Journal of Thermal Sciences 105 (2016) 1e12

3

Fig. 1. (a) Radiative enclosure; (b) Computation domain corresponding to one quarter of the enclosure.

combination (sum) of the view factors of element j with element j* in the one-quarter domain of Fig. 1(b) and with the symmetrical elements to j* in the other quarters of the domain. For a more general presentation of the problem, Eqs. (1) and (2) can be written in dimensionless form according to: J
X Q o; j ¼ εj tj4 þ 1  εj Fjj* Q o; j*

(3)

j*¼1

Q o; j ¼ Q r; j þ

J X

Fjj* Q o; j*

(4)

j*¼1

In the above equations, the dimensionless heat flux and radiosity are defined as Qr;j ¼ qr;j =ðsTd4 Þ and Qo;j ¼ qo;j =ðsTd4 Þ, respectively, where Td is the specified temperature on the design surface. The dimensionless temperature is defined as tj ¼ Tj/Td. In the next developments, the elements on the design surface, on the filament and on the adiabatic walls will be specified by jd, jh and jw, respectively. In the present study, the three dimensional enclosure is divided into finite-sized square elements with dimensions Dx ¼ Dy ¼ Dz, and the view factors are computed with the relatively simple correlations provided in Refs. [28], which lead to deviations less than 0.1% in comparison to the exact relations in Ref. [30]. Fig. 2(a) presents a schematic of the division of the bottom surface of the enclosure into square surface zones. Fig. 2(b) depicts one example of filament heater formed by zone elements placed on the one-

quarter top surface of the enclosure. In the present solution, the shape and position of the heater are not known, but should be found so that the conditions on the design surface are attained. In assembling the filaments, it is required that all the zone elements are contiguous. Moreover, the number of elements that form the filament, Jh, is specified, although different sizes can be tested. In the figure, the filament is formed by 20 elements, that is Jh ¼ 20. One might notice that using the uniform grid to build the filament could lead to a configuration with non-realistic width depending on the size of the square elements that form the grid. In this study, the variables are presented in dimensionless form, so the actual dimensions will depend on the characteristic length L that represents the enclosure. For this reason, a slender filament is placed over the positions of the square heating elements to represent what could be a filament with a more realistic width. In the last part of this study, it will be evaluated the deviations that are introduced when the slender filament replaces the optimum filament with square elements. Except for the design surface and the filament heaters, all the other surfaces that form the enclosure are perfectly adiabatic. Therefore, considering that the net radiative heat flux is uniform on each element that forms the filament, it can be directly determined by the application of energy balance in the enclosure, that is,

J Q Qh ¼  d d Jh

(5)

where Qh is the dimensionless net radiative heat flux on the filament heater, Qh ¼ qh =ðsTd4 Þ, Qd is the dimensionless heat flux

Fig. 2. (a) View of the design surface (shaded area); (b) Filament formed with 20 heating elements. The black line represents the equivalent slender filament.

4

L.D. Lemos et al. / International Journal of Thermal Sciences 105 (2016) 1e12

specified on the design surface, Qd ¼ qd =ðsTd4 Þ, and Jd is the number of design surface elements in the one-quarter of the enclosure. In Eq. (5), the negative sign arises from the sign convention adopted in this study, that is, since the design surface gains energy from the radiative exchanges with the heaters, the net heat flux on the design surface should be negative (Qd < 0). This sign convention is consistent with the relations in Eqs. (2) and (4). Equation (5) also relies on the fact that the areas of the heating and design surfaces are the same, that is, Dx2. Based on Eq. (5), it follows that the net heat flux on the filament heater is known a priori. The problem will be solved according to the implicit formulation. First, it is proposed a certain geometric configuration for the filament heater. Next, it is specified the dimensionless temperature on the design surface, tjd ¼ Tjd =Td ¼ 1. For the heaters and the adiabatic wall, the net radiative fluxes are specified, Qjh ¼ Qh and Qjw ¼ 0, respectively. After the specification of the thermal boundary conditions, Eq. (3) is applied for each design surface (since the temperature is specified), while Eq. (4) is applied for each heating or wall element (for the heat flux is specified), leading to a well-conditioned system of linear equations on the unknown radiosities. After solving for all radiosities, Eq. (4) can be applied for each design surface element to determine the net radiative heat flux, which is then compared to the specified value, Qd. Based on the deviations between the calculated heat flux, Qr;jd , and the specified value, Qd, another geometrical configuration is proposed, and the process is rerun until a satisfactory solution is achieved. In the present study, this calculation process will be formulated as an optimization problem to search for the best geometrical configuration of the filament. Mathematically, this optimization problem can be formulated by means of an objective function that measures the deviation between the specified and the calculated heat flux on the design surface, Qd and Qr;jd , respectively. The objective functions is thus defined as

vffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi u Jd
 uX
2 F ix;1 ; iy;1 ; :::; ix;jh ; iy;jh ; :::; ix;Jh ; iy;Jh ¼ t Qd  Qr;jd

equation:

2m 

vun  vln þ1 q

where vun and vln are the upper and lower bounds, respectively, of variable vn, with n ¼ 1; N. In the present problem, vn represent the values of the integers ix;jh or iy;jh that define the position of the heating element jh. Fig. 3 presents a string with 2 bits, m ¼ 2, for all variables ix;jh and iy;jh , except for ix,1 and iy,1, which are decoded with 4 bits, m ¼ 4. The reason for this choice will be explained in the next section. The total number of bits of the string will be denoted by Lb. In Fig. 3, the number of bits of the string will be therefore Lb ¼ 2  (N þ 2). Conversely, the physical value of each design variable is given by:

  vn ¼ vln þ vun  vln

(6)

where ix;jh and iy;jh are integer numbers that define the spatial position of element jh that forms the filament, as depicted in Fig. 2(b). Considering that the upper surface is divided in Ix and Iy elements in the horizontal and vertical directions in the figure, it follows that 1  ix;jh  Ix and 1  iy;jh  Iy. Therefore, the number of unknown variables is N ¼ 2  Jh, where Jh is the number of heating elements that form the filament, as defined before. The next section presents the methodology to determine the spatial distribution of the heating elements (that is, the values of ix;jh and iy;jh for each heating element jh) that minimizes the objective function defined in Eq. (6). 3. Searching the optimum position of the heating filament 3.1. The generalized extremal optimization method (GEO) Both deterministic and stochastic methods can be applied in the minimization of Eq. (6). One basic characteristic of stochastic methods is the incorporation of some form of randomness to the searching process as an aid to avoid local minima. In this study, the generalized extremal optimization (GEO) algorithm is applied, since it proved effective to solve inverse design problems in Refs. [25,26], but other solutions techniques could be attempted. In the GEO algorithm [29], the N variables of the optimization problem are encoded by only one binary string, and the bits are randomly aligned to form an initial bit configuration, as depicted in Fig. 3. As just seen, N ¼ 2  Jh. The number of bits m to attain a desired precision q for each design variable is given by the following

In 2m  1

(8)

where In is an integer number obtained in the transformation of the variable vn from its binary form to a decimal representation. In the GEO algorithm, one bit is flipped (from 0 to 1 or viceversa) and the objective function value, given by Eq. (6), is calculated and compared to the previous one, then the bit is returned to its original value. This process is repeated until all bits in the string are flipped. Next, the bits are ranked according to the change in the objective function. In minimization problems, the lowest rank (k ¼ 1) is assigned to the bit that leads to the smallest objective function when flipped, while the highest one (k ¼ Lb) is assigned to the bit that leads to the highest objective function. Then, a single bit in the string is randomly selected, and the probability of its mutation (from 0 to 1 or vice versa) is computed according to the following relation:

PðkÞ ¼ kt ;

jd ¼1

(7)

1  k  Lb

(9)

where t is a positive adjustable parameter. Finally, a random number RAN is generated. If RAN  P(k), then the bit is confirmed to mutate. Otherwise, another candidate bit is randomly picked, and the process is repeated until a bit is confirmed to mutate. Equation (9) establishes that the bits with bad fitness (low k) are more likely to mutate than the bits with good fitness (high k). However, even the bits with the best fitness can mutate. This feature of the method permits the solution to escape from local minima.

3.2. Assembling and evolution of the filament In the application of the GEO algorithm, it is necessary to start with a first choice of the unknown parameters to be sought, that is, the values of positions ix;jh and iy;jh of each heating element that forms the filament. This is equivalent to initialize the problem with a random geometrical configuration of the filament. In the present approach, the first element of the filament, defined by (ix,1,iy,1), is randomly chosen to fit in any given position on the top surface, with the constraint that 1  ix;1  Ix and 1  iy;1  Iy . For the particular grid considered in this study, Ix ¼ 15 and Iy ¼ 12. Using relation Eq. (7) and setting the precision q to be 1 (since ix,1 and iy,1 are integers) leads to the 4 bits in Fig. 3 that encodes those variables. After a random generation of bits for these two variables, the first element is placed in a certain position in the grid as shown in Fig. 4(a). The next element (ix,2,iy,2), is again randomly chosen, but within the constraints that ix;1  1  ix;2  ix;1 þ 1 and 1  ix;2  Ix , and iy;1  1  iy;2  iy;1 þ 1 and 1  iy;2  Iy . Fig. 4(b) shows the possible domain where the second element can be randomly placed. Using relation Eq. (7) with precision q equal to 1, each

L.D. Lemos et al. / International Journal of Thermal Sciences 105 (2016) 1e12

5

Fig. 3. String of bits encoding the variables to be optimized, that is, the x and y positions of each heating element of the filament, given by the integers ix;jh and iy;jh .

Fig. 4. (a) The first element of the filament in the grid. (b) The possible locations for the second element (shaded and dark red areas, including the first element element). (c) The second element and the possible locations for the third element (shaded and dark red areas). (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)

variable can be encoded by only 2 bits or four possible configurations: j0j0j corresponds to ix,2 ¼ ix,1-1, j1j0j and j0j1j to ix,2 ¼ ix,1, and j1j1j to ix,2 ¼ ix,1 þ1. The next figure, Fig. 4(c), shows one possible location of the second element, corresponding to ix,2 ¼ ix,1 þ1 (encoded by j1j1j) and iy,2 ¼ iy,1 þ1 (again, j1j1j), and the domain for the random placement of the third element. In this particular example, the second element sits on the border of the symmetry line, so the possible locations for the third element is adjusted accordingly to avoid it to be placed outside the domain. The process is repeated until all the elements of the filament are placed on the top surface. This is shown in Fig. 5(a) for a filament with 20 elements. Once the filament is assembled, the next step in the solution is to let the filament to evolve according the GEO algorithm. As described in Section 3.1, the binary bits that form the string are flipped one by one to perturb the variables towards the optimum solution. To demonstrate how this is accomplished in the present methodology, consider the 7th element of the filament, which is marked in Fig. 5(a). In this example, the horizontal position of the element corresponds to ix,7 ¼ ix,6, represented by j0j1j. When the second bit (the one on the left) that encodes x7 is flipped from 0 to 1, the new position will be ix,7 ¼ ix,6 þ 1, represented by j1j1j. The new position of the seventh element is shown in Fig. 5(b). All the other

bits of the string keep their original values. This means that all the elements of the filament before the 7th element keep their original positions. The vertical position of the 7th element, y7, also retains the same value as the one before this step. The bits of the elements after the 7th element also retain their values, so the relative positions of these elements with respect to the 7th element remain the same. However, since the position of the 7th element is shifted to the right in this step, all the following elements are shifted accordingly. The only exception is the 13th element, which cannot be shifted to the right since it sits in the limit of the domain. In this case, the restriction ix,13  Ix is applied, and the element remains in the same position. This is also shown in Fig. 5(b). According to the GEO algorithm, the heat flux on the design surface is computed for this new configuration of the filament, the objective function is computed according to the methodology outlined in Section 2, and the bit is ranked according to the variation in the objective function. The bit is flipped back to its original value, and so the 7th element and all the following elements return to their previous positions. Then, the next bit of the string is flipped and the process is continued until all bits are ranked. Finally, one of the bits is mutated according to the methodology described in Section 3.1, which leads to a new configuration for the filament heater. The procedure is repeated until the objective function falls below a given tolerance.

Fig. 5. (a) A possible initial configuration of a filament with 20 heating elements (the “seed”). (b) New configuration of the filament after flipping one of the bits that encode the 7th heating element.

6

L.D. Lemos et al. / International Journal of Thermal Sciences 105 (2016) 1e12

One of the advantages of the proposed methodology is the simplicity of implementation. In fact, the standard GEO algorithm is used without any modification in its original form. After flipping each bit, the filament is only slightly modified, which allows a smooth evolution towards the optimum solution. The filament in Fig. 5(b) could be one possible mutation of the configuration in Fig. 5(a) if the bit discussed above was the one chosen to mutate. It should be observed that there is no restriction in the proposed methodology to avoid that two or more elements of the filament occupy the same position (ix, iy) in the grid. In that case, the heat flux is multiplied by the number of heating elements in the position. However, as will be seen in the section of results, none of the optimum solutions presented two or more elements occupying the same position in the grid, an indication that the evolution towards the optimum solution avoids such condition. 4. Results and discussion 4.1. Filaments formed with square elements In this study, the emissivities of the design surface and filament heater are εjd ¼ 0:9 and εjh ¼ 0:9, respectively. The remaining surfaces of the enclosure are adiabatic, so their emissivities do not affect the solution. The aspect ratio of the enclosure base is W/ L ¼ 0.8; the dimensionless height is H/L ¼ 0.2. Due to the geometry symmetry, only one quarter of the enclosure needs to be solved. The adopted number of elements in the x, y and z directions are Ix ¼ 15, Iy ¼ 12, and Iz ¼ 6, respectively. Grid independence tests revealed that the chosen number of elements provides sufficient numerical accuracy in face of the errors that are expected from the inverse analysis. The design surface dimensions, as seen in Fig. 1(a), are taken as Wd/L ¼ 0.6 and Ld/L ¼ 0.8. In this mesh configuration, for the one-quarter domain of the enclosure, the number of elements on the design surface is Jd ¼ 108, as in Fig. 2(a). The specified conditions on the design surface elements are the net radiative heat flux, equal to qd ¼ 3.22  103 W/m2, and the temperature, equal to Td ¼ 673 K, which are representative values found in industrial heating processes. The negative value for the radiative flux follows from the convention that heat out of a surface is positive. In dimensionless form, the specified conditions on the design surface are therefore Qr;jd ¼ Qd ¼ 0:277 and tjd ¼ td ¼ 1:0, respectively. For the adiabatic walls, Qr,jw ¼ 0, while for the filament heater the heat flux is determined from the global energy balance in the enclosure, as given by Eq. (5), and depends on the number of elements that form the filament, Jh. In this study, filaments with Jh ¼ 20, 30, 40 and 50 heating elements are considered, leading to dimensionless heat fluxes of Qr;jh ¼ Qh ¼ 1.50, 1.00, 0.75 and 0.60, respectively. The first case considers that the filament is composed of 20 elements, that is, Jd ¼ 40 and Qr;jh ¼ Qh ¼ 1.50. The methodology described in Section 3.2 is used to make the first assembling of the filament, and then to make it to evolve according to the GEO algorithm. In the application of the algorithm, two parameters are necessary to be chosen: the number of mutations, and the parameter t, which governs the probability of a given bit to mutate according to Eq. (9). The first study considered the number of mutations, which was carried out for a value of t, chosen by experience, equal to 1.0. As described in Sections 3.1 and 3.2, prior to each mutation it is necessary the evaluation of the objective function as many times as the number of bits that forms the string in Fig. 3, which in this case is Lb ¼ 2  (N þ 2) ¼ 84 bits (N ¼ 2  Jd ¼ 40 variables). Fig. 6(a) shows the evolution of the objective function with the number of mutations. As expected, the value of the objective function decreased with the mutations, stabilizing after about 500 mutations or 500  84 ¼ 42,000

evaluations of the objective function. The next step consists of determining the optimum value of the parameter t. Fig. 6(b) shows the behavior of the objective function after 1000 mutations of the objective functions for values of parameter t between 0.50 and 2.5. As seen, the optimum value of parameter t was 1.0, which was adopted for the next calculations. In the proposed methodology, the first assembling of the filament (the “seed” of the solution) is a random process, so one question is whether the final form of the filament is dependent of the seed. The tests carried out in this study showed that this is the case. To quantify the effect of the seed, Table 1 shows the values of the objective function for different seeds after 2500 mutations, corresponding to 2500  84 ¼ 210,000 evaluations of the objective function, and with t ¼ 1.0. As expected from the inspection of Fig. 6(a), increasing the number of mutations did not further decrease the value of the objective function. In the table, the seeds are classified from 1 to 10 according to the value of the objective function, from best to worse. As seen, the objective function after 2500 mutations varied between 0.537 and 0.767 for the different seeds. Fig. 7(a) shows Seed 1 and the final filament configuration that led to the objective function with the smallest value. As seen, the overall shape of the final configuration roughly resembled the seed, although between the two configurations there is a substantial decrease in the objective function from 1.30 to 0.537. Fig. 7(b) presents the complete view of the top surface with four filaments, one in each quarter of the surface. The design surface on the bottom of the enclosure is also presented in the figure for a better visualization of the positions of the filaments in relation to the design surface. The definite way to evaluate the inverse solution is simply computing the radiative heat flux on the design surface for the final filament configuration. Fig. 8(a) shows the dimensionless heat flux on the design surface for the filament configuration shown in Fig. 7(b). As seen, the heat flux ranged from 0.282 to 0.266, thus comprehending the target value of Qd ¼ 0.277, but lacked the required uniformity. The accuracy of the solution can be quantified by computing the deviation between the specified heat flux and the heat flux obtained from the final configuration of the filament heater:

  Q  Q   d r; jd  gjd ¼     Qd

(10)

Based on the calculation of the deviation for each element jd that forms the design surface, the arithmetic mean deviation was found to be 1.81%. However, the maximum deviation was relatively high, reaching a value of 9.81%. One possible way to improve the design is to increase the number elements that form the filament heater, for instance, to 30 elements. In this case, Jd ¼ 60 and the binary string is formed with Lb ¼ 124 bits. As informed above, in this case the dimensionless heat flux on the heater should be Qr;jh ¼ Qh ¼ 1.00 to provide all the heat that is required on the design surface. The procedure to obtain the final configuration of the heating element was the same as the one described for the filament with 20 elements. First, the optimum parameter t was found by computing the objective function for 1000 mutations. As with the 20-filament, the optimum value was t ¼ 1.0. Next, ten seeds were randomly generated to select the one that could lead to the smallest value of the objective function after 2500 mutations or 2500  124 ¼ 310,000 evaluations of the objective function. The corresponding values of the objective functions are shown in Table 1. Overall, the objective function decreased in comparison with the filament with 20 elements, varying between 0.218 and 0.645. The solution for the seed that led to the smallest objective function, F ¼ 0.218, is shown in Fig. 9(a) for one quarter of the domain, and in Fig. 9(b) for the entire view of the

L.D. Lemos et al. / International Journal of Thermal Sciences 105 (2016) 1e12

7

Fig. 6. (a) Variation of the objective function with the number of mutations (t ¼ 1.0). (b) Variation of the objective function with t (after 1000 mutations). Results for the filament with 20 heating elements.

Table 1 Values of the objective function after 2500 mutations for different seeds and for filaments with different number of heating elements. Seed

20 Heating elements

30 Heating elements

40 Heating elements

50 Heating elements

1 2 3 4 5 6 7 8 9 10

0.5374 0.5683 0.5700 0.5783 0.5802 0.5975 0.5975 0.6065 0.7002 0.7672

0.2176 0.2489 0.2902 0.3286 0.3311 0.3373 0.3592 0.6302 0.6332 0.6450

0.1904 0.2105 0.2147 0.2161 0.2162 0.2203 0.2246 0.2347 0.2461 0.3361

0.1186 0.1408 0.1504 0.1795 0.1811 0.1820 0.1861 0.2497 0.2752 0.3433

Fig. 7. (a) The first configuration of the filament, the “seed”, in gray, and the optimum configuration of the filament with 20 heating elements, in y; (b) The full view of the top surface with four filaments. The shaded area represents the design surface.

top surface. It is interesting to observe that in this case one of the ends of the filament touches the symmetry line, which means that it connects with the other filament. As seen in Fig. 9(b), only two filaments would be required in this design solution. The heat flux distribution for this filament is shown in Fig. 8(b). It can be observed a considerable improvement in the uniformity of the heat flux in comparison to the filament with 20 elements. The mean and the maximum relative errors were reduced to 0.730% and 4.62%,

respectively. The final two solutions involved filaments with 40 and 50 elements (Jd ¼ 80 and Jd ¼ 100). In these cases, the required dimensionless heat fluxes on the filaments were Qr;jh ¼ Qh ¼ 0.75 and 0.60, respectively. As with the filaments with 20 and 30 elements, the optimum value of parameter t was 1.0. Table 1 presents the value of the objective function for different seeds after 2500 mutations, corresponding to 410,000 and 510,000 evaluations of the

8

L.D. Lemos et al. / International Journal of Thermal Sciences 105 (2016) 1e12

Fig. 8. Radiative heat flux on the design surface for filaments with (a) 20, (b) 30, (c) 40 and (d) 50 heating elements.

Fig. 9. (a) Optimum configuration of the filament with 30 heating elements; (b) The full view of the top surface with two filaments.

objective function for the filaments with 40 and 50 elements, respectively. As seen in the table, the overall effect of increasing the number of elements was the reduction in the objective function. For the seeds that led to the minimum objective functions, the optimum filaments with 40 and 50 elements are shown in Figs. 10 and 11, respectively. It is interesting to observe that since both ends of the filament with 50 elements touch the symmetry lines, this design would require only one filament. Although its shape is not trivial, it is still quite possible to be manufactured. The heat fluxes are shown in Fig. 8(c) and (d). The mean and the maximum

deviations of the filament with 40 elements were 0.630% and 2.86%. For the filament with 50 elements the errors were 0.398% and 2.43%. Those deviations can be considered quite satisfactory, especially if considered that in this type of inverse problem no exact solution should be expected. Table 2 compares the errors for the different number of heating elements, showing a decrease of the error with the increase in the length of the filament. This trend, however, is expected to invert at some point. In fact, if a very long filament was chosen to fill the entire top surface of the enclosure, the condition would approach uniform heat flux input on the top

L.D. Lemos et al. / International Journal of Thermal Sciences 105 (2016) 1e12

9

Fig. 10. (a) Optimum configuration of the filament with 40 heating elements; (b) The full view of the top surface with four filaments.

Fig. 11. (a) Optimum configuration of the filament with 50 heating elements; (b) The full view of the top surface with a single filament.

Table 2 Average and maximum relative errors in the radiative heat flux on the design surface for filaments with different number of heating elements. Number of heating elements

Average relative error (%)

Maximum relative error (%)

20 30 40 50

1.81 0.730 0.630 0.398

9.81 4.62 2.86 2.43

surface, which is known to lead to non-uniform heating on the bottom surface of the enclosure. The filaments shown in Figs. 7 and 9e11 were the final configurations for the seeds that led to smallest objective functions (Seed 1). The different seeds in Table 1 result in different filament configurations even for a fixed number of heating elements. Although the final solution depends on the initial seed, its accuracy could still be satisfactory after a sufficient number of evolutions. Taking the example of the filament with 40 elements, starting with Seed 2, the maximum and average errors of the inverse solution were 3.99% and 0.71%, respectively; for Seed 3, the errors were 3.61% and 0.72%. Thus, both solutions can still be considered satisfactory. The final geometric configurations starting with these two seeds are presented in Fig. 12(a) and (b). As seen, they are considerably different from each other and from the configuration in Fig. 10 (starting with Seed 1). This illustrates one basic aspect of inverse design: there may be different solutions that approximately satisfy a given

problem. This is another advantage of the method, that is, the designer can make a choice among a set of sufficiently accurate solutions considering other criteria, such as cost and simplicity of manufacturing. 4.2. The slender filament The previous section presented designs of filaments formed with square heating elements. In those figures, slender filaments were drawn over the square elements to represent what could be more realistic filament configurations. Using square elements that were coincident with the mesh highly simplified the search for the optimum configuration, since their locations could be specified by only two integers, ix,1 and iy,1. As previously observed, the dimensions of the enclosure and of the square elements were presented in dimensionless form, which means that, if the characteristic dimension of the enclosure is large, the width of each

10

L.D. Lemos et al. / International Journal of Thermal Sciences 105 (2016) 1e12

Fig. 12. Final configurations of the filament with 40 heating elements starting with (a) Seed 2; and (b) Seed 3 in Table 1.

square element could be several times larger than those of typical electrical filaments or heating tubes. The question that arises is what would be the error introduced if those slender filaments were used to replace the filaments formed with square elements in the post-processing stage. To study this problem, one can consider a section of the slender filament placed over a square element with common center, as seen in Fig. 13(a). The dimensionless width of the element is r ¼ Ds/Dx, so that r  1. The major difference in considering the slender filament regards the view factor with the other elements of the enclosure. To quantify this effect, the view factors between the heating filaments (both the square and the slender section) and one element on the design can be computed using the formula for the configuration shown in Fig. 14, which is available in Ref. [30]. Fig. 13(b) shows a few possible locations of the heating element with relation to one element on the design surface, which is placed in the same x-y position of element 1, but at a distance in z of 6  Dx, the height of the enclosure under study. Table 3 presents the view factors of the square element and the section of the slender element to the design surface element. In the present test case, it was considered that the width of the slender filament was 1/10 of the square element, that is, r ¼ 0.1 As seen in Table 3, the view factors were in general very close to each other, with a maximum relative deviation less than

Fig. 14. Geometrical configuration for the computation of the view factor: rectangle to rectangle in a parallel plane [30].

Fig. 13. (a) Section of the slender filament placed over the square element that forms the computational grid; (b) Possible locations of the heating elements relative to the design surface element, placed in position 1 but at a distance of 6  Dx in the z direction.

L.D. Lemos et al. / International Journal of Thermal Sciences 105 (2016) 1e12

11

Table 3 Comparison of the view factors of the square element (jsq) and the slender heating element (jsl) to the design surface element (jd). The design surface element is placed in the same x-y position of 1 in Fig. 12(b), but at a distance z of 6  Dx. Position of the heating element in Fig. 12(b)

Fjsq jd  103

Fjsl jd  103

Relative deviation (%)

1 2 3 4 5 6 7 8 9

8.682 8.235 7.079 8.235 7.820 6.747 7.079 6.747 5.878

8.721 8.265 7.091 8.271 7.849 6.758 7.108 6.770 5.888

0.4483 0.3682 0.1690 0.4370 0.3608 0.1703 0.4062 0.3400 0.1726

Fig. 15. Heat flux on the design surface: (a) for the filament formed with 50 square heating elements; (b) for the equivalent slender filament.

0.5%. The reason for this proximity is that the perfectly diffuse heating elements are sufficiently distant from the design surface for the differential approximation to be applied. In this approximation, the view factor depends almost solely of the distance between the centers of the elements and incident angles, not the shape of the elements. The heat flux on the design surface was next calculated using the slender filament. Following the results in Table 3, the view factors between the slender filaments and the other elements of the enclosure were taken as the same as the view factors between the square heating elements and those other elements. In this calculation, it was considered that the remaining part of the square element is insulated, as seen in Fig. 13(a). This test was carried out for the heating element formed with 50 elements, as seen in Fig. 11(a) and (b). Keeping r ¼ 0.1, the area of the slender filament was reduced to 12.0% of the original filament (with square elements), so the required dimensionless heat flux on the slender filament was increased to Qr;jh ¼ Qh ¼ 5.0. The heat flux is presented in Fig. 15 together with the heat flux obtained for the filaments made with the square elements, already shown in Fig. 8(d). As seen, the heat fluxes were very close, showing that slender filaments could be designed in the present case by simply placing them over the square elements, whose optimum locations are found from the inverse solution proposed in this study.

design surface. The inverse design was formulated as an optimization problem for the minimization of an objective function that measured the deviation between the specified heat flux and the one computed for a given configuration of the heating elements. The optimization problem was solved with the generalized extreme optimization (GEO) algorithm, a stochastic method. The methodology was applied to solve for filaments with different sizes, ranging from 20 to 50 elements. Overall, increasing the number of elements in the filaments led to more uniformity in the heat flux on the design surface. For the filament with 50 heating elements, the average and the maximum errors in the radiative heat flux on the design surface were about 0.4% and 2.4%, respectively, which can be considered fully satisfactory for this type of problem. In the present solution, the elements of the filament were coincident with the square elements that formed the grid. This greatly simplified the assembling and evolution of the filament, but could lead to nonrealistically wide filaments depending on the characteristic dimension of the enclosure. In the last part of the paper, it was shown that the filament formed with square elements could be replaced with a more realistic slender filament in the postprocessing stage. The main requirement is that the dimensions of the surface elements should be sufficiently small to permit the use of the differential approximation in the computation of the respective view factors.

5. Conclusions

Acknowledgments

This paper presented the application of the inverse analysis to determine the geometrical configuration of filament heaters to achieve uniform temperature and net radiative heat flux on the

This research conducted with support from National Center for Supercomputing (CESUP), at University Federal of Rio Grande do Sul (UFRGS). Author RB also thanks CNPq (Brazil) for research grant

12

L.D. Lemos et al. / International Journal of Thermal Sciences 105 (2016) 1e12

150108/2015-0. Author FHRF thanks CNPq (Brazil) for research grants 309961/2013-0 and 476490/2013-8. References [1] Hansen PC. Rank-deficient and Discrete Ill-posed Problems: Numerical Aspects of Linear Inversion. Philadelphia: SIAM Monographs on Mathematical and Computation; 1998. [2] Ertürk H, Ezekoye OA, Howell JR. Comparison of three regularized solution techniques in three-dimensional inverse radiation problem. J. Quant. Spectrosc. Radiat. Transf. 2000;73:307e16. [3] França F, Ezekoye O, Howell J. Inverse boundary design combining radiation and convection heat transfer. J. Heat Transf. 2001;123:884e91. [4] Sarvari SMH, Mansouri SH, Howell JR. Inverse design of three-dimensional enclosures with transparent and absorbing-emitting media using an optimization technique. Int. Commun. Heat Mass Transf. 2003;30:149e62. [5] Daun K, Howell JR, Morton DP. Optimization of transient settings to provide spatially uniform heating in manufacturing processes involving radiant heating. Numer. Heat. Transf. Part A 2004;46:651e67. [6] Sarvari SMH. Inverse determination of heat source distribution in conductiveradiative media with irregular geometry. J. Quant. Spectrosc. Radiat. Transf. 2005;93:383e95. [7] Daun KJ, Howell JR. Inverse design methods for radiative transfer systems. J. Quant. Spectrosc. Radiat. Transf. 2005;93:43e60. [8] Kowsary F, Pooladvand K, Poushaghaghy A. Regularized variable metric method versus the conjugate gradient method in solution of radiative boundary design problem. J. Quant. Spectrosc. Radiat. Transf. 2007;108: 277e94. [9] Rukolaine SA. Regularization of inverse boundary design radiative heat transfer problems. J. Quant. Spectrosc. Radiat. Transf. 2007;104:171e95. [10] Mossi AC, Vielmo HA, França FHR, Howell JR. Inverse design involving combined radiative and turbulent convective heat transfer. Int. J. Heat Mass Transf. 2008;51:3217e26. [11] Bayat N, Mehraban S, Sarvari SMH. Inverse boundary design of a radiant furnace with diffuse-spectral design surface. Int. Commun. Heat Mass Transf. 2010;37:103e10. [12] Hoffmann RS, Seewald A, Schneider PS, França FHR. Inverse design of thermal systems with spectrally dependent emissivities. Int. J. Heat Mass Transf. 2010;55:931e9. [13] Dehghani M, Hosseini Sarvari SM, Ajam H. Inverse estimation of boundary conditions on radiant enclosures by temperature measurement on a solid object. Int. Commun. Heat Mass Transf. 2011;38:1455e62. [14] Chopade RP, Mishra SC, Mahanta P, Maruyama S, Komiya A. Uniform thermal conditions on 3-D object: optimal power estimation of panel heaters in a 3-D

radiante enclosure. Int. J. Therm. Sci. 2012;51:63e76. [15] Chopade RP, Mishra SC, Mahanta P, Maruyama S. Estimation of power of heaters in a radiant furnace for uniform thermal conditions on 3-D irregular shaped objects. Int. J. Heat Mass Transf. 2012;55:4340e51. [16] Mosavati B, Mosavati M, Kowsari F. Solution of radiative inverse boundary design problem in a combined radiating-free convecting furnace. Int. Commun. Heat Mass Transf. 2013;45:130e6. [17] Moghadassian B, Kowsary F. Inverse boundary design problem of natural convection-radiation in a square enclosure. Int. J. Therm. Sci. 2014;75: 116e26. [18] Sarvari SMH. Optimum placement of heaters in a radiant furnace using the genetic algorithm. In: 13th International Heat Transfer Conference, IHTC-13; 2006 (Sydney, Australia). [19] Ashouri J, Sarvari SMH, Farahat S. Optimum arrangement of heaters in a threedimensional radiant furnace using the genetic algorithm. In: 7th International Conference on Heat Transfer, Fluid Mechanics and Thermodynamics; 2010 (Antalya, Turkey). [20] Farahmand A, Payan S, Hosseini Sarvari SM. Geometric optimization of radiative enclosures using PSO algorithm. Int. J. Therm. Sci. 2012;60:61e9. [21] Gossard D, Lartigue B. Three-dimensional conjugate heat transfer in partitioned enclosures: determination of geometrical and thermal properties by an inverse method. Appl. Therm. Eng. 2013;54:549e58. [22] Rukolaine SA. Shape optimization of radiant enclosures with specular-diffuse surfaces by means of a random search and gradient minimization. J. Quant. Spectrosc. Radiat. Transf. 2015;151:174e91. [23] Pourshaghaghy A, Pooladvand K, Kowasary F, Karimi-Zand K. An inverse radiation boundary design problem for an enclosure filled with an emitting, absorbing and scattering media. Int. Commun. Heat Mass Transf. 2006;33: 381e90. [24] Safavinejad A, Mansouri SH, Sakurai A, Maruyama S. Optimal number and location of heaters in 2-D radiant enclosures composed of specular and diffuse surfaces using micro-genetic algorithm. Appl. Therm. Eng. 2009;29:1075e85. [25] Cassol F, Schneider PS, França FHR, Sousa FL, Silva Neto AJ. Multi-objective optimization as a new approach to illumination design of interior spaces. Build. Environ. 2011;46:331e8. [26] Brittes R, França FHR. A hybrid inverse method for the thermal design of radiative heating systems. Int. J. Heat Mass Transf. 2013;57:48e57. [27] Howell JR, Menguc MP, Siegel R. Thermal Radiation Heat Transfer. sixth ed. CRC Press; 2016. [28] Tucker RJ. Direct exchange areas for calculating radiation transfer in rectangular furnaces. J. Heat Transf. 1986;108:707e10. [29] Sousa FL, Ramos FM, Paglione P, Girardi RM. New stochastic algorithm for design optimization. AIAA J. 2003;41:1808e18. [30] Howell JR. A Catalog of Radiation Heat Transfer Configuration Factors. New York: McGraw-Hill; 1982.