A Dynamic, Distributed Model of Shell-and-Tube Heat Exchangers

ARTICLE pubs.acs.org/IECR

A Dynamic, Distributed Model of Shell-and-Tube Heat Exchangers Undergoing Crude Oil Fouling Francesco Coletti and Sandro Macchietto* Department of Chemical Engineering and Chemical Technology, Imperial College London, South Kensington Campus, London SW7 2AZ, U.K. ABSTRACT: Fouling in refinery preheat trains causes major energy inefficiencies, resulting in increased costs, greenhouse gas emissions, maintenance efforts, and safety hazards. Fouling deposition is not well understood, and current exchanger design and monitoring practices neglect the local effects and dynamics of fouling, in favor of lumped, steady-state, heuristic models (e.g., Tubular Exchanger Manufacturers Association (TEMA) fouling factors). In this paper, a dynamic, distributed model for a multipass shell-and-tube heat exchanger undergoing crude oil fouling is proposed. The model calculates fouling rates as a function of local conditions and time. It accounts for heat exchanger geometry, variation of oil physical properties with temperature, local accumulation of fouling deposits, and their structural change over time (aging). The interaction between fouling growth and fluid dynamics inside the tubes is accounted for by solving a moving boundary problem. Moreover, a procedure to analyze refinery data and support the estimation of a set of model parameters has been established. The model is validated using data from an ExxonMobil refinery and shows excellent agreement (less than 2% error) with primary plant measurements even when it is tested for its predictive capabilities over long periods (i.e., 1 year). It is concluded that the model can be used with confidence to identify and predict the fouling state of exchangers, to assess economic losses due to fouling, to support operating decisions such as planning of cleaning schedules, and to assist in the design and retrofit of heat exchangers.

1. INTRODUCTION In an oil refinery, approximately 6% of the total crude throughput is consumed as energy by the refining process itself.1 About two-thirds of this energy2 is used in the furnace of the crude distillation unit (CDU), which provides the enthalpy required for primary fractionation. An extensive network of heat exchangers; the preheat train (PHT);is used to reduce this requirement by recovering as much energy as possible from the hot product streams of the distillation column (Figure 1). Although a few examples exist3 of refineries using some compact plate and frame design, a typical PHT is made up of as many as 60 shelland-tube heat exchangers.4 Depending on the thermal performance of the PHT, the cold (i.e., ambient temperature) inlet crude oil can be preheated before entering the furnace to a typical range5 of 270-290 C. This temperature is known as the coil inlet temperature (CIT). The efficiency of the PHT is paramount for the whole refinery as it recovers nearly 60-70% of the duty required for distillation.4 For this reason, individual exchangers in the PHT, the exchanger network structure, and their operation are typically optimized extensively through both heuristic and analytical methodologies.6,7 However, over time, layers of unwanted material (fouling) deposit on the heat transfer surfaces, hindering thermal and hydraulic efficiency. Traditional design methodologies, based on empirical Tubular Exchanger Manufacturers Association (TEMA) fouling factors,8 have received several critiques9-11 as they fail to account for the dynamic nature of fouling as well as its dependence on process conditions and composition. Moreover, technologies used for network design, such as pinch analysis,12 do not take into account the progressive deterioration of performance caused by fouling. As a result, r 2011 American Chemical Society

benefits obtained by optimizing a network with such techniques can be greatly reduced and in some cases fouling is actually exacerbated.13 1.1. Impact of Fouling on Crude Distillation Units. Problems related to fouling in the crude distillation unit can be divided into four main areas: 1. operating difficulties 2. economic penalties (summarized in Table 1) 3. environmental impact 4. safety hazards when opening the units for cleaning Several mechanisms among those recognized by Epstein14 are responsible for fouling in PHTs. For example, it is not unusual to find deposits due to biofouling at the cold end of the train whereas corrosion fouling afflicts the whole train. However, fouling is particularly severe at the hot end where there is evidence1 that chemical reactions, triggered by the high temperatures, are responsible for deposition. These deposits have a thermal conductivity;typically between15 0.2 and 1 W m-1 K-1;up to 2 orders of magnitude lower than that of the tube metal wall (ca. 38 W m-1 K-1, depending on the metal used16). As a consequence, heat transfer is impaired, less energy is recovered, and the CIT decreases over time with a typical drop17 of 8-11 C/year. To maintain a constant temperature at the distillation column inlet, the decline in CIT is countered by burning additional fuel in the furnace downstream from the PHT. This not only adversely affects energy costs but also has an environmental impact due to Received: December 15, 2009 Accepted: December 20, 2010 Revised: December 15, 2010 Published: March 10, 2011 4515

dx.doi.org/10.1021/ie901991g | Ind. Eng. Chem. Res. 2011, 50, 4515–4533

Industrial & Engineering Chemistry Research

ARTICLE

the extra release of greenhouse gases to the atmosphere.18 A study in 2006 for the U.S. Department of Energy19 indicates that potential fuel savings up to 55% can be achieved by improving operating practices and capital equipment in refineries. Among the suggested improvements, it was found that fouling mitigation in PHT and fired heater could lead to 15% fuel savings.

Figure 1. Schematic of a typical crude distillation unit (CDU).

Considering that fuel consumption at the atmospheric furnace represents 4% of a refinery total throughput,2 15% of the 2008 global oil consumption of 85 Mbbl/day equates worldwide to a staggering 500 000 bbl/day potential savings, the daily production of a very large refinery. Decreased thermal efficiency also has a direct consequence on refinery throughput. As the furnace duty increases with diminishing efficiency, at some point the furnace typically reaches its maximum load (firing limit). As a consequence, throughput must be reduced, causing major economic losses quantified in 1981 as U.S. $8,500,000 per year for a then typical oil refinery processing 100 000 bbl/day (data from Van Nostrand et al.20 and adjusted for inflation, but not oil prices, to 2009). The effects of fouling on exchangers’ hydraulics can also be relevant. The progressive reduction of the cross-sectional area available to the crude flow causes an increase in pressure drops. If the throughput is to be kept constant, extra pumping power is required, resulting in an increase of electricity consumption. Additionally, the pumps must be oversized to meet this energy duty, adding to capital costs. Although it is the thermal impact of fouling that typically limits throughput, cases exist in which it is the pressure drop in the network that forces throttle back of production. The oil industry has faced the crude oil fouling problem for decades.21 Monitoring the performance of PHTs22,23 as well as managing the cleaning schedule of the units24 is a focus of primary concern for refinery operators. Although several studies are available in the literature on the topic,25-27 crude oil fouling mechanisms remain ill-defined. The number of distinct phenomena involved and the complexity of their interactions (summarized in Figure 2) has so far limited the accuracy of fouling predictions. Predictions of fouling in refineries are typically merely based on current trends; neither crude composition nor process conditions are accounted for and managerial decisions are taken based on simple calculation and past experience, aided at most by simple lumped models.

Table 1. Losses Associated with Fouling in a Typical Oil Refinery Processing 100 000 U.S. bbl/daya type of loss whole refinery

source Van Nostrand et al.20

year refinery size (bbl/day) 1981

100 000

original value U.S. $9,875,000 p.a. due to fouling in

adjusted to 2009 (1000 U.S. $) 22,600 p.a.

the whole refinery whole PHT

Van Nostrand et al.20

1981

100 000

whole PHT

personal communication

2009

210 000

U.S. $2,500,000 p.a.

one branch of the PHT

Liporace and De Oliveira23 2007

360 000

U.S. $21,800,000 p.a.

-

throughput reduction

Bories and Patureaux36

2000

160 000

U.S. $1.5 million over 3 months

1,810 over 3 months

1981 1981

100 000 100 000

U.S. $3,730,000 p.a. U.S. $1,020,000 p.a. without

Panchal and Huangfu4

2000

100 000

U.S. $400,000 p.a.

510 p.a.

Yeap et al.2

2003

100 000

750 tons/year of extra CO2 released

31 p.a.

U.S. $4,785,000 p.a. (48.5% of the total) 11,000 p.a. -

after startup

in total refinery throughput reduction Van Nostrand et al.20 extra fuel burnt in the furnace Van Nostrand et al.20

8,500 p.a. 2,300 p.a.

antifoulantsb

due to drop in CIT

(h22,500 p.a.)c a drop of 1 C in CIT

81

Macchietto et al.

2009

200 000

U.K. £250,000 p.a.

400 p.a.

Baudelet and Krueger82

1998

n/a

1 ton/day of CO2 released equal to

130 p.a.

production loss taking a heat Baudelet and Krueger82

1998

n/a

U.S. $100,000 p.a. U.S. $20,000 per day = U.S. $80,000

24 per day = 100 per cleaning

exchanger out of service

a

per cleaning

cleaning of a single unit

BP17

2006

n/a

U.S. $30,000-50,000 per unit

30-50 per unit

antifoulants

Van Nostrand et al.20

1981

100 000

U.S. $155,000 p.a.

350 p.a.

All values in 2009 U.S.dollars. p.a., per annum. b The process fuel cost used is U.S. $2.80/million Btu. c The carbon tax rate considered is h30/metric ton. 4516

dx.doi.org/10.1021/ie901991g |Ind. Eng. Chem. Res. 2011, 50, 4515–4533


ARTICLE

Figure 2. Phenomena occurring in the tube-side of a shell-and-tube heat exchanger undergoing chemical reaction fouling. Boxes at the top of the figure show the phenomena whose increase reduces fouling rates and has a positive impact on the measured quantities: pressure drop, throughput, and coil inlet temperature. Boxes at the bottom of the figure show the ones that should decrease in order to mitigate fouling. Dashed boxes indicate phenomena not accounted for in the model (however, see Coletti et al.63 for incorporation of roughness). Arrows show the direction of positive impact between the phenomena.

To overcome the limitations given by the fouling factor approach, Kern and Seaton28 followed by Taborek et al.,21 Epstein,29 and Crittenden et al.30 carried out some pioneering work to describe the time dependency of the fouling resistance on process conditions through mathematical relationships. These represent a first attempt to move beyond the fixed fouling factor approach and put the basis for more complex models reviewed in section 1.2. 1.2. Thermal Fouling Models. It has been demonstrated experimentally31,32 that, at certain velocities, the suppression mechanism governed by the shear stress outbalances the deposition mechanism resulting in no fouling. Models based on this threshold concept32 have been used to take into account crude oil fouling when designing or retrofitting single units33-35 and whole PHTs2,36 and to assist in the scheduling of cleaning, thus improving network operability and mitigating costs related to maintenance.37-39 However, all these works focus on the fouling model itself with little attention paid to the underlying heat transfer model. For example, the above are all lumped models that do not take into account local variation in conditions along the exchanger units. As a consequence, large multipass multishell units cannot be properly represented. Indeed, it has been recognized that one of the main limitations of lumped models is the lack of accuracy that results by considering a constant heat transfer coefficient.40 Furthermore, in previous models, steady state energy balance equations have been considered and a quasisteady-state approximation is required.41 In many cases, a thin slab assumption is also used in the calculation of the overall thermal resistance of the fouling layer. Such an assumption is only valid when the maximum thickness of the fouling layer is about 10% of a typical inner tube diameter.2 It is therefore not surprising that all these limitations often lead to large disagreements between simulation results and plant data.

