Transient Heat Diffusion with Temperature-Dependent Conductivity and Time-Dependent Heat Transfer Coefficient

Abstract

Lie point symmetry analysis is performed for an unsteady nonlinear heat diffusion problem modeling thermal energy storage in a medium with a temperature-dependent power law thermal conductivity and subjected to a convective heat transfer to the surrounding environment at the boundary through a variable heat transfer coefficient. Large symmetry groups are admitted even for special choices of the constants appearing in the governing equation. We construct one-dimensional optimal systems for the admitted Lie algebras. Following symmetry reductions, we construct invariant solutions.

1. Introduction

For many years, considerable attention has been paid to the collection, storage, and use of thermal energy to meet various energy demands. The use of solar energy to meet the thermal demands of industries, electronics devices, and residential establishments, and so forth, is fast growing in many countries of the world [1]. Solar energy is provided by the light energy that comes from the sun. An important component of thermal systems designed for such purposes is a thermal energy storage unit. Solar collectors transform short wavelengths into long wavelengths and trap this energy in the form of heat which is transferred and transported into a heat storage vault. The medium in which the energy is stored may be fluid or solid [2]. For instance, in middle- and low-temperature solar energy systems, water and stones are the best and cheapest storing energy medium [3]. The heat energy collected by solar energy collectors increases the temperature of the medium, so the heat energy is stored in the medium. When needed, the heat energy is desorbed for use. The effectiveness of a liquid thermal storage system is determined by how the temperature of the system decays as a result of heat losses to the environment [2, 4]. The thermal energy storage problem in a medium with temperature-dependent thermal conductivity constitutes an unsteady nonlinear heat diffusion problem and the solutions in space and time may reveal the appearance of thermal decay in the system. In order to predict the occurrence of such phenomena, it is necessary to analyze a simplified mathematical model from which insight might be gleaned into an inherently complex physical mechanism.

Meanwhile, the solution of unsteady nonlinear heat diffusion equations in rectangular, cylindrical, and spherical coordinates remains a very important problem of practical relevance in the engineering sciences [5]. Recently, the ideas of hybrid analytical-numerical schemes for solving nonlinear differential equations have experienced a revival (see, e.g., [6]). One such trend is related to the combination of group theoretic approach and Adomian decomposition method [7]. This hybrid analytical-numerical approach is also extremely useful in the validation of purely numerical schemes.

In the present work, we study an unsteady nonlinear heat diffusion problem modeling thermal energy storage in a medium with power law temperature-dependent thermal conductivity and subjected to a convective heat transfer to the surrounding environment at the boundary. The mathematical formulation of the problem is established in Section two. In Section three, we introduce and apply some rudiments of Lie group techniques. In Section four, we construct the one-dimensional optimal systems, and perform reductions by one variable and construct invariant solutions in Section five. Some discussions and conclusions are presented in Section six.

2. Governing Equations

Consider an unsteady thermal storage problem in a body whose surface is subjected to heat transfer by convection to an external environment having a heat transfer coefficient that varies with respect to time. The energy equation in a rectangular, cylindrical, or spherical coordinate system can be used to find the temperature distribution through a region defined in an interval 0<𝑟<𝑎. The unsteady heat conduction problem can be described by the following governing equation [1, 6]:𝜌𝑐_𝑝𝜕𝜃=1𝜕𝑡𝑟^𝑚𝜕𝜕𝑟(𝐾(𝜃)𝑟^𝑚𝜕𝜃𝜕𝑟)−𝑄(𝜃−𝜃_∞) (2.1)with the initial condition𝜃(𝑟,𝑡)=𝜃_𝑖,at𝑡=0, (2.2)and the boundary conditions𝜕𝜃𝜕𝑟=0,at𝑟=0,𝐾(𝜃)𝜕𝜃𝜕𝑟=−ℎ(𝑡)(𝜃−𝜃_∞),at𝑟=𝑎, (2.3) where 𝜃 is the temperature, 𝐾(𝜃)=𝐾₀((𝜃−𝜃_∞)/(𝜃_𝑖−𝜃_∞))^𝑛 is the power law temperature-dependent thermal conductivity, 𝑛 is a constant, 𝑡 is time, 𝜃_𝑖 is the initial temperature of the body, 𝜃_∞ is the temperature of the environment, 𝜌 is the density, 𝑐_𝑝 is the specific heat at a constant pressure, 𝑄 is the heat loss parameter, and ℎ(𝑡)=ℎ₀𝑓(𝑡) is the time-dependent heat transfer coefficient, with ℎ₀ and 𝐾₀ being constants. The geometry of the body is specified by 𝑚=0,1,2 representing rectangular, cylindrical, and spherical coordinates, respectively. Equations (2.1)–(2.3) are made dimensionless by introducing the following quantities: 𝑟𝑟=𝑎,𝐾𝑡=⁰𝑡𝜌𝑐_𝑝𝑎²,𝜃=𝜃−𝜃_∞𝜃_𝑖−𝜃_∞𝑄,𝑠=𝜌𝑐_𝑝,Bi=𝑎ℎ₀𝐾₀. (2.4)

