Dynamic Optimization of the Flow Rate Distribution in Heat Exchanger

Article pubs.acs.org/IECR

Dynamic Optimization of the Flow Rate Distribution in Heat Exchanger Networks for Fouling Mitigation Bruna C. G. Assis,‡ Julia C. Lemos,† Fábio S. Liporace,‡ Sérgio G. Oliveira,‡ Eduardo M. Queiroz,§ Fernando L. P. Pessoa,§ and André L. H. Costa*,† †

Instituto de Química, Rio de Janeiro State University (UERJ), Rua São Francisco Xavier, 524, Maracanã, Rio de Janeiro, RJ, CEP 20550-900, Brazil ‡ Petrobras Research & Development Center (CENPES/PETROBRAS), Cidade Universitária, Avenida Horácio Macedo, 950, Rio de Janeiro, RJ CEP 21949-900, Brazil § Escola de Química, Federal University of Rio de Janeiro (UFRJ), Avenida Athos da Silveira Ramos, 149, Ilha do Fundão, Rio de Janeiro, RJ CEP 21949-900, Brazil ABSTRACT: Heat exchanger networks are structures composed of a set of heat exchangers interconnected in order to reduce utilities consumption. During the network operation, heat exchangers may present a decrease of their thermal effectiveness caused by fouling, which corresponds to the undesirable accumulation of deposits over their thermal surface. In this context, this paper presents a proposal to increase the energy recovered in heat exchanger networks affected by fouling through the optimization of the distribution of the flow rates of the process streams. The problem corresponds to a dynamic optimization problem, because the flow rate optimization affects the surface temperature and velocity, which modifies the fouling rate, thus demanding the simultaneous analysis of the entire time horizon. The objective function is represented by the integral of the utility consumption during the operational time horizon. The main constraints include mass and energy balances, heat exchangers equations (P-NTU method), and fouling rate modeling. The mathematical structure of the problem corresponds to a nonlinear optimization. The utilization of the optimization scheme is illustrated by the analysis of two examples of heat exchanger networks.

1. INTRODUCTION Heat exchanger networks (HENs) are structures composed of a set of interconnected heat exchangers which promote the heat transfer between process streams, therefore allowing a reduction in the utility consumption in process plants. Because of the increase of energy costs, the design of these systems has attracted considerable attention since the 1980s of the last century.1 However, the energy efficiency increase obtained through a HEN may be severely disturbed by fouling. Heat exchanger fouling is the undesirable accumulation of deposits over the heat exchange surface. This phenomenon introduces extra resistances in the thermal circuit, bringing a reduction of the overall heat transfer coefficient. As a consequence, there is a diminution of the heat exchanger effectiveness. The literature presents several approaches for fouling mitigation based on the utilization of computational resources. Considering the operation of a HEN, a subject widely investigated corresponds to the optimization of the cleaning schedule of the heat exchangers. This mitigation technique was explored using mixed-integer nonlinear programming (MINLP),2−4 mixed-integer linear programming (MILP),5−7 and stochastic methods.4,8 More recently, additional aspects were explored together with the cleaning schedule: manipulation of by-passes,8 desalter control,9 and hydraulic behavior.10 Fouling mitigation can also be considered during the design of individual heat exchangers,11,12 HEN synthesis/retrofitting,13,14 or the utilization of heat transfer enhancement devices in the HEN synthesis.15,16 © 2015 American Chemical Society

Another approach for fouling mitigation involves the optimization of the distribution of the flow rates in HENs. This approach was proposed for fouling management in crude preheat trains by Oliveira Filho et al.,17 based on a steady-state modeling of the system using a given set of fouling resistances. The extension of this formulation also considering the hydraulic behavior of HENs through a constrained mathematical programming problem was explored in Assis et al.18 The flow split optimization associated with the cleaning of parallel heat exchangers was studied in Ishiyama et al.19 A previous work involving the operational optimization of the flow rate distribution in HENs, not linked to fouling issues, can be found in Lid et al.20 However, the redistribution of the flow rates along the heat exchangers in a network may have a strong effect on the fouling rate. In these cases, the application of the optimal solution found for a given instant may bring a potential loss of the energy recovery in the future, due to the impact of the modification of the fouling rates; that is, the set of adopted flow rates can determine an increase of fouling in certain heat exchangers which may cause an opposite effect in the energy recovery when considering the entire operational horizon. An additional discussion about the impact of flow rate distribution Received: Revised: Accepted: Published: 6497

May May May May

30, 2014 4, 2015 12, 2015 12, 2015 DOI: 10.1021/acs.iecr.5b00453 Ind. Eng. Chem. Res. 2015, 54, 6497−6507

Article

Industrial & Engineering Chemistry Research on fouling, including thermohydraulic effects, can be found in Ishiyama et al.21 The focus of the current paper is to propose a mathematical programming approach for the optimization of the flow rate distribution considering the dynamic behavior of the HEN associated with the modifications of the fouling rates. The resultant formulation corresponds to a nonlinear optimization problem aiming to minimize the total utility consumption cost. The constraints are composed of mass and energy balances together with heat exchanger equations (P-NTU method) associated with fouling rate models. The rest of this paper is organized as follows. Section 2 presents the formulation of the optimization problem, section 3 discusses the initialization of the optimization procedure in order to avoid nonconvergence problems, section 4 illustrates the application of the proposed optimization scheme through two HEN examples, and section 5 reports the conclusions.

∑

mk , τ −

k ∈ Stin

∑

mk , τ = 0

t ∈ (MX ∪ SP), τ ∈ TI

k ∈ Stout

(1)

∑

mk , τ −

k ∈ Stin

∑

mk , τ + nt , τ = 0

k ∈ Stout

