# Heat Conduction in a Multilayer Sphere

Introduction | Solve the PDE Heat Transfer Model |

Domain and Mesh Generation | Post Processing and Visualization |

Heat Transfer Model | Nomenclature |

Initial and Boundary Conditions | References |

## Introduction

Devices constructed by layers of multiple materials have been used in many industrial applications such as chemical, mechanical, and aerospace engineering. One important feature of multilayered materials lies in the possibility to combine several thermal properties. Performing numerical simulations and obtaining the temperature distribution of such compounds are important because they give insights of how the real system will behave.

This model is based on [Singh et al., 2016] in which a 3D heat transfer transport through a multilayered sphere was analysed. The problem statement is to perform a thermal analysis and find the temperature field distribution T(t,x,y,z) at time t of n concentric hollow spherical layers in thermal contact. A 2D cross section of the geometry is shown for n=3.

Each layer of the sphere has specific material properties, for example, the i th layer will have a thermal conductivity K_{i}, mass density ρ_{i} and specific heat capacity Cp_{i}. For this simulation we will limit ourselves to three layers - three concentric hollow spheres - and we have i =1,2,3. As mentioned in [Singh et al., 2016] two different types of boundary conditions are imposed. One at the external face and a second type at the internal face of the multi layer sphere setup. In this model the external face of the hollow sphere is subjected to a non linear heat source Q_{external}(θ,ϕ) at all times and the internal face of the first layer will be set to a fixed temperature of T_{internal}=0[ºC]. The simulation starts with an initial temperature of T_{initial}=0[ºC] for all layers of the sphere. As heat transfers between layers, the temperature change of each layer will depend, as expected, on its thermal properties.

Refer to the information provided in "Heat Transfer" for a more general theoretical background on heat transfer analysis.

## Domain and Mesh Generation

The following plot shows a 2-dimensional cross sectional cut through the 3 concentric spherical layers. Indicated are the radii r_{i} and the model parameter labels of each layer κ_{i}, ρ_{i} and Cp_{i}.

Next, we combine the faces. This will result in an object with 4 regions, 0<r<r_{0}, r_{0}<r<r_{1},... The multilayered sphere is hollow in the region 0<r<r_{0}. This requirement can be specified during the mesh generation with the ToElementMesh option RegionHoles.

The three layers will be labeled with three markers 1, 2 and 3. For more information on markers and their generation in meshes see the Element Mesh Generation tutorial.

## Heat Transfer Model

A transient heat transfer equation (1) will be used to determine the temperature field distribution T(t,x,y,z):

For this PDE model, two boundary conditions will be applied to the inner face (r=r_{0}) and the external face (r=r_{3}) of the multilayer sphere. On the inner face, a surface boundary condition of constant temperature equal to T_{internal}=0[ºC] will be imposed. On the external face, a nonlinear heat flux source of value Q(θ,ϕ)=q_{0}θ^{2}(π-θ)^{2}ϕ^{2}(π-ϕ)^{2} will be imposed, where θ is the polar angle and ϕ is the azimuthal angle in spherical coordinates.

Also, see this note about how to set up of computationally efficient PDE coefficients.

## Initial and Boundary Conditions

For this PDE model the internal face of layer 1 will be set to a temperature surface boundary condition of T_{internal}=0[ºC]. The external face of the layer 3 will be subjected to a nonuniform time independent heat flux condition Q(θ,ϕ) that depends on the polar angle θ and the azimuthal angle ϕ. The problem description in [Singh et al., 2016] is formulated in spherical coordinates. To apply the boundary conditions, originally defined in spherical coordinates, a coordinate transformation from spherical coordinates to cartesian coordinates is performed.

_{0}, the heat source magnitude, and convert the original boundary conditions to cartesian coordinates.

First, a DirichletCondition is set up with the help of the HeatTemperatureCondition function.

The NeumannValue on the outside of the sphere is created with the help of the HeatFluxValue function. The heat-flux source is applied on the external face of layer 3.

The following plot shows the heat flux boundary condition and the coordinate system used to formulate the problem. The heat flux has it's peak value q_{0} at y=r_{3} and it decreases to 0 in the region y<0.

## Solve the PDE Heat Transfer Model

To analyze the heat conduction between the three layers, we solve the model from t_{initial}=0s to t_{final}=50 s.

## Post Processing and Visualization

We visualize 1D, 2D and 3D cuts of the temperature field distribution T(t,x,y,z) of the multilayered sphere.

### 1 Dimensional Cut: Temperature distribution in radial direction

We plot the transient temperature distribution solutions in the radial direction along θ=π/2 and the azimuthal (ϕ) directions: a) ϕ=0, b) ϕ=π/2, and c) ϕ=π.

The linear section ϕ=0& θ=π/2 corresponds to x>0& y=0, ϕ=π/2 & θ=π/2 corresponds to x=0&y>0, and ϕ=π & θ=π/2 corresponds to x=0& y<0.

### 2 Dimensional Cut: Temperature contour plot distribution

Visualize the temperature field contours in the z=0 cross section of the sphere. The temperature distribution of the whole system increases with time. For each time in t∈{5,10,20,30,45,50} we define a list of nine temperature values to plot as contours.

### 3 Dimensional Cut: Temperature density plot animation

We perform a time animation of the temperature field distribution in a 3D section of the sphere. When the heat source is applied from the region y>0, the temperature starts to increase in the y direction as can be seen in the animation below. Since the relationship between thermal conductivities of each layer is κ_{3}>κ_{2}>κ_{1}, layer 3 (outer layer) transports heat and increases its temperature faster than the internal layer as can be seen in the density plot animation.

A limit in the plot range of the SliceContourPlot3D is used and the reason is that if the center sphere of the slice surface "CenterCutSphere" is exactly the same sphere as the region of the mesh it is plotted over, part of the surface will not be displayed. Therefore, the range of the plot must be smaller than the region of the mesh.

The 3D animations that follow can require an above average amount of disk space. To reduce the disk space this notebook requires the function Rasterize is called on the visualization. The disadvantage of this is that the quality of the animation will not be as sharp as without it. To obtain high quality graphics remove the call to Rasterize.

## Nomenclature

## References

1. Suneet, S., Jain, P. K., & Uddin, R. Analytical Solution for Three-Dimensional, Unsteady Heat Conduction in a Multilayer Sphere. ASME. Journal of Heat Transfer, 2016, 138(10), 101301. doi: 10.1115/1.4033536.