Пример решения задачи о распределении тепла в пластине

Уравнение распространения тепла в пространственной области V имеет вид:

(1)

где u – температурное поле в области V, - плотность теплопроводящей среды, c – теплоемкость, - коэффициент теплопроводности. В случае однородной изотропной среды уравнение (1) можно переписать в виде

(2)

где - коэффициент температуропроводности, - оператор Лапласа.

На части границы области V могут быть заданы условия Дирихле (заданная величина температуры)

(3)

а на остальной части условия Неймана

(4)

где q – тепловой поток, проходящий внутрь области V, n – внешняя нормаль к ее границе.

В PDE ToolBox решается двумерная задача, которая может применяться, например, к пластинке постоянной толщины, теплоизолированной по обеим сторонам, а теплообмен осуществляется через боковые границы.

Задача состоит в том, чтобы найти статическое распределение тепла в области G, изображенной на рис. 2.

Круговые отверстия имеют одинаковый размер диаметра - 0.6 м. Правая и левая границы теплоизолированы. Коэффициент теплопроводности = 225.

Дифференциальное уравнение, описывающее распределение температуры, записывается в виде уравнения Пуассона:

(5)

Условия на левой и правой границах области G имеют вид, соответствующий теплоизоляции:

(6)

Здесь n - вектор нормали к границе.

На верхней и нижней границах области G поддерживается температура равная 580 градусов, поэтому там

(7a)

На окружностях поддерживается температура равная 450 градусов, поэтому там

(7b)

Рис. 2 Область G

Таким образом, задача о распределении тепла в области G полностью описывается уравнением (5) и граничными условиями.

Рассмотрим подробно процесс решения этой задачи в PDE ToolBox. Первый этап - конструирование области. Запустив среду PDE ToolBox (см. рис. 1), вы увидите, что по умолчанию область решения задачи должна лежать в прямоугольнике [-1.5 1.5; -1 1]. Необходимо установить другие границы изменения х и у (в рассматриваемом случае можно установить границы изменения х от -1 до 4, у - от -1 до 3). Для изменения границ отображаемой области по x и y необходимо выполнить команду меню Options ► Axes Limits (Настройки ► Пределы осей), после чего появится окно, в котором указаны границы изменения х и у (рис. 3). Изменим границы х (от -1 до 4) и у (от -1 до 3) (рис. 4).

Рис. 3 Окно Axes Limits Рис. 4 Новые значения границ

Для удобства можно изобразить линии сетки на координатной сетке. Для этого необходимо выполнить команду меню Options ► Grid (Настройки ► Сетка), после чего среда PDE ToolBox будет иметь вид, изображенный на рис. 5.

Рис. 5 Область G

Для конструирования области (в которой решается уравнение) предназначен пункт меню Draw (Рисовать). Для перехода в режим конструирования можно выполнить команду меню Draw ► Draw Mode (Рисовать ► Режим рисования) или просто выбрать один из геометрических примитивов и начать рисовать. Геометрические примитивы можно выбирать в пункте меню Draw или воспользоваться панелью рисования области. В PDE ToolBox есть следующие примитивы для рисования области:

Ø прямоугольник. Для выбора прямоугольника можно выполнить команду меню Draw ► Rectangle/square (Рисовать ► Прямоугольник/квадрат) или щелкнуть по кнопке на панели рисования области; если при прорисовке фигуры удерживать клавишу Ctrl, то вместо прямоугольника получится квадрат;

Ø центрированный прямоугольник. Для выбора этой фигуры можно выполнить команду меню Draw ► Rectangle/square (centered) (Рисовать ► Прямоугольник/квадрат (Центрированный)) или щелкнуть по кнопке на панели рисования области; если удерживать клавишу Ctrl, то вместо прямоугольника получится квадрат. Отличие этого прямоугольника от обычного состоит в том, что при прорисовке прямоугольник растягивается не только в сторону правого нижнего угла, но и левого верхнего;

