**Robust Control of Distributed Parameter Systems with Demonstration in Casting Technology and MATLAB/Simulink/DPS Blockset Software Support**

Cyril Belavý, Gabriel Hulkó and Karol Ondrejkovič

Additional information is available at the end of the chapter

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

## **1. Introduction**

26 Will-be-set-by-IN-TECH

No. 8, pp. 967-978.

pp. 626–638.

1059–1060.

16–23.

No. 7, pp. 262–278.

7-10, Budapest-Hungary, pp. 341–344.

[8] El Aroudi, A. & Leyva, R. (2001). Quasi-periodic route to chaos in a PWM voltage-controlled DC-DC boost converter. *IEEE Trans. on Circuits and Systems*, Vol. 48,

[9] El Aroudi, A.; Debbat, M.; Giral, R.; Olivar, G.; Benadero, L.& Toribio, E. (2005). Bifurcations in DC-DC switching converters: review of methods and applications.

[10] El Guezar, F. & Bouzahir, H. (2008). Chaotic behavior in a switched dynamical system. *Modelling and Simulation in Engineering*, Vol. 2008, Article ID 798395, 6 pages. [11] El Guezar, F.; Acco, P.; Bouzahir, H. & Fournier-Prunaret, D. (2008). Accurate and Fast Event Detection Occurrence in Planar Piecewise Affine Hybrid Systems. *Proceedings of the International Symposium NOLTA (NOn Linear Theory and its Applications)*, September

[12] El Guezar, F. (December 2009). Modélisation et simulation des systèmes dynamiques hybrides affines par morceaux. Exemples en électronique de puissance. *Ph D thesis*,

[13] El Guezar, F.; Bouzahir, H. & Fournier-Prunaret, D. (2011). Event Detection Occurrence For Planar Piecewise Affine Hybrid Systems. *Nonlinear Analysis: Hybrid Systems*, Vol. 5,

[14] Girard, A. (2002). Detection of event occurrence in piece-wise linear hybrid systems. *Proceedings of the 4th International Conference on Recent Advances in Soft Computing*,

[15] Girard, A. (September 2004). Analyse algorithmique des systèmes hybrides. *Ph D thesis*,

[16] Hamill , D. C. & Jeffries, D. J. (1988). Subharmonics and chaos in a controlled switched-mode power converter. *IEEE Trans. Circuits Systs. I*, Vol. 35, No. 8, pp.

[17] Hedayat, C. D.; Hachem, A.; Leduc, Y. & Benbassat, G. (March–April 1997). High-level modeling applied to the second-order charge-pump PLL circuit. *Technical report, Texas*

[19] Tse, C.K. (1994). Chaos from a Buck switching regulator operating in discontinuous mode. *IEEE Transactions on International Journal of Circuit Theory and Application*. Vol. 22,

[20] Tse, C.K. (1994). Flip bifurcation and chaos in three–state boost switching regulators. *IEEE Transactions on Circuits and Systems I: Theory and Applications*, Vol. 41, No. 1, pp.

[22] Van Paemel, M. (July 1994). Analysis of a charge pump PLL: a new model. *IEEE*

[23] Yuan, G. H.; Banerjee, S.; Ott, E. & Yorke, J. A. (1998). Border-collision bifurcations in the buck converter, *IEEE Trans. on Circuits and Systems-I*, Vol. 45, pp. 707–715.

[21] Tse, C.K. (2003). *Complex behavior of switching Power converters*, CRC Press.

*Transactions on Communications*, Vol. 42, No. 7, pp. 2490-2498.

*International Journal of Bifurcation and Chaos*, Vol. 5, pp. 1549–1578.

Institut National des Sciences Appliquées de Toulouse, France.

Nottingham, United Kingdom, December, pp. 19–25.

Institut National Polytechnique de Grenoble, France.

[18] Mira, C. (1987). Chaotic dynamics. *World scientific Publishing*.

*Instrument Technical Journal*. Vol. 14, No. 2.

Most of the dynamical systems analysed in engineering practice have the dynamics, which depends on both position and time. Such systems are classified as distributed parameter systems (DPS). The time-space coupled nature of the DPS is usually mathematically described by partial differential equations (PDE) as infinite-dimensional systems. However, from point of view of implementation of DPS control in technological practice, where a finite number of sensors and actuators for practical sensing and control is at disposal, such infinite-dimensional systems need to be approximated by finite-dimensional systems. There are many dimension reduction methods, which can be used to solve this problem.

In the first mathematical foundations of DPS control, analytical solutions of the underlying PDE have been used (Butkovskij, 1965; Lions, 1971; Wang, 1964). That is the decomposition of dynamics into time and space components based on the eigenfunctions of the PDE. Continuous and approximation theories aimed to control of parabolic systems presents monograph (Lasiecka & Triggiani, 2000). Methodical approach from the view of time-space separation with model reduction is presented in (Li & Qi, 2010). Variety of transfer functions for systems described by PDE are illustrated by means of several examples in (Curtain & Morris, 2009). Well-known reduction methods based on finite difference method (FDM), or finite element method (FEM), spectral method require an accurate nominal PDE model and usually lead to a high-order model, which requires unpractical high-order controller.

An engineering approach for the control of DPS is being developed since the eighties of the last century (Hulkó et al., 1981, 1987, 1998, 2009a, 2009b). In the field of lumped parameters system (LPS) control, where the state/output quantities *x(t)/y(t)* – parameters are given as

© 2012 Belavý et al., licensee InTech. This is an open access chapter distributed under the terms of the Creative Commons Attribution License (http://creativecommons.org/licenses/by/3.0), which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited. © 2012 Belavý et al., licensee InTech. This is a paper distributed under the terms of the Creative Commons Attribution License (http://creativecommons.org/licenses/by/3.0), which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.

finite dimensional vectors, the actuator together with the controlled plant make up a controlled LPS. In this sense the actuators and the controlled plant as a DPS create a controlled lumped-input and distributed-parameter-output system (LDS).

Robust Control of Distributed Parameter Systems with

Demonstration in Casting Technology and MATLAB/Simulink/DPS Blockset Software Support 31

Control Toolbox, Optimization Toolbox, System Identification Toolbox, Real-Time Windows

In general, DPS are systems whose state or output quantities, *X xyzt Y xyzt* , ,, / ,,, are distributed quantities or fields of quantities, where *xyz* , , are spatial coordinates in 3D. These systems are often considered as systems whose dynamics is described by PDE. In the inputoutput relation, PDE define distributed-input/distributed-output systems (DDS) between distributed input, *U xyzt* , , , and distributed output quantities, *Y xyzt* , , , , at initial and boundary conditions given. Distributed parameter systems are very frequently found in various technical and non-technical branches with limited number of manipulated input quantities, or actuators. These lumped input quantities by means of interaction of fields and quantities generate distributed output of real DPS. Representation of such DPS is either in the form of LDS, Fig. 1 a), or in the form of LDS with zero-order hold unit H (HLDS), when discrete-time lumped input quantities are used, Fig. 1. b) , (Hulkó et al., 1981, 1987, 1998).

**Figure 1.** Representation of DPS: a) LDS - lumped-input/distributed-output system, b) HLDS –LDS with block of zero-order hold units H, *U xyzt* ,,, - distributed input quantity, *Y xyzt* ,,, - distributed output quantity, 1, *<sup>i</sup> i n U t* - lumped input quantities, *U k* - vector of discrete-time lumped input

Distributed output of the linear LDS from zero initial conditions either in continuous, or in

 1 1 ,,, *n n*

*Yxk Y xk H xk U k* 

*i i Y xt Y xt xt U t* 

1 1 ,, , *n n*

*i i*

*i ii*

*i ii*

*G* (1)

*<sup>G</sup>* (2)

Target and Simulink Design Optimization.

quantities

**2.1. Dynamics of LDS** 

discrete-time (DT) is in the form:

**2. LDS/HLDS representation of DPS** 

In this chapter the decomposition of dynamics of controlled LDS into time and space components is introduced. Based on this decomposition a methodical framework of control synthesis decomposition into space and time tasks will be presented. In the space domain, approximation problems are solved. In the time domain, synthesis of control is performed by lumped parameter control loops, where robust controllers are used.

The casting technology is a typical case of the DPS. There in order to obtain the desired solidification structure, the casting process requires a specific temperature field of the mould, which is defined on complex-shape 3D definition domain. Modelling, simulation and evaluation of real-time experiments in this area is now widely accepted as an important tool in product design and process development to improve productivity and casting quality. For analysis of the casting process dynamics as DPS, especially temperature fields in the casting mould and control synthesis purposes, the benchmark casting plant with steel mould of complex-shape was designed at Faculty of Mechanical Engineering STU in Bratislava.

The main emphasis of this chapter is to present an engineering approach for the robust control of DPS with demonstration in the casting technology along with software support in the MATLAB & Simulink programming environment. This approach opens a wide space for novel applications of the toolboxes and blocksets of the MATLAB & Simulink software environment. For the software support of modelling, control and design of DPS, given on 1D-3D definition domains, the **Distributed Parameter Systems Blockset for MATLAB & Simulink (DPS Blockset),** a Third-Party Product of The MathWorks Company www.mathworks.com/products/connections/ has been developed at the Institute of Automation, Measurement and Applied Informatics, Faculty of Mechanical Engineering STU, (Hulkó et al. 2003-2010). Also a web portal named **Distributed Parameter Systems Control** - www.dpscontrol.sk has been created for those, who are interested in solving problems of DPS control (Hulkó et al., 2003-2007). This web portal contains application examples from different areas of engineering practice, such as the control of technological and manufacturing processes, mechatronic structures, groundwater remediation, etc. In addition, this web portal offers the demo version of the **DPS Blockset** with the **Tutorial, Show, Demos** and **DPS Wizard** for download, along with the **Interactive Control** service for the interactive solution of model control problems via the Internet.

In this chapter, for the control synthesis purpose, LDS models of temperature fields in the casting mould were created by means of evaluation of real-time experiments. Robust control synthesis based on internal model control (IMC) structure in the time domain has been done. Designed robust controllers were used for the robust control of preheating casting mould in the real-time experiment in accordance to casting technology requirements. Identification, uncertainty analysis of the models, robust control synthesis and experiments were performed with the software support of DPS Blockset and toolboxes of the MATLAB & Simulink, especially the System Identification Toolbox, Control Systems Toolbox, Robust Control Toolbox, Optimization Toolbox, System Identification Toolbox, Real-Time Windows Target and Simulink Design Optimization.

## **2. LDS/HLDS representation of DPS**

30 MATLAB – A Fundamental Tool for Scientific Computing and Engineering Applications – Volume 3

controlled lumped-input and distributed-parameter-output system (LDS).

by lumped parameter control loops, where robust controllers are used.

for the interactive solution of model control problems via the Internet.

Bratislava.

finite dimensional vectors, the actuator together with the controlled plant make up a controlled LPS. In this sense the actuators and the controlled plant as a DPS create a

In this chapter the decomposition of dynamics of controlled LDS into time and space components is introduced. Based on this decomposition a methodical framework of control synthesis decomposition into space and time tasks will be presented. In the space domain, approximation problems are solved. In the time domain, synthesis of control is performed

The casting technology is a typical case of the DPS. There in order to obtain the desired solidification structure, the casting process requires a specific temperature field of the mould, which is defined on complex-shape 3D definition domain. Modelling, simulation and evaluation of real-time experiments in this area is now widely accepted as an important tool in product design and process development to improve productivity and casting quality. For analysis of the casting process dynamics as DPS, especially temperature fields in the casting mould and control synthesis purposes, the benchmark casting plant with steel mould of complex-shape was designed at Faculty of Mechanical Engineering STU in

The main emphasis of this chapter is to present an engineering approach for the robust control of DPS with demonstration in the casting technology along with software support in the MATLAB & Simulink programming environment. This approach opens a wide space for novel applications of the toolboxes and blocksets of the MATLAB & Simulink software environment. For the software support of modelling, control and design of DPS, given on 1D-3D definition domains, the **Distributed Parameter Systems Blockset for MATLAB & Simulink (DPS Blockset),** a Third-Party Product of The MathWorks Company www.mathworks.com/products/connections/ has been developed at the Institute of Automation, Measurement and Applied Informatics, Faculty of Mechanical Engineering STU, (Hulkó et al. 2003-2010). Also a web portal named **Distributed Parameter Systems Control** - www.dpscontrol.sk has been created for those, who are interested in solving problems of DPS control (Hulkó et al., 2003-2007). This web portal contains application examples from different areas of engineering practice, such as the control of technological and manufacturing processes, mechatronic structures, groundwater remediation, etc. In addition, this web portal offers the demo version of the **DPS Blockset** with the **Tutorial, Show, Demos** and **DPS Wizard** for download, along with the **Interactive Control** service

In this chapter, for the control synthesis purpose, LDS models of temperature fields in the casting mould were created by means of evaluation of real-time experiments. Robust control synthesis based on internal model control (IMC) structure in the time domain has been done. Designed robust controllers were used for the robust control of preheating casting mould in the real-time experiment in accordance to casting technology requirements. Identification, uncertainty analysis of the models, robust control synthesis and experiments were performed with the software support of DPS Blockset and toolboxes of the MATLAB & Simulink, especially the System Identification Toolbox, Control Systems Toolbox, Robust In general, DPS are systems whose state or output quantities, *X xyzt Y xyzt* , ,, / ,,, are distributed quantities or fields of quantities, where *xyz* , , are spatial coordinates in 3D. These systems are often considered as systems whose dynamics is described by PDE. In the inputoutput relation, PDE define distributed-input/distributed-output systems (DDS) between distributed input, *U xyzt* , , , and distributed output quantities, *Y xyzt* , , , , at initial and boundary conditions given. Distributed parameter systems are very frequently found in various technical and non-technical branches with limited number of manipulated input quantities, or actuators. These lumped input quantities by means of interaction of fields and quantities generate distributed output of real DPS. Representation of such DPS is either in the form of LDS, Fig. 1 a), or in the form of LDS with zero-order hold unit H (HLDS), when discrete-time lumped input quantities are used, Fig. 1. b) , (Hulkó et al., 1981, 1987, 1998).

**Figure 1.** Representation of DPS: a) LDS - lumped-input/distributed-output system, b) HLDS –LDS with block of zero-order hold units H, *U xyzt* ,,, - distributed input quantity, *Y xyzt* ,,, - distributed output quantity, 1, *<sup>i</sup> i n U t* - lumped input quantities, *U k* - vector of discrete-time lumped input quantities

#### **2.1. Dynamics of LDS**

Distributed output of the linear LDS from zero initial conditions either in continuous, or in discrete-time (DT) is in the form:

$$Y\left(\overline{\boldsymbol{x}},t\right) = \sum\_{i=1}^{n} Y\_i\left(\overline{\boldsymbol{x}},t\right) = \sum\_{i=1}^{n} G\_i\left(\overline{\boldsymbol{x}},t\right) \otimes \mathcal{U}\_i\left(t\right) \tag{1}$$

$$Y\left(\overline{\boldsymbol{x}},k\right) = \sum\_{i=1}^{n} Y\_i\left(\overline{\boldsymbol{x}},k\right) = \sum\_{i=1}^{n} \mathcal{G}H\_i\left(\overline{\boldsymbol{x}},k\right) \oplus \mathcal{U}\_i\left(k\right) \tag{2}$$

where denotes convolution product and denotes convolution sum, *x xyz* , , is position vector in 3D, , *<sup>i</sup> G x t* - distributed parameter impulse response of LDS to the i-th input, , *H xk <sup>i</sup> G* - DT distributed parameter impulse response of LDS with zero-order hold units H (HLDS) to the i-th input, , *Y xt <sup>i</sup>* - distributed parameter output quantity of LDS to the i-th input, , *Y xk <sup>i</sup>* - DT distributed parameter output quantity of HLDS to the i-th input, *U t <sup>i</sup>* - lumped input quantity, *U k <sup>i</sup>* - DT lumped input quantity, (Hulkó et al. 1998).

When *U t <sup>i</sup>* is a unit-step (Heaviside) function, , *Y xt <sup>i</sup>* is in the form of distributed step response function , *<sup>i</sup> H x t* . Similarly, for the unit-step function *U k <sup>i</sup>* : , , *Y xk H xk i i H* . For simplicity in this chapter distributed quantities are considered mostly as continuous scalar quantity fields with unit sampling interval in the time domain. Whereas DT distributed parameter step responses *<sup>i</sup>* , *<sup>i</sup> HH xk* of HLDS can be computed by common analytical or numerical methods, then DT distributed parameter impulse responses can be obtained as

$$\left\{\mathcal{G}\mathcal{H}\_i\left(\overline{\mathbf{x}},k\right) = \mathcal{H}\mathcal{H}\_i\left(\overline{\mathbf{x}},k\right) - \mathcal{H}\mathcal{H}\_i\left(\overline{\mathbf{x}},k-1\right)\right\}\_i\tag{3}$$

Robust Control of Distributed Parameter Systems with

Demonstration in Casting Technology and MATLAB/Simulink/DPS Blockset Software Support 33

**Figure 2.** Distributed parameter feedback control loop: HLDS - LDS with zero-order holds *<sup>i</sup><sup>i</sup> <sup>H</sup>* on the input, CS - control synthesis, TS - control synthesis in time domain, SS - control synthesis in space domain, K - time/space sampling, *V xt* , - disturbance quantity, *Y xt* , - distributed controlled

*W k*

Let us consider a step change of distributed parameter control quantity *W xk W x* , , and *V xt* , 0 . The goal of the control synthesis is to generate a sequence of control inputs *U k* in such manner, that in the steady-state, for *k* , the control error *Exk* , will

First, in the SS blocks, the approximation both of sampled distributed controlled quantity *Yxk* , and reference quantity *W x* , , on the set of reduced steady-state distributed step

> *n n ii i i i <sup>Y</sup> i i Y x k Y x k HR x Y x k Y k HR x*

 1 1 min , ,, , ,

 1 1 min , , , , ,

*n n*

functions in the strictly convex normed linear space of distributed parameter quantities with

*ii i i i <sup>W</sup> i i W x W x HR x W x W HR x* 

*H H* (8)

*H H* (9)

*HHR x* form a finite-dimensional subspace of approximation

*Y k*

min , min , , , *Ex Wx Yx Ex* (7)

*C z* - lumped parameter controllers, *<sup>i</sup> <sup>i</sup>*


*U k* - lumped control


quantity, *Yxk* , - sampled distributed controlled quantity, *<sup>i</sup> <sup>i</sup>*

approach its minimal value *E x* , in the quadratic norm:

*HHR x* , are solved in following form:

controlled quantity, *W xk* , - reference quantity, *<sup>i</sup> <sup>i</sup>*

*E k* - control errors, *<sup>i</sup> <sup>i</sup>*

quantity, *<sup>i</sup> <sup>i</sup>*

responses *<sup>i</sup>* , *<sup>i</sup>*

*i*

*i*

Basis functions *<sup>i</sup>* , *<sup>i</sup>*

quantities

For points *i iii* , , *<sup>i</sup> x xyz* located in surroundings of lumped input quantities *<sup>i</sup> <sup>i</sup> U t* , where partial distributed transient responses *i i* , *<sup>i</sup> H x t* attains maximal amplitudes, partial distributed output quantities are obtained in time-domain and next either continuous *i i* , *<sup>i</sup> Sxs* , or discrete transfer functions *i i* , *<sup>i</sup> SH x z* with sampling period *T* are identified.

$$\left\{ Y\_i \left( \overline{\mathbf{x}}\_i, t \right) = \mathcal{G}\_i \left( \overline{\mathbf{x}}\_i, t \right) \otimes \mathcal{U}\_i \left( t \right) \right\}\_i \quad \rightarrow \quad \left\{ Y\_i \left( \overline{\mathbf{x}}\_i, s \right) = \mathcal{G}\_i \left( \overline{\mathbf{x}}\_i, s \right) \mathcal{U}\_i \left( s \right) \right\}\_i \tag{4}$$

$$\left\{\mathbf{Y}\_{i}\left(\overline{\mathbf{x}}\_{i},k\right) = \operatorname{\mathcal{G}H}\_{i}\left(\overline{\mathbf{x}}\_{i},k\right) \oplus \operatorname{\mathcal{U}}\_{i}\left(k\right)\right\}\_{i} \quad \rightarrow \quad \left\{\mathbf{Y}\_{i}\left(\overline{\mathbf{x}}\_{i},z\right) = \operatorname{\mathcal{G}H}\_{i}\left(\overline{\mathbf{x}}\_{i},z\right)\operatorname{\mathcal{U}}\_{i}\left(z\right)\right\}\_{i} \tag{5}$$

For the space dependency and in the steady-state we can define reduced transient step responses between i-th input quantity at point , , *i iii x xyz* and corresponding partial distributed output quantity in the steady-state:

$$\left\{ \mathcal{H} \text{HR}\_i \left( \overline{\underline{x}}, \infty \right) = \frac{\mathcal{H} \text{H}\_i \left( \overline{\underline{x}}, \infty \right)}{\mathcal{H} \text{H}\_i \left( \overline{\underline{x}}\_i, \infty \right)} \right\}\_i \tag{6}$$

for , <sup>0</sup> *i i <sup>i</sup> HH x* .

Dynamics of LDS/HLDS is decomposed to the time and space components:

*Time Components of Dynamics i i* , *<sup>i</sup> Sxs* , or *i i* , *<sup>i</sup> SH x z* - for given *i* and chosen *<sup>i</sup> x*

*Space Components of Dynamics <sup>i</sup>* , *<sup>i</sup> HHR x* - for given i in ∞

#### **2.2. Feedback control loop based on HLDS dynamics**

Decomposition of dynamics enables also to decompose the control synthesis (CS) to time synthesis (TS) and space synthesis (SS) tasks in the feedback control loop of the distributed parameter system, Fig. 2.

parameter step responses *<sup>i</sup>* , *<sup>i</sup>*

*i i* , *<sup>i</sup>*

for , <sup>0</sup> *i i <sup>i</sup> HH x* .

parameter system, Fig. 2.

where partial distributed transient responses *i i* , *<sup>i</sup>*

*Sxs* , or discrete transfer functions *i i* , *<sup>i</sup>*

distributed output quantity in the steady-state:

*Time Components of Dynamics i i* , *<sup>i</sup>*

*Space Components of Dynamics <sup>i</sup>* , *<sup>i</sup>*

**2.2. Feedback control loop based on HLDS dynamics** 

where denotes convolution product and denotes convolution sum, *x xyz* , , is position vector in 3D, , *<sup>i</sup> G x t* - distributed parameter impulse response of LDS to the i-th input, , *H xk <sup>i</sup> G* - DT distributed parameter impulse response of LDS with zero-order hold units H (HLDS) to the i-th input, , *Y xt <sup>i</sup>* - distributed parameter output quantity of LDS to the i-th input, , *Y xk <sup>i</sup>* - DT distributed parameter output quantity of HLDS to the i-th input, *U t <sup>i</sup>* - lumped input quantity, *U k <sup>i</sup>* - DT lumped input quantity, (Hulkó et al. 1998).

When *U t <sup>i</sup>* is a unit-step (Heaviside) function, , *Y xt <sup>i</sup>* is in the form of distributed step response function , *<sup>i</sup> H x t* . Similarly, for the unit-step function *U k <sup>i</sup>* : , , *Y xk H xk i i H* . For simplicity in this chapter distributed quantities are considered mostly as continuous scalar quantity fields with unit sampling interval in the time domain. Whereas DT distributed

numerical methods, then DT distributed parameter impulse responses can be obtained as

For points *i iii* , , *<sup>i</sup> x xyz* located in surroundings of lumped input quantities *<sup>i</sup> <sup>i</sup>*

partial distributed output quantities are obtained in time-domain and next either continuous

*ii ii i* , , *ii ii i* , , *i i*

*ii ii i* , , *ii ii i* , , *i i*

For the space dependency and in the steady-state we can define reduced transient step responses between i-th input quantity at point , , *i iii x xyz* and corresponding partial

*H x*

*H x*

, , , *i*

*H*

 *H*

*Sxs* , or *i i* , *<sup>i</sup>*

*HHR x* - for given i in ∞

Decomposition of dynamics enables also to decompose the control synthesis (CS) to time synthesis (TS) and space synthesis (SS) tasks in the feedback control loop of the distributed

*i*

*HR x*

Dynamics of LDS/HLDS is decomposed to the time and space components:

*H*

*HH xk* of HLDS can be computed by common analytical or

*iii* , , ,1 *<sup>i</sup> GHH H xk H xk H xk* (3)

*Y x t x t U t Y x s S x sU s G* (4)

*Y x k H x k U k Y x z SH x z U z G* (5)

*i i i*

*SH x z* - for given *i* and chosen *<sup>i</sup> x*

*H x t* attains maximal amplitudes,

(6)

*SH x z* with sampling period *T* are identified.

*U t* ,

**Figure 2.** Distributed parameter feedback control loop: HLDS - LDS with zero-order holds *<sup>i</sup><sup>i</sup> <sup>H</sup>* on the input, CS - control synthesis, TS - control synthesis in time domain, SS - control synthesis in space domain, K - time/space sampling, *V xt* , - disturbance quantity, *Y xt* , - distributed controlled quantity, *Yxk* , - sampled distributed controlled quantity, *<sup>i</sup> <sup>i</sup> Y k* - approximation parameters of controlled quantity, *W xk* , - reference quantity, *<sup>i</sup> <sup>i</sup> W k* - approximation parameters of reference quantity, *<sup>i</sup> <sup>i</sup> E k* - control errors, *<sup>i</sup> <sup>i</sup> C z* - lumped parameter controllers, *<sup>i</sup> <sup>i</sup> U k* - lumped control quantities

Let us consider a step change of distributed parameter control quantity *W xk W x* , , and *V xt* , 0 . The goal of the control synthesis is to generate a sequence of control inputs *U k* in such manner, that in the steady-state, for *k* , the control error *Exk* , will approach its minimal value *E x* , in the quadratic norm:

$$\min \left\| E\left(\overline{\underline{x}}, \infty\right) \right\| = \min \left\| \mathcal{W}\left(\overline{\underline{x}}, \infty\right) - \mathcal{Y}\left(\overline{\underline{x}}, \infty\right) \right\| = \left\| \bar{E}\left(\overline{\underline{x}}, \infty\right) \right\| \tag{7}$$

First, in the SS blocks, the approximation both of sampled distributed controlled quantity *Yxk* , and reference quantity *W x* , , on the set of reduced steady-state distributed step responses *<sup>i</sup>* , *<sup>i</sup> HHR x* , are solved in following form:

$$\min\_{\tilde{Y}\_i} \left\| Y(\overline{\boldsymbol{x}}, k) - \sum\_{i=1}^n Y\_i(\overline{\boldsymbol{x}}\_i, k) \mathcal{H} \mathcal{H} \mathcal{R}\_i(\overline{\boldsymbol{x}}, \boldsymbol{\infty}) \right\| = \left\| Y(\overline{\boldsymbol{x}}, k) - \sum\_{i=1}^n \bar{Y}\_i(k) \mathcal{H} \mathcal{H} \mathcal{R}\_i(\overline{\boldsymbol{x}}, \boldsymbol{\infty}) \right\| \tag{8}$$

$$\min\_{\mathcal{W}\_i} \left\| W(\overline{\boldsymbol{x}}, \boldsymbol{\infty}) - \sum\_{i=1}^n W\_i(\overline{\boldsymbol{x}}\_i, \boldsymbol{\infty}) \mathcal{H} \mathcal{H} \mathcal{R}\_i(\overline{\boldsymbol{x}}, \boldsymbol{\infty}) \right\| = \left\| W(\overline{\boldsymbol{x}}, \boldsymbol{\infty}) - \sum\_{i=1}^n \breve{W}\_i \, \mathcal{H} \mathcal{H} \mathcal{R}\_i(\overline{\boldsymbol{x}}, \boldsymbol{\infty}) \right\| \tag{9}$$

Basis functions *<sup>i</sup>* , *<sup>i</sup> HHR x* form a finite-dimensional subspace of approximation functions in the strictly convex normed linear space of distributed parameter quantities with quadratic norm, where the approximation problem is solved. From approximation theory involves, that solution of the approximation problems (8), (9) is guaranteed as a unique the best approximation in the form 1 , *n i i i Y k HR x <sup>H</sup>* with the vector of optimal approximation parameters *<sup>i</sup> <sup>i</sup> Y k* in task (8) and the best approximation in the form 1 , *n i i i W HR x <sup>H</sup>* with the vector of optimal approximation parameters *<sup>i</sup><sup>i</sup> W* for approximation task (9).

Robust Control of Distributed Parameter Systems with

*Sxs*

*Sxs* , takes

Demonstration in Casting Technology and MATLAB/Simulink/DPS Blockset Software Support 35

In general, a mathematical model for the plant dynamics is the basis for analysis and design of control systems. Also for LDS representation of DPS lumped and distributed models are used. However, in practice, no mathematical model describes exactly a physical process. It is obvious, that although no model represents the process exactly, some of them will do so

The theory of the robust control represents one of the possible approaches to the control system design in the presence of uncertainty. The goal of the robust system design is to retain a good quality of system performance in spite of model inaccuracies and changes. For the design techniques, the following requirements are supposed to be fulfilled: formulation of nominal plant model, different plant uncertainty models and requirements for both,

LDS representation of DPS means decomposition of dynamics to space and time

In distributed parameter control system, according to Fig. 2, single-input, single-output control loops in the block TS are tuned as closed feedback control loops using usual methods. In these loops, as models of the controlled system, transfer functions *i i* , *<sup>i</sup>*

*U k* and *i i* , *<sup>i</sup>*

invariant models, *<sup>i</sup>* . This family in the frequency domain, e.g. for models *i i* , *<sup>i</sup>*

 procedure of dynamics modelling and possible change of parameters in models (4), (5) solution of approximation problem (8), (9), where lumped quantities are obtained

In order to treat uncertainties, it will be further assumed that the dynamic behaviour of a plant is described not by a single linear time invariant model, but by a family of linear time

where , *i i Sxj* is the nominal plant model. Any member of the family *i* fulfils the

where *ai L j* is an additive uncertainty and *ai L* is the bound of additive uncertainty. If

we wish to work with multiplicative uncertainties, we define the relations:

*SH x z* in the z-domain are used. These transfer functions describe the

*i i i i i i ai S Sxj Sxj L* :, , (11)

, , *i i i i ai Sxj Sxj L j* (12)

, *ai ai i i Lj L S* (13)

*Yxk* .

components. Uncertainties may occur in both, time and space components.

**3. Robust control system** 

with greater accuracy than others.

robust stability and performance.

dynamics between sequences *<sup>i</sup> <sup>i</sup>*

and/or *i i* , *<sup>i</sup>*

the following form:

conditions:

**3.1. Sources of uncertainties in the LDS structure** 

In this case, the sources of uncertainties are given by:

Let us formulate a DPS control problem for the distributed reference quantity *W x* , . When *W xk* , is assumed, the space control synthesis is performed in each time step *k*, which gives *<sup>i</sup> <sup>i</sup> W k* parameters. Graphical interpretation of the approximation problem (9) for HLDS defined on 1D space *x* is on Fig. 3.

**Figure 3.** Solution of the approximation problem for the distributed reference quantity *W x* ,

Next, based on the solution of approximation problem, the vector of control error is created:

$$\left\{\overline{E}\left(k\right) = \left\{E\_i\left(k\right)\right\}\_i = \left\{\bar{W}\_i - \bar{Y}\_i\left(k\right)\right\}\_i \tag{10}$$

The control errors vector *<sup>i</sup> <sup>i</sup> Ek E k* enters into the block TS, where the vector of control quantities, *<sup>i</sup> <sup>i</sup> Uk U k* is generated by controllers *<sup>i</sup> <sup>i</sup> C z* in single-parameter control loops. During the control process, for *k* the control task (7) is accomplished.

Finally, we may state as a summary, that in the feedback control of DPS with dynamics represented in the form of HLDS, the control synthesis is performed as:

*Space Tasks of Control Synthesis* – as approximation tasks.

*Time Tasks of Control Synthesis* – on the level of lumped parameter control loops.

## **3. Robust control system**

34 MATLAB – A Fundamental Tool for Scientific Computing and Engineering Applications – Volume 3

best approximation in the form

*Y k*

approximation parameters *<sup>i</sup> <sup>i</sup>*

*W k*

(9) for HLDS defined on 1D space *x* is on Fig. 3.

The control errors vector *<sup>i</sup> <sup>i</sup>*

control quantities, *<sup>i</sup> <sup>i</sup>*

,

approximation task (9).

which gives *<sup>i</sup> <sup>i</sup>*

1

*i*

*i i*

*W HR x*

*n*

quadratic norm, where the approximation problem is solved. From approximation theory involves, that solution of the approximation problems (8), (9) is guaranteed as a unique the

*i i*

*Y k HR x*

,

parameters. Graphical interpretation of the approximation problem

*<sup>H</sup>* with the vector of optimal

*W* for

in task (8) and the best approximation in the form

1

**Figure 3.** Solution of the approximation problem for the distributed reference quantity *W x* ,

Next, based on the solution of approximation problem, the vector of control error is created:

*i ii <sup>i</sup> <sup>i</sup>*

*Uk U k* is generated by controllers *<sup>i</sup> <sup>i</sup>*

Finally, we may state as a summary, that in the feedback control of DPS with dynamics

control loops. During the control process, for *k* the control task (7) is accomplished.

represented in the form of HLDS, the control synthesis is performed as:

*Time Tasks of Control Synthesis* – on the level of lumped parameter control loops.

*Space Tasks of Control Synthesis* – as approximation tasks.

*Ek E k W Y k* (10)

*C z* in single-parameter

*Ek E k* enters into the block TS, where the vector of

*<sup>H</sup>* with the vector of optimal approximation parameters *<sup>i</sup><sup>i</sup>*

Let us formulate a DPS control problem for the distributed reference quantity *W x* , . When *W xk* , is assumed, the space control synthesis is performed in each time step *k*,

*i*

*n*

In general, a mathematical model for the plant dynamics is the basis for analysis and design of control systems. Also for LDS representation of DPS lumped and distributed models are used. However, in practice, no mathematical model describes exactly a physical process. It is obvious, that although no model represents the process exactly, some of them will do so with greater accuracy than others.

The theory of the robust control represents one of the possible approaches to the control system design in the presence of uncertainty. The goal of the robust system design is to retain a good quality of system performance in spite of model inaccuracies and changes. For the design techniques, the following requirements are supposed to be fulfilled: formulation of nominal plant model, different plant uncertainty models and requirements for both, robust stability and performance.

## **3.1. Sources of uncertainties in the LDS structure**

LDS representation of DPS means decomposition of dynamics to space and time components. Uncertainties may occur in both, time and space components.

In distributed parameter control system, according to Fig. 2, single-input, single-output control loops in the block TS are tuned as closed feedback control loops using usual methods. In these loops, as models of the controlled system, transfer functions *i i* , *<sup>i</sup> Sxs* and/or *i i* , *<sup>i</sup> SH x z* in the z-domain are used. These transfer functions describe the dynamics between sequences *<sup>i</sup> <sup>i</sup> U k* and *i i* , *<sup>i</sup> Yxk* .

In this case, the sources of uncertainties are given by:


In order to treat uncertainties, it will be further assumed that the dynamic behaviour of a plant is described not by a single linear time invariant model, but by a family of linear time invariant models, *<sup>i</sup>* . This family in the frequency domain, e.g. for models *i i* , *<sup>i</sup> Sxs* , takes the following form:

$$\Psi\_i = \left\{ S\_i : \left| \begin{array}{c} S\_i \left( \overline{x}\_{i'} \text{joo} \right) \ - & \tilde{S}\_i \left( \overline{x}\_{i'} \text{joo} \right) \end{array} \right| \le \overline{L}\_{ai} \left( \text{oo} \right) \right\} \tag{11}$$

where , *i i Sxj* is the nominal plant model. Any member of the family *i* fulfils the conditions:

$$S\_i\left(\overline{\underline{x}}\_{i'}\right)\text{oo}\right) = \begin{array}{c} \tilde{S}\_i\left(\overline{\underline{x}}\_{i'}\right)\text{oo}\end{array} + \begin{array}{c} L\_{ai}\left(\text{jo}\right) \tag{12}$$

$$\left| L\_{ai} \left( j \alpha \right) \right| \le \left. \overline{L}\_{ai} \left( \alpha \right) \right| \quad \text{ } \forall \text{ } S\_i \in \Psi\_i \tag{13}$$

where *ai L j* is an additive uncertainty and *ai L* is the bound of additive uncertainty. If we wish to work with multiplicative uncertainties, we define the relations:

$$L\_{mi}\left(j\alpha\right) = \frac{L\_{ai}\left(j\alpha\right)}{\tilde{S}\_i\left(\overline{x}\_{i'}j\alpha\right)} \; ; \qquad \overline{L}\_{mi}\left(\alpha\right) = \frac{\overline{L}\_{ai}\left(\alpha\right)}{\left|\overline{S}\_i\left(\overline{x}\_{i'}j\alpha\right)\right|}\tag{14}$$

Robust Control of Distributed Parameter Systems with

*Y k*

*U k* - lumped control quantities, *i i* , *<sup>i</sup>*

*W k*

*e z SH x z Q z z* (16)

, , , *i i Ni i Mi i SH x z SH x z SH x z* (17)

(18)

*W k*



in the form unit-step function

*SH x z* -

Demonstration in Casting Technology and MATLAB/Simulink/DPS Blockset Software Support 37

Structure of the distributed parameter feedback robust control system DPS with HLDS

**Figure 5.** Distributed parameter feedback robust control system: HLDS - LDS with zero-order holds

controlled quantity, *Yxk* , - sampled distributed controlled quantity, *<sup>i</sup> <sup>i</sup>*

*E k* - control errors, *<sup>i</sup> <sup>i</sup>*

models of lumped controlled systems, *<sup>i</sup> <sup>i</sup> Q z* - lumped parameter IMC controllers

are obtained from solution of the following minimization problem:

min min 1 ,

parameters of controlled quantity, *W xk* , - reference quantity, *<sup>i</sup> <sup>i</sup>*

*<sup>H</sup>*<sup>2</sup> optimal IMC controllers *<sup>i</sup> <sup>i</sup> Q z* for inputs *<sup>i</sup> <sup>i</sup>*

*i i*

First, factorize the nominal stable transfer function , *i i SH x z* :

subject to the constraint *Q z <sup>i</sup>* to be stable and causal.

After this, optimal IMC controller *Q z <sup>i</sup>* is given by:

*H* on the input, CS - control synthesis, TS - control synthesis in time domain, SS - control synthesis in space domain, K - time/space sampling, *V xt* , - disturbance quantity, *Y xt* , - distributed

<sup>2</sup> <sup>2</sup>

where *SH x z Ni i* , includes positive zeros or time-delays of the transfer function , *i i SH x z* .

<sup>1</sup> , *Q z SH x z i M i i*

Finally, controller *Q z <sup>i</sup>* is augmented by low-pass filter *<sup>i</sup> F z* with parameter 0 1 *<sup>i</sup>* :

*i i ii i Q z Q z*

dynamics and IMC controllers is on Fig. 5.

*<sup>i</sup><sup>i</sup>*

 <sup>1</sup> *z*

*z*

*z*

of reference quantity, *<sup>i</sup> <sup>i</sup>*

where *L j mi* is a multiplicative uncertainty and *Lmi* is the bound of multiplicative uncertainty.

#### **3.2. Design of IMC robust controllers**

A robust control system for HLDS will be designed using the Internal Model Control (IMC) strategy (Morari & Zafiriou, 1989) with the general structure depicted in Fig. 4. a). It is possible to transform this structure to the classical feedback control loop, Fig. 4. b) and to incorporate it into the TS block of the DPS feedback control system. The relationship between the classical feedback controller *C* and the IMC controller *Q* for the nominal model *P* of the controlled process *P* is as follows, and vice-versa:

$$Q = \frac{C}{1 + \tilde{P}C}; \qquad C = \frac{Q}{1 - \tilde{P}Q} \tag{15}$$

**Figure 4.** a) Internal Model Control structure, b) equivalent classical feedback control loop

It is well known that IMC strategy has the following properties:


The IMC structure thus provides the following benefits with respect to classical feedback: better dynamic response, system stability and robustness. One can search for *Q* instead of *C* without any loss of generality.

Structure of the distributed parameter feedback robust control system DPS with HLDS dynamics and IMC controllers is on Fig. 5.

36 MATLAB – A Fundamental Tool for Scientific Computing and Engineering Applications – Volume 3

*mi mi*

*L j L*

 ; , , *ai ai*

(14)

*C Q Q C PC PQ* (15)

*i i i i*

*Sxj Sxj*

*L j L*

where *L j mi* is a multiplicative uncertainty and *Lmi* is the bound of multiplicative

A robust control system for HLDS will be designed using the Internal Model Control (IMC) strategy (Morari & Zafiriou, 1989) with the general structure depicted in Fig. 4. a). It is possible to transform this structure to the classical feedback control loop, Fig. 4. b) and to incorporate it into the TS block of the DPS feedback control system. The relationship between the classical feedback controller *C* and the IMC controller *Q* for the nominal model

; 1 1

**Figure 4.** a) Internal Model Control structure, b) equivalent classical feedback control loop

stable, then the IMC structure guarantees the closed-loop stability.

1. *Dual stability*: Assume a perfect model, *P P* and if the controller and process are

2. *Perfect control*: Assume a perfect model, *P P* and the closed-loop system is stable, while <sup>1</sup> *Q P* , then there is no output steady-state error for set-point variance and

The IMC structure thus provides the following benefits with respect to classical feedback: better dynamic response, system stability and robustness. One can search for *Q* instead of

It is well known that IMC strategy has the following properties:

disturbances.

*C* without any loss of generality.

uncertainty.

**3.2. Design of IMC robust controllers** 

*P* of the controlled process *P* is as follows, and vice-versa:

**Figure 5.** Distributed parameter feedback robust control system: HLDS - LDS with zero-order holds *<sup>i</sup><sup>i</sup> H* on the input, CS - control synthesis, TS - control synthesis in time domain, SS - control synthesis in space domain, K - time/space sampling, *V xt* , - disturbance quantity, *Y xt* , - distributed controlled quantity, *Yxk* , - sampled distributed controlled quantity, *<sup>i</sup> <sup>i</sup> Y k* - approximation parameters of controlled quantity, *W xk* , - reference quantity, *<sup>i</sup> <sup>i</sup> W k* - approximation parameters of reference quantity, *<sup>i</sup> <sup>i</sup> E k* - control errors, *<sup>i</sup> <sup>i</sup> U k* - lumped control quantities, *i i* , *<sup>i</sup> SH x z* models of lumped controlled systems, *<sup>i</sup> <sup>i</sup> Q z* - lumped parameter IMC controllers

*<sup>H</sup>*<sup>2</sup> optimal IMC controllers *<sup>i</sup> <sup>i</sup> Q z* for inputs *<sup>i</sup> <sup>i</sup> W k* in the form unit-step function <sup>1</sup> *z z z* are obtained from solution of the following minimization problem:

$$\min\_{Q\_i(z)} \left\| e\_i(z) \right\|\_2 = \min\_{Q\_i(z)} \left\| \left( 1 - SH\_i(\overline{x}\_i, z) Q\_i(z) \right) \chi\_i(z) \right\|\_2 \tag{16}$$

subject to the constraint *Q z <sup>i</sup>* to be stable and causal.

First, factorize the nominal stable transfer function , *i i SH x z* :