Neglecting the bar symbol for clarity, the dimensionless boundary value problem (BVP) becomes𝜕𝜃=1𝜕𝑡𝑟^𝑚𝜕𝜕𝑟(𝜃^𝑛𝑟^𝑚𝜕𝜃𝜕𝑟)−𝑠𝜃 (2.5)subject to𝜕𝜃𝜃𝜕𝑟=0,at𝑟=0,_𝑛𝜕𝜃𝜕𝑟=−Bi𝑓(𝑡)𝜃,at𝑟=1, (2.6) where Bi is the Biot number and 𝑠 is the heat loss parameter.

3. Classical Lie Point Symmetry Analysis

In brief, a symmetry of a differential equation is an invertible transformation of the dependent and independent variables that does not change the original differential equation. Symmetries depend continuously on a parameter and form a group; the one-parameter group of transformations. This group can be determined algorithmically. The theory and applications of Lie groups may be obtained in excellent texts such as those of [8–13]. In essence, determining symmetries for the governing equation (2.5) implies seeking transformations of the form𝑟_∗=𝑟+𝜖𝑅(𝑡,𝑟,𝜃)+𝑂(𝜖²𝑡),_∗=𝑡+𝜖𝑇(𝑡,𝑟,𝜃)+𝑂(𝜖²𝜃),_∗=𝜃+𝜖Θ(𝑡,𝑟,𝜃)+𝑂(𝜖²), (3.1) generated by the vector field𝜕Γ=𝑇(𝑡,𝑟,𝜃)𝜕𝜕𝑡+𝑅(𝑡,𝑟,𝜃)𝜕𝜕𝑟+Θ(𝑡,𝑟,𝜃)𝜕𝜃, (3.2)which leaves the governing equation invariant. Note that we seek symmetries that leave a single equation (2.5) invariant rather than the entire boundary value problem, and apply the boundary condition to the obtained invariant solutions. It is well known that the dimension of symmetry algebra admitted by the governing equation may be reduced if one seeks invariance of the entire BVP (see, e.g., [9]). The action of Γ is extended to all the derivatives appearing in the governing equation through second prolongationΓ⁽²⁾=Γ+Θ_[𝑡]𝜕𝜕𝜃_𝑡+Θ_[𝑟]𝜕𝜕𝜃_𝑟+Θ_[𝑟𝑟]𝜕𝜕𝜃_𝑟𝑟, (3.3)whereΘ_[𝑡]=𝐷_𝑡(Θ)−𝜃_𝑟𝐷_𝑡(𝑅)−𝜃_𝑡𝐷_𝑡Θ(𝑇),_[𝑟]=𝐷_𝑟(Θ)−𝜃_𝑟𝐷_𝑟(𝑅)−𝜃_𝑡𝐷_𝑟Θ(𝑇),_[𝑟𝑟]=𝐷_𝑟(Θ_[𝑟])−𝜃_𝑟𝑟𝐷_𝑟(𝑅)−𝜃_𝑟𝑡𝐷_𝑟(𝑇), (3.4)and 𝐷_𝑟 and 𝐷_𝑡 are the operators of total differentiation with respect to 𝑟 and 𝑡, respectively. The operator Γ is a point symmetry of the governing equation (2.5), ifΓ⁽²⁾(Eqn(2.5))|_Eqn(2.5)=0. (3.5)Since the coefficients of Γ do not involve derivatives, we can separate (3.5) with respect to the derivatives of 𝜃 and solve the resulting overdetermined system of linear homogeneous partial differential equations known as the determining equations. Further calculations are omitted at this stage as they were facilitated by a freely available package Dimsym [14], a subprogram of Reduce [15].

The admitted Lie algebra is three-dimensional and spanned by the base vectorsΓ₁=𝜕𝜕𝑡,Γ₂=𝑒^𝑛𝑠𝑡𝜕𝑛𝑠(𝑠𝜃−𝜕𝜕𝜃𝜕𝑡),Γ₃=2𝜃𝑛𝜕𝜕𝜕𝜃+𝑟𝜕𝑟. (3.6)

