The problem I am currently working on involves rapidly cooling liquid aluminum in a cylindrical graphite crucible and getting the radial temperature distribution over time, possibly using finite difference methods.

This is a 1D transient problem in only the radial direction using the following equation and I am only concerned about the liquid portion of aluminum (no source term and constant thermal properties).

$$\frac{1}{a}\frac{dT}{dt} = \frac{1}{r}\frac{d}{dr}\left(r\frac{dT}{dr}\right)$$

The main problem I am having is determining the boundary conditions, especially at the interface between the graphite inner wall with the liquid aluminum (as graphite and liquid aluminum have different thermal diffusivities). The boundary condition at the outer graphite wall involves convection and the BC at the centerline of the cylinder should have no heat flux for symmetry but I am not sure how to get a boundary condition at the interface (I am thinking of equating heat flux of graphite to aluminum).

I was wondering if anyone here has any idea on how to get the boundary condition at the interface or have any credible sources that I can look at. Thank you.

  • 1
    $\begingroup$ "(I am thinking of equating heat flux of graphite to aluminum)." That is correct. It will give you an equation linking $dT/dr$ on either side of the boundary. At the boundary , $T$ is a continuous function of $r$ but $dT/dr$ is discontinuous. $\endgroup$
    – alephzero
    May 22, 2021 at 3:11
  • $\begingroup$ Thank you for answering so quickly. However, I am not completely sure about how to set up this boundary condition in cylindrical coordinates (I'm pretty new to this topic). Could you possibly show me the starting equation at the interface or if you have any recommended sources for this topic? $\endgroup$ May 22, 2021 at 4:06
  • $\begingroup$ if the $\frac{1}{a}$ coefficient in the analysis is a measure related to conductivity or resistance, then you should be able to just adjust the parameter and use different for aluminium and the graphite. Additionally you could also use a different parameter for the aluminium and changing temperature. $\endgroup$
    – NMech
    May 22, 2021 at 9:20

2 Answers 2


Interface condition

From what I gather the $\frac{1}{a}$ coefficient in the equation is a measure related to conductivity or resistance. Therefore all you need to is that you should be able to just adjust the parameter $a$ and use different for aluminum and the graphite depending on the position of the analysis. So the equation will have the following form.

$$\frac{1}{a (r) }\frac{dT}{dt} = \frac{1}{r}\frac{d}{dr}(r\frac{dT}{dr})$$

Additionally you could also use a different parameter for the aluminum and changing temperature (because I expect the aluminum will have a greater temperature difference. In that case you would have a function of the radius and temperature

$$\frac{1}{a (r,T) }\frac{dT}{dt} = \frac{1}{r}\frac{d}{dr}(r\frac{dT}{dr})$$

So the basic formulation does not change. The only thing that changes is how you calculate the $a(r)$.

symmetry boundary condition (center of the crucible)

The gradient needs to be zero at the symmetry boundary, because at the that point you shouldn't have any heat flux across the symmetry boundary. At that point you could use a von Neumann boundary condition.

Basically is a condition where you define the gradient of the property (in this case temperature).

If you have N nodes (starting from 1 and ending at N), where N is the center node, what you need to do (Alephzero can correct me, because I am a bit rusty), is you define another node, and on each update you set the value of the N+1 node equal to N. That will ensure that the property gradient is zero.

the N+1 code is sometimes called a Ghost node. The following image is an example for Ghost node that can be used for a von Neuman condition at node 1.

enter image description here

Figure: Ghost nodes source (Computational ScienceSE)

good resource

Last time I did something similar, I used the following book by Charpa and Canale as a example (actually its my go to book for numerical methods for at least 20 years).

  • $\begingroup$ @alephzero if you have some time, could you please assess the validity of my answer? $\endgroup$
    – NMech
    May 22, 2021 at 9:38
  • $\begingroup$ Thank you for the answer and the resource. I understand what you mean by having the thermal diffusivity as a function of radius but I was just wondering if I could do this problem by separating this into 2 regions. One region will involve only the graphite wall with the convective BC on the outside and the interface BC on the inside, and the other region with only liquid aluminum, with the interface BC and the von Neumann BC at the center, essentially having 2 regions for finite differencing with the equations being connected by the interface boundary condition. $\endgroup$ May 22, 2021 at 14:33
  • $\begingroup$ I am not certain what you mean by separating into 2 regions. You mean to solve two separate problems? Even if you did that (in order to save yourself from the added headache of changing $a(r)$ and just have a constant $a$ in each region), you would still need the interaction between the two regions. I mean, even if the temperature at the external side of the graphene is constant, the T at the interface would still be depended (to some extent) by the temperature of the liquid aluminum. $\endgroup$
    – NMech
    May 22, 2021 at 14:39
  • $\begingroup$ Yes, I meant I wanted to solve this as 2 separate problems. I was hoping to have a heat flux BC at the interface for the graphite section (which should give the temperature at that interface) and using that temperature as the BC for the aluminum section but I'm not sure if this is possible. Does it help that the initial condition for both the graphite and aluminum is the same for this problem? $\endgroup$ May 22, 2021 at 14:48
  • 1
    $\begingroup$ @NMech Your answer looks OK to me, though it's a few decades since I did any serious heat transfer calculations using finite difference methods. Personally I would use finite elements instead. where you don't need things like "ghost points." $\endgroup$
    – alephzero
    May 23, 2021 at 17:06

The cylindrical coordinates don't make any difference.

The physical boundary condition is that the radial heat flux on each side of the boundary is the same. That means $$k_1 \left.\frac{dT}{dr}\right|_1 = k_2 \left.\frac{dT}{dr}\right|_2$$ where $k_1$ and $k_2$ are the thermal conductivities of the materials on each side of the boundary.

You can estimate the gradients using your finite difference mesh. The simplest formula uses "one-sided" diferences at grid point $p$ to give $$\left.\frac{dT}{dr}\right|_p \approx \frac{T_{p+1} - T_p}{h} \approx \frac{T_p - T_{{p-1}}}{h}$$ so (for equally spaced mesh points) the boundary condition at the interface is $$k_1(T_p - T_{p-1}) = k_2(T_{p+1} - T_p).$$

Or you can use a second order one-sided estimate of the gradient $$\left.\frac{dT}{dr}\right|_p \approx \frac{-T_{p+2} + 4T_{p+1} - 3T_p}{2h} \approx \frac{3T_{p} - 4T_{p-1} + T_{p-2}}{2h}$$ so the boundary condition is $$k_1(3T_p + -4T_{p-1} + T_{p-2}) = k_2(-T_{p+2} + 4T_{p+1} - 3T_p).$$

Incidentally, you can also use the one-sided estimate of the gradient to impose the boundary condition $$\left.\frac{dT}{dr}\right|_{r=0} = 0$$ at the origin without using an extra "ghost point" in the mesh.

  • $\begingroup$ Thank you. I will also try this. $\endgroup$ May 24, 2021 at 15:34

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.