Optimization of plate-fin heat exchangers by an improved harmony search algorithm

Applied Thermal Engineering 50 (2013) 877e885

Contents lists available at SciVerse ScienceDirect

Applied Thermal Engineering journal homepage: www.elsevier.com/locate/apthermeng

Optimization of plate-fin heat exchangers by an improved harmony search algorithm Moslem Yousefi a, *, Rasul Enayatifar b, Amer Nordin Darus a, Abdul Hanan Abdullah b a b

Department of Thermo-fluids, Faculty of Mechanical Engineering, Universiti Teknologi Malaysia, 81310 Skudai, Johor, Malaysia Faculty of Computer Science, Universiti Teknologi Malaysia, 81310 Skudai, Johor, Malaysia

h i g h l i g h t s < This is the first application of HS algorithm for plate-fin heat exchanger design. < The original HS algorithm has been modified. < A self-adaptive penalty function is used which eases the constraint handling. < The results show the superiority of this method over GA, PSO and their hybrid.

a r t i c l e i n f o

a b s t r a c t

Article history: Received 18 February 2012 Accepted 30 May 2012 Available online 15 June 2012

This study explores the use of a proposed variant of harmony search algorithm for design optimization of plate-fin heat exchangers. The algorithm deals with a large number of continuous and discrete variables. To handle the constraints in the optimization problem, a self-adaptive penalty function scheme is used. The efficiency and accuracy of the proposed method are demonstrated through an illustrative example taken from previous studies. Numerical results indicate that the presented approach can generate optimum solutions with higher accuracy when compared to Genetic algorithms (GAs), Particle Swarm Optimization (PSO) and GA hybrids with PSO (GAHPSO).
2012 Elsevier Ltd. All rights reserved.

Keywords: Plate-fin heat exchanger Evolutionary computation Harmony search algorithm Constraint handling

1. Introduction Plate-fin heat exchangers (PFHEs), categorized as compact heat exchangers, are commonly used in various industrial areas such as chemical and petrochemical processes, cryogenics, and the automobile and aerospace industries. A drawing of a typical PFHE is presented in Fig. 1. The design of PFHEs is a multifaceted task including thermodynamic and fluid dynamic design, cost estimation, strength calculations and constraint considerations. Moreover, ahead of PFHE design, optimization should take account of the problem requirements [1]. The most common optimization objectives in PFHE design are associated with minimizing capital costs and operating costs. Generally, as the two mentioned objectives are conflicting in nature, a compromise should be achieved for a practical design. In addition to using traditional mathematical methods [2e5], simulated annealing [6] and artificial neural network [7], many * Corresponding author. Tel.: þ60 177290476. E-mail address: yousefi[email protected] (M. Yousefi). 1359-4311/$ e see front matter
2012 Elsevier Ltd. All rights reserved. http://dx.doi.org/10.1016/j.applthermaleng.2012.05.038

researchers have successfully employed evolutionary computation in design optimization of PFHEs [8e17]. These studies have considered various objectives such as minimization of costs, minimization of the number of entropy generation units, minimum weight, minimum total pressure loss and minimum heat transfer area. The results have demonstrated the robustness and applicability of these algorithms in PFHE design. However, due to the continuous improvements in evolutionary computation methods, further studies of the application of these newly introduced methods in PFHE design are needed. Inspired by the improvisation process of musicians, harmony search (HS) algorithm was developed by Geem et al. [18]. It has gained much attention in recent years because (a) it is easier to implement compared to other meta-heuristics and (b) it is adaptable to a wide range of applications given that it can deal with both continuous and discrete variables without additional effort. Recently, the HS algorithm has been effectively used in a wide range of engineering applications. An overview if these applications are presented in Ref. [19]. Moreover, this algorithm has been successfully used in design optimization of thermal systems. A variant of HS algorithm was employed by Fesanghary et al. [20] for

878

M. Yousefi et al. / Applied Thermal Engineering 50 (2013) 877e885

Nomenclature A, AHT Aff bw C Cr Dh f f(X) g(X) G h (W/m2 H HMCR HMS j l L m n N Na, Nb NI NTU

heat exchanger surface area (m2) free flow area (m2) pitch adjusting bandwidth heat capacity rate (W/K) Cmin/Cmax hydraulic diameter (m) friction factor objective function constraint function mass flow velocity (kg/m2 s) K) convective heat transfer coefficient height of fin (m) harmony memory considering rate harmony memory size Colburn factor lance length of the fin (m) heat exchanger length (m) mass flow rate (kg/s) fin frequency (fins per meter) number of decision variables number of fin layers for fluid a and b number of iterations number of transfer units

optimization of shell-and-tube heat exchangers. The presented results have proved the efficiency of this algorithm comparing to the conventional GAs. Therefore, in this study an improved HS approach is presented for solving a PFHE optimization problem. To improve the efficiency of the HS, an additional modification has been applied to the existing improved harmony search [21]. In our proposed algorithm, a Roulette-wheel selection is employed in harmony consideration for increasing the possibility of selection of more powerful individuals. Using this approach, the IHS will be equipped with the advantage of using the global-best individual in the evolution process. In design optimization of a PFHE, various equality and/or inequality constraints such as pressure drops, heat transfer rate, heat exchanger size, etc, have to be considered. However, constraint handling has not received due attention in the previous studies. The evolutionary algorithms, such as genetic algorithms (GAs), particle swarm optimization (PSO), and HS algorithm, cannot deal with the