Extra symmetries may be obtained for the cases (a)𝑚=(3𝑛+4)/(𝑛+2);(b)𝑚=0,𝑛=0;(c)𝑛=0,𝑚=2;(d)𝑚=1,𝑛=−1; and (e)𝑚=0,𝑛=−4/3. One may note that except for case (a), all these cases are realistic. 𝑛=0 renders the governing equation (2.5) linear and we herein omit this case. Extra symmetries, for which 𝑛≠0, are listed in Table 1. Physically the parameters 𝑚 and 𝑛 are not related (since 𝑚 represents the geometry and 𝑛 is exponent of the thermal conductivity). However, it is interesting from symmetry analysis point of view to note that the given relationship between 𝑚 and 𝑛 leads to extra symmetries being admitted.

Table 1 

Extra-admitted symmetries.

4. One-dimensional Optimal Systems of Subalgebras

In this section, we determine nonequivalent subalgebras of the symmetry algebra admitted by (2.5) that is, we construct the one-dimensional optimal system of the symmetry algebra given in (3.6). Reduction of independent variables by one is possible using any linear combination of the admitted base vectors. In order to ensure that a minimal complete set of reductions is obtained from the admitted Lie algebra, an optimal system is constructed (see, e.g., [11, 12]). An optimal system of a Lie algebra is a set of 𝑙 dimensional subalgebras such that every 𝑙 dimensional subalgebra is equivalent to a unique element of the set under some element of the adjoint representation [11]:Ad(exp(𝜖Γ_𝑖))Γ_𝑗=^∞_𝑛=0𝜖^𝑛𝑛!(AdΓ_𝑖)^𝑛Γ_𝑗=Γ_𝑗−𝜖[Γ_𝑖,Γ_𝑗𝜖]+²2[Γ_𝑖,[Γ_𝑖,Γ_𝑗]]−⋯, (4.1)where [Γ_𝑖,Γ_𝑗]=Γ_𝑖Γ_𝑗−Γ_𝑗Γ_𝑖 is the commutator of Γ_𝑖 and Γ_𝑗. To compute the one-dimensional optimal system of the algebra in (3.6), first a commutator table is constructed and given by Table 2.

Table 2 

Commutator table.

The adjoint representation is constructed using formula (4.1). We wish to simplify as much as possible the coefficients 𝑎₁,𝑎₂, and 𝑎₃ by carefully applying the adjoint maps toΓ=𝑎₁Γ₁+𝑎₂Γ₂+𝑎₃Γ₃. (4.2)Starting with a nonzero vector (4.2) with 𝑎₁≠0 and rescaling Γ such that 𝑎₁=1, it follows from Table 3 that acting on Γ by Ad(exp((−𝑎₂/𝑛𝑠)Γ₂)), one obtains Γ₁+𝛼Γ₃. No further simplification is possible. For 𝑎=0 and assuming 𝑎₃≠0 (say, 𝑎₃=1 by rescaling); acting on the remaining vector (4.2) by Ad(exp(𝑐₁Γ₂)) yields 𝑎₂𝑒^{𝑐1𝑛𝑠}Γ₂+Γ₃. However, depending on the sign 𝑎₂, the coefficient of Γ₂ can be assigned to either +1, −1, or 0. Finally, for 𝑎₃=0, we obtain Γ₂. Thus the one-dimensional optimal system is given by {Γ₁+𝛼Γ₃;Γ₃±Γ₂;Γ₃;Γ₂}, (4.3)where 𝛼 is an arbitrary constant.

Table 3 

Adjoint representation for the base vectors given in (3.6).

5. Symmetry Reductions and Invariant Solutions

In this subsection, we reduce the variables of the governing BVP by one. We provide the invariant solution constructed using Γ₂ in (3.6) which satisfies the prescribed boundary condition. Reductions by members of the one-dimensional optimal system are listed in Table 4. In order to find invariants of Γ₂, we have to solve the system𝑑𝜃=𝑠𝜃𝑑𝑡=−1𝑑𝑟0. (5.1)The system (5.1) yields the invariants 𝐽₁=ln𝜃+𝑠𝑡and𝐽₂=𝑟 which give rise to the functional form𝜃=𝑒^−𝑠𝑡𝐺(𝑟). (5.2)

Table 4 

Reductions by elements of the optimal systems.

