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 Ki, mass density ρi and specific heat capacity Cpi. 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 Qexternal(θ,ϕ) at all times and the internal face of the first layer will be set to a fixed temperature of Tinternal=0[ºC]. The simulation starts with an initial temperature of Tinitial=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 ri and the model parameter labels of each layer κi, ρi and Cpi.
Next, we combine the faces. This will result in an object with 4 regions, 0<r<r0, r0<r<r1,... The multilayered sphere is hollow in the region 0<r<r0. 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=r0) and the external face (r=r3) of the multilayer sphere. On the inner face, a surface boundary condition of constant temperature equal to Tinternal=0[ºC] will be imposed. On the external face, a nonlinear heat flux source of value Q(θ,ϕ)=q0θ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 Tinternal=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.
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 q0 at y=r3 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 tinitial=0s to tfinal=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.