Analytic investigations of CNFP-based self-deicing road system on the deicing performance

Cold Regions Science and Technology 103 (2014) 123–132

Contents lists available at ScienceDirect

Cold Regions Science and Technology journal homepage: www.elsevier.com/locate/coldregions

Analytic investigations of CNFP-based self-deicing road system on the deicing performance Hui Li ⁎, Qiangqiang Zhang, Huigang Xiao School of Civil Engineering, Harbin Institute of Technology, Harbin 150090, China

a r t i c l e

i n f o

Article history: Received 17 October 2013 Accepted 4 April 2014 Available online 13 April 2014 Keywords: Deicing Orthogonal expansion Transient analytic solutions Phase change

a b s t r a c t A novel road deicing system consisting of a carbon nano-fiber polymer (CNFP) thermal source, an AlN ceramic insulated encapsulation layer, a multi-wall carbon nanotube (MWCNT)/cement-based thermal conductive layer and a thermal insulated substrate was developed in a previous experimental study. Based on the basic transient heat conduction theory, a mathematical model of composite-media thermal conduction is proposed in this study. Utilizing the orthogonal expansion technique, the non-homogeneous problem is split into the superposition of two steady-state non-homogeneous problems and a homogeneous transient problem; the transient analytic solutions are found in the stage without a phase change. Considering the quasi-steady hypothesis, the other parts of the solutions are determined for the phase change stage. The parameter analyses of the analytic solutions obtained in terms of the time-dependent temperature field reveal the same parameter-dependent influence and curve tendency in the deicing process as in the previous experimental study. To verify the credibility and reliability of the analytic solutions, the results are experimentally confirmed; both the theoretical and experimental approaches present similar trends except for a few slight, acceptable differences. These analytic solutions will be applied for prediction, control and guidance in further deicing studies. © 2014 Elsevier B.V. All rights reserved.

1. Introduction It is well known that there are three states of water: liquid (water), solid (ice and/or snow) and gas (vapor). The transition temperature between the solid phase (ice/snow) to the liquid state (water) is 0 °C, while it is 100 °C for transition from the liquid state (water) to the gas state (water vapor). Due to the phase change, it is difficult to describe these problems exactly, especially for ice/snow melting to water, which is associated with the formation of a moving interface. In the 18th century, research on heat by observation and limited experiments began to be conducted, and in 1822, Fourier published his law, leading to the establishment of the mathematical theory of heat conduction (Fourier, 1822). In 1891, Stefan investigated ice formation and made some progress on this problem. In 1912, F. Neumann presented the exact study of some idealized phase change problems in his notes (Carslaw and Jaeger, 1959; Ozisik, 1993). Since then, this subject has entered a period of rapid development, and many scholars in various fields, such as food thawing and sea ice melting, have studied the ice melting problem associated with heat conduction and phase change with both experiments and numerical simulations (Leng et al., 2005; Perovich et al., 2008; Zeng and Faghri, 1994). These scholars have

⁎ Corresponding author. E-mail addresses: [email protected] (H. Li), [email protected] (Q. Zhang), [email protected] (H. Xiao).

http://dx.doi.org/10.1016/j.coldregions.2014.04.001 0165-232X/© 2014 Elsevier B.V. All rights reserved.

mainly focused on the heat conduction of a single medium and its phase change process, and the only analytic solution yet available is that for a simple and homogeneous scenario. Based on the orthogonal expansion technique (Ozisik, 1993), the transient solution of the composite-media heat conduction problem, with homogeneous boundary conditions and no phase change, can be determined. However, the road deicing process in this study is a coupling problem involving composite-media heat conduction and phase change interference with non-homogeneous conditions, and the existing solution does not cover this case. Therefore, the deicing process of our system is complex; for example, most researchers only utilized experiments and/or numerical simulations to investigate the deicing and snow melting process in certain cases (Feng and Li, 2009; Yang et al., 2012; Yehia et al., 2000; Zhao et al., 2010; Zhao et al., 2011; Zhao et al., 2011), which cannot cover all situations. A carbon nanofiber polymer (CNFP)-based self-deicing road system with high deicing efficiency has been proposed (Li et al., 2013). The whole system consists of four units: the pavements unit embedded with CNFP thermal source, the monitoring unit with various sensors (thermal-couples et al.), the energy supply unit (solar energy harvest facilities and/or national grid) and the computer-based control unit. The efficiency of this system for deicing and field snow-melting has been experimentally investigated under various ambient temperatures, heat flux densities and ice thicknesses. The high efficiency, repeatability, cost and feasibility of the proposed self-deicing road system have been validated by the experiments.

124

H. Li et al. / Cold Regions Science and Technology 103 (2014) 123–132

Fig. 1. Configuration of the CNFP-based deicing system.

During the deicing process, the temperature, ice thickness and other variables collected by the monitoring unit will be feedback to the computer-based control unit, and then the input energy will be optimally determined with minimizing comprehensive index β(β = (T/Tmax) + (Ec/Ec,max)), which is defined as a comprehensive index β and accounts for the trade-off between the optimal time and energy consumption, where T and Ec are the time and energy taken for deicing, respectively, and Tmax and Ec,max are the longest time and maximum energy required for deicing a certain amount of ice, respectively (Li et al., 2013). As is well known, a numerically based control program typically requires too much time and computer capacity. Therefore, a simple and feasible analytic model of the deicing system that considers multiple factors (wind speed, heat flux, ambient temperature and ice thickness) is necessary for prediction with high precision and proper guidance of further deicing investigations and applications. In this study, the analytic solutions of deicing performance are developed, and the time-dependent temperature field is calculated and compared with the experimental results to verify its feasibility and credibility (Li et al., 2013).