t ∈ (PS ∪ PD), τ ∈ TI

∑

∑

mk , τ −

k ∈ (Stin ∩ HSTR )

(2)

mk , τ = 0

k ∈ (Stout ∩ HSTR )

t ∈ HE , τ ∈ TI

∑

mk , τ −

k ∈ (Stin ∩ CSTR )

(3)

∑

mk , τ = 0

k ∈ (Stout ∩ CSTR )


nt , τ − ntspe ,τ = 0

2. OPTIMIZATION PROBLEM FORMULATION 2.1. Network Description. The HEN organization is represented by a digraph.22 The digraph vertices are identified by the index t ∈ VET, representing the network elements: heat exchangers (subset HE), mixers (subset MX), splitters (subset SP), supply units (subset PS,) and demand units (subset PD). The inlet and outlet streams between the network and external equipment are located at the vertices of the supply and demand units. The set of vertices are interlinked by the digraph edges, identified by the index k ∈ STR, which represents the material streams. These streams are divided into hot streams (subset HSTR) and cold streams (subset CSTR). 2.2. Optimization Variables. The problem variables are the network mass flow rates and temperatures, the heat loads and heat transfer coefficients of the heat exchangers, the PNTU dimensionless groups and the variables related to the fouling rate modeling. All sets of variables are established for each time instant τ ∈ TI of the simulation horizon. The mass flow rates and temperatures of the network streams are represented by mk,τ and Tk,τ. The external mass flow rates and temperatures in the demand/supply units are represented by nt,τ and Vt,τ. The heat exchanger variables are the heat load, Qt,τ, the overall and convective heat transfer coefficients, Ut,τ, htube,t,τ and hshell,t,τ, and the corresponding P-NTU dimensionless groups, Pt,τ, NTUt,τ and CRt,τ. The optimization variables associated with the fouling rate model are the fouling rate, dRfdtt,τ, the corresponding values of the fouling rate at the heat exchanger ends, dRfdtt1,τ and dRfdtt2,τ, the fluid flow velocity, vt,τ, the Fanning friction factor, f t,τ, and the surface and film temperatures at the exchanger ends, Tfilm1,t,τ, Tfilm2,t,τ, Ts1,t,τ, Ts2,t,τ. 2.3. Optimization Constraints. The main optimization constraints are composed of mass and energy balances, heat exchanger equations and fouling rate model equations. Because fouling accumulation is much slower than the dynamic behavior of the heat exchangers, the description of the HEN behavior during the time horizon of the optimization is based on a pseudostationary model, in which the evolution of the process variables are described by a sequence of network steady-states. The values of the fouling resistances are updated in each time instant according to a proper fouling rate model. Mass Balances. The mass conservation along the network is represented by

(4)

t ∈ PS , τ ∈ TI

(5)

in

out

where St is the set of the edges k direct to the vertex t, St is the set of edges k direct from the vertex t, and nspe t,τ are the problem specifications of the network inlet flow rates. Energy Balances. The energy conservation along the network is represented by

∑ k ∈ Stin

mk , τ Cpk Tk , τ −

∑ k ∈ Stout

mk , τ Cpk Tk , τ = 0

t ∈ (MX ∪ SP), τ ∈ TI

∑ k ∈ Stin


∑ k ∈ Stout

(6)

mk , τ Cpk Tk , τ + nt , τ CPV t t ,τ = 0

t ∈ (PS ∪ PD), τ ∈ TI

∑ k ∈ (Stin ∩ HSTR )


=0


∑


k ∈ (Stin ∩ CSTR )

=0

(7)

∑ k ∈ (Stout ∩ HSTR )

mk , τ Cpk Tk , τ − Q t , τ (8)

∑ k ∈ (Stout ∩ CSTR )

mk , τ Cpk Tk , τ + Q t , τ


Vt , τ − V tspe ,τ = 0

(9)

t ∈ PS , τ ∈ TI

(10)

where Cpk is the heat capacity of the stream of the edge k, CPt is the heat capacity of the inlet/outlet stream through the supply/ demand vertex t, and Vspe t,τ is the specification of the temperature of the network inlet stream at the supply unit t. This set of constraints must be complemented by an equation which imposes that the splitter outlet temperatures are equal: Tk , τ − Tk ′ , τ = 0

t ∈ SP , k ∈ Stin , k′ ∈ Stout , τ ∈ TI (11)

Heat Exchanger Equations. The heat transfer in heat exchangers is described using the P-NTU method.23 This method is similar to the ε-NTU method, but always directed to one of the streams (without depending on, as in the ε-NTU method, which stream is the minimum and maximum fluids). In this paper, the P-NTU equations are described in relation to the hot stream. 6498

DOI: 10.1021/acs.iecr.5b00453 Ind. Eng. Chem. Res. 2015, 54, 6497−6507

Article

Industrial & Engineering Chemistry Research ⎡

The heat exchanger equations of the P-NTU method can be written using the optimization variables, therefore yielding the following set of constraints:

htube, t , τ −

cm base ⎢ t k , τ htube, t ⎢⎣

+ (1 − ct )mk ′ , τ ⎤ ⎥ base ⎥⎦ mtube, t

0.8

=0

t ∈ HE , k ∈ (CSTR ∩ Stin), k′ ∈ (HSTR ∩ Stin), τ

Pt , τTk , τ + (1 − Pt , τ )Tk ′ , τ − Tk ″ , τ = 0 k″ ∈ (HSTR ∩ Stout ), τ ∈ TI

(18)

∈ TI

t ∈ HE , k ∈ (CSTR ∩ Stin), k′ ∈ (HSTR ∩ Stin), (12)

hshell, t , τ −

