Уравнение теплопроводности в MATLAB

[Constructor]

Здравствуйте. Прошу вашей помощи в решении следующей задачи.

Необходимо смоделировать конкретный физический процесс: нагревание полупроводникового образца (германия) вследствие фото- и электролюминесценции. Для этого я использую MATLAB и уравнение теплопроводности.

1. Общий вид уравнения:
$\frac{\partial u}{\partial t} = \ae \cdot \frac{\partial^2 u}{\partial x^2} + f(t)$,
где $\ae = \frac{k}{c\cdot \rho }$ — коэффициент температуропроводности; k — теплопроводность вещества, c — удельная (изобарная) теплоёмкость, $\rho $ — плотность вещества.
2. Фотолюминесценция
a) В случае фотолюминесценции $f(t) = 0$. Исходная температура образца $T = 298 K$.
b) Граничные условия следующие:
В точке $x = 0$ падает модулированное лазерное излучение (меандр со скважностью 2 и частотой $\upsilon = 80$ Гц). Мощность излучения $P = 40$ мВт.
В точке $x = l$ может быть два варианта: плохой и хороший теплоотвод. Плохой теплоотвод - вакуум, хороший теплоотвод - теплоотводящий держатель.
c) Начальные условия: равномерное значение температуры вдоль всего образца.
3. Электролюминесценция.
a) В случае электролюминесценции $f(t) $ представляет собой меандр со скважностью 2 и длительностью импульса 10 мкс (т.е. это воздействие импульсным электрическим током).
b) Граничные условия следующие:
В точке $x = 0$ имеется некий сток тепла по хладопроводу.
В точке $x = l$ имеется отражение излучения от границы, т.е. производная $\frac{\partial u}{\partial x} = 0$.
c) Начальные условия аналогичные, как и для фотолюминесценции.

Далее, в MATLAB общая запись уравнения имеет вид: $c(x,t,u,\frac{\partial u}{\partial x})\cdot \frac{\partial u}{\partial t} = x^{-m}\cdot \frac{\partial }{\partial x} (x^{m}\cdot b(x,t,u,\frac{\partial u}{\partial x}))+s(x,t,u,\frac{\partial }{\partial x})$.
В нашем случае $m=0$.

У меня получились следующие файлы-функции (самый общий вид):
1. Файл-функция для исходных условий.
function value = initial1(x)
%INITIAL1: MATLAB function M-file that specifies the initial condition for
%a PDE in time and one space dimension.
value = 298;
2. Файл-функция для граничных условий.
function [pl,ql,pr,qr] = bc1(xl,ul,xr,ur,t)
%BC1: MATLAB function M-file that specifies boundary conditions for a PDE
%in time and one space dimension.
pl = ul-298;
ql = 0;
pr = ur-298;
qr = 0;
3. Файл-функция для самого дифференциального уравнения в частных производных.
function [c,b,s]=eqn1(x,t,u,DuDx)
%EQN1: MATLAB function M-file that specifies a PDE in time and one space
%dimension.
c = 1;
b = DuDx;
s = 0;
4. Скрипт для решения дифференциального уравнения и вывода этого решения в виде трёхмерного графика.
%PDE1: MATLAB script M-file that solves and plots
%solutions to the PDE stored in eqn1.m
m = 0;
%NOTE: m=0 specifies no symmetry in the problem. Taking m=1 specifies
%cylindrical symmetry, while m=2 specifies spherical symmetry.

%Define the solution mesh
x = linspace(0,1,10);
t = linspace(0,10,40);
%Solve the PDE
u = pdepe(m,@eqn1,@initial1,@bc1,x,t);
%Plot solution
surf(x,t,u);
title('Surface plot of solution.');
xlabel('Distance x');
ylabel('Time t');

Основная проблема: каким образом соотнести конкретные физические характеристики системы и уравнение теплопроводности. Я совсем запутался.

Благодарю за внимание.

[Munin]

(Оффтоп)

У меня сильное подозрение, что коэффициент температуропроводности обозначается на самом деле "каппа", $\kappa$ или $\varkappa$, а в виде $\ae$ его рисуют на доске из-за скорости и неаккуратности.