Due to the excellent thermal conductivity and limited thickness of the AlN wafer (N175 W/m·K, 0.5 mm) compared to the other materials listed in Table 1, the AlN wafer can be considered to be a perfect thermal conductor, and its participation in thermal conduction modeling can be ignored except to conduct source heat transition perfectly. Therefore, the coordinate system would be established with the bottom of the thermal conductive layer as the origin and the positive direction of heat flux flow as described in Fig. 2. Therefore, this problem can be simplified as a heat conduction problem with two composite layers; the mathematical governing equations of heat-conduction in Cartesian coordinates are given in the following forms:

2. Analytic modeling of a CNFP-based deicing process

where T1(x, t) and T2(x, t) are the temperature distribution of the thermal conductive layer and the ice layer, respectively, both of which are timeand space-dependent functions, and α1 and α2 are the thermal diffusivity of the thermally conductive layer and the ice layer, respectively, which indicate the heat diffusion capacities of the materials. The second kind of boundary condition is the stable heat flux generated by the CNFP-based thermal source and applied at the bottom surface of thermal conductive layer.

2.1. One-dimensional heat conduction problem before phase change The same configuration of CNFP-based deicing system proposed in Li et al. (2013) is employed herein. As shown in Fig. 1, the system consists of four self-governed and mutual-coordinated functional units: a CNFPbased heating source to generate a stable heat flux for system; a MWCNT (multi-walled carbon nano-tube)/cement-based thermal conductive layer; An AlN (Aluminum nitride) ceramic wafer electroinsulated layer and an epoxy thermal insulated substrate. The material parameters of these four functional layers are listed in Table 1. The CNFP-based thermal source and the thermal conductive layer play very critical role in the deicing system. The heat generated by the CNFP flows through the thermally conductive layer into the ice layer one-dimensionally, so the system can be modeled as a one-dimensional heat conduction in a finite laminated slab with non-homogeneous boundary conditions, as shown in Fig. 2. The CNFP-based thermal source at the location of x = 0 will be considered as the second kind of boundary condition. And the temperature and heat flow are continuous at the interface of x = a, based on the assumption of complete thermal contact. Additionally, the system also disperses heat into the ambient by heat convection, which is treated as a convection coefficient h, applied at the upper surface of the ice layer as the third kind of boundary condition.

Ice Water AlN wafer MWCNT/cementbased composites Epoxy

Density

Thermal conductivity

Specific heat

Heat diffusivity

(kg/m3)

(W/m∙K)

(kJ/(kg∙K))

(m2/s × 10‐6)

970 1000 3235 2230

2.22 0.58 N175 2.83

2.1 4.2 0.75 0.88

1.09 0.50 N72.13 1.44

1.1–1.4

0.1–0.13

1400

∂2 T 2 ðx; t Þ 1 ∂T 2 ¼ α 2 ∂t ∂x2

−k1

0.2

∂T 1 ¼ qw ∂x

0 b x b a; t N 0;

a b x b b; t N 0;

x ¼ 0; t N 0

ð1aÞ

ð1bÞ

ð2aÞ

The continuous conditions at the interface between the ice layer and the thermal conductive layer require both heat flux and temperature continuity. k1

∂T 1 ∂T ¼ k2 2 ∂x ∂x

T1 ¼ T2

x ¼ a; t N 0;

x ¼ a; t N 0:

ð2bÞ

ð2cÞ

The third kind of boundary condition on the surface of the ice layer is the heat transfer by means of the convection coefficient h. −k2

Table 1 The material parameters (Li et al., 2013). Parameter

∂2 T 1 ðx; t Þ 1 ∂T 1 ¼ α 1 ∂t ∂x2

∂T 2 ¼ hðT 2 −T f Þ x ¼ b; t N 0 ∂x

ð2dÞ

The deicing system is subjected to the following initial condition, given as a uniform temperature that is coincident with the ambient temperature based on neglecting the slight temperature gradient in the deicing systems: T1 ¼ T f

t ¼ 0; 0 ≤ x ≤ a

ð3aÞ

T2 ¼ T f

t ¼ 0; a ≤ x ≤ b

ð3bÞ

H. Li et al. / Cold Regions Science and Technology 103 (2014) 123–132

125

Fig. 2. The heat conduction schematic of the analytic model of deicing as a one-dimensional non-homogeneous problem in a finite laminated slab.

(2) Similarly, the functions φi(x), i = 1, 2 satisfy the other steadystate heat conduction problem given as follows:

where k1 and k2 are the thermal conductivities of the thermal conductive layer and the ice layer, respectively; Tf is the ambient temperature; and h is the heat convection coefficient for heat transfer between the upper surface of the ice layer and the environment. Based on the orthogonal expansion technique (Evans et al., 2001; Mikhailov and Ozisik, 1984), it is easier to solve a problem with homogeneous boundary conditions than one with non-homogeneous boundary conditions. The time-dependent heat conduction for two heat transfer mediums with non-heat generation and nonhomogeneous conditions can be split mathematically into the superposition of a transient problem with homogeneous conditions and two other steady state cases with one non-homogeneous boundary condition each (Evans et al., 2001). Then, the formula for T i (x, t) can be considered the superposition of three simpler problems in the following form: T i ðx; t Þ ¼ θi ðx; t Þ þ ϕi ðxÞ þ φi ðxÞT f

i ¼ 1; 2:

a b x b b;

ð5bÞ

ð6bÞ

k2

dϕ1 ðxÞ ¼ qw ∂x

x ¼ a;

ð6dÞ

x ¼ a;

dφ2 þ hφ2 ¼ h dx

ð6eÞ

x ¼ b:

ð6fÞ

∂2 θ1 ðx; t Þ 1 ∂θ1 ¼ 0 b x b a; t N 0 α 1 ∂t ∂x2

ð7aÞ

∂2 θ2 ðx; t Þ 1 ∂θ2 ¼ a b x b b; t N 0 α 2 ∂t ∂x2

ð7bÞ

x¼0 ;

ð5cÞ

∂θ1 ¼0 ∂x k1