Ø эллипс. Для его выбора можно выполнить команду Draw ► Ellipse/circle (Рисовать ► Эллипс/круг) или выбрать инструмент на панели рисования области), при нажатой клавише Ctrl вместо эллипса получится круг; кроме обычного эллипса есть центрированный (команда меню Draw ► Ellipse/circle (centered) (Рисовать ► Эллипс/круг (Центрированный)) или кнопка ;

Ø многоугольник (команда меню Draw ► Polygon (Рисование ► Многоугольник) иди кнопка ).

Построенные объекты можно поворачивать с помощью команды Draw ► Rotate… (Рисование ► Повернуть…).

Для решения рассматриваемой задачи необходимо нарисовать прямоугольник, а в нем два круга. Начав с рисования прямоугольника, следует нарисовать его по возможности точнее. Для уточнения координат прямоугольника необходимо дважды щелкнуть по прямоугольнику, что приведет к появлению диалогового окна Object Dialog (Вызов диалогового окна для изменения свойств объекта) - рис. 6.

Рис. 6 Диалоговое окно Object Dialog

Параметры Left (Левый) и Bottom (Нижний) определяют координаты левого нижнего угла прямоугольника, там следует ввести нули (левый нижний угол в данном случае совпадает с началом координат). Параметры Width (Ширина) и Height (Высота) определяют ширину и высоту прямоугольника (для решаемой задачи эти величины равны 3 и 2 соответственно). Параметр Name (Имя) определяет имя объекта. Все прямоугольники по умолчанию начинаются с «R» (от английского слова Rectangle), имя эллипса начинается с буквы «Е», круга - с буквы «С», многоугольника - с буквы «P», квадрата - с букв «SQ».

Следующим этапом будет рисование двух кругов, первый - с центром в точке (1, 1) и радиусом 0.3, второй - с центром в точке (2.2, 1) и радиусом 0.3. Для построения круга используется кнопка при нажатой клавише Ctrl. Затем следует уточнить его координаты: X-Center, Y-Center - координаты центра круга, Radius - радиус круга (рис. 7).

Рис. 7 Диалоговое окно Object Dialog для круга

Рис. 8 Диалоговое окно Object Dialog для круга

Теперь надо определить способ построения области G при помощи операций над множествами точек построенных элементарных объектов. Для этого в PDE ToolBox имеются операции:

- - операция вычитания объектов, имеет наивысший приоритет;

+ - операция объединения объектов;

* - операция пересечения объектов.

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

По умолчанию созданные объекты соединяются с помощью знака «плюс». В рассматриваемом случае необходимо из прямоугольника вычесть два круга, поэтому в области формулы необходимо исправить формулу на R1-C1-C2 (рис. 8).

Рис. 9 Выбор типа решаемой задачи

Рис. 10 Окно PDE Specification - определение типа и коэффициентов дифференциального уравнения в частных производных

На втором этапе решения задачи необходимо ввести уравнение в частных производных. В PDE ToolBox есть возможность задавать тип решаемой задачи. Для этого можно воспользоваться командой Options ► Application (Настройки ► Приложение) - рис. 9 или раскрывающимся списком . Для ввода уравнений служит пункт меню PDE. Первый подпункт этого меню PDE Mode служит для перехода в режим ввода уравнений. Второй подпункт служит для показа подобластей, входящих в область решения. В решаемой задаче есть всего одна подобласть. Третий пункт PDE Specification или кнопка позволяет выбрать тип уравнений и ввести его коэффициенты (рис. 10).

При решении данной задачи следует выбрать в качестве типа задачи Heat Transfer (Задача о распределении тепла) или просто ввести коэффициенты эллиптического уравнения в окне PDE Specification. Установив тип решаемой задачи Heat Transfer, параметры уравнения в частных производных следует определить в окне PDE Specification (рис. 11).

Заметим здесь, что для скалярной функции Т(х, у, z) ee градиентом является вектор grad(T), компоненты которого равны , а дивергенция от градиента div(grad(T)) вычисляется по формуле: (в плоской задаче координата z отсутствует).

Рис. 11 Окно PDE Specification - определение коэффициентов задачи о распределении тепла

Третий этап решения задачи посвящен определению граничных условий. Для ввода граничных условий предназначен пункт меню Boundary (Граничные условия). Для перехода в режим ввода граничных условий можно выполнить команду меню Boundary ► Boundary Mode (Граничные условия ► Режим граничных условий) или воспользоваться кнопкой или же комбинацией клавиш Ctrl+B.

После этого, чтобы определить граничные условия, необходимо щелчком мыши выделить границу и выполнить команду меню Boundary ► Specify Boundary Condition (Граничные условия ► Установить граничные условия) или осуществить двойной щелчок мышью по выделенной границе. (Если необходимо объединить несколько участков границы, то в этом случае щелчок по каждому из участков осуществляется при нажатой клавише Shift).

На верхнем и нижнем участках границы области определены условия Дирихле вида

hT = r ,

и в окне Boundary Condition следует ввести h = 1, r = 580 (рис. 12).

Рис. 12 Ввод условий для верхней и нижней границ области G

Выделив одновременно контуры отверстий, установить температуру 450 градусов (рис. 13).

Рис. 13 Ввод условий на границах отверстий области G

Рис. 14 Ввод условий Неймана для левой и правой границы области G

Осталось ввести условия Неймана вида (6) (или, что то же самое, n • k • grad(T) = 0) на левой и правой границе области G (рис. 14).

Ввод условий задачи завершен.

На четвертом этапе выполнить триангуляцию области - покрыть область сеткой, состоящей из треугольников. Это и есть конечные элементы, на которые разбивается область. Режимом триангуляции управляет пункт меню Mesh (Сетка). Для перехода в режим триангуляции служит команда Mesh ► Mesh Mode (Сетка ► Режим сетки). Переход в режим триангуляции сразу разбивает на крупные треугольники (рис. 15). Этого же можно добиться с помощью команды начала триангуляции - пункт меню Mesh ► Initialize Mesh (Сетка ► Определить (Инициализировать) сетку), комбинация клавиш Ctrl+I или кнопка .

Обычно начальной триангуляции недостаточно для получения решения с удовлетворительной точностью. Для разбиения области на более мелкие треугольники служит команда меню Mesh ► Refine Mesh (Сетка ► Улучшить сетку) или кнопка . Для решаемой задачи критерием достаточной величины треугольников может служить вид отверстий. Двойное применение команды Mesh ► Refine Mesh, дает сетку, необходимую для получения решения с достаточной точностью (рис. 16). Следует учесть, что из-за слишком мелкой сетки время решения задачи значительно увеличивается. Возврат к начальной триангуляции осуществляется с помощью команды меню Mesh ► Initialize Mesh или кнопки .

Рис. 15 Начальная триангуляция области G

На пятом этапе для решения поставленной задачи необходимо выполнить команду меню Solve ► Solve PDE (Решение ► Решить PDE) или щелкнуть по кнопке . После этого найденное решение автоматически изображается на экране в виде цветного контурного графика (рис. 17).

Рис. 16 Уменьшенная сетка в области решения G

Шестой этап решения задачи состоит в обработке графического представления результатов.

Контурный график не всегда достаточно хорошо позволяет проанализировать полученное решение. Пользователь может изменить тип построенного графика, для чего можно выполнить команду меню Plot ► Parameters (График ► Параметры) или щелкнуть по кнопке . В появившемся диалоговом окне (рис. 18) можно установить необходимые параметры графика решения, для построения трехмерного графика достаточно установить флажок Height (3-D plot). На рис. 19 изображен трехмерный график найденного решения задачи.

Рис. 17 Контурный график решения задачи

Рис. 18 Диалоговое окно для установки параметров графика

Рис. 19 Трехмерный график решения задачи

При решении задачи были рассмотрены основные этапы решения дифференциального уравнения в частных производных в графической среде PDE ToolBox.

Теперь рассмотрим, как в PDE ToolBox хранится область решения, сетка и полученное решение. Для экспорта области решения в меню Draw есть команда Export Geometry Description, Set Formula, Labels… (Экспорт описания геометрии области, формулы, меток), с помощью которой (рис. 20) можно экспортировать геометрию области и формулу связи между объектами в MATLAB.

Рис. 20 Диалоговое окно Export, функция Export Geometry Description, Set Formula, Labels...

Как видно из рис. 20, можно экспортировать три объекта gd, sf и ns. Матрица геометрии области хранится в gd, формула связи хранится между объектами в переменной sf. Рассмотрим структуру матрицы gd на примере области G (рис. 8.29).

gd =

3.0000 1.0000 1.0000

4.0000 1.0000 2.2000

0 1.0000 1.0000

3.0000 0.3000 0.3000

3.0000 0 0

0 0 0

0 0 0

0 0 0

2.0000 0 0

2.0000 0 0

Число столбцов матрицы совпадает с количеством объектов, входящих в область решения, каждый столбец отвечает за отдельный объект. Рассмотрим структуру каждого столбца. Первый элемент столбца содержит номер объекта: 1 - круг, 2 - многоугольник, 3 - прямоугольник, 4 - эллипс.

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

Для круга второй и третий элементы содержат координаты центра круга, четвертый - радиус.

Рис. 21 Сетка с пронумерованными узлами и конечными элементами

Эллипс отличается от круга тем, что четвертый и пятый элементы столбца содержат величины полуосей эллипса, а шестой - угол поворота в радианах. В текстовой переменной sf хранится формула:

sf =

R1-C1-C2

Для вывода на экран номеров узлов сетки служит команда Mesh ► Show Node Labels (Сетка ► Показать узлы сетки), а для вывода номеров треугольных конечных элементов - команда Mesh ► Show Triangle Labels (Сетка ► Показать метки треугольных элементов). После нумерации узлов и конечных элементов сетка будет иметь вид, такой как на рис. 21 (Показан фрагмент сетки, выделенный при помощи кнопки ).

Теперь рассмотрим, как хранится сетка. Для этого экспортируем ее в MATLAB с помощью пункта меню Mesh ► Export Mesh (Сетка ► Экспорт сетки) - рис. 22.

Рис. 22 Окно Export, функция Export Mesh

Матрица p содержит координаты всех точек, каждый столбец содержит координаты одной точки (координаты х и у); для области G она имеет вид (фрагмент в редакторе матриц):

Матрица t содержит 4 строки:

В каждом столбце матрицы: первые три элемента содержат номера узлов, образующих треугольник (конечный элемент), четвертый - номер подобласти, которой он принадлежит. Число столбцов матрицы t равно количеству конечных элементов, на которые разделена область, Для области G матрица t будет такой (4-я строка постоянно равна 1, так как область решения состоит из одной подобласти):

В матрице е находится информация об узлах, формирующих границу подобласти. В нашем примере матрица е имеет вид:

Решение представляет собой массив значений u: его элементы соответствуют координатам узлов, определенных матрицей р. Количество элементов массива решений u равно количеству столбцов матрицы р. Для экспорта решения в MATLAB служит команда меню Solve ► Export Solution (Решение ► Экспорт решения). Сопоставляя значения матриц р, t и е, можно анализировать полученное решение в различных точках. Чем мельче сетка, тем более точное решение будет получено. Но, с другой стороны, мелкая сетка требует больше времени решения задачи.

Варианты заданий для самостоятельного решения

Решить статическую задачу распределения тепла в одной из областей в соответствии с номером варианта (рис. 23). Выделенные жирными линиями границы области соответствуют заданию условий Дирихле с указанным значением температуры. На остальных частях границы заданы условия теплоизоляции.

Таблица 1

Материал: сталь (плотность = 7800 кг/м³, теплоемкость c = 465 Дж/(кг×град), коэффициент теплопроводности = 47 Вт/(м×град)). Все данные приведены в системе СИ.

Пластина A:

Пластина B:

Пластина C:

Рис. 23

Поиск по сайту: