Изучаем самостоятельно Toolbox PDE application MATLAB для моделирования динамических процессов с распределенными параметрами

В работе рассмотрены следующие вопросы:

1. Описание техники решения двумерных уравнений с частными производными с использованием MATLAB приложения PDE app.
2. Разбор примера решения задачи нагрева пластины внутренним источником тепла с использованием PDE app.
3. Пример для самостоятельной работы – решение задачи нагрева пластины внешним потоком тепла с использованием MATLAB приложения PDE app.
4. Программные средства решения уравнений с частными производными (на примере уравнения теплопроводности).

Авторы: Крушель Е.Г., Панфилов А.Э., Харитонов И.М.

Введение

В работе приведено описание принципов решения задач с частными производными 2-го порядка в среде одного из приложений системы MATLAB – приложения Partial Differential Equations Toolbox (далее используется принятое в MATLAB сокращенное наименование PDE application и PDE app).

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

Сложность учета пространственной распределенности моделируемых явлений приводит к вынужденному использованию упрощений в ущерб адекватности. В частности, популярным приёмом составления математической модели объекта является следующий. Модель – материальная точка, не имеющая пространственных размеров, но обладающая физическими свойствами (массой, инерционностью, колебательностью и др). Объект в целом имитируется множеством таких точек, реакции которых на внешние и внутренние воздействия одинаковы. Класс таких моделей – модели с сосредоточенными параметрами. Для их описания привлекаются обыкновенные дифференциальные уравнения, и во многих практически важных случаях такой подход позволяет получить работоспособные технические решения.

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

1. Вращающаяся печь цементной промышленности, длина 150м и более. В верхней зоне – процессы высушивания сырья, в средней – разогрев, в нижней – спекание. Эти процессы не сходны.
2. Магистральные каналы ирригационных систем, протяженность – сотни километров. Волновые процессы в канале различны по его длине.
3. Процессы распространения выбросов вредных веществ из источников загрязнения атмосферы с учетом метеорологических условий.
4. Длинные линии электропередачи. Волновые процессы в них также различны по длине.

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

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

Учебной целью работы является предоставление студентам и выпускникам технического вуза помощи в освоении системы компьютерной математики (СКМ) MATLAB для моделирования систем с распределенными параметрами. Представляемый материал может быть встроен в технических вузах в учебные программы дисциплин «Задачи математической физики» и «Системы компьютерной математики» в качестве компьютерной поддержки. После освоения раздела MATLAB, посвященного уравнениям с частными производными (ЧП), выпускник вуза сможет самостоятельно решать задачи моделирования распределенных процессов и систем управления объектами с распределенными параметрами, профессионально используя совершенные математические средства, которые ранее были доступны только квалифицированным математикам.

Для достижения этой цели решены задачи:

1. Освоение решателя Pdesolve() системы Mathcad для решения уравнений и систем уравнений с ЧП.
2. Освоение приложения системы MATLAB PDE toolbox для решения уравнений и систем уравнений с ЧП с постоянными параметрами.
3. Разработка MATLAB программы для решения уравнений с ЧП, параметры которых зависят от пространственных координат и времени.

Форма представления материала – последовательное решение задач с нарастающей сложностью. Выбранный уровень подробности изложения позволяет читателю самостоятельно изучать представленный материал, что делает его доступным для студентов заочной и/или дистанционной формы обучения. Краткие сведения о классификации уравнений с ЧП и о способах их приближенного решения приведены в Приложениях 2 и 3.

Далее используется сокращенное наименование средств для решения уравнений с частными производными, принятое в MATLAB: PDE application или PDE app.

Рассмотренная ниже область применения PDE application MATLAB:

1. Для изучения MATLAB техники решения уравнений с частными производными 2-го порядка (т.е. на плоскости) – как отдельных (скалярных) уравнений, так и систем из двух уравнений с частными производными 2-го порядка. В частности, предусмотрены следующие возможности:
   a. Решение скалярного уравнения с частными производными общего типа; могут быть решены уравнения эллиптического, параболического и гиперболического типов с граничными условиями первого и второго рода (т.е. с условиями Дирихле и Неймана);
   b. Решение систем из двух уравнений общего типа.
2. Для ознакомления с возможностями решения специальных инженерных задач с использованием уравнений с ЧП, в частности:
   a. Решение уравнений структурной механики – распределение напряжений на плоскости;
   b. Решение уравнений структурной механики – распределение деформаций на плоскости;
   c. Решение задач электростатики;
   d. Решение задач магнитостатики;
   e. Изучение электромагнитных явлений в электродвигателях постоянного тока;
   f. Изучение явлений проводимости в средах постоянного тока;
   g. Решение уравнений теплопередачи;
   h. Решение уравнений диффузии.

Позиции п.2 относятся к разделам физики, требуют специальных знаний. Поэтому изложение ограничивается описанием общей техники использования PDE app, освоение которой достаточно для самостоятельного решения любых задач п. 1, 2.

Описание последовательности действий основано на текстах, имеющихся в файлах помощи системы MATLAB. Поясняющие примеры и способы действий разработаны самостоятельно.

 

1. Описание техники решения двумерных уравнений с частными производными с использованием PDE app

Сведения, помещенные ниже, описывают последовательность шагов, необходимых для решения уравнения с частными производными (далее используется сокращение ЧП) – в частности, описание назначения опций меню, предусмотренных для PDE app.

1.1. Пошаговое выполнение расчетов в PDE app

Шаг 1. Черчение. При вызове PDE app из главного меню MATLAB, вкладка Apps, приложение PDE app находится в режиме черчения (Draw Mode), в котором пользователь может использовать 5 типов твердотельных объектов для вычерчивания модели конструктивной твердотельной геометрии, КТГ (соответствующий англоязычный термин Constructive Solid Geometry (CSG) model).

Также можно редактировать формулу, соединяющую элементы геометрии области, на которой нужно найти решение уравнения с ЧП. Твердотельные объекты выбираются либо с помощью пяти кнопок (расположенных слева), либо из меню Draw (рис. 1, 2).

Рис. 1 – Пять кнопок твердотельных объектов, выбираемых иконками

Рис. 2 – Выбор твердотельных объектов из меню Draw

Справа от кнопок твердотельных объектов находятся кнопки, с помощью которых пользователь получает доступ ко всем функциям, необходимым для решения уравнения с ЧП: определить граничные (краевые) условия, сконструировать сетку, которая будет использована для приближенного речения уравнения с ЧП (в PDE app использована триангуляционная сетка, разработанная членом-корреспондентом Академии наук СССР Б.Н. Делоне), а также для визуализации решения.

PDE app используется как графическое средство для вычерчивания областей с 2D-геометрией на координатной плоскости (координаты x,y). Эта операция трактуется как чертеж области, на которой будет осуществляться поиск решения. Для этого используется 4 базовых твердотельных объекта, сетка и операция «привязать к сетке» (snap-to-grid).

Работа с PDE app начинается в режиме черчения (Draw). Пользователю предоставляется возможность выбора типа объекта, комбинирования объектов друг с другом и определения алгебры для соотношений, с помощью которых соединяются комбинируемые объекты (рис. 3).