x ¼ a;

dϕ2 þ hϕ2 ¼ 0 dx

ð7cÞ

x ¼ a; t N 0;

ð7dÞ

x ¼ a; t N 0;

ð7eÞ

ð5dÞ k2

x ¼ a;

x ¼ 0; t N 0;

∂θ1 ∂θ ¼ k2 2 ∂x ∂x

θ1 ¼ θ2 dϕ dϕ k1 1 ¼ k2 2 dx dx

k2

dφ1 dφ ¼ k2 2 dx dx

ð6cÞ

The corresponding boundary and initial conditions are given as

subject to the boundary conditions:

ϕ1 ðxÞ ¼ ϕ2 ðxÞ

x¼0 ;

(3) The functions θi(x, t), i = 1, 2 satisfy the following transient heat conduction problem given in the forms:

2

−k1

a b x b b;

φ1 ðxÞ ¼ φ2 ðxÞ

2

d ϕ2 ðxÞ ¼0 dx2

d2 φ2 ðxÞ ¼0 dx2

k1

(1) The functions ϕi(x), i = 1, 2 satisfy the steady-state heat conduction problem, given as follows: ð5aÞ

ð6aÞ

dφ1 ðxÞ ¼0 ∂x

ð4Þ

0 b x b a;

0 b x b a;

subject to the boundary conditions:

Physically, ϕi(x) indicates the total energy input associated with one boundary condition of heat flux, as shown in Eq. (5c), φi(x) is the total energy dispersion related to heat convection with the environment on the surface of the overlay ice, as shown in Eq. (6f), and θi(x, t) is the transient homogeneous term for homogenous boundary conditions.

d ϕ1 ðxÞ ¼0 dx2

d2 φ1 ðxÞ ¼0 dx2

ð5eÞ

∂θ2 þ hθ2 ¼ 0 ∂x

x ¼ b; t N 0;

ð7fÞ 

t ¼ 0; 0 ≤ x ≤ a:

ð7gÞ



t ¼ 0; a ≤ x ≤ b:

ð7hÞ

θ1 ðx; t Þ ¼ T f −ϕ1 ðxÞ−φ1 ðxÞT f ¼ F 1 ðxÞ x ¼ b:

ð5fÞ

θ2 ðx; t Þ ¼ T f −ϕ2 ðxÞ−φ2 ðxÞT f ¼ F 2 ðxÞ

126

H. Li et al. / Cold Regions Science and Technology 103 (2014) 123–132

As second-order ordinary differential equations with one nonhomogeneous boundary condition, the substitution of boundary conditions Eqs. (5c)–(5f) yields the solution of the functions ϕi(x), i = 1, 2 in the form:  3 2 qw qw qw b 1 1

 þ þ aq − − w 6 k ϕ1 ðxÞ k2 k1 k2 7 h 1 7 x : ð8aÞ ¼6 4 qw 5 1 qw qw b ϕ2 ðxÞ þ − k2 k2 h Similarly, the solutions of equations for functions φi(x), i = 1, 2 can be given as follows:


 φ1 ðxÞ 0 1 x ¼ : ð8bÞ φ2 ðxÞ 0 1 1 Because these equations are eigenvalue problems, it is possible to use the method of variable separation and the orthogonal expansion technique (Mikhailov and Ozisik, 1984), and the solution of equations θi(x, t), i = 1, 2 can be given as   ∞ X λn −λ2 t x 0 b x b a; t N 0; ð9aÞ cn e n cos pffiffiffiffiffiffi θ1 ðx; t Þ ¼ α1 n¼1

θ2 ðx; t Þ ¼

∞ X n¼1

2

−λn t

cn e


    λn λn x þ B2n cos pffiffiffiffiffiffi x a b x b b; t N 0 A2n sin pffiffiffiffiffiffi α2 α2 ð9bÞ

where cn is the coefficient of combination determined from the initial conditions and orthogonality of the eigenfunctions, An2 and Bn2 are the coefficients of the eigenvalue problems, and λn is the eigenvalue determined by the transcendental Eq. (13). We denote that rffiffiffiffiffiffi aλn bλn bh k α2 γ ¼ pffiffiffiffiffiffi ; η ¼ pffiffiffiffiffiffi ;H ¼ ;κ ¼ 1 k2 α1 α2 k2 α 1 2 3 a a
 ð10Þ A2n 6 cosγ sin b η−κ sinγ cos b η 7 ¼4 5; a a B2n cosγ cos η þ κ sinγ sin η b b
pffiffiffiffiffiffi  α1 k1 k sinð2γÞ þ a þ 2 ½A þ B þ C  2α 1 2λn α2  pffiffiffiffiffi 2
a A a  A ¼ 2n ðb−aÞ− 2 sin2η− sin2 η ; b 2 2λn  pffiffiffiffiffi 2
a2 B a  B ¼ 2n ðb−aÞ þ sin2η− sin2 η 2λn b 2 pffiffiffiffiffih i a a C ¼ A2n B2n 2 cos2 η− cos2η : b 2λn Nn ¼

The eigenvalues λn are determined by the following transcendental equation:  a H−η tan 1− η b : κ tanγ ¼  a H tan 1− η þ η b

Then, the temperature distribution Ti(x, t), i = 1, 2 can be written completely in the forms: T 1 ðx; t Þ ¼ θ1 ðx; t Þ þ ϕ1 ðxÞ þ φ1 ðxÞT f      ∞ X λn 1 1 b 1 1 −λ2 t x −qw ¼ cn e n cos pffiffiffiffiffiffi x− − −a − k1 h k2 k1 k2 α1 n¼1 þ T f 0 b x b a; t N 0