Каким образом у вас лазерный свет и электрический ток превращаются в температуру, и почему это называется фотолюминесценцией и электролюминесценцией, я не понял. Ещё я не понял, почему вы пишете об $f(t)$ как о числовой функции, в то время как, судя по вашему уравнению, она должна быть задана по всей области задачи $f(x,t),$ $x\in[0,l].$ И почему вы условие $\frac{\partial u}{\partial x} = 0$ называете отражением излучения от границы.

[Constructor]

Munin,
благодарю за комментарии и вопросы.

Цитата:

Каким образом у вас лазерный свет и электрический ток превращаются в температуру, и почему это называется фотолюминесценцией и электролюминесценцией, я не понял.

В температуру они превращаются, поскольку идёт нагрев образца под действием либо лазерного излучения, либо проходящего через образец электрического тока. Фотолюминесценция и электролюминесценция происходят в дальней инфракрасной области спектра вследствие межзонной накачки, а затем релаксации на мелкие примесные центры.

Цитата:

Ещё я не понял, почему вы пишете об $f(t)$ как о числовой функции, в то время как, судя по вашему уравнению, она должна быть задана по всей области задачи $f(x,t),$ $x\in[0,l].$

Эта функция задает дополнительный источник в случае электролюминесценции. По x она не меняется. По временной шкале представляет собой меандр (величина этого меандра примерно 0,2 А при разности потенциалов 14 В; во втором полупериоде величина равна 0).

Цитата:

И почему вы условие $\frac{\partial u}{\partial x} = 0$ называете отражением излучения от границы.

Отражение теплового излучения. Возможно, я с терминологией напортачил. Лучше сказать "тепловой поток". Отражение идёт вследствие того, что поток через границу равен 0 (в данном случае это определяется граничным условием второго типа, т.е. $u - u_{0} = h\cdot \frac{\partial u}{\partial x}$; Поскольку температура образца $u$ равна температуре окружающей среды $u_{0}$, то и производная равна нулю). На границе в этом случае температура постоянна (условие первого типа).

[Munin]

Constructor в сообщении #538880 писал(а):

В температуру они превращаются, поскольку идёт нагрев образца под действием либо лазерного излучения, либо проходящего через образец электрического тока.

И где физическое описание этого процесса? Начать с того, где происходит этот нагрев: в объёме образца или на его поверхности? Хорошо бы формулы написать с соответствующими коэффициентами.

Constructor в сообщении #538880 писал(а):

Эта функция задает дополнительный источник в случае электролюминесценции. По x она не меняется.

Хорошо, что уточнили. А вы точно знаете, что не меняется?

Constructor в сообщении #538880 писал(а):

Отражение теплового излучения. Возможно, я с терминологией напортачил. Лучше сказать "тепловой поток".

Вообще это условие, да, описывает нечто вроде "отражения". Но это образно, например, как отражение в методе изображений в электростатике, где такое условие накладывается на уравнение Лапласа. Какой физический смысл у этого условия в случае теплопроводности, надо рассматривать отдельно. Кажется, оно соответствует нулевому потоку тепла через границу, абсолютной теплоизоляции или линии симметрии системы.

[Constructor]

Munin

Цитата:

И где физическое описание этого процесса? Начать с того, где происходит этот нагрев: в объёме образца или на его поверхности? Хорошо бы формулы написать с соответствующими коэффициентами.

Нагрев происходит на поверхности образца (в случае фотолюминесценции) и во всём его объёме (в случае электролюминесценции).
1. Фотолюминесценция.
Образец закреплён и подвергается воздействию импульсного лазерного излучения при комнатной температуре.
$\frac{\partial u}{\partial t} = 3.625210 \cdot 10^{-5} \cdot \frac{\partial^2 u}{\partial x^2}$,
т.е. $\ae = 3.625210 \cdot 10^{-5} [\frac{m^{2}}{s}]$ - коэф. температуропроводности для германия.
$\frac{\partial u}{\partial x} (x=0) = - \frac{P(t))}{S\cdot k} $ - поток тепла от лазерного излучения.
$ P(t) = 40\cdot10^{-3} [W], S = \pi \cdot (10^{-2})^{2} [m^{2}], k = 59.9 [\frac{W}{m\cdot K}] $,
$u (x=l) = 290 [K] $ - температура держателя.
$u(x,0)=300 [K]$
2. Электролюминесценция.
$\frac{\partial u}{\partial t} = 3.625210 \cdot 10^{-5} \cdot \frac{\partial^2 u}{\partial x^2}$ +f(t).
$\frac{\partial u}{\partial x} (x=0) = 0$,
$\frac{\partial u}{\partial x} (x=l) = 0$,
$u (x=l) = 295 [K] $,
$u (x=l) = 290 [K] $,
$u(x,0)=300 [K]$,
$f(t) = V(t)$ - меандр.