$$SH\_i\left(\overline{\mathbf{x}}\_{i'}\mathbf{z}\right) = SH\_{Ni}\left(\overline{\mathbf{x}}\_{i'}\mathbf{z}\right) SH\_{Mi}\left(\overline{\mathbf{x}}\_{i'}\mathbf{z}\right) \tag{17}$$

where *SH x z Ni i* , includes positive zeros or time-delays of the transfer function , *i i SH x z* . After this, optimal IMC controller *Q z <sup>i</sup>* is given by:

$$\mathcal{Q}\_i(z) = \mathcal{S} H\_{Mi} \left( \overline{\mathbf{x}}\_{i'} z \right)^{-1} \tag{18}$$

Finally, controller *Q z <sup>i</sup>* is augmented by low-pass filter *<sup>i</sup> F z* with parameter 0 1 *<sup>i</sup>* :

$$F\_i(z) = \frac{1 - \alpha\_i}{z - \alpha\_i} \tag{19}$$

Robust Control of Distributed Parameter Systems with

Demonstration in Casting Technology and MATLAB/Simulink/DPS Blockset Software Support 39

induction motor driven roller vane pumps, a bunch of hoses, a flow divider, a collector with built-in check valves, a plate heat exchanger and an expansion tank. The main cooling circuit is divided into five independent circuits thus enabling the control of heat extraction from the casting via the chills. The coolant flow is controlled by means of frequency converters, since volumetric pumps are used. The main heating circuit is also divided into

Inside of the casting mould are built-in 26 electric heating elements, each with maximal heating power 400 W. Heating elements are grouped to 5 zones and their heating power is actuated by the input voltage range of (0 – 10) V. In the body of the mould is also placed 7 water-cooled copper chills and 11 thermocouples, Fig. 8., Fig 9. Coordinates of measuring

five independent circuits.

**Figure 6.** Scheme of the benchmark casting plant

**Figure 7.** Two-part steel mould in detail

Resulting IMC controller with filter *Q z Fi* is in following form:

$$Q\_{Fi}(z) = Q\_i(z) F\_i(z) = Q\_i(z) \frac{(1 - \alpha\_i)}{z - \alpha\_i} \tag{20}$$

Parameter of the filter *<sup>i</sup>* is the only one tuning parameter to be selected by the user to achieve the appropriate compromise between performance and robustness and to keep the action of the manipulated variable within bounds. It must be chosen with respect to both, robust stability and robust performance condition:

$$\left| F\_i \left( e^{j\alpha T} \right) \right| < \left[ \left| S\_i \left( \overline{\mathbf{x}}\_{i'} \text{jo} \right) \mathbb{Q}\_i \left( e^{j\alpha T} \right) \right| \overline{L}\_{\text{mi}} \left( \alpha \right) \right]^{-1}, \quad 0 \le \alpha \le \frac{\pi}{T} \tag{21}$$

$$\sum\_{i} Q\_{i} \begin{pmatrix} j \alpha \end{pmatrix} \left| \overline{L}\_{ai} \begin{pmatrix} \alpha \end{pmatrix} + \left| 1 - S\_{i} \begin{pmatrix} \overline{x}\_{i'} \end{pmatrix} \right> \mathcal{Q}\_{i} \begin{pmatrix} j \alpha \end{pmatrix} \right| G\_{wi} \begin{pmatrix} \alpha \end{pmatrix} < 1, \quad 0 \le \alpha \le \frac{\pi}{T} \tag{22}$$

where *Gwi* is weighting function

For , , *i i Mi i SH x z SH x z* and low-pass filter (19), robust controller *C z <sup>i</sup>* in equivalent classical feedback control loop takes form:

$$\mathcal{C}\_{i}(\mathbf{z}) = \frac{\operatorname{SH}\_{Mi}\left(\overline{\mathbf{x}}\_{i\prime}, \mathbf{z}\right)^{-1} F\_{i}(\mathbf{z})}{1 - \operatorname{SH}\_{i}\left(\overline{\mathbf{x}}\_{i\prime}, \mathbf{z}\right) \operatorname{SH}\_{Mi}\left(\overline{\mathbf{x}}\_{i\prime}, \mathbf{z}\right)^{-1} F\_{i}(\mathbf{z})} = \frac{1}{\operatorname{SH}\_{Mi}\left(\overline{\mathbf{x}}\_{i\prime}, \mathbf{z}\right)} \cdot \frac{F\_{i}(\mathbf{z})}{1 - F\_{i}(\mathbf{z})}\tag{23}$$

#### **4. Benchmark casting plant**

The casting mould is one of the key components of a casting. It is well known, that the quality of the castings is affected strongly by the surface quality and the distribution of temperature in the mould, which has both time and space dependence. For study of the physical phenomena occurring during the casting solidification, from a DPS control point of view, control system development as well as mathematical model validation, a benchmark of the casting processes was designed (Belavý et al., 2009). At the study of casting processes and design of experimental plant, simulation studies in virtual software environments ProCAST and COMSOL Multiphysics were used.

#### **4.1. Construction of the benchmark casting plant**

Scheme of the benchmark casting plant is depicted in Fig. 6. The core item is the two-part steel mould of a complex-shape mounted in the frame of the ejector mechanism, Fig. 7. This mechanism is hinge-mounted to the main frame, to enable tilting of the mould for optimal filling and metal flow. Further, a hydraulic cooling circuit, which consists of an array of induction motor driven roller vane pumps, a bunch of hoses, a flow divider, a collector with built-in check valves, a plate heat exchanger and an expansion tank. The main cooling circuit is divided into five independent circuits thus enabling the control of heat extraction from the casting via the chills. The coolant flow is controlled by means of frequency converters, since volumetric pumps are used. The main heating circuit is also divided into five independent circuits.

**Figure 6.** Scheme of the benchmark casting plant

38 MATLAB – A Fundamental Tool for Scientific Computing and Engineering Applications – Volume 3

Resulting IMC controller with filter *Q z Fi* is in following form:

robust stability and robust performance condition:

where *Gwi* is weighting function

*i*

**4. Benchmark casting plant** 

*C z*

classical feedback control loop takes form:

ProCAST and COMSOL Multiphysics were used.

**4.1. Construction of the benchmark casting plant** 

*i*

*F z*

*Fi i i i*

, , 0 *j T j T i i <sup>i</sup> i mi Fe Sxj Qe L <sup>T</sup>* 

1 , 1, 0 *Qj L Sxj Qj G i ai i i i wi <sup>T</sup>*

For , , *i i Mi i SH x z SH x z* and low-pass filter (19), robust controller *C z <sup>i</sup>* in equivalent

1

*Mi i i i*

*SH x z F z F z*

*Mi i <sup>i</sup> i i Mi i i*

*SH x z F z SH x z SH x z F z*

1 , 1

, 1 1, ,

The casting mould is one of the key components of a casting. It is well known, that the quality of the castings is affected strongly by the surface quality and the distribution of temperature in the mould, which has both time and space dependence. For study of the physical phenomena occurring during the casting solidification, from a DPS control point of view, control system development as well as mathematical model validation, a benchmark of the casting processes was designed (Belavý et al., 2009). At the study of casting processes and design of experimental plant, simulation studies in virtual software environments

Scheme of the benchmark casting plant is depicted in Fig. 6. The core item is the two-part steel mould of a complex-shape mounted in the frame of the ejector mechanism, Fig. 7. This mechanism is hinge-mounted to the main frame, to enable tilting of the mould for optimal filling and metal flow. Further, a hydraulic cooling circuit, which consists of an array of

*Q z Q zF z Q z*

<sup>1</sup> *<sup>i</sup>*

*z* 

<sup>1</sup> *<sup>i</sup>*

Parameter of the filter *<sup>i</sup>* is the only one tuning parameter to be selected by the user to achieve the appropriate compromise between performance and robustness and to keep the action of the manipulated variable within bounds. It must be chosen with respect to both,

*i*

*i*

1

(21)

(22)

*z*

(19)

(20)

 

(23)

**Figure 7.** Two-part steel mould in detail

Inside of the casting mould are built-in 26 electric heating elements, each with maximal heating power 400 W. Heating elements are grouped to 5 zones and their heating power is actuated by the input voltage range of (0 – 10) V. In the body of the mould is also placed 7 water-cooled copper chills and 11 thermocouples, Fig. 8., Fig 9. Coordinates of measuring

points of thermocouples in *x-y* plane are given in Tab. 1 and *z* - coordinate is -0.05 m. Location of built-in elements has been carefully designed based on simulation studies in software environments ProCAST and COMSOL Multiphysics, in order to have the possibility of preheating the mould in 5 zones achieving desired temperature profile as well as directional solidification of the casting by means of active heat removal. The temperature field in the mould-casting system is possible to estimate through interpolation of data, measured by thermocouples.

Robust Control of Distributed Parameter Systems with

Demonstration in Casting Technology and MATLAB/Simulink/DPS Blockset Software Support 41

**4.2. Measurement and control scheme in MATLAB & Simulink** 

**Figure 10.** Measurement and control scheme in the MATLAB & Simulink

Toolbox and the communication is triggered by *Connect to Target* icon.

The scheme is composed of three main subsystems, namely: SENSING, HEATING and COOLING. Utilizing these subsystems and communication interface between process and computer, it is possible to measure dynamical characteristics of temperature field at zone heating. It is also possible to execute the experiment of controlled preheating of permanent mould before casting operation and controlled cooling during casting solidification. The above mentioned communication interface consists of data acquisition cards Advantech PCI-1710 [Ah], PCI-1710 [Ch] as analog input and Humusoft AD622 [Eh] as digital input.

Communication interface is performed by means of Simulink Real-Time Windows Target

The SENSING subsystem enables the temperature measurement by thermocouples No. 1 to 11, which are permanently located in the bottom part of the steel mould, Fig. 9. These are marked on the scheme as P1 to P11. The subsystem also enables temperature measurements during casting solidification by temporal thermocouples W1 to W9, located in the casting

was setup, Fig. 10.

The measurement and control task of temperature field in the permanent casting mould was performed in the MATLAB & Simulink environment, where a *mould\_exp\_robust.mdl* scheme

**Figure 8.** Location of heating elements (red) and copper chills (blue) in the mould

**Figure 9.** Bottom side of the steel casting mould


**Table 1.** Coordinates of thermocouples in *x-y* plane

## **4.2. Measurement and control scheme in MATLAB & Simulink**

40 MATLAB – A Fundamental Tool for Scientific Computing and Engineering Applications – Volume 3

**Figure 8.** Location of heating elements (red) and copper chills (blue) in the mould

Position 1 2 3 4 5 6 7 8 9 10 11 *x* (m) 0 0 0 0.145 0.090 0.035 0 0 -0.035 -0.145 -0.255 *y* (m) 0.1585 0.0985 0.050 0 0 0 -0.135 -0.075 0 0 0

**Figure 9.** Bottom side of the steel casting mould

**Table 1.** Coordinates of thermocouples in *x-y* plane

measured by thermocouples.

points of thermocouples in *x-y* plane are given in Tab. 1 and *z* - coordinate is -0.05 m. Location of built-in elements has been carefully designed based on simulation studies in software environments ProCAST and COMSOL Multiphysics, in order to have the possibility of preheating the mould in 5 zones achieving desired temperature profile as well as directional solidification of the casting by means of active heat removal. The temperature field in the mould-casting system is possible to estimate through interpolation of data,

The measurement and control task of temperature field in the permanent casting mould was performed in the MATLAB & Simulink environment, where a *mould\_exp\_robust.mdl* scheme was setup, Fig. 10.

**Figure 10.** Measurement and control scheme in the MATLAB & Simulink

The scheme is composed of three main subsystems, namely: SENSING, HEATING and COOLING. Utilizing these subsystems and communication interface between process and computer, it is possible to measure dynamical characteristics of temperature field at zone heating. It is also possible to execute the experiment of controlled preheating of permanent mould before casting operation and controlled cooling during casting solidification. The above mentioned communication interface consists of data acquisition cards Advantech PCI-1710 [Ah], PCI-1710 [Ch] as analog input and Humusoft AD622 [Eh] as digital input.

Communication interface is performed by means of Simulink Real-Time Windows Target Toolbox and the communication is triggered by *Connect to Target* icon.

The SENSING subsystem enables the temperature measurement by thermocouples No. 1 to 11, which are permanently located in the bottom part of the steel mould, Fig. 9. These are marked on the scheme as P1 to P11. The subsystem also enables temperature measurements during casting solidification by temporal thermocouples W1 to W9, located in the casting

domain. Thermocouple sensing junctions are located in the middle of the arm cross section, right above permanent thermocouples "P", except the node of the casting, where one temporal thermocouple points to the center of the node. Sensors S1 to S7 measure the temperature of cooling water in embedded chills. *Display* blocks of the sensors "P" and "W" in the scheme correspond to the real positions of sensors. Measured temperatures are saved to the data file and continuously displayed on the *Scope* block.

Robust Control of Distributed Parameter Systems with

Demonstration in Casting Technology and MATLAB/Simulink/DPS Blockset Software Support 43

**Figure 11.** DPS Blockset on the web portal of The MathWorks Corporation

**Figure 12.** The library of DPS Blockset for MATLAB & Simulink

The HEATING subsystem has two basic operational regimes, which are activated by *Switch2* block. In the manual control regime, it is possible to set heating performance in range (0÷10) Volt and measure the temperature transient characteristics with P1 to P11 thermocouples in given locations of the mould. The second regime activated by *Switch2* block enables to perform controlled preheating of the casting mould to desired temperature profile W defined in 11 points, where thermocouples P1 to P11 are located. The control task is performed by *DPS Robust IMC Control Synthesis* block.

## **5. Distributed Parameter Systems Blockset for MATLAB & Simulink**

For the MATLAB & Simulink based software support of modelling, control and design of Distributed Parameter Systems given on complex 3D domains of definition, the programming environment **Distributed Parameter Systems Blockset for MATLAB & Simulink** (**DPS Blockset**) as Third-Party MathWorks Product has been developed by the Institute of Automation, Measurement and Applied Informatics, Faculty of Mechanical Engineering, Slovak University of Technology in Bratislava, within the program CONNECTIONS of The MathWorks Corporation, Fig. 11. , (Hulkó et al., 2003-2010).

The library of DPS Blockset shows Fig. 12. Blocks **HLDS** and **RHLDS** serve for modelling of distributed parameter systems as lumped-input/distributed-output systems with zero-order hold units. The block **DPS Control Synthesis** provides feedback to distributed parameter controlled systems in control loops with blocks for discrete-time **PID**, **Algebraic**, **State-Space** and **Robust Control**. The block **DPS Input** generates distributed quantities which can be used as distributed control quantities or distributed disturbances, etc. **DPS Display** presents distributed quantities with many options including export to AVI files. The block **DPS Space Synthesis** performs space synthesis as an approximation problem.

The block **DPS Wizard** in step-by-step operation, by means of several model examples with default parameters on 1D-3D definition domains, gives an automatized guide for arrangement and setting distributed parameter control loops. The block **Demos** contains examples oriented to methodology of modelling and control synthesis. The block **Show**  contains motivation examples such as: *Control of temperature field of 3D metal body* (the controlled system was modelled in the virtual software environment COMSOL Multiphysics)**;** *Control of 3D beam of "smart" structure* (the controlled system was modelled in the virtual software environment ANSYS)*; Adaptive control of glass furnace* (the controlled system was modelled by Partial Differential Equations Toolbox of the MATLAB ), and *Groundwater remediation control* (the controlled system was modelled in the virtual software environment MODFLOW)**.** The block **Tutorial** presents methodological framework both for formulation and solution of control tasks for distributed parameter systems.

#### Robust Control of Distributed Parameter Systems with

Demonstration in Casting Technology and MATLAB/Simulink/DPS Blockset Software Support 43

**Figure 11.** DPS Blockset on the web portal of The MathWorks Corporation

42 MATLAB – A Fundamental Tool for Scientific Computing and Engineering Applications – Volume 3

to the data file and continuously displayed on the *Scope* block.

performed by *DPS Robust IMC Control Synthesis* block.

domain. Thermocouple sensing junctions are located in the middle of the arm cross section, right above permanent thermocouples "P", except the node of the casting, where one temporal thermocouple points to the center of the node. Sensors S1 to S7 measure the temperature of cooling water in embedded chills. *Display* blocks of the sensors "P" and "W" in the scheme correspond to the real positions of sensors. Measured temperatures are saved

The HEATING subsystem has two basic operational regimes, which are activated by *Switch2* block. In the manual control regime, it is possible to set heating performance in range (0÷10) Volt and measure the temperature transient characteristics with P1 to P11 thermocouples in given locations of the mould. The second regime activated by *Switch2* block enables to perform controlled preheating of the casting mould to desired temperature profile W defined in 11 points, where thermocouples P1 to P11 are located. The control task is

**5. Distributed Parameter Systems Blockset for MATLAB & Simulink** 

CONNECTIONS of The MathWorks Corporation, Fig. 11. , (Hulkó et al., 2003-2010).

**DPS Space Synthesis** performs space synthesis as an approximation problem.

formulation and solution of control tasks for distributed parameter systems.

For the MATLAB & Simulink based software support of modelling, control and design of Distributed Parameter Systems given on complex 3D domains of definition, the programming environment **Distributed Parameter Systems Blockset for MATLAB & Simulink** (**DPS Blockset**) as Third-Party MathWorks Product has been developed by the Institute of Automation, Measurement and Applied Informatics, Faculty of Mechanical Engineering, Slovak University of Technology in Bratislava, within the program

The library of DPS Blockset shows Fig. 12. Blocks **HLDS** and **RHLDS** serve for modelling of distributed parameter systems as lumped-input/distributed-output systems with zero-order hold units. The block **DPS Control Synthesis** provides feedback to distributed parameter controlled systems in control loops with blocks for discrete-time **PID**, **Algebraic**, **State-Space** and **Robust Control**. The block **DPS Input** generates distributed quantities which can be used as distributed control quantities or distributed disturbances, etc. **DPS Display** presents distributed quantities with many options including export to AVI files. The block

The block **DPS Wizard** in step-by-step operation, by means of several model examples with default parameters on 1D-3D definition domains, gives an automatized guide for arrangement and setting distributed parameter control loops. The block **Demos** contains examples oriented to methodology of modelling and control synthesis. The block **Show**  contains motivation examples such as: *Control of temperature field of 3D metal body* (the controlled system was modelled in the virtual software environment COMSOL Multiphysics)**;** *Control of 3D beam of "smart" structure* (the controlled system was modelled in the virtual software environment ANSYS)*; Adaptive control of glass furnace* (the controlled system was modelled by Partial Differential Equations Toolbox of the MATLAB ), and *Groundwater remediation control* (the controlled system was modelled in the virtual software environment MODFLOW)**.** The block **Tutorial** presents methodological framework both for

**Figure 12.** The library of DPS Blockset for MATLAB & Simulink

## **6. Dynamics of the temperature field in the casting mould**

For control synthesis purpose, LDS/HLDS models of temperature fields in the casting mould have been created by means of evaluation of real-time experiments. The measurement of temperatures fields in the casting mould was performed by MATLAB & Simulink scheme *mould\_exp\_robust.mdl*, Fig. 10.

Robust Control of Distributed Parameter Systems with

*Sxs* , were determined points


y (m)

Temperature ( oC)


x (m)

7

4

Actuating in Zone # 5

2 3

10

1

11

*Sxs* are in Table 3.

(24)

97.87 4578.71 1214.73 3410.50

Demonstration in Casting Technology and MATLAB/Simulink/DPS Blockset Software Support 45


x (m)

7 4

**Figure 14.** Temperature profiles in the steady-state measured by thermocouples No. 1 - 11 after step

Actuating in Zone # 3

2

1

11

10

9


y (m)

Temperature ( oC)

 1,5 , *i ii <sup>i</sup> x xy* , located in each zone closely of lumped input quantities 1,5 *<sup>i</sup> <sup>i</sup> <sup>U</sup>* , where temperatures attain maximal amplitudes by actuating the heating power in each zone separately. These points for actuating of the heating power in each zone are represented by

> Zone No. 1 2 3 4 5 Thermocouple No. 7 9 1 4 11

Identification of measured dynamical characteristics of temperatures was performed in the MATLAB software environment, where graphical user interface (GUI) *ident* from System Identification Toolbox was activated. There after importing the time domain input/output data, from the pop-up menu *Process models* transfer function in the form (24) for identification has been chosen, see Fig. 15., where are also results of identification from actuating in zone #1. Comparison of measured and identified model output in zone #1 is

Continuous transfer functions with structure (24) are converted by means of function *zpk* to zero-pole-gain format (ZPK). Then, for control synthesis purposes, they are transformed to

*SH x z* with sample time T= 10 s.

 

1 1 1

> 49.36 5489.58 1398.20 3201.10

62.56 5340.97 973.90 3195.90

*Sxs* for actuating in each zone

 1 2

*K Tz s Tp s Tp s*

Zone No. 1 2 3 4 5

58.16 4945.00 1958.20 3301.80

change of the heating power separately in Zone #1, #3 and #5


x (m)

4

8

9 6 5


y (m)

150 7

1

11

Actuating in Zone # 1

2 3

10

0 50 100

Temperature ( oC)

positions of thermocouples given in Table 2.

discrete transfer functions 1,5 , *i i <sup>i</sup>*

*K Tp1 Tp2 Tz*

53.40 5881.68 939.66 3616.81

**Table 3.** Identified parameters of transfer functions 1,5 , *i i <sup>i</sup>*

For identification of i-th transfer functions 1,5 , *i i <sup>i</sup>*

**Table 2.** Position of thermocouples for identification of measured temperatures

presented in Fig. 16. Identified parameters of transfer functions 1,5 , *i i <sup>i</sup>*

#### **6.1. Experimental identification of transfer functions**

In the casting mould lumped inputs 1,5 *<sup>i</sup> <sup>i</sup> <sup>U</sup>* are heating elements which act on sub-domains 1,5 *<sup>i</sup> <sup>i</sup>* , (Zone 1 - 5). Distributed output is the temperature field of the casting mould.

Temperatures in the casting mould were measured by 11 thermocouples as a time-response to the step change of heating power, which was activated by the input voltage step from 0 to 2,5 V, for heating elements separately in each zone. Results of measurements for Zone #1, #3 and #5 are depicted in Fig. 13, where are time and *x-y* space dependences of temperatures and Fig. 14 presents temperature profiles in steady-states.

**Figure 13.** Time course of temperatures measured by thermocouples No. 1 - 11 after step change of the heating power separately in Zone #1, #3 and #5

#### Robust Control of Distributed Parameter Systems with Demonstration in Casting Technology and MATLAB/Simulink/DPS Blockset Software Support 45

44 MATLAB – A Fundamental Tool for Scientific Computing and Engineering Applications – Volume 3

For control synthesis purpose, LDS/HLDS models of temperature fields in the casting mould have been created by means of evaluation of real-time experiments. The measurement of temperatures fields in the casting mould was performed by MATLAB & Simulink scheme

In the casting mould lumped inputs 1,5 *<sup>i</sup> <sup>i</sup> <sup>U</sup>* are heating elements which act on sub-domains

Temperatures in the casting mould were measured by 11 thermocouples as a time-response to the step change of heating power, which was activated by the input voltage step from 0 to 2,5 V, for heating elements separately in each zone. Results of measurements for Zone #1, #3 and #5 are depicted in Fig. 13, where are time and *x-y* space dependences of temperatures

**Figure 13.** Time course of temperatures measured by thermocouples No. 1 - 11 after step change of the

y (m) Time (s)


Time (s) x (m)

1 2 3 8 7

<sup>0</sup> <sup>5000</sup> <sup>10000</sup> <sup>15000</sup>

50

50

100

Temperature ( oC)

150

100

Temperature ( oC)

150


<sup>0</sup> <sup>5000</sup> <sup>10000</sup> <sup>15000</sup>

<sup>0</sup> <sup>5000</sup> <sup>10000</sup> <sup>15000</sup>


Time (s) x (m)

1 2 3 8 7

y (m) Time (s)

Temperature ( oC)

Temperature ( oC)


<sup>0</sup> <sup>5000</sup> <sup>10000</sup> <sup>15000</sup>

Actuating in Zone # 5

11 10 9 6 5 4

Actuating in Zone # 3

11 10 9 6 5 4

1,5 *<sup>i</sup> <sup>i</sup>* , (Zone 1 - 5). Distributed output is the temperature field of the casting mould.

**6. Dynamics of the temperature field in the casting mould** 

**6.1. Experimental identification of transfer functions** 

and Fig. 14 presents temperature profiles in steady-states.

heating power separately in Zone #1, #3 and #5

y (m) Time (s)

<sup>0</sup> <sup>5000</sup> <sup>10000</sup> <sup>15000</sup>


Actuating in Zone # 1

11 10 9 6 5 4

0

Time (s) x (m)

1 2 3 8 7

<sup>5000</sup> <sup>10000</sup> <sup>15000</sup>

50

100

Temperature ( oC)

150


50

100

Temperature ( oC)

150

*mould\_exp\_robust.mdl*, Fig. 10.

**Figure 14.** Temperature profiles in the steady-state measured by thermocouples No. 1 - 11 after step change of the heating power separately in Zone #1, #3 and #5

For identification of i-th transfer functions 1,5 , *i i <sup>i</sup> Sxs* , were determined points 1,5 , *i ii <sup>i</sup> x xy* , located in each zone closely of lumped input quantities 1,5 *<sup>i</sup> <sup>i</sup> <sup>U</sup>* , where temperatures attain maximal amplitudes by actuating the heating power in each zone separately. These points for actuating of the heating power in each zone are represented by positions of thermocouples given in Table 2.


**Table 2.** Position of thermocouples for identification of measured temperatures

Identification of measured dynamical characteristics of temperatures was performed in the MATLAB software environment, where graphical user interface (GUI) *ident* from System Identification Toolbox was activated. There after importing the time domain input/output data, from the pop-up menu *Process models* transfer function in the form (24) for identification has been chosen, see Fig. 15., where are also results of identification from actuating in zone #1. Comparison of measured and identified model output in zone #1 is presented in Fig. 16. Identified parameters of transfer functions 1,5 , *i i <sup>i</sup> Sxs* are in Table 3. Continuous transfer functions with structure (24) are converted by means of function *zpk* to zero-pole-gain format (ZPK). Then, for control synthesis purposes, they are transformed to discrete transfer functions 1,5 , *i i <sup>i</sup> SH x z* with sample time T= 10 s.

$$\frac{K\{Tz\,s+1\}}{\left(Tp\_1\,s+1\right)\left(Tp\_2\,s+1\right)}\tag{24}$$


**Table 3.** Identified parameters of transfer functions 1,5 , *i i <sup>i</sup> Sxs* for actuating in each zone


Robust Control of Distributed Parameter Systems with

*Yxk* from each lumped input, 1,5 *<sup>i</sup> <sup>i</sup> <sup>U</sup>* , to

Demonstration in Casting Technology and MATLAB/Simulink/DPS Blockset Software Support 47

For uncertainty analysis, which takes place during the approximation problem solution, scheme from both, DPS Blockset blocks and Simulink blocks was arranged, see Fig. 17.

the corresponding output and approximated quantities corresponding to lumped output

**Figure 17.** Block scheme in the MATLAB & Simulink for the analysis of approximation in the space

uncertainties of their parameters. Analysis was performed in MATLAB environment, where functions from Robust Control Toolbox were used. Using functions *ureal* and *gridureal* were generated families of step responses with defined percentage variability of parameters in the

The gain variability K strongly affects to the value of the transient response in steady-state. The variability of the time constants Tp1 and Tp2 influences dynamics of transient response. Uncertainty region is caused by solution of the approximation task in the space synthesis. It would be appropriate to cover it by variability of nominal transfer function parameters. Cover of the uncertainty region in the time domain by variability of nominal parameters of

the transfer function , *i i Sxs* in Zone #1, #3 and #5 is presented in Fig. 20.

*Sxs* in the structure (24) were also analysed in terms of

*Sxs* . Results for zone 1 are depicted in Fig. 19.

from the block SS. Characteristics with uncertainty regions for actuating in each

There were obtained both, step responses *i i* , *<sup>i</sup>*

 *<sup>i</sup> <sup>i</sup> Y k*

zone are depicted in Fig. 18.

synthesis and step responses in Zones 1 – 5

Identified transfer function 1,5 , *i i <sup>i</sup>*

nominal transfer functions 1,5 , *i i <sup>i</sup>*

**Figure 15.** GUI Process Models menu and results of identification , *i i Sxs* in Zone #1

**Figure 16.** Comparison of measured and identified model output in Zone #1

#### **6.2. Uncertainty analysis in the space and time domain**

The proposed structure of the DPS control systems, Fig. 2, Fig. 5., are significant by decomposition of the control synthesis into the space and time subtasks. In the SS blocks, the approximation both of distributed controlled quantity *Yxk* , and reference quantity *W x* , , formulated as (8), (9) is solved. As the best approximation of controlled quantity *Yxk* , vector of optimal approximation parameters *<sup>i</sup> <sup>i</sup> Y k* is obtained. Dynamics of these lumped quantities is different in compare with lumped quantities, *i i* , *<sup>i</sup> Yxk* , given by transfer functions *i i* , *<sup>i</sup> SH x z* thus is created an uncertainty region in the time domain.

For uncertainty analysis, which takes place during the approximation problem solution, scheme from both, DPS Blockset blocks and Simulink blocks was arranged, see Fig. 17. There were obtained both, step responses *i i* , *<sup>i</sup> Yxk* from each lumped input, 1,5 *<sup>i</sup> <sup>i</sup> <sup>U</sup>* , to the corresponding output and approximated quantities corresponding to lumped output *<sup>i</sup> <sup>i</sup> Y k* from the block SS. Characteristics with uncertainty regions for actuating in each zone are depicted in Fig. 18.

46 MATLAB – A Fundamental Tool for Scientific Computing and Engineering Applications – Volume 3

**Figure 15.** GUI Process Models menu and results of identification , *i i Sxs* in Zone #1

**Figure 16.** Comparison of measured and identified model output in Zone #1

measured model output

The proposed structure of the DPS control systems, Fig. 2, Fig. 5., are significant by decomposition of the control synthesis into the space and time subtasks. In the SS blocks, the approximation both of distributed controlled quantity *Yxk* , and reference quantity *W x* , , formulated as (8), (9) is solved. As the best approximation of controlled quantity



Amplitude

0

1

2

of these lumped quantities is different in compare with lumped quantities, *i i* , *<sup>i</sup>*

*Y k*

*SH x z* thus is created an uncertainty region in the time

is obtained. Dynamics

<sup>0</sup> <sup>5000</sup> <sup>10000</sup> <sup>15000</sup> -3

Measured minus simulated model output

Time (s)

*Yxk* ,

**6.2. Uncertainty analysis in the space and time domain** 

<sup>0</sup> <sup>5000</sup> <sup>10000</sup> <sup>15000</sup> <sup>0</sup>

Time (s)

Measured and simulated model output

*Yxk* , vector of optimal approximation parameters *<sup>i</sup> <sup>i</sup>*

given by transfer functions *i i* , *<sup>i</sup>*

domain.

Amplitude

**Figure 17.** Block scheme in the MATLAB & Simulink for the analysis of approximation in the space synthesis and step responses in Zones 1 – 5

Identified transfer function 1,5 , *i i <sup>i</sup> Sxs* in the structure (24) were also analysed in terms of uncertainties of their parameters. Analysis was performed in MATLAB environment, where functions from Robust Control Toolbox were used. Using functions *ureal* and *gridureal* were generated families of step responses with defined percentage variability of parameters in the nominal transfer functions 1,5 , *i i <sup>i</sup> Sxs* . Results for zone 1 are depicted in Fig. 19.

The gain variability K strongly affects to the value of the transient response in steady-state. The variability of the time constants Tp1 and Tp2 influences dynamics of transient response. Uncertainty region is caused by solution of the approximation task in the space synthesis. It would be appropriate to cover it by variability of nominal transfer function parameters. Cover of the uncertainty region in the time domain by variability of nominal parameters of the transfer function , *i i Sxs* in Zone #1, #3 and #5 is presented in Fig. 20.

Robust Control of Distributed Parameter Systems with

1000

Amplitude

1500

Tp2 (s)

2000

<sup>3000</sup> <sup>4000</sup> <sup>5000</sup> <sup>6000</sup> <sup>7000</sup> <sup>500</sup>

Zone # 5, Ranges of Tp1, Tp2

Ranges Nominal

Tp1 (s)

Zone # 5, Ranges of Tp1, Tp2

<sup>0</sup> <sup>5000</sup> <sup>10000</sup> <sup>15000</sup> <sup>0</sup>

Nominal Ranges Approximation

Time (sec)

Demonstration in Casting Technology and MATLAB/Simulink/DPS Blockset Software Support 49

**Figure 20.** Cover of the uncertainty region in the time domain by variability of the nominal parameters

<sup>0</sup> <sup>5000</sup> <sup>10000</sup> <sup>15000</sup> <sup>0</sup>

Time (sec)

<sup>3000</sup> <sup>4000</sup> <sup>5000</sup> <sup>6000</sup> <sup>7000</sup> <sup>500</sup>

Zone # 3, Ranges of Tp1, Tp2

Tp1 (s)

Zone # 3, Ranges of Tp1, Tp2

Robust control of temperature fields of the casting mould as distributed parameter system was performed with MATLAB & Simulink and DPS Blockset software support. Robust controllers designed by the IMC control strategy were first optimised through their tuning parameters *<sup>i</sup><sup>i</sup>* and then implemented to the control scheme for the real-time control of

In the MATLAB & Simulink environment, by means of the DPS Blockset, distributed parameter system of robust control, *mould\_robust\_DZPK.mdl* was arranged, see Fig. 21. It is DPS feedback control loop, where the *DPS Robust Control Synthesis* block includes both, time and spatial part of the control synthesis, see Fig. 22 a). In this case, the control system consists of five single parameter control loops, each for one zone of the mould, where

*DPS Robust Control Synthesis* contains two blocks named *DPS Space Synthesis*, where approximation of distributed controlled quantity *Yxk* , and reference quantity, *W x* , is executed. The time control synthesis is performed by *Robust controllers based on IMC* block with five discrete controllers given by ZPK transfer function and filters with parameters

Tp1 and Tp2 of the transfer function , *i i Sxs* in Zone #1, #3 and #5

1000

Amplitude

1500

Tp2 (s)

2000

discrete robust controllers based on IMC structure are used.

temperature fields of the casting mould in the benchmark casting plant.

**7. Robust control process** 

<sup>0</sup> <sup>5000</sup> <sup>10000</sup> <sup>15000</sup> <sup>0</sup>

Time (sec)

<sup>3000</sup> <sup>4000</sup> <sup>5000</sup> <sup>6000</sup> <sup>7000</sup> <sup>500</sup>

Zone # 1, Ranges of Tp1, Tp2

Tp1 (s)

Zone # 1, Ranges of Tp1, Tp2

1000

Amplitude

1500

Tp2 (s)

2000

1,5 *<sup>i</sup> <sup>i</sup>* , see Fig. 22 b).

**7.1. Optimization of tuning parameters** 

**Figure 18.** Approximation parameters *<sup>i</sup> <sup>i</sup> Y k* as a result of DPS space synthesis and step responses *i i* , *<sup>i</sup> HH xk* for actuating in Zones 1 – 5

**Figure 19.** Family (Samples) of step responses for the defined percentage variability of parameters in the nominal transfer function <sup>1</sup> , *<sup>i</sup> S xs*

Robust Control of Distributed Parameter Systems with Demonstration in Casting Technology and MATLAB/Simulink/DPS Blockset Software Support 49

**Figure 20.** Cover of the uncertainty region in the time domain by variability of the nominal parameters Tp1 and Tp2 of the transfer function , *i i Sxs* in Zone #1, #3 and #5

## **7. Robust control process**

48 MATLAB – A Fundamental Tool for Scientific Computing and Engineering Applications – Volume 3

Amplitude

<sup>0</sup> <sup>5000</sup> <sup>10000</sup> <sup>15000</sup> <sup>0</sup>

Time (s)

Amplitude

Actuating in Zone # 2

*Y k*

<sup>0</sup> <sup>5000</sup> <sup>10000</sup> <sup>15000</sup> <sup>0</sup>

10% variability of K, Actuating in Zone # 1

0 0.5 1 1.5 2 2.5 3

Time (sec)

20% variability of Tp1, Actuating in Zone # 1

0 0.5 1 1.5 2 2.5 3

Time (sec)

Time (s)

Actuating in Zone # 4

**Figure 19.** Family (Samples) of step responses for the defined percentage variability of parameters in

x 104

 Nominal Samples

x 104

Amplitude

Amplitude

 Nominal Samples

as a result of DPS space synthesis and step responses

10% variability of Tz, Actuating in Zone # 1

0 0.5 1 1.5 2 2.5 3

Time (sec)

20% variability of Tp2, Actuating in Zone # 1

0 0.5 1 1.5 2 2.5 3

Time (sec)

x 104

x 104

 Nominal Samples

 Nominal Samples

<sup>0</sup> <sup>5000</sup> <sup>10000</sup> <sup>15000</sup> <sup>0</sup>

Time (s)

Amplitude

Actuating in Zone # 5

<sup>0</sup> <sup>5000</sup> <sup>10000</sup> <sup>15000</sup> <sup>0</sup>

 Approximation Step Response

Time (s)

Actuating in Zone # 3

**Figure 18.** Approximation parameters *<sup>i</sup> <sup>i</sup>*

*HH xk* for actuating in Zones 1 – 5

<sup>0</sup> <sup>5000</sup> <sup>10000</sup> <sup>15000</sup> <sup>0</sup>

Time (s)

Amplitude

Actuating in Zone # 1

the nominal transfer function <sup>1</sup> , *<sup>i</sup> S xs*

Amplitude

Amplitude

*i i* , *<sup>i</sup>*

Amplitude

Robust control of temperature fields of the casting mould as distributed parameter system was performed with MATLAB & Simulink and DPS Blockset software support. Robust controllers designed by the IMC control strategy were first optimised through their tuning parameters *<sup>i</sup><sup>i</sup>* and then implemented to the control scheme for the real-time control of temperature fields of the casting mould in the benchmark casting plant.

## **7.1. Optimization of tuning parameters**

In the MATLAB & Simulink environment, by means of the DPS Blockset, distributed parameter system of robust control, *mould\_robust\_DZPK.mdl* was arranged, see Fig. 21. It is DPS feedback control loop, where the *DPS Robust Control Synthesis* block includes both, time and spatial part of the control synthesis, see Fig. 22 a). In this case, the control system consists of five single parameter control loops, each for one zone of the mould, where discrete robust controllers based on IMC structure are used.

*DPS Robust Control Synthesis* contains two blocks named *DPS Space Synthesis*, where approximation of distributed controlled quantity *Yxk* , and reference quantity, *W x* , is executed. The time control synthesis is performed by *Robust controllers based on IMC* block with five discrete controllers given by ZPK transfer function and filters with parameters 1,5 *<sup>i</sup> <sup>i</sup>* , see Fig. 22 b).

Robust Control of Distributed Parameter Systems with

Demonstration in Casting Technology and MATLAB/Simulink/DPS Blockset Software Support 51

(25)

The block *DPS Robust Control Synthesis* also contains block *Output Constraint for optimization of alpha*, where optimization of parameters 1,5 *<sup>i</sup> <sup>i</sup>* according to criterion function (25) is

min , ,

Parameters of filters were optimized in the presence of constraints of the criterion function in order to assure nearly aperiodic course of the quadratic norm of the distributed control error with respect to the robust stability and robust performance conditions, see Fig. 23 and optimization progress presents Fig. 24. Control process was simulated for the reference quantities - temperatures on given 11 positions, where thermocouples are embedded. Robust stability was tested for ranges around the nominal parameters Tp1 and Tp2, of the transfer function , *i i Sxs* , see Fig. 20, with MATLAB function *robuststab*, from the Robust Control Toolbox, e.g. in the following is the printout of the robust stability testing for the

*J W xk Yxk*

0

*k*

STABreport = Uncertain System is robustly stable to modeled uncertainty.

**Figure 23.** Minimization of the criterion function in presence of output constraints


*N*

*i*

performed.

control loop in zone 1.

Robust stability testing in Zone #1 StabilityMargin = UpperBound: 1.3333 LowerBound: 1.3296 DestabilizingFrequency: 8.0745e-004

DestbUnc = Tp1: 1.7210 Tp2: 438.5898

**Figure 21.** DPS feedback robust control system in the DPS Blockset for simulation of robust control of the temperature fields in the casting mould

**Figure 22.** Structure of blocks in the scheme mould\_robust\_DZPK.mdl: a) DPS Robust Control Synthesis, b) Robust controllers based on IMC

The block *DPS Robust Control Synthesis* also contains block *Output Constraint for optimization of alpha*, where optimization of parameters 1,5 *<sup>i</sup> <sup>i</sup>* according to criterion function (25) is performed.

$$J = \min\_{a\_i} \sum\_{k=0}^{N} \left\| \left. W(\overline{\boldsymbol{x}}, k) - Y(\overline{\boldsymbol{x}}, k) \right\| \right\| \tag{25}$$

Parameters of filters were optimized in the presence of constraints of the criterion function in order to assure nearly aperiodic course of the quadratic norm of the distributed control error with respect to the robust stability and robust performance conditions, see Fig. 23 and optimization progress presents Fig. 24. Control process was simulated for the reference quantities - temperatures on given 11 positions, where thermocouples are embedded. Robust stability was tested for ranges around the nominal parameters Tp1 and Tp2, of the transfer function , *i i Sxs* , see Fig. 20, with MATLAB function *robuststab*, from the Robust Control Toolbox, e.g. in the following is the printout of the robust stability testing for the control loop in zone 1.

Robust stability testing in Zone #1 StabilityMargin = UpperBound: 1.3333 LowerBound: 1.3296 DestabilizingFrequency: 8.0745e-004

DestbUnc = Tp1: 1.7210

50 MATLAB – A Fundamental Tool for Scientific Computing and Engineering Applications – Volume 3

**Figure 21.** DPS feedback robust control system in the DPS Blockset for simulation of robust control of

11

11

11

5 11

11

11

2 W

> 1 ei

1 Y Discrete robust controllers based on IMC

> 1 ui

**b)**






kcon{5}\*zcon{5}(z) pcon{5}(z) Discrete Zero-Pole5

kcon{4}\*zcon{4}(z) pcon{4}(z) Discrete Zero-Pole4

kcon{3}\*zcon{3}(z) pcon{3}(z) Discrete Zero-Pole3

5 5

kcon{2}\*zcon{2}(z) pcon{2}(z) Discrete Zero-Pole2

kcon{1}\*zcon{1}(z) pcon{1}(z) Discrete Zero-Pole1

11

DPS Space Synthesis

5

DPS Space Synthesis

sfun\_norm S-Function1

Output Constraint for Optimization of alpha 0.9619 0.9769 0.9565 0.9457 0.94 Filter parameters for robust controllers

**Figure 22.** Structure of blocks in the scheme mould\_robust\_DZPK.mdl: a) DPS Robust Control

the temperature fields in the casting mould

5

5 5

e\_approx To Workspace1

eiui Robust controllers based on IMC

Saturation

5 5

1 U

alpha

5

Synthesis, b) Robust controllers based on IMC

**a)**

dps\_norm To Workspace

5

Tp2: 438.5898

STABreport = Uncertain System is robustly stable to modeled uncertainty.


**Figure 23.** Minimization of the criterion function in presence of output constraints

Robust Control of Distributed Parameter Systems with

2 W

0

Time (s) y (m)

1 2 3 8 7

Control quantities

<sup>0</sup> <sup>2000</sup> <sup>4000</sup> <sup>6000</sup> <sup>8000</sup> <sup>10000</sup> <sup>12000</sup> <sup>0</sup>

Time (s)

5000

u1 u2 u3 u4 u5 10000

Controlled quantities

1 Y

Demonstration in Casting Technology and MATLAB/Simulink/DPS Blockset Software Support 53

K\*uvec

K\*uvec

Space Synthesis Y

Space Synthesis W


ui (V)

Controlled quantities

Temperature ( oC)

0

Time (s) x (m)

5000

Temperature ( oC)

**DPS Time and Space Control Synthesis**

Add

11 10 9 6 5 4

**Figure 25.** Structure of the block DPS Robust IMC Control Synthesis in the scheme

eiui

Robust controllers based on IMC

**Figure 26.** Robust control of the temperature field of the casting mould: distributed reference quantities, controlled quantities at given position of the mould, quadratic norm of distributed control

*U k* in Zones 1-5

mould\_exp\_robust.mdl for real-time control


Quadratic norm of Control error

x (m)

**7**

**5**

**6 9 3**

Distributed reference quantities W(x,y)

**2**

**10**

**8**

**4**

y (m)

Amplitude

Temperature ( oC)

**1**

**11**

1 U

error *Exk* , and control quantities 1,5 *<sup>i</sup> <sup>i</sup>*

<sup>0</sup> <sup>2000</sup> <sup>4000</sup> <sup>6000</sup> <sup>8000</sup> <sup>10000</sup> <sup>12000</sup> <sup>0</sup>

Time (s)


**Figure 24.** Optimization progress for tuning parameters *<sup>i</sup> <sup>i</sup>*1,5 of robust controllers


'Tp1' is 100%. Increasing 'Tp1' by 25% leads to a 25% decrease in the margin.

'Tp2' is 54%. Increasing 'Tp2' by 25% leads to a 14% decrease in the margin.

## **7.2. Real-time robust control of preheating in the casting mould**

Real-time robust control of the temperature fields of the casting mould in the benchmark casting plant was performed by means of Simulink designed block scheme *mould\_exp\_robust.mdl*. The structure of the block *DPS Robust IMC Control Synthesis* is depicted in Fig. 25, where are blocks *Space Synthesis Y* and *Space Synthesis W*, and block *Robust controllers based on IMC* with the same structure as in Fig. 22 b).

#### **DPS Time and Space Control Synthesis** 1 U K\*uvec Space Synthesis Y K\*uvec Space Synthesis W eiui Robust controllers based on IMC Add 2 W 1 Y

**Figure 25.** Structure of the block DPS Robust IMC Control Synthesis in the scheme mould\_exp\_robust.mdl for real-time control

52 MATLAB – A Fundamental Tool for Scientific Computing and Engineering Applications – Volume 3

**Figure 24.** Optimization progress for tuning parameters *<sup>i</sup> <sup>i</sup>*1,5 of robust controllers


**7.2. Real-time robust control of preheating in the casting mould** 

*Robust controllers based on IMC* with the same structure as in Fig. 22 b).

 'Tp1' is 100%. Increasing 'Tp1' by 25% leads to a 25% decrease in the margin. 'Tp2' is 54%. Increasing 'Tp2' by 25% leads to a 14% decrease in the margin.

Real-time robust control of the temperature fields of the casting mould in the benchmark casting plant was performed by means of Simulink designed block scheme *mould\_exp\_robust.mdl*. The structure of the block *DPS Robust IMC Control Synthesis* is depicted in Fig. 25, where are blocks *Space Synthesis Y* and *Space Synthesis W*, and block

causing an instability at 0.000807 rad/s.


**Figure 26.** Robust control of the temperature field of the casting mould: distributed reference quantities, controlled quantities at given position of the mould, quadratic norm of distributed control error *Exk* , and control quantities 1,5 *<sup>i</sup> <sup>i</sup> U k* in Zones 1-5

Control process was performed for the reference variable - temperature profile given in 11 positions, where thermocouples are embedded. Results of the real-time robust control process are on Fig. 26. The quality of control both in the time and space domain is given by the quadratic norm of the distributed control error.

Robust Control of Distributed Parameter Systems with

Demonstration in Casting Technology and MATLAB/Simulink/DPS Blockset Software Support 55

Belavý, C. et al. (2009). Uncertainty Analysis and Robust Control of Temperature Fields of Casting Mould as Distributed Parameter Systems, *INCOM ´09: Preprints of the 13th IFAC symposium on Information control problems in manufacturing,* Moscow, June 3-5,

Butkovskij, A. G. (1965). *Optimal Control of Distributed Parameter Systems,* Nauka,

Curtain, R. & Morris K. (2009). Transfer Functions of Distributed Parameter Systems:

Hulkó, G. et al. (1981). On Adaptive Control of Distributed Parameter Systems,

Hulkó, G. et al. (1987). Control of Distributed Parameter Systems by means of Multi-Input and Multi-Distributed-Output Systems*, Proceedings of 10-th World Congress of IFAC*,

Hulkó, G. et al. (1998). *Modeling, Control and Design of Distributed Parameter Systems with Demonstrations in MATLAB,* Publishing House STU, ISBN 80-227-1083-0,

Hulkó, G. et al (2003-2007). Distributed Parameter Systems. Web portal [Online]. Available:

Hulkó, G. et al. (2009a). Engineering Methods and Software Support for Modelling and Design of Discrete-time Control of Distributed Parameter Systems*. European Journal of Control,* Vol. 15, No. Iss. 3-4, *Fundamental Issues in Control,* (May-August 2009), pp. 407-

Hulkó, G. et al. (2009b). Engineering Methods and Software Support for Control of Distributed Parameter Systems*, ASC 2009: Proceedings of the 7th Asian Control Conference,* 

Hulkó, G. et al. (2003-2010). *Distributed Parameter Systems Blockset for MATLAB & Simulink*, www.mathworks.com/products/connections/ - Third-Party Product of The

Morari, M. & Zafirou, E. (1989). *Robust Process Control,* Prentice Hall, Englewood

Lasiecka, I. & Triggiani, R. (2000). *Control Theory for Partial Differential Equations: Continuous and Approximation Theories I. Abstract Parabolic Systems,* Cambridge University Press,

Li, H. X. & Qi, Ch. (2010). Modeling of Distributed Parameter Systems for Application – A Synthesized Review from Time-Space Separation*. Journal of Process Control*, No 20,

Lions, J. L. (1971). *Optimal Control of Systems Governed by Partial Differential Equations,*

MathWorks, Bratislava-Natick, Available from: www.dpscontrol.sk

A tutorial. *Automatica,* Vol. 45, (2009), pp. 1101-1116

*Proceedings of 8-th World Congress of IFAC*, Kyoto, 1981

**9. References** 

2009

Moscow

Munich, 1987

Bratislava

Cliffs

www.dpscontrol.sk

417, ISSN 0947-3580

(2010), pp. 891-901

Honk Kong, China, August 27-29, 2009

ISBN 0-521-43408-4, Cambridge, UK

Springer - Verlag, Berlin - Heidelberg - New York

## **8. Conclusion**

The aim of this chapter was to present the engineering approach for the robust control of DPS, which opens a wide space for novel applications of the toolboxes and blocksets of the MATLAB & Simulink software environment. This approach is based on the general decomposition of controlled DPS dynamics, represented by transient and impulse characteristics, into time and space components. Starting out from this dynamics decomposition a methodical framework was presented for the decomposition of control synthesis into the space and time subtasks. In the space domain an approximation problems were solved, while in the time domain the control synthesis was performed by lumped parameter SISO control loops, where various wellknown methods for design of controller is possible to utilize. The advantage of this approach is the relatively simple LDS model of DPS, which is directly suitable for control purposes and can be easily identified from input-output data by means of classical techniques.

Currently, it is interesting to formulate and solve tasks of control in various engineering branches, including the casting technology, by means of methods and tools of distributed parameter systems. Methodical approach presented in this chapter demonstrates simple possibilities, how to exploit the distributed dynamical characteristics on complex definition domains, obtained by evaluation of measured data for robust IMC control synthesis of DPS with respect to uncertainty of models and the real-time control according to technological requirements.

## **Author details**

Cyril Belavý, Gabriel Hulkó and Karol Ondrejkovič *Institute of Automation, Measurement and Applied Informatics, Faculty of Mechanical Engineering, Center for Control of Distributed Parameter Systems, Slovak University of Technology in Bratislava, Slovak Republic* 

## **Acknowledgement**

This work was supported by the Slovak Scientific Grant Agency VEGA under the contract No. 1/0138/11 and the Slovak Research and Development Agency under contracts No. APVV-0131-10 and No. APVV-0090-10 and also by the European Union with the co-financing of the European Social Fund, grant TAMOP-4.2.1.B-11/2/KMR-2011- 0001.

### **9. References**

54 MATLAB – A Fundamental Tool for Scientific Computing and Engineering Applications – Volume 3

the quadratic norm of the distributed control error.

**8. Conclusion** 

classical techniques.

requirements.

**Author details** 

**Acknowledgement** 

0001.

Cyril Belavý, Gabriel Hulkó and Karol Ondrejkovič

*Institute of Automation, Measurement and Applied Informatics,* 

*Slovak University of Technology in Bratislava, Slovak Republic* 

*Faculty of Mechanical Engineering, Center for Control of Distributed Parameter Systems,* 