T 2 ðx; t Þ ¼ θ2 ðx; t Þ þ ϕ2 ðxÞ þ φ2 ðxÞT f
    ∞ X λn λn −λ2 t x þ B2n cos pffiffiffiffiffiffi x ¼ cn e n A2n sin pffiffiffiffiffiffi α2 α2 n¼1   1 1 b −qw x− − þ T f 0 b x b a; t N 0 k2 h k2

Based on the orthogonality of the eigenvalue problem (Evans et al., 2001; Ozisik, 1993), cn can be determined as
 1 k1 k G1 þ 2 ðG2 þ G3 Þ cn ¼ Nn α1 α2     pffiffiffiffiffiffi

α1 a 1 b 1 1 λn G1 ¼ qw a − − −a − sin pffiffiffiffiffiffi λn k1 h k2 k1 k2 α1  pffiffiffiffiffiffi   pffiffiffiffiffiffi α1 α1 λn þ cos pffiffiffiffiffiffi a − α1 k1 λn k1 λn    pffiffiffiffiffiffi
pffiffiffiffiffiffi
 α2 α2 λn λn b − sin pffiffiffiffiffiffi a G2 ¼ qw A2n sin pffiffiffiffiffiffi ð12Þ α2 α2 λn k 2 λn       1 λn a 1 b λn þ cos pffiffiffiffiffiffi b þ a − − cos pffiffiffiffiffiffi h k2 h k2 α2 α2     pffiffiffiffiffiffi
pffiffiffiffiffiffi
α2 α2 λn λn G2 ¼ qw b − cos pffiffiffiffiffiffi a B2n cos pffiffiffiffiffiffi α2 α2 λn k 2 λn       1 λn 1 b a λn − cos pffiffiffiffiffiffi b þ þ − a : sin pffiffiffiffiffiffi h h k2 k2 α2 α2

ð14aÞ

ð14bÞ

The time- and space-dependent temperature distributions in the ice layer and the thermal conductive layer can be obtained with Eq. (14a) and Eq. (14b). For the deicing procedure described by this problem, the functions are valid for non-phase change states. Therefore, the solutions given in Eq. (14a) and Eq. (14b) are only applied to determine the temperature distributions of the composite mediums and critical points T ðx; t Þjt¼t0 ;x¼a ¼ 0 before the ice begins to melt. The stage of the phase change will be discussed with mathematical transformations. According to the numerical study of our FEM program, after the critical time t0, the temperature gap between the top and bottom surface of the ice layer is normally far less than 10 °C until the phase change is complete. The Stefan number is verified to be smaller than 0.1, as presented in Eq. (15), which implies that the temperature distribution can be considered a quasi-steady procedure with a phase change occurring (Jiji, 2009). Ste ¼

ð11Þ

ð13Þ

cpi ðT b −T a Þ 2:1  103  10 b ¼ 0:062 b 0:1; L 335  103

ð15Þ

where Ste is the Stefan number, Tb and Ta are the temperatures of the top and bottom surfaces of the ice layer, cpi is the specific heat of ice, and L is the latent heat of ice (335 kJ/kg). The heat flux qw is the inflow energy of the ice layer, and the heat flux at the upper surface of the ice layer, denoted qh at the time t0, is the outflow energy. Both of these two kinds of energy flows are applied to the ice layer. However, to study the phase change stage, a quasisteady boundary condition based on energy conservation is used. ∂T 2 ðx; t Þ jx¼b;t¼t 0 ∂x
    ∞ X cn λn −λ2n t λn λn b −B2n sin pffiffiffiffiffiffi b þ qw ; ð16Þ ¼ −k2 A2n cos pffiffiffiffiffiffi pffiffiffiffiffiffi e α2 α2 α2 n¼1

qh ¼ −k2

where qh is the equivalent heat flux applied to the top surface of the ice layer. 2.2. One-dimensional heat conduction problem with phase change When the temperature at the interface between the thermally conductive layer and the ice layer increases up to 0 °C the ice begins a phase change and melts gradually. In the analytic model, the solution of such problems is inherently difficult due to the formation of a moving

H. Li et al. / Cold Regions Science and Technology 103 (2014) 123–132

interface between the melted water and the remaining ice. Therefore, the location of the liquid–solid interface is unknown a priori and must follow as a part of the solution, where the latent heat is absorbed and the temperature field is discontinuous. Based on the quasi-steady simplification, this study focuses on the phase change and the moving interface of the ice layer. The problem is thus given with following mathematical formulas. As shown in the model in Fig. 3, the temperature distribution of the liquid–solid two-phase region is governed by an equation in the forms: ∂T w 2 ðξ; t Þ 2

∂ξ

∂T i 2 ðξ; t Þ ∂ξ2

¼ 0 0 b ξ b sðt Þ; t N 0;

¼0

ð17aÞ

sðt Þ b ξ b δ; t N 0;

ð17bÞ

where δ is the thickness of the initial ice and is equal to b − a, as mentioned in Section 2.1. The melting begins at the boundary ξ = 0, the liquid–solid interface ξ = s(t) moves in the positive ξ direction, and s(t) is the location of the liquid–solid interface, which is not known a priori and is thus determined as a part of the solution. The subscripts w and i refer to water and ice, respectively. Hence, the problem involves three evaluated functions, Tw(ξ, t), Ti(ξ, t) and s(t), but Eq. (17a) and Eq. (17b) only give two governing equations. Therefore, another equation is necessary and is defined by considering the energy conservation at the interface ξ = s(t) in the form: kw

∂T w ðξ; t Þ ∂T ðξ; t Þ dsðt Þ −ki i ¼ ρw L dt ∂ξ ∂ξ

ξ ¼ sðt Þ; t N 0;

ð18aÞ

where L is the latent heat of ice per unit mass (kJ/kg) associated with the phase change (335 kJ/kg), kw and ki are the thermal conductivities of water and ice, respectively, and ρw is the density of water. Additionally, the temperature at the liquid–solid interface is stable at the phase change value Tm given by T w ðξ; t Þ ¼ T i ðξ; t Þ ¼ T m

