Homogenization of time dependent boundary conditions for multi-layer heat conduction problem in cylindrical polar coordinates

International Journal of Heat and Mass Transfer 129 (2019) 721–734

Contents lists available at ScienceDirect

International Journal of Heat and Mass Transfer journal homepage: www.elsevier.com/locate/ijhmt

Homogenization of time dependent boundary conditions for multi-layer heat conduction problem in cylindrical polar coordinates Pranay Biswas a, Suneet Singh a, Hitesh Bindra b,⇑ a b

Department of Energy Science & Engineering, Indian Institute of Technology Bombay, Mumbai 400076, India Department of Mechanical and Nuclear Engineering, Kansas State University, 3002 Rathbone Hall, 1701B Platt St., Manhattan, KS 66506, USA

a r t i c l e

i n f o

Article history: Received 14 May 2018 Received in revised form 6 October 2018 Accepted 7 October 2018

Keywords: Analytical solution Multi-layer heat conduction Time dependent boundary conditions Homogenization Auxiliary function Finite integral transformation

a b s t r a c t In eigenfunction based solutions for heat conduction problems, the inhomogeneous boundary conditions (BCs) and the eigenfunctions (corresponding to homogeneous BCs) are incompatible with each other at the boundaries. Therefore, homogenization of BCs is essential to avoid the abovementioned incompatibility. However, homogenization is usually quite tedious especially for time dependent BCs. The process becomes even more difficult for multi-dimensional multi-layer problems. These difficulties come from the fact that, in general, an auxiliary function (quasi-static part of temperature) needs to be subtracted from the temperature of each layer of the composite. This treatment yields a set of conditions on auxiliary functions. In this work it is demonstrated that these auxiliary functions are not unique and subsequently, two formulations of auxiliary function are proposed here for cylindrical polar (r, h) coordinates. The second formulation leads to relatively simple computational procedure. The computational ease of the proposed methodology arises from the fact that auxiliary functions are defined only in the innermost and the outermost layers of the composite in this approach. The solutions of a few illustrative examples are obtained without homogenization as well as with homogenization (using two distinct approaches). The results without homogenization of BCs exhibits large mismatch near the boundaries, while such mismatch goes way in case of homogenized BCs. Moreover, the novel (second) approach proposed here for homogenization is easy to implement and efficient compared to the first formulation approach. Ó 2018 Published by Elsevier Ltd.

1. Introduction The multi-layer composites such as slabs, cylinders and spheres have attracted considerable attention in a varied range of applications due to added advantage of combining physical, mechanical and thermal properties of different materials. These composites are widely used in automotive, space, chemical and nuclear industries. Real life systems need accurate and efficient determination of heat flux and temperature distributions within each layer of the composite. In order to find temperature and flux distributions by solving multi-layer heat conduction problems, numerical as well as analytical approaches are used [1–3]. The analytical solutions provide better insight into the underlying physics through the mathematical forms. Moreover, these solutions are often used to verify newly developed numerical models. With the advent of softwares capable of symbolic manipulations [4–8], there has been renewed interest in analytical solutions of such problems. ⇑ Corresponding author. E-mail addresses: [email protected] (P. Biswas), [email protected] (S. Singh), [email protected] (H. Bindra). https://doi.org/10.1016/j.ijheatmasstransfer.2018.10.036 0017-9310/Ó 2018 Published by Elsevier Ltd.

The analytical solutions of 1-D multi-layer transient heat conduction problems in Cartesian, cylindrical and spherical coordinates have been developed over several decades. The solutions were developed by using separation of variables (SOV) [1–3], orthogonal expansion [9–11], eigenfunction expansion [12,13], finite integral transform (FIT) [14–18], Green’s functions (GFs) [19–23], Laplace transform (LT) [24,25], Duhamel’s theorem approach [26] and the Galerkin procedure [27]. These 1-D solutions were further extended to multi-dimensional heat conduction in Cartesian coordinate systems [28–33]. However, such multi-dimensional solutions in cylindrical [34–36] and spherical [37] coordinates are neither simple nor straightforward. This is especially true for laminated composites [38–40], composites containing functionally graded materials [41,42] and composites with reinforced as well as fiber materials [37,43,44]. In various scientific, engineering and industrial applications the composites with these materials are widely used and it is vital to determine temperature as well as flux distributions within the composites. In several eigenfunction based approaches, the analytical solutions of the partial differential equations (PDEs) are represented as a linear combination of eigenfunctions. Thus, it is essential to

722

P. Biswas et al. / International Journal of Heat and Mass Transfer 129 (2019) 721–734

Nomenclature a and b

coefficients of Bessel functions used to describe radial eigenfunctions

A; B and C coefficients used to describe boundary conditions (BCs) f

function used to describe time dependency of BCs

F F

initial temperature distributions finite integral transform (FIT) of initial temperature distributions volumetric heat source FIT of volumetric heat source FIT of the first term of modified heat source Bessel’s function of first kind thermal conductivity number of layers in the composite norm in radial direction truncation term in the series for radial direction auxiliary function heat flux magnitude of the applied heat flux at the outermost boundary for example 2 radial coordinate radial eigenfunction inhomogeneous part of ODEs time temperature magnitude of the applied temperature at the outermost boundary for example 1 FIT of temperature function

g G I J k M N P q q00 q000 r R S t T T0 T u and V W X Y Z

v

Greek letters thermal diffusivity C expansion coefficient of temperature h angular coordinate # norm in angular direction k radial directional eigenvalue l derived boundary functions related to the FIT approach without homogenized BCs / and w spatial functions related to auxiliary functions in the first formulation X any quantity of interest x FIT of quantity of interest

a

Abbreviations 1-D, 2-D, 3-D one dimensional, two dimensional, three dimensional BCs boundary conditions FIT finite integral transform GFs Green’s function ICs initial conditions LT Laplace transform ODEs ordinary differential equations PDEs partial differential equations SOV separation of variables

spatial functions to describe the solution of the Euler’s equation in the first formulation

Superscripts c Fourier coefficient associated with cosine function h modified hD reduced modified D reduced s Fourier coefficient associated with sine function

FIT of the first term of modified heat source truncation term in the series for angular direction spatial function to define auxiliary function for the innermost layer in the second formulation Bessel’s function of second kind spatial function to define auxiliary function for the outermost layer in the second formulation

Subscripts i layer or interface in inner surface of the 1st layer m mth term in series for angular eigenfunctions n nth term in series for radial eigenfunctions out outer surface of the Mth layer

homogenize the boundary conditions (BCs) for good convergence of the series solution. In case of time independent BCs, the homogenization process is relatively straightforward and the SOV approach is preferred to find solution for 1-D composite slabs [4,5]. However, numerical calculations for solutions, in case of multi-dimensional composite slabs, become complicated due to the presence of imaginary eigenvalues [6–8]. The analytical solutions for 2-D polar ðr; hÞ [45–48] and 3-D spherical [49] geometries are free from imaginary eigenvalues and therefore facilitate ease of computation. Moreover, the SOV is not applicable to heat conduction problem with time dependent BCs as well as sources and the homogenization process for such BCs becomes difficult. The LT approach is often used to obtain analytical solution of transient heat conduction problem by transforming solution from frequency domain to time domain by using Fourier–Mellin or Bromwich integral. However, the evaluation of this integral is neither straightforward nor trivial. In practice, various numerical methods are developed for the inversion process. The process, in general, deals with evaluation of contour integrals [50] and also requires a very large number of terms in the series for good convergence. The LT approach becomes very inefficient for multi-layer problems because the inversion needs to be carried out for the solution in

each layer of the composite. The approach is even more difficult to apply to multi-dimensional multi-layer heat conductions. The GFs approach and Duhamel’s theorem [26] are also used to solve heat conduction problems with time dependent BCs, without homogenization of BCs. However, implementation of these approaches for the multi-dimensional multi-layer heat conduction problems (especially in spherical and polar geometries with real eigenvalues) is very tedious. The FIT is another widely used eigenfunction based approach. A general formulation of the FIT for multi-region heat conduction problem is presented in [51]. Such analytical solutions for multidimensional transient as well as steady state heat conduction problems are available in literature [17,52–54]. These solutions mainly used time dependent BCs and in general does not require homogenization of BCs. However, this approach produces erroneous solutions at/near boundaries due to mismatch between the BCs and corresponding eigenfunctions [15,16]. In order to remove this type of mismatch a generalized process for homogenization of BCs is essential. For 1-D multi-layer heat conduction, the homogenization of time dependent BCs could be carried out by using approach available in published literature [26]. In this formulation, the