This work was supported by the Slovak Scientific Grant Agency VEGA under the contract No. 1/0138/11 and the Slovak Research and Development Agency under contracts No. APVV-0131-10 and No. APVV-0090-10 and also by the European Union with the co-financing of the European Social Fund, grant TAMOP-4.2.1.B-11/2/KMR-2011-

Control process was performed for the reference variable - temperature profile given in 11 positions, where thermocouples are embedded. Results of the real-time robust control process are on Fig. 26. The quality of control both in the time and space domain is given by

The aim of this chapter was to present the engineering approach for the robust control of DPS, which opens a wide space for novel applications of the toolboxes and blocksets of the MATLAB & Simulink software environment. This approach is based on the general decomposition of controlled DPS dynamics, represented by transient and impulse characteristics, into time and space components. Starting out from this dynamics decomposition a methodical framework was presented for the decomposition of control synthesis into the space and time subtasks. In the space domain an approximation problems were solved, while in the time domain the control synthesis was performed by lumped parameter SISO control loops, where various wellknown methods for design of controller is possible to utilize. The advantage of this approach is the relatively simple LDS model of DPS, which is directly suitable for control purposes and can be easily identified from input-output data by means of

Currently, it is interesting to formulate and solve tasks of control in various engineering branches, including the casting technology, by means of methods and tools of distributed parameter systems. Methodical approach presented in this chapter demonstrates simple possibilities, how to exploit the distributed dynamical characteristics on complex definition domains, obtained by evaluation of measured data for robust IMC control synthesis of DPS with respect to uncertainty of models and the real-time control according to technological


Wang, P. K. C. (1964). *Control of Distributed Parameter Systems* (Advances in Control Systems: Theory and Applications, 1.), Academic Press, New York

**Chapter 3** 