NTUt , τmk , τ Cpk − Ut , τA t = 0 t ∈ HE , k ∈ (HSTR ∩ Stin), τ ∈ TI

⎡ (1 − c )m + c m ⎤0.6 t k ,τ t k′,τ ⎥ =0 base ⎢⎣ ⎥⎦ mshell, t

base ⎢ hshell, t

t ∈ HE , k ∈ (CSTR ∩ Stin), k′ ∈ (HSTR ∩ Stin),

(13)

(19)

τ ∈ TI

CR tmk Cpk − mk ′Cpk ′ = 0

where the parameter ct is equal to 1 if the cold stream flows in the tube-side of the heat exchanger t, otherwise it is equal to 0. This alternative allows the utilization of base case values for the convective coefficients obtained using more complex and accurate methods, as those employed in heat exchanger design softwares (e.g., HTRI), but considering their variations due to the flow rate optimization. Fouling Rate Modeling. The application of the proposed optimization approach demands the availability of a fouling model able to describe the impact of the variation of the operating variables on the fouling rate. The literature contains several alternatives of semiempirical models with this feature.26 Because of the complexity of the fouling phenomenon, each model is usually focused on a particular application. In this context, the current paper explores the utilization of fouling rate models for fouling prediction in crude preheat trains.27 The fouling rate in these systems can be described using the concept of fouling threshold. In these models, the fouling rate is evaluated by a difference between a deposition rate and a suppression rate. If the suppression rate is higher than the deposition rate, there is no fouling, therefore determining a fouling threshold.28 Among the several variant models based on this concept, this paper employed the modified Ebert−Panchal model. Excluding deposit removal effects,28 the expression of this model becomes21,29

t ∈ HE , k ∈ (CSTR ∩ Stin), k′ ∈ (HSTR ∩ Stin), (14)

τ ∈ TI

where Pt,τ is the temperature effectiveness, NTUt,τ is the number of transfer units, CRt,τ is the ratio between the heat capacity flow rates, Ut,τ is the overall heat transfer coefficient, and At is the heat transfer area. The heat exchangers considered in this paper presents countercurrent flow or a multipass configuration with an even number of passes in the tube-side. The corresponding P-NTU expressions are, respectively: Pt , τ −

1 − exp[−NTUt , τ(1 − CR t , τ )] 1 − CR t , τ exp[− NTUt , τ(1 − CR t , τ )]

=0

t ∈ SCC ⊂ HE , τ ∈ TI

(15)

⎧ Pt , τ − 2⎨1 + CR t , τ + (1 + CR t2, τ )0.5 ⎩ ⎪

⎪

−1 1 + exp[−NTUt , τ(1 + CR t2, τ )0.5 ] ⎫ ⎪ ⎬ =0 × ⎪ 1 − exp[−NTUt , τ(1 + CR t2, τ )0.5 ] ⎭

t ∈ SMP ⊂ HE , τ ∈ TI

(16)

⎛ ⎞ ⎛ −Ea ⎞ = max⎜⎜0, αRe−0.66Pr 0.33 exp⎜ ⎟ − γτw⎟⎟ dt ⎝ RTfilm ⎠ ⎝ ⎠

dR f

where SCC is the subset of countercurrent heat exchangers and SMP is the subset of multipass heat exchangers (HE = SCC ∪ SMP). The overall heat transfer coefficient can be calculated from the individual convective heat transfer coefficients:24

where Re and Pr are Reynolds and Prandtl dimensionless groups, R is the universal gas constant, τw is the shear stress, Tfilm is the film temperature and α, γ, and Ea (activation energy) are empirical model parameters. In this case, it is assumed that fouling only occurs due to the crude oil, which flows in the tube-side. The insertion of this fouling rate model in the optimization formulation is represented by the following set of constraints:

D log(De , t /Di , t ) 1 1 ⎛ De , t ⎞ 1 ⎜⎜ ⎟⎟ − e , t − − Ut , τ htube , t , τ ⎝ Di , t ⎠ 2k w , t hshell , t , τ ⎛ De , t ⎞ ⎟⎟ = 0 − R f , t , τ ⎜⎜ ⎝ Di , t ⎠

(20)

t ∈ HE , τ ∈ TI (17)

dRf dtt , τ = 0.5(dRf dt1, t , τ + dRf dt 2, t , τ )

where htube,t,τ and hshell,t,τ are the convective heat transfer coefficients in the tube-side and shell-side, Di,t and De,t are the inner and outer tube diameters, Rf,t,τ is the fouling resistance and kw,t is the thermal conductivity of the tube wall. In eq 17, it is assumed that the fouling is limited to the tube-side, coherent with the fouling rate model discussed in the next subsection. The convective heat transfer coefficients can be evaluated from base case values associated with a power law correction factor applied to flow rate variations:25

t ∈ HE , τ ∈ TI (21)

−1 dRf dt1, t , τ ≥ α(Di , t vt , τρt /μt )−0.66 Prt0.33 exp( −EaR−1Tfilm1, t , τ ))

− γ(ft , τ ρt vt2, τ /2)

dRf dt1, t , τ ≥ 0 6499



(22) (23)


Article

Industrial & Engineering Chemistry Research

Finally, it is necessary to include in the set of the optimization constraints the integration of the fouling rate model:

−1 dRf dt 2, t , τ ≥ α(Di , t vt , τρt /μt )−0.66 Prt0.33 exp( −EaR−1Tfilm2, t ,τ)

− γ(ft , τ ρt vt2, τ /2)

dRf dt 2, t , τ ≥ 0


(24)


R f , t , τ + 1 = R f , t , τ + 0.5(dRf dtt , τ + 1 + dRf dtt , τ )Δt

(25)