PAR Pr Q Re R t T U x X

pitch adjusting rate Prandtl number rate of heat transfer (W) Reynolds number penalty parameter fin thickness (m) temperature  C overall heat transfer coefficient a decision variable a set of each decision variables

Greek symbols 3 effectiveness DP pressure drop (N/m2) m viscosity (N/m2 s) y() constraint violation f ðÞ penalty function r density (kg/m3) Subscripts a, b fluid a and b L lower limit max maximum min minimum U upper limit

constraints directly and constraint handling methods should be employed to help the optimization process. The static penalty function method, which is very simple and is the most popular constraint handling method, has been employed in all the previous related studies. However, this method has the drawback of depending on the proper selection of its parameters [22e24]. To make things worse, the penalty parameter is problem-dependant, and has to be defined by a trial-and-error process that could be time consuming and error-prone. Therefore, in this work, a self-adaptive penalty function scheme [25] is adopted for constraint handling. This scheme not only eases the constraint handling process, as there is no parameter setting involved, but also results in better solutions when compared to the static penalty function method. The rest of the paper is organized as follows. Section 2 presents thermal modeling of the PFHEs. Section 3 presents the HS algorithm, its improved variants, and our proposed Roulette-wheel selection. Section 4 sheds light on constraint handling. Finally, the experimental results and conclusions are presented in Section 5 and Section 6 respectively. 2. Thermal modeling

Fig. 1. A schematic drawing of a typical PFHE.

In thermal modeling, the physical property of fluid is considered to be constant for the sake of simplicity. Other assumptions are: (a) the number of fin layers for the cold side (Nb) is one more than for the hot side (Na); (b) the working condition is steady state; (c) the thermal resistance of the walls is neglected; (d) the heat exchanger works under steady state conditions; (e) the heat transfer coefficient and the area distribution are uniform and constant; and (f) fouling is ignored as its influence is negligibly small for a gas-to-gas heat exchanger. In this study, rectangular offset-strip fins are considered because of their compactness, high heat transfer efficiency and reliability. Fig. 2 gives a schematic view of a typical rectangular offset-strip fin. The 3eNTU method is employed for rating the performance of the heat exchanger in the optimization process. The effectiveness of

M. Yousefi et al. / Applied Thermal Engineering 50 (2013) 877e885

879

h 0:1499 j ¼ 0:6522ðReÞ0:5403 ðaÞ0:1541 ðdÞ ðgÞ0:0678  1 þ 5:269 i0:1 0:456 ðgÞ1:055  105 ðReÞ1:34 ðaÞ0:504 ðdÞ (13) f ¼ 9:6243ðReÞ0:7422 ðaÞ0:1856 ðdÞ

0:3053

 108 ðReÞ4:429 ðaÞ0:920 ðdÞ

3:767

h ðgÞ0:2659  1 þ 7:669 i0:1

ðgÞ0:236

(14)

Fig. 2. A schematic drawing of an offset-strip fin.

where a cross-flow heat exchanger, where both fluids are unmixed, is given by Ref. [26]: 3


  n h i o 1 NTU0:22 exp Cr:NTU0:78  1 ¼ 1  exp Cr

(1)

where Cr ¼ Cmin/Cmax and the NTU is calculated as follows: having neglected the thermal resistance of the walls and the fouling factor.

1 1 1 ¼ þ ðhAÞa ðhAÞb UA

(2)

UA Cmin

(3)

NTU ¼

The heat transfer coefficient is calculated from the Colburn factor j: 23

h ¼ j$G$Cp$Pr

(4)

Here G ¼ m/Aff, where Aff is free-flow cross-sectional area that is calculated using the geometrical parameters.

Aff a ¼ ðHa  ta Þð1  na ta ÞLb Na

(5)

Aff b ¼ ðHb  tb Þð1  nb tb ÞLa Nb

(6)

The heat transfer area for both sides can be calculated similarly.

Aa ¼ La Lb Na ½1 þ 2na ðHa  ta Þ

(7)

Ab ¼ La Lb Na ½1 þ 2nb ðHb  tb Þ

(8)

Then, the total heat transfer area is given by:

AHT ¼ Aa þ Ab

(9)

The heat transfer rate is calculated as follows.

  Q ¼ 3 Cmin Ta;1  Tb;1

(10)

The frictional pressure drop in each side is given by:

DPa ¼

DPb ¼

2fa La G2a

ra Dh;a 2fb Lb G2b

rb Dh;b

s h

a ¼ ;d ¼

t t ; g ¼ ; s ¼ ð1=n  tÞ and h ¼ H  t: lf s

The hydraulic diameter for calculating the Reynolds number is defined as below.

Dh ¼

4s$h$l 2ðs$l þ h$l þ t$hÞ þ t$s

