**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**.

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

**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**

**Figure 2.** The flowchart of the basic PSO algorithm.

48 Optimization Algorithms- Methods and Applications

according to the following steps:

determine the local best location **P***i*

The main procedure of PSO algorithm for solving optimization problems can be carried out

**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

**Step 2**: Fitness evaluation: Evaluate each particle's fitness according to its location and

(*t*) and global best location **P**g(*t*).

numbers *t*max and the stop criterion *ε*, etc. Set the current iteration as *t*=1.

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 *w* is introduced to control the impact of the current velocity on the next velocity. The velocity formula of *i*th particle can be expressed as [28]

$$\mathbf{V}\_{i}(t+1) = \mathbf{w} \cdot \mathbf{V}\_{i}(t) + \mathbf{c}\_{1} \cdot r\_{1} \cdot \left[\mathbf{P}\_{i}(t) - \mathbf{X}\_{i}(t)\right] + \mathbf{c}\_{2} \cdot r\_{2} \cdot \left[\mathbf{P}\_{g}(t) - \mathbf{X}\_{i}(t)\right] \tag{3}$$

where *w* is 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

$$
\delta \mathbf{w} = \mathbf{w}\_{\text{max}} - \frac{t}{t\_{\text{max}}} \cdot \left(\mathbf{w}\_{\text{max}} - \mathbf{w}\_{\text{min}}\right) \tag{4}
$$

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 *i*th particle at *t*+1 iteration is determined by three parameters of *t* iteration, namely **X***i* (*t*), **P***i* (*t*) and **P**g(*t*). The new velocity of *i*th particle can be expressed as [30]

$$\mathbf{V}\_{i}(t+1) = c\_{1} \cdot r\_{1} \cdot \left[\mathbf{P}\_{i}(t) - \mathbf{X}\_{i}(t)\right] + c\_{2} \cdot r\_{2} \cdot \left[\mathbf{P}\_{g}(t) - \mathbf{X}\_{i}(t)\right] \tag{5}$$

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{aligned} \mathbf{P}\_{i} &= \mathbf{X}\_{j}(t+1) \\ \mathbf{P}\_{i} &= \begin{bmatrix} \mathbf{P}\_{i} & F\left(\mathbf{P}\_{i}\right) < F\left[\mathbf{X}\_{i}(t+1)\right] \\ \mathbf{X}\_{i}(t+1), & F\left(\mathbf{P}\_{i}\right) \ge F\left[\mathbf{X}\_{i}(t+1)\right] \end{bmatrix} \\ \mathbf{P}\_{\mathcal{S}}^{\top} &= \arg\min\left\{ F\left(\mathbf{P}\_{i}\right) \middle| i=1,\cdots,\cdots M \right\} \\ \mathbf{P}\_{\mathcal{S}} &= \arg\min\left\{ F\left(\mathbf{P}\_{\mathcal{S}}^{\top}\right), F\left(\mathbf{P}\_{\mathcal{S}}\right) \right\} \end{aligned} \tag{6}$$

After above updates, the following criterions are executed:

(1) If **P**g = **P***j* , which demonstrates that the random location **X***j* is the global best location. In this situation the *j*th particle will not be updated on the basis of Eq. (5) and the algorithm will randomly generate a location **X***j* in 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***g* and **P***j* are updated, respectively.

(2) If **P**g ≠ **P***j* and **P***g* has 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***j* and **P***g* has been updated, which demonstrates that there is a location of *k*th particle (*k* ≠ *j*) meets the requirement **X***k* (*t* + 1)=**P***k* =**P***g*, namely the local best location of *k*th particle is the global best location and it is better than the global best location at the last iteration. In this situation, the *k*th 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***g* and **P***j* are 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**

max ( max min ) max

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

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

In SPSO algorithm, the inertia weight coefficient is set as *w*=0. Hence, the velocity of *i*th particle

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

updated based on the Eq. (5). The whole amending process can be expressed by the following

( ) ( ) ( ) { ( ) }

situation the *j*th particle will not be updated on the basis of Eq. (5) and the algorithm will randomly generate a location **X***j* in 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***g*

1 1 1

*t F Ft Fi M*

ì <é +ù ï ë û = í + ³é +ù ïî ë û

{ ( ) ( )}

'

*F F*

*i ii*

**P PX**

**X PX**

*i ii*

, arg min ,

arg min ,

**P PP**

= =

*g g g*

( ) ( )

1

*F Ft*

LL

, which demonstrates that the random location **X***j* is the global best location. In this

*i* (*t cr t t cr t t* + = × ×é - ù+ × × - ) *i i* ( ) ( ) é ù g ( ) *i* ( ) ë û ë û **V PX P X** 1 1 2 2 1 (5)

at *t*+1 iteration is determined by three parameters of *t* iteration, namely **X***i*

( )

1

*t*

,

= +

*g i*

**P P**

*j j*

**P X**

The new velocity of *i*th particle can be expressed as [30]

a particle in the searching space whose location is **X***j*

'

are updated, respectively.

=

After above updates, the following criterions are executed:

*i*

**P**

applied in engineering fields.

50 Optimization Algorithms- Methods and Applications

**2.3. Stochastic PSO algorithm**

particle swarms.

equations

(1) If **P**g = **P***j*

and **P***j*

*t ww w w t* = -× - (4)

(*t*), **P***i*

(*t* + 1), and other particles' locations are

(*t*) and **P**g(*t*).

(6)

Differential Evolution (DE) algorithm adopts simple differential operation among potential solutions to produce new candidate solution, which is a parallel, direct and stochastic search‐ ing technique. It was first proposed for solving Chebyshev polynomials and global optimiza‐ tion 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

$$\boldsymbol{\beta} = \boldsymbol{\chi} \left[ \mathbf{X}\_{\boldsymbol{\epsilon}\_1}(t) - \mathbf{X}\_{\boldsymbol{\epsilon}\_1}(t) \right] \tag{7}$$

where *χ* is the differential mutation operator which controls the magnification of differential variation **X***r*1 (*t*)−**X***r*2 (*t*) 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 *i*th particle can be expressed as

$$\mathbf{X}\_{i}\left(t+1\right) = \mathbf{X}\_{i}\left(t\right) + \mathbf{V}\_{i}\left(t\right) + \mathcal{J}\left(C - F\_{\text{min}}\right) \tag{8}$$

where *C* represents 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 *i*th 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 *i*th particle is updated according to Eq. (2) and the mutation operation of *i*th 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**

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*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.

**Table 2.** Details of four benchmark functions.

Inverse Geometry Design of Radiative Enclosures Using Particle Swarm Optimization Algorithms http://dx.doi.org/10.5772/62351 53

of DEPSO algorithm, the following judgment should be executed before the location is

(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 *i*th particle is updated

(2) If the new fitness value is worse than before, which indicates this mutation is failed. Then the location of *i*th particle is updated according to Eq. (2) and the mutation operation of *i*th

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.

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*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

**Function Expression Dimension Search Space**

*n* £ 30 [ 100,100]

*n* £ 30 [ 5,5]

*n* £10 [ 5,5]

*n* £ 2 [ 100,100]

*n* -

*n* -

*n* -

*n* -

will continue at the next iteration until the mutation is successful.

**3. Simulation test of PSO algorithms**

updated:

accuracy.

Rosenbrock

Schaffer

Sphere ( ) 2

1

Rastrigin ( ) ( ) 2 2

1 *n i i fx x* = =å

*n*

1

3 1 1

*x x f x*

=

*n*

*i*

( )

**Table 2.** Details of four benchmark functions.

*i fx x x*

[ cos 2 10]

=- + å

p

[100( ) 1 ]

*ii i*

= - +- å

( )

*x x* + - = + é ù + + ë û

1 2 sin 0.5 1 2 1.0 0.001

2 2 2

*i i*

( ) ( )

22 2 1 2 4 2 2 2

*fx x x x* + =

according to Eq. (8).

52 Optimization Algorithms- Methods and Applications

**Figure 3.** The images of (a) Sphere function, (b) Rastrigin function, (c) Rosenbrock function and (d) Schaffer function with dimension *n*=2.

**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 *TS*. 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 require‐ ment, 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\_{\rm obj} = \sum\_{i=1}^{N} \left(\underline{Q}\_{\rm w\_i} - \underline{Q}\_{\rm w}\right)^2 \tag{9}$$

where *N* is the number of the computational node on the design surface, and *Q*wi and *Q*av are the dimensionless radiative heat flux of the *i*th node and average value on the design surface, respectively. The iteration stop criterion of PSO is defined as

$$F\_{\text{obj}} < \mathfrak{c} \tag{10}$$

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

$$
\varepsilon\_{\rm rel} = \left| \frac{\mathcal{Q}\_{\rm w} - \mathcal{Q}\_{\rm w}}{\mathcal{Q}\_{\rm w}} \right| \times 100\% \tag{11}
$$

#### **4.2. Akima cubic interpolation**

**4. Inverse geometry design of two-dimensional radiative enclosures**

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 *TS*. The two side

In order to optimize the geometry shape of the design surface to meet the specified require‐ ment, 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

> obj ( w av ) *i N i F QQ* =

1

*F*obj < e

respectively. The iteration stop criterion of PSO is defined as

= - å 2

where *N* is the number of the computational node on the design surface, and *Q*wi and *Q*av are the dimensionless radiative heat flux of the *i*th node and average value on the design surface,

where *ε* is a small positive value. The smaller the value of *ε* is, the better the homogenization

(10)

(9)

**4.1. Description of the inverse geometry design problem**

54 Optimization Algorithms- Methods and Applications

surfaces AB and CD are cold with temperatures of 0K.

**Figure 4.** Physical model of inverse geometry design.

values which can be expressed as

degree will be.

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]

$$\begin{cases} Y(\mathbf{x}\_k) = \mathbf{y}\_k \\ Y(\mathbf{x}\_{k+1}) = \mathbf{y}\_{k+1} \\ Y'(\mathbf{x}\_k) = \mathbf{g}\_k \\ Y'(\mathbf{x}\_{k+1}) = \mathbf{g}\_{k+1} \end{cases} \tag{12}$$

where *gk* is the slope of the curve at the position *xk*, which can be defined as [33]

$$\mathbf{g}\_{k} = \begin{cases} \frac{I\_{k-1} + I\_k}{2} & I\_{k+1} = I\_k \text{ and } I\_{k-1} = I\_{k-2} \\\\ \frac{\left| I\_{k+1} - I\_k \right| I\_{k-1} + \left| I\_{k-1} - I\_{k-2} \right| I\_k}{\left| I\_{k+1} - I\_k \right| + \left| I\_{k-1} - I\_{k-2} \right|} & \text{else} \end{cases} \tag{13}$$

and the function *lk* can be expressed as [33]

$$I\_k = \frac{y\_{k+1} - y\_k}{x\_{k+1} - x\_k} \tag{14}$$

At the endpoints, the value of the function *lk* can be defined as [33]

$$\begin{cases} l\_0 = 2l\_1 - l\_{2'} & l\_{-1} = 2l\_0 - l\_0 \\ l\_{N\_d} = 2l\_{N\_d - 1} - l\_{N\_d - 2'} & l\_{N\_d + 1} = 2l\_{N\_d} - l\_{N\_d - 1} \end{cases} \text{ at the left endpoint} \tag{15}$$

If the above equations are satisfied, then the cubic polynomial in the subinterval [*xk* , *xk* +1] can be determined as [33]

$$\mathbf{Y}(\mathbf{x}) = \mathbf{C}\_1 + \mathbf{C}\_2(\mathbf{x} - \mathbf{x}\_k) + \mathbf{C}\_3(\mathbf{x} - \mathbf{x}\_k)^2 + \mathbf{C}\_4(\mathbf{x} - \mathbf{x}\_k)^3 \tag{16}$$

where *C*1, *C*2, *C*3, and *C*4 are polynomial coefficients, which can be calculated as [33]

