Решение уравнения теплопроводности в программе

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

Уравнение теплопроводности для двумерной среды в конечных разностях имеет вид:

Для расчета распределения температуры по поверхности пластины используется программа ПР-1. Она содержит цикл по времени с вложенными в него двумя циклами по i и по j, в которых перебираются все элементы пластины и вычисляются их температуры в последующие моменты времени. Результат решения задачи представлен на рис. 1, — на нем разными цветами изображены области с различными температурами. Изотермы (границы разноцветных областей) перпендикулярны теплоизолированным краям пластины и параллельны краям, температура которых поддерживается постоянной.

Рис. 1. Распределение температуры: двумерная среда

Задача 2.

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

Заменяя частные производные их конечно-разностными аппроксимациями, запишите уравнение теплопроводности для неоднородной двумерной среды в конечных разностях:

Для решения задачи используется программа ПР-2. В ней последовательно перебираются узлы двумерной сетки по строкам и по столбцам. Когда переменная naprav равна 0, то элементы перебираются по столбцам сверху вниз и снизу вверх. Когда переменная naprav равна 1, то элементы перебираются по строкам слева направо и справа налево. Нижний край пластины теплоизолирован, это задается циклом:

Все остальные края пластины поддерживаются при постоянной температуре. Результат работы программы представлен на рис. 2.

Рис. 2. Теплопроводность в неоднородной среде

Задача 3.

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

Для решения задачи может быть использован алгоритм АЛ-1.

В программе ПР-3 используются два массива T[i] и TT[i], в которых сохраняются значения температуры элементов стержня в моменты t и t+1 соответственно. Расчет температуры осуществляется по выведенной выше формуле в цикле. Будем считать, что длина стержня l, его коэффициент температуропроводности k при x>0,2l равен 1,8, а при x

Рис. 3. Теплопроводность стержня

Задача 4.

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

Эта нестационарная задача требует решения следующего уравнения теплопроводности в полярных координатах:

В конечных разностях получаем:

Для решения этой задачи используется программа ПР-4. Он содержит два вложенных цикла по i и j, в которых перебираются все узлы двумерной сетки и пересчитываются значения T[i,j] на следующем временном слое. При ее запуске на экране появляется цветное изображение, границы одноцветных областей соответствуют изотермам. Пример результата вычислений приведен на рис. 4.

Рис. 4. Решение уравнения теплопроводности в полярных координатах.

Тексты программ находятся в zip-архиве, файл gl-7.pas.

Метод Finite Volume — реализация на примере теплопроводности

Метод Finite Volume (FVM)

В основе метода лежит разбиение области на непересекающиеся контрольные объемы(элементы), узловые точки, в которых ищется решение.Узловые точки находятся в центрах контрольных объемов.Также, как и для метода конечных разностей, для каждого элемента составляется уравнение, получается система линейных уравнений.Решая ее — находим значения
искомых переменных в узловых точках.Для отдельного элемента уравнение получается путем интегрирования исходного дифф уравнения по элементу и аппроксимации интегралов.

Термин конечный объем в статье будет часто заменятся на Элемент, будем для удобства считать их эквивалентами (элемент в данной статье не имеет ничего общего с методом конечных элементов).

Есть 2 различных способа решения задачи по FVM:
1) грани контрольного объема совпадают с гранями элемента
2) грани контрольного объема проходят через центры граней элементов(на которые разбита область).Искомые переменные хранятся в вершинах этих элементов.Вокруг каждой вершины строится контрольный объем. Для непрямоугольной сетки этот способ имеет еще 2 подвида.

Мы будем использовать способ 1) с контрольными объемами совпадающими с элементами на которые разбита область.

Некоторые плюсы FVM:

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

Метод FVM реализуем на примере уравнения теплопроводности:

Итак основные шаги при реализации FVM:

Перевод дифф уравнения в форму пригодную для FVM — интегрирование по контрольному объему
Составление дискретного аналога, выбор способа перевода производных и других подынтегральных выражений в дискретную форму
Получение уравнения для каждого из контрольных объемов, на которые разбита область.Составление системы линейных уравнений и ее решение.

Дискретизация по времени.

