# Heat equation with const. heat source

1. Oct 29, 2008

### Franck

Hi,

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.

Franck

2. Oct 29, 2008

### nasu

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.

3. Oct 29, 2008

### Franck

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

4. Oct 29, 2008

### Redbelly98

Staff Emeritus
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.

5. Oct 30, 2008

### Franck

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.

6. Oct 30, 2008

### Redbelly98

Staff Emeritus
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?

7. Oct 30, 2008

### Franck

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.
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
8. Oct 30, 2008

### Integral

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

$$U_t(x,t) = kU_{xx}(x,t)$$
Initial condition:
U(x, +0)=0
with U(x,t) =0 as x -> $\infty$
and

$$K U_x(0,t) = \phi_0$$

He has for a solution:

$$U(0,t) = \frac {2 \phi_0 } K \sqrt {\frac {kt} \pi$$

9. Oct 30, 2008

### Redbelly98

Staff Emeritus
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.

10. Oct 31, 2008

### Franck

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
11. Oct 31, 2008

### Redbelly98

Staff Emeritus
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)
and
α = k / (ρ cp)

12. Oct 31, 2008

### Redbelly98

Staff Emeritus
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
Density
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.

EDIT:
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
13. Oct 31, 2008

### Redbelly98

Staff Emeritus
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 \$$

where

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
14. Nov 1, 2008

### Franck

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.

15. Nov 2, 2008

### Redbelly98

Staff Emeritus
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)$$

where
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
and/or
t2 = (c ρ / k) R2 ,
which are vastly different from each other.

(It would be nice if this could be checked somehow with a simulation.)

16. Nov 3, 2008

### Franck

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)
$$\rho$$=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.

17. Nov 3, 2008

### Franck

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. =>
rs=0.39um
R=10um
P=0.03W
and thermal parameters as above

With this I get

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

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

Last edited by a moderator: May 3, 2017
18. Nov 3, 2008

### Redbelly98

Staff Emeritus
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

19. Nov 3, 2008

### Redbelly98

Staff Emeritus
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:

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
20. Nov 10, 2008

### Franck

Dear All,

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

Thanks
Franck