Цитата:

Хорошо, что уточнили. А вы точно знаете, что не меняется?

Уверен. По оси x изменений нет, только вдоль оси t.

Цитата:

Какой физический смысл у этого условия в случае теплопроводности, надо рассматривать отдельно. Кажется, оно соответствует нулевому потоку тепла через границу, абсолютной теплоизоляции или линии симметрии системы.

Это условие ($\frac{\partial u}{\partial x} = 0$) соответствует нулевому тепловому потоку через границу.

[GraNiNi]

Constructor в сообщении #538595 писал(а):

Необходимо смоделировать конкретный физический процесс: нагревание полупроводникового образца (германия) вследствие фото- и электролюминесценции.

По-моему этими терминами Вы запутываете существо проблемы.
Нужно разделить источники поступления тепла (нагрева) и способы потери тепла (охлаждения).
1. Источник тепла - лазерный луч, потеря тепла - фотолюминисценция, конвекция, ИК излучение.
2. Источник тепла - электрический ток, потеря тепла - электролюминисценция, конвекция, ИК излучение.
Так ?

[Munin]

Constructor в сообщении #538595 писал(а):

Основная проблема: каким образом соотнести конкретные физические характеристики системы и уравнение теплопроводности. Я совсем запутался.

Ну, по итогам последнего сообщения:
- для случая 1 вы написали все формулы, связывающие величины уравнения теплопроводности и физические характеристики системы. Я надеюсь, в этом случае вы распутались?
Ах да, осталось ещё рассмотреть случай плохого теплоотвода на границе $x=l.$
- для случая 2 у вас осталось нераскрытым, что такое $f(t) = f_A\cdot V(t).$ Вам задана только форма этой функции по времени, вам нужно уточнить, как вычисляется амплитуда этой функции. У вас должна получиться величина размерности К/с, и надо её получить из характеристик образца как проводника, и пропускаемого электрического тока. Может быть, джоулево тепло, делённое на теплоёмкость?
И забыто граничное условие на $x=0,$ "некий сток тепла по хладопроводу".

GraNiNi в сообщении #539299 писал(а):

По-моему этими терминами Вы запутываете существо проблемы.
Нужно разделить источники поступления тепла (нагрева) и способы потери тепла (охлаждения).
1. Источник тепла - лазерный луч, потеря тепла - фотолюминисценция, конвекция, ИК излучение.
2. Источник тепла - электрический ток, потеря тепла - электролюминисценция, конвекция, ИК излучение.
Так ?

Если рассматривать ещё и охлаждение образца, тогда получится намного более сложное уравнение, типа
$\frac{\partial u}{\partial t} = \kappa \cdot \frac{\partial^2 u}{\partial x^2} + f(t) - g(u),$
где $g(u)$ - функция потерь тепла, зависящая от текущей температуры. Может, численно MATLAB и решит это, но у меня сомнения, что задали нечто столь неподъёмное, скорее, просто слова "электролюминесценция" и "фотолюминесценция" не имеют отношения к сути задачи.

Было произнесено следующее:

Constructor в сообщении #538880 писал(а):

Фотолюминесценция и электролюминесценция происходят в дальней инфракрасной области спектра вследствие межзонной накачки, а затем релаксации на мелкие примесные центры.

У меня складывается впечатление, что люминесценция в данном случае - только следствие нагрева, и процесс, позволяющий измерить полученную функцию температуры $u.$ А потери на люминесценцию на саму теплопередачу не влияют.

[Constructor]

Munin
В данном случае эта задача представляет собой моделирование эксперимента. Такие условия были поставлены коллегами, так сказать.

Цитата:

для случая 1 вы написали все формулы, связывающие величины уравнения теплопроводности и физические характеристики системы. Я надеюсь, в этом случае вы распутались?