Рис. 3 – Примеры объектов и формула для их комбинирования

Для того, чтобы выделять суб-области, нужно щелкнуть по её обозначению (выделится суб-область, после чего её можно перемещать, рис. 4).

Рис. 4 – Положение курсора для выделения суб-области

Можно произвести запись чертежа в файл параметров модели. Если пользователь намерен продолжать расчеты с той же моделью в следующей MATLAB-сессии с PDE app, достаточно вызвать соответствующий файл (например, из командной строки MATLAB или из меню File приложения PDE app). Запись в файл может производиться на всех стадиях работы (при этом все ранее сделанные операции будут сохранены, так что их повторение не потребуется).

Шаг 2. Формирование граничных условий. При нажатии иконки  (рис. 5) произойдёт выделение внешних границ области с граничными условиями, предусмотренными по умолчанию (рис. 6).

Рис. 5 – Расположение иконки формирования граничных условий

Рис. 6 – Выделение внешних границ после нажатия иконки 

Выделяются только внешние границы. Если нужно изменить границы встроенных суб-областей, нужно указать их в формуле комбинирования объектов (см. рис. 3).

Если внешние границы получились неверно (т.е. не совпадают с границами, которые соответствуют решаемой задаче), нужно возвратиться в режим черчения (в этом режиме можно добавлять, удалять или изменять любой твердотельный объект, а также изменять формулу комбинирования объектов, рис. 3). Замечание: формула может корректироваться только в режиме черчения.

Например, если нужно изменить геометрию части объектов, после перехода в режим Draw потребуется указать перечень этих объектов. После перехода в режим Boundary Mode (иконка Boundary) будут показаны границы только тех объектов, которые указаны в формуле комбинирования объектов (рис. 7).

Рис. 7 – Корректировка геометрии выбранной части объектов

Если потребуется, можно устранить часть или все границы (опции Remove Subdomain Border или Remove All Subdomain Borders из меню Boundary).

Например, после удаления всех границ суб-областей (Remove All Subdomain Borders) остаются только внешние границы. Можно пронумеровать каждый элемент границы (опция Show Edge Labels, показать метки границы). Щелчок мыши по любому отрезку границы приведёт к его выделению. Двойной щелчок приведёт к выводу всплывающего окна для ввода условий соответствующей краевой задачи (рис. 8, это же можно сделать с помощью опции Specify Boundary Conditions (определение краевых условий) из меню Boundary).

Рис. 8 – Всплывающее окно для ввода граничных условий на выбранном отрезке границы

Шаг 3. Инициализация триангуляционной сетки. Инициализация производится нажатием иконки  или путем использования опции Initialize Mesh меню Mesh (сетка). Обычно параметры сетки, установленные по умолчанию, генерируют сетку, адекватную задаче. Но если по каким-то причинам она не подходит пользователю, её можно изменить опцией Parameters меню Mesh (рис. 9).

Рис. 9 – Генерация триангуляционной сетки, принятой по умолчанию

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

Рис. 10 – Результаты двух нажатий иконки Refine

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

Шаг 4. Выбор типа уравнения с ЧП. Тип уравнения можно непосредственно выбрать в списке панели элементов окна (рис.11) или в диалоговом окне PDE Specification (спецификация уравнения с ЧП), вызываемого опцией PDE Specification из меню PDE. Замечание: этот шаг можно делать на любой стадии решения задачи (он не зависит от предшествующих операций). Если коэффициенты уравнения с ЧП зависят от выбора материала, то их нужно вводить в режиме PDE путем двойного щелчка по наименованию соответствующей суб-области.

Рис. 11 – Выбор типа уравнения из диалогового окна

Шаг 5. Решение уравнения с ЧП. Для этого нужно нажать иконку  (рис. 12) или выбрать опцию Solve PDE (решить уравнение с ЧП) из меню Solve.

Рис. 12 – Вызов решателя уравнения с ЧП. Здесь – для решения уравнения теплопроводности (Heat Transfer)

Шаг 6. Настройка графического вывода. Если пользователь не намерен запрашивать автоматический вывод графика решения или если он хочет изменить способ представления решения, нужно вызвать диалоговое окно Parameters из меню Plot (графики) до того, как будет запущена процедура решения (нажатием иконки ), рис. 13.

Рис. 13 – Диалоговое окно установки параметров графического вывода

1.2. Возможности, предоставляемые приложением PDE app
В дополнение к описанным шагам (п. 1.1) имеются следующие альтернативы изменений:

• Экспорт сетки или решения в рабочую область MATLAB для последующего анализа.
• Визуализация иных свойств решения.
• Изменение вида и параметров уравнения с ЧП и пересчет решения.
• Изменение сетки и пересчёт решения. Можно либо инициализировать секту (Initialize mesh из меню Mesh), либо уменьшить размеры треугольников сетки (Refine Mesh из того же меню). Из меню Mesh можно также подвигать сетку (Jiggle mesh) или отменить внесенные изменения (Undo Mesh Change).
• Изменить граничные условия (нужно использовать иконку или опцию Boundary Mode из меню Boundary).
• Изменить геометрию модели. Для этого нужно возвратиться в режим черчения (Draw mode) или кликнуть иконку нужного твердотельного объекта. После возврата в режим черчения можно добавлять, изменять или удалять твердотельные объекты, а также изменять формулу комбинирования объектов.

Последовательность шагов, описанная в п. 1.1, является гибкой: предусмотрены комбинации клавиш, позволяющие «проскакивать» некоторые шаги (в этом случае приложение PDE app вставит пропущенные шаги самостоятельно, по умолчанию. В этом проявляется общий стиль системы MATLAB, всемерно содействующий пользователю, даже не имеющим большого опыта работы с информационными технологиями). В частности:

1. Если пользователь до данной стадии работы ещё не определил геометрию модели и вышел из режима черчения (Draw Mode) с пустой моделью, приложение PDE app автоматически создаст L-образную геометрию с граничными условиями по умолчанию и будет выполнять последующие шаги с этой моделью.
2. Если пользователь установил режим черчения (Draw Mode) и нажал иконку  для инициализации сетки, то приложение PDE app вначале автоматически произведёт декомпозицию геометрии согласно текущей формуле комбинирования объектов и установит граничные условия по умолчанию на внешних границах. После этого будет сгенерирована первичная сетка.
3. Если пользователь нажал на иконку  для уточнения размеров треугольников сетки прежде, чем была инициализирована первичная сетка, приложение PDE app автоматически сгенерирует первичную сетку, а также произведёт декомпозицию геометрии, если активен режим черчения (Draw Mode).
4. Если пользователь нажмёт на иконку  для того, чтобы решить уравнение с ЧП, но сетка не создана, приложение PDE app автоматически создаст сетку перед тем, как решать уравнение с ЧП.
5. Если пользователь определит тип графического вывода и запросит вывод графика, приложение PDE автоматически проверит, доступно ли решение рассматриваемого уравнения. Если окажется, что решение недоступно, приложение PDE app автоматически решит уравнение и выведет график с опциями, выбранными пользователем.
   
   
6. Если пользователь не определил тип уравнения с ЧП, то приложение PDE app автоматически решит уравнение, предусмотренное по умолчанию (уравнение Пуассона):
   
   (уравнение соответствует встроенному в приложение PDE app уравнению эллиптического типа с коэффициентами c = 1, a = 0, f = 10).
   
   
