Analysis of the efficiency of an embedded floor heating for snow melting and freeze protection

An algorithm for modelling and numerical simulation of the three dimensional conjugate heat transfer in slab-embedded snow melting and freeze protecting systems is developed. The influence of the climate conditions on the heat transfer in the constructions is modelled according to the correspondent design method of ASHRAE. The models are verified and solved numerically via finite volume method for a road embedded hydronic installation at different piping configurations. The fluid flow parameters and temperature fields in the construction are obtained at steady state climate conditions, hardest for 95 % of the wintertime for the region. An analysis of the efficiency of the modelled configurations is implemented based on the computed snow-free area ratio of the road surfaces.


Introduction
The finite element simulation of the heat transfer in floor constructions with embedded radiant heating allows detailed information about the temperature field of the heat exchanging surfaces at the design and verification tasks [1 -6]. The heat transfer between the floor surfaces and the room environments at indoor installations is usually modelled via Robin boundary conditions at accepted room temperatures and heat transfer coefficients [1]. They are determined according to the well-known methods for design of radiant floor heating [7]. The heat transfer at the slab underside of ground floors of the buildings is computed at Dirichlet conditions using the ground temperature. The temperature drop of the liquid flows is relatively small (up to 5 K). That allows simplifications of the geometrical models, two dimensional and relatively quick numerical solutions.
Outdoor radiant floor systems are used in the recent decades for snow melting and ice protecting of ways, parks and stairs. Such hydronic installations are suitable for utilisation of waste heat of flows with relatively low temperatures [8]. The exchanged heat flow at the upper slab surfaces depends on the climatic factors as the temperature, relative humidity, cloudiness, snow intensity and wind speed [9]. The heat transfer is non-stationary and the temperature drop of the liquid flows can reach 15 -20 K. That results in non-uniform temperature field and local icings of the snow melting surfaces. The mathematical modelling is a useful approach for prediction and estimation of these processes at the design stage [1], but the boundary conditions are more complicated in comparison to the indoor installations. Numerical researches of the conjugate heat transfer in outdoor radiant floor heating systems, including the coupled snow heating and melting, water evaporation, convective and radiative heat transfer between the upper surfaces and the ambient environment have not been published to now. An algorithm for modelling and numerical simulation of these multiphysics processes is proposed in [10]. It was applied for an approximate numerical analysis of the steady state temperature and fluid flow fields at a hydronic system, integrated in a pavement way in order to validate the models [11].
The aim of the present study is a future developing and applying of this algorithm for analysis of the efficiency of floor-embedded snow melting and freeze protecting systems, taking into account the quasi-stationary heat transfer at severe climatic conditions.

Mathematical models
The conceptions for numerical simulations in [10] are expounded below to fully descriptions of the research. The geometrical model has to be three dimensional, covering the solid domain of the floor layers and the liquid domain of the pipe coils ( fig. 1). The pipe walls can be excluded from the geometry -they can be represented by a fluid/solid interface with the thermal resistance of the wall. The geometrical model is discretized by a finite volumes mesh, finer in the pipes and coarse in the solid floor layers.
The system of equations that is solved numerically in order to simulate the conjugate heat transfer in the construction contains: continuity equation, momentum equations, energy equation and k-ɛ turbulence model for the fluid domain; energy equation for the solid domains. The physical properties of the solids and the liquid in the equations above can be used as constants because of the relatively small temperature variations in the constructions.
The boundary conditions of the model are summarized in Table 1. where: = heat flux for snow (ice) heating to the melting temperature and water heating from the melting temperature to the accepted surface temperature (sensible heat flux), Wm -2 ; = heat flux for snow (ice) melting, Wm -2 ; = heat flux, exchanged by radiation and convection with the outer environment,Wm -2 ; = heat flux for water evaporation, Wm -2 ; Ar = snow-free area ratio, dimensionless. It is computed via the equivalent snow-free area Af and the total area of the surface At, m 2 . (2) The sensible heat flux is: where: =specific heat capacity of the ice, Jkg -1 K -1 ; = specific heat capacity of the water, Jkg -1 K -1 ; s = snowfall rate water equivalent with dimensions mmh -1 = lh -1 m -2 ; ta = ambient temperature, ⁰С; tf = liquid film temperature, ⁰С; ts = melting temperature, ⁰С; = water density, kgm -3 ; c1 =1000 mm/m x 3600 s/h = 3.6 x 10 6 . The liquid film temperature is usually accepted as 0.56°C [9]. Information about the snowfall rate can be obtained by the meteorological database [12].
The melting heat flux is given by the equation: where rm is the specific fusion heat of the snow, Jkg -1 . The convective and radiative heat flux, exchanged at the snow-free surface is: where: hc = coefficient of heat transfer by convection, Wm -2 K -1 ; Tf = liquid film temperature, K; TMR = mean radiant temperature of surroundings, K; σ = Stefan-Boltzmann constant (5.67 х 10 -8 Wm -2 K -4 ); ɛs = surface emissivity, dimensionless. The convection heat transfer coefficient is computed by: where: Kair =thermal conductivity of the air, Wm -1 K -1 : L = characteristic length of the slab in the wind direction, m; Pr = Prandtl number for the air; ReL= Reynolds number, obtained by the wind speed near the floor surface and L. The mean radiant temperature is equal to the ambient temperature at snowfall and cloudy sky or is obtained according to [9] at the clear sky.
The heat flux for the water evaporation at the wet surface is given by the equation: where: ρdry air = dry air density, kgm -3 ; hm = mass transfer coefficient, ms -1 ; xa = humidity ratio of ambient air, kgkg -1 ; xf = humidity ratio of the saturated ambient air at tf , kgkg -1 ; r = specific heat of water vaporization, Jkg -1 . The coefficient of the mass transfer is: where Sc= Schmidt number (Sc=0.6).
The computed heat flux according to (1) is used for the sizing of the floor installations at a design stage. The heating systems can be sized to provide enough thermal power for snow melting and freeze protection at the hardest metrological conditions for the region. However, the choice of the power is consistent with а rational cost of the investments and the radiant heating usually is not designed to operate satisfactory for 100 % of the climatic conditions.
A quasi-stationary heat transfer proceeds at the areas with temperatures, near the snow melting temperature if the heat power is not enough to melt the snow in the chilly and windy periods. When the surface temperatures decrease below the melting one, the snow melting stops and a snow cover is formed (Ar=0). Then the surface heat flux is transferred by conduction through the snow layer. That leads to an increasing of the surface temperature. The melting starts again if the melting temperature is reached. When the surface temperature exceeds the melting one, all heat fluxes in (1) participate in the heat transfer. The snow cover is stable if the surface temperature does not reach the melting temperature. These processes can be modeled by the following modification of (3), (4) and (5), using the fluxes as functions of the nodal temperatures t of the upper surface of the construction: (11) (12) where: δsnow=snow cover thickness, m; Ksnow=thermal conductivity of the snow, Wm -1 K -1 .

