Inverse Geometry Design of Radiative Enclosures Using Particle Swarm Optimization Algorithms

Abstract

Three different Particle Swarm Optimization (PSO) algorithms—standard PSO, stochastic PSO (SPSO) and differential evolution PSO (DEPSO)—are applied to solve the inverse geometry design problems of radiative enclosures. The design purpose is to satisfy a uniform distribution of radiative heat flux on the designed surface. The design surface is discretized into a series of control points, the PSO algorithms are used to optimize the locations of these points and the Akima cubic interpolation is utilized to approximate the changing boundary shape. The retrieval results show that PSO algorithms can be successfully applied to solve inverse geometry design problems and SPSO achieves the best performance on computational time. The influences of the number of control points and the radiative properties of the media on the retrieval geometry design results are also investigated.

1. Introduction

Radiative heating devices are encountered in various industrial fields, such as industrial boilers, spacecraft, infrared reflecting ovens, metallurgical equipment, and so on and the design of radiative enclosure has a direct impact on security issues [1]. Inverse design technique is a new method in recent years, whose solving process is to establish an objective function according to the design requirements at first, and then optimize the objective function by some optimization methods, and achieve the pre-specified purpose finally. Inverse design technique has the advantages of simple process, short design circle, good optimization results, etc., and it has got more and more attentions and applications.

Inverse design problems can be divided into two categories according to the design prerequisites [2]. One is inverse boundary design problems, in which the geometry shape of the radiative enclosures is fixed and the boundary conditions are need to be deigned [3–6]. The other is inverse geometry design problems, where the boundary conditions are predetermined and the geometry shape of the design surface needs to be designed [7, 8]. For the reason that the geometry shapes of the boundaries of radiative enclosures are different at each iteration and the grids in the computational domain must be re-meshed in every iterative calculation, the inverse geometry design of radiative enclosures is the most complex inverse radiative problem [2, 9].

During the last few decades, some inverse design techniques have been successfully used for solving inverse geometry design problems. Howell et al. [10, 11] proposed inverse design ideas and applied inverse Monte Carlo techniques, Tikhonov method, truncated singular value decomposition (TSVD), Modified TSVD (MTSVD), artificial neural networks (ANN) and conjugate gradient method (CGM) to solve the inverse design problem of a three-dimensional industrial furnace, which greatly improved the practical design of thermal and environmental systems. Franca and Howell [12] studied a transient inverse design problem that finding the optimal location of a heater on the top surface of a three-dimensional enclosure to produce a prescribed time-dependent temperature distribution on the bottom surface of the enclosure, the TSVD method is used to regularize the ill-conditioned system of linear equations the pre-specified temperature curve is obtained with an error of less than 1.0%. Tan et al. [13] applied meshless method to solve the coupled conductive and radiative heat transfer problem in heating devices, in which a series of nodes are used for discretizing the computational domain to overcome the tedious re-mesh works, and CGM is adopted to optimize the height of the adiabatic diffuse reflection surfaces and the geometry shape of the heating surface to satisfy the required total heat flux on the pre-appointed region of the low temperature heated surface. Sarvari and Mansouri [14] used CGM to minimize the objective function which is expressed as the sum of square residuals between estimated and desired heat fluxes on the design surface to satisfy the specified temperature and heat flux distributions. The radiative heat transfer problem in the two-dimensional irregular enclosure filled with gray participating media with uniform absorption coefficient is solved by discrete transfer method, the effect of optical depth and angular refinement on the inverse design results are also investigated.

However, it can be found that most of the above researches about inverse geometry design problems are solved based on the gradient-based methods. All these methods have the common disadvantages that the computational process of getting the gradient is complicated and the retrieval results are strongly depended on the initial guessed values. On the contrast, intelligent algorithms can abolish the complex calculations about gradient and randomly generate the potential solutions in the search space to overcome the drawback of depending on the initial value. In recent years, some intelligent algorithms have been successfully applied for solving inverse radiative problems, including PSO, genetic algorithm (GA), ant colony optimization (ACO), difference evolution (DE), fruit fly optimization algorithm (FOA), to name a few [15–22]. Compared with conventional techniques, intelligent algorithms can obtain much more potential solutions at each iteration and all searches are executed in parallel, which greatly improved the computational efficiency, especially for solving some high dimensional problems. Intelligent algorithms also have been successfully used to solve inverse design problems. For example, Moparthi et al. [23] solved the coupled radiative and conductive heat transfer problem in one-dimensional planar system based on finite volume method (FVM) and lattice Boltzmann method (LBM) and applied GA to optimize the heater temperature to produce the desired heat flux and temperature distributions on the design surface. The retrieval results show that the temperature or heat flux of the heater surface has a significant effect on the design surface condition, the medium properties and the distance between the two surfaces. Amiri et al. [24] adopted modified discrete ordinate method to solve the radiative transfer equation (RTE), the micro genetic algorithm (MGA) is employed to optimize the objective function which is defined as the sum of the square of the difference between estimated and desired heat fluxes on the design surfaces. The design purpose is finding the best number and locations of the discretized heaters to meet the desired temperature and heat flux distributions on the design surface. Sarvari et al. [25] discretized the design surface into a series of control points and used B-spline to approximate the geometry shape of the boundary, the MGA is employed to optimize the locations of the control points to produce a desired heat flux distribution on the temperature-specified surface. The effects of corresponding parameters on the inverse design results are also investigated and the angular meshes are recommended as $N_{\theta }\times N_{\phi }=10\times 10$ . In addition, the results indicate that optimizing the weights corresponding to the control points can improve the quality of designed shape. However, the MGA has the obvious drawback that the convergence velocity is relatively low.

The PSO algorithm is a kind of biologically inspired algorithm whose search process is similar to foraging of birds and it was proposed in 1995 by Eberhart and Kennedy [26]. The physical model of PSO is very simple and the computational program is easy to be implemented, it also has strong robustness and achieves good performance on computational efficiency and accuracy. In addition, PSO algorithm can well balance the global and local search of particles, which enhance the global convergence of the algorithm. Farahmand et al. [27] investigated the inverse geometry design of two-dimensional radiative enclosures with diffuse gray surfaces based on the PSO and the retrieval results show that PSO algorithm obtains better performance in satisfying the design goal based in terms of computational accuracy and CPU time compared with MGA. However, the standard PSO algorithm also suffers from easily trapping into local optima in solving high dimensional problems. In order to strengthen the applicability of PSO, some improvements have been proposed and widely applied, including Stochastic PSO (SPSO), Differential Evolution PSO (DEPSO), Multi Phase PSO (MPPSO), and so on. However, to the authors’ best knowledge, there are few reports concerning about the applications of improved PSO algorithms for solving inverse geometry design problems of radiative enclosures.