Немного теории или первый шаг в реализации FVM

FVM на стандартной прямоугольной сетке

На рисунке изображен Элемент С и его соседние элементы справа(E), слева(W), сверху(N) и снизу(S).У элемента С есть 4 грани обозначенные буквами e w n s.Именно эти 4 грани и составляют периметр элемента и по ним производится интегрирование.Для каждого элемента в результате получаем дискретный аналог исходного дифф уравнения.

Составим дискретный аналог для элемента С.Для начала нужно разобраться с интегралом (3).Интеграл это ведь по факту сумма.Поэтому мы и заменяем интеграл по всей поверхности элемента, на сумму по 4-м составляющим этой поверхности, тоесть 4 граням элемента.

Уравнение (7) и есть конечное уравнение для элемента С, из него мы на каждом шаге по времени получаем новое значение температуры (Tnew) в элементе С.

Граничные условия на прямоугольной сетке

Мы рассмотрим только 2 вида граничных условий.

Задана температура Tb на границе
Задан поток FluxB на границе, рассмотрим только случай когда FluxB=0, т.е. грань e будет теплоизолирована(Insulated)

Случай 2) самый простой, поскольку получается что грань e не потребуется при дискретизации(т.к. все коэффициенты Flux=0) и можно ее просто пропустить.

Теперь рассмотрим случай 1).Дискретизация грани e будет в целом похожа на ту что уже была описана.Будут только 2 изменения — вместо Te будет известное граничное значение Tb и вместо расстояния DXe будет DXe/2.В остальном можно рассматривать значение Te так, как будто это был бы обычный соседний узел E.Теперь подробнее распишем терм для граничного элемента С.

Пример численных расчетов на прямоугольной сетке

FVM в задачах со сложной геометрией

Здесь как раз проявляется преимущество FVM, где также, как и в методе конечных элементов, можно представлять область с круглыми границами через разбиение на треугольники или любые другие полигоны.Но FVM имеет еще 1 плюс — при переходе от треугольников к полигонам с большим числом сторон не требуется абсолютно ничего менять, конечно если код был написан для произвольного треугольника а не равностороннего.Более того, можно без изменения кода использовать смесь разных элементов — треугольники, полигоны, квадраты и тд.

Рассмотрим общий случай, когда вектор соединяющий центры 2-х элементов не совпадает с вектором нормали к общей грани этих элементов.Вычисление потока flux через грань теперь будет состоять из 2-х частей.В первой будет расcчитываться ортогональная составляющая а во второй так называемая «кросс-диффузия».

На картинке изображены 2 элемента, С — текущий рассматриваемый элемент и F — соседний элемент.Опишем дискретизацию для грани, разделяющей эти 2 элемента.Вектор соединяющий центры элементов — DCF.Вектор e — это единичный вектор по направлению DCF.Вектор Sf — направлен по нормали к грани, его длинна равна длине грани.

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

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

Примеры и проверка результатов

Описание структуры исходников

Гитхаб с исходниками лежит тут
Основная версия в папке heat2PolyV2.То что относится к вычислительной части лежит в heat2PolyV2\Src\FiniteVolume\.

Вначале файла Scene2.cs — параметры которые можно менять для отображения в разных цветовых схемах, масштаб, отображение mesh и т.д.Сами примеры хранятся в heat2PolyV2\bin\Debug\Demos\
Выгрузку из Матлаба сделать просто — нужно открыть pde toolbox, открыть m файл (либо создать самому с нуля), зайти в меню Mesh-Экспорт mesh, нажать ОК; перейти в основной Матлаб, в панельке появятся переменные — матрицы p e t, открыть файл savemymesh.m, выполнить его, появится файл p.out, перенести его в папку Demos.
В исходниках для выбора примера необходимо задать имя файла в строке param.file = «p»;(FormParam.cs).Далее необходимо применить граничные условия — для готовых примеров можно просто раскомментировать соответствующие блоки в MainSolver.cs:

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

Сравнительный анализ аналитического и численного решения задачи для уравнения теплопроводности

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

Для решения прикладных и математических задач в настоящее время существуют различные программы, с помощью которых можно довольно быстро вычислить сложные вычисления. К таким программам можно отнести MATCAD, MACHCAD, MATLAB, MAPLE и другие.