The literature shows that some effort has been put into developing accurate dynamic models of shell-and-tube heat exchangers in order to study the transient response of these units.42-45 Some papers focus on shell-side flow modeling through a threedimensional46,47 or axially dispersed plug flow,48 in an attempt to improve the heat transfer coefficient predictions given by traditional methodologies such as the Bell-Delaware method or flow stream analysis.49 Roetzel and Xuan50 proposed dynamic models for several types of heat exchangers including plate, shell-andtube, and cross-flow heat exchangers. They suggested the use of a simple asymptotic model to represent fouling combined with a detailed, distributed model. However, their work focused primarily on studying the dynamic response to step and ramp disturbances and no results are reported on the fouling behavior. Indeed, only few papers focus on a dynamic model for heat exchangers which incorporates fouling. Fryer and Slater51 developed a model one-dimensionally distributed on the shell-side. They included chemical reaction fouling but applied the model to milk fouling, in which case they could estimate directly from experiments the activation energy of the fouling reaction. However, their model does not account for the heat capacity of the wall and is not for multipass heat exchangers. Georgiadis et al. used the milk fouling model by Toyoda and Fryer,52 coupled it with a distributed heat exchanger model where thermal and transport phenomena are taken into account in detail,53,54 and used it to optimize design and operations.55,56 Even fewer attempts reported in the literature calculate the effects of shell-side fouling because of the complexity related to the calculations of flow patterns and the wall shear stress which, on the shell-side, cannot be inferred by pressure drop calculations. A computationally intensive computational fluid dynamics (CFD) simulation is therefore required. This kind of approach 4517



ARTICLE

Figure 3. Definition of model domains.

has been undertaken by Clarke and Nicolas,57 who showed how CFD work can be directly related to plant data to explore the effects of shell-side fouling. To reduce the computational effort required to simulate the 14-baffle unit considered, the shell volume was modeled as a porous medium. The tube bundle was represented by a reduction in porosity of the shell-side porous volume. It should be noted that, on the tube-side, an arbitrary linear temperature rise along the tube was set and the heat transfer coefficient was considered constant. For this reason, the interactions between tube-side and shell-side were neglected. Nonetheless, to the authors’ knowledge this represents, remarkably, the only attempt to simulate shell-side fouling in refinery shell-and-tube heat exchangers. Some examples of works using neural networks to model fouling in PHTs are also reported in the literature.58,59 However, they require extensive data for “training” the model and appear to have limited extrapolation capabilities compared to either first principles or empirical models. 1.3. Objectives. From the literature survey presented in section 1.2, there seem to exist either lumped models of low accuracy in the description of dynamics and distributed phenomena or dynamic and distributed models which do not incorporate fouling or use simple models. Moreover, the trend in the literature is to fit models not to primary measured data such as temperatures but to fouling resistances calculated using the classic logarithmic mean temperature or ε-number of transfer units (NTU) methods. Such derived fouling resistance embeds all the uncertainties not only in the measurements but also in the assumptions, such as constant density and specific heat capacity, used for its calculation. In this paper, a more comprehensive approach is proposed which aims at the following: 1. Develop a dynamic mathematical model capable of describing tube-side crude oil fouling in shell-and-tube heat exchangers as a function of local conditions throughout the unit. 2. Devise a procedure to systematically analyze plant data and estimate necessary model parameters using primary plant measurements (i.e., temperatures and flow rates) rather than derived fouling resistances. 3. Validate the model with plant measurements and testing it for its predictive capabilities again, against primary quantities directly measured (temperatures). In section 2 the reference system adopted for a generic shelland-tube heat exchanger is described. Model equations are outlined in section 3 for each domain considered, while section

4 deals with some numerical aspects of the simulations. A systematic procedure to analyze refinery data is reported in section 5. Section 6 introduces the case study for an ExxonMobil refinery which is used to validate the model in section 7, where simulation results are compared with plant data. Finally, section 8 discusses some possible applications of the model.

2. SYSTEM DEFINITION The mathematical model proposed in this paper is for a multipass tubular heat exchanger undergoing crude oil fouling on the tube-side. It is acknowledged that fouling can also occur on the shell-side (for example, when the heating fluid is a heavy petroleum fraction or if crude is switched to the shell-side for design reasons). Unless independent pressure drop measurements (typically not available) are used, it is not easy to determine with certainty whether fouling is occurring on the tube-side only. However, visual inspection of exchanger units when they are dismantled for cleaning indicate that very often tube-side fouling is the dominant resistance to heat transfer. This is due to a number of reasons. First, the most foulant liquid is often allocated on the tube-side for ease of cleaning purposes. Second, because of the higher level of turbulence and shear achievable the fouling resistance is smaller than on the tube-side. The exchanger considered for the model is of TEMA type AET (i.e., single pass shell, front end stationary head with removable cover and floating head at the rear end) of length L. However, the model can be used for shell type E with different front and rear head types (e.g., bonnet, floating head, U tube bundle, etc.) but will require some adaptation for different shell types (e.g., two pass with longitudinal baffle, F shell, or divided flow, J shell). Figure 3 depicts the physical system, which is divided into four distinct control volumes (domains): ΩS: shell-side domain, the volume of the exchanger shell outside the tubes ΩW,n: wall domain, between the tube’s inner radius, Ri, and the outer one, Ro ΩL,n: fouling layer domain, defined between the oil/deposit interface of the fouling layer, Rflow, and the inner radius of the tube, Ri ΩT,n: tube-side domain, defined between the center of a tube and Rflow The subscript n is used to indicate the pass number and varies between 1 and the total number of passes on the tube-side, Np. The model is one dimensionally distributed along the z axis for the tube-side and shell-side domains (ΩT and ΩS). In the fouling layer and wall domains (ΩL and ΩW) the model is distributed 4518



ARTICLE

(Figure 4) based on first principles:61 rate of energy accumulation in control volume 1 ¼ net energy transfer by fluid flow 2 - net energy transfer by conduction 3 þ rate of internal heat generation 4 - fnet work transfer to environmentg5

ð1Þ

In the following sections the governing equations of the model are given, together with suitable boundary and initial conditions. 3.1. Shell-Side. The heat balance on the shell as the control volume (domain ΩS defined for z = [0, L]) is as follows. ∂ ∂ ∂ ∂TS ðF cp, S TS Þ ¼ - dirS ðFS cp, S TS uS Þ þ λS ∂t S ∂z ∂z ∂z þ

Figure 4. Reference system. Schematic representation of a one-shell, two-tube pass heat exchanger in the case of first tube pass flow in cocurrent (a) and countercurrent (b) with the shell-side flow. The heat transfer coefficient for the shell-side takes into account deviation from cocurrent/countercurrent flow due to cross-flow.

along both the axial and the radial coordinates in order to capture important phenomena occurring in the radial direction. An example of such phenomena is the aging of the deposits outlined in section 5. Being distributed only axially, the shell-side is assumed to have the same temperature radially, at a given z coordinate in the heat exchanger. Each tube in the same pass is assumed to experience the same temperature difference relative to the shell. The driving force for heat exchange is therefore the same for each tube in the same pass n. As all tubes in each pass behave similarly, the system can be conceptualized as shown in Figure 4, where a single tube per pass is actually modeled and considered representative of all other tubes in the same pass (however, in practice the shell-side flow is not parallel to the tubes, as discussed later). Headers are considered perfectly mixed, and all entrance effects at the tube entrances and exits (i.e., z = 0 and z = L) and heat losses are neglected. The physical properties of the oils flowing in both tube-side and shell-side, such as density F, viscosities μ and ν, thermal conductivity λ, and heat capacity cp, used in the heat balance and the heat transfer relationships are calculated using API relationships60 as a function of temperature and space.

3. GOVERNING EQUATIONS The temperature profiles in each domain considered are given by a thermal balance on a differential control volume Δz

Np 1 X PS, n hS ðTS - TW , n jr ¼ Ro Þ AS n ¼ 1

ð2Þ