(13)
The snow-free area ratio Ar in (1) is defined as additional variable, computed as results of the temperatures for the nodes, belonging to the upper floor surface: Ar can be used for qualitative and quantitative assessment of the efficiency of the embedded radiant heating at the accepted conditions. The accuracy of its calculation depends on the convergence at the iterative solution of the system of equations. The higher the convergence is, the more accurately Ar is.

Object of investigation
The proposed models are applied to investigate the efficiency of a design solution of outdoor hydronic radiant heating. It is being embedded in a roadway near to Sofia. Two variants of  . 2). The construction of the road is consisted of four layers (Table 2).  The liquid is a solution of propylene glycol (35%). The temperature and velocity on the inlet are 45 ºC and 0.45 ms -1 . The temperature drop at one loop is accepted as 15 K. The gauge pressure on the loop outlet is 5 kPa. The system is sized to ensure 95 % satisfactory operation at the conditions in Table 3. More detail information is not possible due to confidential rules. Heat flux for water evaporation Wm -2 38 Total heat flux Wm -2 235

Numerical simulations and results
Initially the models ware solved numerically at a constant heat flux in Table 3 in order to validate them and to precise the finite volumes mesh [11]. Results of numerical solutions using the heat fluxes according to eq. (10 -14) at the same meshes are presented below. The iterative procedures were stopped at achieving of imbalances 0.2 % and 2 % at the mass and energy balances respectively. Insignificant changes of the computed variables and residues are established below these limits. The minimal temperatures on the upper surfaces are almost equal at the investigated coils. The maximal temperatures are in the areas above the first loops of the pipes along the flow path. The difference between them at the investigated variants are 2K ( fig. 3 and 4). The temperature fields in the vertical cross sections are similar. The differences between the obtained at the simulations and the computed at the design stage pressure drops are below 3%. That also proves the adequacy of the models.   The average temperatures are almost equal to the accepted surface temperature at the design stage tf =0.56 °C. However, the snow-free surface ratios are less than one due to the non-uniform temperature fields on the upper surface.
Although the snow-free surface ratio is smaller at the case of the meander coil, the locations with the possible snowing and icing are on the road periphery that is not loaded by car traffic. Such dangerous locations are in the centre of the road at the spiral coil -the icing in these places is more inacceptable.

Conclusions
The proposed algorithm for modelling and numerical simulation of the heat transfer in outdoor embedded hydraulic radiant heating for melting of snow and ice provides relatively accurate and detailed information about the efficiency of the constructions. It allows assessment of the temperature fields and the possibilities for local snowing and icing at different inlet flow parameters and configurations of the pipe coils and slabs. The models can be further developed in order to simulate numerically the non-stationary heat transfer under varying climatic conditions. The algorithm was used to investigate the efficiency of two variants of pipe coils at an outdoor radiant heating. The results of the numerical studies show an advantage of the meandering coil, as the places where it is possible to form ice and snow cover on the heated road section are safer in comparison to the variant with a spirally wound coil.