$$\begin{cases} C\_1 = \mathbf{y}\_k \\ C\_2 = \mathbf{g}\_k \\ C\_3 = \frac{\Im\left(\mu\_k - 2\mathbf{g}\_k - \mathbf{g}\_{k+1}\right)}{\mathbf{x}\_{k+1} - \mathbf{x}\_k} \\ C\_4 = \frac{\mathbf{g}\_k + 2\mathbf{g}\_{k+1} - 2\mathbf{u}\_k}{\left(\mathbf{x}\_{k+1} - \mathbf{x}\_k\right)^2} \end{cases} \tag{17}$$

**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**

=+ - + - + - 2 3

( )

+

1

(17)

*k kk k k kk k k k*

2 2

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

( )

+ + +

1 1 4 2 1

where *C*1, *C*2, *C*3, and *C*4 are polynomial coefficients, which can be calculated as [33]

*u gg C x x*

3 2

ï - - í = - ï

ï + - ï = ï - î

*gg u C x x*

(*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)

*k k*

*C y C g*

ì = ï = ï ï

1 2

3

**Table 4.** Coordinates of control points of three test cases.

56 Optimization Algorithms- Methods and Applications

**Figure 5.** Curve fitting results by means of Akima cubic interpolation.

12 3 4 () ( ) ( ) ( ) *Yx C C x x C x x C x x kk k* (16)

The boundary shape changes with the optimization of the design surface and the computa‐ tional 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(\mathbf{s}, \hat{\mathbf{s}})}{\partial s} = -\beta\_{\mathbf{e}} I(\mathbf{s}, \hat{\mathbf{s}}) + \kappa\_{\mathbf{a}} I\_{\mathbf{b}}(\mathbf{s}) + \frac{\kappa\_{\mathbf{s}}}{4\pi} \int\_{4\pi} I(\hat{\mathbf{s}}\_{i}) \Phi\left(\hat{\mathbf{s}}\_{i}, \hat{\mathbf{s}}\right) d\Omega\_{i} \tag{18}$$

which is an integro-differential type, where *I* is the function of position *s* and 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 **ŝ***i* and scattering direction **ŝ**, which can be defined as *Φ* =1.0 + *a* cos (ŝ*i* ⋅ ŝ). 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 \mathbf{x}} + \beta^{m}\frac{\partial I^{m}}{\partial y} = -\beta\_{\text{e}}I^{m} + \kappa\_{\text{a}}I\_{\text{b}}(\mathbf{s}) + \frac{\kappa\_{\text{s}}}{4\pi} \left(\sum\_{m'=1}^{N\Omega} w^{m'}I\_{\text{w}}^{m'} \Phi^{m',m}\right) \tag{19}$$

where *αm* and *βm* are the direction cosines of the discrete direction *m*. *s* refers to the spatial position, *w* is the quadrature weight. The radiative boundary condition can be directly imposed as follows [34]

$$I\_{\mathbf{w}}^{\rm m} = \mathcal{E}\_{\mathbf{w}} I\_{\mathbf{b}, \mathbf{w}} + \frac{1 - \mathcal{E}\_{\mathbf{w}}}{\pi} \sum\_{\mathbf{n}\_{\nu} \cdot \mathbf{s}^{\rm m} < 0} \varpi^{\rm m} I\_{\mathbf{w}}^{\rm m} \left| \mathbf{n}\_{\mathbf{w}} \cdot \mathbf{s}^{\rm m} \right| \; \quad \mathbf{n}\_{\mathbf{w}} \cdot \mathbf{s}^{\rm m} > 0 \tag{20}$$

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 **s** denotes 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(\mathbf{x}, \mathbf{y})}{\partial(\boldsymbol{\xi}, \boldsymbol{\eta})} \tag{21}$$

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

$$\begin{split} &I \left[ I^{m} \mathcal{I} \left( \boldsymbol{\alpha}^{m} \frac{\partial \boldsymbol{\xi}}{\partial \boldsymbol{\alpha}} + \boldsymbol{\beta}^{m} \frac{\partial \boldsymbol{\xi}}{\partial \boldsymbol{\eta}} \right) \right] \Delta \boldsymbol{\eta} - \left[ I^{m} \mathcal{I} \left( \boldsymbol{\alpha}^{m} \frac{\partial \boldsymbol{\xi}}{\partial \boldsymbol{\alpha}} + \boldsymbol{\beta}^{m} \frac{\partial \boldsymbol{\xi}}{\partial \boldsymbol{\eta}} \right) \right]\_{\mathbf{w}} \Delta \boldsymbol{\eta} + \left[ I^{m} \mathcal{I} \left( \boldsymbol{\alpha}^{m} \frac{\partial \boldsymbol{\eta}}{\partial \boldsymbol{\alpha}} + \boldsymbol{\beta}^{m} \frac{\partial \boldsymbol{\eta}}{\partial \boldsymbol{\eta}} \right) \right]\_{\mathbf{n}} \Delta \boldsymbol{\xi} \\ & - \left[ I^{m} \mathcal{I} \left( \boldsymbol{\alpha}^{m} \frac{\partial \boldsymbol{\eta}}{\partial \boldsymbol{\alpha}} + \boldsymbol{\beta}^{m} \frac{\partial \boldsymbol{\eta}}{\partial \boldsymbol{\eta}} \right) \right]\_{\mathbf{s}} \Delta \boldsymbol{\xi} = \boldsymbol{J} \left[ -\mathcal{B}\_{\mathbf{e}} I^{m} + \kappa\_{\mathbf{a}} \mathbf{I}\_{\mathbf{b}} + \frac{\kappa\_{\mathbf{s}}}{4\pi} \left( \sum\_{m'=1}^{N\Omega} \boldsymbol{w}^{m'} I\_{\mathbf{z}}^{m'} \boldsymbol{\Phi}^{m',m} \right) \right]\_{\mathbf{p}} \Delta \boldsymbol{\xi} \Delta \boldsymbol{\eta} \end{split} \tag{22}$$

where *P* represents the central node of the control volume. The subscripts *e*, *w*, *n*, and *s* represent the eastern, western, northern and southern boundaries around *P*, respectively.

The step scheme is applied to solve the above equations, and Eq.  (24) can be expressed in matrix form

Inverse Geometry Design of Radiative Enclosures Using Particle Swarm Optimization Algorithms http://dx.doi.org/10.5772/62351 59

$$\mathbf{A}\,\Psi = \mathbf{B} \tag{25}$$

where **A** represents the five-diagonal non-symmetric coefficient matrix, **Ψ** represents the vector that consists of the variables *I m* at grid nodes, and **B** represents the vector that consists of the variables *bm* on 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 conver‐ gence rate [35].

**Figure 6.** The schematic diagram of grids in computational domains.

x h

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

 x

k

*s m P*

where *P* represents the central node of the control volume. The subscripts *e*, *w*, *n*, and *s* represent 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

=+ +++ *mm mm m m mm mm m*

where the subscripts *E*, *W*, *N*, and *S* represent the central node of eastern, western, northern,

ì éù éù ü ì ü ï æö æö ¶ ¶ ï ï ¶ ¶ ï =- + D = + D í êú êú ç÷ ç÷ ý í ý êú êú ¶ ¶ ¶ ¶ ïî ëû ëû èø èø ï ï þ î ïþ

max , 0 max , 0

*x y x y*

*x y x y*

ì éù éù ü ì ü ï æö æö ¶ ¶ ï ï ¶ ¶ ï =- + D = + D í êú êú ç÷ ç÷ ý í ý êú êú ¶ ¶ ¶ ¶ ïî ëû ëû èø èø ï ï þ î

*x y x y*

*x y x y*

x h

ìé ùé ù ü ì ü ï æö æö ¶ ¶ ï ï ¶ ¶ ï = + D+ - + D íê úê ú ç÷ ç÷ ý í ý ê úê ú ¶ ¶ ¶ ¶ ïîë ûë û èø èø þ î ï ï þï

max , 0 max , 0

ìé ùé ù ü ì ü ï æö æö ¶ ¶ ï ï ¶ ¶ ï + + D+- + D íê úê ú ç÷ ç÷ ý í ý ê úê ú ¶ ¶ ¶ ¶ ïîë ûë û èø èø ï ï þ î ïþ

*m m m m*

max , 0 max , 0

max , 0 max , 0

p

' 1 4

 h

å s ' ' ',

*P P E E WW NN S S P aI aI a I aI aI b* (23)

*e w*

*n s*

*e w*

*n s*

x

h

ab

ab

x

h

ab

 x

> h

ab

 x

 h

> h

> > x

 h

 x

ïþ

(24)

*w*

W

*ew n*

=

*x y J* (21)

h

ab

x h

 h

 x

(22)

¶ = ¶ (,) (,)

x

é ù æ öæ ö æ ö ¶ ¶ é ù ¶ ¶ é ù ¶ ¶ ê ú ç ÷ç ÷ ç ÷ + D- + D+ + D ê ú ê ú ê ú ¶ ¶ ê ú ¶ ¶ ê ú ¶ ¶ ë û è øè ø è ø ë û ë û

*x y x y x y*

ab

*mm m mm m mm m*

e ab

*I J I J I J*

*N mm m m m m m m*

é ù æ ö ¶ ¶ é ù æ ö -ê ú ç ÷ + D= - + + ê ú ç ÷ F DD ë û è ø ¶ ¶ ê ú ë û è ø

 bk

*I J J I I wI*

and southern control volumes around control volume *P*, and

 x

*a J a J*

 h

*a J a J*

 x

*a J J*

W

é ù æ ö = + ê ú ç ÷ F DD ë û è ø å

=

' 1 4

 h

s ' ' ',

 h

*m mm m mm*

 x

*m mm m mm*

 h

*m mm m m*

 x

*J J*

*m P*

x

*E W*

*N S*

h

ab

x

ab

h

ab

k

*N m P m m m m P P P*

 p

Eq. (24) can be expressed in matrix form

b

k

a

,

*P*

x h

D D

s

p

4

*T b J w I*

+ e

*J*

*P*

*P*

ab

x

h

 b

*x y*

a

ab

 x

58 Optimization Algorithms- Methods and Applications

 h  h

> x

**Figure 7.** The dimensionless radiative heat flux distribution on the bottom boundary for different absorption coeffi‐ cients.

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 *Tg*, 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* × *L y* =1.0×1.0 m. The absorption and scattering coefficients of media are set as *κ*a =2.0 m-1 and *κ*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.

Inverse Geometry Design of Radiative Enclosures Using Particle Swarm Optimization Algorithms http://dx.doi.org/10.5772/62351 61

0K and all the boundaries are assumed as blockbody. The temperature of media is set as *Tg*, 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

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* × *L y* =1.0×1.0 m. The absorption and scattering coefficients of media are set as *κ*a =2.0 m-1 and *κ*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

program for solving radiative problems in irregular enclosures.

**4.4. Inverse design results and discussions**

60 Optimization Algorithms- Methods and Applications

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.

**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 *Nd*=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 *Nd*=3, 5 and 7, whereas the homogenization degree is relatively poor in the case that *Nd*=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 *Nd*=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].

**Table 6.** Comparison of inverse design results for different number of control points.

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

**Figure 9.** (a) Geometry shape of the design surface and (b) dimensionless radiative heat flux distribution on the design

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 *Nd*=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 *Nd*=3,

surface by means of SPSO algorithm for different numbers of control points.

performance than the initial PSO in terms of computational accuracy and efficiency.

62 Optimization Algorithms- Methods and Applications

**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 *ω* =0.5, the extinction coefficients are set as *β* =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 *β* =3.0, the scattering albedos are set as *ω* =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.

**Table 7.** Inverse geometry design results for different extinction coefficients.

Inverse Geometry Design of Radiative Enclosures Using Particle Swarm Optimization Algorithms http://dx.doi.org/10.5772/62351 65

**Table 8.** Inverse geometry design results 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 *ω* =0.5, the extinction coefficients are set as *β* =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 *β* =3.0, the scattering albedos are set as *ω* =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**

64 Optimization Algorithms- Methods and Applications

**Fitness values**

**Dimensionless radiative heat flux**

PSO *β*e=1.0 8.22×10-8 0.2333 0.0329 0.0895

SPSO *β*e=1.0 7.61×10-8 0.2334 0.0295 0.0816