723

P. Biswas et al. / International Journal of Heat and Mass Transfer 129 (2019) 721–734

quasi-static part of temperature (auxiliary function) for each layer of the composite is obtained by setting the Laplacian of auxiliary function to zero (constrained auxiliary function, see Chapter 8 in Ref. [26]). The formulation works for both single as well as multi-layer heat conduction problems. However, its extension to multi-dimensional multi-layer cases is neither straightforward nor trivial. Some work pertinent to constrained auxiliary functions for 2-D cylindrical polar ðr; hÞ coordinates is available in the work by Wang et al. [55]. For 1-D multi-layer problems, an alternate form for auxiliary functions was presented by Antonopoulos et al. [10]. In the alternate approach the auxiliary functions are needed to be defined only for the innermost and the outermost layers of the composite. Such definition of auxiliary function is not applicable to the single layer geometry. It should also be mentioned that the form of function presented in Ref. [10] fails to account for the non-zero heat flux condition (Neumann BCs). The extension of such alternate form of auxiliary functions for multidimensional heat conduction is not well established yet. In order to bridge in this research gap, we have formulated a novel strategy to solve for the single as well as multilayer multidimensional transient heat conduction problems while considering time dependent BCs and sources to provide a reliable framework for the verification and validation of various experimental as well as numerical results. In this paper, analytical solution of multi-layer heat conduction problem in cylindrical polar coordinates ðr; hÞ is found for time dependent BCs and sources. A generalized methodology is proposed to homogenize BCs in the case of single as well as multilayer composites with any kind of BCs. In the current method, an auxiliary function is defined and is subtracted from the original temperature, and an equation is obtained in terms of the modified temperature. The auxiliary function needs to be defined in such a way that the modified temperature equation has homogeneous BCs. This leads to a set of conditions on the auxiliary functions that are not uniquely defined. Therefore, alternative approaches are possible for the homogenization based on different auxiliary functions. Two different forms for auxiliary functions are proposed here and they are termed as ‘‘first formulation” and ‘‘second formulation”. The first formulation deals with evaluation of auxiliary function in each layer of the composite and analogous to the procedure presented in Ref. [55]. On the other hand, the auxiliary functions are defined second formulation in based on form presented in Ref. [10] (although only for 1-D problems), and are applicable to the single as well as multi-layer heat conduction problems with any kind of BCs. In this study it has been shown that both first and second formulation are independently satisfying the same set of conditions on auxiliary functions. However, the second formulation provides added advantage in terms of ease in implementation and facilitates comparatively lower computational cost compare to the first formulation. It is obvious that once homogenization of BCs is done any eigenfunction based approach could be used and the same solution is achieved. In illustrative examples, the multi-layer heat conduction problems are solved by using FIT approach for time dependent Dirichlet as well as Neumann BCs in two different methods, i.e. without and with homogenized BCs. The solutions with homogenized BCs are compared to those without homogenized BCs.

2. Problem description An M-layer annulus composite medium, as illustrated in Fig. 1 is considered here. Initially (t ¼ 0), each layer is at specific temperature T i ðr; h; t ¼ 0Þ ¼ F i ðr; hÞ, in ri1 < r < ri , i ¼ 1; . . . ; M. For times t > 0, energy is generated in each layer at a rate of g i ðr; h; tÞ. The mathematical description of the problem is given here for this

Fig. 1. Schematic representation of an M-layer annulus composite.

2-D multi-layer annulus assuming the interfaces are at perfect thermal contact. Governing equations:


 1 @T i ðr; h; t Þ 1 @ @T i ðr; h; tÞ 1 @ 2 T i ðr; h; t Þ ¼ r þ 2 @t r @r @r r ai @h2 1 þ g i ðr; h; t Þ ki

ð1Þ

BCs and interface conditions: in r-direction  At inner surface of the 1st layer, i ¼ 1

Ain

 @T 1 ðr; h; tÞ  þ Bin T 1 ðr; h; t Þjr0 ¼ C in ðh; tÞ @r

ð2Þ

r0

 At interface between ith and ði þ 1Þth layer, i ¼ 1; . . . ; M  1

T i ðr; h; t Þjri  T iþ1 ðr; h; tÞjri ¼ 0 ki

  @T i ðr; h; t Þ @T iþ1 ðr; h; t Þ  k iþ1   ¼0 @r @r ri ri

ð3Þ ð4Þ

 At outer surface of Mth layer, i ¼ M 

Aout

@T M ðr; h; t Þ  þ Bout T M ðr; h; t ÞjrM ¼ C out ðh; t Þ @r rM

ð5Þ

The variations in temperature and heat flux in the h-direction are periodic and the corresponding BCs are

T i ðr; h; t Þjh¼0 ¼ T i ðr; h; tÞjh¼2p and  @T i ðr; h; tÞ ¼ :  @h h¼2p ICs:

T i ðr; h; t Þjt¼0 ¼ F i ðr; hÞ

 @T i ðr; h; t Þ  @h

h¼0

ð6Þ ð7Þ

The principle objective is to find analytical solution of temperature distribution T i ðr; h; tÞ for the above described heat conduction problem within the composite medium. 3. Solution methodology As the h-directional BCs [Eq. (6)] in the current configuration are periodic and homogeneous in nature, therefore, independent angu-

724

P. Biswas et al. / International Journal of Heat and Mass Transfer 129 (2019) 721–734

lar eigenfunctions turn out to be sinðmhÞ and cosðmhÞ, m 2 I. The PDEs, given in Eqs. (1)–(7) and used to described 2-D heat conduction problem, can be converted to m number of 1-D PDEs with the help of Fourier series expansion using angular eigenfunctions.

In order to solve the reduced 1-D heat conduction equations [in Eqs. (10)–(13)], BCs need to be homogenized to remove mismatch between BCs and corresponding eigenfunctions at/near boundaries.

3.1. Eigenfunction expansion in the h-direction

3.2. Homogenization of BCs

The brief discussion on the Fourier series expansion in the hdirection is given here and details of the expansion technique can be found in Refs. [54,55]. Using the same expansion technique the temperature of each layer of the composite is expressed as given,

nized by defining a reduced modified temperature T hD im ðr; t Þ (dynamic part) for each layer of the composite as follows,

1 1 X X T i ðr; h; t Þ ¼ T cim ðr; t Þ cosðmhÞ þ T sim ðr; t Þ sinðmhÞ: m¼0

ð8Þ

m¼1

Similarly, the heat sources is expanded as follows

g i ðr; h; t Þ ¼

1 X

g cim ðr; tÞ cosðmhÞ þ

m¼0

1 X

g sim ðr; t Þ sinðmhÞ:

ð9Þ

m¼1

The coefficients g cim ðr; tÞ and g sim ðr; tÞ in Eq. (9) are evaluated by applying the orthogonality conditions in the h-direction and are shown below,

Z 2p 1 g cim ðr; t Þ ¼ g i ðr; h; tÞ cosðmhÞdh and #m 0 Z 2p 1 g sim ðr; t Þ ¼ g i ðr; h; tÞ sinðmhÞdh #m 0 

where #m ¼

The BCs of reduced 1-D heat conduction problem are homoge-

D T Dim ðr; tÞ ¼ T hD im ðr; t Þ þ qim ðr; t Þ

where qDim ðr; t Þ is auxiliary function (quasi-static part). The function qDim ðr; t Þ is the representation of two distinct functions qcim ðr; t Þ and qsim ðr; t Þ. The substitution of T Dim ðr; tÞ obtained from Eq. (14) into reduced 1-D heat conduction problem [in Eqs. (10)–(13)] gives rise to reduced modified 1-D heat conduction problem as stated below: Reduced modified 1-D governing equations:

! 1 @T hD 1 @ @T hD m2 1 im ðr; t Þ im ðr; t Þ  2 T hD ¼ r ðr; tÞ þ g hD ðr; t Þ @t r @r @r ki im ai r im 

direction. The same expansion procedure is carried out for boundary functions C in ðh; t Þ, C out ðh; t Þ and ICs F i ðr; hÞ to obtain coefficients C cin;m ðt Þ, C sin;m ðt Þ, C cout;m ðt Þ, C sout;m ðtÞ, F cim ðrÞ and F sim ðr Þ, respectively. The Fourier series expansion of several quantities discussed above are used in Eqs. (1)–(7) yielding two sets of reduced 1-D heat conduction problem equations in terms of reduced temperature T Dim ðr; tÞ.

The reduced 1-D governing equations and corresponding reduced ICs as well as BCs are essentially similar irrespective of the superscript c and s. Therefore, D is used in the superscript to distinguish between them in the general formulation of the reduced 1-D heat conduction problem described in Eqs. (10)–(13).

g hD im ðr; t Þ ¼ ki

ð10Þ

ð11Þ

r0

Aout

 @T DMm ðr; t Þ   @r

  þ Bout T DMm ðr; t Þ

rM

rM

¼ C Dout;m ðt Þ:

ð12Þ

In addition to above conditions, the interface conditions for modified temperature are homogeneous and similar to those given by Eqs. (3) and (4). Such conditions can be used by directly replacing T i ðr; h; t Þ with T Dim ðr; t Þ.

Reduced modified homogeneous BCs:  @T hD ðr; t Þ  hD Ain 1m  þ Bin T 1m ðr; t Þ ¼ 0  @r r0

t¼0

¼ F Dim ðrÞ:

ð17Þ

r0

  @T hD Mm ðr; t Þ Aout   @r

  þ Bout T hD Mm ðr; t Þ

rM

rM

¼0

ð18Þ

The interface conditions can be obtained by simply replacing T i ðr; h; t Þ with T hD im ðr; t Þ in Eqs. (3) and (4).

  T hD im ðr; t Þ

t¼0

¼ F Dim ðr Þ  qDim ðr; 0Þ:

ð13Þ

ð19Þ

The reduced modified 1-D problem for modified temperature described by Eqs. (15)–(19) is solved by any eigenfunction based approach that requires homogeneous BCs, by setting an eigenvalue problem for each layer of the composite. However, such homogenization is possible if and only if the function qDim ðr; t Þ satisfies following set of conditions. The set of conditions on auxiliary function: at two extreme surfaces,

Ain

  @qD1m ðr; tÞ D D   þ Bin q1m ðr; tÞ r0 ¼ C in;m ðt Þ @r r0

 @qDM;m ðr; t Þ Aout   @r

rM

Reduced ICs:

  T Dim ðr; t Þ

ð16Þ

Reduced modified ICs:

Reduced BCs: in r-direction:

 @T D1m ðr; tÞ  D D  þ Bin T 1m ðr; t Þ ¼ C in;m ðtÞ  r0 @r


  @qD ðr; t Þ 1 @ m2 1 @qDim ðr; t Þ r im  2 qDim ðr; tÞ  r @r r ai @r @t

þ g Dim ðr; t Þ

Reduced 1-D governing equations:

Ain

ð15Þ

where modified source term g hD im ðr; t Þ can be expressed as follows,

2p if m ¼ 0 p if m – 0 is normalization factor (norm) in the h-

! 1 @T Dim ðr; t Þ 1 @ @T D ðr; t Þ m2 1 ¼ r im  2 T Dim ðr; t Þ þ g Dim ðr; tÞ @t r @r @r ki ai r

ð14Þ

  þ Bout qDM;m ðr; t Þ

rM

¼ C Dout;m ðt Þ:

ð20Þ

ð21Þ

Additionally, interface conditions are analogous to those described in Eqs. (3) and (4) and are used by replacing T i ðr; h; tÞ with qDim ðr; tÞ. The whole set of conditions is obtained by using

725

P. Biswas et al. / International Journal of Heat and Mass Transfer 129 (2019) 721–734

the fact that the sum of reduced modified temperature and auxiliary function should satisfy each and every equation of reduced 1-D problem [in Eqs. (10)–(13)]. It should be mentioned that although qDim ðr; t Þ is quintessential in obtaining the final solution of the multi-layer heat conduction problem but the evaluation of qDim ðr; tÞ is not straightforward. Both first as well as second formulations are applied to the multi-layer heat condition problem and a detailed description of the methodology is given in the following section. 3.3. Definition and evaluation of auxiliary function A. First formulation: The procedure to evaluate auxiliary function for 2-D multilayer heat conduction was provided by Wang et al. [55] for fiber-reinforced cylindrical composites. The similar method is being used here for cylindrical composites with any kind of time dependent BCs. Each function qDim ðr; t Þ for this formulation satisfies an additional ‘‘pseudo steady state” condition which is defined as follows:


 @qD ðr; t Þ 1 @ m2 r im  2 qDim ðr; t Þ ¼ 0; r @r r @r

r i1 6 r 6 r i :

ð22Þ

This particular choice is made to simplify the modified source term given in Eq. (16). Note that Eq. (22) is solvable if qDim ðr; t Þ is either a separable function, i.e. product of space and time components or linear combinations of separable functions. In this formulation, the auxiliary function qDim ðr; t Þ for each layer of the composite is defined as follows [26],

qDim ðr; tÞ ¼ wim ðrÞC Din;m ðt Þ þ /im ðr ÞC Dout;m ðtÞ;

r i1 6 r 6 r i

ð23Þ

where spatial functions wim ðr Þ and /im ðr Þ are solutions of two independent steady state problems. These steady state problems are constructed by combining the additional condition [in Eq. (22)] and the form in Eq. (23) for each auxiliary function along with the set of conditions on them in Eqs. (20) and (21). The procedure to evaluate spatial functions wim ðrÞ and /im ðr Þ is described here subsequently and is demonstrated only for wi;m ðr Þ. Governing equations for wim ðr Þ: The following ordinary differential equations (ODEs) are constructed by substituting Eq. (23) into Eq. (22),


 1 d dw ðrÞ m2  2 wim ðr Þ ¼ 0; r im r dr dr r

r i1 6 r 6 ri

ð24Þ

½Em 2M2M ½Cm 2M1 ¼ ½B1 2M1 ) ½Cm 2M1 h i ¼ E1 ½B1 2M1 m

ð28Þ

2M2M

C2M;m T and B1 ¼ ½1 h i with 2M  1 elements. The matrix E1 m where Cm ¼ ½C1;m

...

C2;m

2M2M

0

...

0T are

is inverse of

½Em 2M2M which is obtained by putting Eq. (23) into the set of conditions in Eqs. (20) and (21). The matrix ½Em 2M2M has 4M2 elements, of which 4ð2M  1Þ elements are non-zero. In this formulation, the spatial functions wim ðrÞ are unique for each layer of the composite. A similar solution approach is followed to evaluate unique solutions for /im ðr Þ by replacing B1 with B2 , where B2 ¼ ½0 0 . . . 1T . The final solution for qDim ðr; tÞ is obtained by using evaluated spatial functions wim ðr Þ and /im ðr Þ in Eq. (23). The evaluated function qDim ðr; t Þ, for each layer of the composite, is also unique in this formulation and is one among several possible solutions. For an M-layer composite, one needs to evaluate 2M  m number spatial functions [wim ðrÞ and /im ðr Þ, M  m evaluation for each of them]. This is due to fact that functions um ðrÞ and v m ðrÞ in Eq. (27), are dependent of m. Thus, the evolution of qDim ðr; tÞ in the first formulation is computationally expensive. Furthermore, the implementation and code-development for the quantities of interest such as modified heat sources as well as modified ICs are also very time consuming and complex. B. Second formulation: The auxiliary functions prescribed here are based on those presented by Antonopoulos et al. [10] for 1D multi-layer heat conduction problems. These functions are modified here to incorporate Neumann BCs as well. In second formulation, the auxiliary functions are defined as follows,

