Steady-State Heat Conduction in Spherical Shells
Theory and Physics
Difference Between Spherical Shell and Cylinder
Professor, what's the difference between steady-state heat conduction in a spherical shell and in cylindrical coordinates? They both have a "round shape," so they seem similar.
The biggest difference is how the heat transfer area increases. For a cylinder, the area is proportional to the radius ($A = 2\pi r L$), but for a sphere, the area is proportional to the square of the radius ($A = 4\pi r^2$).
If the way the area increases is different, does the temperature distribution also change?
Exactly. For a cylinder, the temperature varies with $\ln r$ (logarithmic function), whereas for a spherical shell, it varies with $1/r$. Intuitively, because the area increases more rapidly towards the outside in a spherical shell, the temperature drop decays more quickly.
| Shape | Area $A(r)$ | Temperature Distribution | Form of Thermal Resistance |
|---|---|---|---|
| Flat Plate | Constant | $T \propto r$ (Linear) | $L/(kA)$ |
| Cylinder | $\propto r$ | $T \propto \ln r$ | $\ln(r_2/r_1)/(2\pi k L)$ |
| Spherical Shell | $\propto r^2$ | $T \propto 1/r$ | $(1/r_1 - 1/r_2)/(4\pi k)$ |
I see, so as the "exponent" of the area increases, the temperature drop becomes more gradual. In practical applications, where are spherical shells used?
There are three typical examples.
- LNG Spherical Tank (approx. 40m diameter) insulation design — Prediction of boil-off gas volume
- Nuclear Reactor Pressure Vessel temperature distribution through the wall thickness — Input for thermal stress evaluation
- Spherical Fuel Pellet (HTGR: TRISO fuel for High-Temperature Gas-cooled Reactors) central temperature estimation
In all cases, spherical symmetry can be assumed, so a 1D analytical solution can provide a good first approximation for design. That's the strength.
Governing Equation and Analytical Solution
So, please tell me the governing equation for a spherical shell. For a cylinder, there was a $1/r$ term, right?
The steady-state, 1D heat conduction equation in spherical coordinates, considering only the radial direction, is as follows.
Compared to the cylinder's $\frac{1}{r}\frac{d}{dr}\left(kr\frac{dT}{dr}\right)=0$, the exponent of $r$ changes from 1 to 2. This reflects the fact that the area is proportional to $r^2$.
If $k$ is constant, how do you integrate it?
When $k$ is constant, integrating $\frac{d}{dr}\left(r^2 \frac{dT}{dr}\right) = 0$ twice gives
Substituting the boundary conditions $T(r_1) = T_1$, $T(r_2) = T_2$ yields the analytical solution for the temperature distribution.
By analogy with the $\ln r$ distribution for a cylinder, it takes the form where $\ln r \to 1/r$.
What shape of graph does a $1/r$ distribution actually produce?
For example, for a steel spherical shell with inner diameter 100mm, outer diameter 200mm, $T_1 = 500$°C, $T_2 = 100$°C:
- $r = 100$ mm: $T = 500$°C (inner surface)
- $r = 133$ mm (one-third of the wall thickness): $T = 300$°C
- $r = 200$ mm: $T = 100$°C (outer surface)
For a flat plate, the temperature at one-third of the thickness would be $(500-100)\times(1-1/3)+100 = 367$°C, so the temperature drop is steeper near the inner surface in a spherical shell. This also means the heat flux near the inner surface is higher.
Thermal Resistance of a Spherical Shell
How do you calculate the heat flow rate through a spherical shell? Can you use the concept of thermal resistance like for a cylinder?
Yes, you can. Substituting $A(r) = 4\pi r^2$ into Fourier's law $q = -kA(r)\frac{dT}{dr}$ and integrating gives
Arranging this into the form $Q = \Delta T / R$, the conductive thermal resistance of a spherical shell is
A characteristic of spherical shells is that the length $L$ does not appear, unlike in cylinders. Because a sphere is a closed shape in all directions, there is no parameter equivalent to the length of a pipe.
If there is convection on the inner and outer surfaces, do we just add the convective thermal resistances?
Yes. The convective thermal resistance for a spherical surface is $R_{\text{conv}} = 1/(4\pi r^2 h)$, so from the inner fluid temperature $T_{\infty,1}$ to the outer fluid temperature $T_{\infty,2}$:
It's a series circuit of convection + conduction + convection. It has exactly the same structure as Ohm's law in electrical circuits.
Multi-Layer Spherical Shell Model
For cases like LNG tanks where insulation material is sandwiched, you calculate it as multiple layers, right?
For an $n$-layer spherical shell where each layer has thermal conductivity $k_i$ and radii $r_0, r_1, \ldots, r_n$, the total conductive thermal resistance is
The overall heat transfer coefficient $U$, including convection, with the reference area taken as the inner surface $A_1 = 4\pi r_0^2$, is
Using this, you can parametrically calculate heat loss when changing the insulation thickness, right?
That's right. Taking an LNG spherical tank (40m diameter) as an example, with a steel shell 30mm ($k = 50$ W/(m·K)) + perlite insulation 500mm ($k = 0.04$ W/(m·K)), the thermal resistance of the insulation accounts for about 2,500 times that of the steel shell. The strength of the thermal resistance method is that it allows quantitative understanding of which layer is dominant.
Case with Internal Heat Generation
What about cases with internal heat generation, like nuclear reactor fuel spheres?
For a solid sphere (radius $R$, surface temperature $T_s$) with uniform volumetric heat generation $\dot{q}$ [W/m³], the governing equation is
Solving with the symmetry condition at the center $dT/dr|_{r=0} = 0$ and $T(R) = T_s$ gives
The center temperature is $T_{\max} = T_s + \dot{q} R^2 / (6k)$. For a cylinder, the denominator was $4k$, so the temperature rise at the center is smaller for a sphere. This is because the surface area is larger, making heat dissipation easier.
So for a High-Temperature Gas-cooled Reactor's TRISO fuel particle (approx. 1mm diameter), you can estimate the center temperature with this formula.
Yes. For a UO₂ kernel ($k \approx 3$ W/(m·K), $\dot{q} \approx 2 \times 10^8$ W/m³, $R = 0.25$ mm), $\Delta T_{\max} = 2\times10^8 \times (0.25\times10^{-3})^2 / (6 \times 3) \approx 0.7$°C. Because the particle is small, the temperature rise is also slight. However, in reality, the contact thermal resistance of the multilayer coating (SiC layer, etc.) needs to be considered.
Spherical Shell Heat Conduction and Earth Science
The theory of steady-state heat conduction in spherical shells dates back to Poisson (1820s) and Lord Kelvin in the 19th century. Kelvin modeled the Earth as a cooling sphere and attempted to estimate the Earth's age from the $1/r$ distribution of steady-state heat conduction (he estimated about 100 million years, but this was a significant underestimate because he was unaware of internal heat generation from radioactive decay). Even today, spherical shell heat conduction is used as an initial approximation in models of Earth's mantle convection.
Derivation of the Spherical Coordinate Laplacian
- Laplacian $\nabla^2 T$ (spherical coordinates): Transforming the Cartesian $\frac{\partial^2 T}{\partial x^2}+\frac{\partial^2 T}{\partial y^2}+\frac{\partial^2 T}{\partial z^2}$ to spherical coordinates $(r,\theta,\phi)$ gives $$\nabla^2 T = \frac{1}{r^2}\frac{\partial}{\partial r}\left(r^2\frac{\partial T}{\partial r}\right) + \frac{1}{r^2\sin\theta}\frac{\partial}{\partial\theta}\left(\sin\theta\frac{\partial T}{\partial\theta}\right) + \frac{1}{r^2\sin^2\theta}\frac{\partial^2 T}{\partial\phi^2}$$ For spherical symmetry ($T$ is a function of $r$ only), the $\theta$ and $\phi$ terms vanish, leaving $\frac{1}{r^2}\frac{d}{dr}\left(r^2\frac{dT}{dr}\right) = 0$.
- Physical meaning of $r^2$: Originates from the area of a sphere at radius $r$, $4\pi r^2$. For a constant heat flow rate $Q$, the heat flux $q'' = Q/(4\pi r^2)$ decreases inversely with $r^2$.
Assumptions and Applicability Limits
- Steady State: No time variation ($\partial T/\partial t = 0$). Temperature field after sufficient time has elapsed for transient effects to be negligible.
- Spherical Symmetry: Temperature is a function of $r$ only (independent of $\theta$, $\phi$). Not applicable if there are local heat sources or non-uniform convection.
- Isotropic & Homogeneous: Thermal conductivity $k$ is constant within each layer. Nonlinear analysis is required if temperature dependence is significant (e.g., steel's $k$ changes from 50→25 W/(m·K) from 0–800°C).
- Fourier's Law: $q = -k \nabla T$. Non-Fourier models are needed for extremely low temperatures (superfluid helium) or femtosecond laser heating.
Dimensional Analysis and Key Parameters
| Variable | SI Unit | Typical Values / Notes |
|---|---|---|
| Temperature $T$ | K or °C | Use absolute temperature [K] when radiation is involved |
| Thermal Conductivity $k$ | W/(m·K) | Steel≈50, Insulation≈0.03–0.05, UO₂≈3 |
| Inner Radius $r_1$ | m | For spherical shell $r_1 > 0$ (solid sphere is a different solution) |
| Spherical Shell Thermal Resistance $R$ | K/W | $(1/r_1 - 1/r_2)/(4\pi k)$ — Does not contain length $L$ |
| Volumetric Heat Generation $\dot{q}$ | W/m³ | Nuclear fuel≈$10^8$, Joule heating≈$10^5$–$10^7$ |
Numerical Methods and Implementation
Axisymmetric FEM Model
When solving spherical shell heat conduction with FEM, do you have to create a full 3D spherical model?
For spherical symmetry, a 2D axisymmetric model is sufficient. You draw a semicircle cross-section, which corresponds to a model rotated around the symmetry axis. The computational cost becomes less than 1/1000th of a 3D model.
| Model Type | Approx. Number of Elements | Degrees of Freedom | Computation Time |
|---|---|---|---|
| 3D Full (Entire Sphere) | 500k– | 500k– | Minutes |
| 3D 1/8 Symmetry | 60k– | 60k– | Seconds–Minutes |
| 2D Axisymmetric | 200– | 200– | Instant |
| 1D Analytical Solution | None | None | Hand Calculation |
What should I be careful about with an axisymmetric model?
There are three main points.
- Element Type Selection: In Ansys Mechanical, use PLANE55 (linear) / PLANE77 (quadratic) with KEYOPT(3)=1 (axisymmetric). In Abaqus, use DCAX4 / DCAX8.
- Boundary Conditions on the Symmetry Axis: An adiabatic condition ($\partial T/\partial r = 0$) is automatically applied on the $r=0$ axis, but some software may require explicit setting.
- Mesh Bias: The temperature gradient is steep near the inner surface, so use a finer mesh there and a coarser mesh towards the outer surface for efficiency. A biased mesh with inner element size 1/3 to 1/5 of the outer size is recommended.
Discretization with Finite Difference Method
If I want to implement it myself with the Finite Difference Method instead of FEM, how should I do it? I'd like to try writing it in Python.
Expanding the steady-state heat conduction equation in spherical coordinates gives
Discretizing this with central differences, for node $i$ (position $r_i$):
Rearranging gives a tridiagonal matrix system of equations.
The influence of $\Delta r / r_i$ becomes larger where $r_i$ is small (near the inner surface), so with a uniform mesh, accuracy degrades near the inner surface. Using a non-uniform mesh or the variable transformation $s = 1/r$ improves accuracy.
The transformation $s = 1/r$ is interesting. What happens after the transformation?
Setting $s = 1/r$, the spherical shell equation becomes $\frac{d^2T}{ds^2} = 0$. That means in $s$-space, the temperature distribution becomes linear. Discretizing with a uniform $s$ grid gives an exact solution with central differences. This transformation is a technique specific to spherical shell heat conduction.
Mesh Design Guidelines
When meshing for FEM, what criteria should I use to decide the element size?
For steady-state heat conduction in a spherical shell, the following guidelines are good.
| Item | Recommended Value | Reason |
|---|---|---|
| Number of layers in thickness direction | 10 or more | To accurately reproduce the $1/r$ distribution |
| Bias ratio on inner side | 1:3 to 1:5 | Temperature gradient is steep at inner surface |
| Element type | Quadratic elements recommended | Linear elements provide a poor approximation for $1/r$ distribution |
| Multi-layer interface | Make nodes coincident | Because temperature gradient is discontinuous at the interface |
| Aspect ratio | 5 or less | Excessively flat elements reduce accuracy |
Related Topics
なった
詳しく
報告