DEPSO *β*e=1.0 5.07×10-8 0.2334 0.0230 0.0849

**Table 7.** Inverse geometry design results for different extinction coefficients.

*β*e=3.0 7.19×10-8 0.1137 0.0496 0.0736

*β*e=5.0 3.94×10-8 0.0735 0.0460 0.1297

*β*e=3.0 5.63×10-8 0.1136 0.0398 0.0723

*β*e=5.0 2.42×10-8 0.0734 0.0301 0.0909

*β*e=3.0 4.16×10-8 0.1136 0.0351 0.0706

*β*e=5.0 2.85×10-8 0.0734 0.0307 0.0889

**Average relative**

**Maximum**

**relative error (%)**

**error (%)**

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 *β* =3.0 and *ω* =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.

**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 *Nd*=3.

#### **Acknowledgements**

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 acknowl‐ edged. A very special acknowledgement is also made to the editors and referees who make important comments to improve this chapter.

#### **Author details**

**Algorithm Scattering**

**5. Conclusions**

points is recommended as *Nd*=3.

important comments to improve this chapter.

**Acknowledgements**

**characteristic**

66 Optimization Algorithms- Methods and Applications

**Fitness values**

**Dimensionless radiative heat flux**

*a* = 0 9.04×10-8 0.1136 0.0403 0.0776 *a* = − 1 7.21×10-8 0.1006 0.0385 0.0941

*a* = 0 8.16×10-8 0.1136 0.0463 0.0825 *a* = − 1 1.60×10-8 0.1007 0.0195 0.0629

*a* = 0 6.94×10-8 0.1136 0.0426 0.0735 *a* = − 1 5.71×10-8 0.1007 0.0413 0.0882

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

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 acknowl‐ edged. A very special acknowledgement is also made to the editors and referees who make

PSO *a* = 1 8.66×10-8 0.1297 0.0354 0.0828

SPSO *a* = 1 7.57×10-8 0.1297 0.0328 0.0982

DEPSO *a* = 1 5.53×10-8 0.1297 0.0378 0.0797

**Table 9.** Inverse geometry design results for different scattering characteristics of media.

**Average relative error (%)**

**Maximum relative error (%)**

> Hong Qi\* , Shuang-Cheng Sun, Zhen-Zong He, Shi-Ting Ruan, Li-Ming Ruan and He-Ping Tan

\*Address all correspondence to<EMAIL_ADDRESS>
School of Energy Science and Engineering, Harbin Institute of Technology, Harbin, China

#### **References**

[23] Moparthi A, Das R, Uppaluri R, Mishra SC Optimization of heat fluxes on the heater and the design surfaces of a radiating-conducting medium. Numerical Heat Transfer, Part A. 2009, 56(10): 846-860.

[11] Daun K, França F, Larsen M, et al. Comparison of methods for inverse design of radiant

[12] Franca FHR, Howell JR Transient inverse design of radiative enclosures for thermal processing of materials. Inverse Problems in Science and Engineering. 2006, 14(4):

[13] Tan JY, Zhao JM, Liu L.H Geometric optimization of a radiation-conduction heating device using meshless method. International Journal of Thermal Sciences. 2011, 50:

[14] Sarvari SMH, Mansouri SH Inverse design for radiative heat source in two-dimensional participating media. Numerical Heat Transfer, Part B. 2004, 46(3): 283-300.

[15] Qi H, Ruan LM, Zhang HC, Wang YM., Tan HP Inverse radiation analysis in a onedimensional participating slab by the stochastic particle swarm optimizer algorithm.

[16] Qi H, Ruan LM, Shi M, An W, Tan HP Application of multi-phase particle swarm optimization technique to inverse radiation problem. Journal of Quantitative Spectro‐

[17] Das R, Mishra SC, Ajith M, Uppaluri R An inverse analysis of a transient 2-D conduc‐ tion-radiation problem using the lattice Boltzmann method and the finite volume method coupled with the genetic algorithm. Journal of Quantitative Spectroscopy and

[18] Mirsephai A, Mohammadzaheri M, Chen L, O'Neill B An artificial intelligence ap‐ proach to inverse heat transfer modeling of an irradiative dryer. International Com‐

[19] Safavinejada A, Mansouria SH, Sakuraib A, Maruyamab S Optimal number and location of heaters in 2-D radiant enclosures composed of specular and diffuse surfaces using micro-genetic algorithm. Applied Thermal Engineering. 2009, 29(5-6): 1075-1085.

[20] Zhang Biao B, Qi Hong H, Ren YT, Sun SC, Ruan LM Application of homogenous continuous ant colony optimization algorithm to inverse problem of one-dimensional coupled radiation and conduction heat transfer. International Journal of Heat and Mass

[21] Ren YT, Qi H, Huang X, Wang W, Ruan LM, Tan HP Application of improved krill herd algorithms to inverse radiation problems. International Journal of Thermal

[22] He ZZ, Qi H, Yao YC, Ruan LM Inverse estimation of the particle size distribution using the fruit fly optimization algorithm. Applied Thermal Engineering. 2015, 88: 306-314.

enclosures. Journal of Heat Transfer. 2005, 128(3): 269-282.

International Journal of Thermal Sciences. 2007, 46: 649-661.

munications in Heat and Mass Transfer. 2012, 39(1): 40-45.

Sciences. 2016. DOI:10.1016/j.ijthermalsci.2015.12.009.

scopy and Radiative Transfer. 2008, 109: 476-493.

Radiative Transfer. 2008, 109(11): 2060-2077.

Transfer, 2013, 66(4): 507-516.

423-436.

68 Optimization Algorithms- Methods and Applications

1820-1831.

## **Shape Optimization of Busemann-Type Biplane Airfoil for Drag Reduction Under Non-Lifting and Lifting Conditions Using Genetic Algorithms**

Yi Tian and Ramesh K. Agarwal

Additional information is available at the end of the chapter

http://dx.doi.org/10.5772/62811

#### **Abstract**

The focus of this chapter is on the shape optimization of the Busemann-type biplane airfoil for drag reduction under both non-lifting and lifting conditions using genetic algorithms. The concept of biplane airfoil was first introduced by Adolf Busemann in 1935. Under design conditions at a specific supersonic flow speed, the Busemann biplane airfoil eliminates all wave drag due to its symmetrical biplane configuration; however, it produces zero lift. Previous research has shown that the original Buse‐ mann biplane airfoil shows a poor performance under off-design conditions. In order to address this problem of zero lift and to improve the off-design-condition perform‐ ance, shape optimization of an asymmetric biplane airfoil is performed. The commer‐ cially available computational fluid dynamics (CFD) solver ANSYS FLUENT is employed for computing the inviscid supersonic flow past the biplane airfoil. A singleobjective genetic algorithm (SOGA) is employed for shape optimization under the nonlifting condition to minimize the drag, and a multi-objective genetic algorithm (MOGA) is used for shape optimization under the lifting condition to maximize both the lift and the lift-to-drag ratio. The results obtained from both SOGA and MOGA show a significant improvement in the design and off-design-condition performance of the optimized Busemann biplane airfoil compared to the original airfoil.

**Keywords:** genetic algorithm, shape optimization, biplane airfoil, computational fluid dynamics, supersonic flow