7. Для разных режимов приложение PDE app автоматически устанавливает параметры, предусмотренные для этого уравнения по умолчанию.

 

2. Разбор примера. Решение задачи нагрева пластины внутренним источником тепла с использованием PDE app

Постановка задачи:

Решить базовое уравнение теплопроводности с источником тепла. Уравнение имеет вид:

или, расшифровывая оператор  , можно записать в привычной форме:

, здесь a = 1

Замечание: ниже физические размерности величин и параметров не указываются. Предполагается, что они записаны в совместимой системе единиц измерения, например:

[Размерность плотности материала]: кг/м³;
[Размерность температуры]: °C;
[Размерность времени]: час;
[Размерность координат чертежа]; м;
[Размерность теплового потока]: ккал/час;
[Размерность коэффициента теплопроводности]: ккал/(м×час×°C);
[Размерность удельной теплоемкости]: ккал/(кг×°C);
[Размерность коэффициента теплопередачи]: ккал/(м²×час×°C).

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

Параметры суб-областей:

1) Квадратная (внешняя) суб-область: коэффициент теплопроводности равен 10, плотность равна 2. Координаты вершин квадрата: (0, 0), (3, 0), (3, 3), (0, 3).

2) Ромбовидная (внутренняя) суб-область: коэффициент теплопроводности равен 2, плотность равна 1. Угол поворота ромба относительно квадратной суб-области (45°), координаты вершин (1.5, 0.5), (2.5, 1.5), (1.5, 2.5), (0.5, 1.5).

На всех внешних границах температура равна 0 (что соответствует граничным условиям, предусмотренным для приложения PDE app по умолчанию – следовательно, при решении задачи их вводить не нужно).

Ниже последовательно описаны этапы решения.

1. Из основного меню MATLAB вкладка Apps вызываем приложение PDE app (рис. 14).

Рис. 14 – Размещение иконки PDE app в составе иконок MATLAB-приложений

2. В появившемся рабочем окне приложения PDE app указываем тип решаемой задачи - Heat Transfer, теплопередача (рис. 15).

Рис. 15 – Установка типа уравнения с ЧП

3. Переходим в режим черчения (позиция Draw Mode меню Draw, рис. 16).

Рис. 16 – Выбор режима черчения (Draw Mode)

4. Выбираем позицию выпадающего меню Axes Limits из горизонтального меню Options и вводим границы координат (X-axes limits для оси абсцисс и Y-axes limits оси ординат, рис. 17). Границы вводим с избытком по сравнению с заданными координатами внешней области – например, (–0.5 3.5) для границ осей x и y. Флажки Auto в окне Axes Limits должны быть отключены (рис. 17).

Рис. 17 – Установка пределов шкал осей абсциссы и ординаты

 5. Нажимаем Apply (применить), Close (закрыть) в окошке, получим чертеж с заданными размерами осей (рис. 18).

Рис. 18 – Результат установки пределов шкал осей

6. Несмотря на то, что мы вводили одинаковые пределы для осей, область получилась неквадратной (приложение PDE app по умолчанию вводит разные масштабы для осей, чтобы чертеж получался «красивым», с соблюдение пропорций «золотого сечения»). Нам же, для наглядности, нужно видеть не искаженные области (квадратную и ромбовидную), поэтому включаем опцию выпадающего меню Axes equal (равные оси) из меню Options (рис. 19).

Рис. 19 – Установка равного масштаба для осей чертежа

Мы видим, что границы оси абсцисс автоматически изменились; также автоматически включился флажок Auto в окошке ввода границ осей, п.4. Не нужно обращать на это внимания: при вводе координат вершин областей ось абсцисс будет автоматически подстраиваться.

7. Выбираем иконку (вычерчивание прямоугольника-квадрата). Зажав клавишу Ctrl, чтобы вычертился квадрат, а не прямоугольник, начинаем рисовать квадрат, начиная с левого нижнего угла, т.е. «на глаз» с точки с абсциссой 0, ординатой 0. Не отпуская клавишу Ctrl, протягиваем вверх указатель мышки примерно до ординаты 3 (квадрат вырисовывается автоматически). В итоге получим затененный квадрат с приблизительно правильными коодинатами вершин и обозначением суб-области SQ1 (square1, квадрат1), рис. 20. Благодаря опции Axes equal (равное оси) квадрат визуально не искажен.

Рис. 20 – Результат вычерчивания квадрата

8. Теперь откорректируем координаты вершин квадрата, которые по постановке задачи равны (0,0), (3,0), (3,3), (0,3) (в п. 7 мы их устанавливали «на глаз»). Устанавливаем курсор на имя области (SQ1) и двойным щелчком левой кнопки мыши вызываем диалоговое окно ввода параметров квадрата. По постановке задачи квадрат имеет ширину сторон 3, слева 0, снизу 0, тогда получится квадрат с требуемыми координатами (0,0), (3,0), (3,3), (0,3). Вводим эти координаты в диалоговое окно (рис. 21), затем нажимаем OK.

Рис. 21 – Уточнение координат вершин квадрата

9. Переходим к вычерчиванию ромбовидной суб-области. Технология вычерчивания:

а) Вычерчиваем прямоугольник примерно посередине области SQ1 с помощью кнопки . Можно нажимать только левую кнопку мыши, без Ctrl, т.к. координаты углов ромба вводятся по очереди, не предполагая равнества его сторон (рис. 22).

Рис. 22 – Результат вычерчивания внутренней суб-области

Мы видим, что появилось новое обозначение суб-области - R1 (rectangle1, прямоугольник1).

б) Из меню Draw позицией Rotate… вызываем диалоговое окно угла вращения (Rotate) и вводим угол поворота 45 (согласно постановке задачи). Суб-область R1 повернется на 45° относительно суб-области SQ1 (рис. 23).

Рис. 23 – Результат вращения прямоугольника

в) Мы ввели координаты вершин ромба неточно. Откорректируем их: установим курсор на обозначение суб-области R1 и двойным щелчком левой кнопки мыши вызовем диалоговое окно последовательного ввода координат вершин. Сначала там будут показаны координаты вершин, которые мы ввели «на глаз» (рис. 24).

Рис. 24 – Начальные значения координат вершин ромба

Красным кружком на рис. 24 показано выпадающее меню координат вершин.

Координаты нужно вводить по направлению обхода ромба, начиная с любой вершины. Важно не допускать нарушения последовательности – иначе приложение PDE app выведет предупреждение о недопустимости пересечения сторон многоугольника (в данном случае – ромба). Не нажимая OK, последовательно вызываем координаты из выпадающего меню и вводим абсциссу и ординату. Например, начнем с вершины (0.5, 1.5), рис. 25.

Рис. 25 – Ввод координат первой вершины ромба

Не нажимая OK, выбираем из выпадающего меню координаты второй вершины по порядку обхода и корректируем ее (рис. 26).