In this chapter, the application of PSO algorithms in solving inverse geometry design problems of two-dimensional radiative enclosures filled with participating media is investigated. The design goal is to satisfy a uniform distribution of radiative heat flux on the designed surface. The discrete ordinate method (DOM) with a body-fitted coordinate system is used to solve the RTE. The standard PSO, SPSO and DEPSO algorithms are applied to optimize the locations of the control points, and Akima cubic interpolation is adopted to obtain the boundary geometry shape through these points. A typical inverse geometry design test is studied to demonstrate the good performance of PSO algorithms and the effects of corresponding parameters are also discussed.

The remainder of this chapter is organized as follows: the theoretical principles of PSO algorithms are introduced in Section 2. The feasibility of PSO algorithms by four famous benchmark functions is verified in Section 3. The inverse geometry design of two-dimensional radiative enclosures and the influences of the number of control points and the radiative properties of media on the inverse design results are investigated in Section 4. The main conclusions of the researches in this chapter are summarized in Section 5.

2. Theoretical overview of PSO algorithms

Bird individuals will communicate with each other to share their information about food when they are foraging, which can help birds find food faster. The advantages of cooperation of bird swarm are much greater than the disadvantages of competition among bird individuals. Based on the features of bird foraging behavior, Kennedy and Eberhart proposed PSO algorithm in 1995 [26]. The solving process of PSO algorithm is similar to the foraging behavior of birds, and the corresponding relationships are shown in Table 1.

Bird foraging behavior	PSO algorithm
Foraging domain of bird individual	Searching space of each particle
Bird individual	Particle
Flight speed of bird	Moving speed of particle
Location of bird	Location of particle, which represents a solution; of optimization problems
Location of food	The best solution of optimization problems

Table 1.

Corresponding relationships between bird foraging and PSO algorithm.

There are two dominant parameters in the PSO algorithm, namely, the speed and the location of particles. The moving speed decides the direction and distance of the particles, and every location of particles can be considered as the potential a solution of optimization problems. PSO adopts a combination of local and global searches and shares the evolutionary information among particle individuals to find the optimal solution.

2.1. Basic PSO algorithm

At the beginning of the optimization of PSO, every location and velocity of particles are randomly generated. During each iteration, there will be two extreme values. One is best location that an individual particle found so far, which is called local best location. Another is the best location that the whole particle swarm found so far, which is called global best location. The velocity and location of each particle are stochastically accelerated according to these two extremes and the evolutionary formula can be expressed as follows [26]

${\mathbf{V}}_i\left(t+1\right)={\mathbf{V}}_i\left(t\right)+c_1\cdot r_1\cdot \left[{\mathbf{P}}_i\left(t\right)-{\mathbf{X}}_i\left(t\right)\right]+c_2\cdot r_2\cdot \left[{\mathbf{P}}_{\mathrm{g}}\left(t\right)-{\mathbf{X}}_i\left(t\right)\right]$E1

where V_i(t) and ${\mathbf{V}}_i\left(t+1\right)$represent the velocity of ith particle at iterations tand t+1, respectively, and V_i ∈ [−V_max, V_max]. $c_1$and $c_2$are two positive constants called acceleration coefficient. $r_1$and $r_2$are two uniform random values in the range of [0, 1]. P_i(t) and P_g(t) indicate the local best location and global best location, respectively. X_i(t) denotes the location of ith particle, which is depending on the search experience of ith particle and surrounding particles. The evolutionary formula of ith particle’s location is defined as [26]

${\mathbf{X}}_i\left(t+1\right)={\mathbf{X}}_i\left(t\right)+{\mathbf{V}}_i\left(t+1\right)$E2

Figure 1.

The evolution of particle’s location.

According to Eq. (1), we can find that the velocity of ith particle consists of three parts: the first part on the right side is the current velocity of ith particle, which can counterpoise the local search and the global search; the second part on the right side denotes the influence of the search memory of ith particle, which makes individual particle has the ability of global search; the third part on the right side indicates the influence of the cooperation among particles. The updating process of the ith particle’ location is shown in Figure 1.

Figure 2.

The flowchart of the basic PSO algorithm.

The main procedure of PSO algorithm for solving optimization problems can be carried out according to the following steps:

1. Step 1: Initialization: Randomly initialize the location and the velocities of particle of every particle in the searching space, input the number of particle swarm, the maximum iteration numbers t_max and the stop criterion ε, etc. Set the current iteration as t=1.

2. Step 2: Fitness evaluation: Evaluate each particle’s fitness according to its location and determine the local best location P_i(t) and global best location P_g(t).

3. Step 3: Updating: Update the velocity and the location of each particle according to Eqs. (1) and (2), respectively. Calculate the new objective function value of each particle and update the local and global best locations P_i(t) and P_g(t).

4. Step 4: Comparison: Compare the objective function value of each newly obtained particle with the corresponding values at the last iteration. If the new objective function value is better than the one in the last generation, then the new location and velocity of this particle is updated. Otherwise, the new location is abandoned.

5. Step 5: Repeating: Check whether one of the following two stop criterion is reached: (1) the objective function value is less than the value of ε, and (2) the iteration number reaches the maximum iteration number. If so, go to the next step; otherwise, go to Step 3.

6. Step 6: Update of iteration number: Update the number of iteration from tto t+1.

7. Step 7: Termination of iteration: Output the global best optima and its corresponding results of optimization problems and then stop the calculation.

The flowchart of the basic PSO algorithm is shown in Figure 2.

However, there are some obvious shortcomings in the basic PSO algorithm, such as slow convergence, easy to fall into local optimum, and even the velocities tend to be infinity in some occasions. Thus, many modifications have been proposed to overcome these drawbacks.