where ρt and μt are the fluid density and viscosity, vt,τ is the fluid velocity in the tube-side, f t,τ is the Fanning friction factor, and dRfdtt,τ is the fouling rate. The indices 1 and 2 correspond to the extremities of the heat exchanger. In fact, the fouling rate varies along the heat transfer surface and the rigorous evaluation of lumped values would demand a numerical integration. To avoid such complexity, an average of the fouling rates evaluated at each extremity of the heat exchanger for an equivalent countercurrent configuration was adopted. Another alternative for evaluation of integral fouling rate values can be found in Ishiyama et al.21 It should be noted that eqs 21−25 compose a relaxation of the original eq 20. However, because the optimization drives the solution toward a reduction of the energy consumption, the optimal solution will always present an activation of one of the constraints of the pairs eqs 22 and 23, and 24 and 25; that is, the converged fouling rate at the solution will not be higher than both lower bounds, thus obeying the original equation. The reason for this representation is to overcome the nondifferentiability present in eq 20 due to the threshold behavior. The fouling rate modeling is complemented by an additional set of constraints representing the evaluation of the Fanning friction factor, fluid flow velocity, film, and surface temperatures: ft , τ = 0.0035 + 0.264(Di , t vt , τρt /μt )−0.42

t ∈ HE , τ ∈ TI *

R f ,t ,τ =

R spe f ,t

(32)

t ∈ HE, τ = 0

(33)

where TI* is equivalent to the set TI without the last element, Δt is the time increment between two instants and Rspe f,t is the initial conditions (τ = 0) of the fouling resistance. Bounds on Variables. Constraints associated with flow velocities, heaters and coolers capacity limits, and pumparounds operational ranges are represented by bounds on temperatures and flow rates: mkLB ≤ mk , τ ≤ mkUB TkLB ≤ Tk , τ ≤ TkUB

k ∈ SPEm , τ ∈ TI k ∈ SPET , τ ∈ TI

(35)

where the superscripts LB and UB represent the lower and upper bounds, and SPEm and SPET are the sets of bounds. These constraints aims to guarantee the operational feasibility of the solution during the entire time horizon. Additional constraints may also be imposed to represent non-negative bounds according to the nature of the corresponding variables. Auxiliary Constraints. Considering that fouling is a slow process (excluding abnormal conditions), the convergence can be facilitated through additional constraints, which limit the variation of the flow rate between consecutive instants: mk , τ (1 − b) ≤ mk , τ + 1 ≤ mk , τ (1 + b)


k ∈ STR , τ ∈ TI * (36)

(26)

where b is assumed equal to 0.05. 2.4. Objective Function. The objective function corresponds to the utility consumption cost during the operational time horizon:

vt , τ = (mk , τ /ρt )/[(πDi2, t /4)Ntp , t ] t ∈ HE , k ∈ (CSTR ∩ Stin), τ ∈ TI

(27)

Tfilm1, t , τ = 0.5(Ts1, t , τ + Tk , τ ) t ∈ HE , k ∈ (CSTR ∩

min fobj =

Stin),

τ ∈ TI

1

=

h tube, t , τ Di , t

(29)

Tk ′ , τ − Ts1, t , τ 1 hshell, t , τ De , t

+

t ∈ HE , k ∈ (CSTR ∩

log(De , t / Di , t ) 2k w

Stin),

+

R f ,t ,τ Di , t

k′ ∈ (HSTR ∩ Stout),

Ts2, t , τ − Tk , τ 1 h tube, t , τ Di , t

sτ = (1 + i)−npτ (30)

τ ∈ TI =

hshell, t , τ De , t

+

log(De , t / Di , t ) 2k w

+

R f ,t ,τ

(38)

3. INITIALIZATION The resultant problem is a nonlinear programming problem (NLP) composed of the objective function in eq 37 subject to the constraints in eqs 1−19 and 21−36. Aiming to reduce nonconvergence problems, an automatic procedure for generation of initial estimates is employed. This procedure is illustrated in Figure 1. This procedure is based on the solution of two linear systems of equations. The initialization of the mass balances involves

Di , t

t ∈ HE , k ∈ (CSTR ∩ Stout), k′ ∈ (HSTR ∩ Stin), τ ∈ TI

(37)

where i is the interest rate and npτ is the number of past periods associated with time instant τ.

Tk ′ , τ − Ts2, t , τ 1

pτ sτ COP , tnt CPt(Vref, t − Vt , τ )

where pτ is the weights of the numerical integration procedure used (Simpson rule), sτ is the present worth factors, COP,t is the utility costs associated with the network outlet streams (COP,t > 0, for a hot utility, and COP,t < 0, for a cold utility; coherent with the signal of the difference Vref,t − Vt,τ in each case, yielding always a positive objective function), and Vref,t is the reference temperature for the evaluation of the utility consumption. The present worth factors are given by30

Tfilm2, t , τ = 0.5(Ts2, t , τ + Tk , τ )

Ts1, t , τ − Tk , τ

∑ ∑ t ∈ PD τ ∈ TI

(28)

t ∈ HE , k ∈ (CSTR ∩ Stout), τ ∈ TI

(34)

(31)

where Ntp,t is the number of tubes per pass, Tfilm1,t,τ and Tfilm2,t,τ are the film temperatures at the exchanger ends, and Ts1,t,τ and Ts2,t,τ are the deposit surface temperatures at the heat exchanger ends, assuming a countercurrent equivalent configuration. 6500


Article


Figure 2. Heat exchanger network of the Example 1: continuous lines, cold streams; dashed lines, hot streams.

Figure 1. Initial estimate generation procedure.

