Heat equation with const. heat source

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 t^0.58 (1D); t^0.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
Franck

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.

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

Redbelly98

Something is weird about the suggested power dependences t^0.58 and t^0.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:

[tex]
\frac{dT}{dx} = \frac{P}{k A}
[/tex]

where P is the heater power, k is the rod's thermal conductivity, and A is its cross-sectional area.

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.

Redbelly98

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:

T_max=T_sink+(P L) / (k A)

However, if it follows the t^0.58 dependence, eventually it would exceed that temperature. So this would be a good test of the t^0.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?

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 T_max
Have a look at the graph.
http://picasaweb.google.com/franck.h...54245291443458
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).

Integral

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

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

Redbelly98

Thanks Integral. The t^1/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.

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

Redbelly98

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 / (ρ c_p)

Redbelly98

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.

Redbelly98

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

"T" is the temperature rise above ambient.

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

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/π²)
α = k / (ρ c_p)
ρ is the density
c_p 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→∞

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.

Redbelly98

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:

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

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

Assuming r_s << 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 T_max, but dimensional analysis suggests time scales on the order of

t1 = (c ρ / k) r_s²
and/or
t2 = (c ρ / k) R² ,
which are vastly different from each other.

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

Franck

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/cm³ (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= 100um² (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.

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 r_s by fitting Redbelly's T_max. =>
r_s=0.39um
R=10um
P=0.03W
and thermal parameters as above

With this I get

t1 = (c ρ / k) r_s²=7ns
and/or
t2 = (c ρ / k) R² =4.8us

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


http://picasaweb.google.com/franck.h...20670106402930