2.2. Standard PSO algorithm

There are two important capabilities in PSO algorithm, namely exploration and exploitation of particles. Exploration is the phenomenon that particles leave the original orbit and search for new space. Exploitation is the phenomenon that particles look for better locations along the original track. In order to better take advantage of these two search way, Shi and Eberhart put forward the standard PSO algorithm on the basis of basic PSO in 1998 [28], in which an inertia weight coefficient wis introduced to control the impact of the current velocity on the next velocity. The velocity formula of ith particle can be expressed as [28]

${\mathbf{V}}_i\left(t+1\right)=w\cdot {\mathbf{V}}_i\left(t\right)+c_1\cdot r_1\cdot \left[{\mathbf{P}}_i\left(t\right)-{\mathbf{X}}_i\left(t\right)\right]+c_2\cdot r_2\cdot \left[{\mathbf{P}}_{\mathrm{g}}\left(t\right)-{\mathbf{X}}_i\left(t\right)\right]$E3

where wis the inertia weight coefficient, which can directly affect the balance between the global and local exploration abilities. At the initial stage of the search process, a big inertia weight coefficient is recommended to improve the global exploration ability in the relatively large space; whereas the inertia weight should be reduced with the iteration number increases to strengthen the local exploitation ability. It is worth pointing that a linearly decreasing inertia weight coefficient can successfully prevent particles from oscillating near the global best location [29]. Therefore, the inertia weight coefficient can be defined as

$w=w_{\text{max}}-\frac{t}{t_{\text{max}}}\cdot \left(w_{\text{max}}-w_{\text{min}}\right)$E4

Comparing Eq. (1) with Eq. (3) we can find that basic PSO algorithm is a special circumstance of standard PSO algorithm that inertia weight is set as w=1. The searching efficiency is significantly improved by introducing inertia weight w. However, the proportional relation of particle velocities at every generation is not all the same, and the standard PSO algorithm can’t be successfully applied for solving some complicated optimization problems. For breaking through the limitation of PSO, many modified techniques have been developed out and widely applied in engineering fields.

2.3. Stochastic PSO algorithm

In order to overcome the drawback that PSO algorithm converges too early and make sure to reach the goal of global convergence, Zeng and Cui proposed SPSO algorithm in 2004 [30, 31], in which a stopped changing particle is utilized to improve the global searching ability of particle swarms.

In SPSO algorithm, the inertia weight coefficient is set as w=0. Hence, the velocity of ith particle at t+1 iteration is determined by three parameters of titeration, namely X_i(t), P_i(t) and P_g(t). The new velocity of ith particle can be expressed as [30]

${\mathbf{V}}_i\left(t+1\right)=c_1\cdot r_1\cdot \left[{\mathbf{P}}_i\left(t\right)-{\mathbf{X}}_i\left(t\right)\right]+c_2\cdot r_2\cdot \left[{\mathbf{P}}_{\mathrm{g}}\left(t\right)-{\mathbf{X}}_i\left(t\right)\right]$E5

According to Eq. (5), it can be found that the local searching ability of SPSO is increased compared with standard PSO. However, the global searching ability is reduced significantly. In order to further strengthen the global searches of SPSO, the algorithm randomly generates a particle in the searching space whose location is X_j(t + 1), and other particles’ locations are updated based on the Eq. (5). The whole amending process can be expressed by the following equations

$\begin{array}{l}{\mathbf{P}}_j={\mathbf{X}}_j\left(t+1\right)\\ {\mathbf{P}}_i=\left\{\begin{array}{cc}\hfill {\mathbf{P}}_i,\hfill & \hfill F\left({\mathbf{P}}_i\right)<F\left[{\mathbf{X}}_i\left(t+1\right)\right]\hfill \\ \hfill {\mathbf{X}}_i\left(t+1\right),\hfill & \hfill F\left({\mathbf{P}}_i\right)\ge F\left[{\mathbf{X}}_i\left(t+1\right)\right]\hfill \end{array}\right.\\ {\mathbf{P}}_g^{\text{'}}=argmin\left\{F\left({\mathbf{P}}_i\right)\left|i=1,\cdots \cdots M\right.\right\}\\ {\mathbf{P}}_g=argmin\left\{F\left({\mathbf{P}}_g^{\text{'}}\right),F\left({\mathbf{P}}_g\right)\right\}\end{array}$E6

After above updates, the following criterions are executed:

(1) If P_g = P_j, which demonstrates that the random location X_jis the global best location. In this situation the jth particle will not be updated on the basis of Eq. (5) and the algorithm will randomly generate a location X_jin the searching space at the next iteration. At the same time, the velocities and locations of other particles are updated according to Eqs. (5) and (2) after P_gand P_jare updated, respectively.

(2) If P_g ≠ P_jand P_ghas not been updated, which indicates that there is no better global best location found, compared with the last iteration, then all the particles’ the velocities and locations are updated based on the Eqs. (5) and (2), respectively.

(3) If P_g ≠ P_jand P_ghas been updated, which demonstrates that there is a location of kth particle (k ≠ j) meets the requirement ${\mathbf{X}}_k\left(t+1\right)={\mathbf{P}}_k={\mathbf{P}}_g$, namely the local best location of kth particle is the global best location and it is better than the global best location at the last iteration. In this situation, the kth particle stops evolving and it is used for storing global optima. Other particles’ velocities and locations of other particles are updated according to Eqs. (5) and (2) after P_gand P_jare updated, respectively.

Through the above analysis we can find that there is at least one particle’s location reaches the global best location at a particular iteration, which indicates that at least one particle is randomly generated at each iteration. Therefore, SPSO algorithm has been proved with strong global search ability.

2.4. Differential evolution PSO algorithm

Differential Evolution (DE) algorithm adopts simple differential operation among potential solutions to produce new candidate solution, which is a parallel, direct and stochastic searching technique. It was first proposed for solving Chebyshev polynomials and global optimization problems over continuous spaces by Storn and Price in 1995 [32]. Taking a cue from DE, the mutation operation is introduced into PSO algorithm to overcome the drawback of trapping in local optima, which is called DEPSO.

In DEPSO algorithm, the differential evolution operator is introduced to increase the diversity of particle swarms which is defined as

