So I have been working on creating a 1-D model of a material with a layer of heat resistant material stuck to it. For this I have been trying to use the implicit finite-differential method. The equations are as follows:

$$ \dot{E}_{in} - \dot{E}_{out} = \dot{E}_{st} $$

So, if there is only one material, an internal node equation would look like (I am assuming all flow in into the node for sign convention, n is the node location, p is the time step, and x is the thickness of a node): $$ \frac{k A}{\Delta x}(T^{p+1}_{n+1}-T^{p+1}_{n}) + \frac{k A}{\Delta x}(T^{p+1}_{n-1}-T^{p+1}_{n}) = \frac{C_p \rho A \Delta x }{\Delta t}(T^{p+1}_{n}-T^{p}_{n}) $$ The surface node equation looks like: $$ {h A}(T_\infty-T^{p+1}_{1}) + \frac{k*A}{\Delta x}(T^{p+1}_{n-1}-T^{p+1}_{n}) = \frac{C_p \rho A \Delta x }{2 \Delta t}(T^{p+1}_{n}-T^{p}_{n}) $$ This seems to work pretty good for me when I use just one material. I give it some initial conditions and run it to steady state then check it against a simple steady state resistance model and there is a match.

I begin to bump into issues when I put in an internal material boundary. This is the equation I was using to simulate it (This calc is for a node in material A adjacent to a node in material B)

$$ \frac{k_b A}{\Delta x_b}(T^{p+1}_{n+1}-T^{p+1}_{n}) + \frac{k_a A}{\Delta x_a}(T^{p+1}_{n-1}-T^{p+1}_{n}) = \frac{C_{pa} \rho_a A \Delta x_a }{\Delta t}(T^{p+1}_{n}-T^{p}_{n}) $$

When I run this, I get a steady-state difference of about 10% from the resistance model. I also created a temperature steady-state modeling by assuming the $\dot{E}_{st}$ term was zero $(T^{p+1}_{1}=T^{p}_{1})$ at SS) then iteratively solved for the temperature at each node. This model was close to the resistive model but still fairly far off the transient model, so that really makes me think there is an issue with my material transition equations.

I also tried using half nodes at the boundary and changing around the node sizes, neither seemed to make a difference.

Any help is appreciated.

  • 2
    $\begingroup$ Also, i just realized, your question is off-topic as it is not about physics, see the Help Center. You are asking why your implementation does not yield the same answer as the steady-state model; this would be more suited for Engineering. $\endgroup$
    – nluigi
    Oct 26, 2015 at 14:33
  • 1
    $\begingroup$ @nluigi: probably more relevant would be Computational Science. As it stands though, I'm not even sure what OP is even asking here. $\endgroup$
    – Kyle Kanos
    Oct 26, 2015 at 14:53
  • $\begingroup$ Really, the question is: What are the energy equations for the boundary layer in a 1-D, 2 material system. All the stuff about modeling is background, perhaps better left out. $\endgroup$
    – chap178
    Oct 26, 2015 at 15:12

1 Answer 1


Your difference equation for the interface should really give an accurate answer for the steady state. It automatically conserves heat.

For the unsteady state, I don't like the right hand side. I would use $$(\rho_a C_{pa}\Delta x_a+\rho_b C_{pb}\Delta x_b)/2$$. However, this doesn't explain the discrepancy for the steady state. The difference between the specially constructed steady state model and the transient model strongly suggests that the transient model had not yet reached steady state. Is this a possibility?

  • $\begingroup$ The change you propose makes sense when the heat capacity is not the same for the two materials - as they might well not be... $\endgroup$
    – Floris
    Mar 29, 2017 at 21:33

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.