structure commonly found in the practice of engineering. This layout can be observed for two individual units or even for two or more branches of identical heat exchangers in series.2−4,10,17 The installation of parallel units is sometimes necessary due to the large magnitude of the thermal service, which may be infeasible to be carried out in a single train of heat exchangers. Another reason to employ parallel heat exchangers is to provide the process with a more flexible flowsheet. In this case, it is possible to clean one of the heat exchangers without interrupting the operation of the whole process. The mass flow rates, temperatures, and heat capacities of the inlet streams are shown in Table 1 and the physical properties of the cold inlet stream are shown in Table 2. The heat capacities and further physical properties of the streams are assumed constant along the HEN.

the solution of eqs 1−5, complemented by specifications for the stream splits. Using the evaluated values of mass flow rates, the heat exchanger variables are then calculated based on the initial fouling conditions. The initialization of the energy balances is conducted by the solution of the linear system formed by eqs 6−12, using the variables already calculated as fixed parameters. Finally, the values of fouling modeling variables are calculated using the flow rates and temperatures determined in the solution of the linear systems. It is important to observe that this structure of the initialization procedure could be proposed due to the utilization of the P-NTU method for describing the heat exchanger behavior. In the P-NTU method, the resultant equations are linear in relation to the temperatures for a given set of flow rates.

Table 1. Inlet Streams Specifications

4. RESULTS The potential of the proposed optimization approach is illustrated using two examples. The first example is a hypothetical HEN composed of two heat exchangers aligned in parallel. This kind of example allows a clear assessment of the features of the optimal solution when compared to the traditional solution that keeps the mass flow rates equally divided into parallel equipment. The second example is a HEN composed of six heat exchangers distributed along two branches. The aim of this analysis is to illustrate the utilization of the proposed approach to a larger HEN with different heat exchangers. In both examples, the GAMS software was used to implement the proposed optimization model, where the NLP problems were solved using the CONOPT solver, and the linear problems used to generate initial estimates were solved using the CPLEX solver. The hardware employed in the resolution of the examples was a computer with a processor Intel Core i7 3.40 GHz and 12 GB of RAM memory. 4.1. Example 1. The process flowsheet of the investigated network on Example 1 is shown in Figure 2. This example investigates the operation of a small network composed of two identical shell-and-tube heat exchangers organized in parallel. Identical heat exchangers designed to operate in parallel with a uniform flow rate distribution is a

supply node

mass flow rate (kg/s)

temperature (°C)

specific heat (J/(kg·K))

1 2

310.3 175.8

157.1 213.8

2378.5 2476.0

Table 2. Physical Properties of the Cold Stream density (kg/m3)

viscosity (Pa·s)

thermal conductivity (W/(m·K))

818.8

1.27 × 10−3

0.10

The data of the heat exchangers are displayed in Table 3, where the base case is associated with a uniform flow rate distribution. The values of the fouling rate model parameters are the same in both heat exchangers: α = 127.7 m2K/J, Ea = 76 kJ/mol, and γ = 3.44 × 10−15 m2K/(J·Pa). The initial fouling resistances in the heat exchangers 5 and 6 are 9.89 × 10−4 m2K/W and 0 m2K/W, respectively. This situation typically originates from a difference in the cleaning record of the heat exchangers.21 The dynamic optimization was conducted considering an operation period of 1 year, for which the time was discretized in 1 week periods. The base case was used to compare the results, which considers the split fractions equal to 0.5 in the two 6501


Article

Industrial & Engineering Chemistry Research Table 3. Heat Exchanger Data parameter

heat exchangers 5 and 6

tube-side flow heat transfer area (m2) number of shell and tube passes total number of tubes inner and outer diameters (mm) tube wall thermal conductivity (W/(m·K)) tube-side heat transfer coefficient at the base case (W/(m2·K)) shell-side heat transfer coefficient at the base case (W/(m2·K))

cold stream 546 1−2 1520 14.85/19.05 55 1365 1540

Figure 3. Mass flow rates of the cold streams for the base and optimized cases.

Figure 4. Heat load recovery.

The resultant problem is composed of 4348 variables and 4975 constraints and was solved in about 1.5 s. The distribution of the crude flow rate in the heat exchangers during the operational time for the base and optimized cases are shown in Figure 3. According to Figure 3, the mass flow rates in the heat exchangers vary during the time in the optimized case, and it can be explained by the relative difference of the fouling between these heat exchangers. The initially clean heat exchanger receives the major part of the flow during the entire

splitters during the entire operation time. Operational bounds were not imposed on the variables. The objective function considers only the hot utility consumption, and it is equivalent to the minimization of total energy consumption due to fouling during the investigated time horizon (COP,t = 1 and i = 0). The reference temperature adopted is 177.9 °C, corresponding to the network outlet temperature at the clean condition and uniform flow distribution. 6502


Article


Figure 5. Fouling rates.

Figure 6. Fouling resistance.

Despite the flow rates in the heat exchangers in the base case are the same, the initially dirty heat exchanger presents a smaller fouling rate. This difference occurs because the presence of the deposits diminishes the deposit surface temperature which decreases the deposition rate of fouling (the positive term in the fouling rate equation). In the optimized case there is another trend, the initially dirty heat exchanger presents a higher fouling rate. This pattern occurs because of the lower value of the flow rate in this exchanger, which severely decreases the suppression rate of fouling (the negative term of the fouling rate equation). In both cases, the fouling rates in all heat exchangers decrease during the operational time. This pattern occurs due to the increase of the fouling resistances, which causes a continuous decrease in the deposit surface temperature and, consequently, a decrease in the fouling rate. Figure 6 illustrates the behavior of fouling resistance in the heat exchangers for the base and optimized cases. The fouling resistances profiles in Figure 6 are coherent with the corresponding patterns of the fouling rates shown in Figure 5. The gain achieved with the dynamic optimization is shown in Table 4, where it is possible to observe that, during the 52

