multilayer materials are composite media composed of several layers . because of the additional benefit of combining various mechanical , physical , and thermal properties of different substances , a construction using multilayer elements is of interest . multilayer materials are used in semicircular fiber insulated heaters , multilayer insulation materials , and nuclear fuel rods . multilayer transient heat conduction finds applications in thermodynamics , fuel cells , and electrochemical reactors . the layered sphere is utilized to investigate the thermal properties of composite media by assuming embedded spherical particles in the composite matrix . for solving the problems of multilayer transient heat conduction , the same methods which are used in solving problems of single layer transient heat conduction are applied . analytical methods are advantageous over numerical methods in two ways : ( 1 ) analytical solutions can be used as benchmark to examine and actually confirm numerical algorithms ; ( 2 ) compared to a discrete numerical solution , the mathematical form of an analytical solution can provide better insight . it should also be mentioned that the analytical methods applied to multilayer transient conduction are analogous to those used in the single - layer transient heat conduction . these analytical methods include green 's function method , the laplace transform , separation of variables , and eigenfunction expansion method . for instance , salt   solved the transient heat conduction problem in a two - dimensional composite slab using an orthogonal eigenfunction expansion technique . mikhailov and ziik , using the orthogonal expansion approach , solved the problem of transient three - dimensional heat conduction in a composite cartesian medium . haji - sheikh and beck   used green 's function method to obtain temperature distribution in a three - dimensional two - layer orthotropic slab . de monte [ 4 , 5 ] applied the eigenfunction expansion method to obtain the transient temperature distribution for the heat conduction in a two - dimensional two - layer isotropic slab with homogenous boundary conditions . lu et al .   and lu and viljanen   combined separation of variables and laplace transforms to solve the transient conduction in the two - dimensional cylindrical and spherical media . singh et al . [ 8 , 9 ] and jain et al . [ 10 , 11 ] used the combination of separation of variables and eigenfunction expansion methods to solve the two - dimensional multilayer transient heat conduction in spherical coordinates . recently , dalir and nourazar   used the eigenfunction expansion method to solve the problem of three - dimensional transient heat conduction in a multilayer cylinder . singh et al . [ 8 , 9 ] and jain et al . [ 10 , 11 ] have studied 2d multilayer transient conduction problems in spherical and cylindrical coordinates . they have obtained analytical solutions for 2d multilayer transient heat conduction in spherical coordinates , in polar coordinates with multiple layers in the radial direction , and in a multilayer annulus . they have used the method of partial solutions to obtain the temperature distributions . in the method of partial solutions , the nonhomogeneous transient problem is split into two subproblems : a nonhomogeneous steady - state subproblem and a homogeneous transient subproblem . then , the eigenfunction expansion method is used to solve the nonhomogeneous steady - state subproblem and the method of separation of variables is used to solve the homogeneous transient subproblem . the literature survey for the exact analytical solution for 3d transient heat conduction in multilayered sphere demonstrates that such a solution has not , so far , been developed . thus , in the present paper , using the eigenfunction expansion method , an analytical triple - series solution for transient heat conduction in the 3d spherical coordinates for radial multilayer domain with spatially nonuniform and time - dependent internal heat sources is obtained . homogenous boundary conditions of the first or second kind can be applied on surfaces of  = constant and  = constant . however , nonhomogeneous boundary conditions of the third kind ( convection )   are used in the r - direction . first , the problem is a boundary - value problem of conduction in spherical ( r-- coordinates ) or part - spherical multilayer geometries . second , volumetric internal heat sources of nonuniform and time - dependent ( r ,  ,  and t - dependent ) types are present . third , on the inner and outer radial boundaries , nonhomogeneous boundary conditions of any kind can be used but , on the boundary surfaces in the  and -directions , only the first or second kind of homogeneous boundary condition can be applied . a n - layer composite spherical slab ( r0  r  rn , 0     ,  and  0     ) is considered . all the layers have perfect thermal contact and are presumed to be isotropic in thermal properties . i and ki are the temperature independent thermal diffusivity and thermal conductivity of the ith layer . at t = 0 , the ith layer is at a specified temperature fi(r ,  ,  ) and time dependent heat sources gi(r ,  ,  , t ) are switched on in each radial layer . for t > 0 , homogenous boundary conditions of first or second kind are applied to the angular surfaces of  = 0 and  =  and azimuthal surfaces of  = 0 and  = . for the inner ( i = 1 , r = r0 ) and the outer ( i = n , r = rn ) radial surfaces , all three kinds of boundary conditions are applicable . the governing differential equation of the 3d transient conduction in a multilayered sphere with heat sources is as follows : ( 1)2tir2 + 2rtir+1r22ti2+cot  r2ti+1r2sin22ti2  + gir,,,tki=1itit , ti = tir,,,t ,  r0rrn , ri1rri ,  1in,0 ,  <,0 ,  <2,t0 . the boundary conditions are as follows.(i)inner surface of 1st layer ( i = 1 ) : ( 2)aint1rr0,,,t+bint1r0,,,t = cin.(ii)outer surface of nth layer ( i = n ) : ( 3)aouttnrrn,,,t+bouttnrn,,,t = cout.(iii) =  surface ( i = 1,2 ,  , n ) : ( 4)tir,=,,t=0  or  tir,=,,t=0,llllllllllllllllllllllllllli=1,  ,n.(iv) = 0 surface ( i = 1,2 ,  , n ) : ( 5)tir,,=0,t=0  or  ti(r,,=0,t)=0,llllllllllllllllllllllllllli=1,  ,n.(v) =  surface ( i = 1,2 ,  , n ) : ( 6)ti(r,,=,t)=0  or  ti(r,,=,t)=0,llllllllllllllllllllllllllli=1,  ,n.(vi)inner interface of the ith layer ( i = 2 ,  , n ) : ( 7)tiri1,,,t = ti1ri1,,,t  i=2,  ,n , kitirri1,,,t = ki1ti1rri1,,,t  i=2,  ,n.(vii)outer interface of the ith layer ( i = 1 , inner surface of 1st layer ( i = 1 ) : ( 2)aint1rr0,,,t+bint1r0,,,t = cin . outer surface of nth layer ( i = n ) : ( 3)aouttnrrn,,,t+bouttnrn,,,t = cout . =  surface ( i = 1,2 ,  , n ) : ( 4)tir,=,,t=0  or  tir,=,,t=0,llllllllllllllllllllllllllli=1,  ,n . inner interface of the ith layer ( i = 2 ,  , n ) : ( 7)tiri1,,,t = ti1ri1,,,t  i=2,  ,n , kitirri1,,,t = ki1ti1rri1,,,t  i=2,  ,n . outer interface of the ith layer ( i = 1 ,  , n  1 ) : ( 8)tiri,,,t = ti+1ri,,,t  i=1, it is worth mentioning that , at r = r0 and r = rn , boundary conditions of first , second , or third kind are applied by appropriate selection of the coefficients in ( 2 ) and ( 3 ) . it should also be mentioned that zero inner radius ( r0 = 0 ) for multilayered sphere is modeled by assigning zero values to bin and cin in ( 2 ) . the eigenfunction expansion method is used to solve the problem . in the eigenfunction expansion method , first , by using the associated eigenvalue problem (  =  ) , the eigenfunctions are attained at every spatial direction of the problem . afterward , the dependent variable and the available nonhomogeneity in the governing differential equation of the problem are separately written as series expansions of the eigenfunctions . in heat conduction problems , the series expansions are then substituted into the differential equation . by performing some mathematical manipulations , an ordinary differential equation ( ode ) the solution of the problem is completed by solving this ode , which is a first order ode in the case of heat conduction problems . as stated before , the method of partial solutions was used by jain and singh   for solving 2d transient heat conduction problems , the reason being that the heat source is independent of time . however , the method of partial solutions can not be used for solving the present 3d transient heat conduction problem because the heat source depends on time . thus , due to time dependence of the heat source , the partial solution of the steady - state subproblem can not include the heat source term and the partial solutions method can not be used . therefore , to the best knowledge of the authors , the most efficient tool for solving the 3d heat conduction problem of the present paper is the eigenfunction expansion method . for the transient problem of present paper , the associated eigenvalue problem is written as follows : ( 11)2i=2i2ir2 + 2rir+1r22i2llllllllllllllilll+cot  r2i+1r2sin22i2llllllllllllll=2i . using the method of separation of variables , ( 11 ) is solved as follows : ( 12)ir,,=ririi,(13)riiiriii+2rriiiriii+1r2riiiriii+cot  r2riiiriii  + 1r2sin2riiiriii=2riiiriii,(14)riri+2rriri+1r2ii+cot  r2ii+1r2sin2ii=2,(15)riri+2rriri+1r2ii+cot  r2ii=1r2sin2ii2,(16)sin2r2riri+2rriri+ii+cot  ii+2r2  = ii=+2,(17)i+im2i=0im=c1sinim+c2cosim,(18)r2riri+2rriri+2r2=iicot  ii+2sin2=+2,(19)r2ri+2rri+(ipn2r2n2)ri=0  ripn(r)=1rc3jn+0.5ipnrllllllllllllllllllllllllllllllllllllll+  c4yn+0.5ipnr , by using the following change of variable : ( 20)=cos12=sin2 the first and second derivatives of  with respect to  are obtained as follows : ( 21)=did=diddd=sindid,i=d2id2=dddid=ddsindid=cosdidsindddid=sin2d2id2cosdid. substituting ( 21 ) in ( 19 ) results in the following : ( 22)i+(cot  )i+22sin2i=0  sin2d2id2cosdid  + cossinsindid+22sin2i=0  (12)d2id22did+n2im212i=0 . if n = n(n + 1 ) , ( 22 ) is the associated legendre equation . its solution is written as follows : ( 23)in=c5pinim()+c6qinim()in()=c5pinim(cos)+c6qinimcosqinimcos0  =  qinim1  =  c6  =  0in=pinimcos. the problem eigenvalues are imnp = ipn + im . it should be stated that the heat fluxes continuity at the interfaces of the radial layers gives the following : ( 24)im=m ,  im=m,ipn=1pn1i . the eigenfunctions ripn(r ) , in( ) , and im( ) in r- , - , and -directions are derived by applying the boundary conditions in each direction in the following equation : ( 25)ripnr=1rc3jnn  +  1  +  0.5ipnr+c4ynn  +  1  +  0.5ipnr,in=pinmcos,m=c1sinm+c2cosm. it is assumed that the solution of the problem is in the form of a triple - series expansion of the derived eigenfunctions as follows : ( 26)tir,,,t=m  =  1n  =  1p  =  1timnptripnrinm. the heat source term is also written as a triple - series expansion of the eigenfunctions such that ( 27)gir,,,t=m  =  1n  =  1p  =  1gimnptripnrinm , where the coefficient gimnp(t ) is obtained by the use of the orthogonality property as follows : ( 28)gimnpt=00ri    1rigi(r,,,t)r2llllllllllllllllllllripnrinmdr  d  d0  00ri    1rir2ripn2rin2lllllllllllllllllllllllllm2dr  d  d01 . substitution of ( 26 ) and ( 27 ) in ( 1 ) results in the following : ( 29)timnp(t)ripn(r)in()im()+2rtimnp(t)ripn(r)in()im()+1r2timnp(t)ripn(r)in()im( )  + cot  r2timnp(t)ripn(r)in()im()+1r2sin2timnptripnrinim  + 1kigimnptripnrinim=1itimnptripnrinim,(30)dtimnptdt+iripnrripnr+2rripnrripnr+1r2inin+cot  r2inin+1r2sin2imim=  imnptimnpt=ikigimnpt . equation ( 30 ) is a first - order nonhomogeneous ode and has the following solution : ( 31)timnpt=ikieimnpt  =  0  =  tgimnpeimnpd+ai1eimnpt . application of the initial condition ( 10 ) on ( 26 ) gives the following : ( 32)tir,,,t=0  = fi(r,, )  = m  =  1n  =  1p  =  1timnp0ripnrinim. the coefficient ai1 of ( 31 ) is found by the use of the orthogonality property for obtained eigenfunctions as follows : ( 33)ai1=timnp(0)=00ri1rifir,,r2ripnrinimdr  d  d00ri    1rir2ripn2rin2im2dr  d  d. the solution of differential equation ( 1 ) having ( 2 ) to ( 9 ) as boundary conditions and ( 10 ) as initial condition is ( 26 ) with ( 31 ) as the coefficients . we consider a three - layer quarter - sphere ( 0  r  r3 ,  0    /2 ,  and  0     ) which is initially ( t = 0 ) at uniform unit temperature . for time t > 0 , thermal convection occurs , from the outer radial surface at r = r3 at ( zero ) ambient temperature . however , the  = 0 ,  = /2 ,  = 0 , and  =  surfaces are at uniform and constant zero temperatures . these boundary conditions lead to the following : ain = 1 , aout = k3 , bin = 0 , bout = h , cin = 0 , and cout = 0 . additionally , the uniformly distributed heat source gi , i = 1 ,  , 3 , is turned on in each layer at t = 0 . the governing differential equation for the 3d transient heat conduction with heat sources in this three - layer quarter - spherical region is as follows : ( 34)2tir2 + 2rtir+1r22ti2+cot  r2ti+1r2sin22ti2+giki  = 1itit , ti = tir,,,t ,  0rr3,ri1rri ,  1i3,02,0,t0 . the boundary conditions have the following forms : ( 35)t1r0,,,t=0,(36)k3t3rr3,,,t+ht3r3,,,t=0,(37)tir,2,,t=0,(38)tir,,0,t=0,(39)ti(r,,,t)=0.(i)inner interface of the ith layer ( i = 2,3 ) : ( 40)tiri    1,,,t = ti    1ri    1,,,t , kitirri    1,,,t = ki1ti    1rri    1,,,t.(ii)outer interface of the ith layer ( i = 1,2 ) : ( 41)tiri,,,t = ti  +  1ri,,,t,(42)kitirri,,,t = ki  +  1ti  +  1rri,,,t.the initial condition is as follows : ( 43)ti(r,,,0)=1 ,  1i3 . according to ( 25 ) , by the use of the eigenfunction expansion method , the following solutions in z- , r- , and -directions are obtained from the associated eigenvalue problem : ( 44)m=c1sinm+c2cosm,in=pinmcos,ripnr=1rc3jnn  +  1  +  0.5ipnr+c4ynn  +  1  +  0.5ipnr . the eigenvalues are imnp = ipn + im . the eigenfunctions ripn(r ) , in( ) , and m( ) in the r , - and -directions are obtained by applying the relevant boundary conditions in each direction . application of the boundary conditions in the -direction due to ( 44 ) gives the following : ( 46)m=c1sinm+c2cosmm(0)=0c2=0  m()=c1sinmm()=0c1sinm=0  c1    0sinm=0=sinm. then the -direction eigenvalues and eigenfunction are obtained as follows : ( 47)m = m ,  m=1,2,  m=sinm. using the -direction boundary condition , ( 37 ) , on  in ( 44 ) results in the following relation : ( 48)in=pinmcosin/2  =  0pinmcos2=0pinm(0)=0 , where pin(0 ) = 0 is only satisfied when n are odd integers ; that is , n = 1 , 3 , 5 ,  . thus the -direction eigenvalues and the eigenfunction are as follows : ( 49)n=1,3,5,  in=pinmcos. applying the r - direction boundary conditions , that is , ( 35 ) and ( 36 ) , gives the following : ( 50)ripnr  = 1rc3jnn  +  1  +  0.5ipnr+c4ynn  +  1  +  0.5ipnr , lllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllli=1,2,3r(0)=0r(0)=finiteyn(n  +  1 )  +  0.5(0 )  =  c4=0  ripn(r)=c31rjn(n  +  1 )  +  0.5(ipnr)kr(r3)+hr(r3)=0  kc31r3jn(n  +  1 )  +  0.5(ipnr3 )    kc312r3r3jnn  +  1  +  0.5(ipnr3 )  +  hc31r3jn(n  +  1 )  +  0.5(ipnr3)=0c3    01r3hk2r3kjn(n  +  1 )  +  0.5(ipnr3)lllllllllllllllll+hk2r3jn(n  +  1 )  +  0.5(ipnr3)=0kjn(n  +  1 )  +  0.5(ipnr3 )  + hk2r3jnn  +  1  +  0.5(ipnr3)=0 . thus the r - direction eigencondition and eigenfunction are derived as ( m = m ) : ( 51)kjn(n  +  1 )  +  0.5ipnr3+hk2r3jn(n  +  1 )  +  0.5ipnr3=0,ripnr=1rjnn  +  1  +  0.5ipnr . according to ( 28 ) to ( 33 ) , the coefficients gimnp(t ) , ai1 , and imnp are obtained as follows : ( 52)gimnpt=00ri    1rigir21rjn  +  0.5ipnrpinmllllllllllllllllllllllsinmdr  d  d0  00ri1rirjn  +  0.52ipnrpinm2lllllllllllllllllllllllllsin2mdr  d  d01=0giri1rir3/2jn  +  0.5ipnrdr0pinmdllllll0sinmd  ri1rirjn  +  0.52ipnrdr0pinm2dlllllllll0sin2mdri1ri1=1mcosm  =  0  =  giinmrjm  +  1inmrr  =  ri    1r  =  rillllll  1mcosm  =  0  =  1mcosm  =  0  =    inmr22jm2inmr+jm+12inmrr  =  ri    1r  =  rillllllllll2l2inmr221=11lgiinmrjm  +  1inmrr  =  ri    1r  =  ri1m11mllllllll11l  inmr22jm2inmr+jm+12inmrr  =  ri    1r  =  rilllllllllll4inmr221gimnp(t)=gimnp = ai1gi,imnp=i  ripn(r)ripn(r)+1rripn(r)ripn(r)+1r2in()in()+im()im()=ijmipnrjmipnr+1rjmipnrjmipnrm2r2m2 . using the form of a triple - series expansion of the obtained eigenfunctions , that is , ( 26 ) , the solution is as follows : ( 53)tir,,,t=m  =  1n  =  1p  =  1timnpt1rjn  +  0.5llllllllllllllllllllllllllllllllllllllipnrpinmsinm , where the coefficient timnp(t ) is attained as follows : ( 54)timnpt=ikieimnptgimnp1imnp(eimnpt1)+ai1eimnpt=igimnpkiimnp+1giikiimnpgimnpeimnpt . therefore , the solution of ( 34 ) having ( 35 ) to ( 42 ) as boundary conditions and ( 43 ) as initial condition is ( 53 ) with ( 54 ) as the coefficients . inner interface of the ith layer ( i = 2,3 ) : ( 40)tiri    1,,,t = ti    1ri    1,,,t , kitirri    1,,,t = ki1ti    1rri    1,,,t . outer interface of the ith layer ( i = 1,2 ) : ( 41)tiri,,,t = ti  +  1ri,,,t,(42)kitirri,,,t = ki  +  1ti  +  1rri,,,t . the exact analytical solution , that is , transient temperature distribution , is derived for the 3d transient heat conduction problem in a multilayered sphere using eigenfunction expansion method . third kind nonhomogeneous boundary conditions are applied in the radial direction but the first or second kind homogenous boundary conditions are used in the angular and azimuthal directions . the heat conduction in a three - layer quarter - sphere is solved as a case study problem and the temperature distribution is found .
an exact analytical solution is obtained for the problem of three - dimensional transient heat conduction in the multilayered sphere . 
 the sphere has multiple layers in the radial direction and , in each layer , time - dependent and spatially nonuniform volumetric internal heat sources are considered . to obtain the temperature distribution , 
 the eigenfunction expansion method is used . 
 an arbitrary combination of homogenous boundary condition of the first or second kind can be applied in the angular and azimuthal directions . 
 nevertheless , solution is valid for nonhomogeneous boundary conditions of the third kind ( convection ) in the radial direction . a case study problem for the three - layer quarter - spherical region 
 is solved and the results are discussed .