ξ ¼ sðt Þ; t N 0

ð18bÞ

Eqs. (17) and (18a) provide three differential equations that determine the temperature distributions in water and ice and the location of the phase change interface. The outer surfaces are subjected the stable heat flux at ξ = 0 and heat convection at ξ = δ in the forms: ∂T w ðξ; t Þ q ¼−w k w ∂ξ

∂T i ðξ; t Þ q ¼− h k i ∂ξ

ξ ¼ 0;

ξ ¼ δ:

Fig. 3. The schematic of the ice phase change.

ð18cÞ

ð18dÞ

127

The solutions of these equations are given as T w ðξ; t Þ ¼ −

T i ðξ; t Þ ¼ −

sðt Þ ¼

qw q ξ þ T m þ w sðt Þ kw kw

t N 0; 0 b ξ b sðt Þ;

qh q ξ þ T m þ h sðt Þ t N 0; sðt Þ b ξ b δ; ki ki

ðqw −qh Þ t: ρw L

ð19aÞ

ð19bÞ

ð19cÞ

According to Eq. (19c), when the location of the moving interface reaches the upper surface of the ice layer (s(t) = b − a), the time required for the phase change of the ice layer can be calculated as tp ¼

ρw L ðb−aÞ: ðqw −qh Þ

ð19dÞ

Then, the time necessary for deicing requirements is given as ta ¼ tp þ t0 ;

ð20Þ

where tp is the time used in the stage of phase change, and ta is the total time used during the deicing process. 3. Results and discussions 3.1. The validation of temperature field between analytic and experimental study Based on the analytic solutions developed, the temperature distributions of the deicing system can be determined. To verify the reliability of the analytic solutions, the parameters are set equal to those used in the experiment: the thickness of the thermally conductive layer was b = 25 mm, the ice layer is 10, 15 and 20 mm thick, the heat flux is 1000, 1400 and 1800 W/m2, and the convection coefficient is 10, 20 and 30 W/(Km2); the other parameters of the materials are listed in Table 1. The time and energy consumed in deicing are important. The time utilizes to finish deicing a certain amount of ice under certain conditions with a fixed heat flux can be found on the curve of the temperature distribution, and then the energy can be calculated with this time. As shown in Fig. 4, with an ambient temperature of −10 °C and a heat flux of 1800 W/m2, for three thicknesses of the ice layer, the time-dependent temperature field of the experimental results and the analytic study at the moving interface present similar trends. The temperature increases dramatically with the solid state of ice at the beginning (b0 °C). There is a sudden change in the temperature gradient when the temperature reaches the melting point Tm = 0 °C, and the curve flattens and remains at 0 °C for a while. This phenomenon implies a phase change and latent heat absorption at the moving interface. The extent of the flattened area signifies the time for the stage of phase change. When the ice melts to water completely, the temperature increases again, similar to that at the beginning stage (N0 °C). Although the bottom surface of the experimental set-up is insulated as perfectly as possible, there is still a small heat flow downwards that dispersed into ambient, which implies that the actual heat flux flowing upwards was less than 1800 W/m2 in the experiment. Therefore, it is not difficult to understand why the results of the experiment had a smaller curve slope and a longer deicing time than in the analytic study. However, from the acceptable difference of the curve slope, the analytic model cannot only be modified by the effective heat flux flowing upwards, but must also consider the effect of extra heat consumption, which is also verified with good agreement of the analytic and experimental studies (Li et al., 2013).

128

H. Li et al. / Cold Regions Science and Technology 103 (2014) 123–132

a)

10

Temperature (oC)

5 0 -5

1000W/m2 1400W/m2 1800W/m2

-10 0

1000

2000

3000

4000

5000

Time (sec)

b) 8000 -10 oC

Time (sec)

7000

-20 oC -30 oC

6000 5000 4000 3000 2000

1000

1200

1400

1600

1800

Heat flux density (W/m2) Fig. 5. The effect of the heat flux on the deicing process for 10 mm-thick ice under a 10 W/(Km2) thermal convection: (a) the time-dependent temperature field at the moving interface under a −10 °C ambient temperature, and (b) the time consumption to deice a certain amount of ice.

Fig. 4. The comparison of the time-dependent temperature field between the experiment and analytic solution at the moving interface between melted water and the remaining ice at an ambient temperature of −10 °C with a heat flux of 1800 W/m2 for (a) 10 mm, (b) 15 mm, and (c) 20 mm thick ice.

3.2. Investigation of parameters influencing deicing time The time-dependent temperature field comparisons between the experimental results and analytic solutions verify the credibility of analytic solutions. To further understand the relationship of the parameters and the deicing time and to supply guidance for future deicing applications and develop an intelligent controlling model, it is necessary to investigate how the parameters influenced the deicing time.

3.2.1. Effect of the heat flux As shown in Fig. 5(a), for 10 mm-thick ice under −10 °C ambient temperature and 10 W/(Km2) thermal convection, the time-dependent temperature of the moving interface increases dramatically at the beginning for 1000, 1400 and 1800 W/m2 heat flux, and then there is a sudden

change of the temperature gradient at 540, 660 and 840 s, respectively, which demonstrates the considerable effect of latent heat on the temperature field in the deicing process and indicates the beginning of the phase change. After approximate 2460, 3180 and 4560 s, respectively, the phase changes are complete, and the temperature increases dramatically again with the same rate as the initial heat stage. The higher heat flux means that the more thermal energy flows into system per unit time and leads to more dramatic change of temperature, the earlier the beginning of the phase change and the less time required for deicing a certain amount of ice. Based on the solution described in Eq. (20), Fig. 5(b) presents the time expenditures of three different ambient temperature cases in the deicing process, which decrease with increasing heat flux as a hyperbolic decreasing function in the first quadrant, which are in agreement with the experimental study paper (Li et al., 2013). It is also found that the time expenditure for the deicing process decreases dramatically with an increasing lower heat flux, but the effect is decreased for continuous increases. 3.2.2. Effect of the ambient temperature As shown in Fig. 6(a), the lower the ambient temperature, the later the phase change begins and the more time is consumed in the deicing process. The curve transitions are similar to those of the previous cases. The duration of the flat stage is slightly different, but the slopes of the three cases' curves are the same, which implies the same temperature gradient due to the same heat flux. For 10 mm-thick ice with a 1000 W/m2 heat flux under 10 W/(Km2) thermal convection, there are approximately 840, 1860 and 2640 s before the phase change and approximately 4860, 5820 and 7320 s for complete deicing for three different ambient temperature cases, −10, −20 and −30 °C, respectively.