qD1m ðr; t Þ ¼ X ðr ÞC Din;m ðt Þ; qDim ðr; t Þ ¼ 0;

r0 6 r 6 r1

ð29Þ

i ¼ 2; . . . ; M  1

qDMm ðr; t Þ ¼ Z ðr ÞC Dout;m ðt Þ;

ð30Þ

r M1 6 r 6 rM

ð31Þ

where

X ðr Þ ¼

ðr 1  r Þ2 ðr 1  r0 Þ½2Ain þ Bin ðr1  r0 Þ

Z ðr Þ ¼

ðr  r M1 Þ2 : ðrM  r M1 Þ½2Aout þ Bout ðr M  r M1 Þ

and

Additionally, the intermediate conditions are also similar to homogenous interface conditions mentioned in Eqs. (3) and (4) and are used during evaluation of wim ðr Þ by replacing T i ðr; h; t Þ with wim ðrÞ. The solutions for wim ðrÞ follow the following form:

Note that the form of functions as given in Eqs. (29)–(31) are defined only for the innermost as well as the outermost layers of the composite and also satisfy the same set of conditions applied on qDim ðr; tÞ in Eqs. (20) and (21), including the interface conditions. Moreover, functions X ðr Þ and Z ðr Þ are represented in a way such that X ðr0 Þ ¼ 1, X ðr1 Þ ¼ 0, Z ðr M1 Þ ¼ 0 and Z ðr M Þ ¼ 1. The form of function proposed here in Eqs. (29)–(31) works with any kind of BCs, however, the form of function presented in Ref. [10] does not work for Neumann BCs case. Moreover, the auxiliary functions prescribed by Antonopoulos et al. [10] in such a way that cannot be applied to a single layer heat conduction problem. The following modification is required while defining the auxiliary functions qDm ðr; tÞ for such heat conduction problem in multi-dimensions and is as follows (i is omitted the subscript),

wim ðr Þ ¼ Cð2i1Þm um ðr Þ þ Cð2iÞm v m ðr Þ

qDm ðr; t Þ ¼ X ðrÞC Din;m ðt Þ þ Z ðr ÞC Dout;m ðtÞ;

with the following imposed constraints

Ain

 dw1m ðrÞ þ Bin w1m ðr Þjr0 ¼ 1 dr r0

Aout

 dwMm ðr Þ þ Bout wMm ðrÞjrM ¼ 0: dr rM



ð25Þ

ð26Þ

ð27Þ

ln r if m ¼ 0 where um ðr Þ ¼ and v m ðr Þ ¼ r m . The coefficients rm if m – 0 Cim ; i ¼ 1; . . . ; 2M (as each layer has two coefficients), in Eq. (27) are found by matrix inversion method shown below,

r0 6 r 6 r1

ð32Þ

2

0Þ where Z ðr Þ ¼ ðr1 r0 Þ½2Aðrr . With this modification, the second out þBout ðr 1 r 0 Þ

formulation could be applied to single as well as multi-layer heat condition problems by providing only two non-zero spatial

726

P. Biswas et al. / International Journal of Heat and Mass Transfer 129 (2019) 721–734

functions [X ðr Þ; r0 6 r 6 r 1 and Z ðr Þ; r M1 6 r 6 r M ] associated with two extreme boundary functions instead of the extreme layers. Just like in the case of first formulation, the spatial functions X ðrÞ and Z ðr Þ used in second formulation are also unique but they are independent of m (eigenvalues in the h-direction). Moreover, only two functions that are independent of m, are sufficient to define the auxiliary functions throughout the entire simulation for this formulation. Therefore, the implementation becomes very simple and the computational cost reduces drastically when mathematical manipulations are performed on the modified sources and modified ICs. Now, the substitution of qDim ðr; tÞ either from first formulation or from second formulation into Eq. (14) does not provide complete solution for the reduced 1-D temperature T Dim ðr; tÞ. Although, the modified sources and modified ICs are known in this case, however, the modified temperature T hD im ðr; t Þ is still unknown. Once the BCs are homogenized, any eigenfunction based approach such as the FIT or orthogonal eigenfunction expansion method could be used to obtain the analytical solutions of T hD im ðr; t Þ. It is noted here that the solution from both these approaches leads to exactly same final solution. The current study considers the FIT approach and the procedural detail particular to this work is discussed. A general procedure to solve heat conduction is demonstrated that applies to both with as well as without homogenized BCs. 3.4. FIT: ODEs for time component

T hD mn ðt Þ

In the FIT approach, the operator

PM

i¼1 ki

R ri r i1

rRimn ðkimn ; r Þdr is

applied on the reduced modified 1-D problem in Eqs. (15)–(19) to formulate an ODE for time component T hD mn ðt Þ corresponding to the temperature field T hD im ðr; t Þ and is as followed:

dT hD mn ðt Þ dt

D þ ai k2imn T hD mn ðt Þ ¼ Smn ðt Þ

ð33Þ

D D D D with IC T Dmn ð0Þ ¼ FhD mn and Smn ðt Þ ¼ Vmn ðt Þ  Imn ðt Þ þ Gmn ðt Þ (see Ref.

[54] for detail of the FIT approach). Note that the quantity Vhmn ðtÞ is zero for first formulation as it obeys the additional condition (pseudo steady state). The quantities in the above equation are found by using the FIT (for any quantity of interest [XDim ðr; t Þ]) and is defined in Eq. (34) as,

x

D mn ðt Þ

¼

Z M X ki i¼1

ai

ri

r i1

rX

D im ðr; t ÞRimn ðkimn ; r Þdr

ð34Þ th

where Rimn ðkimn ; r Þ is radial eigenfunction for i layer of the composite. The remaining quantities are presented in Table 1. It is to be noted that the solution of Eq. (33) is dependent on the number of layer (i) because the product ai k2imn is present in the specified ODE. However, the time component T hD mn ðt Þ is independent of layer. Therefore, to remove such dependency a mathematical relation of the product ai k2imn between different layers is needed. Moreover, the radial eigenfunctions and corresponding eigenvalues ðkimn Þ are to be found by solving the eigenvalue problem discussed in the following subsection. 3.4.1. Finding of radial eigenfunctions and corresponding eigenvalues The details of eigenvalue problem for multi-layer heat conduction have been previously shown by Singh et al. [46,49,54] and Jain et al. [47,48]. A transcendental equation is solved to obtain the radial eigenvalues while, the eigenfunctions for each layer of a cylindrical composite have the following form of solution,

Rimn ðkimn ; r Þ ¼ aimn J m ðkimn r Þ þ bimn Y m ðkimn rÞ

ð35Þ

Table 1 FIT of the required quantities such as temperature, modified temperature, ICs, modified ICs and terms in modified sources. Quantity of interest XD im ðr; t Þ

Transformed quantity xD mn ðt Þ

1. Temperature

TD im ðr; t Þ

TD mn ðt Þ

2. Modified temperature 3. IC

D D T hD im ðr; t Þ ¼ T im ðr; t Þ  qim ðr; t Þ

T hD mn ðt Þ

FD im ðr Þ

FD mn

4. Modified IC

D D F hD im ðr Þ ¼ F im ðr Þ  qim ðr; 0Þ h  D i 2 @qim ðr;t Þ 1 @ ai r @r r @r  mr2 qD im ðr; t Þ

FhD mn

@qD ðr;t Þ im @t

ID mn ðt Þ

ai D g ðr; t Þ k im

GD mn ðt Þ

Name/Specification

5. First term in modified heat source 6. Second term in modified heat source 7. Original heat source

i

VD mn ðt Þ

where aimn and bimn are coefficients. The radial eigenfunctions Rimn ðkimn ; rÞ satisfy the following orthogonality condition [46]:

  0 rRimn ðkimn ; r ÞRimp kimp ; r dr ¼ Nmn r i1

Z M X ki

a i¼1 i

ri

where Nmn ¼

PM

ki i¼1 ai

R ri ri1

if p–n if p ¼ n

ð36Þ

rR2imn ðkimn ; r Þdr is the normalization factor.

This is possible when the following conditions is true:

pffiffiffiffiffiffiffiffiffiffiffiffi

