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.

Load the finite element package and OpenCascadeLink:

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.

Specify the values of the radii for each layer:

Use functions from OpenCascadeLink to generate 4 concentric balls:

Convert the shapes to faces to define internal boundaries in the geometry:

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

Define an object with 4 faces and 4 regions:

Generate the boundary mesh:

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.

Relate layer names to the markers:

Create a full ElementMesh:

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₀) and the external face (r=r₃) 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₀θ²(π-θ)²ϕ²(π-ϕ)² will be imposed, where θ is the polar angle and ϕ is the azimuthal angle in spherical coordinates.

Define the material parameters of the model by layer:

Define a helper function to specify the material parameters for each layer with the use of ElementMarker:

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

Define the heat transfer model parameters for all the geometry:

Set up the model variables vars and parameters pars of the model:

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.

Define transformation rules from spherical to cartesian coordinates:

Specify q₀, 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.

Set the temperature 0[ºC] at all times at r=r₀:

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.

Create the heat flux on the outside where r=r₃, and for 0<=θ<π and 0<=ϕ<=π:

Set up the initial temperature of the multilayer sphere to T_initial=0[ºC]:

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₀ at y=r₃ 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.

Define the initial and final time of integration:

Define the heat transfer equation to solve with the use of the HeatTransferPDEComponent:

Solve the heat transfer PDE model and monitor the time/memory usage:

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) ϕ=π.

We define the times for which the temperature will be plot:

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.

Define the radial cuts of the multilayer sphere for the plots in the directions of {x,y,-y}:

Generate the time legend values:

Plot the temperature field in the three radial cuts:

Show the solutions:

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.

Create six lists of nine temperature contour values for each time t and relate them with the times list:

Construct the cross section of the multilayer sphere:

Color style of contour lines:

Define a function that uses ContourPlot and plots the temperature field at time t:

Display the plots in an array:

Show all contour plots at different times:

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 κ₃>κ₂>κ₁, layer 3 (outer layer) transports heat and increases its temperature faster than the internal layer as can be seen in the density plot animation.

Define the temperature range values and the legend bar format with BarLegend:

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.

Define the plot range limit and the density plot options:

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.

Show the animation:

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.

Top