Рис. 26 – Последовательный ввод координат вершин ромба

Дальше аналогично вводим координаты вершин (2.5, 1.5) и (1.5, 0.5). После этого нажимаем ОК и получаем правильно введенные координаты ромбовидной суб-области внутри квадратной внешней суб-области.

Мы видим также формулу комбинирования суб-областей (Set Formula), в которой записаны обозначения всех суб-областей чертежа – в данном случае SQ1+R1.

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

11. Переходим к вводу параметров уравнения теплопроводности. Сначала проверяем, не произошло ли случайное переключение типа уравнения. Если произошло – исправляем: должно быть выбрано Heat Transfer (рис. 27).

Рис. 27 – Проверка правильности установки типа уравнения с ЧП

12. В последовательности шагов по решению уравнения с ЧП в среде приложения PDE app указано, что после завершения разработки геометрии области нужно определить граничные условия (на всей внешней границе компоновки суб-областей). В рассматриваемом примере этот шаг можно пропустить: по постановке задачи температура на границах предполагается нулевой, а это в приложении PDE app предусмотрено по умолчанию, и вводить те же граничные условия ещё раз не нужно.

13. Переходим к вводу параметров уравнения теплопроводности. По постановке задачи эти параметры различны для суб-областей. Переходим в меню PDE и выбираем опцию PDE mode (режим ввода параметров уравнения с ЧП). В этом режиме суб-области выделяются незаметными белыми линиями, а четко очерченные границы, которые были установлены в режиме черчения (Draw Mode), исчезают (рис. 28).

Рис. 28 – Вид области решения в режиме PDE mode

14. Ввод параметров уравнения с ЧП в режиме PDE mode лучше начинать с внутренних суб-областей (в нашем случае – с ввода параметров ромбовидной области). Однократным щелчком левой кнопки мыши в режиме PDE mode выделяем ромбовидную область (рис. 29).

Рис. 29 ­– Выделение внутренней суб-области в режиме PDE mode

Из меню PDE выбираем опцию PDE specification… (спецификация параметров уравнения с ЧП). В появившемся диалоговом окне устанавливаем тип уравнения Parabolic (уравнение параболического типа, к которому относится уравнение теплопроводности). Вводим параметры уравнения согласно постановке задачи (рис. 30):

• внешняя температура (T_ext) = 0;
• коэффициент конвективного переноса тепла (h) = 0 (а можно любой, т.к. конвекция имеет место только если разность между внешней и граничной температурами не нулевая, но в нашем примере эта разность равна 0, поэтому соответствующее слагаемое в уравнении, записанном в верхней строке окна, останется нулевым при любом значении h);
• мощность внутреннего источника тепла (Q) = 4.

Рис. 30 – Ввод параметров уравнения с ЧП для внутренней суб-области

15. Теперь в режиме PDE mode вводим параметры уравнения с ЧП для квадратной (внешней) суб-области. Щелчком левой кнопки мыши в любой точке квадратной суб-области, но не внутри ромбовидной суб-области, выделим ее (рис. 31).

 
Рис. 31 – Вид чертежа в режиме PDE mode для ввода параметров внешней суб-области

При вводе параметров внешней (квадратной) суб-области параметры внутренней суб-области не изменятся (несмотря на то, что она выделена). Вводим параметры уравнения с ЧП для квадратной суб-области согласно постановке задачи, аналогично описанному в п.14 выше (рис. 32).

Рис. 32 – Ввод параметров уравнения с ЧП для внешней суб-области

Внимание: согласно постановке задачи, квадратная суб-область SQ1 не имеет внутренних источников тепла. В окошке для Q (источник тепла) вводим 0.

16. Из меню Mesh вызываем опцию Initialize Mesh (инициализация сетки) для покрытия области, на которой ищем решение, триангуляционной сеткой Б.Н. Делоне (рис. 33).

Рис. 33 – Результат инициализации триангуляционной сетки

17. Поскольку мы решаем уравнение с ЧП для динамического процесса (т.е. такого, который развивается во времени), нужно ввести начальное значение и перечень моментов времени, для которых мы хотим знать температуру. Из меню Solve (Решить) выбираем опцию Parameters, с помощью которой вызываем окно Solve Parameters (Параметры решения), рис. 34.

Рис. 34 – Параметры решения, установленные в PDE app по умолчанию

Поскольку MATLAB вычисляет очень быстро (в этом примере стационарное состояние температуры достигается всего за 0.1 с), мы не успеем проследить визуально динамику нагрева. Поэтому целесообразно ввести логарифмическую шкалу времени. Соответствующая MATLAB-функция logspace(a, b, n) вводит логарифмическую шкалу (в данном случае – шкалу времени) с декадами от 10^a до 10^b c n точками между ними. Для нашего примера можно выбрать logspace (-2, -1, 10), что обеспечит выдачу 10 точек между логарифмически размещенными отметками времени от 10^-2 (т.е. от 0.01) до 10^-1 (т.е. до 0.1).

Введём логарифмическую шкалу в соответствующее окошко - меню Solve – опция Solve Parameters (рис. 35). Также введём нулевое условие u(t0). Если граничное условие и начальное условие не совпадают, решение будет содержать разрывы непрерывности.

Рис. 35 – Корректировка параметров решения уравнения с ЧП

