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:

<br /> \frac{dT}{dx} = \frac{P}{k A}<br />

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.hertz/Forum#5262954245291443458
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

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

 

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

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

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:

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

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)
\rho=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.hertz/Forum#5264420670106402930

 

[Redbelly98]

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

 

[Redbelly98]

1D simulation

k=1.9 W/cmK (thermal conductivity)
\rho=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).

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 r_s. And for a boundary condition, the temperature rise is zero at r=R.

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

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

For r < r_s, the temperature is ~ (1 - r²):

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

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

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

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

Hmmm. Weird.

 

[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

 

[Enthalpy]

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.

 

[yahya2003]

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