© 2012 Ibrahim, licensee InTech. This is an open access chapter distributed under the terms of the Creative Commons Attribution License (http://creativecommons.org/licenses/by/3.0), which permits unrestricted use,

distribution, and reproduction in any medium, provided the original work is properly cited.

© 2012 Ibrahim, licensee InTech. This is a paper distributed under the terms of the Creative Commons Attribution License (http://creativecommons.org/licenses/by/3.0), which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.

Fouling is generally defined as the deposition and accumulation of unwanted materials such as scale, algae, suspended solids and insoluble salts on the internal or external surfaces of processing equipment including boilers and heat exchangers (Fig 1). Heat exchangers are process equipment in which heat is continuously or semi-continuously transferred from a hot to a cold fluid directly or indirectly through a heat transfer surface that separates the two fluids. Heat exchangers consist primarily of bundles of pipes, tubes or plate coils.

**Fouling in Heat Exchangers** 

Additional information is available at the end of the chapter

Hassan Al-Haj Ibrahim

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

**Figure 1.** Fouling of heat exchangers.

**1. Introduction** 

## **Chapter 3**

## **Fouling in Heat Exchangers**

## Hassan Al-Haj Ibrahim

56 MATLAB – A Fundamental Tool for Scientific Computing and Engineering Applications – Volume 3

Theory and Applications, 1.), Academic Press, New York

Wang, P. K. C. (1964). *Control of Distributed Parameter Systems* (Advances in Control Systems:

Additional information is available at the end of the chapter

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

## **1. Introduction**

Fouling is generally defined as the deposition and accumulation of unwanted materials such as scale, algae, suspended solids and insoluble salts on the internal or external surfaces of processing equipment including boilers and heat exchangers (Fig 1). Heat exchangers are process equipment in which heat is continuously or semi-continuously transferred from a hot to a cold fluid directly or indirectly through a heat transfer surface that separates the two fluids. Heat exchangers consist primarily of bundles of pipes, tubes or plate coils.

**Figure 1.** Fouling of heat exchangers.

© 2012 Ibrahim, licensee InTech. This is an open access chapter distributed under the terms of the Creative Commons Attribution License (http://creativecommons.org/licenses/by/3.0), which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited. © 2012 Ibrahim, licensee InTech. This is a paper distributed under the terms of the Creative Commons Attribution License (http://creativecommons.org/licenses/by/3.0), which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.

Fouling on process equipment surfaces can have a significant, negative impact on the operational efficiency of the unit. On most industries today, a major economic drain may be caused by fouling. The total fouling related costs for major industrialised nations is estimated to exceed US\$4.4 milliard annually. One estimate puts the losses due to fouling of heat exchangers in industrialised nations to be about 0.25% to 30% of their GDP [1, 2]. According to Pritchard and Thackery (Harwell Laboratories), about 15% of the maintenance costs of a process plant can be attributed to heat exchangers and boilers, and of this, half is probably caused by fouling. Costs associated with heat exchanger fouling include production losses due to efficiency deterioration and to loss of production during planned or unplanned shutdowns due to fouling, and maintenance costs resulting from the removal of fouling deposits with chemicals and/or mechanical antifouling devices or the replacement of corroded or plugged equipment. Typically, cleaning costs are in the range of \$40,000 to \$50,000 per heat exchanger per cleaning.

Fouling in Heat Exchangers 59

heat flux is high, as in steam generators, fouling can lead to local hot spots resulting ultimately in mechanical failure of the heat transfer surface. Such effects lead in most cases

Loss of heat transfer and subsequent charge outlet temperature decrease is a result of the low thermal conductivity of the fouling layer or layers which is generally lower than the thermal conductivity of the fluids or conduction wall. As a result of this lower thermal conductivity, the overall thermal resistance to heat transfer is increased and the effectiveness and thermal efficiency of heat exchangers are reduced. A simple way to monitor a heat transfer system is to plot the outlet temperature versus time. In one unit at an oil refinery, in Homs, Syria, fouling led to a feed temperature decrease from 210C to 170C. In order to bring the feed to the required temperature, the heat duty of the furnace may have to be increased with additional fuel required and resulting increased fuel cost. Alternatively, the heat exchanger surface area may have to be increased with consequent additional installation and maintenance costs. The required excess surface area may vary between 10-50%, with an average around 35%, and the additional extra costs involved may add up to a staggering 2.5 to 3.0 times the initial purchase price of the heat

With the onset of fouling and the consequent build up of fouling layer or layers, the cross sectional area of tubes or flow channels is reduced. In addition, increased surface roughness due to fouling will increase frictional resistance to flow. Such effects inevitably lead to an increase in the pressure drop across the heat exchanger, which is required to maintain the flow rate through the exchanger, and may even lead to flow blocks. Experience with pressure drop monitoring has shown, however, that it is not usually as sensitive an indicator of the early onset of fouling when compared to heat transfer data; thus pressure drop is not commonly used for crude preheat monitoring. In situations where significant swings in flow rates are experienced, flow correction can be applied to both pressure drop and to heat

Different fouling deposit structures can lead to under-deposit corrosion of the substrate material such as localised fouling, deposit tubercles and sludge piles. The factors that are most likely to influence the probability of under-deposit corrosion include deposit composition and its porosity and permeability. Even minor components of the deposits can sometimes cause severe corrosion of the underlying metal such as the hot corrosion caused

Fouling is responsible for the emission of many millions of tonnes of carbon dioxide as well as the use and disposal of hazardous cleaning chemicals. Data from oil refineries suggest that crude oil fouling accounts for about 10% of the total CO2 emission of these plants. Wastes generated from the cleaning of heat exchangers may contain hazardous wastes such as lead and chromium, although some refineries which do not produce leaded gasoline and which use non-chrome corrosion inhibitors typically do not generate sludge that contains

these constituents. Oily wastewater is also generated during heat exchanger cleaning.

to production losses and increased maintenance costs.

transfer calculations to normalise the data to a standard flow.

by vanadium in the deposits of fired boilers [5].

exchangers.

Fouling in heat exchangers is not a new problem. In fact, fouling has been recognised for a long time, and research on heat exchanger fouling was conducted as early as 1910 and the first practical application of this research was implemented in the 1920's. Technological progress in prevention, mitigation and removal techniques in industrial fouling was investigated in a study conducted at the Battelle Pacific Northwest Laboratories for the U.S. Department of Energy. Two hundred and thirty one patents relevant to fouling were analysed [3]. Furthermore, great technical advance in the design and manufacture of heat exchangers has in the meantime been achieved. Nonetheless, heat exchanger fouling remains today one of the major unresolved problems in Thermal Science, and prevention or mitigation of the fouling problem is still an ongoing process. Further research on the problem of fouling in heat exchangers and practical methods for predicting the fouling factor, making use in particular of modern digital techniques, are still called for. One significant and clear indication of the relevance and urgency of the problem may be seen in the current international patent activity on fouling (Table 1).


**Table 1.** International Patent Activity [4]

Major detrimental effects of fouling include loss of heat transfer as indicated by charge outlet temperature decrease and pressure drop increase. Other detrimental effects of fouling may also include blocked process pipes, under-deposit corrosion and pollution. Where the heat flux is high, as in steam generators, fouling can lead to local hot spots resulting ultimately in mechanical failure of the heat transfer surface. Such effects lead in most cases to production losses and increased maintenance costs.

58 MATLAB – A Fundamental Tool for Scientific Computing and Engineering Applications – Volume 3

\$50,000 per heat exchanger per cleaning.

**Table 1.** International Patent Activity [4]

the current international patent activity on fouling (Table 1).

Fouling on process equipment surfaces can have a significant, negative impact on the operational efficiency of the unit. On most industries today, a major economic drain may be caused by fouling. The total fouling related costs for major industrialised nations is estimated to exceed US\$4.4 milliard annually. One estimate puts the losses due to fouling of heat exchangers in industrialised nations to be about 0.25% to 30% of their GDP [1, 2]. According to Pritchard and Thackery (Harwell Laboratories), about 15% of the maintenance costs of a process plant can be attributed to heat exchangers and boilers, and of this, half is probably caused by fouling. Costs associated with heat exchanger fouling include production losses due to efficiency deterioration and to loss of production during planned or unplanned shutdowns due to fouling, and maintenance costs resulting from the removal of fouling deposits with chemicals and/or mechanical antifouling devices or the replacement of corroded or plugged equipment. Typically, cleaning costs are in the range of \$40,000 to

Fouling in heat exchangers is not a new problem. In fact, fouling has been recognised for a long time, and research on heat exchanger fouling was conducted as early as 1910 and the first practical application of this research was implemented in the 1920's. Technological progress in prevention, mitigation and removal techniques in industrial fouling was investigated in a study conducted at the Battelle Pacific Northwest Laboratories for the U.S. Department of Energy. Two hundred and thirty one patents relevant to fouling were analysed [3]. Furthermore, great technical advance in the design and manufacture of heat exchangers has in the meantime been achieved. Nonetheless, heat exchanger fouling remains today one of the major unresolved problems in Thermal Science, and prevention or mitigation of the fouling problem is still an ongoing process. Further research on the problem of fouling in heat exchangers and practical methods for predicting the fouling factor, making use in particular of modern digital techniques, are still called for. One significant and clear indication of the relevance and urgency of the problem may be seen in

Country No. of Patents % of Patents U.S.A. 147 63.6 Germany 22 9.5 Japan 21 9.1 Sweden 9 3.9 Switzerland 8 3.5 Other 24 10.4 Total 231 100.0

Major detrimental effects of fouling include loss of heat transfer as indicated by charge outlet temperature decrease and pressure drop increase. Other detrimental effects of fouling may also include blocked process pipes, under-deposit corrosion and pollution. Where the Loss of heat transfer and subsequent charge outlet temperature decrease is a result of the low thermal conductivity of the fouling layer or layers which is generally lower than the thermal conductivity of the fluids or conduction wall. As a result of this lower thermal conductivity, the overall thermal resistance to heat transfer is increased and the effectiveness and thermal efficiency of heat exchangers are reduced. A simple way to monitor a heat transfer system is to plot the outlet temperature versus time. In one unit at an oil refinery, in Homs, Syria, fouling led to a feed temperature decrease from 210C to 170C. In order to bring the feed to the required temperature, the heat duty of the furnace may have to be increased with additional fuel required and resulting increased fuel cost. Alternatively, the heat exchanger surface area may have to be increased with consequent additional installation and maintenance costs. The required excess surface area may vary between 10-50%, with an average around 35%, and the additional extra costs involved may add up to a staggering 2.5 to 3.0 times the initial purchase price of the heat exchangers.

With the onset of fouling and the consequent build up of fouling layer or layers, the cross sectional area of tubes or flow channels is reduced. In addition, increased surface roughness due to fouling will increase frictional resistance to flow. Such effects inevitably lead to an increase in the pressure drop across the heat exchanger, which is required to maintain the flow rate through the exchanger, and may even lead to flow blocks. Experience with pressure drop monitoring has shown, however, that it is not usually as sensitive an indicator of the early onset of fouling when compared to heat transfer data; thus pressure drop is not commonly used for crude preheat monitoring. In situations where significant swings in flow rates are experienced, flow correction can be applied to both pressure drop and to heat transfer calculations to normalise the data to a standard flow.

Different fouling deposit structures can lead to under-deposit corrosion of the substrate material such as localised fouling, deposit tubercles and sludge piles. The factors that are most likely to influence the probability of under-deposit corrosion include deposit composition and its porosity and permeability. Even minor components of the deposits can sometimes cause severe corrosion of the underlying metal such as the hot corrosion caused by vanadium in the deposits of fired boilers [5].

Fouling is responsible for the emission of many millions of tonnes of carbon dioxide as well as the use and disposal of hazardous cleaning chemicals. Data from oil refineries suggest that crude oil fouling accounts for about 10% of the total CO2 emission of these plants. Wastes generated from the cleaning of heat exchangers may contain hazardous wastes such as lead and chromium, although some refineries which do not produce leaded gasoline and which use non-chrome corrosion inhibitors typically do not generate sludge that contains these constituents. Oily wastewater is also generated during heat exchanger cleaning.

The factors that govern fouling in heat exchangers are many and varied. Of such factors some may be related to the feed properties such as its chemical nature, density, viscosity, diffusivity, pour and cloud points, interfacial properties and colloidal stability factors. The chemical nature of the feed in particular can be an important factor affecting to a large degree the rate and extent of fouling. This includes the chemical composition of the feed and the stability of its components and their compatibility with one another and with heat exchanger surfaces as well as the presence in the feed of unsaturated and unstable compounds, inorganic salts and trace elements such as sulphur, nitrogen and oxygen. The feed storage conditions and its exposure to oxygen on storage in particular can in most cases also affect materially the rate and nature of fouling.

Fouling in Heat Exchangers 61

Detailed analysis of deposits from the heat exchanger may provide an excellent clue to fouling mechanisms. It can be used to identify and provide valuable information about such mechanisms. The deposits consist primarily of organic material that is predominantly asphaltenic in nature, with some inorganic deposits, mainly iron salts such as iron sulphide. The inorganic content of the deposits is relatively consistent in most cases at 22-26% [6].

Deposit analysis is performed by taking a sample and extracting any degraded hydrocarbon oil by using a solvent, such as methyl chloride, that is effective at removing hydrocarbon

The remaining material from this extraction will consist of any organic polymers, coke, and inorganic components. The basic analysis of the non-extractable material involves ashing in which organic and volatile inorganic compounds are lost. By this means, volatile inorganics such as chlorides and sulphur compounds which are lost on ashing, may be determined. The detection of iron sulphide or other volatile inorganic materials determines the cause of inorganic fouling. These values can be compared throughout the exchanger train [6]. The non-volatile material or ash will include all oxidised metallic salt–type materials or corrosion products. The presence of iron in the ash may indicate corrosion in tankage in an upstream unit or in the exchanger train itself. This basic analysis indicates if the deposits are

Special techniques and tools such as the use of optical microscopy and solubility in solvents may be used for the analysis of the non-extractable material. Infrared analysis can identify various functional groups present in the deposit which may include nitrogen, carbonyls, and unsaturated paraffinic or aromatic compounds which are polymerisation precursors, identified in feed stream characterisation [6]. The carbon and hydrogen content of the nonextractable deposit can be determined by elemental analysis. If the carbon to hydrogen ratio is very high, it may indicate that the majority of the organic portion of the deposit is coke. The coke may have been particles entrained in the stream or material which has been thermally dehydrogenated in the heat exchangers. The carbon to hydrogen ratio also indicates whether the deposit is more paraffinic or aromatic. This information helps identify

In Table 2 analytical results are shown from deposits obtained from the four chain feed/effluent heat exchangers in which the hot product effluent is used for pre-heating the cold naphtha feedstock for a naphtha hydrotreater plant at the Homs Oil Refinery [7]. This plant is one of the most important units at the Homs Refinery, with an annual capacity of 480,000 tons/yr. It is used to remove impurities such as sulphur, nitrogen, oxygen, halides and trace metal impurities that may deactivate reforming catalysts. Furthermore, the quality of the naphtha fractions is also upgraded by reducing potential gum formation as a result of the conversion of olefins and diolefins into paraffins. The process utilises a catalyst (Hydrobon) in the presence of substantial amounts of hydrogen under high pressures (50 bars) and temperatures (320°C) (Fig. 2). A major fouling problem was encountered early on in the heat exchangers, indicated by an increased pressure drop, decreased flow rate and

oils and low molecular weight polymers that may have been trapped in the deposit.

primarily organic or inorganic.

the polymers formed [6].

lower temperatures at the heat exchangers outlet.

Other factors of equal importance to the feed properties may be related to operating conditions and equipment design, such as feed temperature, bulk fluid velocity or flow rate, heat exchanger geometry, nature of alloy used and wettability of surfaces where fouling occurs. The rate of fouling is feed temperature dependent with different rates of fouling between the feed inlet and outlet sides of the heat exchanger. In a shell and tube heat exchanger, the conventional segment baffle geometry is largely responsible for higher fouling rates. Uneven velocity profiles, back-flows and eddies generated on the shell side of a segmentally-baffled heat exchanger results in higher fouling and shorter run lengths between periodic cleaning and maintenance of tube bundles.

All these and other factors that may affect fouling need to be considered and taken into account in order to be able to prevent fouling if possible or to predict the rate of fouling or fouling factor prior to taking the necessary steps for fouling mitigation, control and removal.

## **2. Fouling mechanisms and stages**

Fouling can be divided into a number of distinctively different mechanisms. Generally speaking, several of these fouling mechanisms occur at the same time and each requires a different prevention technique. Of these different mechanisms some represent different stages in the process of fouling. The chief fouling mechanisms or stages include:


Detailed analysis of deposits from the heat exchanger may provide an excellent clue to fouling mechanisms. It can be used to identify and provide valuable information about such mechanisms. The deposits consist primarily of organic material that is predominantly asphaltenic in nature, with some inorganic deposits, mainly iron salts such as iron sulphide. The inorganic content of the deposits is relatively consistent in most cases at 22-26% [6].

60 MATLAB – A Fundamental Tool for Scientific Computing and Engineering Applications – Volume 3

also affect materially the rate and nature of fouling.

between periodic cleaning and maintenance of tube bundles.

**2. Fouling mechanisms and stages** 

negative total fouling amount.

erosion and removal.

attachment leading to deposit formation.

removal.

The factors that govern fouling in heat exchangers are many and varied. Of such factors some may be related to the feed properties such as its chemical nature, density, viscosity, diffusivity, pour and cloud points, interfacial properties and colloidal stability factors. The chemical nature of the feed in particular can be an important factor affecting to a large degree the rate and extent of fouling. This includes the chemical composition of the feed and the stability of its components and their compatibility with one another and with heat exchanger surfaces as well as the presence in the feed of unsaturated and unstable compounds, inorganic salts and trace elements such as sulphur, nitrogen and oxygen. The feed storage conditions and its exposure to oxygen on storage in particular can in most cases

Other factors of equal importance to the feed properties may be related to operating conditions and equipment design, such as feed temperature, bulk fluid velocity or flow rate, heat exchanger geometry, nature of alloy used and wettability of surfaces where fouling occurs. The rate of fouling is feed temperature dependent with different rates of fouling between the feed inlet and outlet sides of the heat exchanger. In a shell and tube heat exchanger, the conventional segment baffle geometry is largely responsible for higher fouling rates. Uneven velocity profiles, back-flows and eddies generated on the shell side of a segmentally-baffled heat exchanger results in higher fouling and shorter run lengths

All these and other factors that may affect fouling need to be considered and taken into account in order to be able to prevent fouling if possible or to predict the rate of fouling or fouling factor prior to taking the necessary steps for fouling mitigation, control and

Fouling can be divided into a number of distinctively different mechanisms. Generally speaking, several of these fouling mechanisms occur at the same time and each requires a different prevention technique. Of these different mechanisms some represent different

1. Initiation or delay period. This is the clean surface period before dirt accumulation. The accumulation of relatively small amounts of deposit can even lead to improved heat transfer, relative to clean surface, and give an appearance of "negative" fouling rate and

4. Phase separation and deposition involving nucleation or initiation of fouling sites and

5. Growth, aging and hardening and the increase of deposits strength or auto-retardation,

stages in the process of fouling. The chief fouling mechanisms or stages include:

2. Particulate fouling and particle formation, aggregation and flocculation.

3. Mass transport and migration of foulants to the fouling sites.

Deposit analysis is performed by taking a sample and extracting any degraded hydrocarbon oil by using a solvent, such as methyl chloride, that is effective at removing hydrocarbon oils and low molecular weight polymers that may have been trapped in the deposit.

The remaining material from this extraction will consist of any organic polymers, coke, and inorganic components. The basic analysis of the non-extractable material involves ashing in which organic and volatile inorganic compounds are lost. By this means, volatile inorganics such as chlorides and sulphur compounds which are lost on ashing, may be determined. The detection of iron sulphide or other volatile inorganic materials determines the cause of inorganic fouling. These values can be compared throughout the exchanger train [6]. The non-volatile material or ash will include all oxidised metallic salt–type materials or corrosion products. The presence of iron in the ash may indicate corrosion in tankage in an upstream unit or in the exchanger train itself. This basic analysis indicates if the deposits are primarily organic or inorganic.

Special techniques and tools such as the use of optical microscopy and solubility in solvents may be used for the analysis of the non-extractable material. Infrared analysis can identify various functional groups present in the deposit which may include nitrogen, carbonyls, and unsaturated paraffinic or aromatic compounds which are polymerisation precursors, identified in feed stream characterisation [6]. The carbon and hydrogen content of the nonextractable deposit can be determined by elemental analysis. If the carbon to hydrogen ratio is very high, it may indicate that the majority of the organic portion of the deposit is coke. The coke may have been particles entrained in the stream or material which has been thermally dehydrogenated in the heat exchangers. The carbon to hydrogen ratio also indicates whether the deposit is more paraffinic or aromatic. This information helps identify the polymers formed [6].

In Table 2 analytical results are shown from deposits obtained from the four chain feed/effluent heat exchangers in which the hot product effluent is used for pre-heating the cold naphtha feedstock for a naphtha hydrotreater plant at the Homs Oil Refinery [7]. This plant is one of the most important units at the Homs Refinery, with an annual capacity of 480,000 tons/yr. It is used to remove impurities such as sulphur, nitrogen, oxygen, halides and trace metal impurities that may deactivate reforming catalysts. Furthermore, the quality of the naphtha fractions is also upgraded by reducing potential gum formation as a result of the conversion of olefins and diolefins into paraffins. The process utilises a catalyst (Hydrobon) in the presence of substantial amounts of hydrogen under high pressures (50 bars) and temperatures (320°C) (Fig. 2). A major fouling problem was encountered early on in the heat exchangers, indicated by an increased pressure drop, decreased flow rate and lower temperatures at the heat exchangers outlet.

#### Particulate fouling

Particulate fouling, which is the most common form of fouling, can be defined as the process in which particles in the process stream deposit onto heat exchanger surfaces. These particles include particles originally carried by the feed stream before entering the heat exchanger and particles formed in the heat exchanger itself as a result of various reactions, aggregation and flocculation. Particulate fouling increases with particle concentration, and typically particles greater than 1 ppm lead to significant fouling problems.

Fouling in Heat Exchangers 63

exchanger surfaces. These solid particles are for the most part insoluble inorganic particles such as corrosion products (iron sulphide and rust), catalyst particles or fines, dirt, silt and sand particles, and other inorganic salts such as sodium chloride, calcium chloride and magnesium chloride. The feed streams may also contain some organic particles that may

Many streams including cooling water and other product streams from different units or plants may contain solid particles. In particular, streams from such oil refinery units as vacuum units, visbreakers, and cokers may have more particulates and metals than straightrun products due to the heavier nature of the feeds processed. Streams can also be purchased from other refiners. Due to the increased transit time and exposure to oxygen before being fed to the unit these feeds may have higher particulate levels as a result of

Particles in the fluid stream, regardless of whether they are organic or inorganic in nature,

Typically, particles in the fluid stream greater than 1 ptb (pounds per thousand barrels) lead to significant fouling problems in the unit. Their effect on fouling can be avoided however if these particles are removed by solid-liquid filtration, sedimentation, centrifugation or by any of various fluid cleaning devices. The only particles that need to be considered in this regard are those that are not filterable or those particles that are left to proceed to the heat

The amount of filterable solids in the stream, reported in ptb or wt% (weight percent), may be determined by filtration of the unit feed. Filterable solids analysis can evaluate a stream deposition potential by indicating the type of materials that could contribute to fouling if

Table 3 shows the analysis of filterable solids in the naphtha feed stream to the heat exchangers of the hydrotreater unit at the Homs oil refinery. The feedstock for this unit is a blend of light and heavy straight-run naphtha fractions from four different topping units. The resulting blend is left in a blending tank for a sufficient period of time to allow for equilibrium conditions to be established [8]. To evaluate the quantity of particulate solids which are entrained with the naphtha stream before entering the heat exchangers, a number of samples of the naphtha feed were filtered and the amount of entrained particles determined. Two samples of the filterable solids were taken, one sample was taken from the feed entering a macrofilter on the unit boundary and the other from a second macrofilter on the feed pump suction. The nature of the materials entrained was then determined by ashing and analysing these two samples (Table 3). The size distribution of the filterable solid

Examination of the deposit analysis for heat exchanger D (Table 2), where the deposits are a mixture of inorganic (42%) and organic (58%) deposits, indicate particulate and polymerisation fouling. The nature of particulate fouling in D is confirmed by the variation

have been formed during their storage or transport.

polymerisation reactions and corrosion [6].

allowed to pass through to the heat exchanger.

particles was also determined (Table 4).

exchanger.

fall in general into tow classes: basic sediment and filterable solids.


**Table 2.** Analysis of deposits on heat exchanger surfaces [7].

**Figure 2.** Naphtha hydrotreating unit

### **2.1. Particles in the feed stream**

Particles in the fluid feed stream are solid particles which are entrained or contained in the feed stream before entering the heat exchanger and which can settle out upon the heat exchanger surfaces. These solid particles are for the most part insoluble inorganic particles such as corrosion products (iron sulphide and rust), catalyst particles or fines, dirt, silt and sand particles, and other inorganic salts such as sodium chloride, calcium chloride and magnesium chloride. The feed streams may also contain some organic particles that may have been formed during their storage or transport.

62 MATLAB – A Fundamental Tool for Scientific Computing and Engineering Applications – Volume 3

typically particles greater than 1 ppm lead to significant fouling problems.

**Table 2.** Analysis of deposits on heat exchanger surfaces [7].

**Figure 2.** Naphtha hydrotreating unit

**2.1. Particles in the feed stream** 

Particulate fouling, which is the most common form of fouling, can be defined as the process in which particles in the process stream deposit onto heat exchanger surfaces. These particles include particles originally carried by the feed stream before entering the heat exchanger and particles formed in the heat exchanger itself as a result of various reactions, aggregation and flocculation. Particulate fouling increases with particle concentration, and

Heat exchanger A B1 B2 C D Loss at 105°C (wt %) 1.17 1.03 1.05 1.15 1.14 Loss at 550°C (wt %) 79.70 95.10 90.17 94.42 57.17 Loss at 840°C (wt %) 80.00 95.29 90.19 94.48 57.99 Ash (wt %) 20.00 4.71 9.81 5.52 42.01 Chloride (wt %) 170 435 0 664 508 Sulphur (wt %) 17.00 13.50 13.80 10.20 13.00 Ammonium (ppm) 42 1184 43 134 4969 Iron (wt % of ash) 19.30 2.83 1.70 2.80 15.58 Sodium (ppm of ash) 1473 1047 825 3301 914 Calcium (ppm of ash) 459 179 78 377 1431 Magnesium(ppm of ash) 90 41 19 102 1341 Chromium (ppm of ash) 231 107 1166 196 1096 Copper (ppm of ash) 511 319 74 443 126 Nickel (ppm of ash) 378 129 63 90 52

Particles in the fluid feed stream are solid particles which are entrained or contained in the feed stream before entering the heat exchanger and which can settle out upon the heat

Particulate fouling

Many streams including cooling water and other product streams from different units or plants may contain solid particles. In particular, streams from such oil refinery units as vacuum units, visbreakers, and cokers may have more particulates and metals than straightrun products due to the heavier nature of the feeds processed. Streams can also be purchased from other refiners. Due to the increased transit time and exposure to oxygen before being fed to the unit these feeds may have higher particulate levels as a result of polymerisation reactions and corrosion [6].

Particles in the fluid stream, regardless of whether they are organic or inorganic in nature, fall in general into tow classes: basic sediment and filterable solids.

Typically, particles in the fluid stream greater than 1 ptb (pounds per thousand barrels) lead to significant fouling problems in the unit. Their effect on fouling can be avoided however if these particles are removed by solid-liquid filtration, sedimentation, centrifugation or by any of various fluid cleaning devices. The only particles that need to be considered in this regard are those that are not filterable or those particles that are left to proceed to the heat exchanger.

The amount of filterable solids in the stream, reported in ptb or wt% (weight percent), may be determined by filtration of the unit feed. Filterable solids analysis can evaluate a stream deposition potential by indicating the type of materials that could contribute to fouling if allowed to pass through to the heat exchanger.

Table 3 shows the analysis of filterable solids in the naphtha feed stream to the heat exchangers of the hydrotreater unit at the Homs oil refinery. The feedstock for this unit is a blend of light and heavy straight-run naphtha fractions from four different topping units. The resulting blend is left in a blending tank for a sufficient period of time to allow for equilibrium conditions to be established [8]. To evaluate the quantity of particulate solids which are entrained with the naphtha stream before entering the heat exchangers, a number of samples of the naphtha feed were filtered and the amount of entrained particles determined. Two samples of the filterable solids were taken, one sample was taken from the feed entering a macrofilter on the unit boundary and the other from a second macrofilter on the feed pump suction. The nature of the materials entrained was then determined by ashing and analysing these two samples (Table 3). The size distribution of the filterable solid particles was also determined (Table 4).

Examination of the deposit analysis for heat exchanger D (Table 2), where the deposits are a mixture of inorganic (42%) and organic (58%) deposits, indicate particulate and polymerisation fouling. The nature of particulate fouling in D is confirmed by the variation of fouling factor with time, with no induction time or delay period indicated (Fig. 3). The fouling factor curve is linear with saw-tooth shape, where both the fouling factor and the deposition rate increase with time. This means continuous build up of the fouling layer followed by break off periods [9].

Fouling in Heat Exchangers 65

significant effect on the fouling encountered in certain chemical processes. Such contaminants may include oxygen, nitrogen, NH3, H2S, CN, HCN, Hg, unsaturates, organic sulphides and chlorides, and heavy hydrocarbon compounds such as paraffin wax, resins, asphaltenes, and organometallic compounds. Individual metals, which may exist as metal salts in the feed stream, can catalyse different polymerisation reactions. The concentrations of such metals are typically very low, not exceeding few ppms. However, small concentrations of certain metals can have a significant effect on catalysing different foulingrelated polymerisation reactions. Metal detectors on unit feed samples can detect individual

**Corrosion fouling** is fouling deposit formation as a result of the corrosion of the substrate metal of heat transfer surfaces. This type of corrosion should not be confused, however, with the under-deposit corrosion, referred to earlier, which is one of the aftereffects of fouling.

Corrosion fouling is a mechanism which is dependent on several factors such as thermal resistance, surface roughness and composition of the substrate and fluid stream. In particular, impurities present in the fluid stream can greatly contribute to the onset of corrosion. Such impurities include hydrogen sulphide, ammonia and hydrogen chloride. In crude oil, for example, sulphur and nitrogen compounds are two very common contaminants which are mostly decomposed in certain situations to hydrogen sulphide and ammonia respectively. Chlorides which may be found in oil streams are converted to

R-Cl + H2 → HCl + R

The chlorides may enter the refinery as salt with the crude. Chlorides in the oil stream may also be derived from various chemicals used in the oil industry which can contain high levels of chloride. Such chemicals include tertiary oil recovery enhancement chemicals and solvents used to clean tankers, barges, trucks and pipelines. As the crude oil is processed,

metals in the stream at less than 1 ppm.

**Figure 3.** Variation of fouling factor in exchanger D in 2001.

hydrogen chloride by the following reaction.

## **2.2. Particle formation**

Chemical particle formation is the basic mechanism of particle formation in heat exchangers fluid streams, although organic material growth and biological particle formation, or biofouling, may occur in sea water systems and in types of waste treatment systems. Biofouling may be of two kinds: microbial fouling, due to microorganisms (bacteria, algae, and fungi) and their products, and macrobial fouling, due to the growth of macroorganisms such as barnacles, sponges, seaweeds or mussels. On contact with heat-transfer surfaces, these organisms can attach and breed, reducing thereby both flow and heat transfer to an absolute minimum and sometimes completely clogging the fluid passages. Such organisms may also entrap silt or other suspended solids and give rise to deposit corrosion. Corrosion due to biological attachment to heat transfer surfaces is known as microbiologically influenced corrosion. For open recirculating systems, bacteria concentrations of the order of 1 x 105 cells/ml and fungi of 1 x 103 cells/ml may be regarded as limiting values [10].




**Table 4.** Size distribution of the filterable solid particles

**Chemical particle formation** can be the result of either corrosion or decomposition and polymerisation reactions. Trace contaminants present in the fluid stream can have a significant effect on the fouling encountered in certain chemical processes. Such contaminants may include oxygen, nitrogen, NH3, H2S, CN, HCN, Hg, unsaturates, organic sulphides and chlorides, and heavy hydrocarbon compounds such as paraffin wax, resins, asphaltenes, and organometallic compounds. Individual metals, which may exist as metal salts in the feed stream, can catalyse different polymerisation reactions. The concentrations of such metals are typically very low, not exceeding few ppms. However, small concentrations of certain metals can have a significant effect on catalysing different foulingrelated polymerisation reactions. Metal detectors on unit feed samples can detect individual metals in the stream at less than 1 ppm.

**Figure 3.** Variation of fouling factor in exchanger D in 2001.

64 MATLAB – A Fundamental Tool for Scientific Computing and Engineering Applications – Volume 3

followed by break off periods [9].

**2.2. Particle formation** 

of fouling factor with time, with no induction time or delay period indicated (Fig. 3). The fouling factor curve is linear with saw-tooth shape, where both the fouling factor and the deposition rate increase with time. This means continuous build up of the fouling layer

Chemical particle formation is the basic mechanism of particle formation in heat exchangers fluid streams, although organic material growth and biological particle formation, or biofouling, may occur in sea water systems and in types of waste treatment systems. Biofouling may be of two kinds: microbial fouling, due to microorganisms (bacteria, algae, and fungi) and their products, and macrobial fouling, due to the growth of macroorganisms such as barnacles, sponges, seaweeds or mussels. On contact with heat-transfer surfaces, these organisms can attach and breed, reducing thereby both flow and heat transfer to an absolute minimum and sometimes completely clogging the fluid passages. Such organisms may also entrap silt or other suspended solids and give rise to deposit corrosion. Corrosion due to biological attachment to heat transfer surfaces is known as microbiologically influenced corrosion. For open recirculating systems, bacteria concentrations of the order of

1 x 105 cells/ml and fungi of 1 x 103 cells/ml may be regarded as limiting values [10].

Loss at 105°C (wt. %) 10.0 0.1 Loss at 550°C (wt. %) 28.3 25.3 Loss at 840°C (wt. %) 30.4 26.7 Ash (wt. %) 69.5 73.2 Carbon (wt. %) 2.6 6.4 Sulphur (wt. %) 36.9 19.7 Sulphates (wt. %) 55.8 50.7 Chloride (ppm) - 281 Ammonium (ppm) - 52 Iron (wt. % of ash) 45 58 Sodium (ppm of ash) - 9 Calcium (ppm of ash) - 161

Mesh size (μm) < 90 90 125 355 Particle distribution (%) 24 8 36 32

**Chemical particle formation** can be the result of either corrosion or decomposition and polymerisation reactions. Trace contaminants present in the fluid stream can have a

**Table 3.** Analysis of two samples of the filterable solids.

**Table 4.** Size distribution of the filterable solid particles

Feed filter Pump filter

**Corrosion fouling** is fouling deposit formation as a result of the corrosion of the substrate metal of heat transfer surfaces. This type of corrosion should not be confused, however, with the under-deposit corrosion, referred to earlier, which is one of the aftereffects of fouling.

Corrosion fouling is a mechanism which is dependent on several factors such as thermal resistance, surface roughness and composition of the substrate and fluid stream. In particular, impurities present in the fluid stream can greatly contribute to the onset of corrosion. Such impurities include hydrogen sulphide, ammonia and hydrogen chloride. In crude oil, for example, sulphur and nitrogen compounds are two very common contaminants which are mostly decomposed in certain situations to hydrogen sulphide and ammonia respectively. Chlorides which may be found in oil streams are converted to hydrogen chloride by the following reaction.

#### R-Cl + H2 → HCl + R

The chlorides may enter the refinery as salt with the crude. Chlorides in the oil stream may also be derived from various chemicals used in the oil industry which can contain high levels of chloride. Such chemicals include tertiary oil recovery enhancement chemicals and solvents used to clean tankers, barges, trucks and pipelines. As the crude oil is processed,

some of these chemicals and solvents, which are thermally stable and not soluble in water, pass overhead in the main tower of the atmospheric distillation unit along with the naphtha.

Fouling in Heat Exchangers 67

leaving the base metal exposed. The reaction rate is then only limited by the chloride ion

Fe++ + 2Cl → FeCl2

+ 2H+ → H2

As the chloride concentration in water is reduced by removing the source, diluting with additional water or neutralising with a base, the pH will increase. Hydrogen sulphide will begin to react with the exposed iron and start building a new protective layer. This sulphide lattice gets stronger as the pH increases to 6 and above. The corrosion rate falls off to a

Hydrogen chloride will also cause serious fouling problems if ammonia is present in the system. The ammonia reacts with hydrogen chloride to form ammonium chloride which may cause fouling and plugging problems. In the cooler parts of the unit, the ammonium chloride will condense from the vapour phase and solidify and deposit directly and accumulate on the walls. The salt can also break away from the walls and be carried downstream to eventually deposit somewhere else. If free water is present, ammonium chloride will be absorbed directly from the vapour phase into the water and no solid salts will form on the equipment. Another problem associated with ammonium chloride salt deposits is under deposit pitting corrosion as the hygroscopic nature of the salt will result in a wet environment at the wall under the deposit. The chloride ions will react with the iron to form iron chloride causing serious localised corrosion, the reaction rate accelerating in the presence of hydrogen sulphide. The sulphide ion as part of an ammonium sulphide salt will react with the iron chloride to form iron sulphide, thus releasing the chloride ion to start over.

Any excess ammonia available may react with the disulphide ions present in the solution to form ammonium sulphide salts, but only after most of the chloride has been neutralised. While hydrogen sulphide is only slightly soluble in the water, its salt is highly soluble. Therefore, as the pH is raised to 6 and higher the free ammonia present reacts with the small quantity of hydrogen sulphide in solution, making more ammonium sulphide salts. The rich hydrogen sulphide vapour above the water will continuously replace the consumed hydrogen sulphide. The overall sulphide concentration in the water increases making it

The sharp increase in corrosion rate in the 6.8 to 7.3 pH range is related to the concentration of the ammonium chloride and sulphide salts present. In large quantities these salts can become aggressive, especially the sulphide salts. If the pH is raised further the corrosion rate again falls off to a very low value. This is because the sulphide lattice has formed into a very

The iron content in the deposits obtained from the heat exchangers may be an indication of fouling by corrosion. Although polymerisation may account for about 80% of the total

concentration in the solution at low pH. Loss of wall metal takes place very rapidly.

In water the chloride ions react directly with any exposed iron to form FeCl2.

2e-

minimum.

difficult to raise the pH much further.

strong hard film that cannot easily be broken.

In the hydrotreater feed stream, chloride levels as high as 50 wt. ppm have been reported. High levels of chloride were detected with the filterable solids in the naphtha feed stream to the heat exchangers of the hydrotreater unit at the Homs refinery (Table 3) and in the deposits obtained from the heat exchangers (Table 2). Furthermore, the makeup hydrogen from the platforming unit will always contain trace quantities of hydrogen chloride. In order to maintain catalyst performance, modern platforming catalysts require a small, but continuous dosage of chloride, some of which is always stripped and leaves the platforming unit in the net gas stream that supplies the hydrotreater with makeup hydrogen.

In a hydrogen sulphide environment the sulphur reacts with the exposed iron to form iron sulphide compounds. This happens in both the hot and cooler sections of the unit. The sulphur effectively corrodes the plant. However, once reacted, the iron sulphide forms a complex protective scale or lattice on the base metal, which inhibits further corrosion. The sulphide lattice would remain in equilibrium with its surroundings and the corrosion rate would be minimal if no other impurities were present in the system. The presence of other impurities, however, can accelerate corrosion as these impurities interact with the sulphide lattice.

Of the impurities that contribute to corrosion and fouling, hydrogen chloride may be the most important. By itself hydrogen chloride does not cause a problem. It will not foul equipment or corrode the carbon steel in the unit. Chloride corrosion and fouling, however, take place when hydrogen chloride, ammonia, and water all interact in the colder sections of the unit to defeat the protective sulphide lattice. The extent of the damage depends on their concentration and is directly dependent on pH, with the corrosion rate increasing rapidly with pH decrease.

Hydrogen chloride will become corrosive when it comes in contact with free water, i.e. water that is not in the vapour phase or is not saturated in the liquid hydrocarbon. Oil products are almost always saturated with water, and entrained water, even if it is less of a problem, does occur in most cases. Furthermore, continuous water wash at key locations is recommended as part of the solution to minimise the effects of chloride corrosion and fouling and this further contributes to the total water in the system.

Hydrogen chloride is highly soluble in water, and in a free water environment, any hydrogen chloride present in the vapour or hydrocarbon liquid will be quickly absorbed by the water, thus driving the pH down to approximately 1.

If the iron sulphide lattice is intact this chloride competes with the bisulphate ion (SH- ) for the iron ions in the lattice:

$$\text{S-Fe-S-Fe-SH} + \text{Cl}^\cdot \rightleftharpoons \text{Fe-S-S-Fe-Cl} + \text{SH}^\cdot$$

With a high concentration of hydrogen chloride present the reaction shifts to the right. As more and more bisulphate is released from the sulphide lattice, it eventually dissolves leaving the base metal exposed. The reaction rate is then only limited by the chloride ion concentration in the solution at low pH. Loss of wall metal takes place very rapidly.

In water the chloride ions react directly with any exposed iron to form FeCl2.

66 MATLAB – A Fundamental Tool for Scientific Computing and Engineering Applications – Volume 3

unit in the net gas stream that supplies the hydrotreater with makeup hydrogen.

however, can accelerate corrosion as these impurities interact with the sulphide lattice.

with pH decrease.

the iron ions in the lattice:

some of these chemicals and solvents, which are thermally stable and not soluble in water, pass overhead in the main tower of the atmospheric distillation unit along with the naphtha.

In the hydrotreater feed stream, chloride levels as high as 50 wt. ppm have been reported. High levels of chloride were detected with the filterable solids in the naphtha feed stream to the heat exchangers of the hydrotreater unit at the Homs refinery (Table 3) and in the deposits obtained from the heat exchangers (Table 2). Furthermore, the makeup hydrogen from the platforming unit will always contain trace quantities of hydrogen chloride. In order to maintain catalyst performance, modern platforming catalysts require a small, but continuous dosage of chloride, some of which is always stripped and leaves the platforming

In a hydrogen sulphide environment the sulphur reacts with the exposed iron to form iron sulphide compounds. This happens in both the hot and cooler sections of the unit. The sulphur effectively corrodes the plant. However, once reacted, the iron sulphide forms a complex protective scale or lattice on the base metal, which inhibits further corrosion. The sulphide lattice would remain in equilibrium with its surroundings and the corrosion rate would be minimal if no other impurities were present in the system. The presence of other impurities,

Of the impurities that contribute to corrosion and fouling, hydrogen chloride may be the most important. By itself hydrogen chloride does not cause a problem. It will not foul equipment or corrode the carbon steel in the unit. Chloride corrosion and fouling, however, take place when hydrogen chloride, ammonia, and water all interact in the colder sections of the unit to defeat the protective sulphide lattice. The extent of the damage depends on their concentration and is directly dependent on pH, with the corrosion rate increasing rapidly

Hydrogen chloride will become corrosive when it comes in contact with free water, i.e. water that is not in the vapour phase or is not saturated in the liquid hydrocarbon. Oil products are almost always saturated with water, and entrained water, even if it is less of a problem, does occur in most cases. Furthermore, continuous water wash at key locations is recommended as part of the solution to minimise the effects of chloride corrosion and

Hydrogen chloride is highly soluble in water, and in a free water environment, any hydrogen chloride present in the vapour or hydrocarbon liquid will be quickly absorbed by

) for

If the iron sulphide lattice is intact this chloride competes with the bisulphate ion (SH-

S-Fe-S-Fe-SH + Cl- Fe-S-S-Fe-Cl + SH-

With a high concentration of hydrogen chloride present the reaction shifts to the right. As more and more bisulphate is released from the sulphide lattice, it eventually dissolves

fouling and this further contributes to the total water in the system.

the water, thus driving the pH down to approximately 1.

$$\begin{array}{c} \mathrm{Fe^{+} + 2Cl \longrightarrow FeCl\_{2}} \\\\ 2e^{-} + 2H^{+} \longrightarrow H\_{2} \end{array}$$

As the chloride concentration in water is reduced by removing the source, diluting with additional water or neutralising with a base, the pH will increase. Hydrogen sulphide will begin to react with the exposed iron and start building a new protective layer. This sulphide lattice gets stronger as the pH increases to 6 and above. The corrosion rate falls off to a minimum.

Hydrogen chloride will also cause serious fouling problems if ammonia is present in the system. The ammonia reacts with hydrogen chloride to form ammonium chloride which may cause fouling and plugging problems. In the cooler parts of the unit, the ammonium chloride will condense from the vapour phase and solidify and deposit directly and accumulate on the walls. The salt can also break away from the walls and be carried downstream to eventually deposit somewhere else. If free water is present, ammonium chloride will be absorbed directly from the vapour phase into the water and no solid salts will form on the equipment. Another problem associated with ammonium chloride salt deposits is under deposit pitting corrosion as the hygroscopic nature of the salt will result in a wet environment at the wall under the deposit. The chloride ions will react with the iron to form iron chloride causing serious localised corrosion, the reaction rate accelerating in the presence of hydrogen sulphide. The sulphide ion as part of an ammonium sulphide salt will react with the iron chloride to form iron sulphide, thus releasing the chloride ion to start over.

Any excess ammonia available may react with the disulphide ions present in the solution to form ammonium sulphide salts, but only after most of the chloride has been neutralised. While hydrogen sulphide is only slightly soluble in the water, its salt is highly soluble. Therefore, as the pH is raised to 6 and higher the free ammonia present reacts with the small quantity of hydrogen sulphide in solution, making more ammonium sulphide salts. The rich hydrogen sulphide vapour above the water will continuously replace the consumed hydrogen sulphide. The overall sulphide concentration in the water increases making it difficult to raise the pH much further.

The sharp increase in corrosion rate in the 6.8 to 7.3 pH range is related to the concentration of the ammonium chloride and sulphide salts present. In large quantities these salts can become aggressive, especially the sulphide salts. If the pH is raised further the corrosion rate again falls off to a very low value. This is because the sulphide lattice has formed into a very strong hard film that cannot easily be broken.

The iron content in the deposits obtained from the heat exchangers may be an indication of fouling by corrosion. Although polymerisation may account for about 80% of the total

fouling associated with the "A" heat exchanger, in the Homs hydrotreating plant, fouling by corrosion is not negligible, with about 19% of the total fouling may be due to corrosion, as is clearly indicated by the iron content of the deposits obtained from this exchanger (Table 2).

Fouling in Heat Exchangers 69

Some of these catalytically reactive metals are iron, copper, nickel, vanadium, chromium,

In non-free radical polymerisation, polymer formation results from the reaction of two different molecules under the right conditions. One of the reactive molecules may be a radical, or a compound from a free radical-initiated polymerisation step. Basic compounds can react with other compounds or with themselves to form polymers by several different polymerisation mechanisms. In condensation polymerisation, two large radicals or compounds react together to form an even larger compound, but in their reaction also generate a smaller compound, such as water. This new larger compound can continue to react with other reactive species in the feed stream to make higher molecular weight polymers. At some point, the polymer will either get so large in size that it is no longer able to stay entrained or soluble in the fluid stream and deposit, or all the different compounds

Various laboratory tests can provide an indication of a stream's polymerisation potential. These include laboratory simulations and analytical characterisation, to identify specific compounds in the feed which are known to contribute to polymerisation mechanisms. Such polymerisation precursors may include unsaturated hydrocarbons, acidic compounds,

The presence of unsaturated components in the feed stream contributes significantly to polymerisation reactions, particularly at high temperatures. The bromine number is a method of measuring the degree of unsaturation in a feed stream. The unsaturated bonds react with bromine, and the amount of bromine reacted is an indication of the degree of

The neutralisation, or acid, number measures the acidity of the fluid as it is titrated with a base. This number can be an indication of fouling tendency, where the more acidic the feed stream, the greater is its tendency to foul. This is most likely due to the fact that acidic

The basic nitrogen test determines the amount of basic compounds in a sample, assumed to be mostly amines, by titrating with a mixture of organic acids. This method can, however,

A method of determining a sample's oxidative polymerisation potential is to run a potential gums test. This test is a method of determining a sample's oxidative polymerisation potential. In this test, the fluid is subjected to 100% oxygen for four hours, at 100ºC, in a pressurised sample bomb. The measured gum content, as compared to an initial gum value,

Detailed deposit analysis, as mentioned earlier, can also indicate the occurrence of polymerisation. It is apparent from examination of the deposit analysis results shown in Table 2 that most deposits are organic in nature, as the loss reported on heating the deposit

that can react with it are consumed, and no further polymer is formed [6].

compounds, as mentioned above, can promote free radical polymerisation.

will indicate the impact of oxygen on the stream's polymerisation potential.

amines, carbonyls, mercaptans and pyrrole nitrogen.

overestimate the basic nitrogen content.

calcium, and magnesium [6]

unsaturation [12].

**Coking and Polymerisation** are major causes of fouling in heat exchangers. Decomposition of organic products can lead to the formation of very viscous tar or solid coke particles at high temperatures and polymerisation involves the formation of undesirable organic sediments or polymers. The coke particles and polymers formed in the heat exchanger may grow to such a large size that they drop out of solution and deposit on the process equipment. Such deposits can be extremely tenacious and may require burning off the deposit to return the heat exchanger to satisfactory operation.

There are two major polymerisation mechanisms which can occur in the feed stream: free radical and non-free radical polymerisation.

Free radical polymerisation occurs when a free radical is formed and continues to react with other molecules. The free radicals continue to propagate in the feed stream producing longer chain polymers which will continue to be produced as long as free radicals are being formed. Free radical polymerisation is easily initiated in the presence of light and heat and its rate for polymer formation increases exponentially with temperature. A general rule is that for every 10ºC increase in temperature the rate of polymer formation doubles. Free radical polymerisation readily takes place in heat exchanger tubes and storage tanks [6].

The formation of free radicals has been investigated extensively and it is known that numerous types of free radicals can be formed in a feed stream. These include alkyl radicals produced by the breaking of double or unsaturated bonds as well as other types of precursors such as nitrogen and sulphur radicals which arc easily formed at the temperatures found in the heat exchanger train. Organic sulphur, nitrogen and oxygen compounds increase the potential for various polymerisation reactions, depending on the form in which they exist. Acidic compounds can promote free radical polymerisation by initiating free radicals through the formation of a positive ion or cation. Additional polymerisation precursors include carbonyls, mercaptans, and pyrrole nitrogen.

Oxygen may also react with hydrocarbons to form peroxide free radicals, a step that could occur in the storage tank. When the temperature is increased in the heat exchangers, the peroxides start fast polymerisation reactions leading to the formation of polymers which increase in chain length as more hydrocarbons are attached. The oxygen source is typically from air in non-blanketed storage tanks or oxygenated compounds in the feed stream, which become more reactive as the feed stream is heated [6, 11].

At lower temperatures, free radicals may be formed when a ligand is broken from a metal complex or salt. The unshared electrons resulting from this break react with an unsaturated hydrocarbon or oxygen to form a free radical [5]. There are numerous transition metals which, in very low concentrations, can act as a catalyst and initiate polymerisation reactions. Some of these catalytically reactive metals are iron, copper, nickel, vanadium, chromium, calcium, and magnesium [6]

68 MATLAB – A Fundamental Tool for Scientific Computing and Engineering Applications – Volume 3

deposit to return the heat exchanger to satisfactory operation.

radical and non-free radical polymerisation.

(Table 2).

fouling associated with the "A" heat exchanger, in the Homs hydrotreating plant, fouling by corrosion is not negligible, with about 19% of the total fouling may be due to corrosion, as is clearly indicated by the iron content of the deposits obtained from this exchanger

**Coking and Polymerisation** are major causes of fouling in heat exchangers. Decomposition of organic products can lead to the formation of very viscous tar or solid coke particles at high temperatures and polymerisation involves the formation of undesirable organic sediments or polymers. The coke particles and polymers formed in the heat exchanger may grow to such a large size that they drop out of solution and deposit on the process equipment. Such deposits can be extremely tenacious and may require burning off the

There are two major polymerisation mechanisms which can occur in the feed stream: free

Free radical polymerisation occurs when a free radical is formed and continues to react with other molecules. The free radicals continue to propagate in the feed stream producing longer chain polymers which will continue to be produced as long as free radicals are being formed. Free radical polymerisation is easily initiated in the presence of light and heat and its rate for polymer formation increases exponentially with temperature. A general rule is that for every 10ºC increase in temperature the rate of polymer formation doubles. Free radical polymerisation readily takes place in heat exchanger tubes and storage tanks [6].

The formation of free radicals has been investigated extensively and it is known that numerous types of free radicals can be formed in a feed stream. These include alkyl radicals produced by the breaking of double or unsaturated bonds as well as other types of precursors such as nitrogen and sulphur radicals which arc easily formed at the temperatures found in the heat exchanger train. Organic sulphur, nitrogen and oxygen compounds increase the potential for various polymerisation reactions, depending on the form in which they exist. Acidic compounds can promote free radical polymerisation by initiating free radicals through the formation of a positive ion or cation. Additional

Oxygen may also react with hydrocarbons to form peroxide free radicals, a step that could occur in the storage tank. When the temperature is increased in the heat exchangers, the peroxides start fast polymerisation reactions leading to the formation of polymers which increase in chain length as more hydrocarbons are attached. The oxygen source is typically from air in non-blanketed storage tanks or oxygenated compounds in the feed stream,

At lower temperatures, free radicals may be formed when a ligand is broken from a metal complex or salt. The unshared electrons resulting from this break react with an unsaturated hydrocarbon or oxygen to form a free radical [5]. There are numerous transition metals which, in very low concentrations, can act as a catalyst and initiate polymerisation reactions.

polymerisation precursors include carbonyls, mercaptans, and pyrrole nitrogen.

which become more reactive as the feed stream is heated [6, 11].

In non-free radical polymerisation, polymer formation results from the reaction of two different molecules under the right conditions. One of the reactive molecules may be a radical, or a compound from a free radical-initiated polymerisation step. Basic compounds can react with other compounds or with themselves to form polymers by several different polymerisation mechanisms. In condensation polymerisation, two large radicals or compounds react together to form an even larger compound, but in their reaction also generate a smaller compound, such as water. This new larger compound can continue to react with other reactive species in the feed stream to make higher molecular weight polymers. At some point, the polymer will either get so large in size that it is no longer able to stay entrained or soluble in the fluid stream and deposit, or all the different compounds that can react with it are consumed, and no further polymer is formed [6].

Various laboratory tests can provide an indication of a stream's polymerisation potential. These include laboratory simulations and analytical characterisation, to identify specific compounds in the feed which are known to contribute to polymerisation mechanisms. Such polymerisation precursors may include unsaturated hydrocarbons, acidic compounds, amines, carbonyls, mercaptans and pyrrole nitrogen.

The presence of unsaturated components in the feed stream contributes significantly to polymerisation reactions, particularly at high temperatures. The bromine number is a method of measuring the degree of unsaturation in a feed stream. The unsaturated bonds react with bromine, and the amount of bromine reacted is an indication of the degree of unsaturation [12].

The neutralisation, or acid, number measures the acidity of the fluid as it is titrated with a base. This number can be an indication of fouling tendency, where the more acidic the feed stream, the greater is its tendency to foul. This is most likely due to the fact that acidic compounds, as mentioned above, can promote free radical polymerisation.

The basic nitrogen test determines the amount of basic compounds in a sample, assumed to be mostly amines, by titrating with a mixture of organic acids. This method can, however, overestimate the basic nitrogen content.

A method of determining a sample's oxidative polymerisation potential is to run a potential gums test. This test is a method of determining a sample's oxidative polymerisation potential. In this test, the fluid is subjected to 100% oxygen for four hours, at 100ºC, in a pressurised sample bomb. The measured gum content, as compared to an initial gum value, will indicate the impact of oxygen on the stream's polymerisation potential.

Detailed deposit analysis, as mentioned earlier, can also indicate the occurrence of polymerisation. It is apparent from examination of the deposit analysis results shown in Table 2 that most deposits are organic in nature, as the loss reported on heating the deposit samples to 840°C was greater than 80% in both the "B" and "C" heat exchangers, where working temperatures are rather high. Since organic deposits result mainly from polymerisation reactions, the high organic content observed in the deposit analysis could be taken as an indication that the fouling in these two heat exchangers is due mainly to polymerisation, which could take place in the heat exchangers themselves or it could occur prior to the heat exchangers either during storage or in transport. Analysis for metals in the deposits indicates the presence of individual metals in the stream. Although, some of these metals are only found in very low concentrations, this may be sufficient for catalysing different polymerisation mechanisms [7].

Fouling in Heat Exchangers 71

In the diffusion regime, the smaller the particle size, the greater is its propensity to be deposited. Thus, it is precisely the very fine submicron particles that are most difficult to remove by filtration or other means which have the greatest propensity to foul a surface.

The transition from diffusional to inertial control of transport occurs at particle diameter in the order of 1–2 μm. In the inertia regime the particles are sufficiently large that turbulent eddies give some of them a transverse (Free flight) velocity which is not completely dissipated in the viscous sublayer. These particles then possess sufficient momentum to reach the wall. Some of the particles also experience a more gradual movement towards the wall by migration down the turbulence intensity gradient, i.e. by "turbophoresis" [15]. Much work has been done by a large number of investigators on predicting the results of this free

In the inertial regime a more desirable situation prevails. Here the larger particles, which are relatively easy to remove, are those which have the greater propensity to be

In this regime, which starts at particle diameter dp 10–20 μm, the particle velocity towards the wall approaches the friction velocity and the particle stopping distance becomes of the same order as the pipe diameter. The response of such large particles to turbulent fluctuations becomes limited and the transport coefficient therefore levels off. As the particles get still larger they get even more sluggish in their response to turbulent eddies

In the impaction regime, transport-controlled deposition would be virtually independent of

There is considerable experimental evidence to indicate that the effect of surface roughness is usually to enhance the transport of particles to the surface. The enhancement occurs because of the decrease of viscous sublayer thickness and corresponding increase in turbulence level above the roughness elements, because of the smaller stopping distance required for the particles to arrive at the outer asperities of the roughness elements, and because of the additional mechanism of particle interception by those elements along flow

On the other hand, turbulent particle transport may be retarded as a result of deposition of very fine particles which tends to smooth initially rough surfaces. Transport-retardation is, however, far less common than transport-enhancement by surface roughness. The importance of clean and, where feasible also, polished surfaces for mitigating particle

flight excursion or inertial coasting in a turbulent field [16].

and the transport coefficient actually starts to fall gradually [14].

deposition under transport-controlled conditions is thus apparent [15].

*2.4.2. Inertia* 

deposited.

*2.4.3. Impaction* 

particle size.

lines parallel to the macrosurface [15].

## **2.3. Aggregation and flocculation**

Some of the heavy organics, especially asphaltenes, will separate from the oil phase into large particles or aggregates. These aggregates may then remain in the oil by some peptising agents, like resins, which will be adsorbed on their surface and keep them afloat, but the stability of such steric colloids is a function of concentration of the peptising agent in the solution. When this concentration drops to a point at which its adsorbed amount is not high enough to cover the entire surface of heavy organic particles, these particles coalesce together, grow in size and flocculate. Flocculation of asphaltene in paraffinic crude oils is known to be irreversible. Due to their large size and their adsorption affinity to solid surfaces flocculated asphaltenes can cause irreversible deposition. Segments of the separated particles which contain S, N and/or H bonds could also start to flocculate and as a result produce the irreversible heavy organic deposits which may be insoluble in solvents.

Inorganic particles may also act as nuclei on which agglomeration of organic particles proceed until the particles become eventually large enough to drop out.

## **2.4. Transport and migration to the fouling sites**

Starting with submicron particles, three transport mechanisms progressively predominate in turbulent flow as the particle size increases. After Gudmunsson [13], the corresponding regimes are designated simply as diffusion, inertia and impaction, respectively.

### *2.4.1. Diffusion*

In the diffusion regime, suspended colloidal particles i.e., particles smaller than about 1 μm in at least one dimension, move with the fluid and are carried to the wall by the Brownian motion of the fluid molecules and through the viscous sublayer in the case of a turbulent flow. The submicron particles can then be treated like large molecules, so that the transport coefficient becomes equivalent to the conventional mass transfer coefficient, which can be obtained from the relevant empirical correlations or theoretical equations for forced convection mass transfer in the literature [14].

In the diffusion regime, the smaller the particle size, the greater is its propensity to be deposited. Thus, it is precisely the very fine submicron particles that are most difficult to remove by filtration or other means which have the greatest propensity to foul a surface.

## *2.4.2. Inertia*

70 MATLAB – A Fundamental Tool for Scientific Computing and Engineering Applications – Volume 3

different polymerisation mechanisms [7].

**2.3. Aggregation and flocculation** 

which may be insoluble in solvents.

*2.4.1. Diffusion* 

samples to 840°C was greater than 80% in both the "B" and "C" heat exchangers, where working temperatures are rather high. Since organic deposits result mainly from polymerisation reactions, the high organic content observed in the deposit analysis could be taken as an indication that the fouling in these two heat exchangers is due mainly to polymerisation, which could take place in the heat exchangers themselves or it could occur prior to the heat exchangers either during storage or in transport. Analysis for metals in the deposits indicates the presence of individual metals in the stream. Although, some of these metals are only found in very low concentrations, this may be sufficient for catalysing

Some of the heavy organics, especially asphaltenes, will separate from the oil phase into large particles or aggregates. These aggregates may then remain in the oil by some peptising agents, like resins, which will be adsorbed on their surface and keep them afloat, but the stability of such steric colloids is a function of concentration of the peptising agent in the solution. When this concentration drops to a point at which its adsorbed amount is not high enough to cover the entire surface of heavy organic particles, these particles coalesce together, grow in size and flocculate. Flocculation of asphaltene in paraffinic crude oils is known to be irreversible. Due to their large size and their adsorption affinity to solid surfaces flocculated asphaltenes can cause irreversible deposition. Segments of the separated particles which contain S, N and/or H bonds could also start to flocculate and as a result produce the irreversible heavy organic deposits

Inorganic particles may also act as nuclei on which agglomeration of organic particles

Starting with submicron particles, three transport mechanisms progressively predominate in turbulent flow as the particle size increases. After Gudmunsson [13], the corresponding

In the diffusion regime, suspended colloidal particles i.e., particles smaller than about 1 μm in at least one dimension, move with the fluid and are carried to the wall by the Brownian motion of the fluid molecules and through the viscous sublayer in the case of a turbulent flow. The submicron particles can then be treated like large molecules, so that the transport coefficient becomes equivalent to the conventional mass transfer coefficient, which can be obtained from the relevant empirical correlations or theoretical equations for forced

proceed until the particles become eventually large enough to drop out.

regimes are designated simply as diffusion, inertia and impaction, respectively.

**2.4. Transport and migration to the fouling sites** 

convection mass transfer in the literature [14].

The transition from diffusional to inertial control of transport occurs at particle diameter in the order of 1–2 μm. In the inertia regime the particles are sufficiently large that turbulent eddies give some of them a transverse (Free flight) velocity which is not completely dissipated in the viscous sublayer. These particles then possess sufficient momentum to reach the wall. Some of the particles also experience a more gradual movement towards the wall by migration down the turbulence intensity gradient, i.e. by "turbophoresis" [15]. Much work has been done by a large number of investigators on predicting the results of this free flight excursion or inertial coasting in a turbulent field [16].

In the inertial regime a more desirable situation prevails. Here the larger particles, which are relatively easy to remove, are those which have the greater propensity to be deposited.

## *2.4.3. Impaction*

In this regime, which starts at particle diameter dp 10–20 μm, the particle velocity towards the wall approaches the friction velocity and the particle stopping distance becomes of the same order as the pipe diameter. The response of such large particles to turbulent fluctuations becomes limited and the transport coefficient therefore levels off. As the particles get still larger they get even more sluggish in their response to turbulent eddies and the transport coefficient actually starts to fall gradually [14].

In the impaction regime, transport-controlled deposition would be virtually independent of particle size.

There is considerable experimental evidence to indicate that the effect of surface roughness is usually to enhance the transport of particles to the surface. The enhancement occurs because of the decrease of viscous sublayer thickness and corresponding increase in turbulence level above the roughness elements, because of the smaller stopping distance required for the particles to arrive at the outer asperities of the roughness elements, and because of the additional mechanism of particle interception by those elements along flow lines parallel to the macrosurface [15].

On the other hand, turbulent particle transport may be retarded as a result of deposition of very fine particles which tends to smooth initially rough surfaces. Transport-retardation is, however, far less common than transport-enhancement by surface roughness. The importance of clean and, where feasible also, polished surfaces for mitigating particle deposition under transport-controlled conditions is thus apparent [15].

## **2.5. Phase separation**

Separation of solid particles from fluid stream and their eventual deposition onto heat exchanger surfaces may be a result of many physical processes including condensation from gas phase, gravitational settling, crystallisation and electro-kinetic effect.

Fouling in Heat Exchangers 73

and increasing shear stress [18]. When various heavy organic compounds are present in a petroleum fluid, their interactive effects largely determine their collective deposition

Changes in the nature of oil fluids may lead to the precipitation of some heavy hydrocarbons, mainly asphaltenes, exceeding their solubility limits. Aspaltenes precipitation, which may be a major cause of crude unit fouling, is affected by many factors including variations of

The deposition of heavy hydrocarbons is an example of what is known as solidification fouling, another example of which is the solidification of molten ash carried in a furnace

Precipitation fouling can also occur as a result of pressure changes, where the solubility of salts such as calcium sulphate decreases with decreasing pressure. Laboratory tests have further indicated that variations of pressure exerted on a petroleum fluid can cause the

Motion of charged particles in a conduit may lead to the development of electrical potential differences along the conduit. This electrical potential difference could then cause a change in charges of the colloidal particles further down in the conduit, the ultimate result of which is their untimely deposition. The factors influencing this effect are the electrical and thermal characteristics of the conduit, flow regime, flowing oil properties, characteristics of the polar

Deposition and attachment of solid particles on heat exchanger surfaces is a function of several different operating variables which include particle size and concentration, bulk fluid density and bulk fluid velocity through the heat exchanger [4, 19]. Furthermore, the stickiness and attractive or repulsive forces between particles can significantly contribute to the deposition of particles [3]. Organic deposits may also be the result of heavy hydrocarbon particles bound to the metal surfaces by inorganic deposition. Attachment is also a function of the interfacial properties of the fouling material and the roughness and wettability of the surface where the fouling is going to occur. Whereas smooth and nonwetting surfaces may delay fouling, rough surfaces provide "nucleation sites" that encourage the laying down of the initial fouling deposits. Most initially smooth walls would tend to roughen as particle deposition occurred, so that roughness would then have to be taken into account. On the other hand, deposition of very fine particles onto initially rough surfaces can conceivably

Recent studies have shown that particle size and concentration have great impact on every type of particle deposition. The average diameter of particles entrained in the fluid stream may vary widely, between a maximum of over 350 μm and a minimum of less than 90 μm (Table 4). Solid particles which foul heat exchangers range in size from submicron to several

result in filling the roughness cavities, thereby smoothing the surface [15].

especially when one of the interacting heavy organic compounds is asphaltene.

temperature, pressure, composition, flow regime, and wall and electrokinetic effect.

exhaust gas onto a heat exchanger surface.

heavy organics and colloidal particles.

**2.6. Particle deposition** 

deposition of some of its heavy organic contents.

Suspended particles such as sand, silt, clay, and non-oxides may become too large to remain entrained in the flowing fluid stream. If the particles are sufficiently large and/or heavy that gravity controls the deposition process, we then have what is known as sedimentation fouling, which can often be prevented with relative ease by pre-filtration or presedimentation of the offending particles. Sedimentation fouling is strongly affected by fluid velocity, and suspended particles in the process fluids will deposit in low-velocity regions, particularly where the velocity changes quickly, as in heat exchanger water boxes and on the shell side [17]. Wall temperature, on the other hand, may have less effect in general on sedimentation fouling, although a hot wall may cause a deposit to "bake on" and become very hard to remove.

Dissolved inorganic salts in a polydisperse fluid may become supersaturated if any change in temperature, pressure, composition (such as solvent evaporation or degasification or addition of a miscible solvent) or other factors destabilises the fluid. The heavy and/or polar fractions may then separate from the fluid into steric colloids, micelles (= charged groups of molecules), another liquid phase or into a solid precipitate.

The dependence of salt solubility on temperature is often the driving force for precipitation fouling. This temperature dependence may be different for different salts, with salt solubility increasing or decreasing with increasing temperature so that different salts may foul the cooling or heating surfaces depending on their solubility temperature dependence. While for most salts the solubility gets higher with increasing temperatures, there are salts such as calcium sulphate which have retrograde solubility dependence and are therefore less soluble in warm streams. Such salts will crystallise on heat transfer surfaces if the streams encounter a surface at a temperature higher than the saturation temperature of these salts. The calcium sulphate scale is hard and adherent and usually requires vigorous mechanical or chemical treatment to remove it. Other typical scaling problems are calcium and magnesium carbonates and silica deposits.

Crystallisation normally begins at specially-active nucleation sites such as scratches and pits, whereas a scratch-free or a smooth surface can flush salt crystals. Subsequently particle deposit will start and continue to build up as long as the surface in contact with the fluid has a temperature above or below saturation. High fluid velocity, by increasing the attrition, can however reduce the rate of particle deposition and fouling.

The solubility of certain heavy hydrocarbons with high melting points such as paraffin wax and diamondoids depends strongly on temperature. If the temperature is decreased, the heavy hydrocarbons may precipitate in the form of solid crystals. Deposition of paraffin wax in cooled heat exchanger tubes showed an asymptotic behaviour due to decreasing heat flux and increasing shear stress [18]. When various heavy organic compounds are present in a petroleum fluid, their interactive effects largely determine their collective deposition especially when one of the interacting heavy organic compounds is asphaltene.

Changes in the nature of oil fluids may lead to the precipitation of some heavy hydrocarbons, mainly asphaltenes, exceeding their solubility limits. Aspaltenes precipitation, which may be a major cause of crude unit fouling, is affected by many factors including variations of temperature, pressure, composition, flow regime, and wall and electrokinetic effect.

The deposition of heavy hydrocarbons is an example of what is known as solidification fouling, another example of which is the solidification of molten ash carried in a furnace exhaust gas onto a heat exchanger surface.

Precipitation fouling can also occur as a result of pressure changes, where the solubility of salts such as calcium sulphate decreases with decreasing pressure. Laboratory tests have further indicated that variations of pressure exerted on a petroleum fluid can cause the deposition of some of its heavy organic contents.

Motion of charged particles in a conduit may lead to the development of electrical potential differences along the conduit. This electrical potential difference could then cause a change in charges of the colloidal particles further down in the conduit, the ultimate result of which is their untimely deposition. The factors influencing this effect are the electrical and thermal characteristics of the conduit, flow regime, flowing oil properties, characteristics of the polar heavy organics and colloidal particles.

## **2.6. Particle deposition**

72 MATLAB – A Fundamental Tool for Scientific Computing and Engineering Applications – Volume 3

gas phase, gravitational settling, crystallisation and electro-kinetic effect.

molecules), another liquid phase or into a solid precipitate.

and magnesium carbonates and silica deposits.

however reduce the rate of particle deposition and fouling.

Separation of solid particles from fluid stream and their eventual deposition onto heat exchanger surfaces may be a result of many physical processes including condensation from

Suspended particles such as sand, silt, clay, and non-oxides may become too large to remain entrained in the flowing fluid stream. If the particles are sufficiently large and/or heavy that gravity controls the deposition process, we then have what is known as sedimentation fouling, which can often be prevented with relative ease by pre-filtration or presedimentation of the offending particles. Sedimentation fouling is strongly affected by fluid velocity, and suspended particles in the process fluids will deposit in low-velocity regions, particularly where the velocity changes quickly, as in heat exchanger water boxes and on the shell side [17]. Wall temperature, on the other hand, may have less effect in general on sedimentation fouling, although a hot wall may cause a deposit to "bake on" and become

Dissolved inorganic salts in a polydisperse fluid may become supersaturated if any change in temperature, pressure, composition (such as solvent evaporation or degasification or addition of a miscible solvent) or other factors destabilises the fluid. The heavy and/or polar fractions may then separate from the fluid into steric colloids, micelles (= charged groups of

The dependence of salt solubility on temperature is often the driving force for precipitation fouling. This temperature dependence may be different for different salts, with salt solubility increasing or decreasing with increasing temperature so that different salts may foul the cooling or heating surfaces depending on their solubility temperature dependence. While for most salts the solubility gets higher with increasing temperatures, there are salts such as calcium sulphate which have retrograde solubility dependence and are therefore less soluble in warm streams. Such salts will crystallise on heat transfer surfaces if the streams encounter a surface at a temperature higher than the saturation temperature of these salts. The calcium sulphate scale is hard and adherent and usually requires vigorous mechanical or chemical treatment to remove it. Other typical scaling problems are calcium

Crystallisation normally begins at specially-active nucleation sites such as scratches and pits, whereas a scratch-free or a smooth surface can flush salt crystals. Subsequently particle deposit will start and continue to build up as long as the surface in contact with the fluid has a temperature above or below saturation. High fluid velocity, by increasing the attrition, can

The solubility of certain heavy hydrocarbons with high melting points such as paraffin wax and diamondoids depends strongly on temperature. If the temperature is decreased, the heavy hydrocarbons may precipitate in the form of solid crystals. Deposition of paraffin wax in cooled heat exchanger tubes showed an asymptotic behaviour due to decreasing heat flux

**2.5. Phase separation** 

very hard to remove.

Deposition and attachment of solid particles on heat exchanger surfaces is a function of several different operating variables which include particle size and concentration, bulk fluid density and bulk fluid velocity through the heat exchanger [4, 19]. Furthermore, the stickiness and attractive or repulsive forces between particles can significantly contribute to the deposition of particles [3]. Organic deposits may also be the result of heavy hydrocarbon particles bound to the metal surfaces by inorganic deposition. Attachment is also a function of the interfacial properties of the fouling material and the roughness and wettability of the surface where the fouling is going to occur. Whereas smooth and nonwetting surfaces may delay fouling, rough surfaces provide "nucleation sites" that encourage the laying down of the initial fouling deposits. Most initially smooth walls would tend to roughen as particle deposition occurred, so that roughness would then have to be taken into account. On the other hand, deposition of very fine particles onto initially rough surfaces can conceivably result in filling the roughness cavities, thereby smoothing the surface [15].

Recent studies have shown that particle size and concentration have great impact on every type of particle deposition. The average diameter of particles entrained in the fluid stream may vary widely, between a maximum of over 350 μm and a minimum of less than 90 μm (Table 4). Solid particles which foul heat exchangers range in size from submicron to several

hundred microns. Investigation revealed that shell and tube exchangers are generally plugged by particles above 20 microns. On the other hand, plate fin exchangers, having much narrower slots, can be plugged by particles as small as 2 microns [3]. The deposition mechanism for the smaller particles is Brownian diffusion while for the larger particles (10- 100 μm) it is mainly gravitational settling. At areas of minimum flow velocities, the larger particles in the stream deposit first followed by the smaller particles and the fouling layer starts to build up as a consequence.

Fouling in Heat Exchangers 75

However, fouling mitigation and control is a very complex process and anticipating the likely extent of fouling problems to be encountered with different flow streams is a major difficulty faced alike by designers and operators of heat exchangers. In most cases, optimisation of the design and operational conditions is not possible or at least would not be realistic without a comprehensive modelling of the process backed up by practical observations. Modelling, however, is not an easy process, and the different models available

The use of multiple regression analysis (MRA), which is an extension of simple least squares regression analysis on a set of data, is an excellent means of modelling heat exchanger fouling. A dependent variable, such as the heat exchanger outlet temperature, is regressed against a set of independent variables, temperatures, pressures, and flows, which directly impact the dependent variable. Regression analysis results in a model equation of independent variables that combine to yield the dependent response. Variability in data and the interaction between independent variables is taken into account in the model equation which can be used to predict future performance. The impact of a change, such as antifoulant addition, can then be compared to the predicted response from the model to de-

Fouling mitigation and control require scientific considerations in design and construction. In general, high turbulence, absence of stagnant areas, uniform fluid flow and smooth surfaces reduce fouling and the need for frequent cleaning. In addition, designers of heat exchangers must consider the effects of fouling upon heat exchanger performance during the desired operational lifetime of the heat exchangers. The factors that need to be considered in the designs include the extra surface required to ensure that the heat exchangers will meet process specifications up to shutdown for cleaning, the additional pressure drop expected due to fouling, and the choice of appropriate construction materials. The designers must also consider the mechanical arrangements that may be necessary for

Fouling resistances are different in different designs of heat exchangers. More than 35-40% of heat exchangers employed in global heat transfer processes are of the shell and tube type of heat exchangers. In process industries, more than 90% of heat exchangers used are of the shell and tube type [20]. This is primarily due to the robust construction geometry as well as ease of maintenance and upgrades possible with the shell and tube heat exchangers [1]. Well established procedures for their design and manufacture from a wide variety of materials, as well as availability of codes and standards for design and fabrication and many years of satisfactory service make them first choice in most process industries. However, fouling resistance in the shell and tube heat exchangers are usually much greater than in other types of heat exchangers (Table 5). In the shell side in particular lower fluid flow velocities and low-velocity or stagnant regions, for example in the vicinity of baffles, encourage the

in the literature are generally of limited value and application.

termine how effective the treatment programme is.

fouling inspection or fouling removal and cleaning.

**3.1. Plant design and construction** 

## **2.7. Deposit growth, aging and hardening**

Following particle deposition, deposit growth and consolidation or alternatively autoretardation and erosion, re-entrainment or removal may take place.

The rate of deposition growth and the accumulation of particles on heat exchanger surfaces is a function of the nature of the fouling material, the composition of the fluid stream and other variables such as temperature and flow rate.

With time, the surface deposit strength may increase and the deposit hardens through various processes collectively known as aging such as, for example, polymerisation, recrystallisation and dehydration. Some types of particles can bake on the surface and will become more difficult to remove over time. The toughness of the deposits may be further affected by the presence of asphaltenes, which are highly polar compounds, and which could act as glue and mortar in hardening the deposits. Biological deposits, on the other hand, may weaken with time due to contamination of organisms.

## **2.8. Auto-retardation and erosion or removal**

The decline of particle deposition rate is commonly referred to as auto-retardation. This is a desirable but spontaneous process that is in the main not under the control of the designer or operator. Several mechanisms may account for auto-retardation and the progressive decrease in adherence of particles to the surface including the already referred to possibility of slowing down particle transport in cases where very fine particles fill the roughness cavities of a surface.

Depending on the strength of the deposit, erosion occurs immediately after the first deposit has been laid down. In saw-tooth fouling part of the deposit is detached after a critical residence time or once a critical deposit thickness has been reached. The fouling layer then builds up and breaks off again. Sometimes impurities such as sand or other suspended particles in fluid streams may have a scouring action, which will reduce or remove deposits [13].

## **3. Fouling mitigation, control and removal**

In order to prevent or mitigate the impact of fouling problems, various steps can be taken during plant design and construction and also during plant operation and maintenance. However, fouling mitigation and control is a very complex process and anticipating the likely extent of fouling problems to be encountered with different flow streams is a major difficulty faced alike by designers and operators of heat exchangers. In most cases, optimisation of the design and operational conditions is not possible or at least would not be realistic without a comprehensive modelling of the process backed up by practical observations. Modelling, however, is not an easy process, and the different models available in the literature are generally of limited value and application.

The use of multiple regression analysis (MRA), which is an extension of simple least squares regression analysis on a set of data, is an excellent means of modelling heat exchanger fouling. A dependent variable, such as the heat exchanger outlet temperature, is regressed against a set of independent variables, temperatures, pressures, and flows, which directly impact the dependent variable. Regression analysis results in a model equation of independent variables that combine to yield the dependent response. Variability in data and the interaction between independent variables is taken into account in the model equation which can be used to predict future performance. The impact of a change, such as antifoulant addition, can then be compared to the predicted response from the model to determine how effective the treatment programme is.

## **3.1. Plant design and construction**

74 MATLAB – A Fundamental Tool for Scientific Computing and Engineering Applications – Volume 3

starts to build up as a consequence.

**2.7. Deposit growth, aging and hardening** 

other variables such as temperature and flow rate.

**2.8. Auto-retardation and erosion or removal** 

**3. Fouling mitigation, control and removal** 

cavities of a surface.

remove deposits [13].

retardation and erosion, re-entrainment or removal may take place.

hand, may weaken with time due to contamination of organisms.

hundred microns. Investigation revealed that shell and tube exchangers are generally plugged by particles above 20 microns. On the other hand, plate fin exchangers, having much narrower slots, can be plugged by particles as small as 2 microns [3]. The deposition mechanism for the smaller particles is Brownian diffusion while for the larger particles (10- 100 μm) it is mainly gravitational settling. At areas of minimum flow velocities, the larger particles in the stream deposit first followed by the smaller particles and the fouling layer

Following particle deposition, deposit growth and consolidation or alternatively auto-

The rate of deposition growth and the accumulation of particles on heat exchanger surfaces is a function of the nature of the fouling material, the composition of the fluid stream and

With time, the surface deposit strength may increase and the deposit hardens through various processes collectively known as aging such as, for example, polymerisation, recrystallisation and dehydration. Some types of particles can bake on the surface and will become more difficult to remove over time. The toughness of the deposits may be further affected by the presence of asphaltenes, which are highly polar compounds, and which could act as glue and mortar in hardening the deposits. Biological deposits, on the other

The decline of particle deposition rate is commonly referred to as auto-retardation. This is a desirable but spontaneous process that is in the main not under the control of the designer or operator. Several mechanisms may account for auto-retardation and the progressive decrease in adherence of particles to the surface including the already referred to possibility of slowing down particle transport in cases where very fine particles fill the roughness

Depending on the strength of the deposit, erosion occurs immediately after the first deposit has been laid down. In saw-tooth fouling part of the deposit is detached after a critical residence time or once a critical deposit thickness has been reached. The fouling layer then builds up and breaks off again. Sometimes impurities such as sand or other suspended particles in fluid streams may have a scouring action, which will reduce or

In order to prevent or mitigate the impact of fouling problems, various steps can be taken during plant design and construction and also during plant operation and maintenance. Fouling mitigation and control require scientific considerations in design and construction. In general, high turbulence, absence of stagnant areas, uniform fluid flow and smooth surfaces reduce fouling and the need for frequent cleaning. In addition, designers of heat exchangers must consider the effects of fouling upon heat exchanger performance during the desired operational lifetime of the heat exchangers. The factors that need to be considered in the designs include the extra surface required to ensure that the heat exchangers will meet process specifications up to shutdown for cleaning, the additional pressure drop expected due to fouling, and the choice of appropriate construction materials. The designers must also consider the mechanical arrangements that may be necessary for fouling inspection or fouling removal and cleaning.

Fouling resistances are different in different designs of heat exchangers. More than 35-40% of heat exchangers employed in global heat transfer processes are of the shell and tube type of heat exchangers. In process industries, more than 90% of heat exchangers used are of the shell and tube type [20]. This is primarily due to the robust construction geometry as well as ease of maintenance and upgrades possible with the shell and tube heat exchangers [1]. Well established procedures for their design and manufacture from a wide variety of materials, as well as availability of codes and standards for design and fabrication and many years of satisfactory service make them first choice in most process industries. However, fouling resistance in the shell and tube heat exchangers are usually much greater than in other types of heat exchangers (Table 5). In the shell side in particular lower fluid flow velocities and low-velocity or stagnant regions, for example in the vicinity of baffles, encourage the

accumulation of foulants. Furthermore, segmental baffles have the tendency for poor flow distribution if spacing or baffle cut ratio is not in the correct proportions. Fouling resistance in plate heat exchangers, on the other hand, can be much smaller. This may be due to the high degree of turbulence even at low velocities which keeps solids in suspension. Also, in plate exchangers there are no dead spaces where fluids can stagnate and solids deposit. Furthermore, heat transfer surfaces are generally smooth and plates are built with higherquality materials with no corrosion products to which fouling may adhere. Finally, cleaning of plate heat exchangers is a very simple operation and the interval between cleanings is usually smaller [21]. Hence, the fouling factors required in plate heat exchangers are normally 20-25% of those used in shell and tube exchangers [22]. In certain applications, spiral plate exchangers may be chosen for fouling services, where the scrubbing action of the fluids on the curved surfaces minimises fouling. On the other hand, fouling is one of the major problems in compact heat exchangers, particularly with various fin geometries and fine flow passages that cannot be cleaned mechanically [23].

Fouling in Heat Exchangers 77

produced by the corrugated tube result in tube wall temperatures closer to the bulk fluid

Mounting the heat exchanger vertically can minimise the effect of deposition fouling as gravity would tend to pull the particles out of the heat exchanger away from the heat transfer surface even at low velocity levels. Appropriate orientation of heat exchangers may also make cleaning easier [27]. In fluid allocation, it is usually preferred to allocate the most fouling fluid to the tube side as it is easier to clean the tube interiors than the exteriors and the probability of low-velocity or stagnant regions is less on the tube side. Placing the fouling fluid in the tube side tends also to minimise fouling by allowing better velocity control. The use of concurrent flow instead of counterflow is a strategy that may be resorted

Appropriate choice of construction materials for heat transfer surfaces may be necessary to alleviate fouling problems. For example, the use of low-fouling surfaces such as surfaces implanted with ions, very smooth surfaces or surfaces of low surface energy may be an option for some applications. Surface coatings and treatment, ultraviolet, acoustic, electric and radiation treatment, may further help to alleviate fouling problems. Surface treatment by plastics, vitreous enamel, glass, and some polymers can also minimise the accumulation of deposits [13]. Similarly, if biofouling is expected or encountered, the use of non-ferrous high copper alloys, which are poisonous to some organisms, can discourage the settling of these organisms on the heat transfer surfaces. Alloys containing copper in quantities greater than 70% are effective in preventing or minimising biological fouling, and generally 70% to 90% copper and 30% to 10% nickel are used for this purpose. Copper alloys are however prohibited in high-pressure steam power plant heat exchangers, since

temperature of the working fluids.

**Figure 4.** Helixchanger heat exchanger.

to sometimes in order to control solidification fouling [23].


**Table 5.** Fouling risk and effects for different types of heat exchangers [24, 25]

Over the years, there has been much advance in the design and manufacture of shell and tube heat exchangers with resultant improvements in their fouling behaviour in operation. A striking example of a new design is the Helixchanger heat exchanger (Fig. 4) where the conventional segmental baffle plates are replaced by quadrant shaped baffles arranged at an angle to the tube axis creating a uniform velocity helical flow through the tube bundle. Near plug flow conditions are achieved in a Helixchanger heat exchanger with little back-flow and eddies, often responsible for fouling and corrosion. Low fouling characteristics are provided offering much longer exchanger run lengths between scheduled cleaning of tube bundles. Such run lengths are increased by 2 to 3 times those achieved using the conventionally baffled shell and tube heat exchangers. Heat exchanger performance is maintained at a higher level for longer periods of time with consequent savings in total life cycle costs of owning and operating Helixchanger heat exchanger banks [1].

If fouling is expected on the tube side, some engineers recommend using larger diameter tubes (a minimum of 25 mm OD) [26]. The use of corrugated tubes has been shown to be beneficial in minimising the effects of at least two of the common types of fouling mechanisms, *viz. deposition fouling* because of an enhanced level of turbulence generated at lower velocities, and *chemical fouling* because the enhanced heat transfer coefficients produced by the corrugated tube result in tube wall temperatures closer to the bulk fluid temperature of the working fluids.

**Figure 4.** Helixchanger heat exchanger.

76 MATLAB – A Fundamental Tool for Scientific Computing and Engineering Applications – Volume 3

fine flow passages that cannot be cleaned mechanically [23].

Air

poor

**Table 5.** Fouling risk and effects for different types of heat exchangers [24, 25]

cycle costs of owning and operating Helixchanger heat exchanger banks [1].

cooled Lamella Plate

fin

Poor Very poor

Over the years, there has been much advance in the design and manufacture of shell and tube heat exchangers with resultant improvements in their fouling behaviour in operation. A striking example of a new design is the Helixchanger heat exchanger (Fig. 4) where the conventional segmental baffle plates are replaced by quadrant shaped baffles arranged at an angle to the tube axis creating a uniform velocity helical flow through the tube bundle. Near plug flow conditions are achieved in a Helixchanger heat exchanger with little back-flow and eddies, often responsible for fouling and corrosion. Low fouling characteristics are provided offering much longer exchanger run lengths between scheduled cleaning of tube bundles. Such run lengths are increased by 2 to 3 times those achieved using the conventionally baffled shell and tube heat exchangers. Heat exchanger performance is maintained at a higher level for longer periods of time with consequent savings in total life

If fouling is expected on the tube side, some engineers recommend using larger diameter tubes (a minimum of 25 mm OD) [26]. The use of corrugated tubes has been shown to be beneficial in minimising the effects of at least two of the common types of fouling mechanisms, *viz. deposition fouling* because of an enhanced level of turbulence generated at lower velocities, and *chemical fouling* because the enhanced heat transfer coefficients

good Good Poor Fair Poor Fair Fair Fair Very

Coiled tube

Double

pipe Graphite Scraped

Poor Fair Poor Good

surface

good

Plate Spiral plate

Very

effect Poor Good Good Very

Type

Fouling risk

Fouling

Shell and tube

Very poor

accumulation of foulants. Furthermore, segmental baffles have the tendency for poor flow distribution if spacing or baffle cut ratio is not in the correct proportions. Fouling resistance in plate heat exchangers, on the other hand, can be much smaller. This may be due to the high degree of turbulence even at low velocities which keeps solids in suspension. Also, in plate exchangers there are no dead spaces where fluids can stagnate and solids deposit. Furthermore, heat transfer surfaces are generally smooth and plates are built with higherquality materials with no corrosion products to which fouling may adhere. Finally, cleaning of plate heat exchangers is a very simple operation and the interval between cleanings is usually smaller [21]. Hence, the fouling factors required in plate heat exchangers are normally 20-25% of those used in shell and tube exchangers [22]. In certain applications, spiral plate exchangers may be chosen for fouling services, where the scrubbing action of the fluids on the curved surfaces minimises fouling. On the other hand, fouling is one of the major problems in compact heat exchangers, particularly with various fin geometries and

> Mounting the heat exchanger vertically can minimise the effect of deposition fouling as gravity would tend to pull the particles out of the heat exchanger away from the heat transfer surface even at low velocity levels. Appropriate orientation of heat exchangers may also make cleaning easier [27]. In fluid allocation, it is usually preferred to allocate the most fouling fluid to the tube side as it is easier to clean the tube interiors than the exteriors and the probability of low-velocity or stagnant regions is less on the tube side. Placing the fouling fluid in the tube side tends also to minimise fouling by allowing better velocity control. The use of concurrent flow instead of counterflow is a strategy that may be resorted to sometimes in order to control solidification fouling [23].

> Appropriate choice of construction materials for heat transfer surfaces may be necessary to alleviate fouling problems. For example, the use of low-fouling surfaces such as surfaces implanted with ions, very smooth surfaces or surfaces of low surface energy may be an option for some applications. Surface coatings and treatment, ultraviolet, acoustic, electric and radiation treatment, may further help to alleviate fouling problems. Surface treatment by plastics, vitreous enamel, glass, and some polymers can also minimise the accumulation of deposits [13]. Similarly, if biofouling is expected or encountered, the use of non-ferrous high copper alloys, which are poisonous to some organisms, can discourage the settling of these organisms on the heat transfer surfaces. Alloys containing copper in quantities greater than 70% are effective in preventing or minimising biological fouling, and generally 70% to 90% copper and 30% to 10% nickel are used for this purpose. Copper alloys are however prohibited in high-pressure steam power plant heat exchangers, since

the corrosion deposits of copper alloys are transported and deposited in high-pressure steam generators and subsequently block the turbine blades. Environmental protection also limits the use of copper in river, lake, and ocean waters, since copper is poisonous to aquatic life [23].

Fouling in Heat Exchangers 79

most important question to be answered is, whether the cost of filtration is higher than

Antifoulants or chemical fouling inhibitors may be used to reduce fouling in many systems mainly by preventing reactions causing fouling, and minimising or interfering with the different steps of the fouling process such as crystallisation, agglomeration of small insoluble polymeric or coke-like particles, sticking or attachment of particles to tube walls and deposit consolidation [28]. Such antifoulants include antioxidation additives used to inhibit polymerisation reactions, metal coordinators which react with the trace elements and prevent them from functioning as fouling catalysts, corrosion inhibitors and dispersion agents. Other antifoulants may be used to control crystallisation such as distortion and

Various strategies and devices for the continuous mitigation and reduction of fouling have been proposed such as periodical reversal of flow direction for the removal of weakly adherent deposits, intermittent air injection and/or increasing wall shear stress by raising flow velocity or by increasing turbulence level. In order to enhance the removal of the fouling deposits, velocities in tubes should in general be above 2 m/s and about 1 m/s on the shell side [10]. In several patents, tubular heat exchangers equipped with fouling reduction devices mounted inside the exchanger tubes are described [29]. Such fouling reducers may comprise a mobile turbulence generating element that consists of a metallic winding in the form of a solenoid. The solenoid can be held in position by a hanging system in such a manner that the turbulence generating element can be driven in rotation by the liquid that circulates in the exchanger. The mobile components can be made of spring steel to make them unstretchable. Alternatively, an elastic solenoid may be used that extends over the entire length of the tubes and is agitated by the liquid that circulates in the exchanger [19]. In some heat-transfer applications, mechanical mitigation with dynamic scraped surface heat exchangers is an option. In self-cleaning fluidised-bed exchangers, a fluidised bed of particles is used to control fouling on the outside or inside of tubular exchangers. The selfcleaning exchanger consists of a large number of parallel vertical tubes, in which small solid particles are kept in a fluidised condition by the liquid velocities. The particles have a

dispersion agents, sequestering agents and threshold chemicals [28].

slightly abrasive effect on the tube walls, so that they remove the deposits [30].

Finally, different mechanical strategies for continuous on-stream cleaning of the interior surfaces of the tubes have been proposed including such strategies as circulation of cleaning balls such as sponge rubber or grit coated balls and pushing of brushes through tubes. In the sponge rubber ball cleaning system, the balls used for normal operation should have the right surface roughness to gently clean the tubes without scoring the tube surface. To remove heavy deposits, special abrasive balls that have a coating of carborundum are available [31]. In the Amertap System, slightly oversized sponge rubber balls are continuously recirculated through the tubes in order to remove the accumulation of scale or corrosion products. The M.A.N. System provides for on-stream cleaning by passage of

the fouling cost [3].

brushes through the tubes.

Corrosion-type fouling can also be minimised by the choice of a construction material which does not readily corrode or produce voluminous deposits of corrosion products. A wide range of corrosion resistant materials based on stainless steel is now available to the heat exchanger manufacturer. Noncorrosive but expensive materials such as titanium and nickel based alloys may be used sometimes to prevent corrosion. If one of the fluids is more corrosive, it may be convenient to send it through the tube side because the shell can then be built with a lower-quality and cheaper material.

The construction material selected must also be resistant to attack by the cleaning solutions in situations where chemical removal of the fouling deposit is planned. For fluid allocation, it is usually preferred to allocate the most fouling fluid to the tube side as it is easier to clean the tube interiors than the exteriors.

### **3.2. Plant operation and maintenance**

In many cases, even the right design of a heat exchanger will not prevent fouling problems that may not be predictable at the design stage. For the control and mitigation of fouling it is generally necessary to take into account the different plant operational conditions such as temperature range, fluid flow rate and chemical composition, and, where possible, make such changes as are required by the severity and type of the fouling problems. For example, some types of fouling can be minimised by using high flow velocities, with due consideration of the possibility of metal erosion as it may be necessary to restrict the velocity to values consistent with satisfactory tube life.

Several techniques may be used for the control of fouling as part of plant maintenance. Some of these techniques are designed to prevent or mitigate fouling. These include avoidance of feed contact with air or oxygen by nitrogen blanketing, elimination or reduction of unsaturates, prior treatment of feed, the use of anti-foulants and application of mechanical on-line mitigation strategies. Cathodic protection and surface treatment such as passivation of stainless steel will minimise corrosion fouling [23].

Prior treatment of feed includes caustic scrubbing, desalting, filtration or sedimentation of feed. Caustic scrubbing removes sulphur compounds and desalting reduces trace metal contamination, both of which reduce polymerisation fouling [10]. Depending on system parameters, including fluid temperature, viscosity, pressure, solid concentration, particle size distribution, and fluid compatibility with the filter media, a filter can be designed to remove solid particles from the fluid. Filtration, however, can only remove the larger-sized particles leaving the smaller-sized particles in the feed stream. Filters used on the feed line require also regular maintenance. At the filter design stage, the most important question to be answered is, whether the cost of filtration is higher than the fouling cost [3].

78 MATLAB – A Fundamental Tool for Scientific Computing and Engineering Applications – Volume 3

aquatic life [23].

built with a lower-quality and cheaper material.

the tube interiors than the exteriors.

**3.2. Plant operation and maintenance** 

to values consistent with satisfactory tube life.

passivation of stainless steel will minimise corrosion fouling [23].

the corrosion deposits of copper alloys are transported and deposited in high-pressure steam generators and subsequently block the turbine blades. Environmental protection also limits the use of copper in river, lake, and ocean waters, since copper is poisonous to

Corrosion-type fouling can also be minimised by the choice of a construction material which does not readily corrode or produce voluminous deposits of corrosion products. A wide range of corrosion resistant materials based on stainless steel is now available to the heat exchanger manufacturer. Noncorrosive but expensive materials such as titanium and nickel based alloys may be used sometimes to prevent corrosion. If one of the fluids is more corrosive, it may be convenient to send it through the tube side because the shell can then be

The construction material selected must also be resistant to attack by the cleaning solutions in situations where chemical removal of the fouling deposit is planned. For fluid allocation, it is usually preferred to allocate the most fouling fluid to the tube side as it is easier to clean

In many cases, even the right design of a heat exchanger will not prevent fouling problems that may not be predictable at the design stage. For the control and mitigation of fouling it is generally necessary to take into account the different plant operational conditions such as temperature range, fluid flow rate and chemical composition, and, where possible, make such changes as are required by the severity and type of the fouling problems. For example, some types of fouling can be minimised by using high flow velocities, with due consideration of the possibility of metal erosion as it may be necessary to restrict the velocity

Several techniques may be used for the control of fouling as part of plant maintenance. Some of these techniques are designed to prevent or mitigate fouling. These include avoidance of feed contact with air or oxygen by nitrogen blanketing, elimination or reduction of unsaturates, prior treatment of feed, the use of anti-foulants and application of mechanical on-line mitigation strategies. Cathodic protection and surface treatment such as

Prior treatment of feed includes caustic scrubbing, desalting, filtration or sedimentation of feed. Caustic scrubbing removes sulphur compounds and desalting reduces trace metal contamination, both of which reduce polymerisation fouling [10]. Depending on system parameters, including fluid temperature, viscosity, pressure, solid concentration, particle size distribution, and fluid compatibility with the filter media, a filter can be designed to remove solid particles from the fluid. Filtration, however, can only remove the larger-sized particles leaving the smaller-sized particles in the feed stream. Filters used on the feed line require also regular maintenance. At the filter design stage, the Antifoulants or chemical fouling inhibitors may be used to reduce fouling in many systems mainly by preventing reactions causing fouling, and minimising or interfering with the different steps of the fouling process such as crystallisation, agglomeration of small insoluble polymeric or coke-like particles, sticking or attachment of particles to tube walls and deposit consolidation [28]. Such antifoulants include antioxidation additives used to inhibit polymerisation reactions, metal coordinators which react with the trace elements and prevent them from functioning as fouling catalysts, corrosion inhibitors and dispersion agents. Other antifoulants may be used to control crystallisation such as distortion and dispersion agents, sequestering agents and threshold chemicals [28].

Various strategies and devices for the continuous mitigation and reduction of fouling have been proposed such as periodical reversal of flow direction for the removal of weakly adherent deposits, intermittent air injection and/or increasing wall shear stress by raising flow velocity or by increasing turbulence level. In order to enhance the removal of the fouling deposits, velocities in tubes should in general be above 2 m/s and about 1 m/s on the shell side [10]. In several patents, tubular heat exchangers equipped with fouling reduction devices mounted inside the exchanger tubes are described [29]. Such fouling reducers may comprise a mobile turbulence generating element that consists of a metallic winding in the form of a solenoid. The solenoid can be held in position by a hanging system in such a manner that the turbulence generating element can be driven in rotation by the liquid that circulates in the exchanger. The mobile components can be made of spring steel to make them unstretchable. Alternatively, an elastic solenoid may be used that extends over the entire length of the tubes and is agitated by the liquid that circulates in the exchanger [19]. In some heat-transfer applications, mechanical mitigation with dynamic scraped surface heat exchangers is an option. In self-cleaning fluidised-bed exchangers, a fluidised bed of particles is used to control fouling on the outside or inside of tubular exchangers. The selfcleaning exchanger consists of a large number of parallel vertical tubes, in which small solid particles are kept in a fluidised condition by the liquid velocities. The particles have a slightly abrasive effect on the tube walls, so that they remove the deposits [30].

Finally, different mechanical strategies for continuous on-stream cleaning of the interior surfaces of the tubes have been proposed including such strategies as circulation of cleaning balls such as sponge rubber or grit coated balls and pushing of brushes through tubes. In the sponge rubber ball cleaning system, the balls used for normal operation should have the right surface roughness to gently clean the tubes without scoring the tube surface. To remove heavy deposits, special abrasive balls that have a coating of carborundum are available [31]. In the Amertap System, slightly oversized sponge rubber balls are continuously recirculated through the tubes in order to remove the accumulation of scale or corrosion products. The M.A.N. System provides for on-stream cleaning by passage of brushes through the tubes.

Notwithstanding the various control and maintenance techniques which can minimise fouling problems and reduce their severity, fouling may still occur and fouling removal and process equipment cleaning may be necessary. A review of the patent activities related to fouling indicates in fact that most of the work deals with fouling removal and process equipment cleaning techniques. This could also mean that many process equipment manufacturers face problems after they appear rather than proactively prevent them from occurring [3].

Fouling in Heat Exchangers 81

**4. Rate of fouling** 

for description of the amount of fouling.

the flowing fluid remains constant.

Fouling Rate = (deposition term) - (anti-deposition term)

"threshold conditions" in a given heat exchanger.

term and a mitigation term.

falling, accelerating, asymptotic or saw-tooth as the case may be.

The rate of fouling is normally defined as the average deposit surface loading per unit of surface area in a unit of time. Deposit thickness (μm) and porosity (%) are also often used

Depending on the fouling mechanism and conditions, the rate of fouling may be linear,

1. Linear fouling is the type of fouling where the fouling rate can be steady with time with increasing fouling resistance and deposit thickness. This is perhaps the most common type of fouling. It occurs in general where the temperature of the deposit in contact with

Ebert and Panchal [32] have presented a fouling model that expressed the average (linear) fouling rate under given conditions as a result of two competing terms, namely, a deposition

where α, β, γ and δ are parameters determined by regression, τw is the shear stress at the tube wall and Tfilm is the fluid film temperature (average of the local bulk fluid and local wall temperatures). The relationship in Eq. (1) points to the possibility of identifying combinations of temperature and velocity below which the fouling rates will be negligible. Ebert and Panchal [32] present this as the "threshold condition". The model in Eq. (1) suggests that the heat exchanger geometry which affects the surface and film temperatures, velocities and shear stresses can be effectively applied to maintain the conditions below the

2. Falling fouling is the type of fouling where the fouling rate decreases with time, and the deposit thickness does not achieve a constant value, although the fouling rate never drops below a certain minimum value. Falling fouling in general is due to an increase of removal rate with time. Its progress can often be described by two numbers: the initial

3. Accelerating fouling is the type of fouling where the fouling rate increases with time. It is the result of hard and adherent deposit where removal and aging can be ignored. It can develop when fouling increases the surface roughness, or when the deposit surface

exhibits higher chemical propensity to fouling than the pure underlying metal. 4. Asymptotic fouling rate is where rate decreases with time until it becomes negligible after a period of time when the deposition rate becomes equal to the deposit removal rate and the deposit thickness remains constant. In general, this type of fouling occurs

*w*

(1)

*film*

Re Pr exp *<sup>f</sup>*

fouling rate and the fouling rate after a long period of time.

*dR E dt RT* 

There are several different techniques that can be employed for the removal of fouling. All such techniques require, however, costly system shutdown after a longer period of low efficiency heat transfer. The chief techniques normally utilised are either chemical or mechanical cleaning, but other procedures may sometimes be employed for some specific applications such as ultrasonic cleaning, which is a more recent procedure, and abrasive cleaning.

Mechanical cleaning is generally preferred over chemical cleaning because it can be a more environmentally-friendly alternative, whereas chemical cleaning causes environmental problems through the handling, application, storage and disposal of chemicals. However, mechanical cleaning may damage the equipment, particularly tubes, and it does not produce a chemically clean surface. Furthermore, chemical cleaning may be the only alternative if uniform or complete cleaning is required and for cleaning inaccessible areas. The shell side in particular can only be chemically cleaned. The tubes on the other hand can be mechanically cleaned provided that the tube pattern and pitch provide sufficient space and access to the inside of the bundle, and if mechanical cleaning is required for one of the fluids, the usual practice is to put that fluid in the tube side.

For the chemical removal of fouling material, weak acids and special solvents or detergents are normally used. Chlorination may be used for the removal of carbonate deposits. Mechanical techniques for the removal of fouling include scraping and air bumping. Air bumping is a technique that involves the creation of slugs of air, thereby creating localised turbulence as slugs pass through the equipment.

For tightly plugged tubes drilling, generally known as bulleting, may be employed and for lightly plugged tubes roding is employed. Particularly weakly adherent deposits may be mechanically removed by applying high velocity water jets or a mixture of sand and water. Jet cleaning can be used mostly on external surfaces where there is an easy accessibility for passing the high pressure jet [23].

In cases where biofouling occurs it may be removed by either chemical treatment or mechanical brushing processes. In chemical cleaning techniques biocides are employed such as chlorine, chlorine dioxide, bromine, ozone and surfactants. A more usual practice, however, is by continuous or intermittent "shock" chlorination which kills off the responsible organisms. Other cleaning techniques that can be effective in controlling biological fouling include thermal shock treatment by application of heat or deslugging with steam or hot water, and some less well-known techniques like ultraviolet radiation [23].

### **4. Rate of fouling**

80 MATLAB – A Fundamental Tool for Scientific Computing and Engineering Applications – Volume 3

occurring [3].

cleaning.

Notwithstanding the various control and maintenance techniques which can minimise fouling problems and reduce their severity, fouling may still occur and fouling removal and process equipment cleaning may be necessary. A review of the patent activities related to fouling indicates in fact that most of the work deals with fouling removal and process equipment cleaning techniques. This could also mean that many process equipment manufacturers face problems after they appear rather than proactively prevent them from

There are several different techniques that can be employed for the removal of fouling. All such techniques require, however, costly system shutdown after a longer period of low efficiency heat transfer. The chief techniques normally utilised are either chemical or mechanical cleaning, but other procedures may sometimes be employed for some specific applications such as ultrasonic cleaning, which is a more recent procedure, and abrasive

Mechanical cleaning is generally preferred over chemical cleaning because it can be a more environmentally-friendly alternative, whereas chemical cleaning causes environmental problems through the handling, application, storage and disposal of chemicals. However, mechanical cleaning may damage the equipment, particularly tubes, and it does not produce a chemically clean surface. Furthermore, chemical cleaning may be the only alternative if uniform or complete cleaning is required and for cleaning inaccessible areas. The shell side in particular can only be chemically cleaned. The tubes on the other hand can be mechanically cleaned provided that the tube pattern and pitch provide sufficient space and access to the inside of the bundle, and if mechanical cleaning is required for one of the

For the chemical removal of fouling material, weak acids and special solvents or detergents are normally used. Chlorination may be used for the removal of carbonate deposits. Mechanical techniques for the removal of fouling include scraping and air bumping. Air bumping is a technique that involves the creation of slugs of air, thereby creating localised

For tightly plugged tubes drilling, generally known as bulleting, may be employed and for lightly plugged tubes roding is employed. Particularly weakly adherent deposits may be mechanically removed by applying high velocity water jets or a mixture of sand and water. Jet cleaning can be used mostly on external surfaces where there is an easy accessibility for

In cases where biofouling occurs it may be removed by either chemical treatment or mechanical brushing processes. In chemical cleaning techniques biocides are employed such as chlorine, chlorine dioxide, bromine, ozone and surfactants. A more usual practice, however, is by continuous or intermittent "shock" chlorination which kills off the responsible organisms. Other cleaning techniques that can be effective in controlling biological fouling include thermal shock treatment by application of heat or deslugging with steam or hot water, and some less well-known techniques like ultraviolet radiation [23].

fluids, the usual practice is to put that fluid in the tube side.

turbulence as slugs pass through the equipment.

passing the high pressure jet [23].

The rate of fouling is normally defined as the average deposit surface loading per unit of surface area in a unit of time. Deposit thickness (μm) and porosity (%) are also often used for description of the amount of fouling.

Depending on the fouling mechanism and conditions, the rate of fouling may be linear, falling, accelerating, asymptotic or saw-tooth as the case may be.

1. Linear fouling is the type of fouling where the fouling rate can be steady with time with increasing fouling resistance and deposit thickness. This is perhaps the most common type of fouling. It occurs in general where the temperature of the deposit in contact with the flowing fluid remains constant.

Ebert and Panchal [32] have presented a fouling model that expressed the average (linear) fouling rate under given conditions as a result of two competing terms, namely, a deposition term and a mitigation term.

Fouling Rate = (deposition term) - (anti-deposition term)

$$\frac{d\mathbf{R}\_f}{dt} = \alpha \operatorname{Re}^\beta \operatorname{Pr}^\delta \exp\left(\frac{-E}{RT\_{film}}\right) - \gamma \tau\_w \tag{1}$$

where α, β, γ and δ are parameters determined by regression, τw is the shear stress at the tube wall and Tfilm is the fluid film temperature (average of the local bulk fluid and local wall temperatures). The relationship in Eq. (1) points to the possibility of identifying combinations of temperature and velocity below which the fouling rates will be negligible. Ebert and Panchal [32] present this as the "threshold condition". The model in Eq. (1) suggests that the heat exchanger geometry which affects the surface and film temperatures, velocities and shear stresses can be effectively applied to maintain the conditions below the "threshold conditions" in a given heat exchanger.


where the tube surface temperature remains constant while the temperature of the flowing fluid drops as a result of increased resistance of fouling material to heat transfer. Asymptotic fouling may also be the result of soft or poorly adherent suspended solid deposits upon heat transfer surfaces in areas of fast flow where they do not adhere strongly to the surface with the result that the thicker the deposit becomes, the more likely it is to wash off in patches and thus attain some average asymptotic value over a period of time.

Fouling in Heat Exchangers 83

tabulations available that provide typical fouling factors such as TEMA Table RGP-T-2.4 [35], acceptable evaluation of the effects of fouling needs to be judged and evaluated for each particular application. Such tabulations, however, can be used as a guide in the absence

A number of methods empirical or otherwise have been proposed over the years for the prediction of the rate of fouling in heat exchangers or for estimating a fouling factor to be used in heat transfer calculations. With the advent and development of digital computers with their ability to provide rapid means of performing calculations, new and accurate methods became possible. Matlab being a programming language widely used in all scientific fields and in engineering sciences in particular may be used in conjunction with the artificial neural network (ANN) approach to provide an accurate and reliable method for predicting the fouling rate and rate of heat transfer in heat exchangers. The development of high speed digital computers has provided a rapid means of performing the many calculations involved in the ANN method, and has had a stimulating effect on the current

In recent years, the ANN method has been applied in many disciplines of engineering and has produced promising results. The main feature of this method is its ability to learn and generalise the relationships in a data set and to provide quick and satisfactory estimations,

The artificial neural network method is a computational structure inspired by a biological neural system. An ANN consists of very simple and highly interconnected processors called neurons. The neurons are connected to each other by weighted links over which signals can pass. Each neuron receives multiple inputs from other neurons in proportion to their connection weights and generates a single output, which may be propagated to several other neurons [36]. The inputs (X) into a neuron are multiplied by their corresponding connection weights (W) and summed together, a threshold (θ), acting at a bias, is added also to the sum. This sum is transformed through a transfer function (f) to produce a single output (Y), which may be passed on to other neurons. The function of a neuron can be

Where the transfer function (f) of the neuron is the linear activation function, being in the

Y = f (∑*wx* - θ) (2)

f(x) = purelin (x) (3)

f(x) = tansig (x) (4)

expansion of the ANN method which is progressing at an impressive rate.

which make it attractive for many different applications.

of more specific information.

mathematically expressed as

for the output layer and the tansing function

present work given as:

for the hidden layers

The asymptotic fouling resistance increases with increasing particle concentration and decreasing fluid bulk temperature, flow velocity, and particle diameter. The asymptotic fouling model was first described by Kern and Seaton [33]. In this model, the competing fouling mechanisms lead to an asymptotic fouling resistance beyond which no further increase in fouling occurs. The Tubular Heat Exchanger Manufacturers Association (TEMA) standards suggest fouling factors for several fluids based upon the asymptotic values. This approach, however, does not address all fouling phenomena as it does not, for example, address fouling at the "hot" end of a crude oil preheat train, since fouling there does not exhibit the asymptotic behaviour.

5. Saw-tooth fouling occurs where part of the deposit is detached after a critical residence time or once a critical deposit thickness has been reached. The fouling layer then builds up and breaks off again. This periodic variation could be due to pressure pulses, spalling, trapping of air inside the surface deposits during shutdowns or other reasons. It often corresponds to the moments of system shutdowns, startups or other transients in operation.

## **5. Prediction of fouling factor**

The effect of fouling, as has been noted above, is to form an essentially solid deposit of low thermal conductivity upon the heat transfer surface, through which heat must be transferred by conduction. But since the thermal conductivity of the fouling layer and its thickness are not generally known, the only possible solution to the heat transfer problem is by the introduction of a fouling factor in order to take into account the additional resistance to heat transfer and make possible the calculation of the overall heat transfer coefficient. A fouling coefficient also is sometimes specified, which is the reciprocal value of the fouling factor.

In carrying out heat transfer calculations, caution needs to be exerted in selecting fouling factors, particularly where fouling resistances completely dominate the thermal design. The influence of uncertainties inherent in fouling factors is generally greater than that of uncertainties in other design parameters such as fluid properties, flow rates and temperatures [34]. A large fouling factor is sometimes adopted as a safety margin to cover uncertainties in fluid properties and even in process knowledge but the use of an excessively large fouling factor will result in an oversized heat exchanger with two or three times more area than is really necessary. Although there are many experience-based tabulations available that provide typical fouling factors such as TEMA Table RGP-T-2.4 [35], acceptable evaluation of the effects of fouling needs to be judged and evaluated for each particular application. Such tabulations, however, can be used as a guide in the absence of more specific information.

A number of methods empirical or otherwise have been proposed over the years for the prediction of the rate of fouling in heat exchangers or for estimating a fouling factor to be used in heat transfer calculations. With the advent and development of digital computers with their ability to provide rapid means of performing calculations, new and accurate methods became possible. Matlab being a programming language widely used in all scientific fields and in engineering sciences in particular may be used in conjunction with the artificial neural network (ANN) approach to provide an accurate and reliable method for predicting the fouling rate and rate of heat transfer in heat exchangers. The development of high speed digital computers has provided a rapid means of performing the many calculations involved in the ANN method, and has had a stimulating effect on the current expansion of the ANN method which is progressing at an impressive rate.

In recent years, the ANN method has been applied in many disciplines of engineering and has produced promising results. The main feature of this method is its ability to learn and generalise the relationships in a data set and to provide quick and satisfactory estimations, which make it attractive for many different applications.

The artificial neural network method is a computational structure inspired by a biological neural system. An ANN consists of very simple and highly interconnected processors called neurons. The neurons are connected to each other by weighted links over which signals can pass. Each neuron receives multiple inputs from other neurons in proportion to their connection weights and generates a single output, which may be propagated to several other neurons [36]. The inputs (X) into a neuron are multiplied by their corresponding connection weights (W) and summed together, a threshold (θ), acting at a bias, is added also to the sum. This sum is transformed through a transfer function (f) to produce a single output (Y), which may be passed on to other neurons. The function of a neuron can be mathematically expressed as

$$\mathbf{Y} = \mathbf{f} \left( \sum \mathbf{w} \mathbf{x} \cdot \mathbf{-} \boldsymbol{\Theta} \right) \tag{2}$$

Where the transfer function (f) of the neuron is the linear activation function, being in the present work given as:

$$\mathbf{f}(\mathbf{x}) = \text{purelin}\,\mathbf{(x)}\tag{3}$$

for the output layer and the tansing function

$$\mathbf{f}(\mathbf{x}) = \text{tansig}\,\mathbf{(x)}\tag{4}$$

for the hidden layers

82 MATLAB – A Fundamental Tool for Scientific Computing and Engineering Applications – Volume 3

value over a period of time.

exhibit the asymptotic behaviour.

**5. Prediction of fouling factor** 

in operation.

where the tube surface temperature remains constant while the temperature of the flowing fluid drops as a result of increased resistance of fouling material to heat transfer. Asymptotic fouling may also be the result of soft or poorly adherent suspended solid deposits upon heat transfer surfaces in areas of fast flow where they do not adhere strongly to the surface with the result that the thicker the deposit becomes, the more likely it is to wash off in patches and thus attain some average asymptotic

The asymptotic fouling resistance increases with increasing particle concentration and decreasing fluid bulk temperature, flow velocity, and particle diameter. The asymptotic fouling model was first described by Kern and Seaton [33]. In this model, the competing fouling mechanisms lead to an asymptotic fouling resistance beyond which no further increase in fouling occurs. The Tubular Heat Exchanger Manufacturers Association (TEMA) standards suggest fouling factors for several fluids based upon the asymptotic values. This approach, however, does not address all fouling phenomena as it does not, for example, address fouling at the "hot" end of a crude oil preheat train, since fouling there does not

5. Saw-tooth fouling occurs where part of the deposit is detached after a critical residence time or once a critical deposit thickness has been reached. The fouling layer then builds up and breaks off again. This periodic variation could be due to pressure pulses, spalling, trapping of air inside the surface deposits during shutdowns or other reasons. It often corresponds to the moments of system shutdowns, startups or other transients

The effect of fouling, as has been noted above, is to form an essentially solid deposit of low thermal conductivity upon the heat transfer surface, through which heat must be transferred by conduction. But since the thermal conductivity of the fouling layer and its thickness are not generally known, the only possible solution to the heat transfer problem is by the introduction of a fouling factor in order to take into account the additional resistance to heat transfer and make possible the calculation of the overall heat transfer coefficient. A fouling coefficient also is sometimes specified, which is the reciprocal value of the fouling factor.

In carrying out heat transfer calculations, caution needs to be exerted in selecting fouling factors, particularly where fouling resistances completely dominate the thermal design. The influence of uncertainties inherent in fouling factors is generally greater than that of uncertainties in other design parameters such as fluid properties, flow rates and temperatures [34]. A large fouling factor is sometimes adopted as a safety margin to cover uncertainties in fluid properties and even in process knowledge but the use of an excessively large fouling factor will result in an oversized heat exchanger with two or three times more area than is really necessary. Although there are many experience-based Among the various existing kinds of ANNS, the back propagation (BP) learning algorithm has become the most popular in engineering applications [37]. The BP algorithm is designed to solve the problem of determining weight values for a multi-layer ANN with feed forward connections from the input layer to the hidden layers and then to the output layer. The algorithm is an iterative gradient algorithm, designed to minimise the mean square error between the predicted output and the desired output. Fig. 5 shows a scheme of a simple example for BP algorithm, and Fig. 6 is a flow chart for the back propagation (BP) learning algorithm.

Fouling in Heat Exchangers 85

The ANN approach has also been used as an alternative and practical technique to evaluate the rate of heat exchange and heat transfer coefficient for tubular heat exchangers [34], fin tube heat exchangers [37], fluid—particle systems [54], heat rate predictions in humid air—water heat exchangers [55], and in other applications [22] including in particular the prediction of the rate of fouling and fouling factor in a shell-

For predicting the rate of fouling and fouling factor in heat exchangers an ANN model can be developed and the available data set used for training the network and verifying its generalisation capability. The rates of heat transfer are then calculated and the input output pairs presented to the network, and the weights adjusted to minimise the error between the network output and the actual value. Once training is complete, predictions from a new set of data may be done using the already trained network. The proposed algorithm is solved by a Matlab computer programme. After the rate of heat transfer is

In order to test the applicability of the ANN model using the back propagation learning algorithm to predict the rate of heat transfer and fouling factor for heat exchangers, it was applied on a tube-shell heat exchanger used as a preheat device for the naphtha feed to the

Operation data of the heat exchanger are collected for this purpose [8]. A total of 73 readings tabulated in Table 6 are used for the ANN method. Empirical correlations for the rate of heat

In developing an ANN model, the available data set is used for training the network and the same data are used to verify the generalisation capability of the network [58]. The input parameters were liquid naphtha feed quantity (m·1), gas reaction quantity (m·2), inlet temperature (t1d), outlet temperature (t2d), hydrogen purity (N), gas molecular weight (M), naphtha specific heat (Cp1), gas specific heat (Cp2) and the output parameter is heat transfer

The rates of heat transfer (Q) for the same data are calculated using Excel programme (actual values). Input—output pairs are presented to the network, and the weights are

Q = m Cp (t1d-t2d) kJ/hr (5)

and-tube heat exchanger [56].

Q = q1+q2

rate (Q).

q1 = m**·**1× Cp1 × (t1d-t2d) (liquid) q2 = m**·**2× Cp2× (t1d-t2d) (gas)

Cp2 = Cp (hydrogen) + Cp (hydrocarbons) Cp (hydrogen) = 4.19 (6.8+0.0006×t) y

Cp (hydrocarbons) = (0.00428×t+1.5606) (1-y)

Cp1 = 0.0045×t+2.0687

y = 2.016×N/(100×M)

calculated, it can be used to estimate the fouling factor.

transfer (Q) is determined [57] as given below by Eq.5:

reactor in the naphtha hydrotreating unit at the Homs oil refinery.

Results of the application of the ANN approach show that it can be used to develop the best configuration in the training period. In general, this approach may however be time consuming; it is nonetheless feasible due to its ability to learn and generalise the complex data set with a wide range of experimental conditions. For exhaustive and fundamental treatment of the ANN technique, the reader is referred to any of the textbooks [38, 39].

**Figure 5.** Scheme of a simple example of a BP algorithm.

Some of the recent applications for energy systems include modelling of appliances, light and space cooling energy consumption [40], solar radiation estimation [41], modelling gasoline consumption [42], modelling of heat pumps [43], performance prediction of solar domestic water heating systems [44], prediction of energy consumption of a passive solar building [45], prediction of temperature profiles in producing oil wells [46] and developing heating, ventilating and air conditioning systems for automobiles [47]. In addition, various applications of artificial neural networks in energy problems have been presented thematically [48—50].

In the field of heat exchangers, the Ann approach has been widely used, particularly in the design and control of heat exchangers [51, 52], and in the simulation of heat exchanger performance [14] and heat transfer analysis for different systems such as air—water spray cooling [53].

The ANN approach has also been used as an alternative and practical technique to evaluate the rate of heat exchange and heat transfer coefficient for tubular heat exchangers [34], fin tube heat exchangers [37], fluid—particle systems [54], heat rate predictions in humid air—water heat exchangers [55], and in other applications [22] including in particular the prediction of the rate of fouling and fouling factor in a shelland-tube heat exchanger [56].

For predicting the rate of fouling and fouling factor in heat exchangers an ANN model can be developed and the available data set used for training the network and verifying its generalisation capability. The rates of heat transfer are then calculated and the input output pairs presented to the network, and the weights adjusted to minimise the error between the network output and the actual value. Once training is complete, predictions from a new set of data may be done using the already trained network. The proposed algorithm is solved by a Matlab computer programme. After the rate of heat transfer is calculated, it can be used to estimate the fouling factor.

In order to test the applicability of the ANN model using the back propagation learning algorithm to predict the rate of heat transfer and fouling factor for heat exchangers, it was applied on a tube-shell heat exchanger used as a preheat device for the naphtha feed to the reactor in the naphtha hydrotreating unit at the Homs oil refinery.

Operation data of the heat exchanger are collected for this purpose [8]. A total of 73 readings tabulated in Table 6 are used for the ANN method. Empirical correlations for the rate of heat transfer (Q) is determined [57] as given below by Eq.5:

$$\mathbf{Q} = \mathbf{m} \, \mathbf{C} \, \mathbf{p} \text{ (t1d-t2d) kJ/hr} \tag{5}$$

Q = q1+q2 q1 = m**·**1× Cp1 × (t1d-t2d) (liquid) q2 = m**·**2× Cp2× (t1d-t2d) (gas) Cp1 = 0.0045×t+2.0687 Cp2 = Cp (hydrogen) + Cp (hydrocarbons) Cp (hydrogen) = 4.19 (6.8+0.0006×t) y Cp (hydrocarbons) = (0.00428×t+1.5606) (1-y) y = 2.016×N/(100×M)

84 MATLAB – A Fundamental Tool for Scientific Computing and Engineering Applications – Volume 3

**Figure 5.** Scheme of a simple example of a BP algorithm.

thematically [48—50].

cooling [53].

algorithm.

Among the various existing kinds of ANNS, the back propagation (BP) learning algorithm has become the most popular in engineering applications [37]. The BP algorithm is designed to solve the problem of determining weight values for a multi-layer ANN with feed forward connections from the input layer to the hidden layers and then to the output layer. The algorithm is an iterative gradient algorithm, designed to minimise the mean square error between the predicted output and the desired output. Fig. 5 shows a scheme of a simple example for BP algorithm, and Fig. 6 is a flow chart for the back propagation (BP) learning

Results of the application of the ANN approach show that it can be used to develop the best configuration in the training period. In general, this approach may however be time consuming; it is nonetheless feasible due to its ability to learn and generalise the complex data set with a wide range of experimental conditions. For exhaustive and fundamental treatment of the ANN technique, the reader is referred to any of the textbooks [38, 39].

Some of the recent applications for energy systems include modelling of appliances, light and space cooling energy consumption [40], solar radiation estimation [41], modelling gasoline consumption [42], modelling of heat pumps [43], performance prediction of solar domestic water heating systems [44], prediction of energy consumption of a passive solar building [45], prediction of temperature profiles in producing oil wells [46] and developing heating, ventilating and air conditioning systems for automobiles [47]. In addition, various applications of artificial neural networks in energy problems have been presented

In the field of heat exchangers, the Ann approach has been widely used, particularly in the design and control of heat exchangers [51, 52], and in the simulation of heat exchanger performance [14] and heat transfer analysis for different systems such as air—water spray In developing an ANN model, the available data set is used for training the network and the same data are used to verify the generalisation capability of the network [58]. The input parameters were liquid naphtha feed quantity (m·1), gas reaction quantity (m·2), inlet temperature (t1d), outlet temperature (t2d), hydrogen purity (N), gas molecular weight (M), naphtha specific heat (Cp1), gas specific heat (Cp2) and the output parameter is heat transfer rate (Q).

The rates of heat transfer (Q) for the same data are calculated using Excel programme (actual values). Input—output pairs are presented to the network, and the weights are

adjusted to minimise the error between the network output and the actual value (this is done at the training step on a number of data (n = 73)). Once training is complete, predictions from a new set of data may be done using the already trained network. The proposed algorithm in this study was solved by a computer programme developed using the Matlab programming language, and all computations were performed with a personal computer. After the rate of heat transfer is calculated, it can be used to estimate the fouling factor.

Fouling in Heat Exchangers 87

Gs = m**·**/as where as = ID×C×B/Pt

A total of 73 values of input date tabulated in Table 7 are used for the ANN method, which were the rate of heat transfer (Q) obtained from the neural network, log mean temperature difference (LMTD), shell volumetric flow Gs, tube volumetric flow Gt, shell Reynolds Res and Prandtl Prs numbers, tube Reynolds Ret and Prandtl Prt numbers, and the output is the

The Fouling factor (f) for the same data is calculated using Excel programme (actual value). Input—output pairs are presented to the network, and the weights are adjusted to minimise the error between the network output and the actual value (this is done at the training step on a number of data (n=73). Once training is complete, predictions from a new set of data may be done using the already trained network. The proposed algorithm in this study was solved by a computer programme developed using Matlab programming language; also all

The purpose of using the ANN model with the considered BP learning algorithm as a practical approach is to test the ability to predict the rate of heat transfer and fouling factor of heat exchanger. For the heat transfer function a network is designed with eight input parameters, namely liquid naphtha feed quantity (m**·**1), gas reaction quantity (m**·**2), inlet temperature (t1d), outlet temperature (t2d), hydrogen purity (N), gas molecular weight (M), naphtha specific heat (Cp1) and gas specific heat (Cp2), and one output parameter, the rate of heat exchange (Q). Model sensitivity was examined for different numbers of hidden layer nodes in the range of (1-5). There is no general rule for selecting the number of a hidden layer. The choice of hidden layer size is a specific problem and, to some extent, depends on the number and the quality of training patterns. The number of neurons in a neural network must be sufficient for correct modelling of the problem, and also, it should be low to ensure

The number of neurons in a hidden layer drastically affects the outcome of the network training. If too few neurons are included, then the network may not be able to learn properly. On the other band, if too many neurons were included, the network would encourage over fitting. Some researchers mentioned that the upper bound for the required number of neurons in the hidden layer should be greater than twice the number of input units. This rule does not guarantee generalisation of the network [59]. Every stage of an ANN problem requires a little trial and error to establish a suitable and stable network for the problem; so many networks are built by changing their parameters in order to reach a

net\_q=newff(minmax(input),[8,10,1],{'tansig','tansig','purelin'},'trainlm');

computations were performed with a personal computer.

Gt= m**·**/at where at = Nt×at'/z

Res = Gs×De/μ (shell) Prs = Cp× μ /K (shell)

Ret = D×Gt/ μ (tube) Prt = Cp× μ /k (tube)

fouling factor (f).

generalisation.

suitable result in these statements:

net\_q.trainParam.show = 5;

**Figure 6.** Flow chart for the back propagation (BP) learning algorithm.

The same steps used to calculate the rate of heat transfer by neural network (design, training) are used to calculate fouling factor for different data (Input-output pairs). The empirical correlations for the fouling factor are determined [57] as given below in Eq. 6:

$$\text{if = Uc-Ud/Uc\*Ud m2hr°C/k}\tag{6}$$

Uc = ho×hi/ho+hi Ud = Q/A×LMTD ho = Jh×k/De×pr0.33 hi = Jh×k/D×pr0.33 LMTD = (T1-t2)–(T2-t1)/ln(T1-t2)/(T2-t1) Gs = m**·**/as where as = ID×C×B/Pt Res = Gs×De/μ (shell) Prs = Cp× μ /K (shell) Gt= m**·**/at where at = Nt×at'/z Ret = D×Gt/ μ (tube) Prt = Cp× μ /k (tube)

86 MATLAB – A Fundamental Tool for Scientific Computing and Engineering Applications – Volume 3

heat transfer is calculated, it can be used to estimate the fouling factor.

**Figure 6.** Flow chart for the back propagation (BP) learning algorithm.

Uc = ho×hi/ho+hi Ud = Q/A×LMTD ho = Jh×k/De×pr0.33 hi = Jh×k/D×pr0.33

LMTD = (T1-t2)–(T2-t1)/ln(T1-t2)/(T2-t1)

The same steps used to calculate the rate of heat transfer by neural network (design, training) are used to calculate fouling factor for different data (Input-output pairs). The empirical correlations for the fouling factor are determined [57] as given below in Eq. 6:

f = Uc-Ud/Uc×Ud m2hr°C/kJ (6)

adjusted to minimise the error between the network output and the actual value (this is done at the training step on a number of data (n = 73)). Once training is complete, predictions from a new set of data may be done using the already trained network. The proposed algorithm in this study was solved by a computer programme developed using the Matlab programming language, and all computations were performed with a personal computer. After the rate of

> A total of 73 values of input date tabulated in Table 7 are used for the ANN method, which were the rate of heat transfer (Q) obtained from the neural network, log mean temperature difference (LMTD), shell volumetric flow Gs, tube volumetric flow Gt, shell Reynolds Res and Prandtl Prs numbers, tube Reynolds Ret and Prandtl Prt numbers, and the output is the fouling factor (f).

> The Fouling factor (f) for the same data is calculated using Excel programme (actual value). Input—output pairs are presented to the network, and the weights are adjusted to minimise the error between the network output and the actual value (this is done at the training step on a number of data (n=73). Once training is complete, predictions from a new set of data may be done using the already trained network. The proposed algorithm in this study was solved by a computer programme developed using Matlab programming language; also all computations were performed with a personal computer.

> The purpose of using the ANN model with the considered BP learning algorithm as a practical approach is to test the ability to predict the rate of heat transfer and fouling factor of heat exchanger. For the heat transfer function a network is designed with eight input parameters, namely liquid naphtha feed quantity (m**·**1), gas reaction quantity (m**·**2), inlet temperature (t1d), outlet temperature (t2d), hydrogen purity (N), gas molecular weight (M), naphtha specific heat (Cp1) and gas specific heat (Cp2), and one output parameter, the rate of heat exchange (Q). Model sensitivity was examined for different numbers of hidden layer nodes in the range of (1-5). There is no general rule for selecting the number of a hidden layer. The choice of hidden layer size is a specific problem and, to some extent, depends on the number and the quality of training patterns. The number of neurons in a neural network must be sufficient for correct modelling of the problem, and also, it should be low to ensure generalisation.

> The number of neurons in a hidden layer drastically affects the outcome of the network training. If too few neurons are included, then the network may not be able to learn properly. On the other band, if too many neurons were included, the network would encourage over fitting. Some researchers mentioned that the upper bound for the required number of neurons in the hidden layer should be greater than twice the number of input units. This rule does not guarantee generalisation of the network [59]. Every stage of an ANN problem requires a little trial and error to establish a suitable and stable network for the problem; so many networks are built by changing their parameters in order to reach a suitable result in these statements:

```
net_q=newff(minmax(input),[8,10,1],{'tansig','tansig','purelin'},'trainlm'); 
net_q.trainParam.show = 5;
```
Fouling in Heat Exchangers 89

Re Gs LMTD Q

<sup>f</sup>Tube

Pr Tube Re Gt Shell

Pr

0.0011 1.48 21626.83 602030 1.62 21123 757392 66.59 7420754 0.0011 1.47 23438.94 646368 1.62 22744 813172 67.64 7897551 0.0008 1.48 30387.8 845085 1.62 29534 1063172 67.38 10427147 0.0007 1.47 31538.99 361645 1.62 30543 1084005 67.95 10937005 0.0015 1.45 21187.03 541667 1.61 19793 681452 73.76 6782209 0.0010 1.46 29257.62 788462 1.62 27492 991935 72.36 9735495 0.0009 1.46 33515.82 878205 1.62 31264 1104839 73.89 10815713 0.0009 1.46 31068.49 836538 1.63 28929 1052419 73.99 10467559 0.0008 1.45 32702.56 853098 1.61 31140 1073253 70.74 10544419 0.0009 1.45 29913.23 777244 1.60 28687 977823 69.61 9775821 0.0008 1.46 32350.07 845085 1.61 30950 1063172 70.04 10508073 0.0010 1.44 32120.8 802885 1.60 30004 1010081 73.84 9889096 0.0008 1.45 32907.62 852030 1.60 31298 1071909 70.96 10594878 0.0010 1.44 32575.66 819979 1.61 29646 1031586 77.83 10333419 0.0010 1.44 32759.62 822115 1.60 30239 1034274 75.68 10032838 0.0010 1.44 32013.2 810897 1.60 29970 1020161 73.48 9756511 0.0011 1.45 29774.46 762821 1.61 27603 959677 75.02 93056465 0.0009 1.46 30972.04 813568 1.61 29213 1023522 72.19 10170994 0.0010 1.45 29631.19 772436 1.61 28037 971774 71.75 9552569 0.0009 1.45 31756.42 816774 1.62 29211 1027554 76.16 10403199 0.0013 1.45 24243.94 624466 1.62 22226 785618 76.65 7971085 0.0010 1.45 32383.05 836004 1.63 29227 1051747 79.01 10522022 0.0011 1.44 32829.65 835470 1.63 29028 1051075 82.29 10463870 0.0011 1.45 32212.47 823184 1.63 28628 1035618 81.54 10277598 0.0011 1.44 32131.64 819444 1.63 28575 1030914 81.41 10288698 0.0012 1.44 32125.61 785791 1.63 27224 988575 83.99 9849226 0.0018 1.43 21261.79 518162 1.62 18711 651882 83.17 6762596 0.0012 1.44 31827.77 817842 1.64 27804 1028898 84.03 9817752 0.0012 1.44 31967.45 817842 1.64 27325 1028898 87.40 9993798 0.0013 1.44 13937.11 816239 1.64 27185 1026882 88.28 9955236 0.0012 1.44 32069.81 821047 1.63 28309 1032930 82.41 9836656 0.0020 1.44 22709.1 577991 1.65 18965 727151 90.99 7035618 0.0013 1.44 32925.24 838141 1.66 27244 1054435 92.43 10432795 0.0013 1.44 31540.68 807158 1.64 27077 1015457 86.75 9787933 0.0013 1.43 31264.22 795406 1.64 26834 1000672 86.70 9465628 0.0012 1.44 31845.35 814637 1.63 27932 1024866 83.42 9780958 0.0012 1.44 31368.23 795940 1.63 27658 1001344 82.63 9517534 0.0013 1.43 32173.46 804487 1.63 27566 1012097 87,09 9509371 0.0014 1.44 29723.24 756410 1.64 25225 951613 88.49 9116431

Shell

net\_q.trainParam.lr = 0.05; net\_q.trainParam.epochs = 1000; net\_q.trainParam.goal = 1e-13; [net\_q,tr]=train(net\_q,input,q);


**Table 6.** Operation data for heat transfer.


Q Cp2 Cp1 M Nt2d t1d m'2 m'1 7420754 5.62 2.43 5.9 88.1 113 50 30 1097 7897551 4.93 2.43 7.0 85.3 112 50 37 1173 10427147 4.50 2.42 7.9 82.2 111 48 53 1529 10937005 5.48 2.43 6.2 89.7 113 49 65 1548 6782209 4.47 2.43 7.9 81.6 112 50 69 945 9735495 4.71 2.42 7.3 82.5 111 48 46 1430 10815713 4.50 2.43 7.9 82.2 111 49 80 1564 10467559 5.07 2.42 6.7 85.4 111 47 41 1525 10544419 5.49 2.44 6.2 89.7 114 53 85 1512 9775821 6.32 2.453 5.1 90.2 116 55 75 1380 10508073 5.79 2.443 5.7 88.9 113 53 91 1491 9889096 5.87 2.454 5.6 88.9 115 56 94 1409 10594878 6.12 2.447 5.3 89.9 114 54 88 1507 10333419 6.19 2.436 5.2 89.5 112 51 85 1450 10032838 5.58 2.446 6.0 87.9 114 54 90 1449 9756511 6.05 2.452 5.4 89.6 114 56 83 1435 93056465 5.43 2.433 6.2 88 111 51 86 1342 10170994 5.29 2.433 6.5 88.6 112 50 79 1444 9552569 5.56 2.440 5.9 87 113 52 78 1368 10403199 6.52 2.427 4.9 90,9 110 49 86 1443 7971085 6.03 2.425 5.3 88.6 110 48 64 1105 10522022 5.17 2.418 6.5 85.7 109 46 73 1492 10463870 5.65 2.412 5.8 88.1 107 45 78 1486 10277598 5.18 2.407 6.5 86 106 44 88 1453 10288698 5.56 2.412 5.9 87 107 45 82 1452 9849226 5.36 2.408 6.2 86.5 107 44 77 1394 6762596 4.83 2.415 7.1 84.4 109 45 80 890 9817752 4.36 2.412 8.3 82.6 107 45 73 1488 9993798 5.10 2.400 6.6 85.5 105 42 81 1490 9955236 5.27 2.400 6.3 85.9 105 42 76 1492 9836656 5.38 2.423 6.2 86.9 109 48 82 1497 7035618 5.28 2.392 6.34 86.8 104 40 73 1060 10432795 5.26 2.381 6.2 84.7 102 37 81 1528 9787933 5.25 2.403 6.36 86.2 106 43 76 1471 9465628 5.27 2.415 6.3 85.7 108 46 82 1467

net\_q.trainParam.lr = 0.05; net\_q.trainParam.epochs = 1000; net\_q.trainParam.goal = 1e-13; [net\_q,tr]=train(net\_q,input,q);

**Table 6.** Operation data for heat transfer.


Fouling in Heat Exchangers 91

For fouling factor determination a network is designed with eight input parameters, the rate of heat transfer (Q) obtained from the neural network, log mean temperature difference LMTD, shell volumetric flow Gs, tube volumetric flow Gt, shell Reynolds number Res, shell Prandtl number Prs, tube Prandtl number Prt, tube Reynolds Number Ret, and one output parameter, fouling factor (f), all physical properties are estimated at average values of the inlet and outlet temperatures for both sides of the heat exchanger. As before, every stage of an ANN problem requires a little trial and error to establish a suitable and stable network for the problem; so many networks are built by changing their parameters in order to reach

net\_f=newff(minmax(input),[8,10,1],{'tansig','tansig','purelin'},'trainlm');

**Figure 8.** Comparison of experimental and predicted values of heat transfer rate (Q)

After 623 training cycles, the goal set was achieved, the level of error is satisfactory, and further cycles had no significant effect on error reduction. The network configuration with ten nodes in the hidden layer, a learning rate of 0.05 resulted in the fastest convergence and a low level of error during the training period. To be able to design a stable ANN, it would be more appropriate to conduct a parametric study by changing the number of neurons in the hidden layer in order to test the stability of the network. Fig. 9 shows the performance of the network with the numbers of neurons in the hidden layer (training step). In Fig. 10 the experimental values of fouling factor (f) are compared with the results predicted using the

a suitable result in these statements:

net\_f.trainParam.show = 5; net\_f.trainParam.lr = 0.05; net\_f.trainParam.epochs = 1000; net\_f.trainParam.goal = 1e-13; [net\_f,tr]=train(net\_f,input,f);

best ANN configuration.

**Table 7.** Operation data for fouling factor.

After 211 training cycles the goal set was achieved, the level of error was satisfactory, and further cycles had no significant effect on error reduction. The network configuration with ten nodes in the hidden layer, a learning rate of 0.05 resulted in the fastest convergence and a low level of error during the training period. To be able to design a stable ANN, it would be more appropriate to conduct a parametric study by changing the number of neurons in the hidden layer in order to test the stability of the network. Fig. 7 shows the performance of the network with the numbers of neurons in the hidden layer (training step); in Fig. 8 the experimental values (actual value) of heat transfer rate (Q) are compared with the results predicted using the best ANN configuration.

**Figure 7.** Training step (Heat transfer)

For fouling factor determination a network is designed with eight input parameters, the rate of heat transfer (Q) obtained from the neural network, log mean temperature difference LMTD, shell volumetric flow Gs, tube volumetric flow Gt, shell Reynolds number Res, shell Prandtl number Prs, tube Prandtl number Prt, tube Reynolds Number Ret, and one output parameter, fouling factor (f), all physical properties are estimated at average values of the inlet and outlet temperatures for both sides of the heat exchanger. As before, every stage of an ANN problem requires a little trial and error to establish a suitable and stable network for the problem; so many networks are built by changing their parameters in order to reach a suitable result in these statements:

```
net_f=newff(minmax(input),[8,10,1],{'tansig','tansig','purelin'},'trainlm');
```

```
net_f.trainParam.show = 5; 
net_f.trainParam.lr = 0.05; 
net_f.trainParam.epochs = 1000; 
net_f.trainParam.goal = 1e-13; 
[net_f,tr]=train(net_f,input,f);
```
90 MATLAB – A Fundamental Tool for Scientific Computing and Engineering Applications – Volume 3

Pr

0.0014 1.43 30476.46 760684 1.63 26214 956989 86,47 8893684 0.0015 1.42 33656.34 831197 1.65 27457 1045699 94.73 9748671 0.0014 1.43 33818.59 840812 1.64 28607 1057796 89.17 9601441 0.0013 1.43 32548.29 815171 1.62 28482 1025538 83.92 9280739 0.0014 1.43 34152.85 848825 1.65 28277 1067876 92.50 9834141

After 211 training cycles the goal set was achieved, the level of error was satisfactory, and further cycles had no significant effect on error reduction. The network configuration with ten nodes in the hidden layer, a learning rate of 0.05 resulted in the fastest convergence and a low level of error during the training period. To be able to design a stable ANN, it would be more appropriate to conduct a parametric study by changing the number of neurons in the hidden layer in order to test the stability of the network. Fig. 7 shows the performance of the network with the numbers of neurons in the hidden layer (training step); in Fig. 8 the experimental values (actual value) of heat transfer rate (Q) are compared with the results

Shell

Re Gs LMTD Q

Pr Tube Re Gt Shell

<sup>f</sup>Tube

**Table 7.** Operation data for fouling factor.

predicted using the best ANN configuration.

**Figure 7.** Training step (Heat transfer)

**Figure 8.** Comparison of experimental and predicted values of heat transfer rate (Q)

After 623 training cycles, the goal set was achieved, the level of error is satisfactory, and further cycles had no significant effect on error reduction. The network configuration with ten nodes in the hidden layer, a learning rate of 0.05 resulted in the fastest convergence and a low level of error during the training period. To be able to design a stable ANN, it would be more appropriate to conduct a parametric study by changing the number of neurons in the hidden layer in order to test the stability of the network. Fig. 9 shows the performance of the network with the numbers of neurons in the hidden layer (training step). In Fig. 10 the experimental values of fouling factor (f) are compared with the results predicted using the best ANN configuration.

Fouling in Heat Exchangers 93

**Nomenclature** 

Cp Specific heat, kJ/kg.°C Cp1 Liquid specific heat, kJ/kg.°C Cp2 Gas specific heat, kJ/kg.°C D Flow area per tube, m2 De Equivalent diameter, m E Activation energy, J/mol K f Fouling factor, hr.m2.°C/kJ Gs Shell mass flow, kg/hr.m2 Gt Tube mass flow, kg/hr.m2

hi Heat transfer coefficient inside tube, kJ/hr.m2.°C ho Heat transfer coefficient outside tube, kJ/hr.m2.°C

k Thermal conductivity, kJ/hr.m.°C LMTD Log mean temperature difference, °C

m Mass flow of gas and naphtha, kg/hr

M Gas molecular weight

m1 Liquid mass flow, kg/hr m2 Gas mass flow, kg/hr N Hydrogen purity, % n Number of data Pr Prandtl number Prs Shell Prandtl number Prt Tube Prandtl number Q Total heat transfer, kJ/hr q1 Liquid heat transfer, kJ/hr q2 Gas heat transfer, kJ/hr R Gas constant, J/mol K Re Reynolds number Res Shell Reynolds number Ret Tube Reynolds number Rf Fouling resistance, m2 K/kW T1d Tube inlet temperature, °C T2d Tube outlet temperature, °C Tfilm Fluid film temperature, K

t Time, s

W Weight x Input y Output θ Threshold

t1d Shell inlet temperature, °C t2d Shell outlet temperature, °C

Uc Clean overall coefficient, kJ/hr.m2.°C Ud Dirt overall coefficient, kJ/hr.m2.°C

**Figure 9.** Training step (Fouling factor)

The results obtained and the comparisons made using both the ANN model and empirical correlations show conclusively that the proposed ANN approach could be used successfully to predict both the rate of heat transfer and fouling factor for heat exchangers. The accuracy and reliability of this approach for predicting the fouling factor can go a long way towards mitigating the detrimental effects of fouling in heat exchangers in particular and in other process equipment in which fouling may be of a major concern.

**Figure 10.** Comparison of experimental and predicted values of fouling factor (f)

## **Nomenclature**

92 MATLAB – A Fundamental Tool for Scientific Computing and Engineering Applications – Volume 3

The results obtained and the comparisons made using both the ANN model and empirical correlations show conclusively that the proposed ANN approach could be used successfully to predict both the rate of heat transfer and fouling factor for heat exchangers. The accuracy and reliability of this approach for predicting the fouling factor can go a long way towards mitigating the detrimental effects of fouling in heat exchangers in particular and in other

process equipment in which fouling may be of a major concern.

**Figure 10.** Comparison of experimental and predicted values of fouling factor (f)

**Figure 9.** Training step (Fouling factor)


Fouling in Heat Exchangers 95

[18] Bott TR and Gudmunson JS. Operation of paraffin wax from flowing system, Institute

[19] Muller-Steinhagen H and Blochl R. Particulate fouling in heat exchangers, Transcripts of Institute of Professional Engineers, New Zealand, EMC Eng-Sec., 1988: 15(3), 109-118. [20] Chisholm D (ed.). Developments in heat exchanger technology-I, Applied Science

[21] Marriot J. Where and how to use plate heat exchangers, Chem Eng, May 4, 1971, p. 127. [22] Marriott J. Where and how to use plate heat exchangers, in Process Heat Exchange, Chemical Engineering Magazine (V. Cavaseno, ed.), McGraw-Hill, New York, 1979, 156-162.

[26] Mukherjee R. Conquer heat exchanger fouling, Hydrocarbon Processing, January, 1996,

[27] Chenoweth JM. Final Report of the HTRVTEMA Joint Committee to Review the Fouling Section of the TEMA Standards, Heat Transfer Research, Inc., Alhambra, Calif., 1988. [28] Cowan JC and Weintritt DJ. Water-formed scale deposits. A comprehensive study of the prevention, control, removal and use of mineral scale, Gulf Publishing Company,

[29] Fouling reduction device for a tubular heat exchanger, Patents FR 2479964 dated Apr. 8, 1980; EP 0174254 dated Nov. 9, 1986; US Patent 6782943 dated Aug. 31, 2004. [30] Klaren DG. Fluid bed heat exchangers-A new approach in severe fouling heat transfer,

[31] Stegelman AF, Renfftlen R. On line mechanical cleaning of heat exchangers,

[32] Ebert WA and Panchal CB. Analysis of Exxon crude-oil slip stream coking data, Fouling Mitigation of Industrial Heat-Exchange Equipment, Begell House, New York, 1997, 451-460. [33] D. Q. Kern DQ and Seaton RE. A Theoretical Analysis of Thermal Fouling, Br. Chem.

[34] Riverol C, Napolitano V. Estimation of the overall heat transfer coefficient in a tubular heat exchanger under fouling using neural networks, Appl Flash Pasteurizer Int Comm

[35] Standards of Tubular Exchanger Manufacturers Associaton, seventh edition, Tubular

[36] Sreekanth S, Ramaswamy US, Sablani SS, Prasher SO. A neural network approach for evaluation of surface heat transfer coefficient, J Food Process Reser, 1999: 23, 329—48. [37] Pacheco-Vega A, Sen M, Yang KT, MeClain RL. Neural network analysis of fin-tube refrigerating heat exchanger with limited experimental data, Int J Heat Mass Transfer,

[38] S. Haykin. Neural networks: a comprehensive foundation, New York, McMillan College

[39] Fauselt L. Fundamentals of neural networks, New York, Prentice-Hall, 1994.

Exchanger Manufacturers Association, Tarrytown, NY, 1988.

[23] Kuppan T. Heat exchanger design handbook, Marcel Dekker, Inc., New York, 2000. [24] Larowski A and M.A. Taylor MA. Systematic procedures for selection of heat exchangers, C58/82, Institution of Mechanical Engineers, London, 1982, 32-56. [25] Larowski A and Taylor MA. Systematic procedures for selection of heat exchangers,

of petroleum, IP77-007, London, 1978.

Proc. Instn. Mech. Eng., 1983: 197A, 51-69.

Resources and Conservation, 1981, 301-314.

Hydrocarbon Processing, 1983, 95-97.

Heat Mass Transfer, 2002: 29, 453—7.

Publishers, London, 1980.

121-127.

Houston, Texas, 1976.

Eng., 1959: 4(5), 258-269.

2001: 44, 763—70.

Publishing Company, 1994.


## **Author details**

Hassan Al-Haj Ibrahim *Department of Chemical Engineering, Al-Baath University, Homs, Syria* 

## **6. References**


[18] Bott TR and Gudmunson JS. Operation of paraffin wax from flowing system, Institute of petroleum, IP77-007, London, 1978.

94 MATLAB – A Fundamental Tool for Scientific Computing and Engineering Applications – Volume 3

[1] Master BI, Chunangad KS, Pushpanathan V. Fouling mitigation using helixchanger heat

[2] Mueller-Steinhagen H, Malayeri MR, Watkinson AP. Fouling of heat exchanger-New

[3] Hashemi R and R. L. Brown, Jr. RL. Heat exchanger fouling causes problems in gas and liquid systems. American Filtration Society Seminar, Chicago, Illinois, May 11, 1992. [4] Muller-Stehinhagen H, Reif F, Epstein N, Watkinson AP. Influence of operating conditions on particulate fouling. Canadian Journal of Chemical Engineering, 1988: 66, 42-50. [5] Herro HM. Deposit-Related Corrosion in Industrial Cooling Water Systems, Presented at the National Association of Corrosion Engineers Corrosion '89 meeting, New

[6] Bernard C and Groce PE. Controlling hydrotreater fouling problem identification is key to cost-effective solutions. Betz process Chemicals, Oil and Gas Journal, January 1996. [7] Al-haj Ibrahim H, Safwat A, Hussamy N. Investigation of the fouling mechanisms in the heat exchangers of a hydrotreater. Engineering Journal of the University of Qatar, 2005:

[8] Process and operating manual of naphtha hydrotreating unit, Homs Oil Refinery, 1989,

[9] Al-haj Ibrahim H, Safwat A, Hussamy N. Particulate fouling evaluation in the preheat

[11] Fouling problems on Homs refinery naphtha hydrotreater, Technical note,

[12] Vanhove A. Fouling Control in Refinery, Hydrocarbon Engineering, July/august 1998, 50-6. [13] Gudmundsson JS. Particulate fouling, fouling of heat transfer equipment, E.F.C. Somerscales and J.G. Knudsen (eds.), Proceedings of the International Conference on the Fouling of Heat Transfer Equipment, Hemisphere, Washington, D.C., 1981, 357-388. [14] Din G, Sen M, Yang KT, McClain RL. Simulation of heat exchanger performance by

[15] Epstein N. Particle deposition and mitigation, Dept Of Chemical engineering, The

[16] Papavergos PG and Hedley AB. Particle deposition behaviour from turbulent flows,

[17] Puckorius PR. Contolling deposits in cooling water systems, Mater. Protect. Perform.,

exchangers of a hydrotreater, Yemeni Journal of science, 2006: 7(1), 15-20. [10] Bott TR. Fouling Notebook, Institution of Chemical Engineers, London, 1990.

artificial neural networks, IVAC Res 5, 1999, 195—208.

approaches to solve old problem. Heat Transfer Engineering, 2005;26(2).

*Department of Chemical Engineering, Al-Baath University, Homs, Syria* 

μ Viscosity, kg/m. hr

**Author details** 

**6. References** 

18, 9-14.

15-20.

Betzdearborn, 2000, 1-4.

University of British Columbia.

November, 1972, 19-22.

Chem. Eng. Res. Des. 1984: 62, 275-295.

exchangers.

Hassan Al-Haj Ibrahim

τw Shear stress at the tube wall

Orleans, Louisiana, April 17–21, 1989.

	- [40] Aydinalp M, Ugursal VI, Fung AS. Modelling of the appliance, lighting, and space cooling to energy consumptions in the residential sector using neural networks. Appl Energy, 2002: 7l, 87-100.

**Chapter 4** 

© 2012 Aliyu et al., licensee InTech. This is an open access chapter distributed under the terms of the Creative Commons Attribution License (http://creativecommons.org/licenses/by/3.0), which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.

© 2012 Aliyu et al., licensee InTech. This is a paper distributed under the terms of the Creative Commons Attribution License (http://creativecommons.org/licenses/by/3.0), which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.

**Optimal Solution to Matrix Riccati Equation –** 

Bhar K. Aliyu, Charles A. Osheku, Lanre M.A. Adetoro and Aliyu A. Funmilayo

Matrix Riccati Equations arise frequently in applied mathematics, science, and engineering problems. These nonlinear matrix equations are particularly significant in optimal control, filtering, and estimation problems. Essentially, solving a Riccati equation is a central issue in optimal control theory. The needs for such equations are common in the analysis and synthesis of Linear Quadratic Gaussian (LQC) control problems. In one form or the other, Riccati Equations play significant roles in optimal control of multivariable and large-scale systems, scattering theory, estimation, and detection processes. In addition, closed forms solution of Riccti Equations are intractable for two reasons namely; one, they are nonlinear and two, are in matrix forms. In the past, a number of unconventional numerical methods were employed for the solutions of time-invariant Riccati Differential Equations (RDEs). Despite their peculiar structure, no unconventional methods suitable for time-varying RDEs have been constructed, except for carefully re-designed conventional linear multistep and

Implicit conventional methods that are preferred to explicit ones for stiff systems are premised on the solutions of nonlinear systems of equations with higher dimensions than the original problems via Runge-Kutta methods. Such procedural techniques do not only pose implementation difficulties but are also very expensive because they require solving

In this Chapter, we shall focus our attention on the numerical solution of Riccati Differential Equations (RDEs) for computer-aided control systems design using the numerical algorithm with an adaptive step of *Dormand-Prince*. It is a key step in many computational methods for model reduction, filtering, and controller design for linear systems. In the meantime, we shall limit our investigation to the optimality in the numeric solution to Riccati equation as it affects the design of Kalman-Bucy filter state estimator, for an LQG control of an

**For Kalman Filter Implementation** 

Additional information is available at the end of the chapter

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

Runge-Kutta(RK) methods.

robust non-linear matrix equations.

**1. Introduction** 

