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.


2 Answers 2


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.

  • $\begingroup$ Thanks! Could I confirm that the integral limit of T is just the temperatures at the two ends of the metal bar? For example, if one side is 0C and the other side is 100C, the limit would just be from 0 to 100. If that's the case, intuitively, it seems like Q doesn't care about the temperature profile at all. Is there any physical interpretation regarding this? $\endgroup$
    – usr01
    Oct 28, 2022 at 15:12
  • $\begingroup$ "it seems like Q doesn't care about the temperature profile at all" As shown, the integral is affected by $k(T)$ at every individual point. This results in $Q$ depending on all the intermediate temperatures, not just the boundary conditions. $\endgroup$ Oct 28, 2022 at 16:39
  • $\begingroup$ The temperature gradient automatically adjusts to the two end temperatures, the function k(T), and the function A(x), and so Q is intimately related to the temperature profile. You tell me Q and the temperature at one of the ends, and I'll tell you the entire temperature profile. $\endgroup$ Oct 28, 2022 at 20:31

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]$


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

  • $\begingroup$ Thanks a lot. I assume when the system reaches equilibrium state, it can be further simplified to $\dfrac{\partial}{\partial x} \left( A(x) \dfrac{\partial T}{\partial x}k(T) \right)$$= 0$. And we would only need the boundary condition about the temperature and T[x=x1] and T[x=x2] to solve this equation, right? $\endgroup$
    – usr01
    Oct 28, 2022 at 12:54
  • 2
    $\begingroup$ No, you are wrong. The equation you wrote in your previous comment is incorrect. It should read $$\frac{\partial}{\partial x}\left(A(x)k(T)\frac{\partial T}{\partial x}\right)=0$$ $\endgroup$ Oct 28, 2022 at 13:37
  • 1
    $\begingroup$ @ChetMiller is right, I messed up with $K(T)$: it's only under the outermost derivative. I'm editing my answer. Now you can see that his equation is the steady version of the equation above, if you put $\frac{\partial T}{\partial t} = 0$. $\endgroup$
    – basics
    Oct 28, 2022 at 13:53

Your Answer

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

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