18. Переходим к решению уравнения (меню Solve – опция Solve PDE, т.е. решить уравнение с частными производными). По умолчанию в окно графического вывода автоматически выводится распределение температуры на введенной области SQ1+R1, которое достигается в последний момент времени, отведенного на расчет (в данном примере для введенного logspace (-2, -1, 10) – для момента времени 0.1.

19. При желании можно визуализировать динамику изменения температуры с помощью анимации. Для этого используем меню Plot (график) – опция Plot Parameters (рис. 36). Настраивать параметры вывода графика нужно до того, как будет запущена процедура решения - до нажатия опции Solve - Solve PDE.

Рис. 36 – Окно параметров графического вывода результатов

Для получения анимации нужно включить флажок Height (3-D plot) и флажок Animation. Также полезно преобразовать сетку из треугольной в прямоугольную (флажок Plot in x-y grid) – такое преобразование ускоряет расчет для графика. По умолчанию установлена цветовая палитра cool («прохладная»). Графики будут ярче, если из выпадающего меню цветов выбрать hot («горячая»); после выбора нужно нажать кнопку Close.

По умолчанию выводится пять графиков подряд. Можно изменить количество выводов (например, сделать только однократный вывод) через окошко Options диалогового окна (рис. 37).

Рис. 37 – Установка желательного числа выводов графика

Могут быть выведены различные графики с различными опциями, см. рис. 38-40.

Рис. 38 – 3D-график температуры. Для того чтобы снабдить график сеткой, нужно нажать правую кнопку мыши на поле графика и выбрать опцию Griid (сетка)

Рис. 39 – Проекции линий равного уровня температуры на координатную плоскость

Рис. 40 – Поле направлений антиградиента (антиградиент направлен в сторону наиболее сильного убывания температуры)

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

После нахождения решения в командное окно MATLAB выводится сообщение, например:

31 successful steps (31 успешных шагов)
0 failed attempts (0 ошибочных попыток)
64 function evaluations (64 вычислений функций)
1 partial derivatives (1 частная производная)
10 LU decompositions (10 декомпозиций матриц на треугольные)
63 solutions of linear systems (63 решения линейных систем)

 

3. Пример для самостоятельной работы. Решение задачи нагрева пластины внешним потоком тепла с использованием MATLAB-приложения PDE app

Постановка задачи:

Решить задачу нагрева треугольной железной пластины внешним потоком тепла интенсивностью 10 ккал/сек. Координаты трех углов пластины (в метрах, первое число относится к абсциссе, второе – к ординате): (0.0 1.5) (0.6 1.0) (0.7 2.0)

Нагрев происходит внешним потоком тепла через границу пластины по стороне, ограниченной вершинами (0.0 1.5) и (0.6 1.0). Остальные стороны пластины изолированы от внешнего влияния; температура на границах по этим сторонам равна нулю.

Внутренние источники тепла в пластине отсутствуют.

Начальные условия предполагаются нулевыми. 

Теплофизические параметры пластины, необходимые для решения:

• плотность 7.86 кг/м³;
• удельная теплоемкость 0.11 ккал/(кг×°C);
• коэффициент теплопроводности 50.5 ккал/(м×час×°C);
• коэффициент теплопередачи 79.1 ккал/(м²×час×°C);

Анализ условий задачи:

По условию нагрев пластины происходит извне, через ее границу. Следовательно, на границе по стороне, ограниченной вершинами (0.0 1.5) и (0.6 1.0), необходимо задать граничное условие второго рода (условие Неймана) с параметрами: тепловой поток (по условию равный 10 ккал/сек), переведенный в единицу измерения [ккал/час], равен 36000 ккал/час; коэффициент теплопередачи (по условию) 79.1 (ккал/(м²×час×°C). На остальных границах пластины согласно условию задачи имеет место граничное условие первого рода (условие Дирихле).

Виды окна PDE app на разных шагах решения задачи – рис. 41-45, используются как «подсказки» при самостоятельной работе.

Рис. 41 – Примерный вид чертежа рабочей области

Рис. 42 – Ввод граничного условия второго типа (условие Неймана)

Рис. 43 – Ввод параметров уравнения из меню PDE

Рис. 44 – Вывод графика проекции распределения температуры на координатную плоскость (x, y). Стрелками показан антиградиент температуры

Рис. 45 – Два ракурса 3D-графика распределения температуры с показом триангуляционной сетки

 

4. Программные средства решения уравнений с частными производными (на примере уравнения теплопроводности)

Целесообразность освоения элементов программирования для решения уравнений теплопроводности объясняется ограниченными возможностями приложения PDE app, в котором параметры геометрии и теплофизические свойства объекта могут быть только константами. На практике теплофизические характеристики объекта могут зависеть от пространственных координат и времени и, возможно, от искомой переменной. Также возможности моделирования, предоставляемые приложением PDE app, не слишком удобны (например, при изменении геометрических параметров объекта приходится перерисовывать чертеж объекта вручную и начинать решение «с нуля»).

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

Ниже в качестве примера приведена программа для решения уравнения теплопроводности на квадратной области.

Постановка задачи:

Рассматривается задача нагрева квадратной пластины; ширина стороны равна 2 см. На чертеже (рис. 46) оси проходят через центр тяжести пластины. Начальные условия заданы в центре пластины (координаты центра (x=0, y=0); начальные значения температуры равны (–0.5°C) внутри круга радиуса 0.4 см и равны 0°C вне круга. Граничные условия – нулевые (условия Дирихле).

Зоны внутренних источников тепла расположены во втором и четвертом квадранте пластины (рис. 46), значение скорости изменения температуры внутри зон равно 5°C/сек. Границы зон нагрева – окружности заданного радиуса r (см), который можно изменять в пределах: 0.1≤r≤0.45 (см). Значение радиуса вводится в ходе выполнения программы с командной строки MATLAB.

Рис. 46. Чертеж объекта нагрева. Указаны координаты центров зон расположения внутренних источников тепла и зоны ненулевых начальных условий 

Время моделирования можно изменять в пределах от 0.001 сек до достижения установившегося значения температуры (не менее 1 сек.) Нужное значение времени моделирования вводится с командной строки MATLAB.

Уравнение теплопроводности относительно температуры u(x, y, t), для решения которого используется MATLAB-функция parabolic(…), приведено в документации MATLAB в виде:

где - оператор набла, если u – скаляр, то это оператор градиента.

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

Теплофизические параметры объекта:

c(x, y, t) – коэффициент теплопроводности, вт/(м×°C);

d(x, y, t) = ρ(x, y, t) × C(x, y, t)

ρ(x, y, t) – плотность вещества пластины, кг/м³;

C(x, y, t) – удельная теплоемкость, дж/(кг×°C);

a(x, y, t) – коэффициент, учитывающий инерционность, вт/(м³×°C)

f(x ,y, t) – удельный поток тепла, вт/м³.

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

Текст программы с комментариями прикреплены файлами к данной работе.

Форма представления результатов расчета

1. Численные результаты программы - квадратная матрица значений температур u, которые получаются после завершения времени моделирования (T_max). Она может быть просмотрена в рабочей области (Workspace). Размерность матрицы определена в программе-функции heat_fun.m, командой x=linspace(-1,1,61), т.е. в данном расчете размерность 61×61. Трактовка: стороны квадратной пластины делятся на 60 отрезков (иначе – на 61 отсчет по каждой оси), и в каждом узле сетки 61×61 в матрице u содержится значение температуры. На рис. 47 показан способ вызова редактора значений и фрагмент матрицы значений температур.
   
   Рис. 47 – Вызов редактора и фрагмент числовых результатов расчета температур на поверхности пластины
   
   
2. В ходе работы программы графики, иллюстрирующие распределение температур, динамично изменяются. Левый график отображает триангуляционную сетку.
3. На рис. 48 показан вид экрана монитора после завершения расчетов. Левый график с сечениями 3D-графика размещается поверх чертежа триангуляционной сетки.

Рис. 48 – Графический вывод результатов расчета

На рис. 49 представлены результаты расчета распределения температур, достигнутые в разные моменты времени. Иллюстрируется постепенное снижение влияния начальных условий (отрицательных температур внутри круга радиуса 0.4 см в центре пластины). Также иллюстрируется постепенный разогрев области вокруг источников тепла и распространение разогрева по пластине. Верхний ряд рисунков – 3D-графики распространения тепла, нижний ряд – сечения этих графиков (сечения описаны на левом графике рис. 48)

Рис. 49 – Изменения распределений температур во времени (процессы отработки ненулевых начальных условий и разогрева пластины)

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

Рис. 50 – Иллюстрации влияния изменений коэффициента теплопроводности на распределение температуры по пластине: А) – коэффициент теплопроводности постоянный; Б) – коэффициент теплопроводности со временем увеличивается; В) – коэффициент теплопроводности уменьшается.

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

 

Заключение

1. Представлены методические материалы, достаточные как для самостоятельного изучения приложения PDE app, так и для встраивания в учебные программы изучения дисциплин «Уравнения математической физики» и «Системы компьютерной математики» в технических вузах.
2. Подробно описан порядок подготовки исходных данных к решению уравнения или системы уравнений с частными производными.
3. Приведены иллюстративные примеры моделирования процесса нагрева с распределенными параметрами.
4. В прилагаемом файле содержится программная реализация примеров, включенных в основной текст работы.
5. В приложениях содержатся справочные сведения (примеры производственных объектов с распределенными параметрами; краткая классификация уравнений с частными производными; аннотации основных методов приближенного решения).

 

Литература

1. Маркус А.А. Моделирование тепловых процессов в трубчатых вращающихся печах спекания. – Диссертация на соискание ученой степени кандидата технических наук, электронный ресурс http://av.disus.ru/dissertatciya/101123-1-modelirovanie-teplovih-processov-trubchatih-vraschayuschihsya-pechah-spekaniya.php
2. Маковский, Э. Э. Автоматизированные автономные системы трансформации неравномерного стока / Э. Э. Маковский, В. В. Волчкова. –Фрунзе: Илим, 1981. – 380 с.
3. Натальчук, М. Ф. Эксплуатация гидромелиоративных систем / М. Ф. Натальчук, В. И. Ольгаренко, В. А. Сурин. – М.: Колос, 1995. – 320 с.
4. Электронный ресурс http://helpiks.org/2-109648.html
5. Дилигенская А.Н., Данилушкин. Математическое моделирование объектов с распределенными параметрами. – Самара: Самарск. гос. техн. ун-т, 2012.– 65с.
6. Фарлоу С. Уравнения с частными производными для научных работников и инженеров: Пер с англ. – М.: Мир, 1985.– 384с. Электронный ресурс http://padaread.com/?book=170129&pg=1
7. Электронный ресурс http://empheer.16mb.com/data/t1/t1_1.html
8. Казак В.Ф. Конспект лекций по дисциплине «Уравнения математической физики», КТИ (филиал) ВолгГТУ (рукопись).
9. Белов Ю.А. Уравнения с частными производными (учебное пособие).– ФГОУ высшего профессионального образования «Сибирский федеральный университет», г. Красноярск, 2008.– 119с. Электронный ресурс
10. Чупров И.Ф., Канаева Е.А. Уравнения параболического типа и методы их решения: учебное пособие. – Ухта, УГТУ, 2012.– 193с.
11. Введенская Е.В., Горбацевич В.В., Осипенко К.Ю. Уравнения гиперболического типа: учебное пособие по курсу «Уравнения с частными производными», часть 2 – изд-во МАТИ – Российский государственный технологический университет им. К.Э. Циолковского.– М: 2013.– 24с. Электронный ресурс http://osipenko.rstu.ru/urmat2.pdf
12. Костецкая Г.С., Радченко Т.Н. Уравнения эллиптического типа: учебное пособие. – изд-во ФГОУ «Южный федеральный университет». – Ростов на Дону: 2014.– 91с. https://studylib.ru/doc/2059854/g.-s.-kosteckaya--t.-n.-radchenko-uravneniya 
13. Электронный ресурс http://alexlarin.net/Ucheb/milkov.pdf
14. Электронный ресурс http://helpiks.org/5-23013.html
15. Электронный ресурс http://mathmodel.narod.ru/2.2.htm
16. Электронный ресурс http://www.polybook.ru/comma/6.2.pdf
17. Дьяконов В.П. Mathcad 8-12 для студентов. / В.П. Дьяконов.– М.: СОЛОН-Пресс, 2005.– 631с.
18. Дьяконов В.П. MATLAB R2006/2007/2008 + Simulink 5/6/7. Основы применения. / В.П. Дьяконов.– М.: СОЛОН-Пресс, 2008.– 800с.

 

---

ПРИЛОЖЕНИЕ 1

Примеры объектов и технологических процессов с распределенными параметрами

Приведем несколько примеров распространенных производственных объектов, математические модели для которых разработаны в классе моделей с распределенными параметрами.

1. Вращающиеся печи цементных заводов (рис. П1.1), имеющие длину до 200м и диаметр до 7м.

Рис. П1.1 – Вращающиеся печи Себряковского цементного завода (г. Михайловка Волгоградской области)

2. Крупные магистральные и водораспределительные каналы ирригационных систем (рис. П1.2).

Рис. П1.2 – Чумышская плотина на р. Чу (Киргизия) с водозабором в Ат-Башинский магистральный канал (слева) и Георгиевский магистральный канал (справа).

Протяженность магистральных каналов может достигать сотен километров. Например, протяженность Крымской части Северо-Крымского канала составляет 402,6 км (в настоящее время водоподача в канал перекрыта по решению правительства Украины). При транспортировке воды в магистральном канале возникают распределенные по длине волновые процессы, связанные с отражением волн от перегораживающих и водорегулирующих сооружений. Математическое моделирование неустановившегося режима потока в магистральном канале основано на разработке и решении математических зависимостей, реализующих известные законы гидравлики. Исходной системой дифференциальных уравнений неустановившегося течения воды является система уравнений Сен-Венана с частными производными, хорошо описывающая динамические процессы течения воды в длинных призматических руслах, в том числе волновые процессы.

3. Длинные линии электропередачи (рис. П1.3). Электрическая цепь, в которой электрические сопротивления и проводимости, индуктивности и электрические емкости распределены вдоль всей цепи, называют электрической цепью с распределенными параметрами.

Рис. П1.3 – Высоковольтная линия электропередачи

Система линейных проводов, соединяющих генератор с приемником для передачи электрической энергии или сигнала, называется линией. Чтобы учесть изменения тока и напряжения вдоль линии, принимается, что каждый сколь угодно малый элемент линии по длине обладает сопротивлением и индуктивностью, а между проводами – проводимостью и емкостью, т. е. необходимо рассматривать линию как цепь с распределёнными параметрами. Линию, длина которой сравнима с длиной волны или превосходит её, называют длинной. Параметры цепи могут быть распределены неравномерно вдоль цепи. Однако во многих случаях параметры можно полагать распределенными равномерно вдоль цепи, например, для линий передач, в которых сечение проводов, их взаимное расположение и характеристики среды не изменяются по длине линии. Такие линии называют однородными длинными линиями. Эффект непрерывного изменения тока и напряжения вдоль линии возникает вследствие того, что линия обладает распределенными продольными и поперечными параметрами. Математическая модель длинной линии записывается в форме системы с частными производными, называемой телеграфными уравнениями.

В качестве других примеров объектов с распределенными параметрами можно указать: теплообменники, нефтяные и газовые трубопроводы, транспортные системы, пространственное распространение выбросов вредных веществ из источников-загрязнителей и т.п. Все эти объекты имеют очень широкую область промышленного применения, что и обуславливает актуальность освоения вычислительных методов решения уравнений с ЧП.

 

---

ПРИЛОЖНИЕ 2

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

Известно [6-13], что подходы к исследованию и нахождению решений различны для разных типов дифференциальных уравнений с частными производными. Для того чтобы облегчить инженерам-разработчикам математических моделей объектов и процессов с распределенными параметрами выбор подходящего математического аппарата, вводится классификация уравнений с ЧП.

Принципы классификации могут быть различными. Ниже рассмотрены три следующих принципа:

А) в зависимости от математической природы уравнения;
Б) в зависимости от физического или инженерного смысла решаемой задачи;
В) в зависимости от типа краевых (граничных) условий.

Пояснения к п. А. Согласно математической природе уравнения в [6] выделены шесть групп признаков:

1. Порядок уравнения;
2. Число переменных;
3. Линейность;
4. Однородность;
5. Вид коэффициентов;
6. Тип уравнения.

Последние 3 группы признаков относятся к наиболее изученной группе уравнений с ЧП – линейным уравнениям второго порядка.

Рассмотрим каждый из признаков подробнее.

1. Порядок уравнения. Порядком уравнения называется наивысший порядок частных производных, входящих в уравнение. Примеры:

– уравнение второго порядка относительно искомой переменной u(x, t), зависящей от времени t и пространственной координаты x; λ – параметр (константа):
;     (П2.1)

– уравнение первого порядка:
.     (П2.2)

2. Число переменных. Числом переменных называется число независимых переменных, от которых зависит решение. Например, число переменных в (П2.1) и (П2.2) равно двум (время t и пространственная координата x). Искомая переменная u(x, y, z, t) в уравнении (П2.3) зависит от четырех независимых переменных – времени t и пространственных координат трехмерного пространства x, y, z:

     (П2.3)

3. Линейность. Уравнения с частными производными бывают линейными и нелинейными. В левую часть линейных уравнений искомая переменная и все её частные производные входят аддитивно (т.е. в виде взвешенной суммы искомой переменной и ее частных производных, причем весовые коэффициенты не зависят ни от искомой переменной, ни от её частных производных). Например, уравнение второго порядка с ЧП (П2.4) является линейным:

     (П2.4)

В этом уравнении коэффициенты a₁₁, a₁₂, a₂₂ , b₁ и b₂ не зависят от искомой переменной u(x,y) и её производных, но могут зависеть от x и/или y. В (П2.4) F(x,y) – заданная функция, не зависящая от искомой переменной и её производных.

Если хотя бы один признак линейности нарушен, уравнение с ЧП относится к классу нелинейных. Например, таково уравнение (П2.5), в котором нарушена аддитивность (содержится квадратичный элемент):

     (П2.5)

Остальные признаки классификации записаны ниже применительно к уравнению второго порядка вида (П2.4).

4. Однородность. В однородных уравнениях правая часть равна нулю: F(x,y,z)=0 тождественно. В противном случае уравнение вида (П2.4) называется неоднородным.

5. Виды коэффициентов в (П2.4). Если коэффициенты a₁₁, a₁₂, a₂₂, b₁ и b₂ являются константами, то уравнение (П2.4) называется уравнением с постоянными коэффициентами. Если хотя бы один из коэффициентов является функцией независимых переменных (в (П2.4) – функцией x и/или y), то (П2.4) относится к виду уравнений с переменными коэффициентами.

6. Типизация уравнений второго порядка. В зависимости от сочетания значений коэффициентов a₁₁, a₁₂, a₂₂ уравнение (П2.4) относится к одному из трех типов:

• Параболическому: a₁₂² – a₁₁× a₂₂ = 0
• Гиперболическому: a₁₂² – a₁₁× a₂₂ > 0
• Эллиптическому:    a₁₂² – a₁₁× a₂₂ < 0     (П2.6)

Методы решения для перечисленных типов уравнений различны. Так, для уравнений параболического типа используются интегральные преобразования Фурье и метод разделения переменных [10]. Для уравнений гиперболического типа используются метод распространяющихся волн и метод продолжений [11]. Для уравнений эллиптического типа используются метод функции Грина и метод потенциалов [12].

Пояснения к п. Б классификации уравнений с ЧП по физическому или инженерному смыслу решаемой задачи. Согласно этой позиции классификации выделены 4 уравнения, которые называют основными уравнениями математической физики [13]. Как отмечено в [13], «их подробное изучение дает возможность построить теорию широкого круга физических явлений и решить ряд физических и технических задач». Основными уравнениями математической физики являются:

1. Волновое уравнение

     (П2.7)

Здесь c – скорость распространения волны в данной среде.

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

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

     (П2.8)

Здесь a² – коэффициент температуропроводности (физическая величина, характеризующая скорость изменения (выравнивания) температуры вещества в неравновесных тепловых процессах. Численно равна отношению теплопроводности к объёмной теплоёмкости при постоянном давлении, в системе СИ измеряется в м²/с).

Уравнение теплопроводности является уравнением параболического типа и служит основой для изучения процессов распространения тепла в однородном изотропном теле и явлений диффузии.

3. Уравнение Пуассона

     (П2.9)

Здесь f(x,y,z) – заданная функция.

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

4. Уравнение Лапласа (частный случай уравнения Пуассона, в котором f(x,y,z)=0):

     (П2.10)

Уравнение Лапласа является уравнением эллиптического типа. Это уравнение используется при моделировании многих физических явлений: в задачах механики, распределении тепла, электростатики, гидравлики. Большое значение оператор Лапласа имеет в квантовой физике, в частности в уравнении Шрёдингера.

Пояснения к п. В классификации уравнений с ЧП по типу краевых (граничных) условий. Любое из уравнений (П2.7)-(П2.10) имеет бесчисленное множество решений. Для выбора конкретного решения из этого множества должны быть заданы некоторые дополнительные условия, соответствующие смыслу решаемой задачи.

В большинстве случаем такими условиями являются:

1. Начальные условия, относящиеся к какому-то моменту времени T₀ (чаще всего принимаемому за нуль, T₀ = 0), с которого начинается решение задачи. Для этого момента времени нужно задать значение искомой функции: u(x,y,z,T₀) должно быть задано.
2. Граничные условия, т.е. условия на границах рассматриваемого объекта.

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

В представляемой работе большинство примеров относится к уравнению теплопроводности, поэтому перечислим ниже типы граничных условий именно для этого уравнения [12]. Для простоты будем рассматривать случай с единственной пространственной координатой x. Имеют место граничные условия следующих четырех типов:

Граничные условия первого типа (условия Дирихле): тепловой режим задан на всей внешней границе S тела. Если пространственная координата единственная, то эти условия задаются для температуры u(x,t) на концах нагреваемого стержня длиной L для всех моментов времени внутри рассматриваемого интервала 0 ≤ t ≤ T_f , T_f задано:

u(0, t) = α(t);     u(L, t) = β(t)     (П2.11)

Здесь α(t), β(t) – заданные функции (возможно, равные нулю).

Граничные условия первого типа (условия Неймана): предполагается, что через концы стержня во внешнюю среду подается определенный тепловой поток. Тепловой поток Q – это количество тепла, которое проходит через единицу площади за единицу времени. По закону Фурье тепловой поток равен:

     (П2.12)

Здесь k – коэффициент теплопроводности, n_i – внешняя нормаль к поверхности тела. Предположим, что пространственная координата единственная, тогда должно быть задано:

    (П2.13)

Здесь μ(t), ν(t) – заданные функции времени t.

В точке x=0 направление внешней нормали n_i к сечению стержня совпадает с отрицательным направлением оси Ox, а в точке x = L – с положительным направлением, поэтому

     (П2.14)

Используя (П2.12), (П2.13) и (П2.14), получим запись граничных условий второго рода:

     (П2.15)

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

     (П2.16)

Здесь u₀ – температура внешней среды, λ – коэффициент теплообмена. Граничные условия третьего рода записываются так:

     (П2.17)

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

 

---

ПРИЛОЖЕНИЕ 3

Обзор основных способов приближенного решения уравнений с частными производными

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

Приближённые методы решения задач математической физики делятся на две основные группы [14]:

1. Приближённо-аналитические методы – методы, в которых приближённое решение получается в аналитической форме, например, в виде отрезка некоторого ряда. Один из методов этой группы (метод конечных элементов) реализован в системе MATLAB.
2. Численные методы – методы, с помощью которых можно получить таблицу приближённых значений искомого решения в некоторых точках рассматриваемой области. Метод сеток (как один из численных методов) реализован в системе Mathcad.

П3.1. Сведения о группе приближённо-аналитических методов решения краевых задач для уравнений с ЧП. В данную группу входит метод Фурье, в котором решение представляется в виде бесконечного ряда; для приближенного решения используется конечное число элементов ряда. В эту же группу входят вариационные методы, один из которых (метод конечных элементов, МКЭ) выбран в системе MATLAB в качестве математической основы приближенных решений уравнений с ЧП.

Основную идею вариационных методов поясним на примере метода взвешенных невязок (МВН) [15]. Запишем обобщенно уравнение с ЧП с помощью оператора L(Ф)=0, где Ф – искомое решение.

Будем искать приближенное решение  в виде полинома независимой переменной: , коэффициенты которого a₀, a₁,… a_m выбраны произвольно.

Ели подставить  в исходное уравнение, то получится не нуль (поскольку решение неточное), а так называемая невязка: 

Дальше используется подход, аналогичный известному методу наименьших квадратов: найдём коэффициенты a₀, a₁,… a_m так, чтобы обеспечить минимум критерия, построенного на невязках по всей области поиска решения. Такой критерий имеет форму интеграла от невязок во всех точках области поиска решения.

Метод конечных элементов (МКЭ), входящий в группу приближенно-аналитических методов, является в настоящее время самым распространенным методом нахождения решений уравнений с ЧП. По-видимому, это связано с исключительно высоким качеством его реализации в системах компьютерной математики, избавляющих инженера от решения задач оценки скорости сходимости, точности решения, минимизации сложных критериев и другой высококвалифицированной математической работы, для которой математическая подготовка инженера недостаточна.

Идея МКЭ, как и идея взвешенных невязок, состоит в аппроксимации искомой непрерывной функции (в задачах математической физики – температуры, давления, перемещения и пр.) дискретной моделью, которая строится на основе разбиения области на подобласти (конечные элементы). В системе MATLAB форма конечных элементов – треугольники, рис. П3.1).

Рис. П3.1 – Разбиение области решения на конечные элементы в MATLAB

Для создания алгоритмического обеспечения современных систем компьютерной математики (таких как MATLAB, Mathematica) были привлечены самые выдающиеся достижения математической теории. Например, разбиение области поиска решения на конечные элементы автоматически учитывает особенности границ объектов. На рис. П3.1 показаны результаты разбиения на треугольные подобласти, причем в ответственных участках, на которых геометрия подобластей имеет изломы, размеры треугольников автоматически уменьшаются. Более того, авторы MATLAB не рекомендуют изменять разбиение, устанавливаемое автоматически (рис. П3.1 слева), чтобы избежать ненужного увеличения времени расчета и гарантировать точность решения. Но по желанию исследователя размеры конечных элементов можно уменьшить (рис. П3.1 в центре и справа).

На каждом из конечных элементов искомая функция аппроксимируется пробной функцией (обычно полиномом). Пробная функция вне «своего» элемента равна нулю. Эти функции должны гладко сопрягаться по границам соседних элементов. Параметры гладких функций рассчитываются так, чтобы интеграл от ошибок аппроксимации был минимальным. Поскольку пробная функция – полином, для определения коэффициентов получается система линейных уравнений, решение которой находится легко. Кроме того, матрица искомых коэффициентов полинома получается сильно разреженной (т.е. подавляющее число ее элементов – нули). Для разреженных матриц разработаны особые приемы, резко сокращающие затраты времени на вычисления.

На рис. П3.2 показан результат расчета распределения температуры по плоской пластине, нагреваемой извне (граничные условия второго рода). Иллюстрируется гладкое «сшивание» фрагментов решений на отдельных конечных элементах со смежными элементами.

Рис. П3.2 – Результат расчета с использованием МКЭ

П3.2. Сведения о группе численных методов решения краевых задач для уравнений с ЧП. Наиболее распространенными методами этого класса являются метод сеток и метод прямых.

Метод прямых заменяет конечно-разностной аппроксимацией только для одной из независимых координат. В результате для уравнений с единственной пространственной координатой вместо уравнений с частными производными получается система обыкновенных дифференциальных уравнений.

В методе сеток используется конечно-разностная аппроксимация производных (рис. П3.3 для задачи с одной пространственной координатой), в результате чего вместо решения уравнения с частными производными решается система алгебраических уравнений.

Рис. П3.3 – Сетка на плоскости поиска решения

Шаг сетки h (рис. П3.3) выбирается одинаковым по обеим координатным осям x, t. Можно записать следующие формулы для конечно-разностных оценок производных искомой функции u(x, t) в узле сетки k, окруженном 8-ю соседними узлами (рис. П3.3):

     (П3.1)

Возможны и иные варианты замены производных конечными разностями; для решения полученной системы алгебраических уравнений предложены различные методы. Сочетания методов конечно-разностной аппроксимации и методов решения системы алгебраических уравнений порождает многочисленные модификации метода сеток. В частности, используются два типа разностных схем: явная и неявная.

В явной разностной схеме неизвестное значение искомой функции на (n +1)-м шаге стоит только в левой части каждого уравнения системы алгебраических уравнений, а в правой части стоят значения, уже вычисленные на предыдущих шагах. Если же правая часть уравнений содержит неизвестные значения искомой функции, то схема называется неявной.

Основными этапами приближенного решения уравнения с ЧП методом сеток являются следующие:

1. Область поиска решения покрывается сеткой с числом узлов N×M. Для решения нужно определить значение искомой функции в каждом узле, т.е. отыскать N×M неизвестных.
2. Для каждого внутреннего узла (т.е. расположенного не на границе области поиска решения) записывают конечно-разностные значения производных (например, согласно формулам (П3.1)).
3. Полученная система уравнений дополняется соотношениями, записанными для дискретного представления начальных и граничных условий. Общее число уравнений должно получиться равным числу неизвестных N×M.
4. Система уравнений решается одним из выбранных численных методов. Полученное решение считается приближенным решением уравнения с ЧП. Для этого нужно доказать, что при N, M→∞ приближенное решение сходится к истинному.