К прикладным задачам можно отнести условие, при котором допущено что дан тонкий однородный стержень 0 restart:

Задаем уравнение теплопроводности:

>eq:=diff(u(x,t),t)-a^2*(diff(u(x,t),x$2))=0;

где u(x,t) эта температура стержня в точке x в момент времени t, a связано с коэффициентами теплоемкости и теплопроводности.

Уравнение заносится в переменную eq для того, чтобы в дальнейшем можно было сослаться на эту переменную в своих вычислениях. При этом мы используем знак присвоения «:=». Если использовать двоеточие после команды – Maple выполнит ее, но не покажет результат. Это бывает удобно, когда много мелких команд.

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

> InitC:=u(x,0)=phi(x);

> phi:=x->A*x(l-x);

> BoundC:=u(0,t)=0,u(l,t)=0;

Получим решение методом разделения переменных. Для решения уравнения в частных производных имеются команда pdsolve и набор команд пакета PDEtools. Обращение к команде имеет следующий вид: pdsolve(PDE,FUN,OPT), где PDE уравнение, FUN неизвестная функция, а OPT дополнительные параметры.

Ответ может быть представлен специальной структурой PDESolStruc, которая содержит термин &where и состоит из двух полей: функционального представления и перечня обыкновенных дифференциальных уравнений, полученных в результате процедуры разделения переменных. Данная структура предназначена для указания возможных решений.

Решение при этом может быть получено с помощью команды build из пакета PDEtools. В качестве дополнительных параметров OPT могут использоваться: указание на возможную форму решения HINT=HI, параметр INTEGRATE, означающий автоматическое интегрирование получающейся системы обыкновенных дифференциальных уравнений, указание на поиск явного решения build. В качестве HI могут стоять знаки суммы и произведения или алгебраическое выражение. Уравнения могут быть заданы непосредственно в теле команды.

> res:=pdsolve(eq,HINT=X(x)*T(t));

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

> res1:=op(2,res);

> eq1:=op(1,res1[1]);

Хорошо известно, что физически осмысленные решения получаются при выборе отрицательной константы разделения. Теперь нужно сделать замену переменной _c1=-λ. Замену переменной можно провести в ручную, либо воспользовавшись командой subs, структура которой следующая: subs(x=a,EXPR), где x и а заменяемые переменные, а EXPR выражение, в котором происходит замена.

> eq1:=subs(_c[1]=-lambda,eq1);

Введем ограничение переменной λ, которая должна быть больше «0», для этого воспользуемся командой assume()

> assume(lambda>0);

Для аналитического и численного решения обыкновенных дифференциальных уравнений используется команда dsolve, причем во всех случаях используется единый формат команды: dsolve(ODE,VAR,OPT). Для решения задачи Коши в уравнения ODE нужно включить так же начальные условия, а для краевой задачи соответственно краевые условия.

> s1:=dsolve(eq1);

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

> X:=unapply(rhs(s1),x);

Найдем постоянные С1 и С2 из начальных условий:

> e1:=X(0)=0; e1:=_C2=0

> e2:=X(l)=0;

Составим систему из полученных уравнений:

> sistem:=; .

Запишем матрицу коэффициентов системы при помощи linalg[genmatrix]

> with(linalg):

> M:=genmatrix(sistem,[_C1,_C2]);

Вычислим определитель при помощи linarg[det](имя матрицы).

> Delta:=det(M);

Для получения всех решений уравнений, содержащих периодические функции, следует присвоить переменной _EnvAllSolutions значение true. В противном случае, по умолчанию функция solve выдаст только одно решение. При помощи команды solve решим полученное выражение относительно λ.

> _EnvAllSolutions:=true:solve(Delta,lambda);

Заменим Z1=k для удобочитаемости формулы, символ «%» означает предыдущую строку

> subs(_Z2=k,%);

> ev:=unapply(%,k);

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

> X:=’X’; X:=X

Ограничим значение переменной k и заменим переменную λ в уравнении eq1 на функцию ev(k):

> ur:=subs(lambda=ev(k),eq1);