operation time, almost 77% of the split fraction at the beginning that reduces to 70% at the end. Larger time horizons bring additional asymptotically reductions in the split difference at the end of the period. This pattern occurs because the difference between the fouling resistances becomes less important as both resistance values increase with time due to the deposit growth. The heat load recovery in the heat exchangers during the considered period for the base and the optimized cases are shown in Figure 4. The difference between the heat loads of the heat exchangers at the base case is a consequence of the different initial fouling resistance values. In the optimized case, this difference is much larger because of the higher flow rate directed to the initially clean exchanger. In all cases, the heat load recovery in the heat exchangers decreases with time because of the accumulation of deposits. The total heat load recovery illustrates the success of the optimization approach, such that the total heat load recovered in the optimized case is always larger than the corresponding value of the base case. The fouling rates in the heat exchangers for the base and optimized cases can be observed in Figure 5. 6503


Article

Industrial & Engineering Chemistry Research weeks analyzed, the optimization reaches a reduction of 7% on the energy consumption due to fouling.

Table 5. Inlet Streams Specifications

Table 4. Total Energy Consumption case

total energy required (PJ)

average energy consumption rate (kW)

base optimized

0.162 0.151

5149 4801

supply node

mass flow rate (kg/s)

temperature (°C)

specific heat (J/(kg K))

1 2 3 4

200.0 50.0 50.0 50.0

30.0 130.0 151.1 172.3

2100.0 2500.0 2500.0 2500.0

The objective function in this example is also equivalent to the minimization of hot utility consumption due to fouling, where the reference temperature is 99.6 °C, defined from the clean network condition associated with the base case flow rate distribution. The optimization problem contains 11 556 variables and 13 439 constraints and was solved in about 7.5 s. The distribution of the cold stream flow rate between the branches during the operational time for the base and optimized cases is shown in Figure 8. Because the network branches are different, the result of the optimization of the flow rates depends on the fouling rate behavior and the thermal surfaces of the heat exchangers of each branch. The optimal profile of the flow rates is considerably different from the design flow rates. The design split is 60/40 and the optimal split varies from 58/42 to 73/27 during the time horizon. The heat load recovery in the branches for the base and the optimized cases is shown in Figure 9. The ratio between the heat load recoveries in each branch is almost constant in the base case. The larger branch is associated with a higher heat load and both heat loads decrease with time because of the deposit accumulation. In the optimization result, the increase of the flow rate in branch A until week 19 determines an increase of the corresponding heat load, overcoming the deposit accumulation effect. Simultaneously, the flow rate decrease in branch B determines a fast reduction of its heat load. After week 19, there

4.2. Example 2. The process flowsheet of the investigated network in Example 2 is shown in Figure 7. The objective of this example is to demonstrate the utilization of the proposed optimization scheme for more complex networks, composed of multiple heat exchangers with different sizes. This example investigates the operation of a HEN composed of six different shell-and-tube heat exchangers distributed along two parallel branches (A and B). The heat exchangers of branch A are larger than the heat exchangers of branch B. This size difference can be found, for example, in revamps where the additional branch inserted into the HEN is associated with a smaller capacity than the original one. The mass flow rates, temperatures, and heat capacities of the inlet streams are shown in Table 5. The physical properties of the cold inlet stream are the same of Example 1. The heat capacities and further physical properties of the streams are assumed constant along the HEN. The data of the heat exchangers are displayed in Tables 6 and 7. The base case is associated with a flow rate distribution of 60% to branch A and 40% to branch B. The values of the parameters of the fouling rate model are the same in all heat exchangers: α = 127.7 m2K/J, Ea = 60 kJ/mol, and γ = 3.44 × 10−15 m2K/(J Pa). All heat exchangers are clean at the beginning of the horizon. The dynamic optimization was conducted considering a 1 year horizon, discretized in 1 week intervals. Operational bounds were not included in the optimization.

Figure 7. Heat exchanger network of the Example 2: continuous lines, cold streams; dashed lines, hot streams. 6504


Article

Industrial & Engineering Chemistry Research Table 6. Heat Exchanger Data of Branch A parameter

heat exchanger 9

heat exchanger 10

heat exchanger 11

tube-side flow heat transfer area (m2) number of shell and tube passes total number of tubes inner and outer diameters (mm) tube wall thermal conductivity (W/(m ·K)) tube-side heat transfer coefficient at the base case (W/(m2·K)) shell-side heat transfer coefficient at the base case (W/(m2·K))

cold stream 500 1−2 1500 14.85/19.05 55 1039 808

cold stream 535 1−2 1605 14.85/19.05 55 973 803

cold stream 520 1−2 1560 14.85/19.05 55 1000 805

considering the energy consumption during the investigated time horizon. The isolated analysis of each instant would bring suboptimal results for the entire interval. A similar problem where the values of the operational variables are related to the degradation of the process performance can be found in the catalyst deactivation of chemical reactors.31,32 The optimization of the heat exchanger network flow rates corresponds to a nonlinear programming problem, which seeks to minimize the utility costs, subjected to mass and energy balances, heat exchanger equations and fouling rate equations. The problem is formulated using the Ebert-Panchal fouling rate model, but the structure of the optimization formulation is flexible and can accommodate other fouling rate models. The proposed optimization approach can be applied considering the presence of degrees of freedom in the mass balance associated with the existence of branches in parallel, a structure commonly found in heat exchanger networks. The solution of the problem indicates the flow rates and temperatures along the network for each time instant of the investigated horizon. Nonconvergence problems are avoided through an initialization procedure. Numerical results in two examples illustrated that the optimization procedure can reduce energy costs when compared to a uniform flow rate distribution, thus mitigating the fouling impact. This result was observed with identical parallel equipment and nonidentical parallel equipment. In the latter problem, the optimal profile is also affected by the difference in the heat exchangers surface area. Since the application of the proposed approach only depends on the presence of flow rate control loops, this optimization scheme can be easily adopted to existing heat exchanger

