An analysis of hydronic heating pavement to optimize the required energy for anti-icing

Applied Thermal Engineering 144 (2018) 278–290

Contents lists available at ScienceDirect

Applied Thermal Engineering journal homepage: www.elsevier.com/locate/apthermeng

Research Paper

An analysis of hydronic heating pavement to optimize the required energy for anti-icing

T

⁎

Raheb Mirzanamadi , Carl-Eric Hagentoft, Pär Johansson Department of Architecture and Civil Engineering, Chalmers University of Technology, Gothenburg SE-412 96, Sweden

H I GH L IG H T S

G R A P H I C A L A B S T R A C T

of the numerical simula• Separation tion model of the HHP using superposition principle.

of the required energy • Optimization for anti-icing using the HHP system. of the number of hours of • Prediction the slippery condition at the road surface.

A R T I C LE I N FO

A B S T R A C T

Keywords: Equivalent temperature Superposition principle Ice free road Required energy Optimization

The aim of this study is to determine the optimal required energy for anti-icing the road surface using a Hydronic Heating Pavement (HHP). A hybrid 3D numerical model is used to simulate the transient anti-icing operation of the HHP system. Moreover, a superposition principle is used to separate the numerical simulation model into two fundamental sub-models: (i) a model with supplying heat to the HHP system and (ii) a model without any heat supply. The criteria for determining the optimal required energy is to increase the temperature at the road surface to provide an ice-free road. To ensure the anti-icing, the temperature threshold of +0.49 °C is added to the required temperature increase at the road surface. The optimal required energy is calculated using a linear programming optimization method. The climate data are obtained from Östersund, a city in the middle of Sweden with long and cold winters. Furthermore, the maximum heat flux supplied to the HHP system is constrained to be 200 W/m2. The results of optimization reveal that the optimal required energy is 106.6 kWh/ (m2 year). Supplying this amount of energy to the HHP system leads to remain only three hours of the slippery conditions at the road surface.

1. Introduction The most common method for mitigation of the slippery conditions at the road surface is to distribute salt and sand [1]. Generally, there are three main goals for usage of salt and sand, namely: anti-icing, de-icing and anti-compaction of snow [2]. Anti-icing is an action to prevent the wet road from freezing, de-icing is an action to melt the existing ice/

⁎

snow at the road surface and finally the anti-compaction is an action to prevent the snow from compacting into a hard crust [3]. Salting and sanding can result in the pollution of the ground and surface waters [4] as well as the corrosion of vehicles and road infrastructures [5]. Furthermore, there is temperature limitation for usage of the salt and sand. The Scandinavian guidelines operate with a lower limitation of −8 °C at the road surface for distributing the salt and sand [6]. For the surface

Corresponding author. E-mail address: [email protected] (R. Mirzanamadi).

https://doi.org/10.1016/j.applthermaleng.2018.08.053 Received 5 April 2018; Received in revised form 29 June 2018; Accepted 18 August 2018 Available online 20 August 2018 1359-4311/ © 2018 Elsevier Ltd. All rights reserved.

Applied Thermal Engineering 144 (2018) 278–290

R. Mirzanamadi et al.

temperature lower than −8 °C, it is difficult to keep the salt concentration high enough to prevent the ice-formation at the road surface [6]. It should be also noted that the salt and sand agents could be lost from road surface due to wind, traffic and runoff [2]. Moreover, the delay in distribution of salt and sand could lead to traffic congestion and severe traffic accidents at the slippery surface [7]. Hence, it is necessary to use environmental friendly methods to mitigate the slippery condition. Two alternative methods for anti-icing the road surfaces are: ice prevention coating [8] and hydronic heating pavement [9]. Ice prevention coating is a concept to prevent or delay the formation of ice at the surface [8]. The concept is initially introduced in 1805 by explaining the nature of solid-liquid interaction and the equilibrium behavior of a droplet on the surface [10]. In 1996, production of the first synthetic super-hydrophilic surfaces (SHSs) sparked a rapid development in the area of the ice prevention coating [10]. The SHS weakens the interaction of the surface with the liquid water by trapping air in the surface texture, resulting in repelling water droplets from the surface before ice-formation occurs [11]. Further developments in this area leads to introduce the flexible SHS with hierarchical micro-/nanostructure which can trap more air in the surface of texture and can increase the anti-icing capacity of the surface [8,12]. The SHSs are environment friendly idea and can save energy for anti-icing the road surface. However, the idea is still in the stage of laboratory and requires more investigations before being used in a practical situation, especially for road application [8]. Another alternative method for anti-icing the road surface is to use a Hydronic Heating Pavement (HHP). The HHP system is a renewable method which consists of embedded pipes. A fluid as the thermal energy carrier circulated along the pipes. The HHP system can be connected to a seasonal thermal energy storage. In summer, the energy is harvested from the road surface and stored in the storage. In winter, the harvested energy, as the only source of the energy, is pumped back to the pipes to mitigate the slippery conditions at the road surface. Different control systems were previously examined for heating the road and bridge surfaces. Li and et al. [13] examined two different heating systems to keep the surface temperature of bridge higher than the freezing point of water. The first heating system was a simple ON/ OFF system with a constant heating power and the second heating system was a dynamic heating system with a variable heating power. The dynamic heating system was developed based on an inverse heat conduction along the thickness direction of the bridge. Both the ON/ OFF system and the dynamic heating system were able to keep the surface temperature of the bridge higher than 0.5 °C. The required heat power for the ON/OFF system was constant of 206 W/m2, however, the maximum required heat power for the dynamic heating system was approximately 150 W/m2. Furthermore, Pahud [14] used the ambient air temperature to control the heating in the HHP system. When the air temperature was below +4 °C, the heating system was turned on and when the air temperature was below −8 °C, the heating system was turned off. It was assumed that for the air temperature below −8 °C the content of air humidity decreases so much that the risk for the iceformation at the road surface becomes quite low [15]. Moreover, Abbasi [16] simulated the HHP system by considering two control systems: (i) when the surface temperature is below 0 °C and (ii) when the surface temperature is below both dew-point temperature and 0 °C. The results showed that the annual required energy for anti-icing the road surface using the first control system was 10 times higher than that using the second control system. Mirzanamadi et al [17] adopted the second control system for the simulation of the HHP system. However, they added a temperature threshold (+0.1 °C to +1.6 °C) to both the dew-point and freezing temperatures for pre-heating the road surface. The results showed that (i) pre-heating the road surface played a significant role in decreasing the slippery condition at the road surface and (ii) adding a temperature threshold to the dew-point temperature led to a higher reduction in the slippery condition, comparing to adding that to the freezing temperature.

Different heating control systems such as feed-forward [18] and self-regulating heating [19] were previously used for heating the buildings. However, the same heating systems are rarely applied for roads, especially for anti-icing the road surface. The aim of this study is to determine the optimal required energy for anti-icing the road surface using the HHP system. The main idea of the optimization is adapted from a study done by Karlsson and Hagentoft [20] which was about the application of model based predictive control for water-based floor heating of building. A superposition principle is used to separate the numerical simulation model of the HHP system into two fundamental sub-models, namely: (i) a model with supplying heat to the HHP system and (ii) a model without any heat supply. The first sub-model is used to obtain an elementary temperature response and the second sub-model is used to obtain the temperature at the surface of an unheated road. The optimal required energy for anti-icing the road surface is determined using a linear programming optimization method. The obtained results from this study associated with the required energy and the remaining number of hours of the slippery conditions are compared with the same results obtained from another study, which investigated the anti-icing operation of the HHP system with a low fluid temperature [17].