Eqs. (13) and (14) are valid for 120 < Re < 104, 0.134 < a < 0.997, 0.012 < d < 0.048 and 0.041 < g < 0.121. These equations correlate the j and the f factors from the experimental data with 20% accuracy in the laminar, the transition and the turbulence flow regimes. Therefore, there is no need to describe the flow regime for a specific operating condition, which is very useful in most practical applications. 3. Optimization method 3.1. Harmony search In the basic Harmony search algorithm, an n-dimensional real vector called a harmony represents each solution. An initial population of harmony vectors is generated randomly and allocated in harmony memory (HM). The HM is basically similar to the genetic pool in GA. Subsequently, a new harmony is generated based on all of the solutions in the HM using a harmony memory considering rate (HMCR) and a pitch adjustment rule. HMCR is basically the possibility of choosing a value for the new harmony from only the existing values in the HM. Each value, selected by this action, should undergo a pitch adjustment process with a probability of pitch adjustment rate (PAR). Each iteration is completed after updating the HM by comparing it to the newly generated harmony vector. The new vector replaces the worst harmony of the HM if it has a better state. The algorithm is repeated until a certain termination criterion is met (or simply by reaching the allowed maximum number of iterations). The three basic steps of the algorithm are described below in detail. Step 1 e Initialize the problem and algorithm parameters. The optimization problem is described as follows.

Minimise f ðXÞ; X ¼ ½x1 ; .xN  (11)

There are many correlations for the evaluation of the Colburn factor j and the fanning factor f for offset-strip fins. Eqs. (13) and (14) are the correlations presented by Manglik and Bergles [27] and are used in this work.

(16)

where constraints are stated as

gj ðXÞ  0; j ¼ 1; .; m (12)

(15)

(17)

and

xLi  xi  xU i

i ¼ 1; .; N

(18)

Here f(X) is the desired objective function; X is the set of each decision variable; N is the number of decision variables and xi is the L U range of decision variables that are xLi  xU i where xi and xi are the

880

M. Yousefi et al. / Applied Thermal Engineering 50 (2013) 877e885



Table 1 The selected parameters for the HS variants. HMS HMCR PAR

HS IHS GHS The proposed method

5 5 5 5

0.9 0.9 0.9 0.9

bw

PAR PARmin PARmax bw

bwmin bwmax

0.3 N.A. N.A. N.A.

N.A. 1E-04 N.A. 1E-04

N.A. 0.1 0.1 0.1

N.A. 0.99 0.99 0.99

0.01 N.A. N.A. N.A.

x11 x21 «

x12 x22 «

6 6 6 HM ¼ 6 6 HMS1 HMS1 x2 4 x1 xHMS xHMS 1 2

. . « . .

x1N1 x2N1 «

N.A. (UB-LB)/20 N.A. (UB-LB)/20

xHMS1 N1 xHMS N1

3

x1N x2N «

with probability PAR with probabilityð1  PARÞ:

(20)

The pitch adjustment rule is defined as follows.

lower and upper limits. The HS parameters: harmony memory size (HMS), harmony memory considering rate (HMCR), pitch adjusting rate (PAR), bandwidth (bw) and the termination criterion or simply the number of improvisations (NI) are set in this step. The typical values of these parameters are presented in Table 1. Step 2 e Initialize the harmony memory. A HMS number of randomly generated solution vectors are stored in the HM.

2

Yes No

x0i )x0i Hrand



 bw

(21)

where bw is an arbitrary bandwidth that is set to different values for each decision variable and rand () is a random number between 0 and 1. Step 4 e Update the harmony memory. In this step, if the new harmony is better than the worse solution vector of the HM, it is included in the HM, replacing the worse one. Step 5 e Checking termination criterion. The algorithm continues until the termination criterion, that could be simply the number of improvisations (NI) is met.

3.2. The proposed harmony search algorithm Different methods have been proposed for improving the efficiency of the HS algorithm The improved harmony search algorithm (IHS) was presented by Mahdavi et al. [21] and it follows the same steps as the HS algorithm, with the only difference that the values of the PAR and bw are updated in each iteration dynamically as follows.

7 7 7 7 7 HMS1 xN 5 HMS xN

Step 3 e Improvise a new harmony. A new harmony vector, X 0 ¼ ðx01 ; x02 ; .:; x0n Þ; is generated in this step based on the HS rules. In the memory consideration, the value of each decision variable is chosen randomly from all corresponding existing values in the HM. The HMCR that can vary between 0 and 1 is the probability of choosing one value from the existing values stored in the HM, while (1  HMCR) is the rate of randomly generating one value from the possible ranges of values.

PARmax  PARmin t NI

PARðtÞ ¼ PARmin þ

(22)

 0  1 bwmin In B C bwmax B  tC @ A NI bwðtÞ ¼ bwmax e





x0 ˛ x1i ; x2i ; .; xHMS i X0) i 0 xi ˛xi

with probability HMCR with probabilityð1  HMCRÞ (19)

As an example, a typical HMCR of 85% means that the HS algorithm would choose the decision variable from the stored variables in HM with a probability of 85% or from the entire range of variable with a probability of (1e85%). Every decision variable attained from memory consideration should undergo pitch adjustment with the probability of PAR as follows.

(23)

In Eqs. (22) and (23) PAR (t) and bw (t) represent the pitch adjustment rate and bandwidth in iteration t respectively. PARmax is the maximum pitch adjustment rate and PARmin is its minimum. bwmax and bwmin are maximum and minimum limits for bandwidth respectively. In another study, motivated by the concept of swarm intelligence, Omran and Mahdavi [28] presented global-best harmony search algorithm (GHS). They modified the pitch adjustment rule so that, unlike the basic HS, the presented algorithm generates a new harmony by directly using the best harmony vector in the HM. The

