Modelado térmico LED (Cómo resolver la ecuación de calor con una fuente de calor constante)

Tengo un diseño mecánico con LEDs que generan calor. Quiero estimar la temperatura en la unión LED frente al tiempo, pero especialmente en estado estable.

Conociendo la caída de voltaje y la corriente del LED, puedo estimar la potencia disipada a través del calor (parte de la potencia se va a la luz, pero para ser conservador, pensé en estimar la potencia térmica simplemente como I * V).

La temperatura de la unión LED debe mantenerse por debajo de cierta temperatura. El sistema térmico consistirá en la unión del LED a la almohadilla térmica, la placa de circuito y el diseño mecánico. La resistencia térmica se proporciona para la unión a la almohadilla térmica (alrededor de 10 K/W), y el fabricante proporciona un par de diseños de placa de circuito diferentes, cada uno con su propia resistencia térmica. El mejor diseño puede alcanzar 3-4 K/W, pero un diseño menos costoso puede ser una opción. Entonces necesito la resistencia térmica de la carcasa mecánica, que es preexistente y no se puede cambiar, excepto cómo conectar los LED a la carcasa.

La geometría es algo complicada, así que pensé en comenzar con un bloque de aluminio simple (la placa de circuito se monta directamente en el bloque), para asegurarme de que mi modelado sea correcto, y luego pasar a geometrías más complicadas. Los LED están dispuestos en una línea, por lo que supondré que la transferencia de calor es constante a lo largo de la línea LED (llámese eje z). Tampoco quiero considerar la convección o la radiación en este momento, a lo largo del bloque.

Entonces, digamos que adjunto un disipador de calor al extremo del bloque opuesto a los LED, con una resistencia térmica lo suficientemente baja al ambiente para que el extremo del disipador de calor del bloque esté en el ambiente. La resistencia térmica de un bloque de aluminio es L/(kA), donde L es la longitud, k es la conductividad térmica (0,25 W/(mm*K)) y A es el área de la sección transversal. Entonces, si la longitud es de 20 mm (a lo largo del eje x), el área de la sección transversal es de aproximadamente 300 mm ^ 2 (alrededor de 6 mm x 50 mm, con 6 mm a lo largo del eje y y 50 mm a lo largo del eje z), la resistencia térmica es de unos 0,27 K/W. La potencia del LED por unidad de área es de 0,0896 W/mm^2 (potencia total de 26,88 W, distribuida en 300 mm^2).

Nota: Me gusta trabajar en mm para este problema.

Este problema es 1-D y, según la ley de conducción de Fourier, la caída de temperatura en el bloque en estado estacionario debería ser:

$\Delta T = L \cdot P/(kA) = 7.17 K.$

Quiero poder modelar esto y obtener el mismo resultado.

Debido a que la geometría final es más compleja, creo que necesito usar el método de elementos finitos para la solución final, pero para empezar, sería bueno saber cómo comenzar con la ecuación del calor y terminar con el resultado obtenido de Fourier. ley (ya que la primera puede derivarse de la segunda).

Encontré freefem++ ( http://www.freefem.org/ff++/index.htm ), que puede aproximar una solución a la ecuación del calor usando el método de elementos finitos, pero la formulación variacional está más allá de mi comprensión, y ninguno de los ejemplos trata con una fuente de calor constante.

Quiero configurar la ecuación de calor y las condiciones de contorno. Especialmente quiero el resultado de estado estacionario, pero también estoy interesado en la evolución de la temperatura frente al tiempo. La ecuación de calor (3-D) es

$${{\partial u} \over {\partial t}} = \alpha \nabla ^2 u + {Q \over {c \rho}},$$ dónde

$u = u(x,y,z,t)$= temperatura$(K)$

$\alpha = {k \over c \rho}$= difusividad térmica$({mm^2 \over s})$

$c$= calor específico$({J \over {kg \cdot K}})$

$\rho$= densidad de masa$({kg \over mm^3})$

$Q$= potencia por unidad de volumen$({W \over mm^3})$

Dado que el problema de bloque simple es 1-D, la ecuación de calor se convierte en:

$${{\partial u} \over {\partial t}} = \alpha {{\partial ^2 u} \over {\partial x^2}} + {Q \over {c \rho}}$$

o

$$u_t = \alpha u_{xx} + {Q \over {c \rho}}$$

Tenga en cuenta que Q es una función de$(x, t)$. Este documento proporciona una solución analítica para este problema 1-D: http://people.math.gatech.edu/~xchen/teach/pde/heat/Heat-Duhamel.pdf

Para usar el Principio de Duhamel, primero reformulo el problema para que la fuente de calor esté en$L/2$la temperatura a$x=0$y$x=L$se mantiene en cero (ya sea cero o temperatura ambiente, la solución es la misma). Entonces la fuente de calor es la función Dirac-delta en$x=L/2$, y el problema se convierte en:

$$u_t = \alpha u_{xx} + {Q_i \over {c \rho}} \delta(x-{L \over 2}), 0 < x < L, t > 0$$

y las condiciones de contorno son:

$$u(0,t) = 0, u(L, t) = 0, t > 0$$ $$u(x,0) = 0, 0 \le x \le L$$

Aquí,$Q_i = 0.896 W/mm^3$es la fuente de calor. Pude aproximar la solución dada en el documento usando Matlab (ver más abajo), y como$ds \rightarrow 0$, la solución se acerca más a una función triangular centrada en$x=L/2$y aumentando asintóticamente con el tiempo hasta aproximadamente 7 K. Tenga en cuenta que debido a que la fuente de calor está en el centro, dupliqué la potencia y la longitud del bloque, y la solución para el problema anterior es la mitad del bloque ($x \ge L/2$).

Entonces, ¿qué pasa cuando voy a 2-D? Como dije, la geometría real es compleja (pero se puede modelar como un problema en 2D, ya que la geometría a lo largo del eje z es esencialmente constante), por lo que no espero una solución analítica. Quizás mi simple ejemplo 1-D no sea suficiente para demostrar cómo resolver el problema con el método de elementos finitos.

¿Qué pasaría si la fuente de calor entrara por la izquierda y la parte inferior del bloque se mantuviera a temperatura constante?

¿El principio de Duhamel se extiende a 2-D? Si es así, y aproximo el problema homogéneo auxiliar con FEM, ¿cómo transformo esa aproximación en la solución no homogénea?

Alternativamente, ¿cómo formularía el problema no homogéneo usando la formulación variacional, para poder usar el análisis FEM directamente?

Aproximación de Matlab a la solución analítica 1-D

Aquí está mi código de Matlab que se aproxima a la solución usando el Principio de Duhamel. Hice la aproximación utilizando tanto la serie de Fourier como la función de Green.

% 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');