H. Li et al. / Cold Regions Science and Technology 103 (2014) 123–132

a)

a)

20

0 -10 -10 oC

-20

10

5

Temperature (oC)

Temperature (oC)

10

129

0

-5 10mm 15mm 20mm

-20 oC -30 oC

-30 0

2000

4000

6000

-10 0

8000

2000

b)

6000

1800W/m

10000

2 2

7000

1400W/m 1000W/m 2

9000

Time (sec)

6000 5000 4000

8000 7000 6000 1000W/(K*m2) 1400W/(K*m2) 1800W/(K*m2)

5000

3000

4000 2000

8000

b)

8000

Time (sec)

4000

Time (sec)

Time (sec)

-30

-25

-20

-15

10

-10

12

Fig. 6. The effect of the ambient temperature on the deicing process for 10 mm-thick ice with a 10 W/(Km2) thermal convection: (a) the time-dependent temperature field at the moving interface with a 1000 W/m2 heat flux, and (b) the time consumption to deice a certain amount of ice.

As shown in Eq. (17a) and Eq. (17b), the ambient temperature linearly affects the time-dependent temperature field. The energy is used for heating the system and the ice layer before the phase change, which is mainly determined by the specific heat and the temperature change because the specific heat is a constant material parameter, and the temperature change is associated with the ambient temperature as the initial state of the system. Therefore, the time consumption before the phase change stage is linearly affected by the ambient temperature. For the phase change stage, latent heat is uncorrelated with the temperature but is governed by the quality of ice, so the total time consumed in deicing is linearly affected by the ambient temperature under identical other parameters. Fig. 6(b) clearly shows that that the time consumed for deicing for the three heat flux densities decreases linearly with the increase in the ambient temperature, which is also in agreement with the experiment study paper (Li et al., 2013). 3.2.3. Effect of ice thickness To avoid confusion, we emphasize that the whole deicing process can be divided into the heating stage and the phase change stage. The energy cost in the heating stage is mainly dominated by the materials of the system and the ice from the initial state up to the melting point with the specific heat of the materials and the ambient temperature, while the phase change stage is determined by the quantity of ice due to the effect of latent heat. For the unit amount of ice denoted with the ice thickness here, the latent heat is 16 times greater than the specific heat consumed in a 10 °C temperature change. Therefore, the effect of the ice thickness is more significant in the phase change stage than in the stage before the

14

16

18

20

Ice thickness (mm)

Temperature (oC)

Fig. 7. The effect of ice thickness on the deicing process for a −10 °C ambient temperature and 10 W/(Km2) thermal convection: (a) the time-dependent temperature field at the moving interface with a 1000 W/m2 heat flux, and (b) the time consumption to deice a certain amount of ice.

phase change. Fig. 7(a) verifies this finding with approximately 840, 960 and 1080 s to start the phase change and approximately 4560, 6420 and 8160 s to finish deicing for three ice thickness cases, 10, 15 and 20 mm, with a 1000 W/m2 heat flux under a −10 °C ambient temperature and 10 W/(Km2) thermal convection. It is clear that the ice thickness is a linear function of the latent heat, and thus, the time consumed in the deicing process is an approximately linear function of the ice thickness due to the large difference between the specific heat and latent heat of ice. However, the nonlinear aspect becomes obvious with an increase in the ambient temperature. These effects are presented in Eq. (7b), Eq. (14a) and Eq. (14b), which is also verified in the experimental study paper (Li et al., 2013). 3.2.4. Effect of the wind speed Wind at the ice surface enhances heat transfer and carries away heat with the fluid flow, which is characterized by a convection coefficient in the theory of convection heat transfer as the third kind of boundary condition. Therefore, the study of the effect of the wind speed on deicing is, in fact, the effect of the convection coefficient, which can be expressed with the following power equation (Karlekar and Desmond, 1977): h¼

i λ λh n m k N ¼ C Pr Gr Re l u l

ð21aÞ

where Nu is the Nusselt number, which depends on the intensity of the convection heat transfer between ice surface and air; Pr is the Prandtl number, which indicates the properties of the fluid and is approximately 0.7 for air; Re is the Reynolds number, which is used to indicate the

H. Li et al. / Cold Regions Science and Technology 103 (2014) 123–132

fluid states (laminar, turbulence); and Gr is Grashof number, which is used to determine the effect of natural convection on heat transfer. The convection coefficient is a function of many parameters (wind speed, space size, air density, air viscosity, ambient temperature, etc.) (Wen et al., 1987). In the previous experimental study, the wind speed is much lower than 20 m/s, the Reynolds number Remax b 1.38e5 b 5e5. For convection heat transfer over an infinite plate, the natural convection can be ignored due to the strengthening of forced convection for N 104 (Wang, 1980). Then, the flow state present as laminar state, the Nusselt number should be calculated by the format of Eq. (21b) (Eduardo, 2010).

a)

10 5

Temperature (oC)

130

0 -5

10W/(m 2 *K) 20W/(m2 *K) 30W/(m2 *K)

-10 1 1 Lh Nu ¼ ¼ 0:664Re =2 Pr =3 λ

ð21bÞ

2000