where TS(z) denotes the shell-side fluid temperature whereas FS(z), cp,S(z), and λS(z) denote, respectively, its density, specific heat capacity, and thermal conductivity as a function of the axial coordinate. The last term is the “heat source” term (exchanged from the tubes) where Np is the number of tube passes, PS,n is the wetted perimeter from which heat is exchanged, and AS is the cross-sectional area of the shell. Furthermore, the net work from the control volume to its environment (fifth term in eq 1) is neglected. In eq 2, a variable direction term, dirS = (1, is introduced to take into account the internal arrangement with respect to the direction of the flow in the first tube pass; dirS = 1 if the first tube pass is in cocurrent flow (Figure 4a) and dirS = -1 otherwise (Figure 4b). The shell-side heat transfer coefficient hS and pressure drop ΔPS are calculated through the Bell-Delaware method:49 hs ¼ hid Jc Jl Jb Js Jr

ð3Þ

where hid is the heat transfer coefficient for ideal cross-flow, Jc, Jl, Jb, Js, and Jr are correction factors for, respectively, the segmental baffle window, the baffle leakage, the bypass tube bundle to shell, the laminar heat transfer and the nonequal inlet/outlet baffle spacing. Details of the calculations are given by Taborek49 and are not reported here. The only difference from the standard calculation is that the ideal heat transfer coefficient is explicitly considered as a function of the spatial position: hid ¼ hid ðzÞ

ð4Þ

It should be noted that, while in eq 2 the fluid velocity is assumed to be parallel to the tubes, the heat transfer coefficient in eq 3 does take into account the variations in heat transfer produced by a more complex flow path. This is clearly an approximation with respect to describing the shell-side flow field, used to reduce computational effort. 3.2. Tube Wall. For the tube wall (domain ΩW, defined for z = [0, L] and r = [Ri, Ro]), the standard heat conduction equation is used for each pass, n, assuming no variation along the z direction, symmetry with respect to the angular coordinate θ, and no heat sources: FW cp, W ∂TW, n 1 ∂TW, n ∂2 TW , n ¼ þ r ∂r λW ∂t ∂r 2 4519

ð5Þ



ARTICLE

where TW,n denotes the temperature which is a function of spatial coordinates (i.e., TW,n = TW,n(z,r)) while λW is the thermal conductivity, FW is the density, and cp,W is the heat capacity of the metal wall which are assumed constant. The heat flux for each pass, q00W,n(z,r), is calculated as 00

q W , n ¼ - λW

∂TW , n ∂r

ð6Þ

3.3. Fouling Deposit Layer. The deposit layer, domain ΩL (defined for z = [0, L] and r = [Ri,Rflow,n]), is also modeled as a conductive domain neglecting variations in the axial direction and any heat sources, and assuming symmetry with respect to angular coordinate θ: ∂TL, n ∂TL, n 1 ∂ rλL, n ¼ ð7Þ FL cp, L r ∂r ∂t ∂r

In eq 7, TL,n denotes the local layer temperature and λL,n its thermal conductivity, both of which are a function of spatial coordinates z and r (i.e., TL,n(z,r) and λL,n(z,r)). Conversely, the density of the deposit layer, FL, and its heat capacity, cp,L, are here assumed uniform and constant. It should be noted that the model formulation in cylindrical coordinates allows accounting for curvature effects in the heat flux, thus overcoming the thin slab approximation used in other models. Equation 7 also accounts for the variation of the deposit thermal conductivity in the radial direction. The radial heat flux at any point in the deposit layer, q00L,n(z,r), is calculated as 00

qL, n ¼ - λL, n

∂TL, n ∂r

ð8Þ

Equations 7 and 8 are defined in the domain between Ri and the moving boundary coordinate Rflow,n(z): Rflow, n ¼ Ri - δn

ð9Þ

where the deposit thickness, δn, is a function of the axial coordinate (i.e., δn = δn(z)) and is calculated as shown below by eq 29. In order to solve the set of partial differential equations (PDEs) with a moving boundary, eqs 7 and 8 are reformulated in dimensionless form with respect to the radial coordinate, utilizing a coordinate transformation. The following dimensionless radial coordinate is introduced: ~r ¼

r - Ri Rflow, n - Ri

ð10Þ

So that ~r ¼ 0

at

r ¼ Ri ,

at

r ¼ Rflow, n ,~r ¼ 1

ð11Þ

From the definition of the flow radius in eq 9, the new radial coordinate can also be written as ~r ¼

Ri - r δn

ð12Þ

In terms of the new radial coordinate, ~r = [0,1], eq 7 becomes δn 2 FL cp, L

Similarly, eq 8 becomes 00

qL, n ¼

∂2 TL, n ∂~r 2

ð14Þ

3.3.1. Aging Model. Exposure of the fouling layer to wall temperatures over extended periods triggers chemical modifications of the deposit resulting in an alteration of its physical structure from a soft, gel-like material to a harder, coke-like one. This process takes the name of “aging”.14 The aging process for a given layer starts from the time when deposition of that layer occurred. Ultimately this results in a change in the deposit’s physical properties (e.g., viscosity, thermal conductivity, etc.) across its radial direction. Ishiyama et al.62 developed a kinetic model for aging in chemical reaction fouling which was extended to a distributed system by Coletti et al.63 Details of the aging model are given in those publications, but the main equations are briefly summarized here. In the model, aging is expressed via a change in the deposit thermal conductivity, from that of a freshly deposited material (the lower limit), λ0L, to that of a fully coked one with a maximum value, λ¥ L . In the dimensionless domain: λL, n ¼ λ¥L þ ðλ0L - λ¥L Þyn

ð15Þ

where yn(z,~r ) is a “youth” variable with value of 1 when the layer is initially deposited. Its evolution is calculated assuming first order exponential decay dynamics defined by ! dyn Ea ¼ - Aa exp ð16Þ yn dtage Rg TL, n where Ea and Aa are, respectively, the activation energy and preexponential constant characteristic of the aging process whereas Rg is the universal gas constant. A continuous formulation is obtained by introducing an aging time, tage, defined as tage ¼ tð1 -~r Þ

ð17Þ

It should be noted that, while this time transformation is rigorously valid in the case of linear fouling, it has been shown63 that deviation from this behavior only introduces an error comparable to that of other correlations used in heat exchange models (e.g., shell-side heat transfer coefficient). Using eq 17, eq 16 becomes " ! # dyn Ea ¼ ð1 -~r Þ - Aa exp yn ð18Þ dt Rg TL, n Differential eq 18 is expressed in the same time coordinate as the other differential equations and can therefore be integrated simultaneously with them. 3.4. Tube-Side. Assuming no variation in the radial and angular direction and neglecting the net work from the control volume to its environment (fifth term in eq 1), the heat balance equation for the tube-side domain can be written as

∂TL, n ∂λL, n ∂TL, n λL, n δn ∂TL, n ¼ ∂t ∂~r ∂~r Ri -~r δn ∂~r þ λL , n

λL, n ∂TL, n δn ∂~r

∂ ∂ ðFn cp, n Tn Aflow, n Þ ¼ - dirn ðFn cp, n Tn un Aflow, n Þ ∂t ∂z ∂ ∂Tn λn Aflow, n þ þ Pn hn ðTL, n jr ¼ Rflow - Tn Þ ð19Þ ∂z ∂z

ð13Þ 4520



ARTICLE

where Tn(z) denotes the tube-side fluid bulk temperature whereas Fn(z), cp,n(z), and λn(z) are, respectively, its density, specific heat capacity, and thermal conductivity. A variable direction term, dirn = (1, is introduced to take into account the direction of the flow, which depends on the pass: dirn = 1 in the case of odd pass; dirn= -1 in the case of even pass. It is noted that no assumptions have been made on the physical properties (e.g., constant density, thermal conductivity, or heat capacity), which do vary as a function of local conditions. Equation 19 takes into account the variation of the cross-sectional area, Aflow,n(z), which is a function of the flow radius, Rflow,n (z) (defined in eq 9): Aflow, n ¼ πRflow, n

ð20Þ

2

3.4.1. Hydraulics. The hydraulic effect of fouling (velocities and pressure drop changes due to reduction of flow section in pipes) and its interaction with the fouling process itself is paramount. This is captured through the definition of the flow radius, Rflow, in eq 9 and the flow area, Aflow,n, in eq 20. The tube-side velocity, un(z), is a function of this flow radius: un ¼

m_ Aflow, n NT =Np Fn

ð21Þ

The pressure drop inside the tube also reflects the interaction with the growth of the fouling layer inside the tubes and is calculated by dPn Fun 2 ¼ 4Cf dz 4Rflow, n

ð22Þ

where Cf,n is the Fanning friction factor for rough tubes:2 Cf , n ¼ 0:0035 þ 0:264Ren

- 0:42

ð23Þ

In eq 23 Re is the Reynolds number calculated as a function of the axial coordinate (i.e., Re = Re(z)). Another effect of the interaction of the fouling layer with the fluid flow is through the heat transfer coefficient, hn(z), on the tube-side: hn ¼

λn Nun 2Rflow, n

ð24Þ

here to local conditions and applied to a distributed system to calculate the distributed fouling resistance, Rf,n(z): ! dRf , n - Ef - 0:66 - 0:33 ¼ RRen Pr n exp ð27Þ - γτn dt RTf , n where R, γ, and the activation energy Ef are adjustable parameters whereas Tf,n(z) is the film temperature calculated as Tf , n ¼ Tn þ 0:55ðTL, n jr ¼ Rflow - Tn Þ

It is noted that the model in eq 27 was derived for coke deposition in furnaces (i.e., at higher temperature than the usual PHT temperatures) where the fouling mechanism is different from the chemical reaction one believed to occur in the hot end of the PHT. Once the fouling resistance is known through eq 27, the fouling layer thickness can be calculated by dRf , n dδn ¼ λ0L dt dt

ð25Þ

Equation 24 is valid64 (with a standard deviation of 13%) for 0.7 < Pr < 160, Re > 104, and a ratio between the length of the pipe and its diameter of >10. It is noted from eq 24 that hn is enhanced by the progressive reduction in cross-sectional area. The wall shear stress, τn(z), is calculated as a function of the Reynolds number through the following:65 ! Fn u n 2 τ n ¼ Cf , n ð26Þ 2

ð29Þ

Equation 27 gives the fouling resistance for each pass at each axial coordinate, z. This disaggregated information is very useful to increase the accuracy in calculations of temperature profiles. However, in industry a more common measure of fouling is the average resistance. The average fouling resistance per pass, Rhf,n is calculated as follows: ! Z Z 1 1 L δn Rf , n ¼ d~r dz ð30Þ L 0 0 λL , n An overall average fouling resistance, useful for comparison with aggregate models can be defined:

Rf ¼

1 NS

NS X

Np P n¼1

Rf , n

Np L

i¼1

ð31Þ

where, for a multishell unit, NS is the number of shells in the unit. 3.5. Boundary Conditions. At the shell-side inlet, in case of a cocurrent flow with respect to the first tube pass (dirS = 1, Figure 4a) temperature and pressure are given by TS jz ¼ 0 ¼ TS, in PS jz ¼ 0 ¼ PS, in

where Nu(z), the Nusselt number, a function of the Reynolds and Prandtl numbers, is calculated though the Dittus-Boelter relationship:64 Nun ¼ 0:023Ren 0:8 Pr n 0:4

ð28Þ

ð32Þ

whereas in the case of a first tube pass in countercurrent with shell-side flow (dirS = -1, Figure 4b) TS jz ¼ L ¼ TS, in PS jz ¼ L ¼ PS, in

ð33Þ

At the interface between the shell-side domain, ΩS, and the tube wall, ΩW, there is continuity in the radial heat flux: 00

qW , n jr ¼ Ro ¼ - hS ðTS - TW , n jr ¼ Ro Þ

ð34Þ

At the interface between the wall, ΩW, and the fouling layer, ΩL, there is continuity in the heat flux and the temperature: 00

where Cf,n(z) is the Fanning friction factor, given by eq 23. 3.4.2. Fouling Model. For chemical reaction fouling in crude oil systems, a widely accepted model is the Ebert-Panchal one.32 The use of this model, modified by Panchal et al.,66 is extended 4521

00

qW , n jr ¼ Ri ¼ qL, n jr ¼ Ri

ð35Þ

TW , n jr ¼ Ri ¼ TL, n jr ¼ Ri

ð36Þ



ARTICLE

At the moving boundary between the fouling layer, ΩL, and the tube-side domain, ΩT, there is continuity in the heat flux: 00

qL, n jr ¼ Rflow, n ¼ - hS ðTL, n jr ¼ Rflow, n - Tn Þ

ð37Þ

At the inlet of the first tube pass (i.e., z = 0, n = 1): T1 jz ¼ 0 ¼ Tin P1 jz ¼ 0 ¼ Pin

ð38Þ

For subsequent even pass (dirn = -1 and n > 1): T n jz ¼ L ¼ T n - 1 jz ¼ L Pn jz ¼ L ¼ Pn - 1 jz ¼ L

ð39Þ Figure 5. Two-shell, four-pass heat exchanger arrangement. The overall flow arrangement in the unit is cocurrent. The first pass tube-side flow is cocurrent with the shell-side flow in shell A, whereas it is countercurrent in shell B.

whereas for an odd pass (dirn = 1 and n > 1): T n jz ¼ 0 ¼ T n - 1 jz ¼ 0 Pn jz ¼ 0 ¼ Pn - 1 jz ¼ 0

ð40Þ

3.6. Initial Conditions. At time 0 the heat exchanger is assumed to be clean (i.e., zero fouling resistance) in each pass, n:

Rf , n jt ¼ 0 ¼ 0

ð41Þ

The thickness of the fouling layer, δn, in each pass, n, is initialized to 10-7 m for numerical reasons: δn jt ¼ 0 ¼ 10 - 7

ð42Þ

The temperature profiles at time 0 are assumed to be at steady state in all domains: dTW, n dTL, n dTS dTn ¼ ¼ ¼ ¼ 0 ð43Þ dt dt dt dt t¼0

t¼0

t¼0

t¼0

The aging model given by eq 18 introduces a differential equation which needs an initial value of the youth variable: yn j t ¼ 0 ¼ 1

ð44Þ

In summary, the final model is represented by eqs 2-31, together with boundary conditions 32-40 and initial conditions 41-44. This model may be used to simulate the dynamic behavior and performance of a heat exchanger, once all (timevarying) inputs and model parameters are also specified. The inputs include inlet flow rates, temperatures, and pressures of the hot and cold streams, respectively. The model parameters include equipment characteristics (number of shells and passes, number of tubes, diameters, etc.), oil property parameters, and fouling and aging model parameters. Integration of the dynamic equations yields time profiles of all variables in the model, in particular outlet hot and cold stream temperatures, but also other important information, for example fouling deposit layer thickness and age in each part of the exchanger, pressure drops, local heat exchange coefficients, velocities, shear stress at the oil/ deposit interface, etc. This distributed information may also be aggregated to yield commonly used measures of performance such as (average) fouling resistance in each shell and overall.

4. SOLUTION METHOD The model outlined in section 3 comprises a set of partial differential and algebraic equations. It was implemented in the

gPROMS modeling environment67 exploiting its hierarchical structure. This allowed defining a single tube model once and reusing it as many times as needed within a single or multipass exchanger configuration, for the description of a full unit. An Excel interface via CAPE-OPEN68 was used to input the heat exchanger geometry parameters, making it extremely simple to set up the geometries of a unit under investigation. The same gPROMS software was used for the numerical solution of both simulation and parameter estimation steps using standard solution options. The partial differential equations are automatically discretized using a scheme of choice. Here, a second order centered finite difference scheme was selected with all domains discretized in 10 points in the axial direction. It should be noted that aging (discussed in section 5.3) produces steep profiles of the thermal conductivity in the radial direction in regions close to the interface between domains ΩL and ΩT. This may lead to numerical errors in the calculation of the heat flux q00L in eq 14. To ensure accuracy of the solution, a mesh independency test was performed for the most challenging conditions (combination of fast aging and fast fouling rates). For typical model parameters, Figure 6 reports the results of this test, showing that satisfactory accuracy is achieved with 1000 or more uniform radial discretization points. For slow aging conditions (not shown) fewer points are actually required to achieve a satisfactory accuracy. However, the computational load required to solve the set of PDEs is still substantial. A better alternative to achieve both numerical accuracy and fast solution is to discretize the radial direction domain using nonuniform grids. This technique allows increasing the number of grid points in regions where the gradients are larger (i.e., interface between domains ΩL and ΩT). Using an exponential coordinate transformation, the number of grid points required to achieve the same accuracy given by the uniform grid with 1000 elements is reduced to 150, thus increasing the efficiency of the solution. Figure 6 also shows a comparison between the uniform and nonuniform discretization one for different number of nodes. The wall domain ΩW does not present particular numerical challenges; therefore, a uniform discretization in the radial domain of 10 points was satisfactory.

5. MODEL PARAMETERS As previously noted, a number of model parameters must be specified. Some of these can be typically extracted from design data sheets (for example, the equipment geometric 4522



Figure 6. Mesh independency test. In the radial direction, an exponential grid with 150 elements gives the same degree of accuracy as a uniform grid with 1000 points or more.

characteristics). Others should in principle be available (e.g., oil property parameters) or could be measured in the laboratory (e.g., fouling model parameters). The oil physical properties, namely density, F, dynamic and kinematic viscosities, μ and ν, thermal conductivity, λ, and heat capacity, cp, can vary depending on the oil slate processed in the refinery. They are also temperature dependent (i.e., they vary across the heat exchanger). Such variations are sometimes substantial (e.g., see Figure 16) and need to be captured. Ideally, if the oil composition is known at all times, it can be used as an input to thermodynamic packages that can be interfaced with a process simulator to capture such a dependency. However, this is usually not the case and the properties of the fluids on both sides of the exchanger are typically unknown. An alternative is to estimate fluid physical properties and fouling model parameters from actual refinery data, and this section outlines a procedure for doing so. Established relationships60 based on just a few inputs, namely API, mean average boiling point, MeABP, and kinematic viscosity at 100 F, ν100 F, are used here to characterize the oils. These relationships reflect the main temperature dependency of oil density, viscosities, thermal conductivity, and heat capacity. If, as usually happens in practice, even these inputs are not known, these quantities are treated as a set of parameters to be estimated from plant data according to the procedure detailed in this section. This set of parameters, referred to as set A, includes the physical properties (i.e., API, MeABP, and ν100 F) of the fluids on both shell-side and tube-side. Another set of parameters, set B, includes activation energy, Ef, preexponential deposition constant, R, and suppression constant, γ, in the fouling model (eq 27). This set is also estimated from plant data for each shell comprised in the unit under investigation. The parameters are expected to vary from crude to crude but are also specific for the unit under consideration. 5.1. Parameter Estimation and Model Validation Procedure. First the model is fitted by adjusting the values of the parameters to one set of plant data and checked whether it is able to reproduce those data. Then, in order to test the predictive capabilities of the model, the parameters thus obtained are used to predict the outlet temperatures from the unit against a further set of plant data, different from the one used for the estimations.

ARTICLE

It is noted that parameters in set A do not depend on fouling. They can be estimated over a period of time when fouling has not yet started (or is very weak) and the exchanger can be considered clean. By contrast, parameter set B should be estimated over a longer period where fouling is indeed affecting the performance of the unit. The following periods are therefore defined: 1. Period I: Estimation period. This is divided into two subperiods: (a) Clean period. During this time the unit is considered clean. Data in this period are used to estimate values (called A0 ) of parameters in set A. (b) Operation period (immediately following period Ia). With A0 fixed, parameters in set B are estimated (to yield a value B0 ). 2. Period II: Prediction period (immediately following period I). The simulated responses are calculated with parameters A0 and B0 (now fixed) and may be compared with plant measurements. The assessment criteria are detailed in section 7.2. Crudes processed in a refinery change frequently (typically every ca. 3 days). Period Ia must therefore be sufficiently long so that fluid properties are not fitted to a single oil slate (which may not be representative) but are averaged out over the typical crudes processed by the refinery. Conversely, if the period is too long, the assumption that no fouling is present may become invalid. 5.2. Plant Data. The measured inlet temperatures and flow rates on the shell-side and tube-side are used as inputs to the model. The measured outlet temperatures on both shell-side and tube-side are used to fit the model parameters. These data are typically recorded at least on a daily basis, sometimes more frequently. Daily averages were used here. However, it is notoriously difficult to obtain reliable plant data;38,69 thus some judgment is required when analyzing such records. Crittenden et al.70 highlighted how fouling resistances calculated with refinery data after cleaning are nonzero, showing an offset which was attributed not only to the poor quality of plant data but also to the propagation of errors in the calculation of fouling resistances. In refinery practice, heat duties on both sides of the exchanger are calculated using a simple lumped model: Q ¼ mc _ p ΔT

ð45Þ

that assumes constant mass flow rates, fluid physical properties, and heat transfer coefficients. Refineries typically check the accuracy of their measurements by calculating the percentage difference of the heat duty, Q, on the two sides of the heat exchanger, here defined as error j: QS - QT % ð46Þ j ¼ QS Given the many assumptions described above, values of the error j = 10-30% are often found. Furthermore, such errors typically have a significant nonzero average offset and wide variability. Rigorous statistical methodologies exist in the literature to screen data, detect instrument bias, etc.71,72 Here it is assumed that refinery measurements do not suffer from bias and thus a simple procedure is sufficient to screen out grossly inaccurate measurements which would otherwise cause numerical problems in the estimation of model parameters and result in wrong parameter values. The procedure used calculates j from eq 46 (accepting it uses approximate heat balance models) from plant 4523


Industrial & Engineering Chemistry Research temperature and flow rate data, its average over time, jh, and standard deviation, σ. Daily data for which j deviates more than a certain value from the mean are deemed to be grossly erroneous or outliers and are discarded from consideration for parameterfitting purposes. The remaining data are defined as “filtered” and are used in a proper statistically sound estimation with the full detailed model. Such a filtering procedure can be applied in different ways depending on the quantity and quality of available raw data. In the authors’ experience, for the data-rich application considered here a conservative approach works best; i.e., it is better to reject possibly good data than to keep possibly bad ones, as long as a sufficient amount of data remains for the estimation step with the full model. Once highly suspect data have been eliminated, the next step is to estimate parameter sets A and B, using the filtered data set. First, data in period Ia are used to estimate parameter set A. The best estimated values are called set A0 . Then, with A0 fixed the parameter set B is fit to filtered data in period Ib resulting in set B0 . All parameters (sets A0 and B0 ) are then fixed and the full model may be used to simulate period II, in a fully predictive mode. Initial conditions of the exchanger for each period are those predicted by the model at the end of the previous period. Optimal estimates of model parameter sets A and B are obtained with the gPROMS in-built parameter estimation facility, which is based on a sophisticated, state-of-the-art implementation of the maximum likelihood formulation for dynamic systems.73,74 This method attempts to determine values for the uncertain physical model parameters (and measurement variances, if unknown) that maximize the probability that the mathematical model will predict the measurements available. This is achieved by minimizing the following (likelihood) function (subject to the full dynamic model), which assumes a normal probability density function for the variability: N Φ ¼ lnð2πÞ 2 8 2 39 2 = NV NM < NE i X ϑ X X ð~ a a Þ ijk 1 4lnðσ ijk 2 Þ þ ijk 5 þ min ð47Þ ; 2 ϑ :i ¼ 1 j ¼ 1 k ¼ 1 σ ijk 2 where N denotes the total number of measurements, ϑ is the set of parameters to be estimated (subject to given upper and lower bounds, i.e., ϑL e ϑ e ϑU), NE is the number of data sets available (in this case NE = 1), NV is the number of variables measured in the ith data set, NM is the number of measurements of the jth variable in the ith data set, σijk2 is the variance of the kth measurement of variable j in data set i, ~aijk is the kth measured value of variable j in experiment i, and aijk is the kth predicted value of variable j in data set i. Measurements from different experiments are assumed to be independent, but measurements from the same experiment are possibly correlated. The method also provides standard information on the statistical parameter significance and correlation of the parameters calculated, and model adequacy tests. For the estimation, parameters are scaled according to their respective initial guesses in order to avoid numerical problems. Depending on the instrumentation, a variety of measurement variance models may be used, from simple constant variance to full heteroschedastic. In PHT applications, the sensor used for the temperature measurement is typically a thermocouple for which a constant measurement variance of 1-2 C is appropriate.

ARTICLE

The model can be validated by comparing simulation results with plant data using overlay plots, by checking residuals and distribution curves of the residuals, and by performing a number of standard statistical checks.62 5.3. Hydraulics Related Parameters. Hydraulics effects of fouling are captured in eq 29. The foulant thickness is proportional to the thermal fouling resistance through the thermal conductivity of the fouling layer, λL. For a given value of the fouling resistance, Rf, a low λL value corresponds to a small value of thickness δ, and vice versa, affecting the cross-sectional area, the fluid velocity (eq 21), and, ultimately, the shear stress (eq 26) that govern the fouling suppression term in eq 27. The value of λL depends on the chemical composition of the organic deposits and changes over time according to the structural changes due to the exposure at high temperatures.75 Aging of the deposits can change their physical structure from a soft, gel-like substance to a hard coke-like one modifying the thermal conductivity up to 2 orders of magnitude. Few literature data are available for this parameter. To estimate it, pressure drop measurements would be required but, unfortunately, very often pressure data are not available in refineries. Thermal conductivity values15 used in the simulations in this paper were λ0L = 0.2 W m-1 K-1 (initial gel-like -1 -1 K (coked deposit), deposit, λ similar to oil) and λ¥ L =1Wm 62 as reported by Ishiyama et al.

6. CASE STUDY The model described was applied to simulate an industrial unit in the PHT of an ExxonMobil refinery. The unit comprises two identical shells with four tube passes per shell (the geometric parameters and tube inlet pressure are summarized in Table 2). The crude goes through the tube-side of the unit, whereas residuum from the crude distillation column flows on the shell-side (Figure 5). Plant personnel indicated that typical oils processed, although not light, produced severe fouling in this unit. Refinery observation when the unit was dismantled for cleaning indicated that, while some shell-side fouling was detected, the dominant fouling resistance was on the tube-side. Plant data (inlet temperatures and volumetric flow rates) over a period of 1 year after a mechanical cleaning were used for both parameter estimation and a check of simulations. Data points are available with a daily frequency and the inputs are considered to vary in a piecewise linear fashion between data points. The refinery's own calculations of fouling resistance are reported in Figure 7. Each point in this figure is calculated via the classic temperature difference logarithmic mean method, using daily average temperature and flow rate measurements reported in Figure 8 (the scatter in both is typical). An underlying linear trend for the fouling resistance in Figure 7 is clear, although in the period around 60-100 days of operations a deviation is observed before the substantially linear behavior is restored. The sharp drop in tube-side flow rates (Figure 8a) could possibly be responsible for this large temporary increase in the fouling rate. A typical large nonzero initial value of the fouling resistance calculated this way is evident, which should not be there as the exchanger is clean. After just less than a year of operations (347 days) the exchanger underwent another (chemical) cleaning. On this basis, a time horizon of 347 days after the mechanical cleaning was chosen for the study. The filtering procedure outlined in section 5 was used to check the quality of measurements. Figure 9 shows the value of j (the percentage difference between the shell-side heat duty and the tube-side heat duty) calculated via eq 46 for all 347 days. This figure also introduces a 4524



ARTICLE

Table 2. Heat Exchanger Geometric Parameters and Inlet Pressure parameter

symbol

units

value

tube-side inlet pressure

Pin

Pa

4 876 290

tube length outer tube diameter

L do

m mm

6.1 25.40

inner tube diameter

di

mm

19.05

shell diameter

DS

m

1397

tube count/shell

NT

880

pass number/shell

Np

4

Figure 7. Fouling resistance calculated with refinery procedure over a year after mechanical cleaning.

Figure 8. Inlet volumetric flow rates (a) and temperatures (b) for both shell-side and tube-side.

notation used for all the following figures: data are reported as an open symbol if retained by the filtering procedure and as a filled symbol if they are rejected as deemed unreliable. Also, two vertical dashed and dotted lines mark, respectively, the end of period Ia (22 days in this case study) and period I (60 days in this case study). In period I, 30% of the data points were rejected. In the following prediction period II (287 days or 9.5 months), 59 measurements (ca. 17%) were discarded, showing a reasonable quality of the data used to test the simulation model in prediction mode.

7. RESULTS AND DISCUSSION 7.1. Estimation Period. In all estimations, a constant measurement variance of 1 C was assumed for thermocouples. The parameter values obtained from the estimation procedure are reported in Table 3 for both parameter sets (A and B). The estimated crude API corresponds to a light crude, in line with indications by the refinery. The rather large fouling propensity of this oil;usually associated with heavier crudes;is reflected in the fouling parameters estimated in the subsequent step of the estimation procedure (set B). Most parameters within each set are strongly correlated; this was not unexpected, given the empirical nature of the underlying physical properties and fouling models. As a result, statistical tests for individual parameter estimates indicate they have less than 95% statistical confidence. However, the model overall shows an

Figure 9. Heat check, j, calculated via eq 46 from plant data, over a year after mechanical cleaning. Open symbols show data filtered with the procedure developed, whereas filled symbols show data deemed grossly inaccurate. The dashed vertical line indicates the end of period Ia, whereas the dotted line indicates the end of period I.

adequate statistical fit to the plant data, within the variance specified, according to the standard χ2 test (Table 4). These show that, 4525



ARTICLE

Table 3. Parameter Set A and B Estimates parameter set A parameter

parameter set B

optimal estimate

units

parameter

optimal estimate

residuum API

17.4

API

E01A R

0.001 65

residuum ν100 F

28.8

cSt

E01A E

28 491

736.7

C

E01A γ

9.28 10

crude API

37.5

API

E01B R

0.001 65

crude ν100 F

13.9

cSt

E01B E

28 523

350.2

C

E01B γ

K m W-1 s-1 J mol s-1

residuum MeABP

crude MeABP

units 2

-13

N m4 K W-1 s-1 K m 2 W-1 s-1 J mol s-1

9.35 10

-13

N m4 K W-1 s-1

Table 4. Weighted Residuals and Corresponding χ2 Values weighted residual

χ2 value (95%)

period Ia

10.902

21.026

period I

59.651

92.808

while individual parameters in the oil property and fouling/aging models are not to be believed (probably due to attempting to model a range of oil slates for the former and some underlying mismatch in the model physics for the latter), the overall fit is very good. This is attributed to the model representing well the dominant aspects of thermal, hydraulic, and fouling behaviors in an integrated way. Parts a and b of Figure 10 show overlay plots of model simulation results with, respectively, measured crude and residuum outlet temperatures for the estimation period (period I). Each data point is reported with an error bar of (1 C corresponding to the (assumed) variance of the temperature measurements. The model simulations fit well the measured plant data for all the data points for both shell-side and tube-side. In particular, for the tube (crude) side, all the points deemed “reliable” by the filtering procedure are fitted within (1 C. It is noteworthy that several of the data points which were filtered out and not used for parameter estimation turn out to be well simulated as well. 7.2. Prediction Period. To test whether the model can be used in a predictive mode (i.e., to extrapolate its results beyond the estimation period, without reestimating its parameters), the model was used to simulate the exchanger for the rest of the year. Measured inlet temperatures and flow rates were input to the simulation, but all adjustable parameters (sets A and B) were kept fixed. The simulated performance, in terms of output temperatures over the whole time horizon considered (347 days), is shown by the continuous lines in Figure 11a and 11b for, respectively, tube-side and shell-side, together with plant measurements. The exit temperatures on both sides are predicted well over the entire period, even in the presence of rather major excursions, including the period between 60 and 100 days. The contribution of fouling and aging to the overall behavior of the unit is illustrated by the two extra lines plotted in Figure 11a and 11b: (i) The dashed line was obtained using the same parameter set A as before but setting the fouling rate to 0 (i.e., dRf/dt = 0), corresponding to a case of no fouling. (ii) The dashed-dotted line was obtained using the same parameter sets A and B but setting the preexponential term, Aa, in the aging model to 0, corresponding to a case of fouling but with no deposit aging. In particular, from Figure 11a an insight can be gained on the opposite effects of fouling and aging on the crude outlet temperature. While fouling reduces the efficiency of the heat transfer thus

Figure 10. Overlay plot of model simulations with crude outlet temperature (a) and shell-side outlet temperature (b) over the estimation period I.

decreasing the crude outlet temperature, aging acts in the opposite way, by increasing the thermal conductivity of the deposits, thus enhancing the overall heat transfer. Temperatures estimated without fouling or aging in the model do not match the plant data, but they do so when both phenomena are included. A closer analysis of the accuracy of simulation results is made by considering the percentage residuals defined as ε ¼ 4526

~xi - xi % ~xi

ð48Þ



ARTICLE

Figure 11. Overlay plot of model simulations with crude outlet temperature (a) shell-side outlet temperature (b). Dashed line shows simulation run setting Rf = 0 whereas the dashed-dotted line shows simulation results with no ageing (i.e. Aa = 0). White markers show data filtered with the procedure developed whereas colored markers show data deemed inaccurate. The error bars indicate (1 C. The vertical dashed line shows the end of period Ia whereas the dotted line indicates the end of period I and the beginning of period II.

where xi is the simulated value of the variable considered at the ith time step and ~xi is the corresponding plant measurement. Residuals of all data previously deemed “accurate“ on the tube side (Figure 12) lay between -1.2 and þ0.6% (corresponding to absolute error values of -2.5 and þ1.5 C) over the 347 days of operation. Shell-side accuracy is somewhat lower as shown in Figure 12. In this case the residuals of all reliable data lay between -1 and þ2% (-3 and þ6 C). The residuals show some systematic trend, as opposed to purely random scatter, indicating there is possibly some underlying model mismatch (not surprising, for example, in view of likely oil slate changes over 1 year and shell-side assumptions). Even so, analysis of the parity plots (Figure 13) reveals that although the model tends to underoverestimate the shell-side temperature, over 85% of the points are within a (1% error for both sides of the unit.

The effect of aging was investigated in Figure 11a and 11b by simply re-running simulations with the aging parameter Aa set to zero. However, the other fouling parameters used resulted from an estimation with the aging model included. To investigate whether this could mask the effect of aging (and indeed whether it is necessary to include the aging phenomena at all), the parameter estimation procedure was performed again in period Ib with Aa set to zero. The results, reported in Figure 14 in the form of percentage residuals in the output temperatures, show that the error increases significantly. Moreover, a strong trend is introduced in the residuals, highlighting the structural importance of including aging effects. The contribution of each of the two shells, E01A and E01B, to the overall average fouling resistance Rhf, calculated from eq 31, is shown in Figure 15. As expected for a clean exchanger, the value 4527



ARTICLE

Figure 13. Parity plot for tube-side (a) and shell-side (b) outlet temperatures. The dashed lines indicate (1% deviation. Figure 12. Percentage residuals for tube-side (a) and shell-side (b) outlet temperatures. The dashed vertical line shows the end of period Ia, whereas the dotted line indicates the end of period I.

of the fouling resistance increases from an initial zero value. All variability present in the inlet temperatures and flow rates is smoothed out, to give a smooth curve. The rate of increase is high initially but tails off after about 150 days to an approximately constant value. Fouling in shell E01A is slightly higher than that in shell E01B, but the difference (in this case) is not large. Model outputs for Rhf are also compared with refinery calculations. The latter show a zero fouling resistance for the first few days, after which a sharp jump to an offset value of ca. 2 10-3 m2 K W-1 occurs. The model results start from a zero fouling resistance for a clean heat exchanger. Given the high accuracy of the exit temperatures predicted by the model, it is likely that the overall fouling resistance calculated by the refinery merely reflects the gross approximations (e.g., in physical properties) used in its calculation. The difference between model and refinery-calculated fouling resistances in Figure 15 is not considered large by industrial refinery specialists, who acknowledge the approximate nature of their calculations and mostly consider just underlying trends anyway, not absolute values, and would typically subtract from their curve the cleanexchanger offset. The overall trends shown substantially agree, and when the offset is subtracted, numerical values are also close.

Figure 14. Shell-side and tube-side percentage residuals when estimating parameters in period Ib in case of no aging.

In particular, Figure 15 shows that the (perceived) sudden increase in fouling resistance around 80-100 days can be fully explained by the measured changes in inlet flow rates, temperatures, and complex interactions arising. 4528


Industrial & Engineering Chemistry Research In the authors’ view, the reason for the quality of model temperature predictions shown in previous figures is that the model captures the essential physics in an appropriate way. It not only includes fouling and aging (Figure 11a and 11b), it also considers the effects of variation in physical properties over time and space (as a function of temperature) on the heat transfer

Figure 15. Refinery-calculated fouling resistance and model-calculated overall average Rf. The dotted and dashed lines indicate average Rf values calculated for E01A and E01B, respectively.

ARTICLE

coefficient and key interactions between hydraulic and thermal aspects. Figure 16 shows that the variations of physical properties over the length of each shell in the unit and over time are, indeed, significant. One of the main benefits of using a distributed model is that it is possible to track the effects of fouling along the heat exchanger over time. Figure 17 shows the fouling layer thickness along the

Figure 17. Fouling layer thickness along the two shells at the end of period II (347 days after cleaning). Arrows indicate the direction of the crude flow. Shell-side flow direction is positive for shell E01A and negative for shell E010B.

Figure 16. Variation over time of fluid density (a), specific heat capacity (b), dynamic viscosity (c), and thermal conductivity (d) at shell-side and tubeside outlets with respect to inlet conditions, for each shell in unit E01. 4529



Figure 18. Tube-side pressure drops across the two shells E01A and E01B, and total pressure drop across the unit.

four passes of each shell at the end of period II (347 days after the previous cleaning). A first observation is that in the countercurrent passes, owing to a larger film temperature, the increase in fouling thickness along the tubes is much larger compared to the passes in parallel flow. This difference is particularly evident in the colder shell (E01A), where the driving force (i.e., temperature difference) is larger. Figure 17 corroborates the speculations by Polley et al.33 regarding the effect of flow arrangement on fouling. The model predicts a larger thickness of the fouling layer at the tube outlet of shell E01A than at the inlet of shell E01B. This is counterintuitive as the latter shell is the one at higher crude bulk temperature. However, the film temperature (on which fouling depends according to eq 27) is, on average, more than 2.5 C higher at the inlet of shell E01B than the outlet of E01A. It is noted that the fouling deposits accumulated over the 347 days of operations amount to ca. 3 mm, which corresponds to a large reduction (over 40%) in crosssectional area available to the crude flow and has a strong impact on the unit’s hydraulics. It also confirms that the thinlayer assumption often used in fouling models cannot be applied. It is possible to appreciate the effects of fouling on the hydraulics of the exchangers in Figure 18, where the pressure drops are reported for the tube-side of each of the two shells as dotted (E01A) and dashed lines (E01B), respectively. The total pressure drop, the sum of the two, is also shown. The irregular shape of the curves largely depends on the variability of flow rates, but the underlying trends reflect the effects of fouling. As the two shells are geometrically identical, for the first 130 days after the cleaning the pressure drops are the same whereas, when fouling becomes significant, the contribution to the total pressure drops given by E01B becomes more important than that given by E01A. The overall increasing trend of overall pressure drop in Figure 18 shows the large impact of fouling on the unit. The hydraulic model could not be validated against plant data as pressure drop measurements are usually not logged by refineries and are virtually unavailable (they were not in this case study). Nonetheless, the capability of predicting pressure drops represents another important feature of the model that can assist refiners in monitoring the performance of the PHT.

ARTICLE

8. CONCLUSIONS AND FUTURE WORK A novel dynamic, distributed model for a multipass shell-andtube heat exchanger was developed for use in refinery preheat trains, where fouling deposition results in severe energy inefficiency. The model captures the variation of physical properties with temperature, over time and space. It considers several aspects of the fouling phenomenon such as its dependence on exchanger geometry, process conditions, the fouling propensity of the crude, and a deposit aging mechanism. The formulation in cylindrical coordinates makes it possible to overcome the thin slab approximation;often used in the past;accounting for curvature effects in the heat flux. Interactions between the growth of the fouling layer and the crude oil flowing inside the tubes are captured by solving a moving boundary problem. The model enables a detailed analysis to be carried out of the closely coupled spatial and temporal effects of fouling on the thermal and hydraulic performance of a heat exchanger unit. Specific geometric configurations are represented, contributions of distinct shells in the same unit can be differentiated, and critical zones are identified. This goes well beyond current practice, which relies on the use of very aggregate, averaged, and highly empirical fouling factors. A formal validation of the model has been performed for a rather difficult heat exchanger in a refinery site, despite the well-known difficulties in collecting and analyzing industrial data on fouling behavior. The case study has shown that a 2 month estimation period is sufficient to generate parameters yielding an excellent representation of the outlet temperatures during the estimation period itself. Furthermore, the model with parameters thus adjusted could be used in a fully predictive way for another 10 months to predict exit temperatures to an accuracy within 1% on the tube-side and 2% on the shell-side. Primary plant measurements, rather than derived fouling resistances, were used to estimate the model parameters, thus avoiding incorporating into the analysis errors resulting from simplifying assumptions in fouling resistance calculations. A single case study was presented here, and clearly more extensive demonstrations are needed. However, some further industrial exchanger units from different refineries with different crudes and over different temperature ranges have been recently analyzed utilizing the approach described.76 While this experience shows that a careful analysis of the raw data is always required, it also confirms that fitting accuracies and predictive ability comparable to those presented here are obtained in rather diverse situations. In the model presented, fouling is limited to the tube-side. The results presented here and those reported by Coletti and Macchietto76 indicate that (at least for the cases considered) the potential underprediction of shell-side heat transfer resistance does not have a significant effect. This may not always be true. Indeed, a large mismatch between this model and plant measurements can indicate (indirectly) that shell-side fouling may be important. Of course, in the presence of very large shellside fouling rates the model may not be applicable, and for generality, it would be useful to also include shell-side fouling. A validation of the hydraulic part of the model awaits the availability of pressure drop measurements. It is envisaged that the benefits of using the approach proposed will be 2-fold. First, the improved information resulting from the careful combination of the advanced model presented and 4530


Industrial & Engineering Chemistry Research primary data may be used in refineries to assist with operational aspects: • assessment of heat exchanger performance, using historical data • online model-based monitoring of significant events impacting fouling performance (e.g., early detection of problems due a change in crude slate) • assessment of energy losses, thus extra costs, due to fouling77,78 • assisting in operation decisions such as cleaning scheduling The second are of potential use stemming from the ability to incorporate fouling in the design of heat exchangers, to explore, for example the following: • retrofit options to mitigate fouling77,79 • economic impact of a chosen design on refinery PHT costs and CO2 emissions80

’ AUTHOR INFORMATION

ARTICLE

R = radius, m r = radial coordinate, m ~r = dimensionless radial coordinate, dimensionless Re = Reynolds number, dimensionless Rf = fouling thermal resistance, W m-2 K-1 Rhf = average fouling resistance, W m-2 K-1 Rflow = flow radius, m Rg = universal gas constant, J kg-1 mol-1 S = heat transfer surface, m2 T = temperature, K Tf = film temperature, K t = time, s tage = aging time, s u = mean fluid velocity in the tube, m s-1 y = youth variable · V_ = volumetric flow rate, m3 h-1 x = simulated value h~x = measured value z = axial coordinate, m

Corresponding Author

*E-mail: [email protected].

’ ACKNOWLEDGMENT The authors gratefully acknowledge ExxonMobil for the data provided, the Engineering and Physical Sciences Research Council U.K. (EPSRC) for funding (Grant EP/D503051/1), and industrial members of the IHS ESDU Oil Industry Fouling Working Party for useful discussions, encouragement, and feedback. ’ NOTATION A = area, m2 Aa = frequency factor of the aging model, s-1 Aflow = cross-sectional area for flow, m2 a = predicted value of variable ~a = measured value of variable Cf = friction factor, dimensionless cp = specific heat capacity at constant pressure, J kg-1 K-1 DS = shell inner diameter, m d = diameter, m dir = direction of flow Ea = activation energy for aging model, J mol-1 Ef = activation energy for fouling model, J mol-1 h = local heat transfer coefficient, W m-2 K-1 hid = heat transfer coefficient for ideal cross-flow, W m-2 K-1 J = heat transfer coefficient correction factor L = length, m · m_ = mass flow rate, kg s-1 Pr = Prandtl number, dimensionless N = total number of measurements NE = number of data sets NM = number of measurements in data set Np = number of tube passes per shell NS = number of shells per unit NT = number of tubes NV = number of variables in data set Nu = Nusselt number, dimensionless P = pressure, wetted perimeter, m Q = heat duty, W q00 = heat flux, W m-2 q000 = heat source, W m-3 Rg = universal gas constant, J kg-1 mol-1

Subscripts

f = fouling, film i = inner in = inlet L = fouling layer n = pass number o = outer out = outlet S = shell-side T = tube-side W = wall Superscripts

0 = initial state of deposit ¥ = final state of deposit Greek Symbols

R = deposition parameter, m2 K J-1 γ = suppression parameter, m4 K N-1 J-1 δ = foulant thickness, m ε = temperature percentage residuals ζ = filtering procedure parameter, % ϑ = set of parameters to be estimated λ = thermal conductivity, W m-1 K-1 μ = dynamic viscosity, Pa s ν = kinematic viscosity, St F = density, kg m-3 σ = measurement variance, C τ = shear stress, Pa j = heat balance closure parameter, % Ω = model domain

’ REFERENCES (1) ESDU. Heat Exchanger Fouling in the Pre-Heat Train of a Crude Oil Distillation Unit, Data Item 00016; ESDU International PLC: London, 2000. (2) Yeap, B. L.; Wilson, D. I.; Polley, G. T.; Pugh, S. J. Mitigation of Crude Oil Refinery Heat Exchanger Fouling through Retrofits Based on Thermo-Hydraulic Fouling Models. Chem. Eng. Res. Des. 2004, 82A, 53–71. (3) Andersson, E.; Quah, J.; Polley, G. T. Experience in Application of Compabloc Heat Exchangers in Refinery Pre-Heat Trains. In International Conference on Heat Exchanger Fouling and Cleaning VIII, Schladming, Austria, June 14-19, 2009 [Online]; M€uller-Steinhagen, H., 4531


Industrial & Engineering Chemistry Research Malayeri, R., Watkinson, P., Eds.; 2009; pp 39-43. http://heatexchangerfouling.com/papers/papers2009/6_Polley_Compabloc_F.pdf. (4) Panchal, C. B.; Huangfu, E.-P. Effects of Mitigating Fouling on the Energy Efficiency of Crude-Oil Distillation. Heat Transfer Eng. 2000, 21 (3), 3–9. (5) Liebmann, K.; Dhole, V. R. Integrated Crude Distillation Design. Comput. Chem. Eng. 1995, 19 (Suppl. 1), 119–124. (6) Papalexandri, K. P.; Patsiatzis, D. I.; Pistikopoulos, E. N.; Ebbesen, L. Heat Integration Aspects in a Crude Preheat Refinery Section. Comput. Chem. Eng. 1998, 22 (Suppl. 1), S141–S148. (7) Bagajewicz, M.; Soto, J. Rigorous Procedure for the Design of Conventional Atmospheric Crude Fractionation Units. Part II: Heat Exchanger Network. Ind. Eng. Chem. Res. 2000, 40 (2), 627–634. (8) TEMA. Standards of the Tubular Exchangers Manufacturers Association; TEMA: New York, 1968; pp 123-127. (9) Chenoweth, J. M. The TEMA Standards Fouling Section: Another Look. In Fouling Mitigation of Industrial Heat-Exchange Equipment; Panchal, C. B., Ed.; Begell House: San Luis Obispo, CA, 1995; pp 61-71. (10) Rabas, T. J.; Panchal, C. B. Fouling Rates, Not Fouling Resistances. Heat Transfer Eng. 2000, 21 (2), 1–3. (11) Bennett, C. A.; Kistler, R. S.; Lestina, T. G.; King, D. C. Improving Heat Exchanger Designs. Chem. Eng. Prog. 2007, 103 (4), 40–45. (12) Liebmann, K.; Dhole, V. R.; Jobson, M. Integrated Design of a Conventional Crude Oil Distillation Tower Using Pinch Analysis. Chem. Eng. Res. Des. 1998, 76 (3), 335–347. (13) Wilson, D. I.; Polley, G. T.; Pugh, S. J. Mitigation of Crude Oil Preheat Train Fouling by Design. Heat Transfer Eng. 2002, 23 (1), 24–37. (14) Epstein, N. Thinking About Heat Transfer Fouling: A 5 5 Matrix. Heat Transfer Eng. 1983, 4 (1), 43–56. (15) Watkinson, P. Critical Review of Organic Fluid Fouling: Final Report; Argonne National Laboratory: Argonne, IL, 1988. (16) Perry, R. H.; Green, D. W. Perry’s Chemical Engineers’ Handbook, 7th ed.; McGraw-Hill: New York, 1997. (17) BP. Personal communication, London, 2006. (18) M€uller-Steinhagen, H.; Malayeri, M. R.; Watkinson, A. P. Heat Exchanger Fouling: Environmental Impacts. Heat Transfer Eng. 2009, 30 (10), 773–776. (19) DOE. Energy Bandwidth for Petroleum Refining Processes; U.S. Department of Energy, Office of Energy Efficiency and Renewable Energy Industrial Technologies Program, 2006. (20) Van Nostrand, W. L.; Leach, S. H.; Haluska, J. L. Economic Penalties Associated with the Fouling of Refinery Heat Transfer Equipment. In Fouling of Heat Transfer Equipment; Somerscales, E. F. C., Knudsen, J. G., Eds.; Hemisphere: Troy, NY, 1981; pp 619-643. (21) Taborek, J.; Aoki, T.; Ritter, R. B.; Palen, J. W.; Knudsen, J. G. Fouling: The Major Unresolved Problem in Heat Transfer. Chem. Eng. Prog. 1972, 68 (2), 59–67. (22) Jeronimo, M. A. S.; Melo, L. F.; Sousa Braga, A.; Ferreira, P. J. B. F.; Martins, C. Monitoring the Thermal Efficiency of Fouled Heat Exchangers: A Simplified Method. Exp. Therm. Fluid Sci. 1997, 14 (4), 455–463. (23) Liporace, F. S.; De Oliveira, S. G. Real Time Fouling Diagnosis and Heat Exchanger Performance. Heat Transfer Eng. 2007, 28 (3), 193– 201. (24) 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 (14), 3924–3938. (25) Watkinson, A. P.; Wilson, D. I. Chemical Reaction Fouling: A Review. Exp. Therm. Fluid Sci. 1997, 14 (4), 361–374. (26) Watkinson, A. P. Deposition from Crude Oils in Heat Exchangers. Heat Transfer Eng. 2007, 28 (3), 177–184. (27) Bennett, C. A.; Kistler, R. S.; Nangia, K.; Al-Ghawas, W.; Al-Hajji, N.; Al-Jemaz, A. Observation of an Isokinetic Temperature and Compensation Effect for High-Temperature Crude Oil Fouling. Heat Transfer Eng. 2009, 30 (10), 794–804.

ARTICLE

(28) Kern, D. Q.; Seaton, R. E. A Theoretical Analysis of Thermal Surface Fouling. Br. Chem. Eng. 1959, 4 (5), 258–262. (29) Epstein, N. General Thermal Fouling Models. In Fouling Science and Technology; Melo, L. F., Bott, T. R., Bernardo, C. A., Eds.; Kluver Academic Publishers: Alvor, Portugal, 1987; pp 15-30. (30) Crittenden, B. D.; Kolaczkowski, S. T.; Hout, S. A. Modelling Hydrocarbon Fouling. Chem. Eng. Res. Des. 1987, 65 (2), 171–179. (31) Knudsen, J. G.; Dahcheng, L.; Ebert, W. A. The Determination of the Threshold Fouling Curve for a Crude Oil. In Understanding Heat Exchanger Fouling and its Mitigation; Bott, T. R., Melo, L. F., Panchal, C. B., Somerscales, E. F. C., Eds.; Begell House: Lucca, Italy, 1997; pp 265-272. (32) Ebert, W. A.; Panchal, C. B. Analysis of Exxon Crude-Oil-Slip Stream Coking Data. In Fouling Mitigation of Industrial Heat-Exchange Equipment; Panchal, C. B., Ed.; Begell House: San Luis Obispo, CA, 1995; pp 451-460. (33) 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 (7), 763–776. (34) Butterworth, D. Design of Shell-and-Tube Heat Exchangers When the Fouling Depends on Local Temperature and Velocity. Appl. Therm. Eng. 2002, 22 (7), 789–801. (35) Nasr, M. R. J.; Givi, M. M. Modeling of Crude Oil Fouling in Preheat Exchangers of Refinery Distillation Units. Appl. Therm. Eng. 2006, 26 (14-15), 1572–1577. (36) Bories, M.; Patureaux, T. Preheat Train Crude Distillation Fouling Propensity Evaluation by the Ebert and Panchal Model. In ECI Conference on Heat Exchanger Fouling and Cleaning: Fundamentals and Applications; Watkinson, P., M€uller-Steinhagen, H., Malayeri, M. R., Eds.; Berkeley Electronic Press: Santa Fe, NM, 2003; pp 200-210. (37) Smaili, F.; Vassiliadis, V. S.; Wilson, D. I. Mitigation of Fouling in Refinery Heat Exchanger Networks by Optimal Management of Cleaning. Energy Fuels 2001, 15 (5), 1038–1056. (38) Wilson, D. I.; Vassiliadis, V. S. Mitigation of Refinery Fouling by Management of Cleaning. In Understanding Heat Exchanger Fouling and Its Mitigation; Bott, T. R., Melo, L. F., Panchal, C. B., Somerscales, E. F. C., Eds.; Begell House: Lucca, Italy, 1997; pp 299-306. (39) Ishiyama, E. M.; Paterson, W. R.; Wilson, D. I. The Effect of Fouling on Heat Transfer, Pressure Drop and Throughput in Refinery Preheat Trains: Optimisation of Cleaning Schedules. In Fouling and Cleaning of Heat Exchangers VII, Tomar, Portugal, July 1-6, 2007; M€uller-Steinhagen, H.; Malayeri, R.; Watkinson, P., Eds.; 2007; http://services.bepress.com/eci/heatexchanger2007/9/. (40) Gardner, K.; Taborek, J. Mean Temperature Difference: A Reappraisal. AIChE J. 1977, 23 (6), 777–786. (41) Yeap, B. L.; Hewitt, G. F.; Wilson, D. I.; Pugh, S. J. Simulating the Fouling Behaviour in Shell-and-Tube Heat Exchangers. In 4th International Conference on Heat Exchanger Fouling. Fundamental Approaches and Technical Solutions. Publico Publications: Germany, 2002; pp 275-282. (42) Tan, K. S.; Spinner, I. H. Dynamics of a Shell-and-Tube Heat Exchanger with Finite Tube-Wall Heat Capacity and Finite Shell-Side Resistance. Ind. Eng. Chem Fundam. 1978, 17 (4), 353–358. (43) Yin, J.; Jensen, M. K. Analytic Model for Transient Heat Exchanger Response. Int. J. Heat Mass Transfer 2003, 46 (17), 3255– 3264. (44) Roetzel, W.; Xuan, Y. Transient-Response of Parallel and Counterflow Heat-Exchangers. J. Heat Transfer 1992, 114 (2), 510–512. (45) Xuan, Y.; Roetzel, W. Dynamics of Shell-and-Tube Heat Exchangers to Arbitrary Temperature and Step Flow Variations. AIChE J. 1993, 39 (3), 413–421. (46) Pataknar, S. V.; Spalding, D. B. A Calculation Procedure for the Transient and Steady-State Behavior of a Shell-and-Tube Heat Exchangers. In Heat Exchangers;Design and Theory Sourcebook; Afgan, N. H., Schlunder, E. V., Eds.; McGraw-Hill: New York, 1974; pp 155-176. (47) Prithiviraj, M.; Andrews, M. J. Three Dimensional Numerical Simulation of Shell-and-Tube Heat Exchangers. Part I: Foundation and Fluid Mechanics. Numer. Heat Transfer, Part A 1998, 33 (8), 799–816. 4532


Industrial & Engineering Chemistry Research (48) Kukral, R.; Stephan, K. Dynamic Simulation of Shell-and-Tube Heat Exchangers Based on Models of Intermediate Complexity. In Third U.K. National Conference Incorporating 1st European Conference on Thermal Sciences; Institution of Chemical Engineers/Hemisphere: Birmingham, England, 1992; pp 417-423. (49) Taborek, J. Shell-and-Tube Heat Exchangers: Single-Phase Flow. In Heat Exchanger Design Handbook; Hewitt, G. F., Ed.; Begell House: New York, 2002; Vol. 3. (50) Roetzel, W.; Xuan, Y. Dynamic Behaviour of Heat Exchangers; Computational Mechanics Publications 3; WIT Press: Boston, 1999. (51) Fryer, P. J.; Slater, N. K. H. A Direct Simulation Procedure for Chemical Reaction Fouling in Heat Exchangers. Chem. Eng. J. 1985, 31 (2), 97–107. (52) Toyoda, I.; Fryer, P. J. A Computational Model for Reaction and Mass Transfer in Fouling from Whey Protein Solutions. In Fouling Mitigation of Industrial Heat-Exchange Equipment; Panchal, C. B., Ed.; Begell House: San Luis Obispo, CA, 1997; pp 397-406. (53) Georgiadis, M. C.; Rotstein, G. E.; Macchietto, S. Modeling and Simulation of Shell and Tube Heat Exchangers under Milk Fouling. AIChE J. 1998, 44 (4), 959–971. (54) Georgiadis, M. C.; Macchietto, S. Dynamic Modelling and Simulation of Plate Heat Exchangers under Milk Fouling. Chem. Eng. Sci. 2000, 55 (9), 1605–1619. (55) Georgiadis, M. C.; Rotstein, G. E.; Macchietto, S. Optimal Design and Operation of Heat Exchangers under Milk Fouling. AIChE J. 1998, 44 (9), 2099–2110. (56) Georgiadis, M. C.; Rotstein, G. E.; Macchietto, S. Modelling and Simulation of Complex Plate Heat Exchanger Arrangements under Milk Fouling. Comput. Chem. Eng. 1998, 22 (Suppl. 1), S331–S338. (57) Clarke, R. H.; Nicolas, F. CFD Investigation of Maldistribution Effects on Crude Oil Fouling in Shell and Tube Exchangers. Presented at the Second International Conference on Petroleum and Gas Phase Behaviour and Fouling, Copenhagen, Denmark, August, 2000. (58) Radhakrishnan, V. R.; Ramasamy, M.; Zabiri, H.; Do Thanh, V.; Tahir, N. M.; Mukhtar, H.; Hamdi, M. R.; Ramli, N. Heat Exchanger Fouling Model and Preventive Maintenance Scheduling Tool. Appl. Therm. Eng. 2007, 27 (17-18), 2791–2802. (59) Aminian, J.; Shahhosseini, S. Neuro-Based Formulation to Predict Fouling Threshold in Crude Preheaters. Int. Commun. Heat Mass Transfer 2009, 36 (5), 525–531. (60) Riazi, M. R. Characterization and Properties of Petroleum Fractions, 1st ed.; ASTM: Philadelphia, 2005. (61) Bejan, A. Convection Heat Transfer; Wiley & Sons: New York, 1984. (62) Ishiyama, E. M.; Coletti, F.; Macchietto, S.; Paterson, W. R.; Wilson, D. I. Impact of Deposit Ageing on Thermal Fouling: Lumped Parameter Model. AIChE J. 2010, 56 (2), 531–545. (63) Coletti, F.; Ishiyama, E. M.; Paterson, W. R.; Wilson, D. I.; Macchietto, S. Impact of Deposit Ageing and Surface Roughness on Thermal Fouling: Distributed Model. AIChE J. 2010, 56 (12), 3257– 3273. (64) Hewitt, G. F.; Shires, G. L.; Bott, T. R. Process Heat Transfer; CRC Press: London, 1994. (65) Wilkes, J. O. Fluid Mechanics for Chemical Engineers with Microfluidics and CFD, 2nd ed.; Prentice Hall: Upper Saddle River, NJ, 2005. (66) 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 Its Mitigation; Bott, T. R., Ed.; Begell House: Lucca, Italy, 1997; pp 273-281. (67) Process Systems Enterprise. gPROMS. www.psenterprise. com/gproms; 1997-2009. (68) CAPE-OPEN Open Interface Specification: Thermodynamic and Physical Properties. Available at http://www.colan.org/index-3.html. (69) Takemoto, T.; Crittenden, B. D.; Kolaczkowski, S. T. Interpretation of Fouling Data in Industrial Shell and Tube Heat Exchangers. Chem. Eng. Res. Des. 1999, 77 (8), 769–778. (70) Crittenden, B. D.; Kolaczkowski, S. T.; Downey, I. L. Fouling of Crude Oil Preheat Exchanger. Trans. Inst. Chem. Eng. 1992, 70A, 547–557.

ARTICLE

(71) Bagajewicz, M. Data Reconciliation and Instrumentation Upgrade. Overview and Challenges. Presented at Foundations of Computer Aided Process Operations (FOCAPO), Coral Springs, FL, January, 2003. (72) Bagajewicz, M. Process Plant Instrumentation: Design and Upgrade, 1st ed.; CRC Press: Boca Raton, FL, 2000. (73) Parameter Estimation in gPROMS. In gPROMS Advanced User Guide 3.2; Process Systems Enterprise: London, 2009. (74) Bard, Y. Nonlinear Parameter Estimation; Academic Press: London, 1974. (75) Fan, Z.; Watkinson, A. P. Aging of Carbonaceous Deposits from Heavy Hydrocarbon Vapors. Ind. Eng. Chem. Res. 2006, 45 (18), 6104–6110. (76) Coletti, F.; Macchietto, S. Validation of a Novel Model for Shell-and-Tube Heat Exchangers Undergoing Crude Oil Fouling. In 14th International Heat Transfer Conference (IHTC14); ASME: Washington, DC, 2010. (77) Coletti, F.; Macchietto, S. A Heat Exchanger Model to Increase Energy Efficiency in Refinery Pre-Heat Trains. In European Symposium on Computer Aided Process Engineering (ESCAPE) 19; Jezowski, J., Thullie, J., Eds.; Elsevier: Cracow, Poland, 2009; pp 1245-1250. (78) Coletti, F.; Macchietto, S. Predicting Refinery Energy Losses Due to Fouling in Heat Exchangers. In 10th International Symposium on Process Systems Engineering (PSE); de Brito Alves, R. M., Oller do Nascimento, C. A., Chalbaud Biscaia, E., Jr., Eds.; Elsevier: Salvador de Bahia, Brazil, 2009; pp 219-224. (79) Coletti, F.; Macchietto, S.; Polley, G. T. Effects of Fouling on Performance of Retrofitted Heat Exchanger Networks; a ThermoHydraulic Based Analysis. Pierucci, S., Buzzi Ferraris, G., Eds.; In European Symposium on Computer Aided Process Engineering (ESCAPE) 20; Ischia, Naples, Italy, June 6-9, 2010; pp 19-24. (80) Coletti, F.; Macchietto, S. Refinery Pre-Heat Train Network Simulation Undergoing Fouling: Assessment of Energy Efficiency and Carbon Emissions. Heat Transfer Eng. 2011, 32 (3 & 4), 228–236. (81) Macchietto, S.; Hewitt, G. F.; Coletti, F.; Crittenden, B. D.; Dugwell, D. R.; Galindo, A.; Jackson, G.; Kandiyoti, R.; Kazarian, S. G.; Luckham, P. F.; Matar, O. K.; Millan-Agorio, M.; M€uller, E. A.; Paterson, W.; Pugh, S. J.; Richardson, S.M.; Wilson, D. I. Fouling in Crude Oil Preheat Trains: A Systematic Solution to an Old Problem. Heat Transfer Eng. 32 (3-4) 197-215. (82) Baudelet, C. A.; Krueger, A. W. A Practical Approach to Fouling Mitigation in Refinery Units: The Spirelf. Presented at the AIChE Spring National Meeting, New Orleans, March 8-12, 1998.

4533


A Dynamic, Distributed Model of Shell-and-Tube Heat Exchangers

Recommend Documents