Table 7. Heat Exchanger Data of Branch B parameter tube-side flow heat transfer area (m2) number of tube and shell passes total number of tubes inner and outer diameter (mm) tube wall thermal conductivity (W/(m ·K)) tube-side heat transfer coefficient at the base case (W/(m2·K)) shell-side heat transfer coefficient at the base case (W/(m2·K))

heat exchanger 12

heat exchanger 13

heat exchanger 14

cold stream 400 1−2



1500 14.85/19.05

1425 14.85/19.05

1556 14.85/19.05

55

55

55

699

736

674

624

628

622

is a stabilization of the flow split and the heat load of both branches decreases due to fouling at similar rates. The gain achieved with the dynamic optimization is shown in Table 8, where it is possible to observe that, during the 52 weeks analyzed, the optimization reaches a reduction of 1.5% on the energy consumption due to fouling.

5. CONCLUSIONS This paper presented the dynamic optimization of the distribution of the flow rates in heat exchanger networks for fouling mitigation. Because the fouling rate is affected by the flow rate distribution, the problem must be formulated

Figure 8. Mass flow rates of the cold streams for the base and optimized cases. 6505


Article


Figure 9. Heat load recovery.

HSTR = subset of hot streams i = interest rate kw = thermal conductivity of the tube wall, W/(m K) m = stream mass flow rate, kg/s MX = subset of mixers n = network inlet/outlet mass flow rate, kg/s np = number of past periods Ntp = number of tubes per pass NTU = number of transfer units, dimensionless p = weights of the numerical integration procedure P = effectiveness, dimensionless Pr = Prandtl number PD = subset of demand units PS = subset of supply units Q = heat load, W R = universal gas constant, (Pa m3)/(mol K) Re = Reynolds number, dimensionless Rf = fouling resistance, m2K/W s = present worth factor S = subset of streams SCC = subset of countercurrent heat exchangers SMP = subset of multipass heat exchangers SP = subset of splitters STR = set of edges SPEm = set of bounds on flow rates SPET = set of bounds on temperatures t = time, s T = stream temperature, °C TI = set of time instants U = overall heat transfer coefficient, W/(m2K) v = fluid flow velocity, m/s V = network inlet/outlet temperature, °C VET = set of vertices

Table 8. Total Energy Consumption case

total energy required (PJ)

average energy consumption rate (kW)

base optimized

0.153 0.151

4865 4801

networks, thus becoming a viable alternative for fouling mitigation.

■

AUTHOR INFORMATION

Corresponding Author

*E-mail: [email protected]. Notes

The authors declare no competing financial interest.

■

ACKNOWLEDGMENTS The authors thank the National Council for Scientific and Technological Development (CNPq), Coordination for the improvement of Higher Education Personnel (Capes), and Petrobras for the financial support.

■

NOMENCLATURE

List of Symbols

A = heat exchanger area, m2 b = parameter associated with constraint in eq 44 c = parameter which indicates if the cold stream flows in the tube-side (= 1) or shell-side (= 0) Cop = utility costs, $/J Cp = stream heat capacity, J/(kgK) CP = network inlet/outlet stream heat capacity, J/(kgK) CR = ratio between heat capacity flow rates, dimensionless CSTR = subset of cold streams dRfdt = fouling rate, m2K/J De = outer tube diameter, m Di = inner tube diameter, m Ea = activation energy, J/mol f = Fanning friction factor, dimensionless fobj = objective function, $ h = film coefficient, W/(m2K) HE = subset of heat exchangers

Superscripts

base = base case in = edges into a vertex out = edges from a vertex LB = lower bound UB = upper bound spe = specifications 6506


Article

Industrial & Engineering Chemistry Research Subscripts