2. Hybrid 3D numerical simulation model This section presents: (i) the development of hybrid three-dimensional (3D) numerical simulation model and (ii) validation of the hybrid 3D model using a laboratory experimental results. It should be noted that the model and equations for the hybrid 3D numerical simulation model are based on [21].

2.1. Development of the hybrid 3D numerical simulation model The hybrid 3D numerical simulation model of the HHP system is used to simulate the transient operation of the HHP system for antiicing the road surface. As can be seen in Fig. 1, the HHP system is divided into a finite number of sub-sections. For each sub-section, the 3D model of the HHP system is represented by 2D vertical cross sections which are perpendicular to the pipe direction. All 2D vertical cross sections are serially connected to each other through the convective heat transfer in the fluid, circulating along the pipe. The numerical model simultaneously calculates the transient heat flowing out from the pipes using the finite element method and the temperature decline of the fluid along the pipe using a quasi-steady state assumption. In each 2D vertical cross section, the heat flowing out from the pipes is calculated using an equivalent temperature surrounding the pipe, Teq − pipe (K), and an equivalent thermal resistance between the fluid and the embedded materials, R eq − pipe (m K/W). R eq − pipe consists of three serially coupled thermal resistances of: (i) the thermal resistance at inner pipe surface, RPWS (m K/W), (ii) the pipe material resistance, Rp (m K/W), and (iii) the thermal resistance between the outer pipe wall and an additional annulus, Ri . j (m K/W). Hence, the value of R eq − pipe can be calculated as:

R eq − pipe = Rpipe + Ri, j + RPWS

ln Rp =

( )

279

router rinner

2·π·λpipe

ln Ri, j =

(1)

(2)

( ) ri, j

router

2·π·λi, j

(3)

Applied Thermal Engineering 144 (2018) 278–290

R. Mirzanamadi et al.

Fig. 1. A scheme of the hybrid 3D numerical simulation model of a HHP system (a) the 3D model of a HHP system is divided into a finite number of sub-sections (b) the hybrid 3D model of the HHP system is represented by 2D vertical sections.

RPWS =

1 π·λ f ·Nu ⎛ ⎜where ⎝

4 if Re≤2300 ⎞ Nu = ⎧ 1 0.8 ·Pr 3 if Re > 2300 ⎟ ⎨ 0.023·Re ⎩ ⎠

temperature of the Section 2. Applying this strategy for all sections of pipes; i.e. considering the outlet temperature of one section as the inlet temperature of the next section, will lead to obtain the outlet temperature of last section, T ft ,n + 1 (K) . The value of T ft ,n + 1 is equal to the outlet temperature of whole length of the pipe in the HHP system at time t.

(4)

where router (m) is the outer radius of the pipe, rinner (m) is the inner radius of the pipe, λpipe (W/(m K)) is the thermal conductivity of the pipe, ri, j (m) is the radius of the additional annulus surrounding the pipe, λi, j (W/(m K)) is the thermal conductivity of the surrounding materials, λ f (W/(m K)) is the thermal conductivity of fluid, Nu (–) is the Nusselt number, Re (–) is the Reynolds number and Pr (–) is the Prandtl number. Furthermore, if Tf (K) is the temperature of fluid, circulating through the pipe and qi, j (W/m2) is the heat flow between the fluid and the embedding materials, the value of Teq − pipe can be calculated as:

Teq − pipe = Tf −qi, j ·R eq − pipe

2.2. Experimental validation of the hybrid 3D numerical simulation model The hybrid 3D numerical simulation model of the HHP system is validated by an laboratory experimental results, obtained from the literature [22]. The HHP system in the experimental test consisted of a onelayer concrete pavement and the embedded pipes. The material of pipe was high-density polyethylene. The dimention of the experimental sample was 1000 mm × 1000 mm × 300 mm (length × width × thickness). The distance between pipes was 200 mm (center to center) and the embedded depth of pipes was 50 mm (from surface to center). A schematic view of the experimental test is shown in Fig. 2(a). The thermal properties of materials including thermal conductivity, density and specific heat capacity are given in Table 1. The fluid, used in the experimental test, was 20% ethylene glycol-water mixture. The experimental test was done based on the following boundary conditions: (i) the fluid temperature was 39.6 °C, (ii) the fluid flow rate was 7.1 l/min, (iii) the ambient air temperature at the laboratory room was −1.9 °C, (iv) the relative humidity was 40% and (v) the wind speed was 1 m/s. As shown in Fig. 2(a), there is a symmetrical analysis region between pipes. In order to reduce the computational time, the region of AB-C-D-E-F was selected to make the numerical simulation model. Furthermore, as can be seen from see Fig. 2(b), the hybrid 3D model of the HHP system was represented by four 2D vertical cross sections. For details about selecting the number of cross sections, the reader is referred to [21]. For each 2D vertical cross section, the right, left and bottom boundaries were considered to be adiabatic. The surface of the HHP system was exposed to the convective heat transfer, longwave radiation and condensation/evaporation. From literature [22], the average temperature at the surface between two pipes was given at the time of 120 min. Hence, in order to validate the hybrid 3D numerical simulation model, the heating process was solved for 120 min. The average temperature at the surface along the half distance between pipes at time of 120 min, obtained from the experimental test and the hybrid 3D numerical simulation model, is shown in Fig. 3. As can be seen, the experimental and numerical simulation results are matching well with together, so the maximum

(5)

The longitudinal fluid temperature distribution and the transversal heat transfer process along the pipe are calculated using a steady state solution as:

vf ·π·rinner 2·ρf ·cp, f ·

∂Tf (x ) ∂x

+

Tf (x )−Teq − pipe (x ) R eq − pipe

=0

(6)

where x (m) is the longitudinal length coordinate along the fluid motion, ρf (kg/m3) is the density of fluid, vf (m/s) is the average velocity of fluid through a cross section of pipe and cp, f (J/(kg K)) is the specific heat capacity of fluid. As can be seen from Fig. 1(b), the length of pipe between two 2D vertical cross sections with labels of n and n + 1 is Ln . For this section of pipe, the inlet temperature of the fluid at time t is T ft , n and the outlet temperature is T ft , n + 1. Considering the steady state longitudinal fluid temperature distribution along the pipe, the analytical solution to obtain the outlet fluid temperature is calculated as: t t t −(Ln / ln) T ft , n + 1 = Teq − pipe, n + (T f , n−Teq − pipe, n )· e

(7)

where ln (m) is the characteristic length related to the interaction between the convective heat transfer along the pipe and the transversal heat flow. The value of ln can be calculated as:

ln = R eq − pipe ·vf ·π·rinner 2·ρf ·cp, f

(8)

In this study, the inlet temperature of fluid at Section 1, for each time step is a known value. By applying Eqs. (1)–(8), it will be possible obtain the outlet temperature of fluid at the Section 1 at time step t. The outlet temperature of Section 1 can be used as the inlet

T ft ,1 (K) ,

280

Applied Thermal Engineering 144 (2018) 278–290

R. Mirzanamadi et al.

Fig. 2. The schematic view of hydronic heating pavement (a) a 3D model and (b) a hybrid 3D model used for numerical simulation.

the first sub-model, the surface, bottom and initial temperatures of the road are set to be 0 °C. In addition, the fluid flow rate is set to be constant. In Fig. 4(c), as the second sub-model, the bottom and initial temperatures are set to have the same temperatures as the main numerical simulation model, see Fig. 4(a). It is important to note that in this study the variation of the temperature at the road surface is calculated at a point which is located in the middle between two pipes at the surface of the outlet section. This point is called the control point, see Fig. 4. The temperature variation at the control point, Tsurface (t )(K) can be obtained by summing the temheated perature variation at the control point of the heated road, Tsurface (t )(K) , unheated and that at the control point of the unheated road, Tsurface (t )(K) , as:

Table 1 Thermal properties of the concrete pavement, pipe and fluid (20% glycol). Material

Concrete pavement Pipe Fluid (20% glycol)

Thermal conductivity (W/(m K) )

Density (kg/m3 )

Specific heat capacity (J/(kg K) )

1.6

2380

925

0.42 0.49

1100 1025

1465 3852

heated unheated Tsurface (t ) = Tsurface (t ) + Tsurface (t )

(9)

heated Tsurface (t )

The value of is determined based on the value of the supplied heat to the HHP system, Q (t )(W) . Furthermore, the value of unheated Q (t ) is determined based on the value of Tsurface (t ) . It is worth noting heated (t ). In other that the value of Q (t ) is independent of the value of Tsurface words, there is no temperature feedback from the road surface to the control system in order to change the value of Q (t ) . Considering any temperature feedback in the control system will destroy the linear equation assumption and then the principle of superposition. 4. Equivalent temperature and equivalent heat transfer coefficient at the road surface In this study, it is assumed that all snowfalls at the road are immediately plowed off from the road surface. Hence, the HHP system is not used for snow-melting of the road. However, falling snow, before the snow removal is started, will affect the heat balance of the road surface. Accordingly, the latent heat flux of the snow is ignored and only the sensible heat flux from snow is taken into account for calculation of the heat balance at the road surface. The heat balance is obtained based on seven different heat transfer processes, namely: (i) conductive heat from ground and pipes, qcond (W/m2) (ii) convective heat form ambient air, qconv (W/m2) , (iii) longwave radiation, qlw (W/m2) , (iv) shortwave radiation, qsw (W/m2) , (v) latent heat of evaporation/condensation, qevp / con (W/m2) , (vi) sensible heat from rain, qrain (W/m2) and (vii) sensible heat from snow, qsnow (W/m2) . The heat transfers can be calculated as Eqs. (10)–(16) [23]:

Fig. 3. Comparison of the results related to the laboratory experimental test [22] and the hybrid 3D numerical simulation model (at the time of 120 min).

relative error is 2.4%. It should be noted that the relative error is defined as the difference between the results of the numerical simulation model and experimental test divided to the results of the experimental test.

3. Principle of superposition for the separation of the HHP model A superposition principle is mathematically based on the linear equation system [23]. The superposition principle can be used if the heat transfer process at the road surface and the thermal properties of materials are independent of the temperature level. The separation of the numerical simulation model of the HHP system into two sub-models is illustrated in Fig. 4. As can be seen, the surface boundary of the road is exposed to a constant equivalent heat transfer coefficient, heq (W/(m2 K)) , and an equivalent temperature, Teq (K) . In Fig. 4(b), as 281

qcond = −λ·∇T

(10)

qconv = hc ·(Tambient −Tsurface )

(11)

4 4 qlw = ε ·σ·(Tsky −Tsurface )

(12)

qsw = α·I

(13)

Applied Thermal Engineering 144 (2018) 278–290

R. Mirzanamadi et al.

Fig. 4. Separation of the hydronic heating pavement model into two sub-models using the superposition principle (a) the boundary condition at the road surface consists of the equivalent temperature and the equivalent heat transfer coefficient (b) the sub-model with supplying heat perturbation and (c) the sub-model without supplying any heat.

qevp / con = he ·β . (vambient −vs )

(14)

qrain = ṁ rain ·cp − water ·(Tambient −Tsurface )

(15)

qsnow = ṁ snow ·cp − snow·(Tambient −Tsurface )

(16)

heq = hc + hr

Teq =

tslippery =

heq

(20)

∫0

1year

Tsurface < Tdew ⎞ f ·dt ⎜⎛f = 1 if ⎜⎛ ⎟ Tsurface < Tfreezing ⎝ ⎠ ⎝

else

f = 0⎞⎟ ⎠ (21)

5. Calculation of the temperature at the road surface using the simplified and full simulation models The climate data are obtained from Östersund (63.18 N and 14.5 E), a city in middle of Sweden. The city has long and cold winter period which makes it interesting to simulate the anti-icing operation of the HHP system. The annual mean temperature of the ambient air in Östersund is 2.53 °C [25]. The climate data were generated at a 1 h interval and included: dry-bulb/air temperature (°C), relative humidity (%), wind speed (m/s), dew-point temperature (°C), incoming longwave radiation (W/m2 ), short-wave radiation (W/m2 ) and precipitation (mm/h). Considering the climate data, the value of hc is calculated as 13.21W/(m2 K) with the standard deviation of 4.64 W/(m2 K) and the value of hr is calculated as 3.63W/(m2 K) with the standard deviation of 0.24 W/(m2 K) . A two-dimensional (2D) numerical simulation model of the HHP system [26] is applied to calculate the surface temperature of the unheated road, see Fig. 4(c). In the numerical model, the road structure of the HHP system consists of six layers. The three first layers are constructed from asphalt pavement which have 200 mm depth altogether. The next two layers are constructed from crushed aggregate which have 1080 mm depth altogether. The last layer is ground which has 3720 mm depth. The emissivity value of the road surface is 0.89 (–) and the absorptivity value is 0.78 (–). Furthermore, the pipes are embedded in the second layer of the road. The pipe material is polyethylene (PEX) with a thermal conductivity of 0.4 W/(m K). The outer diameter of pipe is 25 mm and its thickness is 2.3 mm. Since it was assumed that no heat is supplied to the HHP system, the surface of pipe wall is set to be

(17)

In the above equation, the value of Tsurface is unknown. By assuming Tsurface = Tambient , Eq. (17) can be written as:

qlw = hr ·(Tsky−Tsurface )

hc ·Tambient + hr ·Tsky + α·I

It should be noted that, the number of hours of the slippery conditions, tslippery (h) is the number of hours during which the surface temperature is below both dew-point temperature, Tdew (K) and water freezing temperature, Tfreezing (K) . The value of tslippery is calculated as:

where λ (W/(m K)) is the thermal conductivity of the road materials, T (K) is the temperature, Tambient (K) is the ambient air temperature, Tsurface (K) is the road surface temperature, hc (W/(m2 K)) is the convective heat transfer coefficient, ε (–) is the emissivity of the surface, σ (W/(m2 K4)) is the Stefan-Boltzmann constant, Tsky (K) is the sky temperature, α (–) is the solar absorptivity of the surface, I (W/m2) is the solar irradiation, he (J/kg) is the latent heat of evaporation of water, β (m/s) is the moisture transfer coefficient, vambient (kg/m3) is the humidity by the volume of the ambient air, vs (kg/m3) is the humidity by the volume of the saturated air at the surface temperature and ṁ (kg/(m2 s)) is the mass rate of snow/rain per square meter of the surface. For more details about the calculation of the different parameters, the reader is referred to [24]. Considering the rule of difference of two squares, Eq. (12) can be written as: 2 2 qlw = ε ·σ·(Tsky + Tsurface )·(Tsky + T surface )·(Tsky−T surface )

(19)

2 2 where hr = ε ·σ·(Tsky + Tambient )·(Tsky + Tambient )

(18)