Table 2 Mean and standard deviation (SD) of the benchmark function optimization results. HS Sphere Rosenbrock Step Rastrigin Ackley

Mean Mean Mean Mean Mean

    

SD SD SD SD SD

0.071377 534.883721 10.054982 11.530824 16.048601

IHS     

0.169067 643.935556 4.526308 3.719364 1.418158

0.099251 504.485368 7.452871 14.840561 17.239631

GHS     

0.222471 252.597506 3.914539 2.901781 0.905531

The proposed HS

0.001031 65.166415 9.89E-04 0.019196 7.970873

    

0.001213 64.582174 9.11E-04 0.055387 7.450370

1.56245E-5 12.493823 8.02956E-6 7.03975E-7 3.82701E-06

    

1.8892E-5 1.768479 4.3204E-6 1.1704E-7 3.2158E-06

Table 3 Mean and standard deviation (SD) with varying HMS (n ¼ 30). HMS/

5

Sphere Rosenbrock Step Rastrigin Ackley

1.56245E-5 12.493823 8.02956E-06 7.03975E-7 3.82701E-06

10     

1.8892E-5 1.768479 4.3204E-06 1.1704E-7 3.2158E-06

5.4275E-03 26.32379303 7.92375E-04 8.6007E-07 5.3490321E-07

20     

8.4952E-03 5.84540037 9.0569E-04 2.4638E-07 2.80034E-07

9.29634E-02 175.800543 6.54657E-03 5.96305E-09 4.834015034

50     

3.94752E-02 14.874671 1.8344E-03 5.1734E-09 0.50274821

7.54645E-02 367.3654873 8.92172E-05 4.728307E-03 7.0987636203

    

9.5471E-02 1.7010405 6.3204E-05 1.63132E-03 1.324762741

M. Yousefi et al. / Applied Thermal Engineering 50 (2013) 877e885

881

6.30194E-3 150.3582218 349.01293 35.6275037 76.2820536

0.99

1.8892E-5 4.3204E-6 1.768479 1.1704E-7 3.2158E-06      1.56245E-5 8.02956E-6 12.493823 7.03975E-7 3.82701E-06

0.9

1.0048E-03 1.987343 35.30185 1.1704E-7 0.141525      0.8

1.3148E-03 4.6904871 399.01293 3.03975E-7 1.3489617 0.09042 2.8351254 94.700832 0.9594321 0.484023      0.0872831 7.83423295 2.8074Eþ03 8.50639237 2.4883206

0.7

5.8223E-3 1.4886214 4.382Eþ05 3.987234 1.984323      9.30194E-3 3.98239027 7.4934Eþ05 22.9233927 16.802136

0.6

0.103931 3.2162238 94.70832 6.2903843 3.9702323 1.3940328 12.3582218 2.80074Eþ3 25. 6275037 78.2820536

     0.5 HMCR/

where xnew and xB are the new harmony and the best harmony vectors respectively. k is an arbitrary integer between 1 and n. Moreover, the same procedure as IHS is employed for dynamically updating PAR values. This method, despite its reported results and great name, has serious deficiencies. It may break the building structures of the solutions, especially when the global optimum of the problem at hand has different numerical values in its various dimensions, and generate infeasible new harmonies. Therefore, GHS has been based on an inaccurate foundation and its reported results do not have due credibility. In this study, inspired by the survival of the fittest concept, a mechanism similar to Roulette-wheel selection is employed in the harmony considering process. In this newly proposed method, instead of randomly selecting among all corresponding harmony memory variables with the same possibility, the better solutions are given a higher possibility of participating in the newly improvised (generated) solution vector. The Roulette-wheel selection points out that the probability of each solution is proportionate to its relative fitness value. This probability is calculated by simply dividing the fitness value of each solution to the sum of all fitness values and subtracting it from 1. This selection method makes use of the “globalbest” concept without interfering in the building structure of the solutions. The efficiency of the proposed algorithm is tested on several well-studied benchmark problems, namely Sphere, Rosenbrock, Step, Rastrigin and Ackley. The description of these functions and their optimum results can be seen in Ref. [21,28]. The attained results are compared with those obtained by HS algorithm, GHS and IHS. The parameters used in different algorithms are presented in Table 1. For a consistent comparison, the parameters of the proposed method are considered similar to those of the existing IHS algorithm [21] and the termination criterion is set to reaching 50,000 iterations. Table 2 presents the results of all the test functions obtained by the HS algorithm its variants and the proposed method. Each algorithm is run 30 times on each problem and the presented results are the average and standard deviations of these runs. It can be seen that the proposed method generates the best results for all test functions and outperforms the HS algorithm and its variants namely, IHS and GHS. In this subsection, to determine the optimum harmony memory size, the proposed algorithm is tested on all test functions considering different values of this parameter, i.e., HMS ¼ 5, 10, 20 and 30. The best results are presented in Table 3 with a bold font. Although for Rastrigin and Ackley a larger HMS results in a better solution, generally the HMS ¼ 5 results in a better performance of the proposed algorithm. Moreover, the attained results for Rastrigin and Ackley with an HMS of 5 are very competitive with those obtained by HMS ¼ 20 and HMS ¼ 10 respectively. Therefore, to avoid additional computational effort, the HMS ¼ 5 is considered as the optimum harmony size. The value of HMCR parameter can vary between 0 and 1 and its typical value in previous studies has been HMCR ¼ 0.9. Table 4 presents the effect of HMCR value on the results of the proposed algorithm. The average and deviation of results in 30 runs of the computer code are presented. Consistent with the previous studies [18,21,28], It is observed that HMCR ¼ 0.9 is a better configuration of the algorithm and the corresponding results, which are presented by a bold font in Table 4, are more accurate in all test functions. In our study for PFHE design, a typical HM generated at the beginning of the optimization process is as follows.

    

