- Email: [email protected]
Optimal design of plate heat exchangers
Applied Thermal Engineering 63 (2014) 33e39
Contents lists available at ScienceDirect
Applied Thermal Engineering journal homepage: www.elsevier.com/locate/apthermeng
Optimal design of plate heat exchangers Fábio A.S. Mota, Mauro A.S.S. Ravagnani*, E.P. Carvalho Chemical Engineering Graduate Studies Program, State University of Maringá, Av. Colombo, 5790 Maringá, PR, Brazil
h i g h l i g h t s The paper presents an algorithm Q1 for the optimization of heat exchange area of plate heat exchangers (PHE). The optimization algorithm is based on the screening method. Optimal configurations are found for each kind of plate with the lowest number of plates and the smallest area. For all kinds of plates considered, their respective local optima are determined. The global optimum is obtained by comparison among the local optima. The differential equations are solved analytically.
a r t i c l e i n f o
a b s t r a c t
Article history: Received 4 March 2013 Accepted 25 September 2013 Available online 23 October 2013
In the current paper, it is presented an algorithm for the optimization of heat exchange area of plate heat exchangers (PHE). The algorithm is based on the screening method. For each kind of plate, subject to certain constraints, optimal configurations are found e in case they exist e which present the lowest number of plates and, consequently, the smallest area. Each of these found configurations have local optima characteristics. For all kinds of plates considered, their respective local optima are determined. Comparison of all obtained local optima gives a global one. The differential equations, which are generated in the model development, are solved analytically. A case study is presented to test the applicability of the developed algorithm. Results show coherency with the literature. 2013 Elsevier Ltd. All rights reserved.
Keywords: Plate heat exchangers Optimization Simulation Screening
1. Introduction The competitive pressures of the global market and the increasing urgency for both energy conservation and reduction of environmental impacts have changed the focus on the industrial processes, so that heat exchangers with high effectiveness can be used. Plate heat exchanger (PHE) is categorized at the lower end of the compactness spectrum. It offers many advantages and unique characteristics applications when compared with other highly compact heat exchangers. This is due to the flexible thermal design (the plates can be simply added or removed to attend different demands of thermal load and processing), ease to clean maintaining extreme hygiene (necessary when food, pharmaceutical or other kind of products are processed), good temperature control (necessary in cryogenic uses and mitigation of thermal degradation of the process fluids) and better performance in heat transfer [1]. Because of their
* Corresponding author. E-mail addresses: [email protected], [email protected] (M.A.S. S. Ravagnani). 1359-4311/$ e see front matter 2013 Elsevier Ltd. All rights reserved. http://dx.doi.org/10.1016/j.applthermaleng.2013.09.046
characteristics, PHE became ideal for diary, pharmaceutical, food and drink industries. Besides is being used in synthetic rubber industry, pulp and paper industry, maritime applications, solar energy etc. [1,2]. Gut and Pinto [3] proposed the screening method to find the optimal set of configurations of a given plate heat exchanger. In this method, the objective function e which implicitly depends on six configuration parameters e is formulated to obtain the minimum number of channels, which are proportional to the heat transfer area. Constraints are given in order to avoid impossible or non-optimal solutions. The initial set of configurations (IS) is given by the maximum number of possible configurations for the problem and is generated by applying the constraint of channel numbers. The reduced setting (RS) is identified by verifying the pressure drop and flow rate, that is, the IS subset which obeys hydraulic constraints. The optimal set (OS) is formed by the configurations with the minimum number of channels of the RS obeying the thermal constraint. Najafi H. and Najafi B. [4] developed a method to obtain an optimal set of the PHE geometric parameters, with the minimum pressure drop and, at the same time, a maximum heat transfer global coefficient considering the multi-objective optimization using genetic algorithms to generate
34
F.A.S. Mota et al. / Applied Thermal Engineering 63 (2014) 33e39
Nomenclature AP ATCP A b B cp ci C Cr De DP E ECC f g Gc Gp IS k kP LP _ M M N Nu NC NP NTU OS P Pr Rf Re
plate effective heat transfer area, m2 total area of the PHE, m2 Eigenvalues and Eigenvectors matrix channel average thickness, m binary vector specific heat of the fluid, J/kg K coefficients defined in Eq. (7) ci coefficients vector heat capacity min. and max. ratio equivalent diameter of the channel port diameter, m exchanger thermal effectiveness, % thermal effectiveness of countercurrent flow Fanning factor gravity acceleration, m/s2 channel mass velocity, kg/s m2 port mass velocity, kg/s m2 initial set of configurations fluids thermal conductivity, W/m K plate thermal conductivity, W/m K plate length, m mass flow, kg/s tri-diagonal matrix number of channels per pass Nusselt number number of channels number of plates number of transfer units optimal set number of passes Prandtl number fouling factor, K/W Reynolds number
the Pareto boundary that gives the set of optimal solutions. On a recent study, Arsenyeva et al. [5] developed an optimization method to obtain the optimal design of a PHE with mixed plates. An objective function is formulated, which is an implicit function of the kind of plate, subject to certain constraints (number of passes for both currents, and the number of plates with different grooves). For the optimization problem formulation, the authors developed a mathematical model, which deals with the multipass PHEs as systems of passage plate packs. For a given kind of plate the computing of single-pass arrangements starts by increasing the number of passes whereas the objective function is decreasing. When it starts to increase, the optimal spot is found by comparing the different sets of passes. The global optimum is found by comparing all local optima. In order to avoid the edge effect, the algorithm was limited to the exchangers with more than 21 channels. To avoid infeasible solutions, the maximum number of 10 passes for a stream and the resulting number of passes among the currents being inferior or equal to 50 was established. For such given process conditions, note that the PHEs with mixed plates usually present less area comparing with the ones that use only one type of plate. The current study presents the development of an optimization algorithm based on the screening method, capable of finding a global optimal set taking into account different kinds of plates. For each kind of plate, a local optima set is found, and by comparing all of them the global optimum is discovered. A literature example
RS si T tP U v WP Yh z
reduced set of configurations channel flow direction parameter, si ¼ 1 or þ1 temperature, K plate thickness, m overall heat transfer coefficient, W/m2 K channel fluid velocity, m/s plate width, m binary parameter for hot fluid location tri-diagonal matrix Eigenvector
Greek symbols a dimensionless variable defined with Eq. (3) b chevron corrugation inclination angle, degrees DP pressure drop, Pa h normalized length dimensionless q a dimensional temperature l tri-diagonal matrix Eigenvalue m viscosity, Pa s r density, kg/m3 F enlargement factor of plate f parameter for feed connection Subscripts cold cold fluid hot hot fluid i ith element in fluid inlet j jth element out fluid outlet Superscripts I side I of PHE II side II of PHE max maximum value min minimum value
served as a case study, and its results prove the functioning and performance of the developed algorithm. 2. Configuration characterization To describe the flow distribution of the hot and cold fluids inside the plate pack, five parameters are used: NC, PI, PII, f, Yh. These parameters were defined in Gut and Pinto [3,6]. 2.1. Number of channels (NC) The space between two adjacent plates is a channel. Then the number of channels of a PHE is the number of plates minus one. The
Fig. 1. PHE feed connections and sides.
F.A.S. Mota et al. / Applied Thermal Engineering 63 (2014) 33e39
end plates are not considered (Fig. 1). The odd-numbered channels belong to side I, and the even-numbered ones belong to side II. The number of channels in each side is NCI and NCII . 2.2. Number of passes (P) A single-pass arrangement occurs when a fluid flows in just one direction in the channels in a PHE. The fluid makes just one pass. A multi-pass arrangement occurs when the stream is split and flows in alternate directions. The number of change direction plus one is the number of passes. PI and PII are the number of passes in each side. 2.3. Hot fluid location (Yh) It is a binary parameter that assigns the fluids to the PHE sides. If Yh ¼ 1, the hot fluid flows in side I and if Yh ¼ 0, the hot fluid flows in side II.
The feed of side I is arbitrarily set at h ¼ 0 as presented in Fig. 1. Thus the parameter f represents relative position of the side II. Fig. 1 illustrates all possibilities of connection. The h parameter is defined as h ¼ x/LP.
dq1 ¼ s1 aI ðq2 q1 Þ first channel (1a) dh dqi si aI ðqi1 2qi þ qiþ1 Þ if i is odd ¼ intermediate channels si aII ðqi1 2qi þ qiþ1 Þ if i is even dh (1b) sNc aI qNc 1 þ qNc sNc aII qNc 1 þ qNc
if NC is odd if NC is even
last channel
where s is the direction of flow in the channels (s ¼ 1 if upward flow and s ¼ 1 if downward flow) and:
aI ¼
Ti Tcold;in Thot;in Tcold;in
AP UN I _ I cI M
aII ¼
p
0q1
AP UNII _ II cII M
(2)
(3)
p
AP is the plate area, U is the overall heat transfer coefficient, N is the _ is the mass flow and cp is the number of channels per pass, M specific heat.
U ¼
(6)
where
2
d1 6 þd2 6 6 0 M ¼ 6 6 « 6 4 0 0 8 < si aI :
si aII
þd1 2d2 þd3
0 þd2 2d3
0 0 þd3
/ /
/ /
0 0
þdNC 1 0
2dNC 1 þdNC
0 0 « 0
3
7 7 7 7 7 7 þdNC 1 5 dNC
if i is odd i ¼ 1; .Nc if i is even
The analytical solution for the system is given by Eq. (7a), where
li and zi are, respectively, the Eigenvalues and Eigenvectors of matrix M. Nc X
ci zi eli h
(7a)
1 tP 2b 2b I II FkI NuI þ kP þ FkII NuI þ Rf þ Rf
It is important to point out that for some particular configurations of PHE (where both streams make a number of pass multiple of two), the considered analytical solution is not valid. Zaleski and Jarzebski [8] showed that among the Eigenvalues always occurs at least a null one. For the cases where multiples null Eigenvalues occur, the proposed analytical solution may not be suitable and alternative solutions must be applied. Zaleski [9] presented an alternative solution for selected configurations Eq. (7b).
qðhÞ ¼
NX c 2
ci zi eli h þ cNc 1 zNc 1 h þ cNc zNc
(7b)
i¼1
(1c)
qi ðTi Þ ¼
dq ¼ M$q dh
i¼1
The model development is based on the proposition of Gut and Pinto [6]. From a thermal energy balance in the exchanger channels, the authors reached the following system of normalized ordinary differential equations:
dqNC ¼ dh
(5)
The empirical parameters a1 and a2 are given in Saunders [7]. Since Eq. (1) are a linear system of differential equations, this system may be written in the matrix form, as follows:
qðhÞ ¼
3. Mathematical model
Nu ¼ a1 ðReÞa2 ðPrÞ1=3
di ¼
2.4. Feed connection (f)
35
(4)
Here, b is the channel average thickness, F is the enlargement factor of the plate, k is the thermal conductivity of the fluid, kP is the thermal conductivity of the plate, tP is the plate thickness, Rf is the fouling factor and Nu is the Nusselt number defined as function of Reynolds and Prandtl numbers:
Considering perfect mixture at the end of each pass, the boundaries conditions are given by: 1. Fluid inlet: the fluid inlet temperature in the channels of the first pass is the fluid feed temperature.
qi ðhÞ ¼ qfluid;in i˛first pass
(8a)
2. Change of pass: the temperature at the beginning of the channels of a determined pass is equal to the arithmetic average of the temperatures in the channels of the previous pass.
qi ðhÞ ¼
1 N
N X
qj ðhÞ i˛new pass
(8b)
j ˛ previous pass
3. Fluid outlet: the outlet temperature of the fluid is equal to the arithmetic average of outlet temperatures of the channels of the last pass.
36
F.A.S. Mota et al. / Applied Thermal Engineering 63 (2014) 33e39
qfluid;out ðhÞ ¼
N X
1 N
qj ðhÞ
(8c)
j ˛ last pass
Applying Eq. (7a) into the boundary condition equations for the fluid inlet and pass change it is possible to create a linear system of NC equations whose variables are ci. After solving the linear system, the outlet temperatures can be determined by the use of the outlet boundary conditions and, consequently, the thermal effectiveness can be determined. In order to illustrate the generation of the linear system, a PHE with 7 thermal plates (or 8 channels) with the cold fluid making two passes and the hot one making one pass is presented in Fig. 2. Applying Eq. (7a), the following analytical solution can be achieved:
2
3
2
z1;1 6 z2;1 7 l h 6 1 6 6 qðhÞ ¼ c1 4 « 7 5e þ c2 4 z8;1
3
2
z1;2 6 z2;2 7 7el2 h þ . þ c8 6 4 « 5 z8;2
3
z1;8 z2;8 7 7el8 h « 5 z8;8
(9)
where A ¼ Eigenvalues and Eigenvectors matrix C ¼ c0i s coefficients vector B ¼ binary vector where
Bi ¼ 0
if Bi ¼ qhot;in
i˛first pass
Bi ¼ 1
if Bi ¼ qcold;in
i˛first pass
Bi ¼ 0 if Bi ¼ qi ðhÞ N1
N P
APHE ¼ f NC ; P I ; P II ; f; Yh ; plate type
(10)
Depending on the type and size, the number of plates that a plate heat exchanger can hold may vary from 3 to 700. An initial estimate may be considered and a constraint for the number of plates must be imposed:
Using the boundary conditions Eqs. (8a) and (8b), to all the channels of the considered PHE, the equations presented in Table 1 are generated. Applying Eq. (7a) to the linear system, the matrix form is achieved:
A$C ¼ B
Fig. 2. PHE streams.
qj ðhÞ i˛intermediate pass
j˛previous
NCmin NC NCmax
The pressure drops of both streams must not exceed allowable ones: max DPhot DPhot
(12a)
max DPcold DPcold
(12b)
The pressure drop can be calculated as follows:
DP ¼
For given process conditions, it is possible to create an optimization problem to find the heat exchanger minimum area, which is an implicit function of the PHE configuration parameters and of the plate type:
(13)
GC ¼
_ M NbWP
(14)
GP ¼
_ 4M pD2P
(15)
4. Optimization
4.1. Optimization problem formulation
2f ðLP þ DP ÞPG2C G2 þ 1; 4 P þ rgðLP þ DP Þ rDe 2r
where LP is the plate length, DP is the port diameter, r is the density, g is the gravity acceleration. The remaining parameters f (Fanning factor), GP (port mass velocity), GC (channel mass velocity) and De (equivalent diameter of the channel) are defined as:
pass
The used optimization method is based on screening, which was studied by Gut and Pinto [3], in which e for a given type of plate e the number of thermal plates is the objective function that has to be minimized. This is subject to certain inequality constraints and is implicitly given by the heat exchanger configuration parameters. In this work, in order to obtain the global optimum, the screening method must be applied to each type of plate, generating local optima. It is done simply with a loop applied to the screening method sweeping all kinds of plates. The global optimum is found by comparing the local optima.
(11)
f ¼
Kp Rem
De ¼
(16)
2b
(17)
F
The parameters Kp and m are given in Saunders [7]. The velocity constraints for the streams are taken into account to avoid the formation of dead spaces and/or bubbles inside the
Table 1 Boundary condition equations. Cold fluid
Hot fluid
q1(h ¼ 0) ¼ [q5(h ¼ 0) þ q7(h ¼ 0)]/2 q3(h ¼ 0) ¼ [q5(h ¼ 0) þ q7(h ¼ 0)]/2 q5(h ¼ 1) ¼ qcold,in ¼ 0 q7(h ¼ 1) ¼ qcold,in ¼ 0
q2(h ¼ 1) ¼ qhot,in ¼ 1 q4(h ¼ 1) ¼ qhot,in ¼ 1 q6(h ¼ 1) ¼ qhot,in ¼ 1 q8(h ¼ 1) ¼ qhot,in ¼ 1
F.A.S. Mota et al. / Applied Thermal Engineering 63 (2014) 33e39
plate pack. According to Kakaç and Liu [10], flow velocities lower than 0.1 m/s are not used in practice.
vhot vmin hot
(18a)
vcold vmin cold
(18b)
The fluid velocity can be calculated as follows:
v ¼
_ M bWP r
(19)
Fluid properties in Eqs. (13)e(19) are calculated at the average fluid temperature. Depending on the required thermal load or outlet temperatures, a thermal effectiveness constraint may be established:
Emin E Emax
(20)
In order to solve this optimization problem, the screening method proposed by Gut and Pinto [3] is adopted. In this method, the constraints are evaluated successively, reducing the number of configurations until the optima set of configurations (OS) is found, in case it exists. The screening process starts with the identification of the initial set of possible configurations (IS). This initial set identification is based on the number of channels constraint. A subset of IS called reduced set of configurations (RS) is created by verifying both flow velocity and pressure drop constraints. The thermal effectiveness constraint is applied to the reduced set in an increasing order of number of channels. The configurations with the lowest number of channels, respecting the thermal constraint, will form the local optima set. The global optimum is found by comparing the set of local optima. Below it is presented the optimization algorithm description. 4.2. Optimization algorithm description In the paper of Gut and Pinto [3] the screening method was programmed in C language to obtain the reduced set of configuration (RS), but the simulation of the systems of differential equations of the reduced set of configurations (RS) was performed numerically using the gPROMS software. Different from the paper of Gut and Pinto [3], in the present paper, the proposed algorithm is fully developed in MatLab, allowing the inclusion of more optimization variables. Also, it is not necessary the use of auxiliary software, like gPROMS. It represents the total independence of additional software to solve the systems of differential equations. As presented in Eq. (10), the types of plates are also optimization variables in this work. In addition, the simulation of the configurations is performed analytically, differently from the work of Gut and Pinto [3]. The developed algorithm is presented as follows: Algorithm 1. Step
1. Input data PHE parameters: LP ; W P ; t P ; DP ; b; F; kP ; b _ Hot fluid characteristics: Thot;in ; M hot ; Rf ;hot _ Cold fluid characteristics: Tcold;in ; M cold ; Rf ;cold Hot fluid average physical proprieties: rhot, mhot, CP,hot, khot Cold fluid average physical proprieties: rcold, mcold, CP,cold, kcold Constraints: max ; vmin ; DP max ; vmin NCmin ; NCmax ; Emin ; Emax ; DPhot hot cold cold
Note that the PHE parameters are vectors that sweep all types of plates considered.
37
Step 2. Generation of the configurations initial set (IS). The vector NC is generated with all possible numbers of channels for each type of plate.
NC ¼ NCmin : NCmax Step 3. For each element of the vector N C all possible number of passes to both sides are computed, which are integer divisors of the number of channels of the corresponding side. Step 4. Verification of the velocity and pressure drop constraints, v and DP, respectively. Step 4.1. Consider Yh ¼ 0 (hot fluid goes through side II whereas I II cold fluid goes through side I). Thus, Pcold ¼ P I and Phot ¼ P II . Step 4.1.1. For each type of plate calculate the cold fluid velocity, vIcold , in descending order of the possible numbers of passes of a given element of NC . If vIcold reach the minimum allowed value, it is not necessary to evaluate configurations with lower numbers of passes. Step 4.1.2. For each type of plate the pressure drop of the cold I fluid must be calculated, DPcold , in ascending order of the I possible numbers of passes of a given element of N C . If DPcold reaches the maximum allowed value, it is not necessary to evaluate configurations with higher numbers of passes. Step 4.1.3. Verify the flow velocity constraint of side II. It is the same as in Step 4.1.1. Step 4.1.4. Verify the pressure drop constraint of side II. It is the same as in Step 4.1.2. Step 4.2. Considerate Yh ¼ 1. It is the same as in Step 4.1. Step 5. Generate the reduced set of configurations (RS). Step 5.1. For Yh ¼ 0, the number of selected passes for sides I and II, for the heat exchanger, are matched, achieving configurations in the format of [NC PI PII Yh]. For each combination of these 4 parameters, 4 different configurations are possible according to the feed arrangement (f parameter). Step 5.2. The same procedure (Step 5.1) is carried out for Yh ¼ 1. Thus it is calculated the remaining configurations that formed the reduced set. Step 6. Calculation of the thermal effectiveness, ECC, considering purely countercurrent flow. If ECC < Emin, than such configurations may be discarded, because due to edge effects and changes of passes, the thermal effectiveness of a PHE is always less than ECC.
ECC ¼
8 > > > < > > > :
1eNTUð1Cr Þ 1Cr eNTUð1Cr Þ
if Cr < 1
NTU NTUþ1
if Cr ¼ 1
(21)
Step 7. Verification of the thermal effectiveness constraint. The selected configurations on Step 5 are simulated in an increasing order of number of channels until the possible optimal set is found. The remaining configurations do not need to be simulated. In the next section, the simulation algorithm for all the reduced set of configurations (RS) is described. Step 8. Achievement of the global optimum. The heat exchanger with the minimal area is determined by comparing all local optima of all types of plates considered. 4.2.1. Simulation algorithm For the development of this algorithm, the boundary conditions equations in the algorithm form given by Gut and Pinto [6] are used. The simulation algorithm is applied for each value of f separately. The algorithm is presented below:
38
F.A.S. Mota et al. / Applied Thermal Engineering 63 (2014) 33e39
Algorithm 2. Step 1. The tri-diagonal matrix coefficients are computed, d(i) ¼ s(i)a(I or II). Step 2. The tri-diagonal matrix is constructed. Step 3. Eigenvalues and Eigenvectors are computed from the tridiagonal matrix. Step 4. The linear system is generated (Eq. (9)). The boundary condition equations in algorithmic form of inlet fluids and change of passes are used. Step 4.1. Generation of the side I portion of Eigenvalues and I I Eigenvectors matrix, A , and of the binary vector, B . Step 4.2. Generation of the side II portion of Eigenvalues and II II Eigenvectors matrix, A , and of the binary vector, B . Step 4.3. Computation of the Eigenvalues and Eigenvectors matrix and computation of the binary vector. Thus the linear system is identified. I
II
I
II
A ¼ A þA B ¼ B þB
Step 5. Determination of the c0i s coefficients by solving the linear system. Step 6. Determination of the output dimensionless temperatures. The boundary conditions of output fluid in algorithmic form is used. Step 7. Computation of the thermal effectiveness. It may be obtained considering any side of the PHE, because the energy conservation is obeyed only if E ¼ EI ¼ EII. Note that for different values of f different values of E can be achieved because NI and NII are dependent on f.
8 I > I ¼ NI max aI ; aII jq q > E > out j in aI NI NII < E ¼ > > > : EII ¼ NIIII max aII ; aIIII jqin qout jII a N N
(22)
Table 3 Thermal effectiveness of RS for b ¼ 60 . #
NC
PI
PII
Yh
ECC
E (f ¼ 1)
E (f ¼ 2)
E (f ¼ 3)
E (f ¼ 4)
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22
40 40 41 41 42 42 43 43 44 44 48 48 48 48 48 48 48 48 48 48 48 48
2 2 3 3 3 3 2 2 2 2 2 2 2 3 3 3 2 2 3 3 4 4
2 2 2 2 3 3 3 3 2 2 2 3 4 2 3 4 2 3 2 3 2 3
0 1 0 1 0 1 0 1 0 1 0 0 0 0 0 0 1 1 1 1 1 1
0.83 0.83 0.86 0.86 0.88 0.88 0.86 0.86 0.84 0.84 0.85 0.87 0.88 0.86 0.89 0.90 0.85 0.86 0.87 0.89 0.88 0.90
0.538 0.538 0.540 0.539 0.534 0.534 0.536 0.538 0.537 0.537 0.537 0.538 0.515 0.535 0.534 0.531 0.537 0.540 0.533 0.534 0.515 0.530
ea ea 0.540 0.539 0.535 0.545 0.536 0.538 0.492 0.492 ea 0.532 0.514 0.540 0.537 0.530 ea 0.535 0.538 0.537 0.514 0.531
0.819 0.819 0.770 0.774 0.774 0.774 0.776 0.772 0.826 0.826 0.833 0.776 0.795 0.772 0.778 0.840 0.833 0.773 0.776 0.778 0.795 0.840
0.691 0.691 0.770 0.774 0.860b 0.860b 0.776 0.772 ea 0.693 ea 0.786 ea 0.782 0.870 0.830 ea 0.782 0.786 0.870 ea 0.830
E (f ¼ 3)
E (f ¼ 4)
a b
Table 4 Thermal effectiveness of RS for b ¼ 45 . # 1 2 3 4 5 6 7 8 a b
5. Results As a case study, an example of optimization was considered for three types of chevron plates with the same dimensions, but with different angles of corrugation pattern, b. Table 2 presents the Table 2 Example data. Plate characteristics
b ¼ [60 45 30 ] F ¼ 1.24
LP ¼ 0.994 m WP ¼ 0.525 m b ¼ 2.85 mm DP ¼ 15 cm
tP ¼ 0.5 mm kP ¼ 17 W/m K
Hot water inlet
Cold water inlet
Tin,hot ¼ 110 C _ M hot ¼ 20118.91 kg/h Rf,hot ¼ 0 cm2 K/W rhot ¼ 965 kg/m3 mhot ¼ 0.315 cP cP,hot ¼ 4.205 J/kg K khot ¼ 0.673 W/m K
Tin,cold ¼ 65 C _ M cold ¼ 23031.38 kg/h Rf,cold ¼ 0 cm2 K/W rcold ¼ 970 kg/m3 mcold ¼ 0.344 cP cP,cold ¼ 4.198 J/kg K kcold ¼ 0.668 W/m K
E >85% 0 DPcold 1 bar
¼ 0.3 m/s vmin hot
vmin ¼ 0.3 m/s cold
NC 40 40 44 44 48 48 48 48
PI 2 2 2 2 2 2 2 3
PII 2 2 2 2 2 3 2 2
Yh 0 1 0 1 0 0 1 1
ECC 0.92 0.92 0.92 0.92 0.93 0.94 0.93 0.94
E (f ¼ 1) 0.539 0.539 0.538 0.538 0.538 0.521 0.538 0.507
E (f ¼ 2) a
e ea ea ea ea 0.507 0.386 0.521
b
0.894 0.894b 0.900 0.900 0.906 0.816 0.906 0.816
0.701 0.701 ea 0.702 ea 0.832 ea 0.832
Proposed analytical solution cannot be applied. Optimal configuration.
example data, which were extracted from Brenner et al. [11], with a modification in the available pressure to 10% of the original value, just to work with a smaller set of configurations. It makes the problem analysis less arduous. The original problem also considered the input and output temperatures as thermal constraints. A minimal thermal effectiveness of 0.85 was considered. The minimal number of channels considered is 40. Usually the plate heat exchangers work with a greater number of channels. With the application of the optimization algorithm until Step 5 for the first element of the vector b, the RS is selected, which is made up of 88 configurations respecting the constraint of the number of plates, pressure drop, and flow velocity (see Table 3). The sixth step of the algorithm discards configurations 1, 2, 9 and 10, because ECC <85%. The seventh step identified two optimal configurations with thermal effectiveness of 86.0% with 42 channels
Table 5 Thermal effectiveness of RS for b ¼ 30 . #
NC
PI
PII
Yh
ECC
E (f ¼ 1)
E (f ¼ 2)
E (f ¼ 3)
E (f ¼ 4)
1 2
48 48
2 2
2 2
0 1
0.94 0.94
0.538 0.538
0.364 ea
0.917b 0.917b
ea ea
Constraints 40 NC 50 0 DPhot 1 bar
Proposed analytical solution cannot be applied. Optimal configuration.
a b
Proposed analytical solution cannot be applied. Optimal configuration.
F.A.S. Mota et al. / Applied Thermal Engineering 63 (2014) 33e39
6. Conclusions
Table 6 Local and global optima.
b
NC
f
PI
PII
Yh
DPhot (bar)
DPcold (bar)
E
60
42 42 40 40 48 48
4 4 3 3 3 3
3 3 2 2 2 2
3 3 2 2 2 2
0 1 0 1 0 1
0.57 0.57 0.44 0.44 0.72 0.72
0.71 0.71 0.54 0.54 0.90 0.90
0.860 0.860 0.894 0.894 0.917 0.917
45 30
39
Values for b ¼ 45 (in bold) correspond to the global optima.
and each stream making 3 passes. The difference between the optimal configurations is the fluids location. One of them considers the hot fluid flowing in side I and the other one considers the cold fluid in the side I. Just for comparison reasons, Table 3 presents all simulations for the reduced set of configurations (RS). Once the optimal set is found the remaining configurations may be discarded. The cells identified by asterisk marks (*) mean that in simulation of such configurations, a pair of null Eigenvalues came up between the Eigenvalues of the tri-diagonal matrix. In this case the proposed analytic solution in Eq. (7a) cannot be applied. An attempt to use Eq. (7b) was done, but the correct solution was not obtained as well. The reduced set of configurations for the next element of the vector b is presented in Table 4. It were selected 32 configurations, amount smaller than that one found for the first element of b. However, it was expected because the smaller the angle the greater the pressure drops. The sixth step of the algorithm does not discard any configuration, as may be observed in Table 4. In this case ECC >85% for all RS configurations. When Step 6 does not eliminate any element of RS, it is because the value for the constraint of minimal number of channels was not adequate. With a smaller number of channels, the thermal constraint could already be obeyed. In this particular case, for example, it was verified that 28 channels with a 2-2 arrangement would be sufficient to meet the thermal constraint. By applying the thermal constraint (Step 7), the local optima set was identified in the two first configurations, with 40 channels, symmetrical arrangements of 2 passes. For the last element of the vector b, the selected RS are also the optima configurations (see Table 5). For this case, however, no simulation would need to be carried out, because no configuration could be the global optimum since the RS configurations have more channels than the ones previously computed. Finally, by comparing all local optima, the global optimum set are identified for b ¼ 45 (Table 6). It is interesting to note that the global optimum set was the one that used the least of the available pressure drop. In order to obtain a better performance for the available pressure drop and, consequently, configurations with fewer plates, plates with different corrugation patterns in the PHE would be a possibility. The configurations with f ¼ 3 will have always equal or better thermal performance because with f ¼ 3 and an even number of passes the heat transfer behavior in the heat exchanger is similar to the countercurrent, that is the most thermal efficient system.
This paper presented an optimization algorithm for the design of plate heat exchangers. It was based on the model proposed by Gut and Pinto [3]. The difference between the present paper and the work of Gut and Pinto [3] is that for a given process and type of plate the proposed algorithm is able to find the configuration of the lowest number of thermal plates (the optimal configuration for a given PHE). The types of plates are considered as optimization variables. The inclusion of this possibility into a systematic procedure is a great novelty in the paper. For each type of plate, a local optima set may be found (the configurations presenting the lowest number of plates). By comparing all the local optima sets, the global optimum set may be found. The system of differential equations was solved analytically, avoiding numerical methods or other approximations. Also, it is important to highlight that no additional software needs to be used, like gPROMS or other ones. This is the great contribution of the present paper when compared with the work of Gut and Pinto [3]. As a case study, an example of the literature was used to show the applicability of the developed algorithm. The example of optimization has considered three types of chevron plates with the same dimensions, but with different angles of corrugation pattern, b. By the optimization example presented, although some configurations cannot be solved by analytical solution proposed, the developed algorithm had great success. The complex Eigenvalues compromised the simulation of some configurations where both streams make a number of pass multiple of two. However, the configuration at f ¼ 3 generally presents best thermal performance for these cases. So, the fact that the proposed analytical solution cannot be applied for such configurations is not a limiting factor for the algorithm use. References [1] L. Wang, B. Sundén, R.M. Manglik, Plate Heat Exchangers: Design, Applications and Performance, WIT Press, Ashurst Lodge, 2007, pp. 27e39. [2] R.K. Shah, D.P. Sekulic, Fundamentals of Heat Exchanger Design, John Wiley & Sons, Inc., New Jersey, 2003, pp. 22e29. [3] J.A.W. Gut, J.M. Pinto, Optimal configuration design for plate heat exchangers, Int. J. Heat Mass Transfer 47 (2004) 4833e4848. [4] H. Najafi, B. Najafi, Multi-objective optimization of a plate and frame heat exchanger via genetic algorithm, Heat Mass Transfer 46 (2010) 639e647. [5] O.P. Arsenyeva, L.L. Tovazhnyansky, P.O. Kapustenko, G.L. Khavin, Optimal design of plate-and-frame heat exchangers for efficient heat recovery in process industries, Energy (2011) 1e11. [6] J.A.W. Gut, J.M. Pinto, Modeling of plate heat exchangers with generalized configurations, Int. J. Heat Mass Transfer 46 (14) (2003) 2571e2585. [7] E.A.D. Saunders, Heat Exchangers: Selection, Design & Construction, Longman S & T, Harlow (UK), 1988. [8] T. Zaleski, A.B. Jarzebski, Remarks on some properties of the equation of heat transfer in multichannel exchangers, Int. J. Heat Mass Transfer 16 (8) (1973) 1527e1530. [9] T. Zaleski, A general mathematical model of parallel-flow, multichannel heat exchangers and analysis of its properties, Chem. Eng. Sci. 39 (7/8) (1984) 1251e1260. [10] S. Kakaç, H. Liu, Heat Exchanger: Selection, Rating and Thermal Design, second ed., CRC Press, Boca Raton, 2002, pp. 373e412. [11] K. Brenner, G. Conrads, J. Dorn, G. Fritz, Z.G. Guo, N. Jahnke, W. Tinz, E. Weib, Plate Heat Exchangers: Installation and Use of Multi-purpose Heat Exchange Systems, Bosch-Druck GmbH, Landshut, 1998.
Contents lists available at ScienceDirect
Applied Thermal Engineering journal homepage: www.elsevier.com/locate/apthermeng
Optimal design of plate heat exchangers Fábio A.S. Mota, Mauro A.S.S. Ravagnani*, E.P. Carvalho Chemical Engineering Graduate Studies Program, State University of Maringá, Av. Colombo, 5790 Maringá, PR, Brazil
h i g h l i g h t s The paper presents an algorithm Q1 for the optimization of heat exchange area of plate heat exchangers (PHE). The optimization algorithm is based on the screening method. Optimal configurations are found for each kind of plate with the lowest number of plates and the smallest area. For all kinds of plates considered, their respective local optima are determined. The global optimum is obtained by comparison among the local optima. The differential equations are solved analytically.
a r t i c l e i n f o
a b s t r a c t
Article history: Received 4 March 2013 Accepted 25 September 2013 Available online 23 October 2013
In the current paper, it is presented an algorithm for the optimization of heat exchange area of plate heat exchangers (PHE). The algorithm is based on the screening method. For each kind of plate, subject to certain constraints, optimal configurations are found e in case they exist e which present the lowest number of plates and, consequently, the smallest area. Each of these found configurations have local optima characteristics. For all kinds of plates considered, their respective local optima are determined. Comparison of all obtained local optima gives a global one. The differential equations, which are generated in the model development, are solved analytically. A case study is presented to test the applicability of the developed algorithm. Results show coherency with the literature. 2013 Elsevier Ltd. All rights reserved.
Keywords: Plate heat exchangers Optimization Simulation Screening
1. Introduction The competitive pressures of the global market and the increasing urgency for both energy conservation and reduction of environmental impacts have changed the focus on the industrial processes, so that heat exchangers with high effectiveness can be used. Plate heat exchanger (PHE) is categorized at the lower end of the compactness spectrum. It offers many advantages and unique characteristics applications when compared with other highly compact heat exchangers. This is due to the flexible thermal design (the plates can be simply added or removed to attend different demands of thermal load and processing), ease to clean maintaining extreme hygiene (necessary when food, pharmaceutical or other kind of products are processed), good temperature control (necessary in cryogenic uses and mitigation of thermal degradation of the process fluids) and better performance in heat transfer [1]. Because of their
* Corresponding author. E-mail addresses: [email protected], [email protected] (M.A.S. S. Ravagnani). 1359-4311/$ e see front matter 2013 Elsevier Ltd. All rights reserved. http://dx.doi.org/10.1016/j.applthermaleng.2013.09.046
characteristics, PHE became ideal for diary, pharmaceutical, food and drink industries. Besides is being used in synthetic rubber industry, pulp and paper industry, maritime applications, solar energy etc. [1,2]. Gut and Pinto [3] proposed the screening method to find the optimal set of configurations of a given plate heat exchanger. In this method, the objective function e which implicitly depends on six configuration parameters e is formulated to obtain the minimum number of channels, which are proportional to the heat transfer area. Constraints are given in order to avoid impossible or non-optimal solutions. The initial set of configurations (IS) is given by the maximum number of possible configurations for the problem and is generated by applying the constraint of channel numbers. The reduced setting (RS) is identified by verifying the pressure drop and flow rate, that is, the IS subset which obeys hydraulic constraints. The optimal set (OS) is formed by the configurations with the minimum number of channels of the RS obeying the thermal constraint. Najafi H. and Najafi B. [4] developed a method to obtain an optimal set of the PHE geometric parameters, with the minimum pressure drop and, at the same time, a maximum heat transfer global coefficient considering the multi-objective optimization using genetic algorithms to generate
34
F.A.S. Mota et al. / Applied Thermal Engineering 63 (2014) 33e39
Nomenclature AP ATCP A b B cp ci C Cr De DP E ECC f g Gc Gp IS k kP LP _ M M N Nu NC NP NTU OS P Pr Rf Re
plate effective heat transfer area, m2 total area of the PHE, m2 Eigenvalues and Eigenvectors matrix channel average thickness, m binary vector specific heat of the fluid, J/kg K coefficients defined in Eq. (7) ci coefficients vector heat capacity min. and max. ratio equivalent diameter of the channel port diameter, m exchanger thermal effectiveness, % thermal effectiveness of countercurrent flow Fanning factor gravity acceleration, m/s2 channel mass velocity, kg/s m2 port mass velocity, kg/s m2 initial set of configurations fluids thermal conductivity, W/m K plate thermal conductivity, W/m K plate length, m mass flow, kg/s tri-diagonal matrix number of channels per pass Nusselt number number of channels number of plates number of transfer units optimal set number of passes Prandtl number fouling factor, K/W Reynolds number
the Pareto boundary that gives the set of optimal solutions. On a recent study, Arsenyeva et al. [5] developed an optimization method to obtain the optimal design of a PHE with mixed plates. An objective function is formulated, which is an implicit function of the kind of plate, subject to certain constraints (number of passes for both currents, and the number of plates with different grooves). For the optimization problem formulation, the authors developed a mathematical model, which deals with the multipass PHEs as systems of passage plate packs. For a given kind of plate the computing of single-pass arrangements starts by increasing the number of passes whereas the objective function is decreasing. When it starts to increase, the optimal spot is found by comparing the different sets of passes. The global optimum is found by comparing all local optima. In order to avoid the edge effect, the algorithm was limited to the exchangers with more than 21 channels. To avoid infeasible solutions, the maximum number of 10 passes for a stream and the resulting number of passes among the currents being inferior or equal to 50 was established. For such given process conditions, note that the PHEs with mixed plates usually present less area comparing with the ones that use only one type of plate. The current study presents the development of an optimization algorithm based on the screening method, capable of finding a global optimal set taking into account different kinds of plates. For each kind of plate, a local optima set is found, and by comparing all of them the global optimum is discovered. A literature example
RS si T tP U v WP Yh z
reduced set of configurations channel flow direction parameter, si ¼ 1 or þ1 temperature, K plate thickness, m overall heat transfer coefficient, W/m2 K channel fluid velocity, m/s plate width, m binary parameter for hot fluid location tri-diagonal matrix Eigenvector
Greek symbols a dimensionless variable defined with Eq. (3) b chevron corrugation inclination angle, degrees DP pressure drop, Pa h normalized length dimensionless q a dimensional temperature l tri-diagonal matrix Eigenvalue m viscosity, Pa s r density, kg/m3 F enlargement factor of plate f parameter for feed connection Subscripts cold cold fluid hot hot fluid i ith element in fluid inlet j jth element out fluid outlet Superscripts I side I of PHE II side II of PHE max maximum value min minimum value
served as a case study, and its results prove the functioning and performance of the developed algorithm. 2. Configuration characterization To describe the flow distribution of the hot and cold fluids inside the plate pack, five parameters are used: NC, PI, PII, f, Yh. These parameters were defined in Gut and Pinto [3,6]. 2.1. Number of channels (NC) The space between two adjacent plates is a channel. Then the number of channels of a PHE is the number of plates minus one. The
Fig. 1. PHE feed connections and sides.
F.A.S. Mota et al. / Applied Thermal Engineering 63 (2014) 33e39
end plates are not considered (Fig. 1). The odd-numbered channels belong to side I, and the even-numbered ones belong to side II. The number of channels in each side is NCI and NCII . 2.2. Number of passes (P) A single-pass arrangement occurs when a fluid flows in just one direction in the channels in a PHE. The fluid makes just one pass. A multi-pass arrangement occurs when the stream is split and flows in alternate directions. The number of change direction plus one is the number of passes. PI and PII are the number of passes in each side. 2.3. Hot fluid location (Yh) It is a binary parameter that assigns the fluids to the PHE sides. If Yh ¼ 1, the hot fluid flows in side I and if Yh ¼ 0, the hot fluid flows in side II.
The feed of side I is arbitrarily set at h ¼ 0 as presented in Fig. 1. Thus the parameter f represents relative position of the side II. Fig. 1 illustrates all possibilities of connection. The h parameter is defined as h ¼ x/LP.
dq1 ¼ s1 aI ðq2 q1 Þ first channel (1a) dh dqi si aI ðqi1 2qi þ qiþ1 Þ if i is odd ¼ intermediate channels si aII ðqi1 2qi þ qiþ1 Þ if i is even dh (1b) sNc aI qNc 1 þ qNc sNc aII qNc 1 þ qNc
if NC is odd if NC is even
last channel
where s is the direction of flow in the channels (s ¼ 1 if upward flow and s ¼ 1 if downward flow) and:
aI ¼
Ti Tcold;in Thot;in Tcold;in
AP UN I _ I cI M
aII ¼
p
0q1
AP UNII _ II cII M
(2)
(3)
p
AP is the plate area, U is the overall heat transfer coefficient, N is the _ is the mass flow and cp is the number of channels per pass, M specific heat.
U ¼
(6)
where
2
d1 6 þd2 6 6 0 M ¼ 6 6 « 6 4 0 0 8 < si aI :
si aII
þd1 2d2 þd3
0 þd2 2d3
0 0 þd3
/ /
/ /
0 0
þdNC 1 0
2dNC 1 þdNC
0 0 « 0
3
7 7 7 7 7 7 þdNC 1 5 dNC
if i is odd i ¼ 1; .Nc if i is even
The analytical solution for the system is given by Eq. (7a), where
li and zi are, respectively, the Eigenvalues and Eigenvectors of matrix M. Nc X
ci zi eli h
(7a)
1 tP 2b 2b I II FkI NuI þ kP þ FkII NuI þ Rf þ Rf
It is important to point out that for some particular configurations of PHE (where both streams make a number of pass multiple of two), the considered analytical solution is not valid. Zaleski and Jarzebski [8] showed that among the Eigenvalues always occurs at least a null one. For the cases where multiples null Eigenvalues occur, the proposed analytical solution may not be suitable and alternative solutions must be applied. Zaleski [9] presented an alternative solution for selected configurations Eq. (7b).
qðhÞ ¼
NX c 2
ci zi eli h þ cNc 1 zNc 1 h þ cNc zNc
(7b)
i¼1
(1c)
qi ðTi Þ ¼
dq ¼ M$q dh
i¼1
The model development is based on the proposition of Gut and Pinto [6]. From a thermal energy balance in the exchanger channels, the authors reached the following system of normalized ordinary differential equations:
dqNC ¼ dh
(5)
The empirical parameters a1 and a2 are given in Saunders [7]. Since Eq. (1) are a linear system of differential equations, this system may be written in the matrix form, as follows:
qðhÞ ¼
3. Mathematical model
Nu ¼ a1 ðReÞa2 ðPrÞ1=3
di ¼
2.4. Feed connection (f)
35
(4)
Here, b is the channel average thickness, F is the enlargement factor of the plate, k is the thermal conductivity of the fluid, kP is the thermal conductivity of the plate, tP is the plate thickness, Rf is the fouling factor and Nu is the Nusselt number defined as function of Reynolds and Prandtl numbers:
Considering perfect mixture at the end of each pass, the boundaries conditions are given by: 1. Fluid inlet: the fluid inlet temperature in the channels of the first pass is the fluid feed temperature.
qi ðhÞ ¼ qfluid;in i˛first pass
(8a)
2. Change of pass: the temperature at the beginning of the channels of a determined pass is equal to the arithmetic average of the temperatures in the channels of the previous pass.
qi ðhÞ ¼
1 N
N X
qj ðhÞ i˛new pass
(8b)
j ˛ previous pass
3. Fluid outlet: the outlet temperature of the fluid is equal to the arithmetic average of outlet temperatures of the channels of the last pass.
36
F.A.S. Mota et al. / Applied Thermal Engineering 63 (2014) 33e39
qfluid;out ðhÞ ¼
N X
1 N
qj ðhÞ
(8c)
j ˛ last pass
Applying Eq. (7a) into the boundary condition equations for the fluid inlet and pass change it is possible to create a linear system of NC equations whose variables are ci. After solving the linear system, the outlet temperatures can be determined by the use of the outlet boundary conditions and, consequently, the thermal effectiveness can be determined. In order to illustrate the generation of the linear system, a PHE with 7 thermal plates (or 8 channels) with the cold fluid making two passes and the hot one making one pass is presented in Fig. 2. Applying Eq. (7a), the following analytical solution can be achieved:
2
3
2
z1;1 6 z2;1 7 l h 6 1 6 6 qðhÞ ¼ c1 4 « 7 5e þ c2 4 z8;1
3
2
z1;2 6 z2;2 7 7el2 h þ . þ c8 6 4 « 5 z8;2
3
z1;8 z2;8 7 7el8 h « 5 z8;8
(9)
where A ¼ Eigenvalues and Eigenvectors matrix C ¼ c0i s coefficients vector B ¼ binary vector where
Bi ¼ 0
if Bi ¼ qhot;in
i˛first pass
Bi ¼ 1
if Bi ¼ qcold;in
i˛first pass
Bi ¼ 0 if Bi ¼ qi ðhÞ N1
N P
APHE ¼ f NC ; P I ; P II ; f; Yh ; plate type
(10)
Depending on the type and size, the number of plates that a plate heat exchanger can hold may vary from 3 to 700. An initial estimate may be considered and a constraint for the number of plates must be imposed:
Using the boundary conditions Eqs. (8a) and (8b), to all the channels of the considered PHE, the equations presented in Table 1 are generated. Applying Eq. (7a) to the linear system, the matrix form is achieved:
A$C ¼ B
Fig. 2. PHE streams.
qj ðhÞ i˛intermediate pass
j˛previous
NCmin NC NCmax
The pressure drops of both streams must not exceed allowable ones: max DPhot DPhot
(12a)
max DPcold DPcold
(12b)
The pressure drop can be calculated as follows:
DP ¼
For given process conditions, it is possible to create an optimization problem to find the heat exchanger minimum area, which is an implicit function of the PHE configuration parameters and of the plate type:
(13)
GC ¼
_ M NbWP
(14)
GP ¼
_ 4M pD2P
(15)
4. Optimization
4.1. Optimization problem formulation
2f ðLP þ DP ÞPG2C G2 þ 1; 4 P þ rgðLP þ DP Þ rDe 2r
where LP is the plate length, DP is the port diameter, r is the density, g is the gravity acceleration. The remaining parameters f (Fanning factor), GP (port mass velocity), GC (channel mass velocity) and De (equivalent diameter of the channel) are defined as:
pass
The used optimization method is based on screening, which was studied by Gut and Pinto [3], in which e for a given type of plate e the number of thermal plates is the objective function that has to be minimized. This is subject to certain inequality constraints and is implicitly given by the heat exchanger configuration parameters. In this work, in order to obtain the global optimum, the screening method must be applied to each type of plate, generating local optima. It is done simply with a loop applied to the screening method sweeping all kinds of plates. The global optimum is found by comparing the local optima.
(11)
f ¼
Kp Rem
De ¼
(16)
2b
(17)
F
The parameters Kp and m are given in Saunders [7]. The velocity constraints for the streams are taken into account to avoid the formation of dead spaces and/or bubbles inside the
Table 1 Boundary condition equations. Cold fluid
Hot fluid
q1(h ¼ 0) ¼ [q5(h ¼ 0) þ q7(h ¼ 0)]/2 q3(h ¼ 0) ¼ [q5(h ¼ 0) þ q7(h ¼ 0)]/2 q5(h ¼ 1) ¼ qcold,in ¼ 0 q7(h ¼ 1) ¼ qcold,in ¼ 0
q2(h ¼ 1) ¼ qhot,in ¼ 1 q4(h ¼ 1) ¼ qhot,in ¼ 1 q6(h ¼ 1) ¼ qhot,in ¼ 1 q8(h ¼ 1) ¼ qhot,in ¼ 1
F.A.S. Mota et al. / Applied Thermal Engineering 63 (2014) 33e39
plate pack. According to Kakaç and Liu [10], flow velocities lower than 0.1 m/s are not used in practice.
vhot vmin hot
(18a)
vcold vmin cold
(18b)
The fluid velocity can be calculated as follows:
v ¼
_ M bWP r
(19)
Fluid properties in Eqs. (13)e(19) are calculated at the average fluid temperature. Depending on the required thermal load or outlet temperatures, a thermal effectiveness constraint may be established:
Emin E Emax
(20)
In order to solve this optimization problem, the screening method proposed by Gut and Pinto [3] is adopted. In this method, the constraints are evaluated successively, reducing the number of configurations until the optima set of configurations (OS) is found, in case it exists. The screening process starts with the identification of the initial set of possible configurations (IS). This initial set identification is based on the number of channels constraint. A subset of IS called reduced set of configurations (RS) is created by verifying both flow velocity and pressure drop constraints. The thermal effectiveness constraint is applied to the reduced set in an increasing order of number of channels. The configurations with the lowest number of channels, respecting the thermal constraint, will form the local optima set. The global optimum is found by comparing the set of local optima. Below it is presented the optimization algorithm description. 4.2. Optimization algorithm description In the paper of Gut and Pinto [3] the screening method was programmed in C language to obtain the reduced set of configuration (RS), but the simulation of the systems of differential equations of the reduced set of configurations (RS) was performed numerically using the gPROMS software. Different from the paper of Gut and Pinto [3], in the present paper, the proposed algorithm is fully developed in MatLab, allowing the inclusion of more optimization variables. Also, it is not necessary the use of auxiliary software, like gPROMS. It represents the total independence of additional software to solve the systems of differential equations. As presented in Eq. (10), the types of plates are also optimization variables in this work. In addition, the simulation of the configurations is performed analytically, differently from the work of Gut and Pinto [3]. The developed algorithm is presented as follows: Algorithm 1. Step
1. Input data PHE parameters: LP ; W P ; t P ; DP ; b; F; kP ; b _ Hot fluid characteristics: Thot;in ; M hot ; Rf ;hot _ Cold fluid characteristics: Tcold;in ; M cold ; Rf ;cold Hot fluid average physical proprieties: rhot, mhot, CP,hot, khot Cold fluid average physical proprieties: rcold, mcold, CP,cold, kcold Constraints: max ; vmin ; DP max ; vmin NCmin ; NCmax ; Emin ; Emax ; DPhot hot cold cold
Note that the PHE parameters are vectors that sweep all types of plates considered.
37
Step 2. Generation of the configurations initial set (IS). The vector NC is generated with all possible numbers of channels for each type of plate.
NC ¼ NCmin : NCmax Step 3. For each element of the vector N C all possible number of passes to both sides are computed, which are integer divisors of the number of channels of the corresponding side. Step 4. Verification of the velocity and pressure drop constraints, v and DP, respectively. Step 4.1. Consider Yh ¼ 0 (hot fluid goes through side II whereas I II cold fluid goes through side I). Thus, Pcold ¼ P I and Phot ¼ P II . Step 4.1.1. For each type of plate calculate the cold fluid velocity, vIcold , in descending order of the possible numbers of passes of a given element of NC . If vIcold reach the minimum allowed value, it is not necessary to evaluate configurations with lower numbers of passes. Step 4.1.2. For each type of plate the pressure drop of the cold I fluid must be calculated, DPcold , in ascending order of the I possible numbers of passes of a given element of N C . If DPcold reaches the maximum allowed value, it is not necessary to evaluate configurations with higher numbers of passes. Step 4.1.3. Verify the flow velocity constraint of side II. It is the same as in Step 4.1.1. Step 4.1.4. Verify the pressure drop constraint of side II. It is the same as in Step 4.1.2. Step 4.2. Considerate Yh ¼ 1. It is the same as in Step 4.1. Step 5. Generate the reduced set of configurations (RS). Step 5.1. For Yh ¼ 0, the number of selected passes for sides I and II, for the heat exchanger, are matched, achieving configurations in the format of [NC PI PII Yh]. For each combination of these 4 parameters, 4 different configurations are possible according to the feed arrangement (f parameter). Step 5.2. The same procedure (Step 5.1) is carried out for Yh ¼ 1. Thus it is calculated the remaining configurations that formed the reduced set. Step 6. Calculation of the thermal effectiveness, ECC, considering purely countercurrent flow. If ECC < Emin, than such configurations may be discarded, because due to edge effects and changes of passes, the thermal effectiveness of a PHE is always less than ECC.
ECC ¼
8 > > > < > > > :
1eNTUð1Cr Þ 1Cr eNTUð1Cr Þ
if Cr < 1
NTU NTUþ1
if Cr ¼ 1
(21)
Step 7. Verification of the thermal effectiveness constraint. The selected configurations on Step 5 are simulated in an increasing order of number of channels until the possible optimal set is found. The remaining configurations do not need to be simulated. In the next section, the simulation algorithm for all the reduced set of configurations (RS) is described. Step 8. Achievement of the global optimum. The heat exchanger with the minimal area is determined by comparing all local optima of all types of plates considered. 4.2.1. Simulation algorithm For the development of this algorithm, the boundary conditions equations in the algorithm form given by Gut and Pinto [6] are used. The simulation algorithm is applied for each value of f separately. The algorithm is presented below:
38
F.A.S. Mota et al. / Applied Thermal Engineering 63 (2014) 33e39
Algorithm 2. Step 1. The tri-diagonal matrix coefficients are computed, d(i) ¼ s(i)a(I or II). Step 2. The tri-diagonal matrix is constructed. Step 3. Eigenvalues and Eigenvectors are computed from the tridiagonal matrix. Step 4. The linear system is generated (Eq. (9)). The boundary condition equations in algorithmic form of inlet fluids and change of passes are used. Step 4.1. Generation of the side I portion of Eigenvalues and I I Eigenvectors matrix, A , and of the binary vector, B . Step 4.2. Generation of the side II portion of Eigenvalues and II II Eigenvectors matrix, A , and of the binary vector, B . Step 4.3. Computation of the Eigenvalues and Eigenvectors matrix and computation of the binary vector. Thus the linear system is identified. I
II
I
II
A ¼ A þA B ¼ B þB
Step 5. Determination of the c0i s coefficients by solving the linear system. Step 6. Determination of the output dimensionless temperatures. The boundary conditions of output fluid in algorithmic form is used. Step 7. Computation of the thermal effectiveness. It may be obtained considering any side of the PHE, because the energy conservation is obeyed only if E ¼ EI ¼ EII. Note that for different values of f different values of E can be achieved because NI and NII are dependent on f.
8 I > I ¼ NI max aI ; aII jq q > E > out j in aI NI NII < E ¼ > > > : EII ¼ NIIII max aII ; aIIII jqin qout jII a N N
(22)
Table 3 Thermal effectiveness of RS for b ¼ 60 . #
NC
PI
PII
Yh
ECC
E (f ¼ 1)
E (f ¼ 2)
E (f ¼ 3)
E (f ¼ 4)
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22
40 40 41 41 42 42 43 43 44 44 48 48 48 48 48 48 48 48 48 48 48 48
2 2 3 3 3 3 2 2 2 2 2 2 2 3 3 3 2 2 3 3 4 4
2 2 2 2 3 3 3 3 2 2 2 3 4 2 3 4 2 3 2 3 2 3
0 1 0 1 0 1 0 1 0 1 0 0 0 0 0 0 1 1 1 1 1 1
0.83 0.83 0.86 0.86 0.88 0.88 0.86 0.86 0.84 0.84 0.85 0.87 0.88 0.86 0.89 0.90 0.85 0.86 0.87 0.89 0.88 0.90
0.538 0.538 0.540 0.539 0.534 0.534 0.536 0.538 0.537 0.537 0.537 0.538 0.515 0.535 0.534 0.531 0.537 0.540 0.533 0.534 0.515 0.530
ea ea 0.540 0.539 0.535 0.545 0.536 0.538 0.492 0.492 ea 0.532 0.514 0.540 0.537 0.530 ea 0.535 0.538 0.537 0.514 0.531
0.819 0.819 0.770 0.774 0.774 0.774 0.776 0.772 0.826 0.826 0.833 0.776 0.795 0.772 0.778 0.840 0.833 0.773 0.776 0.778 0.795 0.840
0.691 0.691 0.770 0.774 0.860b 0.860b 0.776 0.772 ea 0.693 ea 0.786 ea 0.782 0.870 0.830 ea 0.782 0.786 0.870 ea 0.830
E (f ¼ 3)
E (f ¼ 4)
a b
Table 4 Thermal effectiveness of RS for b ¼ 45 . # 1 2 3 4 5 6 7 8 a b
5. Results As a case study, an example of optimization was considered for three types of chevron plates with the same dimensions, but with different angles of corrugation pattern, b. Table 2 presents the Table 2 Example data. Plate characteristics
b ¼ [60 45 30 ] F ¼ 1.24
LP ¼ 0.994 m WP ¼ 0.525 m b ¼ 2.85 mm DP ¼ 15 cm
tP ¼ 0.5 mm kP ¼ 17 W/m K
Hot water inlet
Cold water inlet
Tin,hot ¼ 110 C _ M hot ¼ 20118.91 kg/h Rf,hot ¼ 0 cm2 K/W rhot ¼ 965 kg/m3 mhot ¼ 0.315 cP cP,hot ¼ 4.205 J/kg K khot ¼ 0.673 W/m K
Tin,cold ¼ 65 C _ M cold ¼ 23031.38 kg/h Rf,cold ¼ 0 cm2 K/W rcold ¼ 970 kg/m3 mcold ¼ 0.344 cP cP,cold ¼ 4.198 J/kg K kcold ¼ 0.668 W/m K
E >85% 0 DPcold 1 bar
¼ 0.3 m/s vmin hot
vmin ¼ 0.3 m/s cold
NC 40 40 44 44 48 48 48 48
PI 2 2 2 2 2 2 2 3
PII 2 2 2 2 2 3 2 2
Yh 0 1 0 1 0 0 1 1
ECC 0.92 0.92 0.92 0.92 0.93 0.94 0.93 0.94
E (f ¼ 1) 0.539 0.539 0.538 0.538 0.538 0.521 0.538 0.507
E (f ¼ 2) a
e ea ea ea ea 0.507 0.386 0.521
b
0.894 0.894b 0.900 0.900 0.906 0.816 0.906 0.816
0.701 0.701 ea 0.702 ea 0.832 ea 0.832
Proposed analytical solution cannot be applied. Optimal configuration.
example data, which were extracted from Brenner et al. [11], with a modification in the available pressure to 10% of the original value, just to work with a smaller set of configurations. It makes the problem analysis less arduous. The original problem also considered the input and output temperatures as thermal constraints. A minimal thermal effectiveness of 0.85 was considered. The minimal number of channels considered is 40. Usually the plate heat exchangers work with a greater number of channels. With the application of the optimization algorithm until Step 5 for the first element of the vector b, the RS is selected, which is made up of 88 configurations respecting the constraint of the number of plates, pressure drop, and flow velocity (see Table 3). The sixth step of the algorithm discards configurations 1, 2, 9 and 10, because ECC <85%. The seventh step identified two optimal configurations with thermal effectiveness of 86.0% with 42 channels
Table 5 Thermal effectiveness of RS for b ¼ 30 . #
NC
PI
PII
Yh
ECC
E (f ¼ 1)
E (f ¼ 2)
E (f ¼ 3)
E (f ¼ 4)
1 2
48 48
2 2
2 2
0 1
0.94 0.94
0.538 0.538
0.364 ea
0.917b 0.917b
ea ea
Constraints 40 NC 50 0 DPhot 1 bar
Proposed analytical solution cannot be applied. Optimal configuration.
a b
Proposed analytical solution cannot be applied. Optimal configuration.
F.A.S. Mota et al. / Applied Thermal Engineering 63 (2014) 33e39
6. Conclusions
Table 6 Local and global optima.
b
NC
f
PI
PII
Yh
DPhot (bar)
DPcold (bar)
E
60
42 42 40 40 48 48
4 4 3 3 3 3
3 3 2 2 2 2
3 3 2 2 2 2
0 1 0 1 0 1
0.57 0.57 0.44 0.44 0.72 0.72
0.71 0.71 0.54 0.54 0.90 0.90
0.860 0.860 0.894 0.894 0.917 0.917
45 30
39
Values for b ¼ 45 (in bold) correspond to the global optima.
and each stream making 3 passes. The difference between the optimal configurations is the fluids location. One of them considers the hot fluid flowing in side I and the other one considers the cold fluid in the side I. Just for comparison reasons, Table 3 presents all simulations for the reduced set of configurations (RS). Once the optimal set is found the remaining configurations may be discarded. The cells identified by asterisk marks (*) mean that in simulation of such configurations, a pair of null Eigenvalues came up between the Eigenvalues of the tri-diagonal matrix. In this case the proposed analytic solution in Eq. (7a) cannot be applied. An attempt to use Eq. (7b) was done, but the correct solution was not obtained as well. The reduced set of configurations for the next element of the vector b is presented in Table 4. It were selected 32 configurations, amount smaller than that one found for the first element of b. However, it was expected because the smaller the angle the greater the pressure drops. The sixth step of the algorithm does not discard any configuration, as may be observed in Table 4. In this case ECC >85% for all RS configurations. When Step 6 does not eliminate any element of RS, it is because the value for the constraint of minimal number of channels was not adequate. With a smaller number of channels, the thermal constraint could already be obeyed. In this particular case, for example, it was verified that 28 channels with a 2-2 arrangement would be sufficient to meet the thermal constraint. By applying the thermal constraint (Step 7), the local optima set was identified in the two first configurations, with 40 channels, symmetrical arrangements of 2 passes. For the last element of the vector b, the selected RS are also the optima configurations (see Table 5). For this case, however, no simulation would need to be carried out, because no configuration could be the global optimum since the RS configurations have more channels than the ones previously computed. Finally, by comparing all local optima, the global optimum set are identified for b ¼ 45 (Table 6). It is interesting to note that the global optimum set was the one that used the least of the available pressure drop. In order to obtain a better performance for the available pressure drop and, consequently, configurations with fewer plates, plates with different corrugation patterns in the PHE would be a possibility. The configurations with f ¼ 3 will have always equal or better thermal performance because with f ¼ 3 and an even number of passes the heat transfer behavior in the heat exchanger is similar to the countercurrent, that is the most thermal efficient system.
This paper presented an optimization algorithm for the design of plate heat exchangers. It was based on the model proposed by Gut and Pinto [3]. The difference between the present paper and the work of Gut and Pinto [3] is that for a given process and type of plate the proposed algorithm is able to find the configuration of the lowest number of thermal plates (the optimal configuration for a given PHE). The types of plates are considered as optimization variables. The inclusion of this possibility into a systematic procedure is a great novelty in the paper. For each type of plate, a local optima set may be found (the configurations presenting the lowest number of plates). By comparing all the local optima sets, the global optimum set may be found. The system of differential equations was solved analytically, avoiding numerical methods or other approximations. Also, it is important to highlight that no additional software needs to be used, like gPROMS or other ones. This is the great contribution of the present paper when compared with the work of Gut and Pinto [3]. As a case study, an example of the literature was used to show the applicability of the developed algorithm. The example of optimization has considered three types of chevron plates with the same dimensions, but with different angles of corrugation pattern, b. By the optimization example presented, although some configurations cannot be solved by analytical solution proposed, the developed algorithm had great success. The complex Eigenvalues compromised the simulation of some configurations where both streams make a number of pass multiple of two. However, the configuration at f ¼ 3 generally presents best thermal performance for these cases. So, the fact that the proposed analytical solution cannot be applied for such configurations is not a limiting factor for the algorithm use. References [1] L. Wang, B. Sundén, R.M. Manglik, Plate Heat Exchangers: Design, Applications and Performance, WIT Press, Ashurst Lodge, 2007, pp. 27e39. [2] R.K. Shah, D.P. Sekulic, Fundamentals of Heat Exchanger Design, John Wiley & Sons, Inc., New Jersey, 2003, pp. 22e29. [3] J.A.W. Gut, J.M. Pinto, Optimal configuration design for plate heat exchangers, Int. J. Heat Mass Transfer 47 (2004) 4833e4848. [4] H. Najafi, B. Najafi, Multi-objective optimization of a plate and frame heat exchanger via genetic algorithm, Heat Mass Transfer 46 (2010) 639e647. [5] O.P. Arsenyeva, L.L. Tovazhnyansky, P.O. Kapustenko, G.L. Khavin, Optimal design of plate-and-frame heat exchangers for efficient heat recovery in process industries, Energy (2011) 1e11. [6] J.A.W. Gut, J.M. Pinto, Modeling of plate heat exchangers with generalized configurations, Int. J. Heat Mass Transfer 46 (14) (2003) 2571e2585. [7] E.A.D. Saunders, Heat Exchangers: Selection, Design & Construction, Longman S & T, Harlow (UK), 1988. [8] T. Zaleski, A.B. Jarzebski, Remarks on some properties of the equation of heat transfer in multichannel exchangers, Int. J. Heat Mass Transfer 16 (8) (1973) 1527e1530. [9] T. Zaleski, A general mathematical model of parallel-flow, multichannel heat exchangers and analysis of its properties, Chem. Eng. Sci. 39 (7/8) (1984) 1251e1260. [10] S. Kakaç, H. Liu, Heat Exchanger: Selection, Rating and Thermal Design, second ed., CRC Press, Boca Raton, 2002, pp. 373e412. [11] K. Brenner, G. Conrads, J. Dorn, G. Fritz, Z.G. Guo, N. Jahnke, W. Tinz, E. Weib, Plate Heat Exchangers: Installation and Use of Multi-purpose Heat Exchange Systems, Bosch-Druck GmbH, Landshut, 1998.