(W/(m2

K)) is the heat transfer coefficient associated with the where hr longwave radiation. In order to simplify the calculation of equivalent temperature at the road surface, Teq (K) , it is assumed that the values of qevp / con , qrain and qsnow are negligible. Hence, the value of Teq is obtained only based on qconv , qlw and qsw . Furthermore, in order to apply the superposition principle, annual mean values of the hc and hr are calculated. In order to obtain a better approximation of the surface temperature when the heating system is turned on, the mean values are obtained only for the period during which the road surface is slippery. By considering hc as the annual mean value of the convective heat transfer coefficient and hr as the annual mean value of the heat transfer coefficient of the longwave radiation during which the road surface is slippery, the value of heq (W/(m2 K)) and Teq (K) is calculated as [23]: 282

Applied Thermal Engineering 144 (2018) 278–290

R. Mirzanamadi et al.

adiabatic. The numerical simulation model is created based on the finite element method using the implicit time-stepping. The numerical model is developed in COMSOL Multiphysics 5.2. For more detail about the 2D numerical simulation model of the HHP system, the reader is referred to [26]. eq (K) If the surface temperature induced by Teq and heq is called Tsurface and the surface temperature induced by the all heat fluxes of Eqs. full (K) , the temperature error at the road surface, (11)–(16) is called Tsurface

ME = +0.35 °C, see Table 2. Furthermore, from Table 2 it is seen that for January, February and December as well as for when the road surface is slippery the value of ME is above 0 °C. For this period, full full eq eq tslippery > tslippery > tslippery . In addition, for April, tslippery however ME < 0 °C. It should be noted that for this period, the slippery conditions occur on the 29th and 30th of April when ME = +0.29 °C. Moreover, from March to November, the mean temperature at the surface of the road exposed to the equivalent temperature is 9.32 °C with the standard deviation of 10.37 °C. For this period, the mean temperature at the surface of the road corresponding to the full simulation model is 8.83 °C with the standard deviation of 10.08 °C. During this period, ME < 0 °C full eq < tslippery and consequently tslippery .

full eq eboundary (°C), between Tsurface and Tsurface for each time step will be obtained as: full eq eboundary = Tsurface −Tsurface

(22)

It should be noted that the time step used in the numerical simulation model is 1 h in line with the climate data. Furthermore, the (arithmetic) mean error ME (°C), and the mean absolute error, ME (°C), will be obtained as [27]:

6. Elementary temperature response due to the unit step-change in the heat supply

n

ME =

The elementary temperature response is calculated using the hybrid 3D numerical simulation model of the HHP system, see Section 2. In order to calculate the elementary temperature response, the boundary condition of the road surface is exposed to the equivalent temperature and the equivalent heat transfer coefficient, see Fig. 4(b). Furthermore, the heat supply, Q (t ) , is set to follow Heaviside’s unit step-change function, H (t ) . The Equation of H (t ) is as [20]:

∑i = 1 eboundary, i (23)

n n

MAE =

∑i = 1 |eboundary, i | n

(24)

Considering a standard normal distribution, the lower and upper endpoints of the temperature error at a 95% confidence interval will be calculated as [27]:

H (t ) =

S. D. Lower Endpoint = ME −1.96· n

(25)

S. D. Upper Endpoint = ME −1.96· n

(26)

