Boundary conditions at interface for1D transient heat conduction problem in cylindrical coordinates

Question

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.

Answer 1

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.

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).

Answer 2

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.