Здесь я распутался. В принципе, все размерности совпадают. Только сам MATLAB странно себя ведёт: в случае фотолюминесценции, когда исходное значение температуры образца, а также значение температуры в точке $x=0$ и $x=l$ имеют одно и то же значение, то даже в случае, когда через грань $x=0$ течет тепловой поток, в конечном итоге солвер pdepe выдает константу (как по x, так и по t). А этого быть не должно.
Плохой теплоотвод в данном случае рассматривать нет необходимости: процесс происходит при комнатной температуре, а крепление для образца обладает хорошей теплопроводностью.

Цитата:

для случая 2 у вас осталось нераскрытым, что такое $f(t) = f_A\cdot V(t).$ Вам задана только форма этой функции по времени, вам нужно уточнить, как вычисляется амплитуда этой функции. У вас должна получиться величина размерности К/с, и надо её получить из характеристик образца как проводника, и пропускаемого электрического тока. Может быть, джоулево тепло, делённое на теплоёмкость?

В случае электролюминесценции, через образец пропускается электрический ток, который модулируется во времени. По поводу амплитуды стоит уточнить. Я понимаю, что размерность должна быть К/с. Думаю, Ваше предложение в этом случае достаточно закономерно. Я рассмотрю эту возможность.
Да, и в этом случае образец находится в гелии. Поэтому при $x=l$ температура равна 5 К. А на границе $x=0$ температура равна либо 5 К (гелиевая среда), либо 0 К (вакуум). При этом, как я понимаю, граничные условия являются условиями первого типа, и производные на границах равны 0. Сток тепла по хладопроводу в данном случае выражается просто температурой этого хладопровода, поскольку его размеры и теплоемкость очень велики по сравнению с аналогичными параметрами самого образца.

Цитата:

Если рассматривать ещё и охлаждение образца...

Охлаждение образца, думаю, здесь рассматривать не надо. Сначала необходимо разобраться с нагревом. Хотя, охлаждение образца в рассматриваемых условиях эксперимента пренебрежимо мало и практически не влияет на изменение температуры $u$.

GraNiNi

Цитата:

По-моему этими терминами Вы запутываете существо проблемы.
Нужно разделить источники поступления тепла (нагрева) и способы потери тепла (охлаждения).
1. Источник тепла - лазерный луч, потеря тепла - фотолюминисценция, конвекция, ИК излучение.
2. Источник тепла - электрический ток, потеря тепла - электролюминисценция, конвекция, ИК излучение.
Так ?

Охлаждение образца в данном случае пренебрежимо мало. Не думаю, что стоит сразу так усложнять систему. В рамках первой итерации нужно смоделировать простой процесс, а затем добавлять новые параметры и детали.

[Munin]

Constructor в сообщении #540692 писал(а):

В данном случае эта задача представляет собой моделирование эксперимента. Такие условия были поставлены коллегами, так сказать.

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

В любом случае, вам нужно поправить "нагревание вследствие ...люминесценции" на "нагревание в ходе ...люминесценции" или даже "для".

Constructor в сообщении #540692 писал(а):

в конечном итоге солвер pdepe выдает константу (как по x, так и по t). А этого быть не должно.

А почему вы думаете, что не должно? В принципе, в уравнении теплопроводности скорость распространения тепла бесконечная, так что может быть, сколько втекает, столько и вытекает. Проверьте, поставив плохой теплоотвод, например. Или аналитически прикинув искомый вид решения.

Constructor в сообщении #540692 писал(а):

Хотя, охлаждение образца в рассматриваемых условиях эксперимента пренебрежимо мало

Это как раз и означает, что его рассматривать не надо :-)

Я думаю, вы переходите от стадии "соотнести физические характеристики и что запихивать в солвер" к стадии "мучаться с солвером, пытаясь получить от него что-то осмысленно выглядящее" :-) После неё светит третья "получать от солвера искомые результаты", так что всё оптимистично.

[ewert]

Constructor в сообщении #538880 писал(а):

Отражение идёт вследствие того, что поток через границу равен 0 (в данном случае это определяется граничным условием второго типа, т.е. $u - u_{0} = h\cdot \frac{\partial u}{\partial x}$; Поскольку температура образца $u$ равна температуре окружающей среды $u_{0}$, то и производная равна нулю). На границе в этом случае температура постоянна (условие первого типа).

Тут какая-то сплошная путаница. $u - u_{0} = h\cdot \frac{\partial u}{\partial x}$ -- это условие третьего типа. И если вдруг производная уверенно равна нулю, то это вовсе не от того, что какая-то там температура равна температуре, а просто потому, что вот теплоизоляция -- и всё тут. И эта теплоизоляция противоречит постоянству температуры.