(24)

Sphere Step Rosenbrock Rastrigin Ackley

j ¼ 1; 2; .; n

Table 4 Mean and standard deviation (SD) with varying HCMR.

Xnew ðjÞ ¼ XB ðkÞ;

5.8223E-3 3.2162238 35.30185 6.2903843 3.9702323

presented pitch adjustment rule, which eliminates the difficulties of selecting bw values, is described as follows.

882

M. Yousefi et al. / Applied Thermal Engineering 50 (2013) 877e885 Table 5 Pseudo code for the proposed method. Input: Plate-fin heat exchanger variables Output: Minimum objective function 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39

2

0:9568 6 0:8634 6 HM ¼ 6 6 0:6697 4 0:4280 0:4318

STEP 1: Initialize the general parameters of Harmony search HMS ) 0 WHILE number of population has not reached DO Increment HMS j)1 WHILE values of all variables have not assigned DO HM (HMS, j) ) Randomly assign the value of each variable in its range Increment j END WHILE Fitness (HMS) ) Evaluate HM (HMS, 1: j) considering desired objective Adding penalty according to Scheme 2 END WHILE STEP 2: Sort all population according to their corresponding Fitnesses Pop ) Select one of the population with the worst Fitness STEP 3: WHILE stopping condition is not met DO k)1 WHILE all variables in Pop are not considered DO IF Rand (0 1)  HMCR THEN PopðkÞ)Using Roulette  Wheel for determining PopðkÞ˛fVariable1k ; Variable2k ; .; VariableHMS g k ELSE Pop (k) ) Determine Pop (k) randomly in its range END IF Update PAR, BW according to Eqs. (22) and (23) respectively IF Rand (0,1)  PAR THEN PopðkÞ)PopðkÞ  randð0; 1Þ  BW ELSE Pop (k) will not change END IF Increment k END WHILE Fit ) Evaluate Pop corresponding to their objectives IF Fit  Fitness (HMS) THEN HM (HMS,:) ¼ Pop END IF Increment loop counter END WHILE

0:88 0:99 0:82 0:11 0:10

0:0021 0:0041 0:0024 0:0041 0:0032

0:00013 0:0019 0:0017 0:00020 0:00010

3 124 197 7 7 58 7 7 102 5 101

For pitch adjusting, discrete values (number of hot side layers) are initially treated as continuous ones, then the achieved values are rounded up. The pseudo code of the proposed method is presented in Table 5.

4. Constraint handling In the present optimization problem, a modified IHS is used for constraint minimization. The problem can generally be described as Eqs. (16)e(18). To handle the constraints in the optimization algorithm, a penalty function method is employed. This method converts the unconstrained problem to a constrained one by adding an additional value corresponding to the level of constraint violation to the original objective value.

Minimise f ðXÞ þ

m  X

f gj ðXÞ



f is a predefined function. The static penalty function, defined in Scheme 1, is the most common penalty function scheme for constrained problems. This scheme is the only method that has been employed in PFHE design so far. Throughout this paper, this static penalty function is referred as Scheme 1.where R is a penalty parameter that has a relatively large value comparing to f(x) and yj(x) is the constraint violation for the jth constraint which is the difference between the actual value of any constraint and its specified limit. The difficulty of using this scheme is determining how to decide the proper amount of R for controlling the severity of the penalty associated with the constraint violation. Excessively large amounts of penalties eliminate the infeasible solutions quickly, but decrease search efficiency, while insufficient penalties usually result in infeasible final solutions. The penalty parameter, R, has to be set by a trial-and-error process. This process should be repeated for any desired objective function because R is problemdependant. To avoid a time-consuming, error-prone and challenging process of selecting penalty parameters, an adaptive penalty scheme presented by Barbosa and Lemonge [25] is employed in this study, and is referred as Scheme 2 in this work. This scheme makes use of the information from all the solutions

(25)

j¼1

Subject to xLi  xi  xU i

i ¼ 1; .; N:

Scheme 1.

M. Yousefi et al. / Applied Thermal Engineering 50 (2013) 877e885

883

Table 7 Variation ranges of design parameters [8]. Parameter

Min

Max

Hot flow length (La) (m) Cold flow length (Lb) (m) Fin height (H) (mm) Fin thickness (t) (mm) Fin frequency (n) (m1) Fin offset length (lf) (mm) Number of hot side layers (Na)

0.1 0.1 2 0.1 100 1 1

1 1 10 0.2 1000 10 200

Scheme 2.

