I am trying to solve a 1D transient heat conduction problem using the finite volume method (FVM), with a fully implicit scheme, in polar coordinates. The governing equation is written as:

$ \frac{\partial \rho C_p T}{\partial t} = \frac{1}{r}\frac{\partial}{\partial r}\left(kr\frac{\partial T}{\partial r}\right)$.

Essentially, the problem is heat conduction in an infinitely long cylinder. The cylinder is made of two solid concentric cylinders (they have the same thermal properties). The body is heated by convection. The schematic heat flow diagram is shown in the figure below.

enter image description here


Using the standard notation of P, as in the figure below, could you please help me to discretize the equation near the interface and at the boundary?

enter image description here

Could you please kindly confirm that the discretization for the first control volume (CV) and the middle CVs is correct?

1st CV: $(\rho C_p \frac{\Delta r r_P}{\Delta t} + \frac{k_w r_w}{\delta r_w})T_P = \frac{k_w r_w}{\delta r_w}T_W + \rho C_p \frac{\Delta r r_P}{\Delta t} T_P^o$,


Middle CVs: $(\rho C_p \frac{\Delta r r_P}{\Delta t} + \frac{k_e r_e}{\delta r_e} + \frac{k_w r_w}{\delta r_w})T_P = \frac{k_e r_e}{\delta r_e}T_E + \frac{k_w r_w}{\delta r_w}T_W + \rho C_p \frac{\Delta r r_P}{\Delta t} T_P^o$,

where $r_P$ is the distance from the center of the cylinder to the node P, $r_e$ is the distance from the center of the cylinder to the $e$ face of the CF, $r_w$ is the distance from the center of the cylinder to the $w$ face of the CF, $\Delta r$ is the width of the CV, $\delta r_e$ is the distance from the node P to the node E, and $\delta r_w$ is the distance from the node P to the node W.

  • $\begingroup$ Make sure to apply symmetry conditions to simplify your problem, i.e. $\partial_rT\left(0\right)=0$ where $r=0$ is the center of the cylinder. $\endgroup$
    – nluigi
    Commented Jul 18, 2017 at 8:05

2 Answers 2


I am assuming you want heat flux continuity at the interface:

$$T_i^1=T_i^2 \quad k_1\left.\frac{\partial T}{\partial r}\right|_i^1=k_2\left.\frac{\partial T}{\partial r}\right|_i^2$$

where $i$ indicates the interface and $1,2$ is the material respectively left, right of the interaces.

Furthermore, you want a heat transfer model at the boundary:

$$-k_2\left.\frac{\partial T}{\partial r}\right|_b = h\left(T_b-T_a\right)$$

where $T_b$ is the temperature at the boundary and $h$ is some heat transfer coefficient which depends on the fluid flow cooling the boundary.

Furthermore, since you are using FVM you are using a staggered grid (i.e. domain boundaries are located in between grid nodes) and I will assume the use of ghost nodes to manipulate the values at the interface and boundary.

Discretization of the temperature and its derivative at the boundary is done by:

$$T_i^\sigma = \frac{T_{i-}^\sigma+T_{i+}^\sigma}{2} \quad \left.\frac{\partial T}{\partial r}\right|_i^\sigma = \frac{T_{i+}^\sigma-T_{i-}^\sigma}{\Delta r}$$

where $-/+$ indicates the (ghost) nodes just respectively left/right from the interface. Applying it to the above conditions at the interface we arrive at the following system of equations:

$$T_{i+}^1 - T_{i-}^2 = -T_{i-}^1 + T_{i+}^2 $$ $$k_1T_{i+}^1 + k_2T_{i-}^2 = k_1T_{i-}^1 + k_2T_{i+}^2 $$

where $T_{i+}^1$ and $T_{i-}^2$ are unknown temperatures at the ghost nodes in material 1 to the right of the interface and material 2 to the left of the interface respectively. Solving for these temperatures yields:

