I have a mechanical design with LEDs that generate heat. I want to estimate the temperature at the LED junction vs. time, but especially at steady state.
Knowing the LED voltage drop and current, I can estimate the power dissipated via heat (some of the power goes to light, but to be conservative, I thought to estimate heat power as simply I*V).
The LED junction temperature needs to be held below a certain temperature. The thermal system will consist of the LED junction to thermal pad, the circuit board, and the mechanical design. The thermal resistance is given for the junction to thermal pad (about 10 K/W), and the manufacturer provides a couple different circuit board designs, each with their own thermal resistance. The best design can achieve 3-4 K/W, but a less expensive design may be an option. So I need the thermal resistance of the mechanical housing, which is pre-existing, and cannot be changed except for how to attach the LEDs to the housing.
The geometry is somewhat complicated, so I thought I'd start with a simple aluminum block (the circuit board mounts directly to the block), to make sure my modeling is correct, and then move to more complicated geometries. The LEDs are arranged in a line, so I'm going to assume the heat transfer is constant along the LED line (call it the z-axis). I also don't want to consider convection or radiation right now, along the length of the block.
So let's say I attach a heat sink to the end of the block opposite the LEDs, with a low enough thermal resistance to ambient so the heat sink end of the block is at ambient. The thermal resistance of an aluminum block is L/(kA), where L is the length, k is the thermal conductivity (0.25 W/(mm*K)), and A is the cross-sectional area. So if the length is 20 mm (along x-axis), the cross-sectional area is about 300 mm^2 (about 6 mm x 50 mm, with 6 mm along y-axis, and 50 mm along z-axis), the thermal resistance is about 0.27 K/W. The LED power per unit area is 0.0896 W/mm^2 (total power of 26.88 W, spread across the 300 mm^2).
Note: I like to work in mm for this problem.
This problem is 1-D, and from Fourier's conduction law, the temperature drop across the block at steady state should therefore be:
$\Delta T = L \cdot P/(kA) = 7.17 K.$
I want to be able to model this and come up with the same result.
Because the final geometry is more complex, I believe I need to use the Finite Element Method for the final solution, but for starters, it'd be nice to know how to start with the heat equation and end up with the result obtained from Fourier's law (since the former can be derived from the latter).
I found freefem++ (http://www.freefem.org/ff++/index.htm), which can approximate a solution to the heat equation using the Finite Element Method, but variational formulation is beyond my understanding, and none of the examples deal with a constant heat source.
I want to setup the heat equation and boundary conditions. I especially want the steady-state result, but I'm also interested in the temperature evolution vs. time. The (3-D) heat equation is
$${{\partial u} \over {\partial t}} = \alpha \nabla ^2 u + {Q \over {c \rho}},$$ where
$u = u(x,y,z,t)$ = temperature $(K)$
$\alpha = {k \over c \rho}$ = thermal diffusivity $({mm^2 \over s})$
$c$ = specific heat $({J \over {kg \cdot K}})$
$\rho$ = mass density $({kg \over mm^3})$
$Q$ = power per unit volume $({W \over mm^3})$
Since the simple block problem is 1-D, the heat equation becomes:
$${{\partial u} \over {\partial t}} = \alpha {{\partial ^2 u} \over {\partial x^2}} + {Q \over {c \rho}}$$
or
$$u_t = \alpha u_{xx} + {Q \over {c \rho}}$$
Note, Q is a function of $(x, t)$. This paper provides an analytical solution for this 1-D problem: http://people.math.gatech.edu/~xchen/teach/pde/heat/Heat-Duhamel.pdf
In order to use Duhamel's Principle, I first reformulate the problem so the heat source is at $L/2$ the temperature at $x=0$ and $x=L$ is held to zero (whether zero or ambient temperature, the solution is the same). Then the heat source is the Dirac-delta function at $x=L/2$, and the problem becomes:
$$u_t = \alpha u_{xx} + {Q_i \over {c \rho}} \delta(x-{L \over 2}), 0 < x < L, t > 0$$
and the boundary conditions are:
$$u(0,t) = 0, u(L, t) = 0, t > 0$$ $$u(x,0) = 0, 0 \le x \le L$$
Here, $Q_i = 0.896 W/mm^3$ is the heat source. I was able to approximate the solution given in the paper using Matlab (see below), and as $ds \rightarrow 0$, the solution becomes closer to a triangle function centered at $x=L/2$ and increasing asymptotically over time to about 7 K. Note that because the heat source is in the center, I doubled the power and block length, and the solution for the above problem is one half of the block ($x \ge L/2$).
So now what happens when I go to 2-D? As I said, the real geometry is complex (but can be modeled as a 2-D problem, since the geometry along the z-axis is essentially constant), so I don't expect an analytical solution. Perhaps my simple 1-D example is not sufficient to demonstrate how to solve the problem with the Finite Element Method.
What if the heat source came in from the left and the bottom of the block was held at constant temperature?
Does the Duhamel Principle extend to 2-D? If so, and I approximate the auxiliary homogeneous problem with FEM, then how do I transform that approximation into the nonhomogenous solution?
Alternately, how would I formulate the nonhomogenous problem using variational formulation, in order to be able to use do FEM analysis directly?
Matlab approximation to 1-D analytical solution
Here is my Matlab code that approximates the solution using Duhamel's Principle. I did the approximation using both Fourier series and Green's function.
% Approximate the analytical solution of the heat equation with a heat
% source in the center of a block.
% System parameters.
H = 6; % the block height (mm)
L = 40; % the block length (mm)
W = 50; % the block width (mm)
kAl = 0.25; % Aluminum thermal conductivity (W/(mm*K))
c = 897; % Aluminum specific heat capacity (J/(kg * K)).
rho = 2.7E-6; % Density (kg/mm^3).
alpha = kAl / (c * rho); % Thermal diffusivity (mm^2/s).
Qi = 2 * 27 / 300; % Input power per unit volume length (?).
dx = 0.2;
dt = .2;
x = 0:dx:L;
tmax = 10;
t = 0:dt:tmax;
% Approximate heat equation using Fourier series and Duhamel's Principle.
ds = 0.1;
N = 200;
n = 1:N;
b = 2*Qi*sin(n*pi/2)/(c*rho*L);
% As N goes to infinity, the solution
% approximates a triangle function centered on L/2. Because we can't go to
% infinity, there will always be a sharp spike at x = L/2.
u = zeros(length(x), length(t));
for xi = 1:length(x)
for ti = 1:length(t)
tc = t(ti);
for ni = 1:length(n)
s = 0:ds:tc;
sint = 0;
for si = 1:length(s)
sint = sint + b(ni)*exp(-alpha*(n(ni)*pi/L)^2*(tc-s(si)))*ds;
end
u(xi, ti) = u(xi, ti) + sin(n(ni)*pi*x(xi)/L) * sint;
end
end
end
figure;
mesh(t, x, u);
ylabel('x (mm)');
xlabel('t (s)');
zlabel('Temperature (deg C)');
title('Approximation to heat equation solution with constant heat source at L/2, using Fourier series');
% Approximate solution using Green's function. Note that as ds -> zero,
% the solution approximates a triangle function centered at L/2, and
% increasing asymptotically over time.
u = zeros(length(x), length(t));
N = 40;
n = -N:N;
ds = 0.01;
for xi = 1:length(x)
for ti = 1:length(t)
tc = t(ti);
if tc == 0
continue;
end
s = 0:ds:(tc-ds);
for si = 1:length(s)
nint = 0;
for ni = 1:length(n)
nint = nint + exp(-(x(xi)-2*n(ni)*L-L/2)^2/(4*alpha*(tc-s(si)))) - ...
exp(-(x(xi)-2*n(ni)*L+L/2)^2/(4*alpha*(tc-s(si))));
end
u(xi, ti) = u(xi, ti) + ...
(Qi/(c*rho)) * nint * ds / sqrt(4*pi*alpha*(tc-s(si)));
end
end
end
figure;
mesh(t, x, u);
ylabel('x (mm)');
xlabel('t (s)');
zlabel('Temperature (deg C)');
title('Approximation to heat equation solution with constant heat source at L/2, using Green''s function');