vectors to define different penalties for different constraints. f is described as follows. Where hf ðXÞi is the average of the objective function values of all the solution vectors in the current iteration. For calculating this value, the violation penalties are not considered and the actual fitness value is used. hvl ðXÞi is the violation of the lth constraint averaged over all the solutions existing in the current iteration. It means that for all solutions in the HM, which are 5 in our study, only the violation of the lth constraint is considered and averaged. Scheme 2 is employed for constraint handling in this work unless mentioned otherwise.

For case study A, the total heat transfer area, AHT, is calculated using Eq. (9). In case study B, the total pressure drop is calculated as follows.

f ðxÞ ¼

DPa DPb þ Pa;max Pb;max

(29)

The pressure drop and heat capacity constraints should be considered in the optimization process. Moreover, additional constraints are implemented. In a nutshell, the optimization problem at hand is a large-scale, combinatorial problem that deals with both continuous and discrete variables. 5.3. Algorithm parameters

5. Computational results In this study, PARmin and PARmax are considered to be 0.1 and 0.99 respectively. bwmin and bwmax are defined as follows.

5.1. Illustrative example An illustrative example from the author’s previous work [8] is considered to demonstrate the application of the adaptive penalty function scheme and proposed IHS. A gas-to-air single pass crossflow heat exchanger having a heat duty of 1069.8 kW needs to be designed to achieve minimum heat transfer area and minimum total pressure drop separately. The maximum flow length and noflow length, Lc, of the heat exchanger are 1 m and 1.5 m respectively. The operating conditions are listed in Table 6. In this study, a total number of seven parameters, namely the hot flow length (La), the cold flow length (Lb), the number of hot side layers (Na), the fin frequency (n), the fin thickness (t), the fin height (H), and the fin strip length (lf) are considered as the optimization variables. All variables except the number of hot side layers are continuous. The thickness of the plate, tp, is considered to be constant at 0.5 mm and is not to be optimized. The variation ranges of the variables are shown in Table 7.

5.2. Objective function In this work, the optimization targets two single-objective functions. The first is minimization of the heat transfer area, which is mainly associated with the capital cost of the heat exchanger, and the other is minimization of the total pressure drop, which represents the operating cost.

Table 6 Operating parameters selected for the case study [8]. Parameters

Hot side

Cold side

Mass flow rate, m (kg/s) Inlet temperature, T ( C) Specific heat, Cp (J/kg K) Density, r (kg/m3) Dynamic viscosity, m (N s/m2) Prandtl number, Pr Maximum pressure drop, DP (N/m2)

1.66 900 1122 0.6296 401E-7 0.731 9.50

2 200 1073 0.9638 336E-7 0.694 8.00

  bwðiÞmax ¼ 0:5  xU  xL i U i L  bwðiÞmin ¼ 0:01  xi  xi

(30)

i ¼ 1,.7 The HMCR and HMS are set to 0.9 and 5 respectively, while the number of iterations (NI) is 5000. 5.4. Minimum heat transfer area (case study A) Considering the pressure drop restrictions and required heat duty, the optimization problem is to find the PFHE geometrical design parameters that minimize the heat transfer area. Fig. 3 presents the evolution process of the IHS method. A significant decrease in the target function is observed in the early stages (the first 10 iterations). After a certain number of iterations (more than 50), the changes in the fitness function become relatively small. Finally, the minimum heat transfer area of the PFHE is found to be 109.62 m2 after 4237 iterations. Table 8 shows the features of the PFHE preliminary design and the optimum configuration found by the proposed algorithm. 5.5. Minimum total pressure drop (case study B) Fig. 4 presents the iteration process of the proposed algorithm to find the minimum total pressure drop. A considerable decline in the target function is observed in the early stage (the first 100 iterations). Then, the changes in the objective function become minute until the final result is achieved after 3461 iterations. The minimum total pressure drop is found to be 0.054. The optimized configuration of the PFHE is presented in Table 9. 5.6. A comparison between Schemes 1 and 2 for constraint handling The optimization results using Scheme 1 for constraint handling are presented denoting * in Tables 8 and 9. It can be observed that Scheme 2 can provide better results in terms of accuracy compared

884

M. Yousefi et al. / Applied Thermal Engineering 50 (2013) 877e885 Table 9 Optimized results for minimum pressure drop.

Design variables

Constrained condition Objective

Fig. 3. Minimizing heat transfer area by IHS.

Side a length (La) (m) Side b length (Lb) (m) Fin height (H) (mm) Fin thickness (t) (mm) Fin frequency (n) (m1) Lance length (lfa) (mm) Number of hot side layers (Na) Constrained DPa at hot side (kPa) Condition DPb at clod side (kPa) No-flow length, LC (m) Objective Heat transfer area (m2) Design variables

Hybrid GA [8]

Present approach

Side a length (La) (m) Side b length (Lb) (m) Fin height (H) (mm) Fin thickness (t) (mm) Fin frequency (n) (m1) Lance length (lfa) (mm) Number of hot side layers (Na) DPa at hot side (kPa) DPb at clod side (kPa) No-flow length, LC (m) Total pressure drop

1.00 1.00 10 0.1 241 10 71 0.29 0.21 1.5 0.056

1.00 1.00 10 0.1 211 10 71 0.28 0.19 1.5 0.054 (0.056*)

Table 10 The effect of the penalty function parameter on the results in Scheme 1.

Table 8 Preliminary design and IHS results for minimum heat transfer area. Variables

Variables