© 2016 The Author(s). Licensee InTech. This chapter is distributed under the terms of the Creative Commons Attribution License (http://creativecommons.org/licenses/by/3.0), which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.

#### **1. Introduction**

For decades, the speed of commercial aircrafts has been bounded by the sound barrier. Even the most successful supersonic transport (SST) plane, the Concorde, could only be deployed in very few routes due to government regulation, low efficiency, and excessive noise genera‐ tion. Since the retirement of Concorde in 2003, the desire for developing a replacement for Concorde space still remains. In order to accomplish that, one of the biggest design challeng‐ es is to eliminate, or at least greatly reduce, the strong bow shock wave generated during supersonic flight, which can cause high wave drag and substantial noise. At supersonic speed, a bow shock is generated ahead of the airplane. This shock wave generates a substantially high wave drag which needs to be overcome by the engine by providing a much higher thrust compared to that for a conventional subsonic/transonic airplane, which results in higher fuel consumption and low propulsive efficiency. In order to address the high-wave-drag prob‐ lem at supersonic speed, a biplane concept was proposed by Adolf Busemann in 1935 [1] which can potentially avoid the formation of bow shock and thus does not create sonic boom.

During the period from 1935 to 1960, substantial research was conducted on the Busemann biplane concept. Moeckel [2] and Licher [3] performed the theoretical analysis of the optimized lifting supersonic biplanes. Tan [4] took a further step and derived analytical expressions for the drag and lift of a three-dimensional supersonic biplane with a finite span of rectangular plan form. Some experimental results were obtained by Ferri [5] using a wind tunnel, and comparisons were made between the experimental and analytical results. During the past 10 years, considerable interest has been shown in supersonic biplane airfoils. Igra and Arad [6] tested and analyzed the effects of different parameters on the drag coefficient of the Busemann airfoil under various flow conditions. Kusunose et al. [7] proposed a concept for the nextgeneration supersonic transport using the Busemann biplane design. A series of studies using both computational fluid dynamics (CFD) and wind-tunnel experimental methods have been completed by Kusunose's research group [7–16].

Recently, Hu et al. [17] employed a multi-point adjoint-based aerodynamic design and optimization method to improve the baseline Busemann biplane airfoil's off-design perform‐ ance and alleviate the flow-hysteresis problem. They also addressed the problem of minimiz‐ ing its drag by shape optimization. This chapter addresses the aerodynamic shape optimization of two-dimensional symmetric Busemann biplane airfoil under both design and off-design conditions for reducing its wave drag, and the alleviation of the flow-hysteresis and choked-flow effects using a single-objective genetic algorithm. The symmetric Busemann biplane airfoil generates zero lift; this problem is addressed by introducing asymmetry in the shape of two profiles of the biplane. The asymmetric configuration is shape optimized for maximum lift and minimum drag by employing a multi-objective genetic algorithm. The flow field is computed using the commercial CFD software ANSYS FLUENT. Body-fitted H-grids around the airfoils are generated using the ICEM software. Random airfoil shapes with constraints in a given generation of the genetic algorithm are obtained using the Bezier curves (third-order polynomials).

### **2. Flow field simulation of the standard diamond-shaped airfoil and the baseline Busemann biplane airfoil under design and off-design conditions**

The flow fields of both the standard diamond-shaped airfoil and the Busemann biplane airfoil are computed at zero angle of attack (the zero-lift condition). For meaningful comparison, the total thickness of the two airfoils is set to be equal. For this specific case under consideration, the thickness-to-chord ratio of the diamond-shaped airfoil is *t*/*c* = 0.1, while the thickness-tocho rd ratio of the Busemann airfoil is *t*/*c* = 0.05 for both its upper and lower components. The distance between the upper and lower components of the Busemann airfoil is set to be one half of the chord length to obtain the theoretical minimum drag under the design condition at Mach number *M* = 1.7. The angle of attack is set to be 0 as we first consider the non-lifting design and off-design conditions by varying the Mach number below 1.7 and greater than 1.7.

#### **2.1. Mesh generation**

**1. Introduction**

72 Optimization Algorithms- Methods and Applications

For decades, the speed of commercial aircrafts has been bounded by the sound barrier. Even the most successful supersonic transport (SST) plane, the Concorde, could only be deployed in very few routes due to government regulation, low efficiency, and excessive noise genera‐ tion. Since the retirement of Concorde in 2003, the desire for developing a replacement for Concorde space still remains. In order to accomplish that, one of the biggest design challeng‐ es is to eliminate, or at least greatly reduce, the strong bow shock wave generated during supersonic flight, which can cause high wave drag and substantial noise. At supersonic speed, a bow shock is generated ahead of the airplane. This shock wave generates a substantially high wave drag which needs to be overcome by the engine by providing a much higher thrust compared to that for a conventional subsonic/transonic airplane, which results in higher fuel consumption and low propulsive efficiency. In order to address the high-wave-drag prob‐ lem at supersonic speed, a biplane concept was proposed by Adolf Busemann in 1935 [1] which can potentially avoid the formation of bow shock and thus does not create sonic boom.

During the period from 1935 to 1960, substantial research was conducted on the Busemann biplane concept. Moeckel [2] and Licher [3] performed the theoretical analysis of the optimized lifting supersonic biplanes. Tan [4] took a further step and derived analytical expressions for the drag and lift of a three-dimensional supersonic biplane with a finite span of rectangular plan form. Some experimental results were obtained by Ferri [5] using a wind tunnel, and comparisons were made between the experimental and analytical results. During the past 10 years, considerable interest has been shown in supersonic biplane airfoils. Igra and Arad [6] tested and analyzed the effects of different parameters on the drag coefficient of the Busemann airfoil under various flow conditions. Kusunose et al. [7] proposed a concept for the nextgeneration supersonic transport using the Busemann biplane design. A series of studies using both computational fluid dynamics (CFD) and wind-tunnel experimental methods have been

Recently, Hu et al. [17] employed a multi-point adjoint-based aerodynamic design and optimization method to improve the baseline Busemann biplane airfoil's off-design perform‐ ance and alleviate the flow-hysteresis problem. They also addressed the problem of minimiz‐ ing its drag by shape optimization. This chapter addresses the aerodynamic shape optimization of two-dimensional symmetric Busemann biplane airfoil under both design and off-design conditions for reducing its wave drag, and the alleviation of the flow-hysteresis and choked-flow effects using a single-objective genetic algorithm. The symmetric Busemann biplane airfoil generates zero lift; this problem is addressed by introducing asymmetry in the shape of two profiles of the biplane. The asymmetric configuration is shape optimized for maximum lift and minimum drag by employing a multi-objective genetic algorithm. The flow field is computed using the commercial CFD software ANSYS FLUENT. Body-fitted H-grids around the airfoils are generated using the ICEM software. Random airfoil shapes with constraints in a given generation of the genetic algorithm are obtained using the Bezier curves

completed by Kusunose's research group [7–16].

(third-order polynomials).

The commercial meshing software ANSYS ICEM CFD is used to generate the mesh for computing the flow field. The far-field boundary is set at 21 chord length by a 20.5-chordlength rectangle from the center of the airfoil. **Figure 1** shows the H-mesh configuration generated for the simulations. There are 64 nodes around the two components both horizon‐ tally and vertically. In between the two components, the grid has a dimension of 32 nodes (vertically) × 64 nodes (horizontally). An ICEM replay script file is created to automatically generate the mesh for the flow past different airfoil shapes once it is called.

**Figure 1.** (a) H-mesh generated in the computational domain and (b) zoomed-in view of the mesh near the biplane airfoil.

#### **2.2. Steady-state flow field simulation at** *M* **< 1.7**

As has been mentioned above, the baseline Busemann biplane airfoil has a much lower drag under its design condition at Mach number *M* = 1.7. However, it shows a poor aerodynamic performance under off-design conditions, which may cause much higher drag compared to the standard diamond-shaped airfoil. The configurations of both the standard diamondshaped airfoil and the baseline Busemann biplane airfoil are shown in **Figure 2**.

**Figure 2.** Configurations of standard diamond-shaped airfoil (left) and baseline Busemann biplane airfoil (right).

**Figure 3** shows the drag coefficient *cd* of both the standard diamond-shaped airfoil and the baseline Busemann biplane airfoil under the zero-lift condition over a range of Mach numbers (0.3 ≤ *M* ≤ 3.3). The simulations are performed using ANSYS FLUENT; the flow field is initialized with an impulsive uniform flow.

**Figure 3.** Comparison of *Cd* for two different airfoils under non-lifting condition.

**2.2. Steady-state flow field simulation at** *M* **< 1.7**

74 Optimization Algorithms- Methods and Applications

initialized with an impulsive uniform flow.

**Figure 3.** Comparison of *Cd* for two different airfoils under non-lifting condition.

As has been mentioned above, the baseline Busemann biplane airfoil has a much lower drag under its design condition at Mach number *M* = 1.7. However, it shows a poor aerodynamic performance under off-design conditions, which may cause much higher drag compared to the standard diamond-shaped airfoil. The configurations of both the standard diamond-

shaped airfoil and the baseline Busemann biplane airfoil are shown in **Figure 2**.

**Figure 2.** Configurations of standard diamond-shaped airfoil (left) and baseline Busemann biplane airfoil (right).

**Figure 3** shows the drag coefficient *cd* of both the standard diamond-shaped airfoil and the baseline Busemann biplane airfoil under the zero-lift condition over a range of Mach numbers (0.3 ≤ *M* ≤ 3.3). The simulations are performed using ANSYS FLUENT; the flow field is

As shown in **Figure 3**, the baseline Busemann biplane airfoil has a higher drag compared to the standard diamond-shaped airfoil when the Mach number is low (0.3 ≤ *M* ≤ 1.6). In the range of 1.7 ≤ *M* ≤ 2.7, however, the drag generated by the baseline Busemann biplane airfoil is smaller than that generated by the standard diamond-shaped airfoil, especially at *M* = 1.7 which is the design condition for the Busemann airfoil. **Figure 3** demonstrates the advantage of the Busemann biplane airfoil as it produces much lower drag near its design condition (*M* = 1.7) due to the wave-reduction and wave-cancellation effects. The simula‐ tions shown above demonstrate that the baseline Busemann biplane airfoil shows a very good performance under its design conditions (*M* = 1.7), while it has a much higher drag under off-design conditions due to the choked-flow phenomenon at lower Mach numbers (*M* < 1.7). **Figure 4(a)** through **4(d)** shows the flow field around the baseline Busemann bi‐ plane airfoil under two different off-design conditions, both of which result in high drag but due to different reasons. At *M* = 0.8, as shown in **Figure 4(a)**, the flow between the two airfoil components reaches Mach 1 at the mid-chord of the airfoil, and then increases fur‐ ther downstream to supersonic flow. After reaching the trailing edge of the airfoil, howev‐ er, the flow speed drops to subsonic and hence forms a vertical shock wave. The supersonic flow in the rear part of the airfoil (after the mid-point) creates a low-pressure region, as shown in **Figure 4(b)**, which leads to a higher wave drag (*cd* = 0.148) at this offdesign condition of *M* = 0.8. At *M* = 1.6, however, as shown in **Figure 4(c)**, the flow field is very different from that at *M* = 0.8 under subsonic conditions. A bow shock is formed in front of the leading edge of the airfoil. The flow speed drops to subsonic and a high-pres‐ sure region, as shown in **Figure 4(d)**, is created which again generates a high drag (*cd* = 0.0926).

**Figure 4.** Velocity and pressure contours around the baseline Busemann airfoil: (a) velocity contours at *M* = 0.8, (b) pressure contours at *M* = 0.8, (c) velocity contours at *M* = 1.6, (d) pressure contours at *M* = 1.6.

#### **2.3. Flow field simulation during acceleration and deceleration**

Despite the high drag being generated under off-design conditions, there is an even worse problem caused by the flow-hysteresis phenomenon during acceleration and choked-flow phenomenon during deceleration that need to be addressed under off-design conditions. To demonstrate these phenomena, flow field simulations are conducted under both acceleration and deceleration using the previous simulation results shown in **Figure 4** as the initial condition. **Figure 5** shows two separated *cd* curves during acceleration and deceleration of the biplane airfoil. As shown in **Figure 5**, the change in acceleration and deceleration *cd* curves occurs in the range *M* = 1.6 to *M* = 2.1, where the blue dashed line represents the *cd* of the Busemann biplane airfoil during acceleration and the red solid line represents the *cd* of the Busemann biplane airfoil during deceleration. This separation between the two *cd* curves is caused by the flow-hysteresis phenomenon during acceleration and the choked-flow phe‐ nomenon during deceleration. Therefore, in order to minimize the difference in the drag coefficients of the Busemann biplane airfoil during acceleration and deceleration as well as to significantly decrease the drag, shape optimization of the baseline Busemann airfoil is conducted using single-objective genetic algorithm (SOGA) to eliminate the flow-hysteresis phenomenon during acceleration and the choked-flow phenomenon during deceleration; it is discussed in the next section.

**Figure 5.** *cd* plots for diamond-shaped airfoil and baseline Busemann airfoil during acceleration and deceleration under non-lifting conditions.

#### **2.4. Flow field of the baseline Busemann biplane airfoil during acceleration**

**2.3. Flow field simulation during acceleration and deceleration**

76 Optimization Algorithms- Methods and Applications

discussed in the next section.

non-lifting conditions.

Despite the high drag being generated under off-design conditions, there is an even worse problem caused by the flow-hysteresis phenomenon during acceleration and choked-flow phenomenon during deceleration that need to be addressed under off-design conditions. To demonstrate these phenomena, flow field simulations are conducted under both acceleration and deceleration using the previous simulation results shown in **Figure 4** as the initial condition. **Figure 5** shows two separated *cd* curves during acceleration and deceleration of the biplane airfoil. As shown in **Figure 5**, the change in acceleration and deceleration *cd* curves occurs in the range *M* = 1.6 to *M* = 2.1, where the blue dashed line represents the *cd* of the Busemann biplane airfoil during acceleration and the red solid line represents the *cd* of the Busemann biplane airfoil during deceleration. This separation between the two *cd* curves is caused by the flow-hysteresis phenomenon during acceleration and the choked-flow phe‐ nomenon during deceleration. Therefore, in order to minimize the difference in the drag coefficients of the Busemann biplane airfoil during acceleration and deceleration as well as to significantly decrease the drag, shape optimization of the baseline Busemann airfoil is conducted using single-objective genetic algorithm (SOGA) to eliminate the flow-hysteresis phenomenon during acceleration and the choked-flow phenomenon during deceleration; it is

**Figure 5.** *cd* plots for diamond-shaped airfoil and baseline Busemann airfoil during acceleration and deceleration under

In this section, the flow-hysteresis phenomenon during acceleration is examined for the baseline Busemann biplane airfoil under the non-lifting condition. The pressure coefficient contours for flow past the baseline Busemann biplane airfoil during acceleration are shown in **Figure 6**. To simulate the non-lifting condition, the angle of attack is set to 0. The resulting nonlifting flow field around the baseline Busemann airfoil is shown in **Figure 6** at various supersonic Mach numbers ranging from *M* = 1.7 to *M* = 2.13. As shown in **Figure 6**, the bow shock exists in front of the Busemann airfoil and does not disappear until the Mach number reaches 2.13. It can also be noticed that there is a subsonic region behind the bow shock between the upper and lower components of the airfoil where the pressure coefficients are high. The presence of this bow shock in front of the airfoil results in a substantial increase in the drag compared to that at the design condition. However, when the Mach number increases from 2.12 to 2.13, the bow shock is swallowed into the upper and lower components of the Busemann airfoil and is replaced by two oblique shock waves, and the subsonic region between the two airfoil components finally disappears as shown in **Figure 6(d)**. The drag coefficient of the airfoil also decreases dramatically, and the flow past the airfoil develops into a state similar to that under the design condition. **Figure 6** illustrates the poor performance of the baseline Busemann airfoil under off-design conditions. During acceleration, the design condition cannot be achieved at Mach 1.7. And the drag coefficient at Mach 1.7 is much higher (*cd* = 0.0923) compared to that of the standard diamond-shaped airfoil (*cd* = 0.0292). Furthermore, due to the flow-hysteresis phenomenon, the drag coefficients *cd* during acceleration and deceleration are different as shown in **Figure 6**.

**Figure 6.** *cp*-contours around the baseline Busemann biplane airfoil with zero lift during acceleration. (a) *M*∞ = 1.7, *cd* = 0.0923; (b) *M*∞ = 2.0, *cd* = 0.0883; (c) *M*∞ = 2.12, *cd* = 0.0860; (d) *M*∞ = 2.13, *cd* = 0.010264.

#### **2.5. Flow field of the baseline Busemann biplane airfoil during deceleration**

In this section, the choked-flow phenomenon of the baseline Busemann biplane airfoil during deceleration is examined. The contours of pressure coefficients for flow past the baseline Busemann airfoil during deceleration are shown in **Figure 7**. To simulate the non-lifting condition, the angle of attack is set to 0 as before. The resulting non-lifting flow field around the baseline Busemann airfoil is shown in **Figure 7** at various supersonic Mach numbers ranging from *M* = 1.8 to *M* = 1.62. As shown in **Figure 7**, a different flow field appears within the small range near the design Mach number (*M* = 1.7); a relatively high drag coefficient still occurs as the Mach number further decreases during deceleration. A strong bow shock is formed in front of the airfoil when the Mach number drops from 1.63 to 1.62, while the drag coefficient increases dramatically from 0.005594 to 0.0926, which is substantially higher than that of the standard diamond-shaped airfoil (*cd* = 0.03158). The flow between the two compo‐ nents of the airfoil is choked at the location of the maximum thickness of the Busemann airfoil, and a subsonic region is formed. This is also a clear indication of the poor performance of the baseline Busemann biplane airfoil under off-design conditions as the drag coefficient of the baseline Busemann airfoil is much higher than that of the standard diamond-shaped airfoil for *M* ≤ 1.62.

**Figure 7.** *cp*-contours around the baseline Busemann biplane airfoil with zero lift during deceleration: (a) *M*∞ = 1.8, *cd* = 0.003638; (b) *M*∞ = 1.7, *cd* = 0.002183; (c) *M*∞ = 1.65, *cd* = 0.004144; (d) *M*∞ = 1.64, *cd* = 0.004778; (e) *M*∞ = 1.63, *cd* = 0.005594; (f) *M*∞ = 1.62, *cd* = 0.09261.

In conclusion, the baseline Busemann biplane airfoil produces a substantially higher drag in the low Mach number range (below the design Mach number of *M* = 1.7). Additionally, it is necessary to accelerate the baseline Busemann biplane airfoil to a much higher Mach number (*M* = 2.13) to reach the shock-wave-swallowing state, while dramatically producing higher drag and largely decreasing its efficiency, and then decelerate to a lower velocity to achieve the design condition at Mach number *M* = 1.7. As a result, the baseline Busemann biplane airfoil design needs to be modified and optimized to avoid or at least reduce the high drag coefficient caused by the flow-hysteresis phenomenon during acceleration and the choked-flow phe‐ nomenon during deceleration.

### **3. Optimization of Busemann-type biplane airfoil**

**2.5. Flow field of the baseline Busemann biplane airfoil during deceleration**

78 Optimization Algorithms- Methods and Applications

*M* ≤ 1.62.

(f) *M*∞ = 1.62, *cd* = 0.09261.

In this section, the choked-flow phenomenon of the baseline Busemann biplane airfoil during deceleration is examined. The contours of pressure coefficients for flow past the baseline Busemann airfoil during deceleration are shown in **Figure 7**. To simulate the non-lifting condition, the angle of attack is set to 0 as before. The resulting non-lifting flow field around the baseline Busemann airfoil is shown in **Figure 7** at various supersonic Mach numbers ranging from *M* = 1.8 to *M* = 1.62. As shown in **Figure 7**, a different flow field appears within the small range near the design Mach number (*M* = 1.7); a relatively high drag coefficient still occurs as the Mach number further decreases during deceleration. A strong bow shock is formed in front of the airfoil when the Mach number drops from 1.63 to 1.62, while the drag coefficient increases dramatically from 0.005594 to 0.0926, which is substantially higher than that of the standard diamond-shaped airfoil (*cd* = 0.03158). The flow between the two compo‐ nents of the airfoil is choked at the location of the maximum thickness of the Busemann airfoil, and a subsonic region is formed. This is also a clear indication of the poor performance of the baseline Busemann biplane airfoil under off-design conditions as the drag coefficient of the baseline Busemann airfoil is much higher than that of the standard diamond-shaped airfoil for

**Figure 7.** *cp*-contours around the baseline Busemann biplane airfoil with zero lift during deceleration: (a) *M*∞ = 1.8, *cd* = 0.003638; (b) *M*∞ = 1.7, *cd* = 0.002183; (c) *M*∞ = 1.65, *cd* = 0.004144; (d) *M*∞ = 1.64, *cd* = 0.004778; (e) *M*∞ = 1.63, *cd* = 0.005594; In this chapter, the shape optimization procedure for the baseline Busemann biplane airfoil using the genetic algorithms (GAs) and the results of the optimized Busemann-type airfoil under both non-lifting and lifting conditions are presented. The optimization process is established by coupling a single-objective genetic algorithm (SOGA)- or a multi-objective genetic algorithm (MOGA)-based optimization method with the mesh generation software ANSYS-ICEM and the CFD solver ANSYS FLUENT as shown in **Figure 8**.

**Figure 8.** Schematic of information flow in the optimization process.

Every individual (airfoil) in each generation of SOGA/MOGA is represented by a set of control points, which randomly generate the airfoil shape by using the Bezier curves. The mesh generation software ICEM is used to generate a two-dimensional structured mesh around the airfoil as an input to the CFD solver FLUENT, which is then used to calculate the supersonic inviscid flow fields for specific flow conditions. Based on the fitness values of all airfoil shapes in a given generation, SOGA/MOGA is applied to create the next generation of airfoils, and this whole process is repeated until the optimal fitness value is obtained. The airfoil shape that corresponds to the optimal fitness value is the final shape of the optimized airfoil [18].

#### **3.1. Overview of genetic algorithm**

Genetic algorithms (GAs) are a class of stochastic optimization algorithms inspired by the biological evolution [19]. They efficiently exploit historical information to speculate on new offspring with improved performance [20]. For the Busemann biplane airfoil in particular, genetic algorithms are used to generate new shapes that produce much lower drag by minimizing the flow-hysteresis phenomenon during acceleration and the choked-flow phenomenon during deceleration. Generally, The GA employs the following steps to complete the optimization process:

#### *3.1.1. Application of single-objective genetic algorithm (SOGA)*

For the non-lifting condition, when the lift coefficient *cl* = 0, a single-objective genetic algorithm (SOGA) is employed for the shape optimization of the Busemann biplane airfoil. The single objective to be achieved is to minimize the drag coefficient *cd*. Multiple design points (Mach numbers) are used during the optimization process due to the flow-hysteresis and chokedflow phenomena caused by the baseline Busemann biplane airfoil during acceleration and deceleration, respectively. For the multi-point optimization, the fitness function employed is a weighted average of the drag coefficients *cd*, which can be written in the following form:

$$I = \frac{\sum\_{i=1}^{n} w\_i I\_i}{\sum\_{i=1}^{n} w\_i} \tag{1}$$

where *I* is the weighted average drag coefficient, and "*i*" denotes a design point (related to Mach number). For the optimization process under the non-lifting condition, an evenly weighted average of drag coefficients *cd* at different Mach numbers is employed. Therefore, Eq. (1) reduces to

$$I = \frac{\sum\_{i=1}^{n} I\_i}{n} \tag{2}$$

As the wave drag of the Busemann biplane airfoil is much lower when the flow is unchoked, it is highly desirable that the strong bow shock wave in front of the airfoil is swallowed into the area between the upper and lower components of the airfoil before the flow speed approaches the design Mach number (*M* = 1.7). Although a higher weight can be assigned to the most important design Mach number (*M* = 1.7) to produce a lower drag at that design point, previous research conducted by Hu et al. [17] has shown that a slightly higher drag coefficient *cd* is obtained at lower Mach numbers during acceleration and deceleration, as given in **Tables 1** and **2**, respectively, if such an uneven weighting method is used.

**Table 1.** Comparison of *cd* at under zero-lift conditions during acceleration (1 count = 0.0001) [17].

genetic algorithms are used to generate new shapes that produce much lower drag by minimizing the flow-hysteresis phenomenon during acceleration and the choked-flow phenomenon during deceleration. Generally, The GA employs the following steps to complete

**4.** Reproduction: pairs of individuals are picked to produce the offspring, which is often done by roulette wheel sampling. A crossover function is then used to produce the

**6.** Check for convergence: if the current generation has converged, the individual having the best fitness will be returned. Otherwise, the process will be repeated starting from Step 2 until the predefined tolerance criteria for acceptable change in fitness between generations

(SOGA) is employed for the shape optimization of the Busemann biplane airfoil. The single objective to be achieved is to minimize the drag coefficient *cd*. Multiple design points (Mach numbers) are used during the optimization process due to the flow-hysteresis and chokedflow phenomena caused by the baseline Busemann biplane airfoil during acceleration and deceleration, respectively. For the multi-point optimization, the fitness function employed is a weighted average of the drag coefficients *cd*, which can be written in the following form:

> 1 1

*w I*

where *I* is the weighted average drag coefficient, and "*i*" denotes a design point (related to Mach number). For the optimization process under the non-lifting condition, an evenly weighted average of drag coefficients *cd* at different Mach numbers is employed. Therefore,

> 1 *n i i I*

As the wave drag of the Busemann biplane airfoil is much lower when the flow is unchoked, it is highly desirable that the strong bow shock wave in front of the airfoil is swallowed into

*I n*

*n i i i n i i*

*I w* = =

= å

= 0, a single-objective genetic algorithm

å (1)

= = å (2)

**1.** Initialization: randomly generates a group of individuals.

*3.1.1. Application of single-objective genetic algorithm (SOGA)*

For the non-lifting condition, when the lift coefficient *cl*

**2.** Evaluation: evaluates the fitness of each individual generated.

**3.**  Mutation: randomly modifies some small percentage of the population.

the optimization process:

80 Optimization Algorithms- Methods and Applications

offspring.

are met.

Eq. (1) reduces to

**Table 2.** Comparison of *cd* at under zero-lift conditions during deceleration (1 count = 0.0001) [17].

**Table 3.** GA parameters for shape optimization of Busemann biplane airfoil under non-lifting conditions.

Here, we have a total of seven design points in the Mach numbers ranging from *M* = 1.1 to *M* = 1.7, for both acceleration and deceleration. The reason for picking these seven design points is that this range of Mach numbers, *M* = 1.1 to *M* = 1.7, provides the critical region before the Mach number increases to the design condition *M* = 1.7; it is the region in which we want to maintain the drag coefficient as low as possible. Both acceleration and deceleration scenarios are considered in order to reduce the flow-hysteresis phenomenon and the choked-flow phenomenon, respectively.

A code in the MATLAB package is developed and utilized in the optimization process of the airfoil as well as for ICEM meshing and FLUENT flow field calculations. All SOGA parameters are defined based on the GA methodology, as shown in **Table 3**.

#### *3.1.2. Airfoil parameterization*

The shape of the airfoil is generated by using the Bezier curves (third-order polynomials). Bezier curves are frequently used in computer graphics to obtain curves that appear to be reasonably smooth at all scales. One of the main reasons that Bezier curves are used in computer graphics is that they can be efficiently constructed; each Bezier curve is simply defined by a set of control points [21].

For the non-lifting case, with *cl* = 0, since the upper and lower components of the airfoil are symmetric with respect to the horizontal axis, only the upper component of the airfoil needs to be defined; the lower component is the mirror image with respect to the horizontal axis. For the upper component of the airfoil, there are two lines that need to be drawn to define its shape. Since the thickness distribution for the entire airfoil remains the same as that of the baseline Busemann biplane airfoil, the *y*-coordinates of the lower line are defined by Eqs. (3) and (4). Thus, the upper line (which is a straight horizontal line for the baseline Busemann airfoil) is the only line that needs to be generated by the Bezier curves in order to define the shape of the whole biplane airfoil.

$$\mathbf{y}\_{\text{low}} = \mathbf{y}\_{\text{up}} - \mathbf{0}.05 - \mathbf{0}.1 \times \mathbf{x}\_{\text{up}} (-\mathbf{0}.\\$ \le \mathbf{x}\_{\text{up}} \le \mathbf{0}) \tag{3}$$

$$\text{l.y}\_{\text{low}} = \text{y}\_{\text{up}} - 0.05 + 0.1 \times \text{x}\_{\text{up}} (0 \le \text{x}\_{\text{up}} \le 0.5) \tag{4}$$

In Eqs. (3) and (4), the subscripts "low" and "up" correspond to the lower and upper lines of the airfoil; the origin is at the center of the upper line. Two Bezier curves are used to generate the shape-defining upper line. Each Bezier curve is defined by a set of four control points. Each control point is constrained by a specified range of *x*- and *y*-coordinates. **Figure 9** shows a randomly generated Busemann-type biplane airfoil shape using Bezier curves.

**Figure 9.** Randomly generated Busemann-type biplane airfoil shape using Bezier curves.

#### *3.1.3. Optimization results*

A code in the MATLAB package is developed and utilized in the optimization process of the airfoil as well as for ICEM meshing and FLUENT flow field calculations. All SOGA parameters

The shape of the airfoil is generated by using the Bezier curves (third-order polynomials). Bezier curves are frequently used in computer graphics to obtain curves that appear to be reasonably smooth at all scales. One of the main reasons that Bezier curves are used in computer graphics is that they can be efficiently constructed; each Bezier curve is simply

symmetric with respect to the horizontal axis, only the upper component of the airfoil needs to be defined; the lower component is the mirror image with respect to the horizontal axis. For the upper component of the airfoil, there are two lines that need to be drawn to define its shape. Since the thickness distribution for the entire airfoil remains the same as that of the baseline Busemann biplane airfoil, the *y*-coordinates of the lower line are defined by Eqs. (3) and (4). Thus, the upper line (which is a straight horizontal line for the baseline Busemann airfoil) is the only line that needs to be generated by the Bezier curves in order to define the shape of the

In Eqs. (3) and (4), the subscripts "low" and "up" correspond to the lower and upper lines of the airfoil; the origin is at the center of the upper line. Two Bezier curves are used to generate the shape-defining upper line. Each Bezier curve is defined by a set of four control points. Each control point is constrained by a specified range of *x*- and *y*-coordinates. **Figure 9** shows a

randomly generated Busemann-type biplane airfoil shape using Bezier curves.

**Figure 9.** Randomly generated Busemann-type biplane airfoil shape using Bezier curves.

= 0, since the upper and lower components of the airfoil are

0.05 0.1 0. )5 0 ( *low up up up yy x x* = - -´ -£ £ (3)

0.05 0.1 0 .5 0( ) *low up up up yy x x* = - +´ £ £ (4)

are defined based on the GA methodology, as shown in **Table 3**.

*3.1.2. Airfoil parameterization*

defined by a set of control points [21].

82 Optimization Algorithms- Methods and Applications

For the non-lifting case, with *cl*

whole biplane airfoil.

After implementing SOGA for 20 generations with eight individuals in each generation, an optimal shape result for symmetric Busemann-type biplane airfoil under the non-lifting condition with a minimum drag is obtained. **Figure 10** shows the geometry of the original Busemann biplane airfoil (red) and the optimized Busemann biplane airfoil (blue) under the non-lifting condition.

**Figure 10.** Geometry of both the original and optimized Busemann biplane airfoil under non-lifting conditions.

**Table 4.** Comparison of *cd* for the original and optimized Busemann biplane airfoil under non-lifting conditions during acceleration (1 count = 0.0001).

**Table 5.** Comparison of *cd* for the original and optimized Busemann biplane airfoil under non-lifting conditions during deceleration (1 count = 0.0001).

The drag coefficients for the seven design points are compared in **Tables 4** and **5** for both the original and optimized Busemann biplane airfoil under the non-lifting condition during acceleration and deceleration, respectively. As shown in **Tables 4** and **5**, the baseline Buse‐ mann biplane airfoil is choked at all Mach numbers within the optimization range, while the optimized Busemann biplane airfoil un-chokes at *M* = 1.5 during acceleration and chokes at *M* = 1.3 during deceleration. Even under choked conditions during both acceleration and deceleration, the optimized Busemann biplane airfoil has significantly lower drag compared to the baseline Busemann biplane airfoil. The only point where the optimized Busemann biplane airfoil has a higher drag compared to the original airfoil is at *M* = 1.7 during deceler‐ ation.

**Figures 11** and **12**, respectively, show the change in the pressure coefficient *cp* around the original Busemann biplane airfoil as the Mach number increases and decreases within the design-point range. The corresponding *cp* contours for the optimized Busemann biplane airfoil are shown in **Figures 13** and **14** to illustrate the wave-cancellation effect.

**Figure 11.** *cp* contours around the original Busemann biplane airfoil during acceleration: (a) *cp* contours at *M* = 1.1; (b) *cp* contours at *M* = 1.3; (c) *cp* contours at *M* = 1.5; (d) *cp* contours at *M* = 1.7.

**Figure 12.** *cp* contours around the original Busemann biplane airfoil during deceleration: (a) *cp* contours at *M* = 1.7; (b) *cp* contours at *M* = 1.5; (c) *cp* contours at *M* = 1.3; (d) *cp* contours at *M* = 1.1.

Shape Optimization of Busemann-Type Biplane Airfoil for Drag Reduction Under Non-Lifting and Lifting Conditions Using Genetic Algorithms http://dx.doi.org/10.5772/62811 85

deceleration, the optimized Busemann biplane airfoil has significantly lower drag compared to the baseline Busemann biplane airfoil. The only point where the optimized Busemann biplane airfoil has a higher drag compared to the original airfoil is at *M* = 1.7 during deceler‐

**Figures 11** and **12**, respectively, show the change in the pressure coefficient *cp* around the original Busemann biplane airfoil as the Mach number increases and decreases within the design-point range. The corresponding *cp* contours for the optimized Busemann biplane airfoil

**Figure 11.** *cp* contours around the original Busemann biplane airfoil during acceleration: (a) *cp* contours at *M* = 1.1; (b) *cp*

**Figure 12.** *cp* contours around the original Busemann biplane airfoil during deceleration: (a) *cp* contours at *M* = 1.7; (b) *cp*

are shown in **Figures 13** and **14** to illustrate the wave-cancellation effect.

contours at *M* = 1.3; (c) *cp* contours at *M* = 1.5; (d) *cp* contours at *M* = 1.7.

contours at *M* = 1.5; (c) *cp* contours at *M* = 1.3; (d) *cp* contours at *M* = 1.1.

ation.

84 Optimization Algorithms- Methods and Applications

**Figure 13.** *cp* contours around the optimized Busemann biplane airfoil during acceleration: (a) *cp* contours at *M* = 1.1; (b) *cp* contours at *M* = 1.3; (c) *cp* contours at *M* = 1.5; (d) *cp* contours at *M* = 1.7.

**Figure 14.** *cp* contours around the optimized Busemann biplane airfoil during deceleration: (a) *cp* contours at *M* = 1.7; (b) *cp* contours at *M* = 1.5; (c) *cp* contours at *M* = 1.3; (d) *cp* contours at *M* = 1.1.

**Figure 15** shows the comparison of the drag coefficients for the standard diamond-shaped airfoil, the baseline Busemann biplane airfoil, and the optimized Busemann biplane airfoil under the non-lifting condition. As shown in this figure, the separation between the accelera‐ tion and deceleration curves for *cd* still exists for the optimized Busemann biplane airfoil, which means that the flow-hysteresis and the choked-flow effects are not totally eliminated. How‐ ever, as clearly shown in **Figure 15**, the flow-hysteresis area has been significantly reduced, and the drag increase during deceleration due to the choked-flow phenomenon is much smaller than that for the original Busemann biplane airfoil. The drag of the optimized Buse‐ mann biplane airfoil in the subsonic region is also smaller than that of the original Busemann biplane airfoil, although it is slightly higher than that of the standard diamond-shaped airfoil with 0.6 < *M* < 0.85. Under both subsonic and supersonic conditions, the optimized Busemann biplane airfoil has been able to significantly reduce the wave drag compared to the original Busemann biplane airfoil. At the design condition Mach number of 1.7, however, the drag coefficient of the optimized Busemann biplane airfoil (*cd* = 0.01038) is much higher than that of the original Busemann biplane airfoil *cd* = 0.002182. This is due to the fact that for the shape optimization, our focus has been on reducing the flow-hysteresis and choked-flow effects, and we chose to assign equal weights to all Mach numbers used as the multiple design points. To address this problem, we could have put more weight on the design condition (*M* = 1.7) during the optimization process.

**Figure 15.** *cd* plot for different airfoils under nonlifting condition.

Next, we examine the details of the flow field during acceleration and deceleration for the optimized Busemann biplane airfoil and compare them with the original baseline Busemann biplane airfoil and the optimization results obtained by Hu et al. [17] using an adjoint-based optimization technique. **Figures 16** and **17** show the pressure coefficient contours of the optimized Busemann biplane airfoil under acceleration and deceleration, respectively. **Figures 18** and **19** show the pressure coefficient contours of the optimized Busemann biplane airfoil using the adjoint-based technique [17] under acceleration and deceleration, respectively.  The swallowing of the bow shock wave happens when the Mach number increases from 1.49 to 1.50 in our GA optimization as shown in **Figure 16(f)**; it happens when the Mach number increases from 1.52 to 1.53 in the adjoint-based optimization as shown in **Figure 18(e)**, and it happens when the Mach number increases from 2.12 to 2.13 for the original Busemann biplane airfoil as shown in **Figure 6(d)**. The drag coefficient decreases from 0.03556 to 0.01316 in the present GA-based optimization and from 0.03336 to 0.01221 in the adjointbased optimization [17].

smaller than that for the original Busemann biplane airfoil. The drag of the optimized Buse‐ mann biplane airfoil in the subsonic region is also smaller than that of the original Busemann biplane airfoil, although it is slightly higher than that of the standard diamond-shaped airfoil with 0.6 < *M* < 0.85. Under both subsonic and supersonic conditions, the optimized Busemann biplane airfoil has been able to significantly reduce the wave drag compared to the original Busemann biplane airfoil. At the design condition Mach number of 1.7, however, the drag coefficient of the optimized Busemann biplane airfoil (*cd* = 0.01038) is much higher than that of the original Busemann biplane airfoil *cd* = 0.002182. This is due to the fact that for the shape optimization, our focus has been on reducing the flow-hysteresis and choked-flow effects, and we chose to assign equal weights to all Mach numbers used as the multiple design points. To address this problem, we could have put more weight on the design condition (*M* = 1.7) during

Next, we examine the details of the flow field during acceleration and deceleration for the optimized Busemann biplane airfoil and compare them with the original baseline Busemann biplane airfoil and the optimization results obtained by Hu et al. [17] using an adjoint-based optimization technique. **Figures 16** and **17** show the pressure coefficient contours of the optimized Busemann biplane airfoil under acceleration and deceleration, respectively. **Figures 18** and **19** show the pressure coefficient contours of the optimized Busemann biplane airfoil using the adjoint-based technique [17] under acceleration and deceleration, respectively.  The swallowing of the bow shock wave happens when the Mach number increases from 1.49 to 1.50 in our GA optimization as shown in **Figure 16(f)**; it happens when the Mach number increases from 1.52 to 1.53 in the adjoint-based optimization as shown in

the optimization process.

86 Optimization Algorithms- Methods and Applications

**Figure 15.** *cd* plot for different airfoils under nonlifting condition.

**Figure 16.** *cp* contours of the GA-optimized Busemann biplane airfoil with zero lift during acceleration: (a) *M* = 1.2, *cd* = 0.05064; (b) *M* = 1.3, *cd* = 0.04275; (c) *M* = 1.4, *cd* = 0.03889; (d) *M* = 1.48, *cd* = 0.03601; (e) *M* = 1.49, *cd* = 0.03556; (f) *M* = 1.50, *cd* = 0.01316.

**Figure 17.** *cp* contours of the GA-optimized Busemann biplane airfoil with zero lift during deceleration: (a) *M* = 1.6, *cd* = 0.01138; (b) *M* = 1.5, *cd* = 0.01316; (c) *M* = 1.4, *cd* = 0.01629; (d) *M* = 1.37, *cd* = 0.02420; (e) *M* = 1.36, *cd* = 0.04045; (f) *M* = 1.3, *cd* = 0.04277.

**Figure 18.** *cp* contours of the adjoint-based-optimized Busemann biplane airfoil with zero lift during acceleration [17]: (a) *M*∞ = 1.3, *cd* = 0.04190; (b) *M*∞ = 1.4, *cd* = 0.03761; (c) *M*∞ = 1.5, *cd* = 0.03318; (d) *M*∞ = 1.52, *cd* = 0.03336; (e) *M*∞ = 1.53, *cd* = 0.0221; (f) *M*∞ = 1.6, *cd* = 0.01125.

Shape Optimization of Busemann-Type Biplane Airfoil for Drag Reduction Under Non-Lifting and Lifting Conditions Using Genetic Algorithms http://dx.doi.org/10.5772/62811 89

**Figure 18.** *cp* contours of the adjoint-based-optimized Busemann biplane airfoil with zero lift during acceleration [17]: (a) *M*∞ = 1.3, *cd* = 0.04190; (b) *M*∞ = 1.4, *cd* = 0.03761; (c) *M*∞ = 1.5, *cd* = 0.03318; (d) *M*∞ = 1.52, *cd* = 0.03336; (e) *M*∞ = 1.53, *cd* =

0.0221; (f) *M*∞ = 1.6, *cd* = 0.01125.

88 Optimization Algorithms- Methods and Applications

**Figure 19.** *cp* contours of the adjoint-based-optimized Busemann biplane airfoil with zero lift during deceleration [17]: (a) *M*∞ = 1.6, *cd* = 0.01125; (b) *M*∞ = 1.5, *cd* = 0.01273; (c) *M*∞ = 1.4, *cd* = 0.01526; (d) *M*∞ = 1.38, = 0.01582; (e) *M*∞ = 1.37, *cd* = 0.03886; (f) *M*∞ = 1.3, *cd* = 0.04191.

In conclusion, as shown in **Figure 20**, the drag coefficient of the GA-optimized Busemann biplane airfoil is significantly reduced when compared to the original Busemann biplane airfoil, and it matches with the adjoint-based optimization result obtained by Hu et al. [17].

**Figure 20.** *cd* plot for different airfoils under non-lifting conditions.

#### **3.2. Shape optimization of Busemann biplane airfoil under lifting conditions**

#### *3.2.1. Application of multi-objective genetic algorithm (MOGA)*

Under the lifting condition with the lift coefficient *cl* ≠ 0, a multi-objective genetic algorithm (MOGA) is employed for shape optimization of the original Busemann biplane airfoil. The two objectives to be achieved are to minimize the drag coefficient *cd* while maximizing the lift coefficient *cl* . Similar to that for the non-lifting condition described in Section 3.2, a total of seven design points ranging from *M* = 1.1 to *M* = 1.7 are used during the optimization process. For the fitness functions for the lifting case, we use the sum of evenly weighted average of both *cd* and *cl* . The GA parameters used for the lifting case are listed in **Table 6**.

**Table 6.** GA parameters for shape optimization of Busemann biplane airfoil under lifting conditions.

#### *3.2.2. Airfoil parameterization*

Similar to that for the non-lifting case, the random shapes of the airfoil are generated by using the Bezier curves with control points. For the lifting case with *cl* ≠ 0, as the upper and lower components of the biplane airfoil are not symmetric, both the upper and lower components need to be defined separately. The thickness distribution for both the upper and lower components is still maintained the same as for the non-lifting case. Now, a total of four Bezier curves are needed to define the shape of the airfoil. **Figure 21** shows a randomly generated Busemann-type biplane airfoil shape using the Bezier curves for lifting conditions.

#### *3.2.3. Optimization results*

**Figure 20.** *cd* plot for different airfoils under non-lifting conditions.

90 Optimization Algorithms- Methods and Applications

coefficient *cl*

**GA parameters Description**

Crossover rate 0.7 Mutation rate 0.1

Generation size 8 individuals per generation

Number of design variables 28 in total, 14 (7 for *cd* & 7 for *cl*

Selection type Roulette Wheel Selection

*cd* and *cl*

*3.2.1. Application of multi-objective genetic algorithm (MOGA)*

**3.2. Shape optimization of Busemann biplane airfoil under lifting conditions**

Under the lifting condition with the lift coefficient *cl* ≠ 0, a multi-objective genetic algorithm (MOGA) is employed for shape optimization of the original Busemann biplane airfoil. The two objectives to be achieved are to minimize the drag coefficient *cd* while maximizing the lift

seven design points ranging from *M* = 1.1 to *M* = 1.7 are used during the optimization process. For the fitness functions for the lifting case, we use the sum of evenly weighted average of both

. The GA parameters used for the lifting case are listed in **Table 6**.

Error of mutation constant 0.8, which determines how much mutation affects the curves as generations go on

**Table 6.** GA parameters for shape optimization of Busemann biplane airfoil under lifting conditions.

Number of generations Maximum of 50 generations if convergence not obtained

. Similar to that for the non-lifting condition described in Section 3.2, a total of

) for acceleration and 14 (7 for *cd* & 7 for *cl*

) for deceleration

After implementing MOGA for 20 generations with eight individuals in each generation, an optimal shape for asymmetric Busemann-type airfoil under lifting conditions with maximum *cl* and minimum *cd* is obtained. **Figure 22** shows the geometry of the original Busemann biplane airfoil (red) and the optimized Busemann biplane airfoil (blue) under lifting conditions.

**Figure 22.** Geometry of both the original and optimized Busemann biplane airfoil under lifting condition.

The drag coefficients for the seven design points are compared in **Tables 7** and **8** for both the original and optimized Busemann biplane airfoil under lifting condition. As shown in

**Tables 7** and **8**, the baseline Busemann biplane airfoil is choked at all Mach numbers within the optimization range, while the optimized Busemann biplane airfoil is unchoked at *M* = 1.6 during acceleration and is choked at *M* = 1.3 during deceleration. Even under choked conditions during both acceleration and deceleration, the optimized Busemann biplane air‐ foil has significantly lower drag compared to the baseline Busemann biplane airfoil. Similar to the non-lifting condition, the only point where the optimized Busemann biplane airfoil has a higher drag compared to the original airfoil is at *M* = 1.7 during deceleration.

**Table 7.** Comparison of *cd* for the original and optimized Busemann biplane airfoil under lifting conditions during acceleration (1 count = 0.0001).

**Table 8.** Comparison of *cd* for the original and optimized Busemann biplane airfoil under lifting conditions during deceleration (1 count = 0.0001).

**Figures 23** and **24** show a change in the pressure coefficient *cp* around the optimized Busemann biplane airfoil under lifting condition as the Mach number increases and decreases within the design-point range. As can be seen from these figures, the bow shock wave in front of the airfoil disappears at *M* = 1.6 during acceleration and is not generated until the flow speed drops down to *M* = 1.3 during deceleration.

**Figure 23.** *cp* contours for the optimized Busemann biplane airfoil under lifting conditions during acceleration: (a) *cp* contours at *M* = 1.1; (b) *cp* contours at *M* = 1.3; (c) *cp* contours at *M* = 1.5; (d) *cp* contours at *M* = 1.7.

Shape Optimization of Busemann-Type Biplane Airfoil for Drag Reduction Under Non-Lifting and Lifting Conditions Using Genetic Algorithms http://dx.doi.org/10.5772/62811 93

**Figure 24.** *cp* contours for the optimized Busemann biplane airfoil under lifting condition during deceleration: (a) *cp* contours at *M* = 1.7; (b) *cp* contours at *M* = 1.5; (c) *cd* contours at *M* = 1.3; (d) *cd* contours at *M* = 1.1.

**Figure 25.** *cd* plot for different airfoils under lifting conditions.

**Tables 7** and **8**, the baseline Busemann biplane airfoil is choked at all Mach numbers within the optimization range, while the optimized Busemann biplane airfoil is unchoked at *M* = 1.6 during acceleration and is choked at *M* = 1.3 during deceleration. Even under choked conditions during both acceleration and deceleration, the optimized Busemann biplane air‐ foil has significantly lower drag compared to the baseline Busemann biplane airfoil. Similar to the non-lifting condition, the only point where the optimized Busemann biplane airfoil

has a higher drag compared to the original airfoil is at *M* = 1.7 during deceleration.

**Mach number 1.1 1.2 1.3 1.4 1.5 1.6 1.7** Baseline 1097 1040 999 967 943 926 923 Optimized 683 633 529 461 418 148 136

**Table 7.** Comparison of *cd* for the original and optimized Busemann biplane airfoil under lifting conditions during

**Table 8.** Comparison of *cd* for the original and optimized Busemann biplane airfoil under lifting conditions during

**Figures 23** and **24** show a change in the pressure coefficient *cp* around the optimized Busemann biplane airfoil under lifting condition as the Mach number increases and decreases within the design-point range. As can be seen from these figures, the bow shock wave in front of the airfoil disappears at *M* = 1.6 during acceleration and is not generated until the flow speed drops down

**Figure 23.** *cp* contours for the optimized Busemann biplane airfoil under lifting conditions during acceleration: (a) *cp*

contours at *M* = 1.1; (b) *cp* contours at *M* = 1.3; (c) *cp* contours at *M* = 1.5; (d) *cp* contours at *M* = 1.7.

**Mach number 1.7 1.6 1.5 1.4 1.3 1.2 1.1** Baseline 32 926 940<PHONE_NUMBER> 1098 Optimized 136 148 173 224 529 633 682

acceleration (1 count = 0.0001).

92 Optimization Algorithms- Methods and Applications

deceleration (1 count = 0.0001).

to *M* = 1.3 during deceleration.

**Figure 25** shows the comparison of the drag coefficients for the standard diamond-shaped airfoil, the baseline Busemann biplane airfoil and the optimized Busemann biplane airfoil under lifting condition. As shown in the figure, the separation between the acceleration and deceleration *cd* lines still exists for the optimized Busemann biplane airfoil, which means that the flow-hysteresis and the choked-flow effects are not totally eliminated. However, as clearly shown in **Figure 25**, the flow-hysteresis area has been significantly reduced, and the drag increase during deceleration due to the choked-flow phenomenon is much smaller than that for the original Busemann biplane airfoil. The drag of the optimized Busemann biplane airfoil in the subsonic region is smaller than that of the original Busemann biplane airfoil, although it is slightly higher than that of the standard diamond-shaped airfoil with 0.3 < *M* < 0.85. Under

**Figure 26.** *cp* contours of the GA-optimized Busemann biplane airfoil with lift during acceleration: (a) *M* = 1.3, *cd* = 0.05286; (b) *M* = 1.4, *cd* = 0.04609; (c) *M* = 1.5, *cd* = 0.04179; (d) *M* = 1.51, *cd* = 0.04116; (e) *M* = 1.52, *cd* = 0.01665; (f) *M* = 1.6, *cd* = 0.01477.

both subsonic and supersonic conditions, the optimized Busemann biplane airfoil has signif‐ icantly reduced the wave drag compared to the original Busemann biplane airfoil. Similar to the non-lifting condition, the drag coefficient of the optimized Busemann biplane airfoil at *M* = 1.7 (*cd* = 0.01362) is higher than that of the original Busemann biplane airfoil (*cd* = 0.002182).

Next, we examine the details of the flow field during acceleration and deceleration for the optimized Busemann biplane airfoil under lifting conditions and compare them with the flow field of the original Busemann biplane airfoil. **Figures 26** and **27** show the pressure coefficient contours of the optimized Busemann biplane airfoil using GA during acceleration and deceleration, respectively.  The swallowing of the bow shock wave happens when the Mach number increases from 1.49 to 1.50 as shown in **Figure 16(f)** and from 2.12 to 2.13 for the original Busemann biplane airfoil as shown in **Figure 6(d)**. The corresponding drag coefficient decreases from 0.04116 to 0.01665 for the optimized Busemann biplane airfoil.

During deceleration, the choked-flow effect still exists; however, it is shifted to a lower Mach number of 1.38 compared to 1.6 for the original Busemann biplane airfoil.

In conclusion, the drag coefficient of the optimized Busemann biplane airfoil under lifting conditions is significantly reduced compared to the original Busemann biplane airfoil as shown Shape Optimization of Busemann-Type Biplane Airfoil for Drag Reduction Under Non-Lifting and Lifting Conditions Using Genetic Algorithms http://dx.doi.org/10.5772/62811 95

**Figure 27.** *cp* contours of the GA-optimized Busemann biplane airfoil with lift during deceleration: (a) *M* = 1.6, *cd* = 0.01477; (b) *M* = 1.5, *cd* = 0.01733; (c) *M* = 1.4, *cd* = 0.02243; (d) *M* = 1.39, *cd* = 0.02219; (e) *M* = 1.38, *cd* = 0.04716; (f) *M* = 1.3, = 0.05289.

in **Figures 25**–**27**. **Table 9** gives the lift coefficient of the optimized Busemann airfoil under lifting conditions at different design points.

**Table 9.** *cl* for the optimized Busemann biplane airfoil under lifting conditions at different design points.

#### **4. Conclusions**

both subsonic and supersonic conditions, the optimized Busemann biplane airfoil has signif‐ icantly reduced the wave drag compared to the original Busemann biplane airfoil. Similar to the non-lifting condition, the drag coefficient of the optimized Busemann biplane airfoil at *M* = 1.7 (*cd* = 0.01362) is higher than that of the original Busemann biplane airfoil (*cd* = 0.002182).

**Figure 26.** *cp* contours of the GA-optimized Busemann biplane airfoil with lift during acceleration: (a) *M* = 1.3, *cd* = 0.05286; (b) *M* = 1.4, *cd* = 0.04609; (c) *M* = 1.5, *cd* = 0.04179; (d) *M* = 1.51, *cd* = 0.04116; (e) *M* = 1.52, *cd* = 0.01665; (f) *M* = 1.6,

*cd* = 0.01477.

94 Optimization Algorithms- Methods and Applications

Next, we examine the details of the flow field during acceleration and deceleration for the optimized Busemann biplane airfoil under lifting conditions and compare them with the flow field of the original Busemann biplane airfoil. **Figures 26** and **27** show the pressure coefficient contours of the optimized Busemann biplane airfoil using GA during acceleration and deceleration, respectively.  The swallowing of the bow shock wave happens when the Mach number increases from 1.49 to 1.50 as shown in **Figure 16(f)** and from 2.12 to 2.13 for the original Busemann biplane airfoil as shown in **Figure 6(d)**. The corresponding drag coefficient decreases from 0.04116 to 0.01665 for the optimized Busemann biplane airfoil.

During deceleration, the choked-flow effect still exists; however, it is shifted to a lower Mach

In conclusion, the drag coefficient of the optimized Busemann biplane airfoil under lifting conditions is significantly reduced compared to the original Busemann biplane airfoil as shown

number of 1.38 compared to 1.6 for the original Busemann biplane airfoil.

In this chapter, numerical simulations of the flow past the standard diamond-shaped airfoil and the baseline Busemann-type biplane airfoil have been conducted. An impulsive uniform flow, a flow during acceleration, and a flow during deceleration are simulated. The original Busemann biplane airfoil shows a poor performance under off-design conditions due to the flow-hysteresis phenomenon during acceleration and the choked-flow phenomenon during deceleration. For shape optimization, a single-objective genetic algorithm (SOGA) and a multiobjective genetic algorithm (MOGA) are employed to optimize the shape of the Busemanntype biplane airfoil under non-lifting and lifting conditions, respectively, to improve its performance under off-design conditions. The commercially available CFD solver FLUENT is employed to calculate the flow field on an unstructured mesh generated using the ICEM grid generation software. A second-order accurate steady-density-based solver in FLUENT is employed to compute the supersonic flow field. The optimization results for the non-lifting case show a significant improvement in reducing the drag coefficient under off-design conditions for the optimized Busemann-type biplane airfoil shape compared to the original shape. The flow-hysteresis phenomenon during acceleration and the choked-flow phenom‐ enon during deceleration are both alleviated significantly for the optimized shape. For the lifting case, the optimized Busemann biplane airfoil is able to significantly reduce the drag coefficient under off-design conditions while generating lift at the same time.

#### **Author details**

Yi Tian and Ramesh K. Agarwal\*

\*Address all correspondence to<EMAIL_ADDRESS>
Department of Mechanical Engineering and Materials Science, Washington University in St. Louis, St. Louis, MO, USA

#### **References**

Transport Aircraft," AIAA Paper 2006-654, 44rd AIAA Aerospace Sciences Meeting and Exhibit, American Institute of Aeronautics and Astronautics, Reston, VA, 2006.

[8] Maruyama, D., Matsushima, K., Kusunose, K., and Nakahashi, K., "Aerodynamic Design of Biplane Airfoils for Low Wave Drag Supersonic Flight," AIAA Paper 2006-3323, 24th Applied Aerodynamics Conference, 2006; American Institute of Aeronautics and Astronautics, Reston, VA.

objective genetic algorithm (MOGA) are employed to optimize the shape of the Busemanntype biplane airfoil under non-lifting and lifting conditions, respectively, to improve its performance under off-design conditions. The commercially available CFD solver FLUENT is employed to calculate the flow field on an unstructured mesh generated using the ICEM grid generation software. A second-order accurate steady-density-based solver in FLUENT is employed to compute the supersonic flow field. The optimization results for the non-lifting case show a significant improvement in reducing the drag coefficient under off-design conditions for the optimized Busemann-type biplane airfoil shape compared to the original shape. The flow-hysteresis phenomenon during acceleration and the choked-flow phenom‐ enon during deceleration are both alleviated significantly for the optimized shape. For the lifting case, the optimized Busemann biplane airfoil is able to significantly reduce the drag

Department of Mechanical Engineering and Materials Science, Washington University in St.

[1] Busemann, A., *Aerodynamic Lift at Supersonic Speeds*, 12th ed., No. 6, Luftfahrtforschung,

[2] Moeckel, W. E., *Theoretical Aerodynamic Coefficients of the Two-Dimensional Supersonic*

[3] Licher, R., *Optimum Two-Dimensional Multiplanes in Supersonic Flow*, Douglas Aircraft

[4] Tan, H. S., *The Aerodynamics of Supersonic Biplanes of Finite Span*, WADC Tech. Rept.

[5] Ferri, A., *Elements of Aerodynamics of Supersonic Flow*, Macmillan Company, New York,

[6] Igra, D., and Arad, E., "A Parametric Study of the Busemann Biplane Phenomena,"

[7] Kusunose, K., Matsushima, K., Goto, Y., Yamashita, H., Yonezawa, M., Maruyama, D., and Nakano, T., "A Fundamental Study for the Development of Boomless Supersonic

52-276, 1950. Wright Air Development Center. Dayton, OH.

Shock Waves, 16, 3, 2007, 269–273, doi:10.1007/s00193-006-0070-x.

coefficient under off-design conditions while generating lift at the same time.

**Author details**

Yi Tian and Ramesh K. Agarwal\*

96 Optimization Algorithms- Methods and Applications

1935, pp. 210–220.

Louis, St. Louis, MO, USA

**References**

1949.

\*Address all correspondence to<EMAIL_ADDRESS>
*Biplane*, NACA Tech. Rept. 1316, 1947.

Co. Tech Dept. SM-18688, 1955.

**Chapter 5**

### **Performance Analysis of the Differential Evolution and Particle Swarm Optimization Algorithms in Cooperative Wireless Communications**

Arif Basgumus, Mustafa Namdar, Gunes Yilmaz and Ahmet Altuncu

Additional information is available at the end of the chapter

http://dx.doi.org/10.5772/62453

#### **Abstract**

[19] He, Y., "Shape Optimization of Airfoils Without and With Ground Effect Using a Multi-Objective Genetic Algorithm," M.S. Thesis, School of Engineering and Applied Science,

[20] "Some Single- and Multiobjective Optimization Techniques," in Unsupervised Classification – Similarity Measures, Classical and Metaheuristic Approaches, and Applications, Chapter 2, 17–21, Springer-Verlag, Berlin, Heidelberg, Bandyopadhyay,

[21] Available from: https://www.math.ubc.ca/~cass/gfx/bezier.html (accessed 9 Dec 2014).

Washington University in St. Louis, 2014.

98 Optimization Algorithms- Methods and Applications

S., and Saha, S., 2013. doi: 10.1007/978-3-642-32451-2.

In this study, we evaluate the performance of differential evolution (DE) and particle swarm optimization (PSO) algorithms in free-space optical (FSO) and mobile radio communications systems. In particular, we obtain the optimal transmission distances for multiple-relay nodes in FSO communication systems and optimal relay locations in mobile radio communications systems for the cooperative-diversity networks, using both algorithms. We investigate the performance comparison of DE and PSO algorithms for the parallel decode-and-forward (DF) relaying. Then, we analyze the cost functions. Furthermore, we present the execution time and the stability of the DE and PSO algorithms.

**Keywords:** free-space optical communications, cooperative-diversity networks, opti‐ mal distance, differential evolution algorithm, particle swarm optimization algorithm

#### **1. Introduction**

The aim of the optimization was to provide the best-suited solution to a problem under the given constraints. The optimization algorithms have recently been much attention and gained significant importance in plenty of engineering problems [1–8]. In this study, we evaluate the performance analysis of differential evolution (DE) and particle swarm optimization (PSO) algorithms both in free-space optical (FSO) and in mobile radio communications systems.

© 2016 The Author(s). Licensee InTech. This chapter is distributed under the terms of the Creative Commons Attribution License (http://creativecommons.org/licenses/by/3.0), which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.

FSO communications have been proposed as a solution for the various applications including fiber backup, and backhaul for wireless communications networks [9]. Despite the fact that the usage of FSO communications is widespread in major applications for wireless commu‐ nications, the performance limitations for long-range links due to the atmospheric turbulenceinduced fading have had profound impacts in FSO communications systems. The method used for relay-assisted FSO transmission links is one of the fading mitigation technique has attracted significant attentions recently in FSO communications networks [9–11]. In [9], the authors consider relay-assisted FSO communications and investigate the outage performance under serial and parallel scheme with amplify-and-forward (AF), and decode-and-forward (DF) relaying models. The authors in [10] consider a cooperative FSO communications via optical AF relay and investigate the bit error probability performance. Bit error rate (BER) analysis of cooperative systems in FSO networks is presented in [11]. The outage performance analysis of FSO communications is presented in both [12] and [13]. Kashani et al. [14] consider the diversity gain analysis and determine the optimal relay locations for both the serial and the parallel relaying schemes. Although cooperative transmissions have greatly been considered in the above manuscripts, to the best of the authors' knowledge, there has not been any notable research for the relay-assisted FSO communications systems using the optimization algo‐ rithms. To fill the research gap, in this paper, we analyze the performance of both DE and PSO algorithms in terms of the transmission distances when applied for the parallel DF relaying in FSO systems. Moreover, we investigate the performance comparison of these two algorithms in respect to the execution time, cost, and stability analysis.

In our study, as a second part of this paper, we focus on dual-hop cooperative-diversity network to study the impact of the relay location between the source and the destination. Cooperative-diversity relay networks provide a significant performance increment in the radio frequency power transmission and the spatial diversity. They are also shown to be a promising solution to mitigate the signal fading arising from the multipath propagation in wireless communications [15, 16]. In the cooperative-diversity networks, relay terminals are employed between the transmitter and receiver nodes, over multiple communications routes, in which two main protocols are used as follows: (i) amplify-and-forward and (ii) DF [15–21].

Most of the previous publications have studied the cooperative-diversity performance over different fading channels [15–26]. In [15], the authors analyzed the cooperative-diversity network using AF cooperation protocol, operating over independent, but not necessarily identically distributed Nakagami-m fading channels. The paper in [17] addressed the multibranch adaptive DF scheme for cooperative-diversity networks. The best relay selection scheme for cooperative-diversity network is studied in [18]. Furthermore, [19] investigated the advantage of the diversity over direct transmission and conventional non-cooperative relaying scheme. In all these papers, analytical framework for performance analysis of BER and the outage probability is provided [15, 17–19]. As far as we know, both DE and PSO algorithms have not been applied for obtaining the optimal location of the relaying terminal over Nakagami-m fading channel.