1. Limited time only! Sign up for a free 30min personal tutor trial with Chegg Tutors
    Dismiss Notice
Dismiss Notice
Join Physics Forums Today!
The friendliest, high quality science and math community on the planet! Everyone who loves science is here!

Heat equation with const. heat source

  1. Oct 29, 2008 #1

    I have the following problem and would appreciate any help you might be able to give.

    In reality, I have a 2D/3D problem (self heating of a semiconductor device), but I'm absolutely happy if I would understand the 1D simplification:

    I have a rod of length L (semi infinite would be even better). One end is at a fixed temperature (heat sink) and the other end is attached to a electric heater with constant power. The entire rod is at the temperature of the heat sink at time t=0. Then the heater is switched on. I'm interested in the temperature T as function of time at the place of the heater.

    I ran thermal simulations of the problem which suggest that T is proportional to t0.58 (1D); t0.38 (real 2D) for the initial heating. But why?

    I guess there is an analytical solution for it where you can see the relation ship, but I wasn't able to find any.

    Thanks for your help
  2. jcsd
  3. Oct 29, 2008 #2
    I don't know how you did the simulations but something is missing: how much heat is transferred from the heater to the rod. The temperature of the heater is not enough. You can have a very hot heater but giving very little heat to the rod.
  4. Oct 29, 2008 #3
    Dear Nasu

    The heater represents Joule heating of a electronic device, dissipating constant power (e.g. 1W). In my simulation tool I can just put this as a heat load (neglecting any volume of the heater).
  5. Oct 29, 2008 #4


    User Avatar
    Staff Emeritus
    Science Advisor
    Homework Helper

    Something is weird about the suggested power dependences t0.58 and t0.38

    Once steady state is reached, the 1-D temperature profile should be a straight line with a constant slope, such that one end is at the heatsink temperature, and the slope of the line obeys the equation:

    \frac{dT}{dx} = \frac{P}{k A}

    where P is the heater power, k is the rod's thermal conductivity, and A is its cross-sectional area.
  6. Oct 30, 2008 #5
    I got the power dependence when I plotted the 1D and 3D simulation (Temperature vs. time) in a log-log-plot. I get a nearly perfect straight line. It might be just a coincidence, or it might be an approximation of the analytical solution.
  7. Oct 30, 2008 #6


    User Avatar
    Staff Emeritus
    Science Advisor
    Homework Helper

    Hmm, that's suprising. Perhaps you could let the simulation run longer, and see if the 1D case approaches the solution I mentioned in post #4 for long times. According to that, the temperature at the heat source should approach, but never exceed:

    Tmax=Tsink+(P L) / (k A)

    However, if it follows the t0.58 dependence, eventually it would exceed that temperature. So this would be a good test of the t0.58 curve vs. what I gave in post #4.

    Are you looking at the temperature at a specific location, or some kind of average over the rod's length?
  8. Oct 30, 2008 #7
    Dear Redbelly

    I'm interested in the temperature evolution right at the location of the heater.

    When I let the simulation run for longer it converges as expected to your Tmax
    Have a look at the graph.
    http://picasaweb.google.com/franck.hertz/Forum#5262954245291443458 [Broken]
    To me, the beginning and the middle part look pretty straight (in log-log) plot, i.e. it follows a power dependence, but why? I guess it must be an approximation of the analytical solution (which I unfortunately don't have).
    Last edited by a moderator: May 3, 2017
  9. Oct 30, 2008 #8


    User Avatar
    Staff Emeritus
    Science Advisor
    Gold Member

    There is an analytic solution to your problem in Operational Mathematics by Churchill. page 145

    [tex] U_t(x,t) = kU_{xx}(x,t) [/tex]
    Initial condition:
    U(x, +0)=0
    with U(x,t) =0 as x -> [itex] \infty [/itex]

    [tex] K U_x(0,t) = \phi_0 [/tex]

    He has for a solution:

    [tex] U(0,t) = \frac {2 \phi_0 } K \sqrt {\frac {kt} \pi [/tex]
  10. Oct 30, 2008 #9


    User Avatar
    Staff Emeritus
    Science Advisor
    Homework Helper

    Thanks Integral. The t1/2 behavior is evident in the first 2000 ns of Franck's graph.

    For a limited time, there would be little difference between the finite and semi-infinite cases. But at some point, the temperature at x=L in a semi-infinite bar would begin to rise appreciably. At that time and later, the two solutions would differ.

    I've found an old Heat Transfer book of mine and will see what it has in the way of conduction along finite-length bars. I'll repost if I find anything useful.
  11. Oct 31, 2008 #10
    Dear Integral

    very interesting and thanks for the reference. I'll checked with our library, we have the book. So I'll go and get it.......

    ... hmmm, someone else borrowed it, so I have to wait for a couple of days.....
    Last edited: Oct 31, 2008
  12. Oct 31, 2008 #11


    User Avatar
    Staff Emeritus
    Science Advisor
    Homework Helper

    FYI, no luck in finding this specific problem in my book (Schaum's Outline Series: Heat Transfer by D.R. Pitts & L.E. Sissom). But the solutions to similar problems are quite complicated, involving the error function erf(z) and exponential functions of x, t, and √t, where
    z = x / √(4 α t)
    α = k / (ρ cp)
  13. Oct 31, 2008 #12


    User Avatar
    Staff Emeritus
    Science Advisor
    Homework Helper

    Okay, I think I have solved it.

    Here is a plot

    I have pretty much set all parameters=1. Franck, if you can tell me the following I could try to reproduce your curve quantitatively as well as qualitatively:

    Length of bar
    Cross-sectional area of bar
    Thermal conductivity
    Heat capacity
    The power of the heat source (1W?)

    My solution is an infinite series of decaying exponentials (in time) and cosine terms (in x). I'll post it shortly, but wanted to get this graph posted 1st.

    I double-checked and for short times this agrees with Integral's solution for the semi-infinite bar. With my parameters this is
    T = 2 (t/π)1/2
    and is indicated by the orange line in the graph.
    Last edited: Oct 31, 2008
  14. Oct 31, 2008 #13


    User Avatar
    Staff Emeritus
    Science Advisor
    Homework Helper

    Here is my solution for a finite length bar, 1D conduction.

    "T" is the temperature rise above ambient.

    \ T(x,t) = m x + b - C \sum e^{-\lambda_n^2 \ \alpha t} \ cos(\lambda_n x) \ / \ n^2 \


    n = 1, 3, 5, ... is the summation index
    λn = π n / (2L)
    mx+b = (L-x) P / (kA) [ This is T(x) as t→∞ ]
    P is the heater power
    L is the rod length
    A is the rod's cross-sectional area
    k is the thermal conductivity
    C = PL / (kA) * (8/π2)
    α = k / (ρ cp)
    ρ is the density
    cp is the heat capacity

    The solution seems to check out. Putting in "1" for all the inputs (L, A, k, etc.) I plotted the results in Excel and got ageement with the conditions:

    T=0 at t=0 for all x
    T=0 at x=L for all t
    dT/dx = -1 at x=0 for t>0
    T → (1-x) as t→∞
    Last edited: Oct 31, 2008
  15. Nov 1, 2008 #14
    Dear Redbelly

    Unfortunately I need to wait till Monday before I can tell you the exact thermal constant I have used. your solution definitely looks looks great. Thanks very much.
  16. Nov 2, 2008 #15


    User Avatar
    Staff Emeritus
    Science Advisor
    Homework Helper

    Okay. I'm working on this at home, and in the morning will leave for work at 9 am Eastern USA time. If I can, I'll try to plug in your numbers and post before going to work. Otherwise, it will have to wait until evening, in which case I think you would see my response Tuesday morning your time (Europe?)

    3D case

    I have calculated the steady state (i.e., t→∞) temperature rise for the 3D case. What I got was:

    T(r=0,t \rightarrow \infty) \equiv T_{max} = \frac{P}{8 \pi k} \left( \frac{3}{r_s} - \frac{1}{R} \right)

    P = input power
    k = thermal conductivity
    rs = radius of small spherical volume, where the heat is generated
    R = radius of object (temperature rise is zero at r=R)

    Assuming rs << R, we can neglect the 1/R term. This means the final temperature is largely determined by the size of the volume in which heat is being generated, and is nearly independent of how large the whole object is.

    I don't know how long it would take to get close to Tmax, but dimensional analysis suggests time scales on the order of

    t1 = (c ρ / k) rs2
    t2 = (c ρ / k) R2 ,
    which are vastly different from each other.

    (It would be nice if this could be checked somehow with a simulation.)
  17. Nov 3, 2008 #16
    Code (Text):
    Dear Redbelly

    Yes, you are right, I'm in Europe (UK).

    I checked the input parameters of my 1D simulation. The thermal parameters are of GaN:

    k=1.9 W/cmK (thermal conductivity)
    [tex]\rho[/tex]=6.1 g/cm3 (density)
    c=1.49 J/gK (heat capacity; not of GaN. GaN would be 0.49, it was a typo)
    P=0.03W power
    A= 100um2 (Cross section)
    L=10um (length).

    Regarding a 3D simulation: I need to check if I can easily construct a sphere. All my thermal models are blocks so for, but I'll try so that you can compare it with your calculation.
  18. Nov 3, 2008 #17
    I managed to get a 3D sphere working in the simulation, but I don't know how accurate it is a it is a simple model (i.e. low mesh density)

    I assumed a point heater at the center. This would lead to an infinite temperature, which is not sensible and the simulation tool does some sort of averaging. Therefore extracted rs by fitting Redbelly's Tmax. =>
    and thermal parameters as above

    With this I get

    t1 = (c ρ / k) rs2=7ns
    t2 = (c ρ / k) R2 =4.8us

    It seem that t1 is too short and t2 too long.

    http://picasaweb.google.com/franck.hertz/Forum#5264420670106402930 [Broken]
    Last edited by a moderator: May 3, 2017
  19. Nov 3, 2008 #18


    User Avatar
    Staff Emeritus
    Science Advisor
    Homework Helper

    Running out of time this morning to post the graph, but using your 1D numbers:

    I get 15.8 C for the max temperature rise as t→∞
    At 4.78 us, the temp rise is 14.7 C
  20. Nov 3, 2008 #19


    User Avatar
    Staff Emeritus
    Science Advisor
    Homework Helper

    1D simulation

    Using those parameters, (c=1.49 J/g-K, not 0.49) I get:

    Time units are seconds.
    My temperature is about half of what you got before, not sure why that is so. The time scales match up pretty well:

    3D steady-state

    I'm posting the complete equation for the 3D steady state case that I mentioned in Post #15.

    As you've noticed, we can't use a point-source for the heater like we're doing for the 1D case. So I'm assuming the heat is uniformly generated in a "source" sphere of radius rs. And for a boundary condition, the temperature rise is zero at r=R.

    For r > rs, the temperature has a 1/r dependence:

    T = P/(4 pi k) * (1/r - 1/R)

    For r < rs, the temperature is ~ (1 - r2):

    T = T(rs) + P/(8 pi k) * (1/rs) * [1 - (r/rs)2]

    I was surprised to find that (***) roughly 1/3 of the total temperature rise occurs within r < rs. I imagine that this temperature change sets in on a fast time scale (~alpha*rs2 ?), and the remaining 2/3 of the rise happens more slowly as the heat spreads into the surrounding medium.

    *** Assuming R >> rs, so that 1/R is negligible compared to 1/rs.

    Hmmm. Weird.
    Last edited: Nov 3, 2008
  21. Nov 10, 2008 #20
    Dear All,

    Thank you very much for your help. I got some answers to my questions and borrowed the math book Integral referenced.

  22. Nov 11, 2008 #21
    Hello everybody!

    Fun, as I solved this in a 1-D case for my own needs. I couldn't find the solution in a book.

    The general solution is a function of x/sqrt(t) which is multiplied by sqrt(t). In other words, the depth as well as the surface temperature increase as sqrt(t). So the energy increases as t (consistent with a constant power density) and the temperature gradient is constant (same remark).

    Then, this function of a single variable involves the integral of erfc, the complementary error function. The advantage is that several software (Mathcad or Maple) have this integral built-in.

    I first got a complicated solution by convolving erfc functions over time to meet the boundary condition of constant power density. Then, after losing one week, I found a solution on the Web: the guy used a Laplace transform with fractional powers... Not simple, but he got farther to a simpler expression. Finally, I did it again with simpler functions and got the simpler solution more easily.

    I'll have a look in my archives to put the solution here. Please be patient!

    Be aware that I had to assume constant heat capacity and conductivity. This precludes any change of state and needs the absolute temperature to vary relatively little.

    I didn't look for 2D nor 3D solutions.
  23. Nov 27, 2008 #22
    Hello there
    If anyone can help me please, I am doing modelling on a small container exposed to a fire from distance x meter, now I am trying to calculate the heat conducted through the container shell I am trying to use finite difference explicit method in one dimension heat conduction,, the container is half filled with fuel and the other half with the fuel vapour.

    If anyone can help please what should I consider the boundaries from the both sides (outside and inside the container)
    I tried to consider the outside derivative boundary as :

    dT/dx=the radiation equation

    And from inside as:
    dT/dx= radiation +convection

    Is what I am doing correct
    Because I am still confused about the first and the last nodes equations how can I get them from these boundaries if the boundaries are correct in the first place.
    I really appreciate the help
    Thanks alot

    Attached Files:

Share this great discussion with others via Reddit, Google+, Twitter, or Facebook