{10 tt >≤ 00

(27)

Fig. 6(a) shows the variation of the elementary temperature response, uq (°C/W), induced by a unit step-change in the heat supply. The value of uq corresponds to the variation of temperature at the control point, see Fig. 4. The numerical simulation model is run for twoweeks (336 h). For further time, it is assumed that the temperature variation reaches steady-state. To improve the precision of the result, the time step of the simulation is set to be 1 min. Fig. 6(b) shows the results of uq for the first 30 min of the simulation. As can be seen, it takes a few minutes that the supplied heat conducts from the embedded pipes to reach the road surface. Altogether, the variation of uq versus time consists of three parts: (i) an approximately five minutes initial delay, (ii) a sharp increase for whent < 50h and (iii) a gradual increase until the value of uq becomes constant at t = 336h . Applying a unit step-change in the supplied heat will lead to a variation in inlet and outlet temperatures of the fluid. If Tin (K) is the inlet temperature and Tout (K) is the outlet temperature, the variation of the fluid temperature during the time span of two-weeks will be as Fig. 7(a). Furthermore, the variation of Tin and Tout for the first 30 min will be as Fig. 7(b). As can be seen, the fluid temperature immediately responses to the unit step-change in the supplied heat. It is understood that the temperature variation has to start from the initial temperature of the fluid, here 0 °C. However, in the numerical simulation model of the HHP system developed by [21], the time which is required to transport the heat from the supply to the outlet section is not taken into account. In general, the induced time-shift in the inlet and outlet temperatures due to the heat transport from the supply to the outlet section is in the order of 30 s to a few minutes [28]. By considering a constant fluid flow rate, Vḟ (l/min) , a constant specific heat capacity, cp (J/(kg K)) and a constant density, ρ (kg/m3) , the difference between the inlet and outlet temperatures induced by the unit step-change of H (t ) will be constant, see Eq. (28).

where S. D . (−) is the standard division associate with ME and n is the degree of freedom (n = 8760 for a year). full eq The variation of Tsurface and Tsurface for January is shown in Fig. 5. Moreover, Table 2 shows the results of ME and MAE , S. D ., Lower Endpoint and Upper Endpoints as well as the number of hours of the eq slippery conditions. The value of tslippery (h) is the number of hours of the slippery conditions when the road surface is exposed to Teq and heq . In full (h) is the number of hours of the slippery addition, the value of tslippery conditions when the road surface is exposed to the all heat fluxes of Eqs. (11)–(16). The results are calculated separately for a year, different months and for when the road surface is slippery (either full eq tslippery > 0hortslippery > 0h ). As seen from Fig. 5, for January the surface temperature corresponding to the equivalent temperature and the constant heat transfer coefficients is colder than that corresponding to the full simulation full eq ≤ Tsurface model (Tsurface ). The statement is confirmed with

Tin−Tout =

H (t ) Vḟ ·cp·ρ

(28)

In this study, Vḟ = 8 l/min , cp = 3852 J/(kg K) , ρ = 1025 kg/m3 . By considering H (t ) = 1W when the heating system is turned on, the value of Tin−Tout will be 1.9 × 10−3 °C.

Fig. 5. The temperature variation at the road surface during January corresponding to the simplified and the full simulation models. 283

Applied Thermal Engineering 144 (2018) 278–290

R. Mirzanamadi et al.

Table 2 Difference between the surface temperatures corresponding to the simplified and full simulation models. (The data are obtained from the numerical simulation model for the climate data of Östersund). (Standard deviation is related to the mean equivalent temperature). Time

ME (°C)

MAE (°C)

S. D.(°C)

Lower endpoint (°C)

Upper endpoint (°C)

eq tslippery (h)

full tslippery (h)

Year Jan. Feb. Mar. Apr. May June July Aug. Sept. Oct. Nov. Dec.

−0.28 0.35 0.20 −0.42 −0.26 −0.62 −0.80 −0.87 −0.77 −0.42 −0.26 −0.02 0.59

0.62 0.37 0.32 0.46 0.65 0.92 1.14 1.10 0.89 0.48 0.28 0.21 0.61

0.91 0.38 0.40 0.55 0.83 1.20 1.34 1.10 1.01 0.48 0.19 0.28 0.57

−0.30 0.32 0.17 −0.46 −0.32 −0.70 −0.90 −0.95 −0.84 −0.45 −0.28 −0.04 0.55

−0.26 0.38 0.23 −0.38 −0.20 −0.53 −0.71 −0.79 −0.70 −0.38 −0.25 0.00 0.63

2073 653 428 60 4 0 0 0 0 0 35 276 617

2021 637 405 78 3 0 0 0 0 0 41 280 577

0.33

0.37

2073

2021

full eq > 0 or tslippery > 0) For when the road surface is slippery (either tslippery

Year

0.35

0.41

0.47

7. Optimal required energy for anti-icing the road surface

Fig. 9(a) illustrates the assembling of two sequences of uT associated with Q101 (W) , the heat pulse between 100 h and 101 h. Furthermore, Fig. 9(b) illustrates the temperature response corresponding to the heat pulse of Q101. As can be seen, there is a delay in the temperature response to reach the maximum value. The delay is due to the time that it takes the heat conducts from the embedded pipes to reach the road surface. The maximum value of uTSP is 0.0011 °C/W which occurs 94 min after beginning of the heat pulse. The maximum value is located between 1 × tp and 2 × tp .

This section presents the results of: (i) the temperature responses at the road surface induced by supplying a heat pulse to the HHP system (ii) the required temperature increase at the road surface and (iii) the optimal required energy for anti-icing the road surface 7.1. Temperature responses at the road surface induced by supplying a heat pulse The supplied heat to the HHP system, Q (t ), is set to follow a stepwise function. As shown in Fig. 8, the supplied heat is composed by a sequence of heat pulses. The magnitude of heat pulses is set to be 1 W. The length of each heat pulse is tp (s) . The value of tp is set to be 60 min in line with the climate data intervals. It should be noted that the unit of tp is set to be minute in consistence with the resolution of the elementary temperature response, see Fig. 6. The temperature response corresponding to a single heat pulse, uTSP (°C/W), can be obtained by assembling two sequences of the elementary temperature responses, uq (°C/W). As presented in Eq. (29), the value of uTSP , induced by the heat pulse of Qn (W) , can be obtained by assembling the elementary temperature responses at tn − 1 = tp·(n−1) and tn = tp·n .

u (t −tp·(n−1))−uq (t −tp·n) tp > 0, n > 0 uTSP (tn ) = ⎧ q ⎨ 0 else ⎩

7.2. Required temperature increase at the road surface for anti-icing the road surface In order to calculate the optimal required energy, it is necessary to calculate the required temperature increase at the road surface for antiicing the road surface. The required temperature increases, T + (°C), can be obtained as:

T+

unheated ⎧ ⎧Tsurface < Tdew unheated ⎪Tdew−Tsurface + Tthreshold if unheated = ⎨Tsurface < Tfreezing ⎨ ⎩ ⎪ 0 ∘C else ⎩

(30)

unheated Tsurface

where Tdew (°C) is the dew-point temperature and (°C) is the unheated surface temperature when no heat is injected to the HHP system, see Fig. 4(c). Moreover, Tthreshold (°C) is a threshold temperature

(29)

a

b

Fig. 6. Elementary temperature response induced by a step change in the supplied heat (a) the results associated with two weeks simulation and (b) the results associate with the first 30 min of the simulation. 284

Applied Thermal Engineering 144 (2018) 278–290

R. Mirzanamadi et al.

a

b

Fig. 7. Variation of the inlet and outlet temperatures induced by a unit step-change in the supplied heat (a) the results associated with two weeks simulation and (b) the results associate with the first 30 min of the simulation.

A × Q ≥ T+

(31)

1

⋯ ⋯ 0 ⎤ ⎡ uTSP (t1) 0 ⎢ ⋮ 0 ⋮ ⎥ ⋮ ⋱ ⎢ i ⋱ ⋮ ⎥ A = ⎢ uTSP (t1) ⋮ 0 ⎥ ⋮ ⋱ ⋮ ⎥ ⎢ ⋮ ⋮ n·tp · n t ⋯ ⎢ n·tp ⋯uTSP (tm) u p (tn ) ⎥ TSP ⎣uTSP (t1) ⎦n·tp× n ⎡ Q1 ⎤ ⎢⋮⎥ Q = ⎢ Qi ⎥ T + = ⎢ ⎥ ⎢⋮⎥ ⎢Qn ⎦ ⎥n ⎣

Fig. 8. The continuous heat supply, Q (t ) , which is composed by a sequence of heat pulses with the length of tp .

+

⎡ T1 ⎤ ⎢ ⋮ ⎥ ⎢ +⎥ ⎢ Ti ⎥ ⎢ ⋮ ⎥ ⎢T + ⎥ ⎣ n·tp ⎦n·tp× 1

In Eq. (31), A is a matrix of the temperature responses associated with all sequences of the heat pulses. The matrix of A is a triangular matrix with n·tp row and n column. Considering a constant value for tp , the value of A only depends on the geometry and the properties of the i (tn ) is the momentary value of the temperature reHHP system. uTSP sponse corresponding to Qn at t = i·tp . Furthermore, Ti+ is the required temperature increase at t = i . Eq. (31) can be exactly solved for A × Q = T +. However, the results of the solution include some negative elements in Q . It understood that negative elements in Q means cooling the road surface and positive elements means heating the road surface. Since, in this study, the HHP system is simulated only for heating the road surface the values of all elements have to be equal to or more than 0 W, Qi ≥ 0W .

to ensure that the temperature increase at the road surface will keep the road surface ice-free. It should be noted that Tthreshold ≥ 0 °C and T + ≥ 0 °C. The minimum value of T + is obtained if Tthreshold = 0 °C. Furthermore, the statement of T + > 0 °C is true if the supplied heat to the HHP system is already conducted to the road surface, otherwise T + = 0 °C. As it is presented in Eq. (31), the temperature increase at the road surface can be obtained by multiplying the temperature response, corresponding to each heat pulse, by the magnitude of the heat pulse. The temperature increase at the road surface is desired to be equal to or more than T +.

a

b

Fig. 9. The temperature response corresponding to Q101, the heat pulse between 100 h and 101 h (a) the heat pulse and the elementary temperature responses at the time range from 98 h to 110 h and (b) the temperature response for the time range from 80 h to 240 h. 285

Applied Thermal Engineering 144 (2018) 278–290

R. Mirzanamadi et al.

Fig. 10. Optimization results associated with anti-icing the road surface during January using the simplified simulation model (a) the optimal heat fluxes supplied to the HHP system (b) the surface and dew-point temperatures. analytical unheated Tsurface = A × Qoptimal + Tsurface

7.3. Calculation of the optimal required energy for anti-icing the road surface

(where

The annual required energy for anti-icing the road surface, Er (kWh/(m2 year) , is calculated as:

Er =

1year

∫i=1

qi ·dt

where

qi =

Qi c·L

where qi (W/m2) is the heat flux supplied to the HHP system at time i, c (m) is the distance between the pipes and L (m) is the pipe length. In this study, c = 100mm and L = 50m . By considering qmax (W/m2 ) as the maximum possible heat flux which can be supplied to the HHP system and Er as the objective function of the optimization, the optimal required energy for anti-icing the road surface will be calculated as:

analytical < 0 °C, the value of lustrated in Fig. 10(b). As shown, for when Tsurface analytical Tsurface is equal to or more than the value of Tdew , see Eq. (30). Table 3 presents the maximum and mean heat fluxes, the number of hours during which the heating system is turned on as well as the optimal required energy for anti-icing the road surface during a year and different months. As can be seen, the maximum heat flux for when the heating system is turned on flow is 955.8 W/m2 which occurs in December. Moreover, the minimum heat flux for when the heating system is turned on is 143.75 W/m2 which occurs in October. The annual mean heat flux is 67.86 W/m2 . By considering 1197 h for when the heating system is turned on, the annual optimal required energy for anti-icing the road surface will be 81.23 kWh/(m2 year) . Approximately 90% of the required energy is used for anti-icing the road surface from the fist of December to the end of February. Comparing the results related to the different months reveals that the maximum required energy is 30.52 kWh/(m2 month) which occurs in December. Moreover, the minimum required energy is 0.51 kWh/(m2 month) which occurs in April. Interestingly, the maximum value of mean heat flux is 170.84 W/m2 which occurs in April as well. For this period, the heating system is turned on only for 3 h.

(33)

A linear programming optimization was used to solve Eq. (33). The function was solved using Gurobi Optimizer 7.5 [29], interfaced MATLAB R2016.b environment. The program was run in a computer with 64 GB RAM and a processor of Intel (R) Core (TM) i7-6900 K CPU @ 3.2 GHz. It is worth noting that Gurobi Optimizer solved Eq. (33) in less than five minutes, while, solving the same problem using MATLAB Optimization toolbox required much longer time (the solver was still busy after 7 h without any solution). The calculated optimal heat flux supplied to the HHP system, qoptimal (W/m2 ), can be used to obtain the surface temperature of the heated analytical road. The analytical solution to obtain the surface temperature, Tsurface (°C), can be written as:

Table 3 Required energy for anti-icing the road surface. (The boundary conditions is the same as Fig. 4(a) and the climate data are obtained from Östersund). Time

Year Jan. Feb. Ma. Apr. May June July Aug. Sept. Oct. Nov. Dec.

(34)

By considering Tthreshold = 0 °C and setting no constrain to qmax , the optimization results associated with the optimal heat fluxes supplied to the HHP system and the variation of the surface temperature during January will be as Fig. 10. As can be seen in Fig. 10(a), all the heat fluxes are above 0 W/m2 . The maximum heat flux during January is less analytical and Tdew are ilthan 600 W/m2 . Furthermore, the variation of Tsurface

(32)

A × Q ≥ T+ min Er such that ⎧ 2 ⎨ ⎩ 0 W/m ≤ Q/(c·L) ≤ qmax

Qoptimal = qoptimal × c × L)

Max. heat flux when the heating

Mean heat flux when the heating system is on (W/m2 )

Number of hours during which the heating system is on (h)

Optimal required energy(kWh/m2)

system is on (W/m2 ) 955.80 590.82 312.96 167.04 226.38 – – – – – 143.75 432.68 955.80

67.86 69.71 60.56 38.57 170.84 – – – – – 32.75 46.61 86.45

1197 355 248 37 3 0 0 0 0 0 26 175 353

81.23 24.75 15.02 1.43 0.51 0 0 0 0 0 0.85 8.16 30.52

286

Applied Thermal Engineering 144 (2018) 278–290

R. Mirzanamadi et al.

Table 4 Temperature difference at the road surface between the analytical calculation and the numerical simulation model as well as the maximum inlet/outlet temperatures of fluid (the surface boundary consists of Teq and heq ). (Standard deviation is related to the mean equivalent temperature). Time

ME (°C)

MAE (°C)

S.D. (°C)

Lower endpoint

Upper endpoint

analytical tslippery (h)

simulation (h) tslippery

Max. Tin (°C)

Max. Tout (°C)

Year Jan. Feb. Mar. Apr. May June July Aug. Sept. Oct. Nov. Dec.

−0.04 −0.17 −0.06 −0.06 −0.01 0.00 0.00 0.00 0.00 0.00 0.00 −0.01 −0.11

0.11 0.21 0.12 0.10 0.09 0.10 0.10 0.09 0.08 0.06 0.04 0.07 0.24

0.19 0.34 0.16 0.10 0.11 0.13 0.12 0.11 0.10 0.08 0.06 0.12 0.40

−0.04 −0.17 −0.06 −0.06 −0.01 0.00 0.00 0.00 −0.01 0.00 0.00 −0.01 −0.11

−0.04 −0.17 −0.06 −0.06 −0.01 0.00 0.00 0.00 0.00 0.00 0.00 −0.01 −0.11

0 0 0 0 0 0 0 0 0 0 0 0 0

121 10 23 1 0 0 0 0 0 0 13 37 37

23.69 11.11 7.41 1.89 7.45 – – – – – 2.94 11.31 23.69

15.99 6.75 4.99 1.23 5.61 – – – – – 2.22 7.2 15.99

0.08

0.12

0

121

–

–

simulation > 0h ) For when the road surface is slippery (tslippery

year

0.10

0.10

0.11

8. Application of the optimal heat fluxes to the full simulation model of the HHP system

obtained from the numerical simulation model is colder than that obtained from Eq. (34). The colder surface resulted in the occurrence of the slippery conditions at the road surface. As mentioned in Section 7.3, the optimal heat fluxes were initially calculated by considering Tthreshold = 0 °C. Hence, a small variation in the surface temperature can result in the occurrence of the slippery conditions at the road surface. As it is presented in Table 4, for when the heating system is turned on, the annual maximum inlet temperature of fluid is 23.69 °C and the annual maximum outlet temperatures is 15.99 °C. It is worth noting that the maximum values of inlet and outlet temperatures are influenced by both the supplied heat fluxes to the HHP system and the ambient conditions. For example, the mentioned maximum inlet and outlet temperatures are related to the supplied heat flux of 824.84 W/m2 , occurring on the 2nd of December. However, the maximum heat flux supplied to the HHP system during a year is 955.8 W/m2 , occurring on the 6th of December. The inlet and outlet temperatures corresponding to the maximum heat flux are respectively 22.14 °C and 13.07 °C.

In this section, the calculated optimal heat flux, qoptimal , are imported to the numerical simulation model of the HHP system. The numerical simulation model is run for two cases: (i) when the road surface is exposed to the simplified boundary conditions, consisting of Teq and heq as well as (ii) when the road surface is exposed to the all heat fluxes of Eqs. (11)–(16); i.e. a full simulation model is run. For the second case, a temperature threshold is used in Eq. (30) and a maximum value is considered for the heat fluxes supplied to the HHP system in Eq. (33). 8.1. Simulation of the HHP system using the calculated optimal heat fluxes and considering the simplified boundary conditions at the road surface The numerical simulation model of the HHP system is used to obsimulation tain: (i) the surface temperature at each time step, Tsurface (°C), (ii) the simulation remaining number of hours of the slippery conditions, tslippery (h) , (iii) the maximum inlet temperature of fluid, Tin (°C), and (iv) the maximum outlet temperature of fluid, Tout (°C). In order to determine the difference between the surface temperature obtained from the numerical simulation model and that obtained analytical simulation from Eq. (34), The values of Tsurface (°C) and Tsurface are compared simulation with each other. The temperature error, emethod (°C), between Tsurface analytical and Tsurface for each time step is calculated as: analytical simulation emethod = Tsurface −Tsurface

8.2. Full simulation model of the HHP system using the calculated optimal heat fluxes In order to reduce the risk for the occurrence of the slippery conditions due to the temperature variation at the road surface, it is important to solve the optimization solution by considering a temperature threshold in Eq. (30). The value of Tthreshold can be calculated based on the investigated temperature errors in this study. Two types of the errors were investigated: (i) the error corresponding to the unheated surface temperature, see Table 2 and (ii) the error corresponding to the heated surface temperature, see Table 4. The value of Tthreshold is obtained by summing the values of ME , Lower Endpoint and Upper Endpoint presented in Tables 2 and 4. Furthermore, in order to obtain a proper Tthreshold for anti-icing the road surface, only the errors corresponding to the period during which the road surface is slippery is taken into account. The three different values for Tthreshold will be obtained as:

(35)

Table 4 presents the results of ME , MAE , S. D ., Lower Endpoint and Upper Endpoints associate with emethod for a year and different months, the remaining number of hours of the slippery conditions as well as the maximum values of inlet and outlet temperatures. The values ofME , MAE , Lower Endpoint and Upper Endpoint are obtained from Eqs. (23)–(26) by replacing eboundary with emethod . As can be seen in Table 4, the annual value of ME as well as the values of ME from November to April are below 0 °C. ME < 0 °C means that the calculated surface temperature using the numerical simulation analytical simulation > Tsurface . model is warmer than that obtained from Eq. (34), Tsurface In addition, the number of hours of the slippery conditions at the road simulation surface which are calculated based on Tsurface is 121 h. Considering analytical tslippery = 0h , a misinterpretation will occur as: by increasing the temperature at the road surface, the remaining number of hours of the slippery conditions will increase. To correct the misinterpretation, the value of ME is recalculated for the period during which the road surface simulation > 0h ). The results showed that the value of ME is is slippery (tslippery above 0 °C for this period. This means that the surface temperature

corresponding to ME : 0.35 °C + 0.1 °C = 0.45 °C •T corresponding to Lower •T 0.33 °C + 0.08 °C = 0.41 °C corresponding to Upper •T threshold threshold

Endpoint

threshold

Endpoint

0.37 °C + 0.12 °C = 0.49 °C

The value of 0.49 °C is selected to ensure that the variation of temperature at the road surface does not significantly affect the antiicing operation of the HHP system. Fig. 11 illustrates the process of importing the calculated optimal 287

Applied Thermal Engineering 144 (2018) 278–290

R. Mirzanamadi et al.

Fig. 11. Calculation of optimal heat flux supplied to the HHP system for keeping the road surface ice-free at the control point (a) considering a simplified simulation model and (b) considering the full simulation model.

flux supplied to the HHP system would be 970.44 W/m2 . From Table 5, the annual mean heat flux for when the heating system is turned on is 71.59 W/m2 . The maximum mean heat flux is 136.03 W/m2 which occurs during April for when the heating system is turned on only for 3 h. This short period of heating produces 0.68 kWh/(m2 month) energy, which is the minimum value for the required energy. The total number of hours during which the heating system is turned on is 1489 h. This time of heating produces 106.6 kWh/(m2 year) energy. This value is 31% more than that presented in Table 3. Comparing the results associated with the number of hours of the slippery conditions presented in Tables 3 and 5 shows that considering the temperature threshold leads to decrease the remining number of hours of the slippery condition from 121 h for when Tthreshold = 0 °C to only 3 h for when Tthreshold = 0.49 °C. Even for these 3 h, the mean difference between the surface temperature and the dew-point temperature is 0.05 °C with the standard division of 0.02 °C.

heat fluxes from the model (a), at which the road surface is exposed to the simplified boundary conditions, to the model (b), at which the road surface is exposed to the all heat fluxes of Eqs. (11)–(16); i.e. the full simulation model is run. In this section, in order to reduce the maximum inlet temperature of fluid, the maximum heat flux supplied to the HHP system, qmax , is constrained to be 200 W/m2 , see Eq. (33). Table 5 presents the maximum and mean heat fluxes supplied to the HHP system for a year and different months, the number of hours during which the heating system is turned on, the required energy for anti-icing the road surface, the remaining number of hours of the slippery conditions as well as the maximum inlet and outlet temperatures for when the heating system in turned on. As can be seen from Table 5, considering qmax = 200 W/m2 for calculation of the optimal heat fluxes supplied to the HHP system results in the maximum inlet temperature of 11.96 °C and the maximum outlet temperature of 10.06 °C. If there was no constrain for qmax , the maximum inlet temperature would be 27.96 °C and the maximum outlet temperature would be 19.11 °C. Furthermore, the maximum heat

Table 5 Maximum and mean heat fluxes supplied to the HHP system, required energy for anti-icing, temperature difference at the road surface between the analytical calculation and the numerical simulation as well as the maximum inlet/outlet temperatures of fluid (the surface is exposed to the all heat fluxes of Eqs. (11)–(16); i.e. a full simulation model is run). (Standard deviation is related to mean equivalent temperature). Time

Year Jan. Feb. Mar. Apr. May June July Aug. Sept. Oct. Nov. Dec.

Max. heat flux when the

Mean heat flux when the

heating system is on (W/m2 )

heating system is on (W/m2 )

Number of hours during which the heating system is on (h)

(kWh/m2)

200 200 200 173.39 200 – – – – – 172.88 200 200

71.59 68.86 65.15 54.40 136.03 – – – – – 54.64 55.09 90.73

1489 454 298 56 5 0 0 0 0 0 30 226 420

106.60 31.26 19.41 3.05 0.68 0.00 0.00 0.00 0.00 0.00 1.64 12.45 38.10

288

Optimal required energy

simulation (h) tslippery

Max. Tin (°C)

Max. Tout (°C)

3 2 0 0 0 0 0 0 0 0 1 0 0

11.96 7.57 7.25 4.11 8.63 – – – – – 6.26 7.28 11.96

10.06 5.67 5.35 2.52 6.73 – – – – – 4.62 5.64 10.06

Applied Thermal Engineering 144 (2018) 278–290

R. Mirzanamadi et al.

9. Conclusions

Acknowledgements

In this study, the optimal required energy was calculated for antiicing the road surface using the HHP system. The climate data obtained from Östersund, a city in the middle of Sweden. A hybrid 3D numerical model was applied to simulate the transient ant-icing operation of the HHP system. The boundary condition at the road surface was simplified by calculating the constant equivalent heat transfer coefficient and the equivalent temperature. The results showed that the mean difference between the surface temperature corresponding to the simplified boundary conditions at the road surface and that corresponding to the full simulation model, exposed to all heat fluxes of Eqs. (11)–(16), was 0.35 °C with the standard deviation of 0.47 °C for when the road surface is slippery. The superposition principle was used to separate the numerical simulation model to two sub-models: (i) a model with supplying heat to the HHP system and (ii) a model without any heat supply. The results of the first model showed that it takes approximately five minutes to transfer heat from the embedded pipes to the road surface. Moreover, the results of the second model showed that there are annually 2073 h risk for the occurrence of the slippery condition at the road surface. The required temperature increase at the road surface for anti-icing was obtained based on the difference between the dew-point temperature and the surface temperature of the unheated road. Furthermore, the temperature threshold of 0.49 °C was added to the required temperature increase. A linear programming optimization was used to calculate the optimal required energy. The maximum heat flux supplied to the HHP system was set to be 200 W/m2 . The results showed that 106.6kWh/(m2 year) energy was required for anti-icing the road surface. Applying this amount of energy to the full simulation model of the HHP system resulted in remaining only 3 h of risk for the slippery condition at the road surface. The main benefit of anti-icing the road surface based on a predictive heating control system, as described in this study, is that the future heat demand can be predicted before supplying heat to the HHP system. Predicting heat demand makes it possible to supplying heat in advance before the ice is formed at the road surface, resulting in counteracting the delay in the heat conduction from the embedded pipes to the road surface. Counteracting the delay time can be result in the lower energy consumption for anti-icing the road surface. The results obtained from this study was compared with the results obtained from a 2D numerical simulation model of the HHP system [17]. The obtained results from the 2D numerical simulation model was: 154.2 kWh/(m2 year) for the required energy and 2 h for the remaining number of hours of the slippery condition at the road surface. The required energy obtained from the optimization solution is approximately 30% less than that obtained from the 2D model, while the remining number of hours of the slippery conditions for both methods is approximately the same. This study focused only on the road part of the HHP system. Further studies are necessary to investigate the interaction between the road and storages parts of the HHP system. The anti-icing operation of the HHP system is required to be experimentally validated by a filed data. Moreover, the presented results in this study is only valid for when the snowfalls are removed from the road surface. Hence, a further study is needed to obtain the amount of the required energy for melting the existing snow at the road surface. The aim of installing the HHP system is to use solar energy for mitigating the slippery conditions at the road surface, so it is possible to consider the HHP system as a sustainable solution. However, the sustainability of the HHP system should be analyzed by considering additional parameters such as the costs associated with the initial investment, operation and maintenance as well as the impacts of the HHP system on reducing the environment pollutions and improving the traffic safety.

This work was financially supported by Norwegian Public Road Administration and Chalmers University of Technology. The authors are grateful to the reviewers for their valuable comments on the paper. The first author would like to acknowledge Pepe Tan from Chalmers University of Technology and Henrik Karlsson from Rise Research Institute of Sweden for their interest in the work and fruitful discussions. Appendix A. Supplementary material Supplementary data associated with this article can be found, in the online version, at https://doi.org/10.1016/j.applthermaleng.2018.08. 053. References [1] A.K. Arvidsson, The winter model – a new way to calculate socio-economic costs depending on winter maintenance strategy, Cold Reg. Sci. Technol. 136 (2017) 30–36, https://doi.org/10.1016/j.coldregions.2017.01.005. [2] H.R. Vignisdottir, G.K. Booto, R.A. Bohne, B. Ebrahimi, H. Brattebø, H. Wallbaum, Life Cycle assessment of Anti-and De-icing Operations in Norway, CIB World Build. Congr., Tampere, 2016, pp. 441–454. [3] J. Wåhlin, A. Klein-Paste, A salty safety solution, Phys. World 30 (2017) 27–30. [4] E. Asensio, V.J. Ferreira, G. Gil, T. García-Armingol, A.M. López-Sabirón, G. Ferreira, accumulation of de-icing salt and leaching in Spanish soils surrounding roadways, Int. J. Environ. Res. Public Health 14 (2017), https://doi.org/10.3390/ ijerph14121498. [5] L. Anuar, A. Amrin, R. Mohammad, A. Ourdjini, Vehicle accelerated corrosion test procedures for automotive in Malaysia, in: MATEC Wev Conf 2017;90:№. 01040. http://doi.org/10.1051/matecconf/20179001040. [6] H. Norem, Selection of strategies for winter maintenance of roads based on climatic parameters, J. Cold Reg. Eng. 23 (2009) 113–135, https://doi.org/10.1061/(ASCE) 0887-381X(2009) 23:4(113). [7] L. Gáspár, Z. Bencze, Salting route optimization in Hungary, Transp. Res. Proc. 14 (2016) 2421–2430, https://doi.org/10.1016/j.trpro.2016.05.285. [8] S. Zhang, J. Huang, Y. Cheng, H. Yang, Z. Chen, Y. Lai, Bioinspired surfaces with superwettability for anti-icing and ice-phobic application: concept, mechanism, and design, Small 13 (2017) 1–20, https://doi.org/10.1002/smll.201701867. [9] H. Wang, A. Jasim, X. Chen, Energy harvesting technologies in roadway and bridge for different applications – a comprehensive review, Appl. Energy 212 (2018) 1083–1094, https://doi.org/10.1016/j.apenergy.2017.12.125. [10] M.J. Kreder, J. Alvarenga, P. Kim, J. Aizenberg, Design of anti-icing surfaces: Smooth, textured or slippery? Nat. Rev. Mater. 1 (2016), https://doi.org/10.1038/ natrevmats.2015.3. [11] J. Lv, Y. Song, L. Jiang, J. Wang, Bio-inspired strategies for anti-icing, ACS Nano 8 (2014) 3152–3169, https://doi.org/10.1021/nn406522n. [12] Y. Shen, G. Wang, J. Tao, C. Zhu, S. Liu, M. Jin, et al., Anti-icing performance of superhydrophobic texture surfaces depending on reference environments, Adv. Mater. Interfaces 4 (2017) 1–7, https://doi.org/10.1002/admi.201700836. [13] K. Li, N. Hong, Dynamic heat load calculation of a bridge anti-icing system, Appl. Therm. Eng. 128 (2018) 198–203, https://doi.org/10.1016/j.applthermaleng.2017. 09.024. [14] D. Pahud, Simulation Tool for the System Design of Bridge Heatig for Ice Prevention with Solar Heat Stored in a Seasonal Ground Duct Store, User Manual, Lugano, Switzerland, 2008. [15] D. Pahud, Serso, Stockage Saisonnier Solaire pour le Dégivrage d’un Pont, Rapport Final (in French), Bern, Switzerland, 2007. [16] M. Abbasi, Non-skid Winter Road, Investigation of Deicing System by Considering Different Road Profiles, Chalmers University of Technology, Gothenburg, Sweden, 2013. [17] R. Mirzanamadi, C.-E. Hagentoft, P. Johansson, Hydronic heating pavement with low temperature: the effect of pre-heating and fluid temperature on anti-icing performance, in: 9th Int. Cold Clim. HVAC Conf., Kiruna, Sweden, 2018. [18] C.E. Hagentoft, A.S. Kalagasidis, Effect smart solutions for district heating networks based on energy storage in buildings. Impact on indoor temperatures, Energy Proc. 78 (2015) 2244–2249, https://doi.org/10.1016/j.egypro.2015.11.346. [19] H. Karlsson, Self-regulating floor heating systems in low energy buildings, in: 8th Symp. Build. Phys. Nortic Ctries., Copenhagen, 2008. [20] H. Karlsson, C.-E. Hagentoft, Application of model based predictive control for water-based floor heating in low energy residential buildings, Build. Environ. 46 (2011) 556–569, https://doi.org/10.1016/j.buildenv.2010.08.014. [21] R. Mirzanamadi, C.-E. Hagentoft, P. Johansson, Harvesting solar energy and heating road surfaces using hydronic heating pavement with low temperature-hybrid 3D numerical simulation model, Energies (2018) (to be submitted for publications). [22] H. Wang, Z. Chen, Study of critical free-area ratio during the snow-melting process on pavement using low-temperature heating fluids, Energy Convers. Manage. 50 (2009) 157–165, https://doi.org/10.1016/j.enconman.2008.08.019. [23] C.-E. Hagentoft, Introduction to Building Physics, 1:7. Lund, Sweden:

289

Applied Thermal Engineering 144 (2018) 278–290

R. Mirzanamadi et al.

Technol. 145 (2018) 106–118, https://doi.org/10.1016/j.coldregions.2017.10.006. [27] A.A. Taylor, L.M. Leslie, A single-station approach to model output statistics temperature forecast error assessment, Weather Forecast 20 (2005) 1006–1020, https://doi.org/10.1175/WAF893.1. [28] H. Karlsson, Thermal Modeling of Water-Based Floor Heating System- Supply Temperature Optimization and Self Requlation Effects, Chalmers University of Technology, 2010. [29] Gurobi Optimization Inc. Gurobi Optimizer, 2017.

Studentlitteratur; 2001. [24] R. Mirzanamadi, Thesis for the Degree of Licentiate of Engineering Ice Free Roads Using Hydronic Heating Pavement with Low Temperature, Chalmers University of Technology, Gothenburg, Sweden, 2017. [25] Meteotest. Meteonorm: Meteonorm, Global Meteorological Database. Handbook part II: Theory, version 6.1. Bern, Switzerland, 2010. [26] R. Mirzanamadi, C.-E. Hagentoft, P. Johansson, J. Johnsson, Anti-icing of road surfaces using hydronic heating pavement with low temperature, Cold Reg. Sci.

290