a1 k21mn ¼ ai k2imn ) kimn ¼ k1mn a1 =ai :

ð37Þ

hD 3.4.2. Evaluation of time component T hD mn ðt Þ and inverse FIT for T im ðr; t Þ Using the relation between the eigenvalues given Eq. (37), the Eq. (33) is rewritten as follows:

dT hD D mn ðt Þ þ a1 k21mn T hD mn ðt Þ ¼ Smn ðt Þ dt

ð38Þ

hD with T hD mn ð0Þ ¼ Fmn . This ODE is independent of layer. The solution of Eq. (38) is as follows,

  Z t 2 2 hD T hD SDmn ðsÞea1 k1mn s ds ea1 k1mn t : mn ðt Þ ¼ Fmn þ

ð39Þ

0

The time component T hD mn ðt Þ is a representation of two distinct hs parts, i.e. T hc mn ðt Þ and T mn ðt Þ. Now, using the inverse FIT provided in Ref. [54], the reduced

modified temperature T hD im ðr; t Þ can be written as follows,

T hD im ðr; t Þ ¼

1 X T hD ðtÞRimn ðkimn ; r Þ mn

n¼1

N mn

:

ð40Þ

It is noted that once homogenization of BCs is performed, the effects from inhomogeneous boundaries are automatically included in the modified sources and therefor the solution given by Eq. (40) is free from boundary effects. However, this treatment does not create any mismatch between the reduced 1-D temperature T Dim ðr; tÞ and original inhomogeneous BCs [as given in Eqs. (11) and (12)], as the auxiliary function qDim ðr; tÞ is added to T hD im ðr; t Þ to obtain T Dim ðr; tÞ [as defined in Eq. (14)]. The solution in Eq. (40) is also can be used to obtain reduced 1D temperature T Dim ðr; t Þ for without homogenized BCs because the solution in Eq. (39) obtained for time component T hD mn ðt Þ holds true for such BCs as well (direct application of the FIT as reported by Singh et al. in Ref. [54]). In that case the following replacements are made: D T hD mn ðt Þ ¼ T mn ðt Þ;

D FhD mn ¼ Fmn

and SDmn ðtÞ

¼ kM r M lDout;mn ðt Þ  k1 r 0 lDin;mn ðt Þ þ GDmn ðt Þ:

727

P. Biswas et al. / International Journal of Heat and Mass Transfer 129 (2019) 721–734

The quantities T hmn ðt Þ and Fhmn are defined in Table 1, whereas, the functions lDin;mn ðtÞ and lDout;mn ðt Þ are constructed from inhomogeneous BCs as mentioned in Eqs. (11) and (12), respectively (see Ref. [54] for detail derivations) as follows,

lDin;mn ðtÞ ¼

8 D  C ðtÞ dR 1mn ðk1mn ;r Þ > <  in;m  B dr

if Ain ¼ 0

> : C Din;m ðtÞ R

if Ain ¼ 0

in

Ain

lDout;mn ðtÞ ¼

r0

1mn ðk1mn ; r 0 Þ

8 D  C ðt Þ dRMmn ðkMmn ;rÞ > <  out;m  dr B out

> : C Dout;m ðtÞ R Aout

rM

Mmn ðkMmn ; r M Þ

if Aout ¼ 0

:

if Aout ¼ 0

It should be mentioned that the mathematical manipulation of lDin;mn ðtÞ and lDout;mn ðtÞ is costly, because each of them is associated with m  n number of terms. Evidently, this solution for reduced 1D temperature T him ðr; tÞ is significantly different from the solution provided in Eq. (40).


  h h T 3 ðr ¼ r 3 ; h; t Þ ¼ T 0 1 þ p  f ðtÞ ; 2 2

0 6 h 6 2p:

ð42Þ

The details of time varying function f ðtÞ are given later in this section. Similar polynomial expression in the h-direction has been previously used by Singh et al. [54] for heat flux condition. The inner surface of the innermost layer (r ¼ r 0 ) is isothermal at zero temperature. This leads to Ain ¼ 0; Bin ¼ 1; Aout ¼ 0; Bout ¼ 1; C in ðh; t Þ ¼ 0 and C out ðh; t Þ ¼ T 3 ðr ¼ r3 ; h; tÞ in the respective BCs. The layers are free from any volumetric heat generations, i.e., g i ðr; h; tÞ ¼ 0. In this illustration of Dirichlet BCs, the radial eigenvalues are found by using procedure given in Refs. [46–49,54]. It is pointed out that, in the results that follow, r, t and T i ðr; h; tÞ are in the units of r0 , r 20 =a1 and T 0 , respectively. The analysis is carried out for the following cases: Case 1: f ðtÞ ¼ 1 and Case 2: f ðt Þ ¼ 1 þ sinðt Þ

3.5. Final solution of temperature The eigenfunctions in Eq. (35) and corresponding eigenvalues are used along with the auxiliary functions as well as the solutions of the reduced modified 1-D temperature [in Eq. (40)] to obtain the final solution of the temperature distribution ½T i ðr; h; t Þ. T i ðr; h; t Þ is expressed by using Eq. (8) as below,

" # 1 X 1 X T hc c mn ðt ÞRimn ðkimn ; r Þ T i ðr; h; t Þ ¼ þ qim ðr; tÞ cosðmhÞ N mn m¼0 n¼1 " # 1 X 1 X T hs mn ðt ÞRimn ðkimn ; r Þ þ qsim ðr; tÞ sinðmhÞ: þ Nmn m¼1 n¼1

t > 0, temperature of the outer surface of the outermost layer (r ¼ r 3 ) is assumed to be of the following mathematical form,

ð41Þ

The FIT solution in Eq. (41) is obtained with homogenized BCs and holds true for any auxiliary function. In order to find solution for T i ðr; h; tÞ by using the FIT without homogenized BCs, the conditions are to be taken as follows: c hs s T hc mn ðt Þ ¼ T mn ðt Þ and T mn ðt Þ ¼ T mn ðt Þ (the time components of the reduced temperature); qcim ðr; t Þ ¼ qsim ðr; tÞ ¼ 0 (i.e. no auxiliary function).

4. Illustrative examples and results In the illustrative examples, the radial eigenvalues and eigenfunctions are found by employing the method given by Singh et al. [46,49,54] and Jain et al. [47,48]. For a particular problem, the infinite series solution for the temperature T i ðr; h; t Þ is truncated at n ¼ P and m ¼ W. In this case, the total number of eigenvalues is ðW þ 1ÞP. For the FIT results with homogenize BCs, both first as well as second formulations are used and are compared with the FIT results without homogenized BCs. In this work, the following non-dimensional values are used: r 1 =r0 ¼ 2, r2 =r0 ¼ 4, r3 =r0 ¼ 6, k2 =k1 ¼ 2 and k3 =k1 ¼ 4, a2 =a1 ¼ 4 and a3 =a1 ¼ 9. Although, the current approach is applicable to any combination of BCs (Dirichlet, Neumann and Robin), however, the present study shows results only for Dirichlet and Neumann BCs. 4.1. Specified temperature (Dirichlet BCs) A three-layer ðM ¼ 3Þ 2-D cylindrical composite (r 0 6 r 6 r 3 , 0 6 h 6 2p, see Fig. 2) is considered in this example. The initial temperature is zero everywhere within the composite. For time