$\beta =\chi \left[{\mathbf{X}}_{r_1}\left(t\right)-{\mathbf{X}}_{r_2}\left(t\right)\right]$E7

where χis the differential mutation operator which controls the magnification of differential variation ${\mathbf{X}}_{r_1}\left(t\right)-{\mathbf{X}}_{r_2}\left(t\right)$and it is usually set as a constant in the interval [0, 2]. $r_1$and $r_2$are two integers which are randomly chosen in the interval [1, M] and they are not equal to the index i. Hence, the location evolution of ith particle can be expressed as

${\mathbf{X}}_i\left(t+1\right)={\mathbf{X}}_i\left(t\right)+{\mathbf{V}}_i\left(t\right)+\beta \left(C-F_{\text{min}}\right)$E8

where Crepresents a preset constant value which satisfies C ≤ F_min, F_min indicates the minimum objective function value at the current iteration. The differential evolution term can force the particle to change if only C ≠ F_min, which can effectively prevent PSO algorithm from falling into local optima.

However, the location obtained by mutation operation maybe a worse result which will cause a bad influence on the search of other particles. In order to make sure the rapidity and stability of DEPSO algorithm, the following judgment should be executed before the location is updated:

(1) If the fitness value of the new location is better than the fitness value of the earlier location, which demonstrates the mutation is successful. Then the location of ith particle is updated according to Eq. (8).

(2) If the new fitness value is worse than before, which indicates this mutation is failed. Then the location of ith particle is updated according to Eq. (2) and the mutation operation of ith will continue at the next iteration until the mutation is successful.

In addition, there is a significant difference between DEPSO and basic PSO that the velocities of particles are not limited in the searching process, which can increase the convergence rate.

3. Simulation test of PSO algorithms

Function	Expression	Dimension	Search Space
Sphere	$f_1\left(x\right)={\displaystyle \sum }_{i=1}^n{x}_i^2$	$n\le 30$	${\left[-100, 100\right]}^n$
Rastrigin	$f_2\left(x\right)={\displaystyle \sum }_{i=1}^n[x_i^2-\cos \left(2\pi x_i\right)+10]$	$n\le 30$	${\left[-5, 5\right]}^n$
Rosenbrock	$f_3\left(x\right)={\displaystyle \sum }_{i=1}^n[100{(x_{i+1}-x_i^2)}^2+{\left(x_i-1\right)}^2]$	$n\le 10$	${\left[-5, 5\right]}^n$
Schaffer	$f_4\left(x\right)=\frac{{\text{sin}}^2\sqrt{x_1^2+x_2^2}-0.5}{{\left[1.0+0.001\left(x_1^2+x_2^2\right)\right]}^2}+\frac{1}{2}$	$n\le 2$	${\left[-100, 100\right]}^n$

Table 2.

Details of four benchmark functions.

In order to test the performance of the above PSO algorithms, four benchmark optimization functions are used for verification whose details are shown in Table 2. The parameters in PSO algorithms are set as follows: the number of particle swarm population is set as $M=50$, the maximum velocity is set as $V_{\text{max}}=3.0$, the inertia weight coefficient is chose according to Eq. (4). The acceleration constants are set as $c_1=1.2$and $c_2=0.8$in SPSO and DEPSO algorithms, whereas $c_1=2.0$and $c_2=2.0$in standard PSO algorithm, respectively. The objective functions are set as the same as the test functions and the iteration stops upon the fitness value reaches the maximum iteration number t_max=1000 or the objective function value is less than 10^-30. The test results of the four functions with dimension n=2 are shown in Figure 3. Table 3 lists the average objective function values and the corresponding mean squared errors for these four benchmark functions with 1000 independent runs. As shown, both the SPSO and DEPSO algorithms achieve better performance than the standard PSO algorithm in terms of computational accuracy.

Figure 3.

The images of (a) Sphere function, (b) Rastrigin function, (c) Rosenbrock function and (d) Schaffer function with dimensionn=2.

Function	Dimension	PSO	SPSO	DEPSO
Sphere	10	5.94×10^-28±3.46×10^-27	8.33×10^-31±1.67×10^-31	8.25×10^-31±1.91×10^-31
	20	2.82×10^-13±7.76×10^-13	4.80×10^-17±3.63×10^-16	2.98×10^-25±1.54×10^-24
	30	5.27×10^-8±8.6910^-8	1.23×10^-6±3.97×10^-6	7.30×10^-15±2.34×10^-14
Rastrigin	10	4.82×10^-28±1.82×10^-27	8.35×10^-31±1.62×10^-31	6.87×10^-31±2.16×10^-31
	20	2.68×10^-13±5.52×10^-13	7.22×10^-17±9.23×10^-16	2.03×10^-25±1.46×10^-24
	30	5.39×10^-8±7.56×10^-8	1.48×10^-6±4.34×10^-6	7.08×10^-15±6.66×10^-14
Rosenbrock	3	3.36×10^-8±8.39×10^-8	3.53×10^-17±1.06×10^-16	2.95×10^-16±9.21×10^-16
	5	5.62×10^-2±6.02×10^-2	4.12×10^-2±1.99×10^-2	1.37×10^-2±2.85×10^-2
	10	3.18×10^0±1.56×10^0	4.62×10^0±3.45×10^0	5.78×10^0±2.61×10^0
Schaffer	2	5.11×10^-2±5.61×10^-2	9.87×10^-4±2.92×10^-3	4.66×10^-4±2.60×10^-3

Table 3.

The retrieval results of three test functions with 1000 independent runs.

4. Inverse geometry design of two-dimensional radiative enclosures

4.1. Description of the inverse geometry design problem

Considering a radiative equilibrium problem in two-dimensional irregular enclosures filled with participating media whose schematic diagram is shown in Figure 4. The curve EF represents the design surface and the design purpose is to produce a uniform distribution of radiative heat flux on the design surface. The bottom surface AD is the heating surface which can considered as the radiative heat source and its temperature is fixed as T_S. The two side surfaces AB and CD are cold with temperatures of 0K.

Figure 4.

Physical model of inverse geometry design.

In order to optimize the geometry shape of the design surface to meet the specified requirement, the objective function (being equal to the fitness function in PSO algorithms) is defined as the square residuals between the estimated and average dimensionless radiative heat flux values which can be expressed as