3000

4000

5000

6000

ð22Þ

where u is the wind speed and l is the length of the wind flow and is on the same order as the test model (0.1 m). As shown in Fig. 8, the effect of wind speed on the deicing efficiency is modeled with the convection coefficient in theory, and the convection coefficient increases as a power function of the wind speed, which actually means the higher the wind speed is, the more heat is dispersed. Thus, the influence of the wind speed on deicing is equivalent to the convection coefficient. As denoted in Fig. 9(a), for 10 mm-thick ice with a 1000 W/m2 heat flux, there are approximately 840, 900 and 960 s before the phase change begins and approximately 4560, 4920 and 5340 s before complete deicing for the three thermal convection states 10, 20 and 30 W/(K ∗ m2), respectively. The time-dependent temperature field at the moving interface is not sensitive to the thermal convection before the phase change stage due to the overlapping of the ice layer on the surface, but it is significantly affected by the length of the flat part of the curve associated with the phase change stage because the convection coefficient is one of most important parameters for determining the outflow heat flux at the upper surface of the ice layer at the time t0. The equivalent heat flux directly determines the time consumed at the phase change stage. The more dramatic the thermal convection, the more time is consumed in the deicing process due to the increase of the thermal convection with the environment. In particular, the time of the deicing process increases slightly with an increase in the convection coefficient at lower levels, but it is more

12 9 6 3

b)

5000 -10 oC

Temperature (oC)

u1=2 λ λ 1=3 1=2 Nu ¼  0:664  P r Re ¼ 3:778 l l l

Wind speed (m/s)

1000

Time (sec)

Thus, the convection coefficient can be simplified to a function of both the Prandtl number Pr and the Reynolds number Re and determined by the following formula: h¼

0

-20 oC

4000

-30 oC

3000 2000 1000 0

10

15

20

25

30

Convection coefficient (W/K*m2) Fig. 9. The effect of thermal convection on the deicing process for 10 mm-thick ice with a 1000 W/m2 heat flux: (a) the time-dependent temperature field at the moving interface under a −10 °C ambient temperature, and (b) the time consumption to deice a certain amount of ice.

sensitive for higher thermal convection, as shown in Fig. 9(b), which is similar to the experimental study paper (Li et al., 2013). 3.3. The prediction and control model for future deicing studies 3.3.1. The optimal evaluation model of deicing time and energy consumption Based on the principle of least time for deicing a certain amount of ice with the least energy consumption, the essential intention of the analytic investigation is to propose the mathematical model for time and energy consumed in deicing for future intelligent control and prediction of the optimal heat flux input under certain conditions. It is known that the time expenditure of deicing is a monotonic hyperbolic-decreasing function of the heat flux and is an approximately linear function of the ice thickness, ambient temperature and convection coefficient. Based on the verified relationship with the parameters, the mathematical model of the deicing time and energy are developed as a function of the proposed parameters (ice thickness, ambient temperature, heat flux and wind speed): t a ¼ t ðqw ; b−a; T f ; hÞ ¼ t p þ t 0

ð23aÞ

h i Ec ¼ Ec ðqw ; b−a; T f ; hÞ ¼ qw  t p þ t 0

ð23bÞ

0 0

10

20

30

40

h (W/(K*m2)) Fig. 8. The relationship between the wind speed and the convection coefficient.

where ta is the time required for deicing and Ec is the energy consumed in the deicing process. Utilizing the derivation of Eq. (23b), the corresponding heat flux leads to an expected minimum value of the energy expended on deicing, which is expressed in Eq. (24). Therefore, as a promising

H. Li et al. / Cold Regions Science and Technology 103 (2014) 123–132

131

development, this parameter could be employed for the prediction of further deicing applications and the optimum control of the heat flux input considering the comprehensive efficiency of time and energy.

the energy EP can also be estimated by the production of time Tp and heat flux qw (EP = tp × qw).

h i ∂Ec ¼ 0⇔Ec; min ¼ qw; prior  t p;prior þ t 0;prior ∂qw

4. Conclusions ð24Þ

3.3.2. The road ice or snow free control model Usually, to ensure the smooth traffic mostly, the anti-icing scenario of road is needed to keep ice or snow free throughout the chilled or snowy day, which therefore requires that the road would have been heated beyond the freezing point (°C) before icing or snowing. To conduct this purpose, a mathematical model should be established for PC control unit to predict the accurate and optimal time that the deicing system should be running in advance for a certain hours under certain ambient temperature conditions. The governing equation of the model is theoretically given as Eq. (25) ~ Eq. (26), which is actually the simplified format of Eq. (1) ~ Eq. (3) with only a single layer thermal conduction media (road ~ MWCNT/Cement-based thermal conductive layer), based on the similar processes to solve Eq. (1) ~ Eq. (3). The analytic solution of Eq. (25) ~ Eq. (26) can be yield as Eq. (27), which is involved in the consideration of multiple weather parameters (ambient temperature, wind speed, snow density and thickness, weather conditions and input power). And the optimal time and performance power for system running in advance could be predicted accurately based on the optimal efficiency both time and energy consumption through the calculation of Eq. (27). ∂2 T 1 ðx; t Þ 1 ∂T 1 ¼ 0 b x b a; t N 0 α 1 ∂t ∂x2

−k1

∂T 1 ¼ qw ∂x

−k1

∂T 1 ¼ hðT 1 −T f Þ x ¼ a; t N 0 ∂x

T1 ¼ T f

x ¼ 0; t N 0

t ¼ 0; 0 ≤ x ≤ a

ð25Þ

ð26aÞ

ð26bÞ

ð26cÞ