The time-dependent heat transfer coefficient may be represented by 𝑓(𝑡)=𝑒^{−𝑛𝑠𝑡}. This choice of 𝑓(𝑡) renders the boundary condition invariant under Γ₂ and it is also realistic (note that the form of 𝑓(𝑡) is obtained by substituting (5.2) into (2.6)) (see also [6, 16]). Substituting this expression for 𝑓(𝑡) and the functional form (5.2) into the governing equation (2.5), one obtains𝑛𝐺⁻¹(𝐺^′)²+𝑚𝑟𝐺^′+𝐺^′′=0, (5.3)and the boundary conditions (2.6) transform to

𝑑𝐺𝑑𝑟=0,𝑟=0,𝑑𝐺𝑑𝑟=−Bi𝐺^1−𝑛,𝑟=1. (5.4)

Note that the trivial solution to (5.3) is given by a constant. Four cases arise for the nontrivial solution of (5.3) subject to different choices of m and n (see also [17, page 365]). Case (a)𝑛=−1,𝑚≠1,𝐺=𝑐₂exp{𝑐₁𝑟^1−𝑚}. (5.5)Case (b)𝑛=−1,𝑚=1,𝐺=𝑐₂𝑟^𝑐1. (5.6)Case (c)𝑛≠−1,𝑚=1,𝐺=±[𝑐₁ln𝑟+𝑐₂]^−(𝑛+1). (5.7)Case (d)𝑛≠−1,𝑚≠1,𝑐𝐺=±[¹(𝑛+1)𝑟^1−𝑚+(𝑚−1)(𝑛+1)𝑐₂]𝑚−1^1/(𝑛+1). (5.8)

We consider cases (a) and (d) only as examples. For case (a) and in terms of original variables we obtain𝑐𝜃=¹(𝑚−1)Bi𝑒^𝑐1exp{−𝑠𝑡+𝑐₁𝑟^1−𝑚}. (5.9)

The invariant solution (5.9) satisfies the prescribed boundary condition (2.6) at 𝑟=0 only if 𝑚=2 (spherical geometry) and for any constant 𝑐₁<0. One may rewrite (5.9) as𝑐𝜃=¹Bi𝑒^𝑐1exp{−𝑠𝑡+𝑐₁𝑟⁻¹}. (5.10)

Without loss of generality we let 𝑐₂=0 in case (d) and in terms of original variables, the invariant solution satisfying the boundary conditions (2.6) for 𝑚=0 (rectangular geometry) is given by𝜃=(−Bi)^1/𝑛(𝑛+1)^1/𝑛𝑒^−𝑠𝑡𝑟^1/(𝑛+1), (5.11)where values of 𝑛 must be chosen such that the singularity at 𝑟=0 is avoided. Furthermore, we obtain real solutions for −1<𝑛<0 and 0<𝑛<1. Invariant solution (5.11) is depicted in Figures 1 and 2.

Figure 1 

Graphical representation of the invariant solution (5.11). Parameters used 𝑛=0.5,𝑠=1, and 𝑟=0.55. The temperature profile is given for varying Biot number.

Figure 2 

Graphical representation of the invariant solution (5.11). Parameters used 𝑛=−0.33,𝑠=1, and 𝑡=1. Temperature is plotted against spatial variable.

Symmetry analysis may lead to extra solutions, if we use the linear combinations of the admitted symmetries or elements of the optimal systems as listed in Table 4. For example, the Γ₃-invariant is given by𝜃=𝑟^(2/𝑛)[𝑚𝑛+2+𝑛+𝑐₁𝑛𝑠𝑒^𝑛𝑠𝑡]𝑛𝑠^−1/𝑛. (5.12)

6. Some Discussions and Conclusions

We have determined some examples of group invariant solutions which satisfy the realistic boundary conditions (it is a well-known fact that more often symmetries do not lead to solutions which satisfy the boundary conditions). In this manuscript, Lie group analysis resulted in some exotic admitted point symmetries. Furthermore, reduction by one variable of the governing equation has been performed using members of the optimal system.

The invariant solution (5.11) shows thermal decay due to heat losses by convection to the surrounding environment. Figure 1, depicts an increase in the Biot number due to the resistance of the medium surface heat losses which leads to an increase in the medium temperature, and hence enhancing its energy storage capabilities. The medium temperature decreases with time. In Figure 2, temperature is much lower at the device surface than at 𝑟=0, and this is due to heat loss to the surrounding.

We have given some exact (invariant) solutions to nonlinear heat diffusion equations with temperature-dependent conductivity and time-dependent heat transfer coefficient.

Acknowledgments

Raseelo J. Moitsheki wishes to thank the National Research Foundation of South Africa under Thuthuka program, for the generous financial support. The author is also grateful to the anonymous reviewers for their useful comments.