Preliminary Hybrid GA [8] Present design [8] approach 0.3 0.3 2.49 0.10 782 3.2 167

0.21 0.23 5.9 0.10 1000 2.1 91

0.20 0.21 6.7 0.10 999.8 1.7 87

9.34 6.90 1 142.75

9.50 8.00 1.18 112.69

9.44 7.96 1.26 109.62 (112.45*)

to Scheme 1. To better understand the difficulty of parameter setting in Scheme 1, the optimization results for case study A using different values of penalty parameter, R, is presented in Table 10. After a trial-and-error process, the best penalty factor value, R, for this problem is found to be 100. In Table 10 any constraint violation or infeasible solution is marked with *. It can be observed that selecting light penalties could result in the infeasible solution having constraint violations, while selecting heavy penalties, R ¼ 200 in this case, degrades the search diversity and hence results in poor solutions. Choosing excessive penalties (R ¼ 1e3 and R ¼ 1e6) degrades constraint handling into a reject strategy where all infeasible solutions are eliminated immediately. These results, presented in Table 11, are the minimum amount achieved in 10 executions of the algorithm. However, it

Penalty parameter (R) DPa DPb Lc Q Objective (AHT)

1

25

50

75

100

200

1000

1e6

9.95* 8.7* 2.4* 1070 89.65*

9.61* 8.00 1.63* 1069.9 104.1*

9.45 8.08* 1.6* 1068.9* 104.84*

9.33 7.92 1.54* 1069.9 107.7*

9.39 7.84 1.15 1069.9 112.50

9.5 7.7 1.13 1069.9 113.23

9.38 7.38 0.9 1069.9 121.2

9.38 7.38 0.9 1069.9 121.2

should be noted that, because the optimization algorithm is naturally random-based, even with a harsh penalty, it is possible to achieve an optimum or near-optimum solution with numerous executions. This means that the success rate (the number of times an optimum solution is achieved in 100 executions) is decreased with an excessive penalty factor. In the case of R ¼ 1e6, an optimum of 113.23 is achieved once in 500 executions of the code. 5.7. A comparison of the proposed algorithm with GA, PSO and GA hybrid with PSO (GAHPSO) In this subsection, the proposed method is compared with GA, PSO and Hybrid GA algorithms in terms of computational time and accuracy. The parameters of these algorithms are presented in Ref. [8]. To present the effects of the adaptive penalty function scheme used in this work, the results obtained by the proposed algorithm using a static penalty function (Scheme 1) is also presented in Table 11. All algorithms are programmed in MATLAB and run on an Intel Core i5 CPU. The mentioned CPU time is the average of 100 executions of the computer code. The aforementioned results indicate that the proposed optimization method, an improved harmony search combined with the adaptive-penalty function scheme, can converge to the optimum results in both cases with higher accuracy in an approximately similar time. The proposed algorithm is superior to GA and PSO in term of accuracy, even when using the same penalty function

Table 11 A comparison of the proposed algorithm with GA, PSO and GA hybrid with PSO (GAHPSO). Parameters

Case Heat transfer study A area (m2) CPU time (s) Case Total study B pressure drop CPU time(s) Fig. 4. Evolution process of the objective of minimum total pressure drop.

GA

PSO

The GAHPSO The proposed IHS proposed IHS Scheme 2 Scheme 1

132.5 123.5 112.7

112.45

109.62

4.40 3.59 3.40 0.070 0.063 0.056

3.17 0.056

3.19 0.054

4.30

3.05

3.07

3.49

3.26

M. Yousefi et al. / Applied Thermal Engineering 50 (2013) 877e885

scheme. In this case, the achieved optimum results are as accurate as the GAHPSO results. 6. Conclusions This study presents the successful application of a proposed variant of harmony search algorithm that works with a selfadaptive penalty function scheme for optimal design of plate-fin heat exchangers. To improve the convergence rate and accuracy of the IHS, inspired by the survival of the fittest concept, a Roulettewheel selection technique was employed in harmony consideration to increase the possibility of selection of more powerful individuals. The numerical results on benchmark problems prove the efficiency of the proposed improvement. HS algorithm is simple in concept, easy to implement, and has few parameters. Moreover, it can deal with both continuous and discrete variables and does not need any information about derivatives. Therefore, harmony search algorithm is a promising algorithm for solving engineering problems, especially in thermal system design where the problems are usually non-convex and have a large amount of discrete variables or discontinuity in the objective function. Furthermore, applying a self-adaptive penalty function scheme for constraint handling eliminates the trial-anderror process of selecting penalty parameters, and results in better solutions when compared to a static penalty function scheme. The simplicity of implementing this scheme and its efficiency make it favorable for any constrained optimization problem. Acknowledgements This work was partially financed by International Doctoral Fellowship (IDF) provided by University Technology Malaysia (UTM) and the Ministry of Higher Education, Malaysia. The authors are also grateful to the anonymous referees for their constructive comments and suggestions that help them elevate the quality of the work to its present level. References [1] W.M. Kays, London L Compact Heat Exchangers, McGraw-Hill, New York, 1984. [2] J.M. Reneaume, H. Pingaud, N. Niclout, Optimization of plate fin heat exchangers: a continuous formulation, Chemical Engineering Research and Design 78 (6) (2000) 849e859. [3] M.R.J. Nasr, G.T. Polley, An algorithm for cost comparison of optimized shelland-tube heat exchangers with tube inserts and plain tubes, Chemical Engineering and Technology 23 (3) (2000) 267e272. [4] K. Muralikrishna, U.V. Shenoy, Heat exchanger design targets for minimum area and cost, Chemical Engineering Research and Design 78 (2) (2000) 161e167. [5] J.M. Reneaume, N. Niclout, MINLP optimization of plate fin heat exchangers, Chemical and Biochemical Engineering Quarterly 17 (2003) 65e76. [6] J.M. Reneaume, N. Niclout, Plate Fin Heat Exchanger Design Using Simulated Annealing, Elsevier, City, 2001.

