Тепловое моделирование светодиодов (как решить уравнение теплопроводности с постоянным источником тепла)

У меня есть механическая конструкция со светодиодами, которые выделяют тепло. Я хочу оценить температуру на переходе светодиода в зависимости от времени, но особенно в установившемся режиме.

Зная падение напряжения и ток светодиода, я могу оценить мощность, рассеиваемую за счет тепла (некоторая часть мощности уходит на свет, но, чтобы быть консервативным, я решил оценить тепловую мощность как просто I * V).

Температура перехода светодиода должна поддерживаться ниже определенной температуры. Тепловая система будет состоять из соединения светодиода с термопрокладкой, печатной платы и механической конструкции. Тепловое сопротивление указано для соединения с термопрокладкой (около 10 К/Вт), и производитель предлагает несколько различных конструкций печатных плат, каждая из которых имеет свое собственное тепловое сопротивление. Наилучшая конструкция может достигать 3-4 кВт/Вт, но возможна и менее дорогая конструкция. Поэтому мне нужно тепловое сопротивление механического корпуса, которое уже существует и не может быть изменено, кроме как прикрепить светодиоды к корпусу.

Геометрия несколько сложна, поэтому я решил начать с простого алюминиевого блока (печатная плата крепится непосредственно к блоку), чтобы убедиться, что мое моделирование правильное, а затем перейти к более сложной геометрии. Светодиоды расположены в линию, поэтому я буду предполагать, что теплопередача постоянна вдоль линии светодиодов (назовем ее осью Z). Я также не хочу сейчас рассматривать конвекцию или излучение по длине блока.

Итак, скажем, я прикрепляю радиатор к концу блока напротив светодиодов с достаточно низким тепловым сопротивлением окружающей среде, чтобы конец блока с радиатором находился на температуре окружающей среды. Термическое сопротивление алюминиевого блока составляет L/(кА), где L – длина, k – теплопроводность (0,25 Вт/(мм*K)), а A – площадь поперечного сечения. Таким образом, если длина равна 20 мм (по оси X), площадь поперечного сечения составляет около 300 мм ^ 2 (около 6 мм x 50 мм, 6 мм по оси Y и 50 мм по оси Z), тепловое сопротивление составляет около 0,27 К/Вт. Мощность светодиода на единицу площади составляет 0,0896 Вт/мм^2 (общая мощность 26,88 Вт, распределенная по 300 мм^2).

Примечание: мне нравится работать в мм для этой задачи.

Эта задача является одномерной, и, согласно закону проводимости Фурье, падение температуры на блоке в установившемся режиме должно быть:

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

Я хочу иметь возможность смоделировать это и получить тот же результат.

Поскольку окончательная геометрия более сложная, я считаю, что мне нужно использовать метод конечных элементов для окончательного решения, но для начала было бы неплохо знать, как начать с уравнения теплопроводности и закончить с результатом, полученным из уравнения Фурье. закон (поскольку первое может быть получено из второго).

Я нашел freefem++ ( http://www.freefem.org/ff++/index.htm ), который может аппроксимировать решение уравнения теплопроводности с помощью метода конечных элементов, но вариационная формулировка находится за пределами моего понимания, и ни один из примеров не касается с постоянным источником тепла.

Я хочу настроить уравнение теплопроводности и граничные условия. Мне особенно нужен стационарный результат, но меня также интересует изменение температуры в зависимости от времени. (3-D) уравнение теплопроводности:

$${{\partial u} \over {\partial t}} = \alpha \nabla ^2 u + {Q \over {c \rho}},$$ куда

$u = u(x,y,z,t)$= температура$(K)$

$\alpha = {k \over c \rho}$= температуропроводность$({mm^2 \over s})$

$c$= удельная теплоемкость$({J \over {kg \cdot K}})$

$\rho$= массовая плотность$({kg \over mm^3})$

$Q$= мощность на единицу объема$({W \over mm^3})$

Поскольку простая блочная задача является одномерной, уравнение теплопроводности принимает вид:

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

или же

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

Обратите внимание, Q является функцией$(x, t)$. В этой статье представлено аналитическое решение этой одномерной задачи: http://people.math.gatech.edu/~xchen/teach/pde/heat/Heat-Duhamel.pdf .

Чтобы использовать принцип Дюамеля, я сначала переформулирую задачу так, чтобы источник тепла находился на$L/2$температура при$x=0$а также$x=L$удерживается равным нулю (будь то ноль или температура окружающей среды, решение остается тем же). Тогда источником тепла является дельта-функция Дирака при$x=L/2$, и проблема становится:

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

и граничные условия:

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

Здесь,$Q_i = 0.896 W/mm^3$является источником тепла. Я смог аппроксимировать решение, приведенное в статье, с помощью Matlab (см. Ниже), и как$ds \rightarrow 0$, решение становится ближе к функции треугольника с центром в$x=L/2$и увеличивается асимптотически со временем примерно до 7 К. Обратите внимание, что, поскольку источник тепла находится в центре, я удвоил мощность и длину блока, и решение вышеуказанной задачи составляет половину блока ($x \ge L/2$).

Итак, что происходит, когда я перехожу в 2-D? Как я уже сказал, реальная геометрия сложна (но может быть смоделирована как двумерная задача, поскольку геометрия вдоль оси z практически постоянна), поэтому я не ожидаю аналитического решения. Возможно, моего простого одномерного примера недостаточно, чтобы продемонстрировать, как решить задачу с помощью метода конечных элементов.

Что, если бы источник тепла поступал слева, а нижняя часть блока поддерживалась при постоянной температуре?

Распространяется ли принцип Дюамеля на 2D? Если да, и я аппроксимирую вспомогательную однородную задачу с помощью МКЭ, то как мне преобразовать это приближение в неоднородное решение?

В качестве альтернативы, как бы я сформулировал неоднородную задачу, используя вариационную постановку, чтобы иметь возможность напрямую использовать анализ FEM?

Приближение Matlab к одномерному аналитическому решению

Вот мой код Matlab, который приближает решение с использованием принципа Дюамеля. Я сделал приближение, используя как ряд Фурье, так и функцию Грина.

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