The solutions for case 1 are obtained by using the FIT approach, with and without homogenized BCs. The temporal temperature distributions within the composite are shown in Fig. 3 for (W = 5, P = 5) at different angular directions h = 0, p/4, p/2 and p during time t = 0.1, 0.5, 1, 2 and steady state (ss as t ? 1). The solutions without homogenized BCs are indicated by dotted lines. The approach leads to zero temperature on the outer boundary (r = r3), however, the applied temperature is non-zero irrespective of angular direction. It is due to fact that the eigenfunctions are zero at boundaries for Dirichlet case. The non-zero boundary temperature and the zero eigenfunctions at outer boundary lead to large mismatch which in turn results in spurious oscillations in solutions (dotted fluctuating lines). On the contrary, the solutions with homogenized BCs for both first formulation (solid lines) as well as second formulation (dashed lines) have good match with the boundary temperature as see in Fig. 3 and these solutions are free from spurious oscillations. This is due to the fact that the process of homogenization of BCs is carried out prior to application of the integral transformation. The homogenization of BCs removes the mismatch at boundary between the eigenfunction and the temperature via defining modified temperature and auxiliary function. The BCs for modified problem are same as that used for eigenvalue problem and the inhomogeneity is being taken care by auxiliary function through modified source term. However, the Fig. 3 shows that solutions with those two formulations have slight difference when a few number of terms are considered in the series. Fig. 3 also shows that the temporal distribution of temperature at a particular location with as well as without homogenized BCs has significant difference when Dirichlet BCs are considered. Fig. 4 shows temperature distribution at h = p for different combinations of W and P. As the number of terms increases, the solutions without homogenized BCs show severe oscillations at/near the boundary. Moreover, the oscillations persist near the boundary even with increasing number of terms in the series (dotted lines). As one goes away from the boundary these solutions approaches the correct values, especially if high number of terms are used in the series. On the other hand, the solutions for both first as well as second formulations (with homogenized BCs), are found to produce almost identical results even for a small number of terms. The solutions for first formulation with terms (W = 10, P = 10) are considered as the reference solutions for time t  0.5. It is due to fact that the errors in solutions, having terms (W = 5, P = 5), with the reference solutions are less than 1% in the entire domain (for r0 = 1  r  r3 = 6 and 0  h  2p) of the composite. It

728

P. Biswas et al. / International Journal of Heat and Mass Transfer 129 (2019) 721–734

Fig. 2. 2-D heat conduction in a 3-layer cylindrical composite subjected to time dependent temperature which is a function of angle h and time t.

Fig. 3. Temporal temperature distributions at different angular directions (h = 0, p/4, p/2 and p) using the FIT with as well as without (dotted lines) homogenized BCs during different time instants when constant temperature is applied at the outer surface of the outermost layer. In the case of homogenized BCs, results for first and second formulations are shown by solid and dashed lines, respectively.

729

P. Biswas et al. / International Journal of Heat and Mass Transfer 129 (2019) 721–734

Fig. 4. Effect of number of terms in series summation on temperature distributions at h = p with as well as without homogenized BCs during different time instants when constant temperature is applied at the outer surface of the outermost layer.

is also noted the maximum error occurred at angle h = p where applied temperature is maximum. Accuracy of the second formulation is illustrated via estimating of error in the obtained solutions in comparison with the reference solution (obtained by using first formulation) at location (r = r2, h = p) on the second interface when time t  0.5. It is found that the accuracy is strongly dependent on the number of terms in the r-direction, i.e. on P. The convergence in results is given in Table 2 for the second formulation with different combinations of terms {W, P} for t = 0.5, 1, 2 and also for the steady state (ss). The error in the solution is about 5.51–7.37% depending on the location when (W = 5, P = 3) terms are used in the series. This error reduced to 2.32–3.11% and then to 0.54–0.72% for (W = 5, P = 5) and (W = 5, P = 7) terms in the series, respectively. Since, the FIT approach without homogenized BCs fails to obtain solutions at boundaries for Dirichlet BCs (see plots in Fig. 4) and the comparisons of the results between the cases without and with homogenized BCs could not be shown. Such comparisons are also not shown even away from boundaries because a reasonably large number of terms are required in the series for FIT without homog-

