- Email: [email protected]
Estimation of thermophysical properties of fouling using inverse problem and its impact on heat transfer efficiency
Available online at www.sciencedirect.com
Solar Energy 83 (2009) 1619–1628 www.elsevier.com/locate/solener
Estimation of thermophysical properties of fouling using inverse problem and its impact on heat transfer efficiency A. Adili *, C. Kerkeni, S. Ben Nasralla Centre de Recherches et des Technologies de l’Energie, Laboratoire d’Energe´tique et des Proce´de´s Thermiques, Route touristique de Soliman, B.P. 95, Hammam-Lif 2050, Tunisia Received 11 August 2008; received in revised form 6 April 2009; accepted 24 May 2009 Available online 28 June 2009 Communicated by: Associate Editor Estrada-Gasca
Abstract At high temperature, the circulation of fluid in heat exchangers provides a tendency for fouling accumulation to take place on the internal surface of tubes. In brief, the deposits on heat exchanger tubes are caused by the presence of inorganic salts, of small quantities of organic materials and products of corrosion in the water. From thermophysical point of view, the deposited fouling has harmful effects on the heat exchanger efficiency. Indeed, it increases the thermal resistance which can raise the energy consumption. This study shows an experimental and a theoretical process of estimation of thermophysical properties of the fouling deposited on a section of a heat exchanger and its effects on the heat transfer efficiency. The estimation method is based on the Gauss–Newton algorithm that minimizes the ordinary least squares function comparing a measured temperature and a theoretical one. The temperature response is measured on the rear face of a bi-layer system composed of a section of a heat exchanger and the fouling deposited on during and after a finite width pulse heat flux on its front face. The theoretical temperature, that is a function of the unknown thermophysical properties of the bi-layer system, is calculated by the resolution of the one-dimensional linear inverse conduction problem, and by the use of the quadrupole formalism. The results of the estimation procedure show, on the one hand the efficiency and the stability of the optimization algorithm to estimate the thermophysical properties of the fouling. On the other hand they underline the necessity of the maintenance of fluid circulating tubes at high temperature. Ó 2009 Elsevier Ltd. All rights reserved. Keywords: Inverse problem; Thermophysical properties; Fouling; Thermal quadrupole formalism
1. Introduction Thermophysical characterization of materials is extremely important in the control current for today fabrication processes and in the prediction of structure failure due to thermal stresses. Accuracy in the estimation of the thermal properties can be improved if the experiments are designed carefully (Garcia, 1999). For these reasons, many researchers are interested in the determination of these parameters.
*
Corresponding author. E-mail address: [email protected] (A. Adili).
0038-092X/$ - see front matter Ó 2009 Elsevier Ltd. All rights reserved. doi:10.1016/j.solener.2009.05.014
In many industrial systems, in which there is a water circulation at high temperature, the problem of fouling deposited on the internal surfaces of tubes is found. The fouling formed has a considerable effect on the considered systems efficiency. Indeed, it can modify its thermophysical properties. For these reasons, the determination of the thermophysical properties of the deposited composite material has become very important to understand the performance of the considered system. Many researchers are interested in the determination of thermophysical properties of one or multilayer materials without application of their methodologies and results in industrial or current application. Albouchi et al. (2005)
1620
A. Adili et al. / Solar Energy 83 (2009) 1619–1628
Nomenclature thermal diffusivity (m2 s1) heat capacity (J kg1 K1) layer thickness (m) measurement errors collector overall heat removal transfer factor global solar radiation per unit area incident on the collector (W m2) H Hessian matrix h heat transfer coefficient (W m2 K1) Q(t) crenel excitation (W m2) p Laplace parameter (s1) thermal contact resistance (W1 K m2) Rc thermal resistance (W1 K) RTh S surface (m2) t time (s) heating time (s) tc ambient temperature (K) Ta Tcrenel calculated temperature (K) Tmeasured measured temperature (K) Ti, To inlet, outlet fluid temperature (K) mean fluid temperature (K) Tm collector overall energy loss coefficient UL (W m2 K) a Cp ei eb FR G
used the photothermal crenel technique to determine the effective thermophysical properties of a glass powder. Cheheb et al. (2008) measured thermal radiative and conductive properties of semitransparent materials using a photothermal crenel method. Faugeroux et al. (2004) and Mzali et al. (2003/2004) theoretically used the flash method to estimate thermal properties of samples which are assumed to be opaque and homogeneous. In this work, on the one hand we present a theoretical and an experimental study allowing the identification of the thermophysical properties of a bi-layer system composed of a heat exchanger section and the fouling deposited on. In brief, the fouling deposited on the heat exchanger tubes is caused by the presence of inorganic salt, small quantities of organic materials and products of corrosion in the water. The experimental bench uses a photothermal method using a crenel heat excitation. The system under investigation is submitted to a crenel heat flux on its front face. The temperature response, during and after irradiation, is measured at the opposite face using a thermocouple. In this paper, the identification of thermophysical properties is performed by an iterative procedure based on minimizing the ordinary least squares function comparing the measured temperature on the rear face of the bi-layer system during a crenel heating flux to the response of a theoretical model with the Gauss–Newton method. The optimization algorithm developed to solve the problem is
V Xij x y z Z
surrounding air speed (m s1) sensitivity coefficients space coordinate (m) space coordinate (m) space coordinate (m) relative sensitivity coefficients
Greek symbols b vector of estimated parameters g thermal efficiency Laplace temperature on the front face of the hf first layer (K s) Laplace temperature on the rear face of the sechr ond layer (K s) k thermal conductivity (W m1 K1) q density (Kg m3) standard deviation of the measurement errors rn (K) standard deviation of the estimated parameters rbi (K) W heat flux density (W m2)
efficient and the unknown parameters are estimated with a good accuracy. On the other hand, in order to point out the effect of the fouling formation, a set of measures of the thermal efficiency and the heat loss coefficient of a flat plate solar collector is realized. It has been shown that the deposited fouling decreases the thermal efficiency and increases the collector heat loss coefficient. 2. System description and mathematical model A heat exchanger can be defined as any device that transfers heat from one medium to another (Kraus, 2003). They are widely used in space heating, refrigeration, air conditioning, power plants, chemical plants, petrochemical plants, petroleum refineries, and natural gas processing. When a heat exchanger is placed in service, the heat transfer surfaces are, presumably, clean. With time, in some services and process industries, the apparatus may undergo a decline in its ability to transfer heat. This is due to the accumulation of heat insulating substances on either or both of the heat transfer surfaces (Kraus, 2003). (Fig. 1) shows a fouling deposited on the internal surface of the heat exchanger. To determine the thermal properties of this fouling accumulation of thickness e1, a section of the system is studied (Fig. 2):
A. Adili et al. / Solar Energy 83 (2009) 1619–1628
1621
Fig. 4. One-dimensional boundary condition. Fig. 1. System under investigation.
e1 and e2. Their interface is characterized by an imperfect contact (thermal contact resistance Rc). The thermal properties and densities of both layers are assumed to be uniform and constant. The convective and radiative heat transfer on the two faces with the uniform environment are expressed by two heat transfer coefficients h1 and h2 (Albouchi et al., 2005). The transient temperature distribution in the sample can be obtained by solving the one-dimensional heat equation for each layer: ki Fig. 2. Section of the heat exchanger.
This sample is subjected to a constant crenel heat flux density Q(t) (W m2) at t0 = 0 on the upper face as shown in Fig. 3: W 0 6 t 6 tc ð1Þ QðtÞ ¼ 0 t > tc where tc is the crenel duration and W is its amplitude as shown in Fig. 3. In the one-dimensional experimental design shown in Fig. 4, the sides of the sample were insulated while an imposed heat flux was applied across the entire top surface. 3. Formulation of the direct problem 3.1. Energy equation The model assumes one-dimensional heat flux through a two-layer sample constituted by two materials of thickness
@ 2 T i ðx; tÞ @T i ðx; tÞ ¼ qi C pi 2 @x @t
i ¼ 1; 2
ð2Þ
where Ti is the temperature of the layer i. Coupled to initial and boundary conditions: T i ¼ 0 at t ¼ 0 @T 1 ðe1 ; tÞ ¼ QðtÞ h1 T 1 ðe1 ; tÞ at x ¼ e1 k1 @x @T 1 ð0; tÞ 1 k1 ¼ ðT 1 ð0; tÞ T 2 ð0; tÞÞ at x ¼ 0 @x Rc @T 1 ð0; tÞ @T 2 ð0; tÞ ¼ k2 at x ¼ 0 k1 @x @x @T 2 ðe2 ; tÞ ¼ h2 T 2 ðe2 ; tÞ at x ¼ e2 k2 @x
ð3Þ ð4Þ ð5Þ ð6Þ ð7Þ
3.2. The thermal quadrupole formalism To solve the system (2)–(7), the thermal quadrupole formalism is used. The entire system can be described in Laplace space as: " # hf 1 0 A1 B 1 1 R c A 2 B 2 ¼ wð1expðptc ÞÞ C 2 D2 h1 1 C 1 D1 0 1 p 1 0 hr A B hr ¼ ð8Þ h2 1 0 C D 0 Here hf and hr are the Laplace transforms of the front and the rear face temperatures of the sample, respectively. The coefficients Ai, Bi, Ci and Di depend on the Laplace parameter p, on the thickness ei of the layer i, and on the thermophysical properties of the material. They are given by:
Fig. 3. Principle of the crenel transient method.
Ai ¼ Di ¼ coshðai ei Þ; C i ¼ ki ai sinhðai ei Þ qffiffiffi Bi ¼ ki1ai sinhðai ei Þ; ai ¼ api
ð9Þ
1622
A. Adili et al. / Solar Energy 83 (2009) 1619–1628
In the Laplace space, the rear face temperature is given by: hr ðpÞ ¼
W ½1 expðptc Þ pC
ð10Þ
where W is the density of the crenel heating flux (Fig. 3). With dimensional parameters, the rear face temperature is a function of several dimensional parameters given by the following expression: hr p ¼ f ðp; tc ; a1 ; k1 ; q1 ; C p1 ; e1 ; h1 ; a1 ; k2 ; q2 ; C p2 ; e2 ; h2 ; Rc ; WÞ ð11Þ Due to the large number of parameters encountered in the mathematical model, this study is presented in dimensionless space with dimensionless parameters (Albouchi et al., 2005). Taking into account that we work at ambient temperature, we have used the approximation of equal heat transfer coefficient on the two sample faces (h1 = h2). The rear face temperature in Laplace space for a crenel heating excitation is given by: hr ¼
b3 ð1 expðs21 tc Þ b21 s21 ðd þ u þ /Þ
ð12Þ
where: d ¼ ½s1 chðs2 Þshðs1 Þ þ b4 b5 s21 shðs2 Þshðs1 Þ ð13Þ þ b5 s1 chðs1 Þshðs2 Þ 1 1 chðs1 Þshðs2 Þ þ b4 chðs1 Þchðs2 Þ þ shðs1 Þchðs2 Þ u ¼ b22 b 5 s1 s1 ð14Þ / ¼ b2 ½2chðs1 Þchðs2 Þ þ b4 b5 s1 chðs1 Þshðs2 Þ 1 shðs1 Þshðs2 Þ þ b4 s1 shðs1 Þchðs2 Þ þ b5 þ b5 The dimensionless parameters are defined by: qffiffiffi qffiffiffiffi 8 > s1 ¼ bp1 ; b1 ¼ ae21 ; s2 ¼ ap2 e2 ; tc ¼ b1 tc > > 1 > < b2 ¼ hek11 ; b3 ¼ q1 Cwp1 e1 > > 12 > > k q C : b4 ¼ Rec1k1 ; b5 ¼ k21 q21 Cp2 p1
ð15Þ
ð16Þ
The variation of the reduced temperature Tcrenel(t, b) with time in the usual space domain is calculated using the numerical algorithm proposed by Graver–Stehfest of hr as shown in the following equation (Maillet et al., 2000): n Lnð2Þ X iLnð2Þ T crenel ðt; bÞ ¼ ð17Þ V i hr t t i¼1 4. Sensitivity analysis A sensitivity coefficient, Xij, is defined as the effect that a change in a particular parameter b has on the variable state. Mathematically, the sensitivity coefficients are defined as the first derivative of the measured variable with
respect to the parameters of the model (Popescu et al., 2004):
@T crenel ðt; bÞ
ð18Þ Xb ¼
@b bj – b where b = [b1, b2, b3, b4, b5]. bj are all parameters other than b that remain constant. The larger Xb, the more sensitive Tcrenel(t, b) is to b and the easier the estimation of this parameter. In addition, viewing the sensitivity coefficients can allow insight into the adequacy of the mathematical model and that of the experimental design. The sensitivity matrix is given by the following expression: 2 @T 3 1crenel 1crenel @T@b @b1 5 6 7 6 7 6 7 ð19Þ X ij ¼ 6 7 4 5 @T ncrenel ncrenel @T@b @b1 5 In this work, the convenient finite difference method was implemented. This choice was supported by the fact that the finite difference accuracy obtained in computing these coefficients is generally acknowledged to be sufficient for the experimental optimization. The sensitivity coefficients are given by: X ij ¼
T icrenel ðb1 ; ... ; bj þ db; . .. ; bp Þ T icrenel ðb1 ; ... ; bj ; ... ; bp Þ db ð20Þ
where db = (0, . . . , 0, dbj, . . . , 0) is a small variation of the bj parameter. In this study, whenever the finite difference method is used to compute sensitivity coefficients, several perturbation sizes, db, were tested to verify that reliable results were obtained. Generally, the perturbation sizes used were between 104 and 106 using double precision. When performing a sensitivity study, it is meaningful to examine the reduced sensitivity Zbj which are obtained by multiplying the original coefficients by the parameter referred to (Raynaud, 1999). These coefficients have the same unit as the state variable (temperature). Using the finite difference approximation, reduced sensitivity coefficients are written as: Z bj ¼ bj X ij ¼ bj
T icrenel ðt; bj þ dbÞ T icrenel ðt; bj Þ dbj
ð21Þ
The simultaneous identification of many parameters is possible if different Zbj are uncorrelated in the course of time. In general, these sensitivity coefficients must be large and uncorrelated with each other (Faugeroux et al., 2004). According to (Fig. 5), we note that in the time interval [0, 190 s] the reduced sensitivity coefficients do not have the same forms and their maxima are reached at different times. We can also deduce that these sensitivity coefficients are uncorrelated in this interval, which is a desirable
A. Adili et al. / Solar Energy 83 (2009) 1619–1628 12
1623
Ζ β1 Ζ β2
10
Ζ β3
reduced sensitivity
8
Ζ β4 Ζ β5
6 4 2 0 -2 -4 -6 0
50
100
150
200
250
t(s) Fig. 5. Time variation of reduced sensitivity coefficients of the reduced model.
condition for the estimation of the unknown thermophysical parameters. However, for t > 190 s the sensitivity coefficients of b1, b2 and b4 are linearly dependent, which indicates that the latest measurements are not useful for estimating of these parameters simultaneously. On the other hand, it can be seen that the reduced model is very sensitive to b3 and b5. So price values of these parameters are of capital importance. Besides, according to (Figs. 5 and 6), we remark that the sensitivity of the model to b4 is the smallest one, which shows that the estimation of this parameter is not easy. In order to verify the approximation of equal heat transfer coefficient on the two faces of the sample (h1 = h2) taken in consideration in the reduction of the mathematical model, the time variations of the reduced sensitivity (Z h1 and Z h2 ) of these two parameters are presented in Fig. 7. (Fig. 7) shows that Z h1 and Z h2 have the same shape and the same influence on the model.
-3
1
x 10
0
reduced sensitivity
-1 -2 -3 -4 -5 -6 -7
The comparison of the two curves shows a symmetrical effect on the model, that means it is possible to permute the heat transfer coefficient on the two faces without modifying the mathematical model on the face rear. Besides these two parameters are correlated, for this reason; we have chosen the same heat transfer coefficient on the two sample faces. 5. Inverse parameter estimation method There are many methods to minimize linear or non-linear functions, among them we mention exhaustive searches, simplex exploration, gradient methods (conjugate gradient, Gauss–Newton, Levenberg–Marquart, etc.), iterative minimizations, and adjoint method (Jarny et al., 1991) which are more or less sophisticated. The choice which can lead to discussion that never ends depend on the number of parameters to be estimated and to some extent on the structure of the criterion. The Gauss–Newton method is one of the simplest methods. It is an iterative procedure based on the minimization of a cost function. This gradient method has been successfully used to solve inverse problems of parameter estimation. The unknown parameters are estimated from adjusting the theoretical temperature history obtained from a mathematical model to the measured temperature history. This can be achieved from the minimization of the gap between the measured and the calculated temperature (Albouchi et al., 2005; Cheheb et al., 2008; Mzali et al., 2003/2004) as shown in Fig. 8. The gap between predictions of the model and the physical reality, called cost function, is given by the following expression: SðbÞ ¼
-8 -9 0
Fig. 7. Reduced sensitivity of the model to the heat transfer coefficients on the two faces.
n X
ðT measured T crenel ðt; bÞÞ
2
ð22Þ
i¼1
50
100
150
200
250
t(s) Fig. 6. Time variation of the reduced model to of Zb4.
where: n is the number of measurement, Tmeasured is a vector containing the measured temperatures, and Tcrenel(t, b) is a vector containing the calculated temperatures obtained
1624
A. Adili et al. / Solar Energy 83 (2009) 1619–1628
error, eb, of estimation is given by the following expression (Jarny and Maillet, 1999): eb ¼ b b
ð28Þ
b: is an estimator of the exact parameter b. eb: is the result of many kind of error sources: eb ¼ e1 þ e2 þ e3 þ e4 þ e5
Fig. 8. Schematic of parameter estimation inverse methodology.
from the solution of the direct model by using the current available estimate for the unknown parameter vector b. The cost function S(b) is minimized by differentiating Eq. (22) with respect to each of the unknown vector bj = (b1, b2, b3, b4, b5) and then setting the resulting expression equal to zero yielding the following set of algebraic equation: n X
2
i¼1
@T crenel ðt; bÞ ½T crenel ðt; bÞ T measured ¼ 0 @bj
ð23Þ
This means that: 2X T ½T crenel ðt; bÞ T measured ¼ 0
ð24Þ
Due to the non-linearity of this system of equations, the vector of estimated temperature is expanded in a Taylor’s series keeping only the first order terms as written in Eq. (25): T crenel ðt; bkþ1 Þ ¼ T crenel ðt; bk Þ þ X T ðbkþ1 bk Þ K
ð25Þ
K
where Tcrenel(t, b ) and X are evaluated at the iteration k. Eq. (25) is substituted to Eq. (24) and we obtain the following iterative procedure: bkþ1 ¼ bk þ Db
ð26Þ
where: Db ¼ ½X T X 1 X T ½T measured T crenel The iterative procedure starts with an initial guess b0, and at each step the vector b is modified until: jbkþ1 bki j i 6d jbki j þ n
for i ¼ 1; 2; . . . 5
ð27Þ
The resolution error e1 is directly related to the method used to solve the direct problem. It is also related to the approximations made by the computer. The error e2 is due to the assumption of constant thermophysical properties according to the temperature. If this assumption is not justified, a skew is then present. Uncertainty on the known parameters generates an error e3 on the estimation of the parameters model. It may be that the model is more or less sensitive to these parameters. The identification is more precise if the sensitivities to the supposed known parameters are low and the uncertainties on the supposed known parameters are minimal. The fourth term e4 presents an error of the thermal metrology. It is due to measurement devices and the signal conversion. The dimension of the sensor can perturb measurement, it is then necessary to control its dimension and its fixing mode. With that the non-linearity of the response of the temperature sensors is added. The measured temperatures are affected of an additive and a random noise, which represents the measurement noise e5 that is deliver. The resulting error, eb, that is equal to the sum of the five sources of errors described above is usually assumed to be gaussian, zero mean and with constant variance r2n Albouchi et al., 2005. It is demonstrated that the uncertainty on the estimated parameter b is directly related to uncertainty in the errors of measurement via the inverse of the Hessian matrix H = XTX as shown in the following equation: 2 2 3 rb1 r2b1;5 6 7 r2b2 7 6 6 7 1 2 T ð30Þ covðbÞ ¼ 6 7 6 7 ¼ ½X X rn 6 7 5 4 r2b5 r2b1;5 The standard deviation on the estimated parameters is given by: 1 1=2
where d is a small number that must be chosen by the investigator (typically 106) and n < 1010 prevents overflow if bk1 ¼ 0. In order to estimate the accuracy of the results, there are many kinds of errors that leads to a standard deviation of the estimated parameters b have to be determined. 6. The various sources of measurement errors During a parametric inversion, the estimated values are not exact but are rather sullied with a certain error. The
ð29Þ
rbi ¼ rn ½diagðX T X Þ
;
i ¼ 1; . . . ; 5
ð31Þ
Eq. (30) shows that the inverse of XTX can amplify the measurement noise. We understand well that maximizing the determinant XTX will minimize the standard deviation of each estimated parameter. According to statistical rules, the exact parameters bj are included in a confidence interval. When the parameters are estimated with 99% confidence level (Mzali et al., 2003/ 2004), this interval is given by the following equation: b 2:576rb 6 b 6 b þ 2:576rb
ð32Þ
A. Adili et al. / Solar Energy 83 (2009) 1619–1628
Hence, the interest to evaluate the error on the estimated parameter vector values, resulting from these modeling errors is evident. 7. Experimental study 7.1. Experimental setup The experimental apparatus is schematically shown in Fig. 9. It involves a stabilized power, a heat source, a sample to be characterized, a thermocouple, a data acquisition system and a computer. The investigated sample, composed of two layers, is a section of a heat exchanger with deposited fouling. The first layer is composed of the fouling of thickness e1 = 0.1 mm, the second layer is the copper of the heat exchanger of thickness e2 = 1 mm. In order to put under the conditions of the one-dimensional heat transfer, the sides of the sample were insulated while an imposed crenel heat flux was applied across the entire top surface using a halogen lamp for 15 s. The experimental temperature at the rear face of the sample is measured by a K-type thermocouple located on the central section of the sample. Transient experimental temperatures are read and recorded with a data acquisition unit (Agilent 34970A) which allows transferring data to a computer via an RS-232 interface. 7.2. Results and discussion The measurement is performed for 250 s with a constant heat flux density, Q = 1000 W m2, for 15 s. The sampling interval is set as 0.25 s throughout the entire temperature recording. The initial temperature is measured as Ta = 300 K. Input parameters for the estimation are summarized in Table 1. Due to the relation between the uncertainty on the estimated parameters bi and the measurement errors via the inverse of the Hessian matrix H as shown in Eq. (31), the measurement errors fluctuation was recorded for a finite duration before turning on the heating source (Fig. 10). A Gaussian distribution was found with a zero mean and a standard deviation equal to 9.62 103 K.
1625
Table 1 Input parameters for the estimation. 1000 W m2 0.1 mm 1 mm 300 K 15 s
Power density (Q) Fouling thickness (e1) Heat exchanger thickness (e2) Reference temperature (Ta) Crenel duration (tc)
Fig. 10. Recorded noise measurement.
Using the measured and calculated temperature and the standard deviation of the measurement errors, the inverse parameter estimation is conducted to determine the unknown parameters. The results of the estimation of the unknown thermophysical properties are shown in Table 2; the measurement uncertainty of each parameter is calculated by considering the global effect of all kind of measurement errors described above. In this study, thermophysical parameters of the second layer of copper p areffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi fixed. Its thermal diffusivity a2 and thermal effusivity k2 q2 C p2 are fixed at known values of 1.16 104 m2 s1 and 1371.898 106 W m2 K1 s1/2, respectively. Using the definition given in Eq. (16), the unknown dimensional parameters of the fouling can be calculated. These parameters are the thermal diffusivity (a1), the thermal conductivity (k1), the volumetric heat capacity (q1Cp1), the global heat transfer coefficient (h), and the contact resistance Rc; their values are given in Table 3. The relative uncertainties on calculated parameters are obtained using the following expressions:
Table 2 Estimated paramaters.
Fig. 9. Experimental device.
b1 b2 b3 b4 b5
Estimated values
Standard deviation, rbi
Relative uncertainty, %
87 9.44 104 4.83 0.36 19.1884322
2.44 102 1.57 106 5.23 102 9.76 103 0.23
0.0723 0.428 2.79 6.98 3.09
1626
A. Adili et al. / Solar Energy 83 (2009) 1619–1628
Table 3 Calculated dimensional parameters. Parameters
Values
2 1
Relative uncertainty, % 6
0.87 10 2070 1.8 17 2 105
a1 (m s ) q1Cp1 (KJ/m3 K) k1 (W m1 K1) h (W m2 K1) Rc (W m2 K1)
2.072 4.790 6.862 8.29 14.842
Da1 Db1 De1 ¼ þ2 a1 b1 e1 Dq1 C p1 Db3 DW De1 þ ¼ þ q1 C p1 b3 e1 W Dk Da1 Dq1 C p1 þ ¼ a1 q1 C p1 k Dh Db2 Dk De1 þ ¼ þ b2 e1 h k DRc Db4 Dk De1 þ ¼ þ Rc b4 e1 k
Fig. 12. Residual between calculated and measured temperatures.
characteristics that the measured errors: zero mean and a standard deviation equal to 9.75 103 K. We can thus consider that these estimations are completely acceptable.
The thickness of the deposited fouling and the heat flux density are known with an accuracy of 1%. According to (Table 3), we can observe that the estimation of the thermophysical properties of the deposited fouling is good. They are calculated with an accuracy of 2.072% for the thermal diffusivity, 4.79% for the volumetric heat capacity, and 6.862% for the thermal conductivity. On the other hand, the thermal contact resistance Rc is estimated with a large uncertainty of 14.842%. The quality of the estimation is analyzed, by comparing the experimental response and the calculated temperature using the estimated values. (Fig. 11) presents a comparison between the measurements and the optimal model using the estimated parameters. The curve shows a good agreement between the measured and calculated temperatures. To qualify this evaluation, the residual of the estimation is represented in Fig. 12. We remark that it is random, centered on zero and it does not present fluctuations. Moreover, It can be observed that this residual has the same 1.2 Calculated Temperature Measured Temperature
1
T/Tmax
0.8 0.6
7.3. Effect of the fouling deposited on the heat transfer Whatever the exact nature of the deposits, an additional thermal resistance to heat transfer is introduced and the operational capability of the heat exchanger is correspondingly reduced. As noted above the consequence of fouling is to form an essentially solid deposit upon the surface, through which heat must be transferred by conduction. If we knew both the thickness and the thermal conductivity of the fouling, we could treat the heat transfer problem simply as another conduction resistance in series with the wall (Fig. 13). In general, we know neither of these quantities and the only possible technique is to introduce the additional resistance as fouling factors in computing the overall heat transfer coefficient. But the present study permits us to know all thermophysical properties of the deposited layer which allows us to know its additional thermal resistance. As well, even though the thickness is unknown, the used inverse method permits us to estimate its value with high accuracy. Then the global thermal resistance is given by the following relation: e1 e2 þ ¼ 0:581 W K1 ð33Þ RTh ¼ k1 S k2 S where S = 102 cm2 is the area of the investigated sample. Let’s recall that, without circulating fluid and before the deposition of the fouling, the thermal resistance is equal to
0.4 0.2 0 -0.2
0
50
100
150
200
250
t(s) Fig. 11. Comparison between measured and calculated temperature.
Fig. 13. Global thermal resistance.
A. Adili et al. / Solar Energy 83 (2009) 1619–1628
1627
The heat loss coefficient (UL) and the thermal efficiency of the collector are determined by the two following empirical relations (Duffie and Beckman, 1991; Kerkni et al., 1990): U L ¼ u1 þ u2 ðT m T a Þ g ¼ F R ðg0 u1 T st
Fig. 14. Flat plate collector.
the one of metal (copper) with which the heat exchanger is made of. We remark thus, during the operational lifetime the thermal resistance increases until 95.7% of its initial value. In order to evaluate the effect of the additional thermal resistance on the heat transfer efficiency, experimental measurements were carried out on a flat plate solar collector before and after its descaling. The method consists in determining the heat losses from the collector indoor in a controlled environment, and the collector thermal efficiency outdoor in natural working conditions (Fig. 14). The heat losses from the collector are determined by taking the following measurements: - The temperature difference of the fluid between the outlet (To) and the inlet (Ti) of the collector. - The ambiant air temperature (Ta). - The heat flow rate of the transfer fluid. - The surrounding air speed (V). The instantaneous thermal efficiency of the collector is determined in quasi steady state conditions. The measurement procedure follows the recommendation for Tunisian solar collector test method (INNORPI, 1991) according to the following conditions: - The collector should be located in such a way that shadow will not fall on the collector at any time during the measurement period. - The incidence angle of direct solar radiation at the plane collector should be less than 40°. - The total solar irradiation on the plane collector during the measurement should be grater than 800 ± 50 W m2. The fluid flow rate should be equal to 0.02 kg s1 ± 10% per square meter of the collector aperture area. - The surrounding air speed should be in a range of 1– 5 m s1. - The fluid temperature at the inlet of the collector should be stable by ±0.1 K.
ð34Þ
u2 T 2st Þ
ð35Þ
where Tm = (Ti + To)/2 is the mean fluid temperature. And Tst = (Tm Ta)/G The values of u1, u2 and FR are experimentally derived constants and g0 is the collector transmissivity–absorptivity product (Lunde, 1981). The operating conditions are given in Table 4: The representation of Tm, G, and V are given in the nomenclature. The results of measurement are presented in Table 5. Results show clearly that the thermal efficiency and the heat loss transfer coefficient depend on the state of the fluid tubes. The maintenance of the fluid circulating tubes can increase the thermal efficiency by 10% and decrease the collector heat loss coefficient by 25%. 9. Conclusion The determination of thermophysical properties from an inverse method is an attractive technique both from the experimental and methodological point of view, because of its accuracy and short time for estimation of parameters. Therefore, an inverse method for the determination of the thermophysical properties of the fouling deposited on the internal surface of fluid circulating tubes at high temperature has been described. The method consists in minimizing a quadratic temperature residual criterion between a measured temperature and a calculated one which has been developed using the quadrupole formalism. Estimated parameters are obtained with high accuracy. The evaluation of confidence intervals resulting from many kinds of errors has been described. These errors have been analyzed and also used to estimate the uncertainties of estimated thermophysical parameters.
Table 4 Operating conditions. Operating conditions
Value
Tm Ta G V
50 K 800 W m2 <4 m s1
Table 5 Thermal efficiency and heat loss coefficient of a tested flat plate collector before and after scaling its fluid conduit. Fluid conduit
Thermal efficiency (g)
Heat loss coefficient (W m2 K)
With fouling Cleaned
0.465 0.520
6.51 5.17
1628
A. Adili et al. / Solar Energy 83 (2009) 1619–1628
To highlight the effect of the characterized fouling on the thermal efficiency, experimental measurements on a flat plate solar collector have been carried out. Results show a spectacular effect on the thermal efficiency. Indeed, a deposited fouling of thickness of an egg shell of 1mm decreases the thermal efficiency by 10%. References Garcia Sandrine, 1999. Experimental design optimization and thermophysical parameter estimation of composite materials using genetic algorithms. Dissertation submitted to the Faculty of the Virginia Polytechnic Institute and State University. Albouchi, F., Fetoui, M., Rigollet, F., Sassi, M., Ben Nasrallah, S., 2005. Optimal design and measurement of the effective thermal conductivity of a powder using a crenel heat excitation. International Journal of Thermal Sciences 44, 1090–1097. Cheheb, Z., Albouchi, F., Ben Nasrallah, S., 2008. Measurement of thermal radiative and conductive properties of Semitransparent materials using a photothermal crenel method. Journal of Quantitative Spectroscopy & Radiative Transfer 109, 620–635. Faugeroux, O., Claudet, B., Be´net, S., Serra, J.J., Boisson, D., 2004. Caracte´risation thermophysique de reveˆtements par me´thode photothermique impulsionnelle en face avant. International Journal of Thermal Sciences 43, 383–401. Mzali, F., Sassi, L., Jemni, A., Ben Nasrallah, S., 2003/2004. Parameter estimation of orthotropic solids with uncertainty in the sensor position:
use of Levenberg–Marquardt and conjugate gradient method. Journal of High Temperatures – High Pressures 35/36, 281–288. Allan D. Kraus, 2003. Heat Transfer Hand Book, Chapter 11. Heat Exchangers, pp. 798.. Maillet, D., Andre´, S., Batsale, J.-C., Degiovanni, A., Moyne, C., 2000. Thermal Quadrupoles, Solving the Heat Equation Through Integral Transforms. Wiley, New York. Popescu, Bogdan, Scudeller, Yves, Brousse, Thierry, Garnier, Bertrand, 2004. Thermal characterization of dielectric thin films using an improved genetic algorithm. Superlattices and Microstructures 35, 239–252. Raynaud, M., 1999. Strategy for experimental design and the estimation of parameters. High Temperatures High Pressures 31, 1–15. Jarny, Y., Ozisik, M.N., Bardon, J.P., 1991. A general optimization method using adjoint equation for solving multidimensional inverse heat conduction. International Journal of Heat and Mass Transfer 34 (11), 2911–2919. Jarny, Y., Maillet, D., 1999. Proble`mes inverses et estimation de grandeurs en thermique. In: Cours C1A, Ecole d’Hiver METTI ’99. Presses Universitaires de Perpignan, pp. 1–50. Recommendation for Tunisian solar collector test method (INNORPI), 1991. Tunisie. Duffie, J.A., Beckman, W., 1991. Solar Engineering of Thermal Processes. Wiley J. & Sons Inc., New York. Kerkni, C., Hamdi, K., Belguith, A., 1990. Aging test of flat plate solar collectors made in Tunisia. 4e`me Congre`s Arabe Internationale sur l’Energie Solaire, Amman, Jordan, pp. 135–141. Lunde, P.J., 1981. New heat transfer factors for flat plate solar collectors. Solar Energy 27 (2), 109–113.
Solar Energy 83 (2009) 1619–1628 www.elsevier.com/locate/solener
Estimation of thermophysical properties of fouling using inverse problem and its impact on heat transfer efficiency A. Adili *, C. Kerkeni, S. Ben Nasralla Centre de Recherches et des Technologies de l’Energie, Laboratoire d’Energe´tique et des Proce´de´s Thermiques, Route touristique de Soliman, B.P. 95, Hammam-Lif 2050, Tunisia Received 11 August 2008; received in revised form 6 April 2009; accepted 24 May 2009 Available online 28 June 2009 Communicated by: Associate Editor Estrada-Gasca
Abstract At high temperature, the circulation of fluid in heat exchangers provides a tendency for fouling accumulation to take place on the internal surface of tubes. In brief, the deposits on heat exchanger tubes are caused by the presence of inorganic salts, of small quantities of organic materials and products of corrosion in the water. From thermophysical point of view, the deposited fouling has harmful effects on the heat exchanger efficiency. Indeed, it increases the thermal resistance which can raise the energy consumption. This study shows an experimental and a theoretical process of estimation of thermophysical properties of the fouling deposited on a section of a heat exchanger and its effects on the heat transfer efficiency. The estimation method is based on the Gauss–Newton algorithm that minimizes the ordinary least squares function comparing a measured temperature and a theoretical one. The temperature response is measured on the rear face of a bi-layer system composed of a section of a heat exchanger and the fouling deposited on during and after a finite width pulse heat flux on its front face. The theoretical temperature, that is a function of the unknown thermophysical properties of the bi-layer system, is calculated by the resolution of the one-dimensional linear inverse conduction problem, and by the use of the quadrupole formalism. The results of the estimation procedure show, on the one hand the efficiency and the stability of the optimization algorithm to estimate the thermophysical properties of the fouling. On the other hand they underline the necessity of the maintenance of fluid circulating tubes at high temperature. Ó 2009 Elsevier Ltd. All rights reserved. Keywords: Inverse problem; Thermophysical properties; Fouling; Thermal quadrupole formalism
1. Introduction Thermophysical characterization of materials is extremely important in the control current for today fabrication processes and in the prediction of structure failure due to thermal stresses. Accuracy in the estimation of the thermal properties can be improved if the experiments are designed carefully (Garcia, 1999). For these reasons, many researchers are interested in the determination of these parameters.
*
Corresponding author. E-mail address: [email protected] (A. Adili).
0038-092X/$ - see front matter Ó 2009 Elsevier Ltd. All rights reserved. doi:10.1016/j.solener.2009.05.014
In many industrial systems, in which there is a water circulation at high temperature, the problem of fouling deposited on the internal surfaces of tubes is found. The fouling formed has a considerable effect on the considered systems efficiency. Indeed, it can modify its thermophysical properties. For these reasons, the determination of the thermophysical properties of the deposited composite material has become very important to understand the performance of the considered system. Many researchers are interested in the determination of thermophysical properties of one or multilayer materials without application of their methodologies and results in industrial or current application. Albouchi et al. (2005)
1620
A. Adili et al. / Solar Energy 83 (2009) 1619–1628
Nomenclature thermal diffusivity (m2 s1) heat capacity (J kg1 K1) layer thickness (m) measurement errors collector overall heat removal transfer factor global solar radiation per unit area incident on the collector (W m2) H Hessian matrix h heat transfer coefficient (W m2 K1) Q(t) crenel excitation (W m2) p Laplace parameter (s1) thermal contact resistance (W1 K m2) Rc thermal resistance (W1 K) RTh S surface (m2) t time (s) heating time (s) tc ambient temperature (K) Ta Tcrenel calculated temperature (K) Tmeasured measured temperature (K) Ti, To inlet, outlet fluid temperature (K) mean fluid temperature (K) Tm collector overall energy loss coefficient UL (W m2 K) a Cp ei eb FR G
used the photothermal crenel technique to determine the effective thermophysical properties of a glass powder. Cheheb et al. (2008) measured thermal radiative and conductive properties of semitransparent materials using a photothermal crenel method. Faugeroux et al. (2004) and Mzali et al. (2003/2004) theoretically used the flash method to estimate thermal properties of samples which are assumed to be opaque and homogeneous. In this work, on the one hand we present a theoretical and an experimental study allowing the identification of the thermophysical properties of a bi-layer system composed of a heat exchanger section and the fouling deposited on. In brief, the fouling deposited on the heat exchanger tubes is caused by the presence of inorganic salt, small quantities of organic materials and products of corrosion in the water. The experimental bench uses a photothermal method using a crenel heat excitation. The system under investigation is submitted to a crenel heat flux on its front face. The temperature response, during and after irradiation, is measured at the opposite face using a thermocouple. In this paper, the identification of thermophysical properties is performed by an iterative procedure based on minimizing the ordinary least squares function comparing the measured temperature on the rear face of the bi-layer system during a crenel heating flux to the response of a theoretical model with the Gauss–Newton method. The optimization algorithm developed to solve the problem is
V Xij x y z Z
surrounding air speed (m s1) sensitivity coefficients space coordinate (m) space coordinate (m) space coordinate (m) relative sensitivity coefficients
Greek symbols b vector of estimated parameters g thermal efficiency Laplace temperature on the front face of the hf first layer (K s) Laplace temperature on the rear face of the sechr ond layer (K s) k thermal conductivity (W m1 K1) q density (Kg m3) standard deviation of the measurement errors rn (K) standard deviation of the estimated parameters rbi (K) W heat flux density (W m2)
efficient and the unknown parameters are estimated with a good accuracy. On the other hand, in order to point out the effect of the fouling formation, a set of measures of the thermal efficiency and the heat loss coefficient of a flat plate solar collector is realized. It has been shown that the deposited fouling decreases the thermal efficiency and increases the collector heat loss coefficient. 2. System description and mathematical model A heat exchanger can be defined as any device that transfers heat from one medium to another (Kraus, 2003). They are widely used in space heating, refrigeration, air conditioning, power plants, chemical plants, petrochemical plants, petroleum refineries, and natural gas processing. When a heat exchanger is placed in service, the heat transfer surfaces are, presumably, clean. With time, in some services and process industries, the apparatus may undergo a decline in its ability to transfer heat. This is due to the accumulation of heat insulating substances on either or both of the heat transfer surfaces (Kraus, 2003). (Fig. 1) shows a fouling deposited on the internal surface of the heat exchanger. To determine the thermal properties of this fouling accumulation of thickness e1, a section of the system is studied (Fig. 2):
A. Adili et al. / Solar Energy 83 (2009) 1619–1628
1621
Fig. 4. One-dimensional boundary condition. Fig. 1. System under investigation.
e1 and e2. Their interface is characterized by an imperfect contact (thermal contact resistance Rc). The thermal properties and densities of both layers are assumed to be uniform and constant. The convective and radiative heat transfer on the two faces with the uniform environment are expressed by two heat transfer coefficients h1 and h2 (Albouchi et al., 2005). The transient temperature distribution in the sample can be obtained by solving the one-dimensional heat equation for each layer: ki Fig. 2. Section of the heat exchanger.
This sample is subjected to a constant crenel heat flux density Q(t) (W m2) at t0 = 0 on the upper face as shown in Fig. 3: W 0 6 t 6 tc ð1Þ QðtÞ ¼ 0 t > tc where tc is the crenel duration and W is its amplitude as shown in Fig. 3. In the one-dimensional experimental design shown in Fig. 4, the sides of the sample were insulated while an imposed heat flux was applied across the entire top surface. 3. Formulation of the direct problem 3.1. Energy equation The model assumes one-dimensional heat flux through a two-layer sample constituted by two materials of thickness
@ 2 T i ðx; tÞ @T i ðx; tÞ ¼ qi C pi 2 @x @t
i ¼ 1; 2
ð2Þ
where Ti is the temperature of the layer i. Coupled to initial and boundary conditions: T i ¼ 0 at t ¼ 0 @T 1 ðe1 ; tÞ ¼ QðtÞ h1 T 1 ðe1 ; tÞ at x ¼ e1 k1 @x @T 1 ð0; tÞ 1 k1 ¼ ðT 1 ð0; tÞ T 2 ð0; tÞÞ at x ¼ 0 @x Rc @T 1 ð0; tÞ @T 2 ð0; tÞ ¼ k2 at x ¼ 0 k1 @x @x @T 2 ðe2 ; tÞ ¼ h2 T 2 ðe2 ; tÞ at x ¼ e2 k2 @x
ð3Þ ð4Þ ð5Þ ð6Þ ð7Þ
3.2. The thermal quadrupole formalism To solve the system (2)–(7), the thermal quadrupole formalism is used. The entire system can be described in Laplace space as: " # hf 1 0 A1 B 1 1 R c A 2 B 2 ¼ wð1expðptc ÞÞ C 2 D2 h1 1 C 1 D1 0 1 p 1 0 hr A B hr ¼ ð8Þ h2 1 0 C D 0 Here hf and hr are the Laplace transforms of the front and the rear face temperatures of the sample, respectively. The coefficients Ai, Bi, Ci and Di depend on the Laplace parameter p, on the thickness ei of the layer i, and on the thermophysical properties of the material. They are given by:
Fig. 3. Principle of the crenel transient method.
Ai ¼ Di ¼ coshðai ei Þ; C i ¼ ki ai sinhðai ei Þ qffiffiffi Bi ¼ ki1ai sinhðai ei Þ; ai ¼ api
ð9Þ
1622
A. Adili et al. / Solar Energy 83 (2009) 1619–1628
In the Laplace space, the rear face temperature is given by: hr ðpÞ ¼
W ½1 expðptc Þ pC
ð10Þ
where W is the density of the crenel heating flux (Fig. 3). With dimensional parameters, the rear face temperature is a function of several dimensional parameters given by the following expression: hr p ¼ f ðp; tc ; a1 ; k1 ; q1 ; C p1 ; e1 ; h1 ; a1 ; k2 ; q2 ; C p2 ; e2 ; h2 ; Rc ; WÞ ð11Þ Due to the large number of parameters encountered in the mathematical model, this study is presented in dimensionless space with dimensionless parameters (Albouchi et al., 2005). Taking into account that we work at ambient temperature, we have used the approximation of equal heat transfer coefficient on the two sample faces (h1 = h2). The rear face temperature in Laplace space for a crenel heating excitation is given by: hr ¼
b3 ð1 expðs21 tc Þ b21 s21 ðd þ u þ /Þ
ð12Þ
where: d ¼ ½s1 chðs2 Þshðs1 Þ þ b4 b5 s21 shðs2 Þshðs1 Þ ð13Þ þ b5 s1 chðs1 Þshðs2 Þ 1 1 chðs1 Þshðs2 Þ þ b4 chðs1 Þchðs2 Þ þ shðs1 Þchðs2 Þ u ¼ b22 b 5 s1 s1 ð14Þ / ¼ b2 ½2chðs1 Þchðs2 Þ þ b4 b5 s1 chðs1 Þshðs2 Þ 1 shðs1 Þshðs2 Þ þ b4 s1 shðs1 Þchðs2 Þ þ b5 þ b5 The dimensionless parameters are defined by: qffiffiffi qffiffiffiffi 8 > s1 ¼ bp1 ; b1 ¼ ae21 ; s2 ¼ ap2 e2 ; tc ¼ b1 tc > > 1 > < b2 ¼ hek11 ; b3 ¼ q1 Cwp1 e1 > > 12 > > k q C : b4 ¼ Rec1k1 ; b5 ¼ k21 q21 Cp2 p1
ð15Þ
ð16Þ
The variation of the reduced temperature Tcrenel(t, b) with time in the usual space domain is calculated using the numerical algorithm proposed by Graver–Stehfest of hr as shown in the following equation (Maillet et al., 2000): n Lnð2Þ X iLnð2Þ T crenel ðt; bÞ ¼ ð17Þ V i hr t t i¼1 4. Sensitivity analysis A sensitivity coefficient, Xij, is defined as the effect that a change in a particular parameter b has on the variable state. Mathematically, the sensitivity coefficients are defined as the first derivative of the measured variable with
respect to the parameters of the model (Popescu et al., 2004):
@T crenel ðt; bÞ
ð18Þ Xb ¼
@b bj – b where b = [b1, b2, b3, b4, b5]. bj are all parameters other than b that remain constant. The larger Xb, the more sensitive Tcrenel(t, b) is to b and the easier the estimation of this parameter. In addition, viewing the sensitivity coefficients can allow insight into the adequacy of the mathematical model and that of the experimental design. The sensitivity matrix is given by the following expression: 2 @T 3 1crenel 1crenel @T@b @b1 5 6 7 6 7 6 7 ð19Þ X ij ¼ 6 7 4 5 @T ncrenel ncrenel @T@b @b1 5 In this work, the convenient finite difference method was implemented. This choice was supported by the fact that the finite difference accuracy obtained in computing these coefficients is generally acknowledged to be sufficient for the experimental optimization. The sensitivity coefficients are given by: X ij ¼
T icrenel ðb1 ; ... ; bj þ db; . .. ; bp Þ T icrenel ðb1 ; ... ; bj ; ... ; bp Þ db ð20Þ
where db = (0, . . . , 0, dbj, . . . , 0) is a small variation of the bj parameter. In this study, whenever the finite difference method is used to compute sensitivity coefficients, several perturbation sizes, db, were tested to verify that reliable results were obtained. Generally, the perturbation sizes used were between 104 and 106 using double precision. When performing a sensitivity study, it is meaningful to examine the reduced sensitivity Zbj which are obtained by multiplying the original coefficients by the parameter referred to (Raynaud, 1999). These coefficients have the same unit as the state variable (temperature). Using the finite difference approximation, reduced sensitivity coefficients are written as: Z bj ¼ bj X ij ¼ bj
T icrenel ðt; bj þ dbÞ T icrenel ðt; bj Þ dbj
ð21Þ
The simultaneous identification of many parameters is possible if different Zbj are uncorrelated in the course of time. In general, these sensitivity coefficients must be large and uncorrelated with each other (Faugeroux et al., 2004). According to (Fig. 5), we note that in the time interval [0, 190 s] the reduced sensitivity coefficients do not have the same forms and their maxima are reached at different times. We can also deduce that these sensitivity coefficients are uncorrelated in this interval, which is a desirable
A. Adili et al. / Solar Energy 83 (2009) 1619–1628 12
1623
Ζ β1 Ζ β2
10
Ζ β3
reduced sensitivity
8
Ζ β4 Ζ β5
6 4 2 0 -2 -4 -6 0
50
100
150
200
250
t(s) Fig. 5. Time variation of reduced sensitivity coefficients of the reduced model.
condition for the estimation of the unknown thermophysical parameters. However, for t > 190 s the sensitivity coefficients of b1, b2 and b4 are linearly dependent, which indicates that the latest measurements are not useful for estimating of these parameters simultaneously. On the other hand, it can be seen that the reduced model is very sensitive to b3 and b5. So price values of these parameters are of capital importance. Besides, according to (Figs. 5 and 6), we remark that the sensitivity of the model to b4 is the smallest one, which shows that the estimation of this parameter is not easy. In order to verify the approximation of equal heat transfer coefficient on the two faces of the sample (h1 = h2) taken in consideration in the reduction of the mathematical model, the time variations of the reduced sensitivity (Z h1 and Z h2 ) of these two parameters are presented in Fig. 7. (Fig. 7) shows that Z h1 and Z h2 have the same shape and the same influence on the model.
-3
1
x 10
0
reduced sensitivity
-1 -2 -3 -4 -5 -6 -7
The comparison of the two curves shows a symmetrical effect on the model, that means it is possible to permute the heat transfer coefficient on the two faces without modifying the mathematical model on the face rear. Besides these two parameters are correlated, for this reason; we have chosen the same heat transfer coefficient on the two sample faces. 5. Inverse parameter estimation method There are many methods to minimize linear or non-linear functions, among them we mention exhaustive searches, simplex exploration, gradient methods (conjugate gradient, Gauss–Newton, Levenberg–Marquart, etc.), iterative minimizations, and adjoint method (Jarny et al., 1991) which are more or less sophisticated. The choice which can lead to discussion that never ends depend on the number of parameters to be estimated and to some extent on the structure of the criterion. The Gauss–Newton method is one of the simplest methods. It is an iterative procedure based on the minimization of a cost function. This gradient method has been successfully used to solve inverse problems of parameter estimation. The unknown parameters are estimated from adjusting the theoretical temperature history obtained from a mathematical model to the measured temperature history. This can be achieved from the minimization of the gap between the measured and the calculated temperature (Albouchi et al., 2005; Cheheb et al., 2008; Mzali et al., 2003/2004) as shown in Fig. 8. The gap between predictions of the model and the physical reality, called cost function, is given by the following expression: SðbÞ ¼
-8 -9 0
Fig. 7. Reduced sensitivity of the model to the heat transfer coefficients on the two faces.
n X
ðT measured T crenel ðt; bÞÞ
2
ð22Þ
i¼1
50
100
150
200
250
t(s) Fig. 6. Time variation of the reduced model to of Zb4.
where: n is the number of measurement, Tmeasured is a vector containing the measured temperatures, and Tcrenel(t, b) is a vector containing the calculated temperatures obtained
1624
A. Adili et al. / Solar Energy 83 (2009) 1619–1628
error, eb, of estimation is given by the following expression (Jarny and Maillet, 1999): eb ¼ b b
ð28Þ
b: is an estimator of the exact parameter b. eb: is the result of many kind of error sources: eb ¼ e1 þ e2 þ e3 þ e4 þ e5
Fig. 8. Schematic of parameter estimation inverse methodology.
from the solution of the direct model by using the current available estimate for the unknown parameter vector b. The cost function S(b) is minimized by differentiating Eq. (22) with respect to each of the unknown vector bj = (b1, b2, b3, b4, b5) and then setting the resulting expression equal to zero yielding the following set of algebraic equation: n X
2
i¼1
@T crenel ðt; bÞ ½T crenel ðt; bÞ T measured ¼ 0 @bj
ð23Þ
This means that: 2X T ½T crenel ðt; bÞ T measured ¼ 0
ð24Þ
Due to the non-linearity of this system of equations, the vector of estimated temperature is expanded in a Taylor’s series keeping only the first order terms as written in Eq. (25): T crenel ðt; bkþ1 Þ ¼ T crenel ðt; bk Þ þ X T ðbkþ1 bk Þ K
ð25Þ
K
where Tcrenel(t, b ) and X are evaluated at the iteration k. Eq. (25) is substituted to Eq. (24) and we obtain the following iterative procedure: bkþ1 ¼ bk þ Db
ð26Þ
where: Db ¼ ½X T X 1 X T ½T measured T crenel The iterative procedure starts with an initial guess b0, and at each step the vector b is modified until: jbkþ1 bki j i 6d jbki j þ n
for i ¼ 1; 2; . . . 5
ð27Þ
The resolution error e1 is directly related to the method used to solve the direct problem. It is also related to the approximations made by the computer. The error e2 is due to the assumption of constant thermophysical properties according to the temperature. If this assumption is not justified, a skew is then present. Uncertainty on the known parameters generates an error e3 on the estimation of the parameters model. It may be that the model is more or less sensitive to these parameters. The identification is more precise if the sensitivities to the supposed known parameters are low and the uncertainties on the supposed known parameters are minimal. The fourth term e4 presents an error of the thermal metrology. It is due to measurement devices and the signal conversion. The dimension of the sensor can perturb measurement, it is then necessary to control its dimension and its fixing mode. With that the non-linearity of the response of the temperature sensors is added. The measured temperatures are affected of an additive and a random noise, which represents the measurement noise e5 that is deliver. The resulting error, eb, that is equal to the sum of the five sources of errors described above is usually assumed to be gaussian, zero mean and with constant variance r2n Albouchi et al., 2005. It is demonstrated that the uncertainty on the estimated parameter b is directly related to uncertainty in the errors of measurement via the inverse of the Hessian matrix H = XTX as shown in the following equation: 2 2 3 rb1 r2b1;5 6 7 r2b2 7 6 6 7 1 2 T ð30Þ covðbÞ ¼ 6 7 6 7 ¼ ½X X rn 6 7 5 4 r2b5 r2b1;5 The standard deviation on the estimated parameters is given by: 1 1=2
where d is a small number that must be chosen by the investigator (typically 106) and n < 1010 prevents overflow if bk1 ¼ 0. In order to estimate the accuracy of the results, there are many kinds of errors that leads to a standard deviation of the estimated parameters b have to be determined. 6. The various sources of measurement errors During a parametric inversion, the estimated values are not exact but are rather sullied with a certain error. The
ð29Þ
rbi ¼ rn ½diagðX T X Þ
;
i ¼ 1; . . . ; 5
ð31Þ
Eq. (30) shows that the inverse of XTX can amplify the measurement noise. We understand well that maximizing the determinant XTX will minimize the standard deviation of each estimated parameter. According to statistical rules, the exact parameters bj are included in a confidence interval. When the parameters are estimated with 99% confidence level (Mzali et al., 2003/ 2004), this interval is given by the following equation: b 2:576rb 6 b 6 b þ 2:576rb
ð32Þ
A. Adili et al. / Solar Energy 83 (2009) 1619–1628
Hence, the interest to evaluate the error on the estimated parameter vector values, resulting from these modeling errors is evident. 7. Experimental study 7.1. Experimental setup The experimental apparatus is schematically shown in Fig. 9. It involves a stabilized power, a heat source, a sample to be characterized, a thermocouple, a data acquisition system and a computer. The investigated sample, composed of two layers, is a section of a heat exchanger with deposited fouling. The first layer is composed of the fouling of thickness e1 = 0.1 mm, the second layer is the copper of the heat exchanger of thickness e2 = 1 mm. In order to put under the conditions of the one-dimensional heat transfer, the sides of the sample were insulated while an imposed crenel heat flux was applied across the entire top surface using a halogen lamp for 15 s. The experimental temperature at the rear face of the sample is measured by a K-type thermocouple located on the central section of the sample. Transient experimental temperatures are read and recorded with a data acquisition unit (Agilent 34970A) which allows transferring data to a computer via an RS-232 interface. 7.2. Results and discussion The measurement is performed for 250 s with a constant heat flux density, Q = 1000 W m2, for 15 s. The sampling interval is set as 0.25 s throughout the entire temperature recording. The initial temperature is measured as Ta = 300 K. Input parameters for the estimation are summarized in Table 1. Due to the relation between the uncertainty on the estimated parameters bi and the measurement errors via the inverse of the Hessian matrix H as shown in Eq. (31), the measurement errors fluctuation was recorded for a finite duration before turning on the heating source (Fig. 10). A Gaussian distribution was found with a zero mean and a standard deviation equal to 9.62 103 K.
1625
Table 1 Input parameters for the estimation. 1000 W m2 0.1 mm 1 mm 300 K 15 s
Power density (Q) Fouling thickness (e1) Heat exchanger thickness (e2) Reference temperature (Ta) Crenel duration (tc)
Fig. 10. Recorded noise measurement.
Using the measured and calculated temperature and the standard deviation of the measurement errors, the inverse parameter estimation is conducted to determine the unknown parameters. The results of the estimation of the unknown thermophysical properties are shown in Table 2; the measurement uncertainty of each parameter is calculated by considering the global effect of all kind of measurement errors described above. In this study, thermophysical parameters of the second layer of copper p areffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi fixed. Its thermal diffusivity a2 and thermal effusivity k2 q2 C p2 are fixed at known values of 1.16 104 m2 s1 and 1371.898 106 W m2 K1 s1/2, respectively. Using the definition given in Eq. (16), the unknown dimensional parameters of the fouling can be calculated. These parameters are the thermal diffusivity (a1), the thermal conductivity (k1), the volumetric heat capacity (q1Cp1), the global heat transfer coefficient (h), and the contact resistance Rc; their values are given in Table 3. The relative uncertainties on calculated parameters are obtained using the following expressions:
Table 2 Estimated paramaters.
Fig. 9. Experimental device.
b1 b2 b3 b4 b5
Estimated values
Standard deviation, rbi
Relative uncertainty, %
87 9.44 104 4.83 0.36 19.1884322
2.44 102 1.57 106 5.23 102 9.76 103 0.23
0.0723 0.428 2.79 6.98 3.09
1626
A. Adili et al. / Solar Energy 83 (2009) 1619–1628
Table 3 Calculated dimensional parameters. Parameters
Values
2 1
Relative uncertainty, % 6
0.87 10 2070 1.8 17 2 105
a1 (m s ) q1Cp1 (KJ/m3 K) k1 (W m1 K1) h (W m2 K1) Rc (W m2 K1)
2.072 4.790 6.862 8.29 14.842
Da1 Db1 De1 ¼ þ2 a1 b1 e1 Dq1 C p1 Db3 DW De1 þ ¼ þ q1 C p1 b3 e1 W Dk Da1 Dq1 C p1 þ ¼ a1 q1 C p1 k Dh Db2 Dk De1 þ ¼ þ b2 e1 h k DRc Db4 Dk De1 þ ¼ þ Rc b4 e1 k
Fig. 12. Residual between calculated and measured temperatures.
characteristics that the measured errors: zero mean and a standard deviation equal to 9.75 103 K. We can thus consider that these estimations are completely acceptable.
The thickness of the deposited fouling and the heat flux density are known with an accuracy of 1%. According to (Table 3), we can observe that the estimation of the thermophysical properties of the deposited fouling is good. They are calculated with an accuracy of 2.072% for the thermal diffusivity, 4.79% for the volumetric heat capacity, and 6.862% for the thermal conductivity. On the other hand, the thermal contact resistance Rc is estimated with a large uncertainty of 14.842%. The quality of the estimation is analyzed, by comparing the experimental response and the calculated temperature using the estimated values. (Fig. 11) presents a comparison between the measurements and the optimal model using the estimated parameters. The curve shows a good agreement between the measured and calculated temperatures. To qualify this evaluation, the residual of the estimation is represented in Fig. 12. We remark that it is random, centered on zero and it does not present fluctuations. Moreover, It can be observed that this residual has the same 1.2 Calculated Temperature Measured Temperature
1
T/Tmax
0.8 0.6
7.3. Effect of the fouling deposited on the heat transfer Whatever the exact nature of the deposits, an additional thermal resistance to heat transfer is introduced and the operational capability of the heat exchanger is correspondingly reduced. As noted above the consequence of fouling is to form an essentially solid deposit upon the surface, through which heat must be transferred by conduction. If we knew both the thickness and the thermal conductivity of the fouling, we could treat the heat transfer problem simply as another conduction resistance in series with the wall (Fig. 13). In general, we know neither of these quantities and the only possible technique is to introduce the additional resistance as fouling factors in computing the overall heat transfer coefficient. But the present study permits us to know all thermophysical properties of the deposited layer which allows us to know its additional thermal resistance. As well, even though the thickness is unknown, the used inverse method permits us to estimate its value with high accuracy. Then the global thermal resistance is given by the following relation: e1 e2 þ ¼ 0:581 W K1 ð33Þ RTh ¼ k1 S k2 S where S = 102 cm2 is the area of the investigated sample. Let’s recall that, without circulating fluid and before the deposition of the fouling, the thermal resistance is equal to
0.4 0.2 0 -0.2
0
50
100
150
200
250
t(s) Fig. 11. Comparison between measured and calculated temperature.
Fig. 13. Global thermal resistance.
A. Adili et al. / Solar Energy 83 (2009) 1619–1628
1627
The heat loss coefficient (UL) and the thermal efficiency of the collector are determined by the two following empirical relations (Duffie and Beckman, 1991; Kerkni et al., 1990): U L ¼ u1 þ u2 ðT m T a Þ g ¼ F R ðg0 u1 T st
Fig. 14. Flat plate collector.
the one of metal (copper) with which the heat exchanger is made of. We remark thus, during the operational lifetime the thermal resistance increases until 95.7% of its initial value. In order to evaluate the effect of the additional thermal resistance on the heat transfer efficiency, experimental measurements were carried out on a flat plate solar collector before and after its descaling. The method consists in determining the heat losses from the collector indoor in a controlled environment, and the collector thermal efficiency outdoor in natural working conditions (Fig. 14). The heat losses from the collector are determined by taking the following measurements: - The temperature difference of the fluid between the outlet (To) and the inlet (Ti) of the collector. - The ambiant air temperature (Ta). - The heat flow rate of the transfer fluid. - The surrounding air speed (V). The instantaneous thermal efficiency of the collector is determined in quasi steady state conditions. The measurement procedure follows the recommendation for Tunisian solar collector test method (INNORPI, 1991) according to the following conditions: - The collector should be located in such a way that shadow will not fall on the collector at any time during the measurement period. - The incidence angle of direct solar radiation at the plane collector should be less than 40°. - The total solar irradiation on the plane collector during the measurement should be grater than 800 ± 50 W m2. The fluid flow rate should be equal to 0.02 kg s1 ± 10% per square meter of the collector aperture area. - The surrounding air speed should be in a range of 1– 5 m s1. - The fluid temperature at the inlet of the collector should be stable by ±0.1 K.
ð34Þ
u2 T 2st Þ
ð35Þ
where Tm = (Ti + To)/2 is the mean fluid temperature. And Tst = (Tm Ta)/G The values of u1, u2 and FR are experimentally derived constants and g0 is the collector transmissivity–absorptivity product (Lunde, 1981). The operating conditions are given in Table 4: The representation of Tm, G, and V are given in the nomenclature. The results of measurement are presented in Table 5. Results show clearly that the thermal efficiency and the heat loss transfer coefficient depend on the state of the fluid tubes. The maintenance of the fluid circulating tubes can increase the thermal efficiency by 10% and decrease the collector heat loss coefficient by 25%. 9. Conclusion The determination of thermophysical properties from an inverse method is an attractive technique both from the experimental and methodological point of view, because of its accuracy and short time for estimation of parameters. Therefore, an inverse method for the determination of the thermophysical properties of the fouling deposited on the internal surface of fluid circulating tubes at high temperature has been described. The method consists in minimizing a quadratic temperature residual criterion between a measured temperature and a calculated one which has been developed using the quadrupole formalism. Estimated parameters are obtained with high accuracy. The evaluation of confidence intervals resulting from many kinds of errors has been described. These errors have been analyzed and also used to estimate the uncertainties of estimated thermophysical parameters.
Table 4 Operating conditions. Operating conditions
Value
Tm Ta G V
50 K 800 W m2 <4 m s1
Table 5 Thermal efficiency and heat loss coefficient of a tested flat plate collector before and after scaling its fluid conduit. Fluid conduit
Thermal efficiency (g)
Heat loss coefficient (W m2 K)
With fouling Cleaned
0.465 0.520
6.51 5.17
1628
A. Adili et al. / Solar Energy 83 (2009) 1619–1628
To highlight the effect of the characterized fouling on the thermal efficiency, experimental measurements on a flat plate solar collector have been carried out. Results show a spectacular effect on the thermal efficiency. Indeed, a deposited fouling of thickness of an egg shell of 1mm decreases the thermal efficiency by 10%. References Garcia Sandrine, 1999. Experimental design optimization and thermophysical parameter estimation of composite materials using genetic algorithms. Dissertation submitted to the Faculty of the Virginia Polytechnic Institute and State University. Albouchi, F., Fetoui, M., Rigollet, F., Sassi, M., Ben Nasrallah, S., 2005. Optimal design and measurement of the effective thermal conductivity of a powder using a crenel heat excitation. International Journal of Thermal Sciences 44, 1090–1097. Cheheb, Z., Albouchi, F., Ben Nasrallah, S., 2008. Measurement of thermal radiative and conductive properties of Semitransparent materials using a photothermal crenel method. Journal of Quantitative Spectroscopy & Radiative Transfer 109, 620–635. Faugeroux, O., Claudet, B., Be´net, S., Serra, J.J., Boisson, D., 2004. Caracte´risation thermophysique de reveˆtements par me´thode photothermique impulsionnelle en face avant. International Journal of Thermal Sciences 43, 383–401. Mzali, F., Sassi, L., Jemni, A., Ben Nasrallah, S., 2003/2004. Parameter estimation of orthotropic solids with uncertainty in the sensor position:
use of Levenberg–Marquardt and conjugate gradient method. Journal of High Temperatures – High Pressures 35/36, 281–288. Allan D. Kraus, 2003. Heat Transfer Hand Book, Chapter 11. Heat Exchangers, pp. 798.. Maillet, D., Andre´, S., Batsale, J.-C., Degiovanni, A., Moyne, C., 2000. Thermal Quadrupoles, Solving the Heat Equation Through Integral Transforms. Wiley, New York. Popescu, Bogdan, Scudeller, Yves, Brousse, Thierry, Garnier, Bertrand, 2004. Thermal characterization of dielectric thin films using an improved genetic algorithm. Superlattices and Microstructures 35, 239–252. Raynaud, M., 1999. Strategy for experimental design and the estimation of parameters. High Temperatures High Pressures 31, 1–15. Jarny, Y., Ozisik, M.N., Bardon, J.P., 1991. A general optimization method using adjoint equation for solving multidimensional inverse heat conduction. International Journal of Heat and Mass Transfer 34 (11), 2911–2919. Jarny, Y., Maillet, D., 1999. Proble`mes inverses et estimation de grandeurs en thermique. In: Cours C1A, Ecole d’Hiver METTI ’99. Presses Universitaires de Perpignan, pp. 1–50. Recommendation for Tunisian solar collector test method (INNORPI), 1991. Tunisie. Duffie, J.A., Beckman, W., 1991. Solar Engineering of Thermal Processes. Wiley J. & Sons Inc., New York. Kerkni, C., Hamdi, K., Belguith, A., 1990. Aging test of flat plate solar collectors made in Tunisia. 4e`me Congre`s Arabe Internationale sur l’Energie Solaire, Amman, Jordan, pp. 135–141. Lunde, P.J., 1981. New heat transfer factors for flat plate solar collectors. Solar Energy 27 (2), 109–113.