$F_{\mathrm{o}\mathrm{b}\mathrm{j}}={\displaystyle \sum _{i=1}^N{\left(Q_{{\mathrm{w}}_i}-Q_{\mathrm{a}\mathrm{v}}\right)}^2}$E9

where Nis the number of the computational node on the design surface, and Q_wi and Q_av are the dimensionless radiative heat flux of the ith node and average value on the design surface, respectively. The iteration stop criterion of PSO is defined as

$F_{\mathrm{o}\mathrm{b}\mathrm{j}}<\epsilon$E10

where εis a small positive value. The smaller the value of εis, the better the homogenization degree will be.

To evaluate the optimization results, the relative error is defined as

${\epsilon }_{\mathrm{r}\mathrm{e}\mathrm{l}}=\left|\frac{Q_{{\mathrm{w}}_i}-Q_{\mathrm{a}\mathrm{v}}}{Q_{\mathrm{a}\mathrm{v}}}\right|\times 100\%$E11

4.2. Akima cubic interpolation

The design surface is discretized into a series of control points, and Akima cubic interpolation is used to approximate the geometry shape of the surface in optimization process. Akima interpolation is formulated by a cubic polynomial between two control points. First, the prerequisite is introduced as [33]

$\left\{\begin{array}{l}Y\left(x_k\right)=y_k\\ Y\left(x_{k+1}\right)=y_{k+1}\\ Y^{\text{'}}\left(x_k\right)=g_k\\ Y^{\text{'}}\left(x_{k+1}\right)=g_{k+1}\end{array}\right.$E12

where g_kis the slope of the curve at the position x_k, which can be defined as [33]

$g_k=\left\{\begin{array}{l}\frac{l_{k-1}+l_k}{2}{l}_{k+1}=l_k\mathrm{and}{l}_{k-1}=l_{k-2}\\ \frac{\left|l_{k+1}-l_k\right|l_{k-1}+\left|l_{k-1}-l_{k-2}\right|l_k}{\left|l_{k+1}-l_k\right|+\left|l_{k-1}-l_{k-2}\right|}\mathit{else}\end{array}\right.$E13

and the function l_kcan be expressed as [33]

$l_k=\frac{y_{k+1}-y_k}{x_{k+1}-x_k}$E14

At the endpoints, the value of the function l_kcan be defined as [33]

$\left\{\begin{array}{l}{l}_0=2l_1-l_2,l_{-1}=2l_0-l_0\mathit{at}\mathit{the}\mathit{left}\mathit{endpoint}\\ l_{N_d}=2l_{N_d-1}-l_{N_d-2},l_{N_d+1}=2l_{N_d}-l_{N_d-1}\mathit{at}\mathit{the}\mathit{right}\mathit{endpoint}\end{array}\right.$E15

If the above equations are satisfied, then the cubic polynomial in the subinterval [$x_k$, $x_{k+1}$] can be determined as [33]

$Y\left(x\right)=C_1+C_2\left(x-x_k\right)+C_3{\left(x-x_k\right)}^2+C_4{\left(x-x_k\right)}^3$E16

where C₁, C₂, C₃, and C₄ are polynomial coefficients, which can be calculated as [33]

$\left\{\begin{array}{l}{C}_1=y_k\\ C_2=g_k\\ C_3=\frac{3\left(u_k-2g_k-{g_k}_{+1}\right)}{x_{k+1}-x_k}\\ C_4=\frac{g_k+2g_{k+1}-2u_k}{{\left(x_{k+1}-x_k\right)}^2}\end{array}\right.$E17

(x, y)	(x1, y1)	(x2, y2)	(x3, y3)	(x4, y4)	(x5, y5)	(x6, y6)
Case 1	(0.0, 1.0)	(0.2, 1.2)	(0.4, 1.5)	(0.6, 1.5)	(0.8, 1.2)	(1.0, 1.0)
Case 2	(0.0, 1.0)	(0.2, 1.5)	(0.4, 1.3)	(0.6, 1.3)	(0.8, 1.5)	(1.0, 1.0)
Case 3	(0.0, 1.0)	(0.28, 1.2)	(0.41, 1.4)	(0.67, 1.1)	(0.88, 1.3)	(1.0, 1.0)

Table 4.

Coordinates of control points of three test cases.

Figure 5.

Curve fitting results by means of Akima cubic interpolation.

In order to test the performance of the Akima cubic interpolation, three interpolation cases are applied, in which six specified points are used as control points whose coordinates are shown in Table 4. The curves in Figure 5 indicate the Akima cubic interpolation can successfully applied for geometry shape fitting.

4.3. Discrete ordinate method with a body-fitted coordinate

The boundary shape changes with the optimization of the design surface and the computational domain must be re-meshed, which greatly increase the solving difficulty of the inverse geometry design problems. In addition, the forward radiative heat transfer problem cannot be precisely solved by the normal numerical method. For the purpose of fitting the irregular boundary shape, the DOM with a body-fitted coordinate system is adopted to solve the RTE. For the participating media, the forward can be written as [34]

$\frac{\partial I\left(s,\widehat{\mathbf{s}}\right)}{\partial s}=-{\beta }_{\mathrm{e}}I\left(s,\widehat{\mathbf{s}}\right)+{\kappa }_{\mathrm{a}}{I}_{\mathrm{b}}\left(s\right)+\frac{{\kappa }_{\mathrm{s}}}{4\pi }{\displaystyle {\int }_{4\pi }I\left({\widehat{\mathbf{s}}}_i\right)}\Phi \left({\widehat{\mathbf{s}}}_i,\widehat{\mathbf{s}}\right)d{\Omega }_i$E18

which is an integro-differential type, where Iis the function of position sand direction ŝ. β_e, κ_a, and κ_s are the extinction, absorption, and scattering coefficients of media, respectively. The Φ(ŝ_i, ŝ) is the scattering phase function between incoming direction ŝ_iand scattering direction ŝ, which can be defined as $\Phi =1.0+acos\left({\stackrel{}{\mathrm{ŝ}}}_i\cdot \stackrel{}{\mathrm{ŝ}}\right)$. The coefficient “a” are different for different scattering characteristics of media and values are a=0, a=1 and a=-1 for the isotropic scattering, the forward scattering and the backward scattering of media. The forward radiative heat transfer problem in the irregular enclosures is solved by using the computationally feasible DOM with a body-fitted coordinate system in this research. The 2D RTE discretized by the DOM can be expressed as [34, 35]

${\alpha }^m\frac{\partial I^m}{\partial x}+{\beta }^m\frac{\partial I^m}{\partial y}=-{\beta }_{\mathrm{e}}{I}^m+{\kappa }_{\mathrm{a}}{I}_{\mathrm{b}}\left(s\right)+\frac{{\kappa }_{\mathrm{s}}}{4\pi }\left({\displaystyle \sum _{m\text{'}=1}^{N\Omega }{w}^{m\text{'}}{I}_{\mathrm{w}}^{m\text{'}}{\Phi }^{m\text{'},m}}\right)$E19

where α^mand β^mare the direction cosines of the discrete direction m. srefers to the spatial position, wis the quadrature weight. The radiative boundary condition can be directly imposed as follows [34]

$I_{\mathrm{w}}^m={\epsilon }_{\mathrm{w}}{I}_{\mathrm{b},\mathrm{w}}+\frac{1-{\epsilon }_{\mathrm{w}}}{\pi }{\displaystyle \sum _{{\mathbf{n}}_w\cdot {\mathbf{s}}^{m\text{'}}<0}{w}^{m\text{'}}}{I}_{\mathrm{w}}^{m\text{'}}\left|{\mathbf{n}}_{\mathrm{w}}\cdot {\mathbf{s}}^{m\text{'}}\right|,{\mathbf{n}}_{\mathrm{w}}\cdot {\mathbf{s}}^{m\text{'}}>0$E20

where ε_w is the emissivity of boundaries, I_b,w is radiative intensity of blackbody boundaries, n_w represents the unit normal vector on the boundary, and sdenotes the direction of radiative transfer.

The Jacobian matrix is used for coordinate transformation to solve the radiative heat transfer in irregular enclosures

$J=\frac{\partial \left(x,y\right)}{\partial \left(\xi ,\eta \right)}$E21

The spatially discretized RTE with a body-fitted coordinate system can be expressed as

$\begin{array}{l}{\left[I^{m}J\left({\alpha }^m\frac{\partial \xi }{\partial x}+{\beta }^m\frac{\partial \xi }{\partial y}\right)\right]}_e\Delta \eta -{\left[I^{m}J\left({\alpha }^m\frac{\partial \xi }{\partial x}+{\beta }^m\frac{\partial \xi }{\partial y}\right)\right]}_w\Delta \eta +{\left[I^{m}J\left({\alpha }^m\frac{\partial \eta }{\partial x}+{\beta }^m\frac{\partial \eta }{\partial y}\right)\right]}_n\Delta \xi \\ -{\left[I^{m}J\left({\alpha }^m\frac{\partial \eta }{\partial x}+{\beta }^m\frac{\partial \eta }{\partial y}\right)\right]}_s\Delta \xi =J{\left[-{\beta }_{\mathrm{e}}{I}^m+{\kappa }_{\mathrm{a}}{I}_{\mathrm{b}}+\frac{{\kappa }_{\mathrm{s}}}{4\pi }\left({\displaystyle \sum _{m\text{'}=1}^{N\Omega }{w}^{m\text{'}}{I}_w^{m\text{'}}{\Phi }^{m\text{'},m}}\right)\right]}_P\mathit{\Delta \xi \Delta \eta }\end{array}$E22

where Prepresents the central node of the control volume. The subscripts e, w, n, and srepresent the eastern, western, northern and southern boundaries around P, respectively.

The step scheme is applied to solve the above equations, and Eq. (22) can be expressed as

$a_P^m{I}_P^m=a_E^m{I}_E^m+a_W^m{I}_W^m+a_N^m{I}_N^m+a_S^m{I}_S^m+b_P^m$E23

where the subscripts E, W, N, and Srepresent the central node of eastern, western, northern, and southern control volumes around control volume P, and

$\begin{array}{l}{a}_E^m=max\left\{-{\left[J\left({\alpha }^m\frac{\partial \xi }{\partial x}+{\beta }^m\frac{\partial \xi }{\partial y}\right)\right]}_e\Delta \eta ,0\right\}{a}_W^m=max\left\{{\left[J\left({\alpha }^m\frac{\partial \xi }{\partial x}+{\beta }^m\frac{\partial \xi }{\partial y}\right)\right]}_w\Delta \eta ,0\right\}\\ a_N^m=max\left\{-{\left[J\left({\alpha }^m\frac{\partial \eta }{\partial x}+{\beta }^m\frac{\partial \eta }{\partial y}\right)\right]}_n\Delta \xi ,0\right\}{a}_S^m=max\left\{{\left[J\left({\alpha }^m\frac{\partial \eta }{\partial x}+{\beta }^m\frac{\partial \eta }{\partial y}\right)\right]}_s\Delta \xi ,0\right\}\\ a_P^m=max\left\{{\left[J\left({\alpha }^m\frac{\partial \xi }{\partial x}+{\beta }^m\frac{\partial \xi }{\partial y}\right)\right]}_e\Delta \eta ,0\right\}+max\left\{-{\left[J\left({\alpha }^m\frac{\partial \xi }{\partial x}+{\beta }^m\frac{\partial \xi }{\partial y}\right)\right]}_w\Delta \eta ,0\right\}\\ +max\left\{{\left[J\left({\alpha }^m\frac{\partial \eta }{\partial x}+{\beta }^m\frac{\partial \eta }{\partial y}\right)\right]}_n\Delta \xi ,0\right\}+max\left\{-{\left[J\left({\alpha }^m\frac{\partial \eta }{\partial x}+{\beta }^m\frac{\partial \eta }{\partial y}\right)\right]}_s\Delta \xi ,0\right\}\\ +J_P{\beta }_{\mathrm{e},P}\mathit{\Delta \xi \Delta \eta }\\ b_P^m=J_P{\left[{\kappa }_{\mathrm{a}}\frac{\sigma T_P^4}{\pi }+\frac{{\kappa }_{\mathrm{s}}}{4\pi }\left({\displaystyle \sum _{m\text{'}=1}^{N\Omega }{w}^{m\text{'}}{I}_P^{m\text{'}}{\Phi }^{m\text{'},m}}\right)\right]}_P\mathit{\Delta \xi \Delta \eta }\end{array}$E24

Eq. (24) can be expressed in matrix form

$\mathbf{A}\mathbf{\Psi }=\mathbf{B}$E25

where Arepresents the five-diagonal non-symmetric coefficient matrix, Ψrepresents the vector that consists of the variables I^mat grid nodes, and Brepresents the vector that consists of the variables b^mon the right side of Eq. (23). The conjugate gradients stabilized (CGSTAB) method is adopted to solve the final discretized RTE because of its stability and fast convergence rate [35].

Figure 6.

The schematic diagram of grids in computational domains.

Figure 7.

The dimensionless radiative heat flux distribution on the bottom boundary for different absorption coefficients.

Consider a non-radiative equilibrium problem in the two-dimensional irregular enclosure filled with participating media. The four boundaries are cold surfaces whose temperature is 0K and all the boundaries are assumed as blockbody. The temperature of media is set as T_g, which can be considered as the heat source. The computational grids are meshed as 10 × 10 and the schematic in the computational domain is shown in Figure 6. The radiative heat transfer problem for different absorption coefficients in the radiative enclosure is solved by the body-fitted DOM and the retrieval results are compared with the results obtained by FVM in Ref. [36] which is shown in Figure 7. The curves show that the retrieval results achieve good agreements with FVM, which verified the accuracy and reliability of the computational program for solving radiative problems in irregular enclosures.

4.4. Inverse design results and discussions

The inverse geometry design model in Section 4.1 is considered and the standard PSO algorithm is abbreviated as PSO if there is no special instruction. The initial shape of the enclosure is rectangular, whose size is set as $L_x\times L_y=1.0\times 1.0\mathrm{m}$. The absorption and scattering coefficients of media are set as ${\kappa }_{\mathrm{a}}=2.0$m^-1 and ${\kappa }_{\mathrm{s}}=0.5$m^-1, respectively. The parameters in the PSO algorithms are set as the same as the ones in test cases of Section 3. The stopping criteria are set as follows: (1) the objective function value is less than 10^-7, (2) the iteration times reach a maximum of 1000 and (3) the fitness value remains unchanged in consequent 100 iterations. The initial and final geometry shape of the design surfaces and their dimensionless radiative heat flux distributions are shown in Figure 8(a) and 8(b), respectively. The curves demonstrate that a uniform distribution of radiative heat flux on the design surface is obtained by means of PSO algorithms. Figure 8(c) shows the relative error distributions of the dimensionless radiative heat flux on the design surface. The average relative errors are 0.0310%, 0.0234% and 0.0279%, and the maximum relative errors are 0.0728%, 0.0540% and 0.0807% for PSO, SPSO and DEPSO algorithms, respectively. Overall, the retrieved results clearly reveal that the specified requirement of producing a uniform radiative heat flux distribution on the design surface can be obtained by PSO algorithm.

Figure 8.

(a) Initial and final geometry shape of the design surface, (b) dimensionless radiative heat flux distribution on the design surface, and (c) relative error distributions of dimensionless radiative heat flux on the design surface, of a two-dimensional radiative enclosure.

Algorithms	CPU; time (s)	Fitness values	Average relative; error (%)	Maximum relative; error (%)
PSO	40037	7.28×10^-8	0.0313	0.0737
SPSO	31572	4.05×10^-8	0.0244	0.0518
DEPSO	36944	5.62×10^-8	0.0275	0.0669

Table 5.

Comparison of inverse design results by PSO, SPSO and DEPSO algorithms.

In view of the random characteristic of intelligent algorithms, all the tests are repeated 50 trials to compare the performance of PSO algorithms. Table 5 shows the comparison of the results obtained by PSO, SPSO and DEPSO algorithms. It can be found that all the PSO algorithms have reached the special design requirement and both SPSO and DEPSO achieve better performance than the initial PSO in terms of computational accuracy and efficiency.

Figure 9.

(a) Geometry shape of the design surface and (b) dimensionless radiative heat flux distribution on the design surface by means of SPSO algorithm for different numbers of control points.

In order to enhance the computational efficiency of inverse geometry design problems, the effects of corresponding parameters are investigated in this study. For the fact that the design surface is discretized into a series of control points in the inverse design process, the number of control point has a direct impact on the inverse geometry design results. The radiative physical parameters of media are kept as the same as the above typical example and the numbers of control point are set as N_d=1, 3, 5, and 7, respectively. The SPSO algorithm is adopted as the inverse design method and the initial and optimized geometry shape and dimensionless radiative heat flux on the design surface are shown in Figure 9, respectively. The curves show that the special design requirement is satisfied under the conditions that N_d=3, 5 and 7, whereas the homogenization degree is relatively poor in the case that N_d=1. Table 6 lists the iteration numbers, fitness values and relatives errors for different numbers of control point. It can be found that both the iteration number and the relative error are smallest when N_d=3, which is because too many control points will decrease the sensitivity of radiative heat flux to the shape changing of design surface and one control point cannot provide enough necessary information for the geometry shape of the boundary [2, 9].

Control point; numbers	Iteration; numbers	Fitness values	Average relative; error (%)	Maximum relative; error (%)
1	164	3.89×10^-8	0.2711	0.5597
3	19	5.36×10^-8	0.0280	0.0796
5	26	6.72×10^-8	0.0309	0.0875
7	37	7.63×10^-8	0.0321	0.1024

Table 6.

Comparison of inverse design results for different number of control points.

Figure 10.

Geometry shape of the design surface for different extinction coefficients of media.

Figure 11.

Geometry shape of the design surface for different scattering albedo of media.

The physical properties of the media have an important influence on the energy transfer and then affect the optimization results of radiative enclosures. The effects of the extinction coefficient and the scattering albedo on the inverse design results are studied here. The scattering albedo of media is fixed as $\omega =0.5$, the extinction coefficients are set as $\beta =1.0$, 3.0 and 5.0, respectively. The initial guessed and final optimized geometry shapes of the design surface are shown in Figure 10. It can be seen that three PSO algorithms achieve similar boundaries and the geometry shapes are significantly different for different extinction coefficients. The extinction coefficient of media is fixed as $\beta =3.0$, the scattering albedos are set as $\omega =0.1$, 0.5 and 0.9, respectively. The optimized geometry shapes of the design surface by means of SPSO are shown in Figure 11. It can be seen that the designed boundaries are very close to each other for different scattering albedos. Tables 7 and 8 list the detailed inverse design results for different extinction coefficients and scattering albedos, respectively. As shown, both the optimized geometry the retrieval results for different extinction coefficients are significantly different, and both the height of the final designed surface and the dimensionless radiative heat flux decrease with the extinction coefficient increases. However, the retrieval results for different scattering albedo are close. Therefore, the scattering albedo of the media has little influence on the inverse geometry design results when the extinction coefficient is determined.

Algorithm	Extinction; coefficient	Fitness; values	Dimensionless; radiative heat flux	Average relative; error (%)	Maximum; relative error (%)
PSO	${\beta }_{\mathrm{e}}=1.0$	8.22×10^-8	0.2333	0.0329	0.0895
	${\beta }_{\mathrm{e}}=3.0$	7.19×10^-8	0.1137	0.0496	0.0736
	${\beta }_{\mathrm{e}}=5.0$	3.94×10^-8	0.0735	0.0460	0.1297
SPSO	${\beta }_{\mathrm{e}}=1.0$	7.61×10^-8	0.2334	0.0295	0.0816
	${\beta }_{\mathrm{e}}=3.0$	5.63×10^-8	0.1136	0.0398	0.0723
	${\beta }_{\mathrm{e}}=5.0$	2.42×10^-8	0.0734	0.0301	0.0909
DEPSO	${\beta }_{\mathrm{e}}=1.0$	5.07×10^-8	0.2334	0.0230	0.0849
	${\beta }_{\mathrm{e}}=3.0$	4.16×10^-8	0.1136	0.0351	0.0706
	${\beta }_{\mathrm{e}}=5.0$	2.85×10^-8	0.0734	0.0307	0.0889

Table 7.

Inverse geometry design results for different extinction coefficients.

Algorithm	Scattering; albedo	Fitness; values	Dimensionless; radiative heat flux	Average relative; error (%)	Maximum; relative error (%)
PSO	$\omega =0.1$	7.93×10^-8	0.1137	0.0459	0.0867
	$\omega =0.5$	7.19×10^-8	0.1137	0.0496	0.0736
	$\omega =0.9$	7.23×10^-8	0.1136	0.0355	0.636
SPSO	$\omega =0.1$	3.02×10^-8	0.1136	0.0297	0.0804
	$\omega =0.5$	5.63×10^-8	0.1136	0.0398	0.0723
	$\omega =0.9$	4.49×10^-8	0.1136	0.0322	0.0665
DEPSO	$\omega =0.1$	5.07×10^-8	0.1136	0.0352	0.0349
	$\omega =0.5$	4.16×10^-8	0.1136	0.0351	0.0706
	$\omega =0.9$	3.49×10^-8	0.1136	0.0358	0.0710

Table 8.

Inverse geometry design results for different scattering albedo of media.

The scattering of media will affect the original transfer direction of radiative heat or energy, so the scattering characteristic of media also has influence on the inverse geometry design of radiative enclosures. The extinction coefficient and scattering albedo of the media are set as $\beta =3.0$and $\omega =0.5$, respectively. The initial and final geometry shapes of the design surface by using SPSO algorithm for three kinds of scattering characteristics of media are shown in Figure 12. It can be found that the optimal geometry shapes of the design surface are close. Table 9 plots the inverse geometry design results for different scattering characteristics. As shown, the values of dimensionless radiative heat flux on the design surface and relative errors are significantly different. Therefore, the scattering characteristics of media mainly affect the radiative heat flux on the boundaries and have no obvious effect on the final geometry shape of inverse radiative design problems.

Figure 12.

Geometry shape of the design surface scattering characteristics of media.

Algorithm	Scattering; characteristic	Fitness; values	Dimensionless; radiative heat flux	Average relative; error (%)	Maximum; relative error (%)
PSO	a = 1	8.66×10^-8	0.1297	0.0354	0.0828
	a = 0	9.04×10^-8	0.1136	0.0403	0.0776
	a = − 1	7.21×10^-8	0.1006	0.0385	0.0941
SPSO	a = 1	7.57×10^-8	0.1297	0.0328	0.0982
	a = 0	8.16×10^-8	0.1136	0.0463	0.0825
	a = − 1	1.60×10^-8	0.1007	0.0195	0.0629
DEPSO	a = 1	5.53×10^-8	0.1297	0.0378	0.0797
	a = 0	6.94×10^-8	0.1136	0.0426	0.0735
	a = − 1	5.71×10^-8	0.1007	0.0413	0.0882

Table 9.

Inverse geometry design results for different scattering characteristics of media.

5. Conclusions

In this chapter, the basic theoretical principles of PSO algorithm is introduced in detail and three kinds of PSO algorithms—standard PSO, SPSO and DEPSO—are applied for solving the inverse geometry design problem of a two-dimensional radiative enclosure filled with participating media. The design purpose is to produce a uniform distribution of radiative heat flux on the designed surface. The design surface is discretized into a series of control points and the Akima cubic interpolation is used to approximate the geometry shape of the boundary. The radiative heat transfer problem in the irregular enclosures is solved by the DOM with a body-fitted coordinate system. The pre-required radiative heat flux distribution is satisfied by optimizing the positions of control points based on the PSO algorithms. The retrieval results show that PSO algorithms can be successfully applied to solve inverse geometry design problems and SPSO achieves the best performance on computational time. Meanwhile, the scattering albedo and scattering characteristics of media have little effect on the geometry shape of the design surface. To improve the computational efficiency, the number of control points is recommended as N_d=3.

Acknowledgments

The supports of this work by the National Natural Science Foundation of China (Nos. 51476043 and 51576053), the Major National Scientific Instruments and Equipment Development Special Foundation of China (No. 51327803), and the Foundation for Innovative Research Groups of the National Natural Science Foundation of China (No. 51421063) are gratefully acknowledged. A very special acknowledgement is also made to the editors and referees who make important comments to improve this chapter.