[Constructor]

Munin

Цитата:

Я думаю, вы переходите от стадии "соотнести физические характеристики и что запихивать в солвер" к стадии "мучаться с солвером, пытаясь получить от него что-то осмысленно выглядящее" После неё светит третья "получать от солвера искомые результаты", так что всё оптимистично.

Я уже чёрт знает сколько времени мучаюсь с солвером. Конечный вид решения я себе представляю и понимаю, каким он должен быть. Только вот солвер выдает такие решения, что хочется ноутбук в окно выкинуть :) Хотя, компьютер лишь совершает те операции, на которые его запрограммировал человек. Поэтому проблема, видимо, во мне :) Пока ещё не всё встало на свои места.

ewert

Цитата:

Тут какая-то сплошная путаница. $u - u_{0} = h\cdot \frac{\partial u}{\partial x}$ -- это условие третьего типа. И если вдруг производная уверенно равна нулю, то это вовсе не от того, что какая-то там температура равна температуре, а просто потому, что вот теплоизоляция -- и всё тут. И эта теплоизоляция противоречит постоянству температуры.

Я понимаю, что это условие третьего типа. Просто я имел ввиду саму производную как условие второго типа.
То есть, как я понимаю, теплоизоляция - это когда $\frac{\partial u}{\partial x} = 0$, однако, при этом на границе может наблюдаться скачок температур $u - u_{0} \neq 0$?

[Munin]

Constructor в сообщении #541345 писал(а):

Я уже чёрт знает сколько времени мучаюсь с солвером. Конечный вид решения я себе представляю и понимаю, каким он должен быть. Только вот солвер выдает такие решения, что хочется ноутбук в окно выкинуть :) Хотя, компьютер лишь совершает те операции, на которые его запрограммировал человек. Поэтому проблема, видимо, во мне :) Пока ещё не всё встало на свои места.

Попробуйте дать солверу самую простую задачку, с тем чтобы потом усложнить её условия до реальных. Добейтесь, чтобы он её решал. Потом усложняйте условия небольшими шагами. Заодно приглядывайтесь к документации солвера, что там сказано про специфику его численных методов, и как с ней управляться. Уравнение теплопроводности обычно более смирно себя ведёт при уменьшении шага по времени, до порядка $\Delta t\sim\Delta x^2/\kappa$ (подробнее см. литературу по численным методам).

[Constructor]

Munin

Цитата:

Попробуйте дать солверу самую простую задачку, с тем чтобы потом усложнить её условия до реальных. Добейтесь, чтобы он её решал. Потом усложняйте условия небольшими шагами. Заодно приглядывайтесь к документации солвера, что там сказано про специфику его численных методов, и как с ней управляться. Уравнение теплопроводности обычно более смирно себя ведёт при уменьшении шага по времени, до порядка $\Delta t\sim\Delta x^2/\kappa$ (подробнее см. литературу по численным методам).

Я экспериментирую. Только пока неудачно. Простые примеры солвер решает нормально. А сложные запутывает до невозможности.
Попробую решить эту проблему. Благодарю за помощь.

[Munin]

Какой самый простейший пример, который не удаётся решить?

[Constructor]

Munin

На данный момент ситуация следующая: у меня получилось правильно записать граничные условия. Однако теперь появилась проблема в устойчивости и сходимости. MATLAB порой выдаёт довольно странные варианты решения. Простые примеры солвер тоже решает.

Допустим, простой пример:
Исходное уравнение: $\frac{\partial u}{\partial t} = \frac{\partial^2 u}{\partial x^2};$
Начальное условие: $u (x,0) = 293;$
Граничное условие 1: $u (0,t) = 293;$
Граничное условие 2:$u (x,t) = 293$, при этом $x = 1$ ;
Поток через грань 1: $\frac{\partial u}{\partial x} (x = 0) = - sin(2 \cdot pi \cdot t)$
Однако, даже здесь сложно подобрать область устойчивости. Соотношение $\tau \leqslant 0.5 \cdot h^{2}$ не всегда действует. Аналогично и в более сложном примере.
Можно вместо синуса использовать меандр. Но при этом появляются провалы, которые объяснить сложно. Видимо, это связано с той же устойчивостью. Если менять шаг по t и x, то изменения происходят. Однако выявить закономерность не представляется возможным.