885

[7] R. Jia, B. Sunden, Optimal design of compact heat exchangers by an artificial neural network method, ASME Conference Proceedings 1 (2003) 655e664. [8] M. Yousefi, R. Enayatifar, A.N. Darus, Optimal design of plate-fin heat exchangers by a hybrid evolutionary algorithm, International Communications in Heat and Mass Transfer 39 (2) (2012) 258e263. [9] M. Yousefi, A.N. Darus, H. Mohammadi, An imperialist competitive algorithm for optimal design of plate-fin heat exchangers, International Journal of Heat and Mass Transfer 55 (11e12) (2012) 3178e3185. [10] M. Yousefi, A.N. Darus, Optimal design of plate-fin heat exchangers by particle swarm optimization, Proceedings of the SPIE (2011) (8350) 83500Re83500R-6. [11] P. Ahmadi, H. Hajabdollahi, I. Dincer, Cost and entropy generation minimization of a cross-flow plate fin heat exchanger using multi-objective genetic algorithm, Journal of Heat Transfer 133 (2) (2011) 021801e021810. [12] S. Sanaye, H. Hajabdollahi, Thermal-economic multi-objective optimization of plate fin heat exchanger using genetic algorithm, Applied Energy 87 (6) (2010) 1893e1902. [13] R.V. Rao, V.K. Patel, Thermodynamic optimization of cross flow plate-fin heat exchanger using a particle swarm optimization algorithm, International Journal of Thermal Sciences 49 (9) (2010) 1712e1721. [14] H. Peng, X. Ling, E. Wu, An improved particle swarm algorithm for optimal design of plate-fin heat exchangers, Industrial and Engineering Chemistry Research 49 (13) (2010) 6144e6149. [15] M. Mishra, P.K. Das, S. Sarangi, Second law based optimisation of crossflow plate-fin heat exchanger design using genetic algorithm, Applied Thermal Engineering 29 (14e15) (2009) 2983e2989. [16] G.N. Xie, B. Sunden, Q.W. Wang, Optimization of compact heat exchangers by a genetic algorithm, Applied Thermal Engineering 28 (8e9) (2008) 895e906. [17] H. Peng, X. Ling, Optimal design approach for the plate-fin heat exchangers using neural networks cooperated with genetic algorithms, Applied Thermal Engineering 28 (5e6) (2008) 642e650. [18] Z.W. Geem, J.H. Kim, G.V. Loganathan, A new heuristic optimization algorithm: harmony search, Simulation 76 (2) (2001) 60e68. [19] G. Ingram, T. Zhang, Overview of Applications and Developments in the Harmony Search Algorithm Music-Inspired Harmony Search Algorithm, Springer, Berlin/Heidelberg, City, 2009. [20] M. Fesanghary, E. Damangir, I. Soleimani, Design optimization of shell and tube heat exchangers using global sensitivity analysis and harmony search algorithm, Applied Thermal Engineering 29 (5e6) (2009) 1026e1031. [21] M. Mahdavi, M. Fesanghary, E. Damangir, An improved harmony search algorithm for solving optimization problems, Applied Mathematics and Computation 188 (2) (2007) 1567e1579. [22] A. Ponsich, C. Azzaro-Pantel, S. Domenech, L. Pibouleau, Constraint handling strategies in genetic algorithms application to optimal batch plant design, Chemical Engineering and Processing: Process Intensification 47 (3) (2008) 420e434. [23] V.S. Summanwar, V.K. Jayaraman, B.D. Kulkarni, H.S. Kusumakar, K. Gupta, J. Rajesh, Solution of constrained optimization problems by multi-objective genetic algorithm, Computers and Chemical Engineering 26 (10) (2002) 1481e1492. [24] C.A. Coello Coello, Theoretical and numerical constraint-handling techniques used with evolutionary algorithms: a survey of the state of the art, Computer Methods in Applied Mechanics and Engineering 191 (11e12) (2002) 1245e1287. [25] H.J.C. Barbosa, A.C.C. Lemonge, A new adaptive penalty scheme for genetic algorithms, Information Sciences 156 (3e4) (2003) 215e251. [26] F.P. Incropera, D.P. Dewitt, T.L. Bergman, A.S. Lavine, Fundamentals of Heat and Mass Transfer, John Wiley & Sons Inc, 2010. [27] R.M. Manglik, A.E. Bergles, Heat transfer and pressure drop correlations for the rectangular offset strip fin compact heat exchanger, Experimental Thermal and Fluid Science 10 (2) (1995) 171e180. [28] M.G.H. Omran, M. Mahdavi, Global-best harmony search, Applied Mathematics and Computation 198 (2) (2008) 643e656.