Thermal conduction with changing shape and temperatured dependent thermal conductivity

Question

I asked a simpler version of this question earlier and some of the comments are very helpful, especially the youtube tutorials that talk about the case where the cross-section area and the thermal conductivity changes linearly as a function of x. Let's say in a more realistic scenario, I have a metal bar with the same material all the way. The cross-section area changes linearly as a function of x with the relation of A=10+x, while the thermal conductivity depends on the local temperature instead of changing linearly with x, and the relationship is k=10+T, where T is the local temperature.

My question is, for such a complicated case, is it possible to write a differential equation as given in the previous videos and solve it numerically or do we have to rely on some simulation such as finite elements? If a differential equation can be written, could you give me some hints? Thanks a lot.

Answer 1

Well, the total rate of flow of heat Q is constant. So we have $$-k(T)A(x)\frac{dT}{dx}=Q$$This separates into $$-k(T)dT=Q\frac{dx}{A(x)}$$Both sides can be integrated.

Answer 2

3D equation. You can write a PDE for conduction in solids with continuous properties. This equation comes from the energy balance equation and reads

$\dfrac{\partial e}{\partial t} = - \nabla \cdot \mathbf{q}$

being $e$ the energy that (under some assumptions) can be written as $e = \rho c T$ and $\mathbf{q}$ the heat flux vector, that can be written as $\mathbf{q} = - k \nabla T$ if Fourier's law holds. To complete the problem, you need initial and boundary conditions.

1D equation. Let's focus now on the 1D dimensional equation we can derive for elongated solid elements, like beams or rods.

• assuming constant density $\rho$ and heat capacity $c$ (that we thus take outside the time derivative),
• after integration on the section of the solid element

$\rho c A \dfrac{\partial T}{\partial t} = - \dfrac{\partial (A q)}{\partial x} - Q_{lateral}$,

being $Q_{lateral}$ the heat flux exchanged through the lateral surface of the element. Assuming that the element lateral surface is insulated, $Q_{lateral} = 0$. The 1D expression of Fourier's law reads

$q = - k \dfrac{\partial T}{\partial x}$,

and thus the differential equation reads

$\rho c A(x) \dfrac{\partial T}{\partial t}(x,t) = \dfrac{\partial}{\partial x} \left( A(x) k(T) \dfrac{\partial}{\partial x} \left( T(x,t) \right) \right)$$\quad, \quad t \in [0, T_{fin}], \quad x \in [x_1, x_2]$

with:

• the initial conditions $T(x,0) = T_0(x)$, $x \in [x_1, x_2]$
• the boundary conditions in $x_1$, $x_2$ for $t \in [0, T_{fin}]$, as an example:
  • prescribed temperature (Dirichlet or essential boundary conditions): $T(x_1, t) = T_1(t)$
  • prescribed heat flux (Neumann or natural boundary conditions): $Aq(x_2, t) = A(x_2) k(T) \dfrac{\partial T}{\partial x}\bigg|_{x_2} = Aq_2(t)$.

Numerical solution. Despite this equation looks quite complex, it can be solved quite in a easy way using its weak formulation and finite element approach. We can get the weak formulation as

$\displaystyle \int_{\Omega} w(x) \rho c A \dfrac{\partial T}{\partial t} + \int_{\Omega} A k \dfrac{\partial w}{\partial x}\dfrac{\partial T}{\partial x} = w(x_2) Q_2(t)$, $\forall w(x) \in W$.

(This is just a hint about the numerical method, without explicitly referring to mathematical details, especially about the weak formulation and the test space $W$ the test function $w(x)$ belongs to, and the treatment of the essential boundary conditions)