(15) Wang, Y.; Smith, R. Retrofit of a heat-exchanger network by considering heat-transfer enhancement and fouling. Ind. Eng. Chem. Res. 2013, 52, 8527. (16) Pan, M.; Bulatov, I.; Smith, R. Exploiting tube inserts to intensify heat transfer for the retrofit of heat exchanger networks considering fouling mitigation. Ind. Eng. Chem. Res. 2013, 52, 2925. (17) Oliveira Filho, L. O.; Liporace, F. S.; Queiroz, E. M.; Costa, A. L. H. Investigation of an alternative operating procedure for fouling management in refinery crude preheat trains. Appl. Therm. Eng. 2009, 29, 3073. (18) Assis, B. C. G.; Gonçalves, C. O.; Liporace, F. S.; Oliveira, S. G.; Queiroz, E. M.; Pessoa, F. L. P.; Costa, A. L. H. Constrained thermohydraulic optimization of the flow rate distribution in crude preheat trains. Chem. Eng. Res. Des. 2013, 91, 1517. (19) Ishiyama, E. M.; Paterson, W. R.; Wilson, D. I. Platform for techno-economic analysis of fouling mitigation. Energy Fuels 2009, 23, 1323. (20) Lid, T.; Strand, S.; Skogestad, S. On-line optimization of a crude unit heat exchanger network. Symposium Chemical Process Control 6, Tucson, Arizona, USA, 2001. (21) Ishiyama, E. M.; Paterson, W. R.; Wilson, D. I. Thermohydraulic channeling in parallel heat exchangers subject to fouling. Chem. Eng. Sci. 2008, 63, 3400. (22) Mah, R. S. H. Chemical Process Structures and Information Flows; Butterworth Publishers: Markham, 1990. (23) Shah, R. K.; Sekulic, D. P. Fundamentals of Heat Exchanger Design; John Wiley and Sons: Hoboken, 2003. (24) Incropera, F. P.; De Witt, D. P. Fundamentals of Heat and Mass Transfer; John Wiley & Sons: New York, 2002. (25) Mukherjee, R. Effectively design shell-and-tube heat exchangers. Chem. Eng. Prog. 1998, 94, 21. (26) Bott, T. R. Fouling of Heat Exchangers; Elsevier Science: Amsterdam, 1995. (27) Panchal, C. B.; Huangfu, E. Effects of mitigating fouling on the energy efficiency of crude-oil distillation. Heat Transfer Eng. 2000, 21, 3. (28) Wilson, D. I.; Polley, G. T.; Pugh, S. J. Ten Years of Ebert, Panchal and the ‘Threshold Fouling’ Concept. ECI Symp. Ser. 2005, RP2, 25. (29) Panchal, C. B.; Kuru, W. C.; Liao, C. F.; Ebert, W. A.; Palen, J. W. Threshold Conditions for Crude Oil Fouling. In Understanding Heat Exchanger Fouling and Mitigation; Bott, T. R., Melo, L. F., Panchal, C. B., Somerscales, E. F. C., Eds.; Begell House: New York, 1999; pp 273−279. (30) Couper, J. R. Process Engineering Economics; Marcel-Dekker: New York, 2003. (31) Lovik, I.; Hillestad, M.; Hertzberg, T. Long term dynamic optimization of a catalytic reactor system. Comput. Chem. Eng. 1998, 22 (Suppl.), S707. (32) Zahedi, G.; Elkamel, A.; Lohi, A. Dynamic optimization strategies of a heterogeneous reactor for CO2 conversion to methanol. Energy Fuels 2007, 21, 2977.

1 = heat exchanger extremity 2 = heat exchanger extremity film = film temperature ref = reference temperature s = surface temperature tube = heat exchanger tube-side shell = heat exchanger shell-side k = edge index k′ = alias of index k k″ = alias of index k t = vertex index τ = time index Greek Letters

α = empirical parameter of the fouling rate equation, m2K/J γ = empirical parameter of the fouling rate equation, m2K/(J Pa) μ = fluid viscosity, Pa·s ρ = fluid density, kg/m3 τw = shear stress, Pa

■

REFERENCES

(1) Furman, K. C.; Sahinidis, N. V. A critical review and annotated bibliography for heat exchanger network synthesis in the 20th century. Ind. Eng. Chem. Res. 2002, 41, 2335. (2) Smaïli, F.; Vassiliadis, V. S.; Wilson, D. I. Mitigation of fouling in refinery heat exchanger networks by optimal management of cleaning. Energy Fuels 2001, 15, 1038. (3) Smaïli, F.; Vassiliadis, V. S.; Wilson, D. I. Optimization of cleaning schedules in heat exchanger networks subject to fouling. Chem. Eng. Commun. 2002, 189, 1517. (4) Smaïli, F.; Vassiliadis, V. S.; Wilson, D. I. Long-term scheduling of cleaning of heat exchanger networksComparison of outer approximation-based solutions with a backtracking threshold accepting algorithm. Chem. Eng. Res. Des. 2002, 80, 561. (5) Lavaja, J. H.; Bagajewicz, M. J. On a new MILP model for the planning of heat exchanger network cleaning. Ind. Eng. Chem. Res. 2004, 43, 3924. (6) Lavaja, J. H.; Bagajewicz, M. J. On a new MILP model for the planning of heat exchanger network cleaning. Part II: Throughput loss considerations. Ind. Eng. Chem. Res. 2005, 44, 8046. (7) Lavaja, J. H.; Bagajewicz, M. J. On a new MILP model for the planning of heat exchanger network cleaning. Part III: Multiperiod cleaning under uncertainty with financial risk management. Ind. Eng. Chem. Res. 2005, 44, 8136. (8) Rodriguez, C.; Smith, R. Optimization of operating conditions for mitigating fouling in heat exchanger networks. Chem. Eng. Res. Des. 2007, 85, 839. (9) Ishiyama, E. M.; Heins, A. V.; Paterson, W. R.; Spinelli, L.; Wilson, D. I. Scheduling cleaning in a crude oil preheat train subject to fouling: Incorporating desalter control. Appl. Therm. Eng. 2010, 30, 1852. (10) Ishiyama, E. M.; Paterson, W. R.; Wilson, D. I. The effect of fouling on heat transfer, pressure drop, and throughput in refinery preheat trains: Optimization of cleaning schedules. Heat Transfer Eng. 2009, 30, 805. (11) Polley, G. T.; Wilson, D. I.; Yeap, B. L.; Pugh, S. J. Use of crude oil fouling threshold data in heat exchanger design. Appl. Therm. Eng. 2002, 22, 763. (12) Butterworth, D. Design of shell-and-tube heat exchangers when the fouling depends on local temperature and velocity. Appl. Therm. Eng. 2002, 22, 789. (13) Wilson, D. I.; Polley, G. T.; Pugh, S. J. Mitigation of crude oil preheat train fouling by design. Heat Transfer Eng. 2002, 23, 24. (14) Yeap, B. L.; Wilson, D. I.; Polley, G. T.; Pugh, S. J. Retrofitting crude oil refinery heat exchanger networks to minimize fouling while maximizing heat recovery. Heat Transfer Eng. 2005, 26, 23. 6507


Dynamic Optimization of the Flow Rate Distribution in Heat Exchanger

Recommend Documents