CHEMICAL ENGINEERING TRANSACTIONS VOL. 70, 2018 A publication of The Italian Association of Chemical Engineering Online at www.aidic.it/cet Guest Editors: Timothy G. Walmsley, Petar S. Varbanov, Rongxin Su, Jiří J. Klemeš Copyright © 2018, AIDIC Servizi S.r.l. ISBN 978-88-95608-67-9; ISSN 2283-9216 Mathematical Modelling of the Thermal and Hydraulic Behaviour of Plate Heat Exchanger in the Fouling Conditions Petro O. Kapustenkoa,*, Oleksiy V. Demirskiya, Vladimir I. Tovazhnyanskyia, Jiří J. Klemešb, Oleksandr I. Matsegorac, Pavlo Y. Arsenyevc, Olga P. Arsenyevad aNational Technical University “Kharkiv Polytechnic Institute”, Dep. ITPA, 21 Frunze st., 61002 Kharkiv, Ukraine bSustainable Process Integration Laboratory – SPIL, NETME Centre, Faculty of Mechanical Engineering, Brno University of Technology - VUT Brno, Technická 2896/2, 616 69 Brno, Czech Republic cAO Spivdruzhnist-T LLC, 2 Kaplunovskiy per., 19, 61002 Kharkiv, Ukraine dUniversity of Paderborn, Chair of Fluid Process Engineering, Paderborn, Germany sodrut@gmail.com The mathematical model of Plate Heat Exchanger (PHE) subjected to fouling is proposed. It is represented by the system of ordinary differential equations. The model is accounting for the distribution of process parameters along the PHE channel that allows predicting fouling development in time at different locations along the channel length. The development of the fouling deposit is accounted with the fouling model presented by the equation in dimensionless form. The relative influence of different terms is characterized by empirical coefficients which can be identified with the data of monitoring the PHE thermal and hydraulic performance. The model allows also the prediction of pressure drop variation in PHE with the development of fouling deposition layer and respective reduction in channels cross-section area. The application of the model and its accuracy is demonstrated with two case studies considering the monitoring of PHEs thermal and hydraulic performance in the industry at sugar factory and in District Heating system. 1. Introduction One of the important prerequisites for the sustainable development of mankind is the efficient use of resources, among which energy plays the leading part. It is generated nowadays mostly with the combustion of fossil fuels, not only exhausting their limited reserves but also emitting harmful combustion products into the environment. For efficient use of energy, the increasing of heat recuperation is of primary importance. It is possible with the process integration technique supplemented by the use of efficient heat transfer equipment (Klemeš et al., 2013). Plate heat exchanger (PHE) is one of the modern efficient types of compact heat exchangers (Klemeš et al., 2015). The compact construction consisting of corrugated plates, stamped from a thin metal sheet, and intensified heat transfer in their channels allows a significant reduction of heat transfer area for the same duty compare to conventional shell and tube heat exchangers. The channels of PHEs are of a small hydraulic diameter and have an intricate geometry which promotes turbulence enhancing heat transfer coefficients. In such conditions, the correct accounting for fouling became of primary importance. The same amount of fouling deposit taken according to recommendations for shell and tube heat exchanger can lead to excessive heat transfer area and lower velocities in channels. It can be followed by a higher reduction in heat transfer performance and more severe blockage of narrow channels. The correct accounting of fouling in PHE design requires the use of the sufficiently detailed mathematical models with the application of accurate enough models of fouling formation. As it is presented in a paper by Kapustenko et al. (2017), for the modelling of PHEs thermal performance, it is possible to account for local parameters of fouling development in PHE channels with one-dimensional model. However, besides PHE thermal behaviour the correct prediction of pressure drops and its time development is very important for estimation of its operational performance and maintenance schedule. The thermo-hydraulic fouling modelling DOI: 10.3303/CET1870019 Please cite this article as: Kapustenko P.O., Demirskiy O.V., Tovazhnyanskyi V.I., Klemes J.J., Matsegora O.I., Arsenyev P.Y., Arsenyeva O.P., 2018, Mathematical modelling of the thermal and hydraulic behaviour of plate heat exchanger in the fouling conditions , Chemical Engineering Transactions, 70, 109-114 DOI:10.3303/CET1870019 109 for crude oil refinery heat exchanger was proposed by Yeap et al. (2004). There is presented an analysis of some fouling models of reaction and transport type based on the equation proposed for fouling deposition term by Epstein, that according to the paper of Epstein (2011) is written as follows:  s 1 E R T1 1 d f b 1 m 2 w C k k k e               (1) where λf is thermal conductivity of fouling deposit, J/(m.K); μ is dynamic viscosity of liquid, Pa.s; τw is wall shear stress, Pa; E is activation energy, J/mol; R = 8.314 J/mol is universal gas constant; Ts is temperature of the fouling layer surface, K; km is mass transfer coefficient, m/s; k1 and k2 are dimensional constants; Cb is concentration of fouling precursors, mol/m3 The fouling models based on Eq(1) are demonstrating good results in correlating the experimental data on fouling at intensified heat transfer surfaces, as shown in the paper of Yang and Crittenden (2012) for crude oil fouling in tubes with inserts and in the paper of Arsenyeva et al. (2013a) for water fouling in PHE channels. However, the presentation of the model in these papers in dimensional form requires obtaining empirical coefficients of complicated dimensional units and jeopardises the generalisation capacity of the model. The aim of the present paper is to modernize fouling model to dimensionless form and to develop on its base the mathematical model of fouling in PHEs accounting for local parameters of their thermal and hydraulic performance. 2. Mathematical modelling The core issue of PHE modelling under fouling formation conditions is the selection of reliable and sufficiently accurate fouling model. As it was shown in the paper by Kapustenko et al. (2017) the good results on thermal modelling for water fouling in PHEs can be achieved with Kern and Seaton approach where fouling rate is described as a difference between fouling deposition term φd and fouling removal term φrm: f d rm /       (2) In that paper the deposition term was presented in a form derived from Eq(1). In this study, the Eq(1) is presented in dimensionless form by the introduction of the set of dimensionless complexes with the use of following assumptions. 1. The analogy between heat and mass transfer holds the place and it can be written for the diffusional Nusselt number:   1/3 D2 m e 2 D2 2 Nu k d / D Nu Pr Pr    (3) where Nu2 is Nusselt number for heat transfer in cooling media; Pr2 = cp2.μ2/λ2 is Prandtl number; PrD =D.ρ2/μ2 is diffusional Prandtl number; ρ2 is liquid density, kg/m3; λ2 is thermal conductivity of liquid, J/(m.K); cp2 is specific heat capacity of liquid, J/(kg.K); de is channel equivalent diameter, m; D is coefficient of diffusion, m/s. 2. The link between the coefficient of diffusion and dynamic viscosity for specific substances can be determined with modification of Stocks-Einstein equation (Einstein, 1906) in the following form:  s B 2 mD T k / r     (4) where kB = 1.38048.10-23 J/K is Boltzmann constant; Ts is temperature of the deposit surface, K; μ is dynamic viscosity, Pa.s; rm is Van der Vaals radius of molecule, m, that is introduced as scale for molecule radius; χ is parameter depending on the solution substances nature which is accounting for radius of solute molecule and deviations from Stocks-Einstein equation for specific solution content. Such character of the relation between D, μ and T was confirmed experimentally by Wilke and Chang (1955) for different solutions in their study establishing the dependence of proportionality factor from solute properties. The same result is also following from more recent works, as e.g. Karunanithi and Bogeshwaran (2016). It can be concluded that for the specific solute the parameter χ can be regarded as not changing with concentration and temperature of the solution. As the radius of the molecule for different species is rather uncertain it is taken the radius for water molecule rm = 1.36.10-10 m according to Zhang and Xu (1995). It is introduced as scaling factor and the differences in rm are accounted in parameter χ. Using Eq(3) and Eq(4), after rearranging the component terms, Eq(1) can be written in the dimensionless form as follows: 110    1 2 1 3 3 d e 2 2 D D 2 2 R R s d / c K Pr / Nu c K exp E / R T               (5) Here KD and KR are dimensionless complexes of variables which value can change with temperature and flow characteristics. These complexes are expressed as follows:  2D 2 m s 2 BK r / T k ;      R w 2 eK / d g .    (6) Another two dimensionless complexes cD and cR are:   2 3 D 1 b f c k / C ;       R 2 b fc k g / C .   (7) These complexes are dependent only from constants k1 and k2 in Eq(1) and variables which are not changing for the same solution. For experiments with the same media and fouling substance cD and cR can be regarded as constants and determined by correlating the tests data. The removal term in dimensionless form can be expressed using dimensional parameter B: *2 rm e 2 2 f w e 2 2 rm f 2 e d / B d / c Re Pr / d                  (8) Here Re* is Reynolds number calculated with the velocity of wall shear stress: * w 2 e 2 2 Re / d / ;      rm 2 p 2c B / c  (9) Counting that for the same liquid the variation of cp2 and λ2 are rather small, the parameter crm can be regarded as constant when correlating the experimental data on fouling for specific media. The mathematical model of PHE thermal behaviour under conditions of fouling was presented in a paper by Kapustenko et al. (2017). However, there is used dimensional fouling model and dimensionless fouling model presented here require validation. It also requires an development to include calculation of the PHE hydraulic behaviour. The basic differential equations of thermal design for countercurrent flows are: 2 2 p2 T / q / (g c )    x (10) 1 2 2 p2 1 p1 T / ( T / ) g c / (g c )       x x (11) where T1 and T2 are temperatures of hot and cold streams, K; g1 and g2 are the mass flow rates of hot and cold streams through one channel, kg/s; cp1 and cp2 are specific heat capacities of hot and cold streams, J/(kg.K); Π is the channel perimeter, m; x is the distance from cold stream entrance, m; q is the specific heat flux, W/m2.   -1 1 2 f w w 1 2 q = 1 / h +1 / h + R + / × (T - T )  (12) Here δw is the thickness of the plate metal, m; λw is the heat conductivity of the plate metal, W/(m.K); h1 and h2 are film heat transfer coefficients for hot and cold streams, respectively, W/(m2.K). The temperature at the outer surface of the deposited layer is:   -1 s 1 2 f w w 2 1 2 2 T = 1 / h + 1 / h + R + / × h × (T - T ) + T    (13) The coefficients h1 and h2 are calculated according to correlations presented by Kapustenko et al. (2011) for pressure drop and heat transfer at the main corrugated field depending on channel geometry and fluid thermo- physical properties. In general form: h h (W ,T ,T , , ,d )  j j j j s ej (14) Here γ=2.b/S is the corrugation doubled height to pitch ratio; β is corrugation angle, degrees. With fouling formation, the free cross-section area of the channel (fch, m2) is becoming smaller inducing the increase of flow velocity that can be determined as follows: 111 2 2 2 g W ( )      ch f f (15) As it is discussed in a paper by Yeap et al. (2004), the deposited fouling layer is leading to increasing the pressure loss in the channel not only causing higher flow velocity but also creating a roughness on the flow boundary. For the majority of streams in heat exchangers, the thermo-physical properties are not much dependent on the change of pressure inside channels (especially for liquids). For this reason, the thermal part of the mathematical model can be solved without considering pressure loss in channel thus enabling to receive information about fouling layer development prior to consequent calculation of pressure losses. The pressure drop in the PHE can be considered as a sum of pressure drop at the main corrugated field ΔPmf, pressure drops at inlet and outlet distribution zones ΔPzin and ΔPzout, the pressure drop in port and collectors ΔPpc (Arsenyeva et al., 2013b). The total pressure drop with expressions for corresponding constituent parts is as follows: 22 2 2 L 2 22 2 2 2 2 2 2 2 0 WW W W P 1.3 2 d 2 2 2                     p pin out DZin DZout e dx (16) where W2in, W2out and W2p are velocities at channel inlet, outlet and ports of PHE, m/s; ζ2 is the friction factor in PHE channel determined depending on channel geometry by Eq (17) proposed by Arsenyeva et al. (2011) with the term A accounting for roughness created by fouling: 1 3 12 12 16 2 2 12 2 37530 1 8 A Re Re                           p p ; 16 1 0.9 7 3 A 4 ln 5 2.7 Re d                      f e p p p (17)             2 2 2.63 0.01 1 exp 0.157 ; 2 / 3; 3 exp / 180 ; 5 1 ; 10 4 0.061 0.69 / 180 1 1 0.9                                   p p p p p tg (18) The coefficient of local hydraulic resistance in inlet distribution zone is determined according to the paper by Arsenyeva et al. (2013), namely ζ2 = 38. For outlet distribution zone to value 38 the correction for fouling roughness is taken as for friction factor at the end of the main corrugated field and also the velocity W2out is calculated for cross-section area reduced by fouling. The Eqs(2)-(18) with correlations for temperature dependences of thermal and physical properties of streams and geometrical relations for PHE represent the system of ordinary differential equations with nonlinear right parts. The numerical solution of this system with finite difference method is implemented for PC using Mathcad software. To check the validity of the developed model in the following section two case studies are considered. 3. Model validation and discussion of the results For illustration of model application and its validation in this section are presented the comparisons of calculations results with data of tests in the industry. 3.1 Case study 1. PHE at the sugar factory The PHE for heating the juice directed to evaporation by the heat from condensate after evaporation effect presented by Demirskiy et al. (2018) is considered. The PHE is of the M15M type produced by AlfaLaval. It consists of 151 plates forming 75 channels for thin juice. The heat transfer area of one plate is 0.62 m 2, the length of corrugated field is 1.38 m, cross section area of channel equals 0.00176 m2, the height of corrugations is 4 mm, the corrugations inclination angle β = 35º, the doubled height to pitch ratio γ = 0.56. The results of monitoring operating parameters during 13 d of operation are presented in Table 1. By these data the temperatures of heated juice and cooled condensate after PHE are calculated and compared to experimental values in Table 1. 112 Table 1: The results of monitoring parameters and calculation by the model of PHE heating thin juice Time θ, h 144 216 264 312 The flowrate of thin juice G2, m3/h 260 270 277 265 Inlet temperature of thin juice t21, ºC 101 100.5 102 101.7 The outlet temperature of thin juice: experimental t22exp, ºC calculated t22clc, ºC 105 104.9 106 105.9 107 106.8 106 105.9 Inlet temperature of condensate t11, ºC 123.5 123.5 123.5 123.5 Outlet temperature of condensate: experimental t12exp, ºC calculated t12clc, ºC 102.8 102.6 104.8 104.4 106.1 105.8 104.8 104.6 The calculated results are obtained for the following values of dimensionless empirical parameters in fouling model, which were determined by least squares method on experimental data: cD = 2.29.106; cR = 0.126; crm = 0.45.10-15. The discrepancies between the calculated and experimental temperatures are not bigger than 0.3 ºC that confirms the validity of the thermal modelling. After rearranging the piping at PHE to make two connections for thin juice, the hydraulic parameters of operation were monitored for 93 d during the operation the next year. The comparison of experimental pressure drops and calculated is presented in Figure 1a. The calculations are made for the average value of thin juice flow rate G2 = 290 m3/h, which was actually changed ± 7 % during the operation. Counting for uncertainties in industrial operation conditions the results of calculations and model accuracy can be regarded as acceptable for predicting the hydraulic behaviour of PHE under the fouling conditions of thin juice heating at the sugar factory. (a) (b) Figure 1: (a) Pressure drop for Case study 1; where solid line -calculated; dots - tests. (b) Overall heat transfer coefficient for Case study 2; where solid dots – calculation; empty dots – tests 3.2 Case study 2. PHE at District Heating system Nowadays the District Heating systems are one of the biggest users of PHEs, which sake for efficiency and compactness practically completely overtaken new installations in this field from traditional shell and tube units. Demirskiy et al. (2018) have checked the validity of PHE thermal model employing a dimensional model of fouling by experimental data of Chernyshov (2002). These experiments performed with AlfaLaval M10B PHE for heating fresh water coming with temperature in the range 7.9 - 9.5 ºC to about 60 ºC by hot water from the boiler with temperature varied from 74 to 98 ºC. The inlet temperature of hot water was rising eventually to ensure the required temperature of the heated fresh water with fouling development. The empirical values of dimensionless parameters in presented here fouling model were determined by least squares method using experimental data for one set of tests with flow velocity 0.57 m/s. These values are: cD = 8.18.105; cR = 0.045; crm = 0.17.10-16. The comparison of the data for all experiments with overall heat transfer coefficients calculated by the model is presented in Figure 1b. The discrepancies between calculated and experimental results are not exceeding ±7 %. It confirms the thermal model validity and its ability to predict PHE fouling behaviour in an investigated range of flow velocities and temperatures. The absence of data on pressure losses does not permit to check the hydraulic part of the modelling in this case. 113 4. Conclusions The mathematical model of PHEs thermal and hydraulic performance under fouling deposition on heat transfer surface is presented. The integral part of it is fouling model in dimensionless form with empirical dimensionless coefficients. It simplifies the use of the model and adds for prospects of its generalisation that should be checked in further research with different fouling conditions. There is a good agreement of the modelling results with obtained data on temperature program and pressure loss in heat exchangers. The calculation by the model enables more accurate design of PHEs to improve their performance and to implement fouling mitigation effects. It allows significant increase of the time between cleaning of heat transfer surface. The developed software can be used for the optimal calculation of PHEs subjected to fouling with an accounting of their capital cost and operational expenses including cost of cleaning. Acknowledgements This research has been supported by the EU project “Sustainable Process Integration Laboratory – SPIL”, project No. CZ.02.1.01/0.0/0.0/15_003/0000456 funded by EU “CZ Operational Programme Research, Development and Education”, Priority 1: Strengthening capacity for quality research in a collaboration agreement with National Technical University “Kharkiv Polytechnical Institute” and AO Spivdruzhnist-T LLC. Olga Arsenyeva is grateful to the Alexander von Humboldt Foundation for the financial support. References Arsenyeva O.P., Crittenden B., Yang M., Kapustenko P.O., 2013a. Accounting for the thermal resistance of cooling water fouling in plate heat exchangers. Applied Thermal Engineering, 61(1), 53-59. Arsenyeva O., Kapustenko P., Tovazhnyanskyy L., Khavin G., 2013b. The influence of plate corrugations geometry on plate heat exchanger performance in specified process conditions. Energy, 57, 201-207. Arsenyeva O.P., Tovazhnyanskyy L.L., Kapustenko P.O., Khavin G.L., 2011. The Generalized Correlation for Friction Factor in Criss-cross Flow Channels of Plate Heat Exchangers, Chemical Engineering Transactions, 25, 399-404. Chernyshov D.V., 2002. Prognosis of scaling in plate water heaters to increase reliability of their work. Thesis for Candidate of Technical Sciences. Tula State University, Tula, Russian Federation. (In Russian). Demirskiy O.V., Kapustenko P.O., Arsenyeva O.P., Matsegora O.I., Pugach Y. A., 2018. Prediction of fouling tendency in PHE by data of on-site monitoring. Case study at sugar factory. Applied Thermal Engineering, 128, 1074-1081. Einstein A., 1906. Investigation on the Theory of the Brownian Movement. Annalen der Physik (4) 19, 371-381 (In German), English translation in Dover Publications Inc., New York, 1956. Epstein N., 2011. Comments on “Relate Crude Oil Fouling Research to Field Fouling Observations by Joshi et al.,”. In Proc. International Conference on Heat Exchanger Fouling and Cleaning - 2011, 62-64. Kapustenko P., Arsenyeva O., Dolgonosova O., 2011. The heat and momentum transfers relation in channels of plate heat exchangers, Chemical Engineering Transactions, 25, 357-362. Kapustenko P., Arsenyeva O., Matsegora O., Kusakov S , Tovazhnianskyi V., 2017. The mathematical modelling of fouling formation along PHE heat transfer surface, Chemical Engineering Transactions, 61, 247-252. Karunanithi B.,Bogeshwaran, K., 2016. Liquid Diffusion-Measurement and Correlation of Diffusion Coefficient in Acetic acid-Carbon tetra chloride System. International Journal of ChemTech Research. 9(08), 465-478. Klemeš J.J., Arsenyeva O., Kapustenko P., Tovazhnyanskyy L., 2015. Compact Heat Exchangers for Energy Transfer Intensification: Low Grade Heat and Fouling Mitigation. CRC Press, Boca Raton, FL, USA. Klemeš J.J., Varbanov P.S., Kapustenko P.O., 2013. New developments in heat integration and intensification, including total site, waste-to-energy, supply chains and fundamental concepts. Applied Thermal Engineering, 61, 1–6. Wilke C.R., Chang P., 1955. Correlation of diffusion coefficients in dilute solutions. AIChE Journal, 1(2), 264- 270. Yang M, Crittenden B., 2012. Fouling thresholds in bare tubes and tubes fitted with inserts, Applied Energy, 89, 67–73. Yeap B.L., Wilson D.I., Polley G.T., Pugh S.J., 2004. Mitigation of crude oil refinery heat exchanger fouling through retrofits based on thermo-hydraulic fouling models. Chemical Engineering Research and Design, 82(1), 53-71. Zhang Y., Xu Z., 1995. Atomic radii of noble gas elements in condensed phases. American Mineralogist, 80 (7-8), 670-675. 114