$$T_{i+}^1 = \frac{k_1-k_2}{k_1+k_2}T_{i-}^1 + \frac{2k_2}{k_1+k_2}T_{i+}^2$$ $$T_{i-}^2 = \frac{2k_1}{k_1+k_2}T_{i-}^1 - \frac{k_1-k_2}{k_1+k_2}T_{i+}^2$$

Note that if $k_1=k_2$ (i.e. the two materials are the same) then $T_{i+}^1=T_{i+}^2$ and $T_{i-}^2=T_{i-}^1$ which indicates the material interface has disappeared.

The same method can be applied to get a value at the ghost node near the boundary: $$T_{b+}^2 = \frac{1-\mathrm{Nu}/2}{1+\mathrm{Nu/2}}T_{b-}^2 + \frac{\mathrm{Nu}}{1+\mathrm{Nu}/2}T_a$$

where $\mathrm{Nu}=\frac{h\Delta r}{k_2}$ is the Nusselt number which is the relative strength of convection to conduction.

How to get to this equation will be left as an exercise for the reader.

  • $\begingroup$ Hi, nluigi! Thanks for posting your method - the Kirchhoff Approximation (if I am not mistaken) or similar to it - of discretizing temperature at an interface between two materials. I saw it was described by Voller V. R. and Swaminathan C. R. (1993) in "Treatment of discontinuous thermal conductivity in control-volume solutions of phase-change problems", Numerical Heat Transfer, Part B Fundamentals, 24(2), 161-180. Your explanation, however, is more elegant and clear. Thanks a lot again! $\endgroup$ Commented Jul 20, 2017 at 4:24

I found the problem in my discretization.

When I integrated the PDE over a control volume and time step, I obtained the following equation:

$\rho C_p (T_P -T_P^o) \Delta V = \left(k_e A_e \frac{T_E - T_P}{\delta r_e} - k_w A_w \frac{T_P - T_W}{\delta r_w} \right)\Delta t$.

My mistake was to equate $\Delta V$ and $A$ to $\Delta r$ and $1$,respectively; as that was 1D in the Cartesian coordinates. However, in 1D of the cylindrical coordinates, $\Delta V = \pi (r_e^2 -r_w^2)$ and $A = 2 \pi r$.

Hence, the general solution in 1D cylinrical coordinates can be written as:

$\rho C_P (T_P - T_P^o)(r_e^2 - r_w^2) = \left(2 k_e r_e \frac{T_E - T_P}{\delta r_e} - 2 k_e r_e \frac{T_E - T_P}{\delta r_e} \right) \Delta t$.

For the 1st CV: $\left(\rho C_P \frac{r_e^2}{\Delta t} + \frac{2 k_e r_e}{\delta r_e} \right) T_P = \frac{2 k_e r_e}{\delta r_e} T_E + \rho C_P \frac{r_e^2}{\Delta t}T_P^o$

For the middle CVs: $\left(\rho C_P \frac{r_e^2 - r_w^2}{\Delta t} + \frac{2 k_e r_e}{\delta r_e} + \frac{2 k_w r_w}{\delta r_w} \right) T_P = \frac{2 k_e r_e}{\delta r_e} T_E + \frac{2 k_w r_w}{\delta r_w} T_W + \rho C_P \frac{r_e^2 - r_w^2}{\Delta t}T_P^o$

For the last CV at the boundary: $\left(\rho C_P \frac{r^2 - r_w^2}{\Delta t} + \frac{2 k_w r_w}{\delta r_w} - ( -2 rh) \right) T_P = \frac{2 k_w r_w}{\delta r_w} T_W + \rho C_P \frac{r^2 - r_w^2}{\Delta t}T_P^o + 2rhT_a$

To discretize the equation at an interface use a harmonic mean, thouroughly described by Patankar (1980).

Patankar, Suhas V., Numerical heat transfer and fluid flow, Series in Computational Methods in Mechanics and Thermal Sciences. Washington - New York - London: Hemisphere Publishing Corporation; New York etc.: McGraw-Hill Book Company. XIII, 197 p. \$ 37.50 (1980). ZBL0521.76003.


Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.