enization approach to match the solutions. As the approaches are related to eigenfunctions, the maximum number of terms are required to match the initial temperature distributions. In view of the above discussion, the problem for case 2 [f ðt Þ ¼ 1 þ sinðt Þ is solved using the FIT approach, with homogenized as well as without homogenized BCs. The temporal variation of temperature profile for case 1 and case 2 are shown in Fig. 5 at midpoint of each layer in direction h = p with homogenized BCs. Fig. 5 shows that both formulations (first and second) produce quite close results for constant as well as time dependent BCs. The temperature profile for time dependent BCs varies sinusoidally around the profile for constant BCs.

4.2. Specified heat flux (Neumann BCs) Fig. 6 represents the composite configuration that was used for a problem in which the heat flux was specified at the outer surface of the outermost layer and the corresponding solution was found

Table 2 Error in temperature obtained by using second formulation compared to the reference solution at location (r = r2 = 4, h = p) for different number of terms and different time instants. Combinations of terms {W, P} (W + 1)  P total number of terms (W = 5, P = 3); (5 + 1)  3 = 18 (W = 5, P = 5); (5 + 1)  5 = 30 (W = 5, P = 7); (5 + 1)  7 = 42

Error (r = r2, h = p, t, W, P) in solution for the second formulation (%) t = 0.5

t=1

t=2

ss (t ? 1)

7.37 3.11 0.72

6.11 2.57 0.60

5.61 2.36 0.55

5.51 2.32 0.54

730

P. Biswas et al. / International Journal of Heat and Mass Transfer 129 (2019) 721–734

Fig. 5. Comparison of temperature profiles at midpoint of each layer and h = p with homogenized BCs when constant [non-oscillating plot indicates f (t) = 1] and time dependent [oscillating plot indicates f (t) = 1 + sin (t)] temperatures are applied on the outer surface of the outermost layer. For first and second formulations the results are shown by solid and dashed lines, respectively.

Fig. 6. 2-D heat conduction in a 3-layer cylindrical composite subjected to time dependent flux which is a function of angle h and time t.

by Singh et al. [54] using the FIT without homogenized BCs. The same form of heat flux condition is used in the current illustration. It is to be mentioned that the temperature T i ðr; h; tÞ is represented in the units of q000 r0 =k1 while remaining variables have the same units as mentioned in the previous example. The analysis is carried out for the following cases: Case 1: f ðt Þ ¼ 1 and Case 2: f ðt Þ ¼ 1 þ sinðt Þ The temperature distributions, for the case of constant heat flux at the outer boundary of the outermost layer, are shown in Fig. 7 at

h = 0, p/4, p/2 and 3p/2. The FIT results without homogenized BCs (dotted lines) have error at/near boundaries (in enclosed region) compared to the FIT results with homogenized BCs. The error is due to mismatch between the non-zero applied heat flux (BCs) on the outer boundary and the zero slopes, at outer boundary, of the associated eigenfunctions corresponding to the homogeneous BCs. It is also can be seen from Fig. 7 that solutions with as well as without homogenized BCs match well everywhere except in the regions where applied heat flux is non-zero (for 0 < h < p) (see plots for h = p/4 and h = p/2 in Fig. 7). Moreover, the spurious oscillations are not present as seen in the case of Dirichlet BCs

P. Biswas et al. / International Journal of Heat and Mass Transfer 129 (2019) 721–734

731

Fig. 7. Temporal temperature distributions at different angular directions (h = 0, p/4, p/2 and 3p/2) using the FIT with as well as without (dotted lines) homogenized BCs when constant heat flux is applied at the outer surface of the outermost layer. In the case of homogenized BCs, results for first and second formulations are shown by solid and dashed lines, respectively.

(Fig. 3 and Fig. 4). The significantly better results and lack of spurious oscillations in case of Neumann BCs is due to the fact that there is a mismatch in the slopes of the eigenfunctions and the temperature at the boundary. The mismatch in slopes leads to relatively small mismatch in temperature. Moreover, with increasing number of terms, the abovementioned mismatch reduces. However, in case of Dirichlet BCs, the mismatch is in temperature itself which results in large oscillations and obviously does not change with number of terms at the boundary. Furthermore, the results for both first and second formulations are found to be identical for a small number of terms in the series. Fig. 8 shows radial temperature distribution at h = p/2 with different number of terms in the series solution. The FIT without homogenized BCs requires significantly higher number of terms compared to the case of homogenized BCs for reasonable convergence (dotted lines in enclosed region in plots). It is observed that the FIT results using Neumann BCs (without homogenized BCs) are identical with previously reported results [54]. The improvement in the FIT solution comes for the process homogenization of BCs is carried out prior to application of the integral transformation. For illustrative example with Neumann BCs, the solutions for first formulation with terms (W = 10, P = 10) are considered as the reference solutions. Because, the differences between solutions obtained by using homogenization of BCs (for both first as well as second formulations) in cases (W = 5, P = 5) and (W = 10, P = 10) are less than 0.2% at outer surface of the outermost layer (r = r3) (see Fig. 8) for h = p/2. The convergence in the FIT results without homogenized BCs for different combinations of terms {W, P} at r = r3 and h = p/2 during t = 5, 10, 15 and also for the steady state

(ss) are given in Table 3. Accuracy of the solutions keep increasing with increasing number of terms in the FIT series. The error of the solution is around 2.66–3.9% with the reference solution when (W = 5, P = 5) terms are considered in the series and it reduces to about 1.35–1.98% for (W = 10, P = 10) terms in the series. Table 3 shows that the number of eigenvalues in the h- and r-directions increases rapidly to obtain better accuracy. The effect of number of terms in r-direction is more prominent than terms in hdirection. One can see first and third row as well as second and fourth row to see change in error due to number of terms in rdirection. In this particular example, the error is about 2.66–3.9% for (W = 5, P = 5) and reduces to 1.39–2.05% for (W = 5, P = 10), whereas, the error is about 2.62–3.85% for (W = 10, P = 5) and reduces to 1.35–1.98% for (W = 10, P = 10). On the other hand, error due to h-direction is not that significant. As seen from the Table 3, the error is about 2.66–3.9% for (W = 5, P = 5) and reduces to 2.62– 3.85% for (W = 10, P = 5), whereas, the error is about 1.39–2.05% for (W = 5, P = 10) and reduces to 1.35–1.98% for (W = 10, P = 10). While improving the accuracy, the computational cost also increases resulting in a trade-off between them. The computational cost is drastically reduced because of the homogenization of BCs carried out prior to the application of integral transformation. In many real systems the outermost surfaces are subjected to time-varying heat flux. Therefore, the heat conduction problems are extended to show the effect of time dependent Neumann BCs. The problem is solved using the FIT for both with as well as without homogenized BCs. The relative error between the results without and with homogenized BCs is reasonably low for this specific problem. The temporal variations of the temperature pro-

732

P. Biswas et al. / International Journal of Heat and Mass Transfer 129 (2019) 721–734

Fig. 8. Effect of number of terms in series summation on temperature distributions at h = p/2, with as well as without homogenized BCs, during different time instants when constant heat flux is applied at the outer surface of the outermost layer.

Table 3 Error in temperature obtained by using FIT without homogenized BCs compared to the reference solution at location (r = r3 = 6 and h = p/2) for different number of terms and different time instants. Combinations of terms {W, P} (W + 1)  P total number of terms (W = 5, P = 5); (5 + 1)  5 = 30 (W = 10, P = 5); (10 + 1)  5 = 55 (W = 5, P = 10); (5 + 1)  10 = 60 (W = 10, P = 10); (10 + 1)  10 = 110

Error (r = r3, h = p/2, t, W, P) in the FIT without homogenized BCs (%) t=5

t = 10

t = 15

ss (t ? 1)

3.9 3.85 2.05 1.98

3.14 3.10 1.65 1.60

2.89 2.85 1.50 1.45

2.66 2.62 1.39 1.35

Fig. 9. Comparison of temperature profiles at midpoint of each layer and h = p/2 with homogenized BCs when constant [non-oscillating plot indicates f (t) = 1] and time dependent [oscillating plot indicates f (t) = 1 + sin (t)] heat fluxes are applied on the outer surface of the outermost layer. For first and second formulations the results are shown by solid and dashed lines, respectively.

file are shown in Fig. 9 at midpoint of each layer (i. e. r = 1.5, 3 and 5) only with homogenized BCs. Like in the case of Dirichlet BCs, in this example using time dependent Neumann BCs the temperature profile oscillates sinusoidally around the profile obtained for

constant BCs. It is seen from the comparative plots provided in Fig. 9, the mismatch between the results for first and second formulations is negligible using both constant as well as time dependent Neumann BCs cases.

P. Biswas et al. / International Journal of Heat and Mass Transfer 129 (2019) 721–734

5. Conclusion The homogenization of BCs is required in order to avoid spurious oscillations at the boundaries due to mismatch of BCs and the corresponding eigenfunctions. In present work two alternate auxiliary functions are proposed for 2D cylindrical polar coordinates which are subtracted from the temperature resulting in a problem for modified temperature which has homogenous BCs. It is seen through illustrative examples that the solutions without homogenization for Dirichlet BCs lead to severe oscillations in the solution making homogenization essential. With Neumann BCs the solution without homogenization does not have oscillations but requires relatively large number of terms when compared with solution with homogenization. Furthermore, the second formulation proposed here requires auxiliary functions only in two extreme layers (innermost and outermost). The auxiliary function in the second formulation is also very simple and does not need any evaluation. Therefore, second formulation can be considered as better than the first formulation. Through solution of illustrative example with Dirichlet BCs it is seen that the solution obtained using the second formulation requires only slightly more number of terms for same level of accuracy. However, it is pointed out here that since auxiliary function is defined only in two layers in case of second formulation (compared to this being defined in each layer for first formulation) the computational cost of first formulation will increase rapidly with increasing number of layers. In case of Neumann BCs the solutions for both formulations are nearly identical irrespective of number of terms used. Conflict of interest The authors declared that there is no conflict of interest. References [1] V. Vodicka, Eindimensionale Warmeleitung in Geschichteten Koorpern, Math. Nachrichten 14 (1) (1955) 47–55. [2] C.W. Tittle, Boundary value problems in composite media: quasi-orthogonal functions, J. Appl. Phys. 36 (4) (1965) 1486–1488. [3] G.P. Mulholland, M.H. Cobble, Diffusion through composite media, Int. J. Heat Mass Transf. 15 (1) (Jan. 1972) 147–160. [4] F. de Monte, Transient heat conduction in one-dimensional composite slab. A ‘natural’ analytic approach, Int. J. Heat Mass Transf. 43 (19) (2000) 3607–3619. [5] F. de Monte, An analytic approach to the unsteady heat conduction processes in one-dimensional composite media, Int. J. Heat Mass Transf. 45 (6) (2002) 1333–1343. [6] F. de Monte, Unsteady heat conduction in two-dimensional two slab-shaped regions. Exact closed-form solution and results, Int. J. Heat Mass Transf. 46 (8) (2003) 1455–1469. [7] F. de Monte, Transverse eigenproblem of steady-state heat conduction for multi-dimensional two-layered slabs with automatic computation of eigenvalues, Int. J. Heat Mass Transf. 47 (2) (2004) 191–201. [8] F. de Monte, Multi-layer transient heat conduction using transition time scales, Int. J. Therm. Sci. 45 (9) (2006) 882–892. [9] M.D. Mikhailov, M.N. Özisßik, N.L. Vulchanov, Diffusion in composite layers with automatic solution of the eigenvalue problem, Int. J. Heat Mass Transf. 26 (8) (1983) 1131–1141. [10] K.A. Antonopoulos, C. Tzivanidis, Analytical solution of boundary value problems of heat conduction in composite regions with arbitrary convection boundary conditions, Acta Mech. 118 (1–4) (1996) 65–78. [11] Y. Sun, I.S. Wichman, On transient heat conduction in a one-dimensional composite slab, Int. J. Heat Mass Transf. 47 (6–7) (2004) 1555–1559. [12] R.M. Shvets, O.I. Yatskiv, Extension of the method of eigenfunctions to the boundary-value problems of mechanical diffusion for multilayer bodies with interlayers, J. Math. Sci. 107 (1) (2001) 3691–3696. [13] C.P. Naveira-Cotta, R.M. Cotta, H.R.B. Orlande, O. Fudym, Eigenfunction expansions for transient diffusion in heterogeneous media, Int. J. Heat Mass Transf. 52 (21–22) (2009) 5029–5039. [14] G.A. Surkov, The use of finite integral transforms to solve problems of unsteady heat conduction in hollow cylinders with moving internal boundaries, J. Eng. Phys. 8 (3) (1965) 261–264. [15] R.M. Cotta, Hybrid numerical/analytical approach to nonlinear diffusion problems, Numer. Heat Transf. Part B Fundam. 17 (2) (1990) 217–226.

733

[16] P.R. Johnston, Diffusion in composite media: Solution with simple eigenvalues and eigenfunctions, Math. Comput. Model. 15 (10) (1991) 115–123. [17] R.M. Cotta, Integral Transforms in Computational Heat and Fluid Flow, CRC Press, Boca Raton, 1993. [18] R.M. Cotta, M.D. Mikhailov, Integral transform method, Appl. Math. Model. 17 (3) (Mar. 1993) 156–161. [19] S.C. Huang, Y.P. Chang, Heat conduction in unsteady, periodic, and steady states in laminated composites, J. Heat Transfer 102 (4) (1980) 742. [20] J.V. Beck, K.D. Cole, A. Haji-Sheikh, B. Litkouhi, Heat Conduction Using Green’s Function, first ed., Taylor & Francis, Washington, DC, 1992. [21] Z.-G. Feng, E.E. Michaelides, The use of modified Green’s functions in unsteady heat transfer, Int. J. Heat Mass Transf. 40 (12) (1997) 2997–3002. [22] A. Haji-Sheikh, J.V. Beck, Temperature solution in multi-dimensional multilayer bodies, Int. J. Heat Mass Transf. 45 (9) (2002) 1865–1877. [23] K.D. Cole, J.V. Beck, A. Haji-Sheikh, B. Litkouhi, Heat Conduction Using Green’s Functions, second ed., Taylor & Francis, 2010. [24] F.J. Rizzo, D.J. Shippy, A method of solution for certain problems of transient heat conduction, AIAA J. 8 (11) (1970) 2004–2009. [25] Z.W. Zhou, Analytical solution for transient heat conduction in hollow cylinders containing well-stirred fluid with uniform heat sink, Int. J. Heat Mass Transf. 38 (15) (1995) 2915–2919. [26] M.N. Özisik, Heat Conduction, second ed., John Wiley & Sons Inc, New York, 1993. [27] A. Haji-Sheikh, J.V. Beck, Green’s function partitioning in Galerkin-based integral solution of the diffusion equation, J. Heat Transfer 112 (1) (1990) 28. [28] H. Salt, Transient conduction in a two-dimensional composite slab—I. Theoretical development of temperature modes, Int. J. Heat Mass Transf. 26 (11) (1983) 1611–1616. [29] H. Salt, Transient conduction in a two-dimensional composite slab—II. Physical interpretation of temperature modes, Int. J. Heat Mass Transf. 26 (11) (1983) 1617–1623. [30] M.D. Mikhailov, M.N. Özisßik, Transient conduction in a three-dimensional composite slab, Int. J. Heat Mass Transf. 29 (2) (1986) 340–342. [31] I. Dülk, T. Kovácsházy, Steady-state heat conduction in multilayer bodies: An analytical solution and simplification of the eigenvalue problem, Int. J. Heat Mass Transf. 67 (2013) 784–797. [32] I. Dülk, T. Kovácsházy, A method for computing the analytical solution of the steady-state heat equation in multilayered media, J. Heat Transfer 136 (9) (2014) 91303. [33] I. Dülk, Temperature solution in composites by domain decomposition, Compos. Struct. 147 (2016) 54–64. [34] M. Bakhtiari, K. Daneshjou, R. Alibakhshi, A new approach to thermal analysis of a multilayered cylindrical structure with imperfect bonds and internal heat source, J. Heat Transfer 138 (12) (2016) 122801. [35] M. Norouzi, H. Rahmani, A.K. Birjandi, A.A. Joneidi, A general exact analytical solution for anisotropic non-axisymmetric heat conduction in composite cylindrical shells, Int. J. Heat Mass Transf. 93 (2016) 41–56. [36] B. Yang, S. Liu, Closed-form analytical solutions of transient heat conduction in hollow composite cylinders with any number of layers, Int. J. Heat Mass Transf. 108 (2017) 907–917. [37] A. Amiri Delouei, M. Norouzi, Exact analytical solution for unsteady heat conduction in fiber-reinforced spherical composites under the general boundary conditions, J. Heat Transfer 137 (10) (2015) 101701. [38] M.H. Kayhani, M. Norouzi, A. Amiri Delouei, A general analytical solution for heat conduction in cylindrical multilayer composite laminates, Int. J. Therm. Sci. 52 (1) (2012) 73–82. [39] A. Amiri Delouei, M.H. Kayhani, M. Norouzi, Exact analytical solution of unsteady axi-symmetric conductive heat transfer in cylindrical orthotropic composite laminates, Int. J. Heat Mass Transf. 55 (15–16) (2012) 4427–4436. [40] M. Norouzi, A. Amiri Delouei, M. Seilsepour, A general exact solution for heat conduction in multilayer spherical composite laminates, Compos. Struct. 106 (2013) 288–295. [41] H.M. Wang, An effective approach for transient thermal analysis in a functionally graded hollow cylinder, Int. J. Heat Mass Transf. 67 (Dec. 2013) 499–505. [42] K. Daneshjou, M. Bakhtiari, R. Alibakhshi, M. Fakoor, Transient thermal analysis in 2D orthotropic FG hollow cylinder with heat source, Int. J. Heat Mass Transf. 89 (2015) 977–984. [43] M. Norouzi, H. Rahmani, On exact solutions for anisotropic heat conduction in composite conical shells, Int. J. Therm. Sci. 94 (2015) 110–125. [44] M. Norouzi, H. Rahmani, An exact analysis for transient anisotropic heat conduction in truncated composite conical shells, Appl. Therm. Eng. 124 (2017) 422–431. [45] E.J. Corrêa, R.M. Cotta, H.R.B. Orlande, On the reduction of computational costs in eigenfunction expansions of multidimensional diffusion problems, Int. J. Numer. Methods Heat Fluid Flow 7 (7) (1997) 675–695. [46] S. Singh, P.K. Jain, Rizwan-uddin, Analytical solution to transient heat conduction in polar coordinates with multiple layers in radial direction, Int. J. Therm. Sci. 47 (3) (2008) 261–273. [47] P.K. Jain, S. Singh, Rizwan-uddin, Analytical solution to transient asymmetric heat conduction in a multilayer annulus, J. Heat Transfer 131 (1) (2009) 11304. [48] P.K. Jain, S. Singh, Rizwan-uddin, An exact analytical solution for twodimensional, unsteady, multilayer heat conduction in spherical coordinates, Int. J. Heat Mass Transf. 53 (9–10) (2010) 2133–2142. [49] S. Singh, P.K. Jain, Rizwan-uddin, Analytical solution for three-dimensional, unsteady heat conduction in a multilayer sphere, J. Heat Transfer 138 (10) (2016) 101301.

734

P. Biswas et al. / International Journal of Heat and Mass Transfer 129 (2019) 721–734

[50] M. Li, A.C.K. Lai, Analytical solution to heat conduction in finite hollow composite cylinders with a general boundary condition, Int. J. Heat Mass Transf. 60 (May 2013) 549–556. [51] Y. Yener, M.N. Özisik, On the solution of unsteady heat conduction in multiregion media with time-dependent heat transfer coefficient, in: Proceedings of the 5th International Heat Transfer Conference, 1974, pp. 188–192. [52] R.M. Cotta, Benchmark results in computational heat and fluid flow: The integral transform method, Int. J. Heat Mass Transf. 37 (SUPPL. 1) (1994) 381–393.

[53] L.A. Sphaier, R.M. Cotta, Integral transform analysis of multidimensional eigenvalue problems with irregular domains, Numer. Heat Transf. Part B Fundam. 38 (2) (2000) 157–175. [54] S. Singh, P.K. Jain, Rizwan-uddin, Finite integral transform method to solve asymmetric heat conduction in a multilayer annulus with time-dependent boundary conditions, Nucl. Eng. Des. 241 (1) (2011) 144–154. [55] H.M. Wang, C.B. Liu, Analytical solution of two-dimensional transient heat conduction in fiber-reinforced cylindrical composites, Int. J. Therm. Sci. 69 (2013) 43–52.