where α1 and k1 are the thermal diffusivity and thermal conductivity of the pavement concrete respectively, Tf is the environmental temperature; and V is the wind velocity flowing over the pavement surface, h1 is heat convection coefficient for heat transfer between the upper surface of the pavement and the environment, l is width of the pavement, a is thickness of the pavement, qw is heat flux density to power the thermal source, ! rffiffiffiffi pffiffiffiffi   ∞  X −4qw α −λ2n t λn 1 l d p ffiffiffi ffi x−0:265 − T 1 ðx; t Þ ¼ cos e x −qw k V k πkλn α n¼1 þ T f 0 b x b d; t N 0

ð27Þ

where λn is eigenvalue determined by the transcendental equation: qffiffiffiffiffi l pffiffinffi ¼ −0:265kλ tan dλ n αV . α According to the experimental results of outdoor deicing or snowthawing applications of heating road system (Li et al., 2013; Yang et al., 2012; Yehia et al., 2000; Zhao et al., 2010; Zhao et al., 2011), the ice or snow free road normally can be realized after the temperature of road surface rising up to and stabilizing at freezing point (0 °C), and it would be better if there is a 1 ~ 2 °C surplus. The predicted time tp for running system in advance could be calculated by Eq. (27) with Tp = 0 °C, and

Given the thermal and electrical properties of the novel material CNFP, a deicing approach was employed in a previous experimental study paper. Based on the mathematical splitting operation, the orthogonal expansion technique and a quasi-steady assumption at the phase change stage, analytic solutions are successfully developed to present the deicing performance and to predict the time-dependent field of a deicing system initially at a uniform temperature equal to the ambient temperature, which is lower than the phase change point, and gradually heated with the input of a stable heat flux under certain wind speed conditions. The whole process can be divided into three stages according to the increase in temperature: the original heating stage (b0 °C), the phase change stage (0 °C) and the continuous heating stage (N0 °C). The critical time denoting state change as a sudden change on the curve can be obtained, such as the time that deicing starts and that deicing is complete. Utilizing the parametric analysis of the analytic solutions developed, the effects of heat flux, ambient temperature, ice thickness and wind speed on the deicing of time consumption are discussed. It is found that the time consumption increases linearly with an increase in ice thickness, decreases linearly with ambient temperature and decreases as a hyperbolic function of the heat-flux density. For wind speed, time consumption is a power function convection coefficient, and the time consumption linearly increases with an increase in the convection coefficient at a small scale (b 30 W/(m2K)) but increases nonlinearly at a larger scale (N 30 W/(m2K)). Similar parameter-dependent relationships were found in the experimental study, which verified the feasibility and credibility of the analytic solutions. Furthermore, the analytic solutions offer a reasonably accurate prediction for the optimal control and operation of deicing applications, in addition to many other applications in food thawing, soil freezing and thawing and others. Acknowledgments This study was financially supported by NCET grant no. 2011BAK02B02, NSFC grant no. 50808059 and HIT NSRIF grant no. 2010019. References Carslaw, H.S., Jaeger, J.C., 1959. Conduction of heat in solids, 2nd ed. Clarendon Press, Oxford. Eduardo, C., 2010. Heat transfer in process engineering. McGraw-Hill, New York pp. 43–44. Evans, G., Blackledge, J., Yardley, P., 2001. Analytic methods for partial differential equations, 2nd ed. Springer, New York pp. 67–75. Feng, J., Li, H., 2009. Development of self-heating concrete using carbon nano-fiber paper. Proc. of SPIE, vol. 7292 (CA: San Diego). Fourier, J.B., 1822. Theorie Analytique de la Chaleur (English trans. By A. Freeman, 1955). Dover Publications, New York. Jiji, L.M., 2009. Heat conduction, 3rd ed. Scientific Publishing Services Pvt. Ltd., Chennai, India. Karlekar, B.V., Desmond, R.M., 1977. Engineering heat transfer. West Publishing Company, St. Paul, MA. Leng, M., Ching, W.H., Leung, D.Y., Lam, G.C.K., 2005. Analytic study of heat transfer with moving phase-change interface in thawing of frozen food. J. Phys. D. Appl. Phys. 38, 477–482. Li, H., Zhang, Q., Xiao, H., 2013. Self-deicing road system with a CNFP high-efficiency thermal source and MWCNT/cement-based high-thermal conductive composites. J. Cold Reg. Eng. 86, 22–35. Mikhailov, M.D., Ozisik, M.N., 1984. Unified analysis and solutions of heat and mass diffusion. Wiley, New York. Ozisik, M.N., 1993. Heat conduction, 2nd ed. Wiley, New York. Perovich, D.K., Richter-Menge, J.A., Jones, K.F., Light, B., 2008. Sunlight, water, and ice: extreme Arctic sea ice melt during the summer of 2007. Geophys. Res. Lett. 35, L11501. Wang, B., 1980. Engineering heat and mass transfer. Science PressI, Beijing, China.

132

H. Li et al. / Cold Regions Science and Technology 103 (2014) 123–132

Wen, Z., Chen, H., Dai, H., 1987. Heat transfer. Shanghai Jiaotong University Press, Shanghai, China. Yang, T., Yang, Z., Singla, M., Song, G., Li, Q., 2012. Experimental study on carbon fiber tape-based deicing technology. J. Cold Reg. Eng. 26, 55–70. Yehia, S., Tuan, C.Y., Ferdon, D., Bing, C., 2000. Conductive concrete overlay for bridge deck deicing: mixture proportioning, optimization and properties. ACI Mater. J. 97, 172–181.

Zeng, X., Faghri, A., 1994. Experimental and numerical study of microwave thawing heat transfer for food materials. ASME J. Heat Transf. 116, 446–455. Zhao, H., Wang, S., Wu, Z., Che, G., 2010. Concrete slab installed with carbon fiber heating wire for bridge deck deicing. J. Transp. Eng. 136, 500–509. Zhao, H., Wu, Z., Wang, S., Zheng, J., Che, G., 2011. Concrete pavement deicing with carbon fiber heating wires. Cold Reg. Sci. Technol. 65, 413–420.