Femlab 2.3. Руководство пользователя (перевод с английского с редакторской правкой В.Е.Шмелева):
1.5. Моделирование в системе FEMLAB

В оглавление книги \ К следующему разделу \ К предыдущему разделу

1.5.2. Мультифизика: термоэлектрические эффекты

В этом примере моделируется термоэлектрическое нагревание в медной пластине. Материал нагревается, когда через него протекает электрический ток, преодолевая электрическое сопротивление. С другой стороны, удельное электрическое сопротивление материала изменяется при изменении температуры. В этом примере сначала рассматривается модель постоянного электрического поля в проводящей среде в прикладном режиме 2D/ Physics Modes/ Conductive media DC. Затем к созданной модели добавляется прикладной режим моделирования теплопередачи 2D/ Physics Modes/ Heat transfer/ Time dependent. Объединение этих двух прикладных режимов даёт мультифизический режим.

Двумерная модель электрического поля в проводящей среде

Пусть имеется квадратная медная пластина 1 x 1 м с круглым отверстием радиуса 0.1 м в центре. К двум противоположным сторонам медного квадрата от внешнего источника подано электрическое напряжение 0.1 В. Сторонние источники тока в системе отсутствуют. Требуется рассчитать и отобразить графически распределение потенциала в системе, а также распределение напряжённости электрического поля, плотности тока и объёмной плотности мощности тепловых потерь.

Навигатор моделей

Сразу после запуска GUI-приложения femlab будет развёрнуто окно Навигатора моделей (закладка New). В группе Dimension выбираем радиокнопку 2D, делаем двойной щелчок мышью по группе прикладных режимов Physics Modes, выбираем там режим Conductive media DC (рис. 1.68). Нажатие кнопки OK приведёт к закрытию Навигатора моделей и переходу GUI-приложения femlab в требуемый прикладной режим.

Рис. 1.68. Выбор прикладного режима моделирования постоянного электрического поля в проводящей среде

Параметры и настройки

Прежде чем приступать к прорисовке геометрии, нужно настроить параметры координатной сетки. Для этого выполним команду меню Options/ Axes/Grid Settings. Эта команда развёртывает диалоговое окно настройки параметров координатной сетки. Вид этого окна с установленными значениями параметров показан на рис. 1.69.

Рис. 1.69. Окно настройки параметров координатной сетки (закладка Axis)

Для удобства прорисовки геометрии шаг координатной сетки зададим равным 0.1 м. Вид диалогового окна Axes/Grid Settings с установленными значениями параметров показан на рис. 1.70. Установив нужные значения параметров в обеих закладках, можно нажимать кнопку OK.

Рис. 1.70. Окно настройки параметров координатной сетки (закладка Grid)

Группа команд меню Options позволяет также создавать новые переменные модели. Для удобства моделирования создадим несколько констант. Выполним команду меню Options/ Add/Edit Constants, в результате чего развернётся диалоговое окно создания и редактирования констант.

Пусть нужно определить следующие константы модели:

r0 - удельное электрическое сопротивление меди при температуре T0, r0=1.754e-8 Ом м;
T0 - температура, при которой известны электрофизические свойства меди, T0=293 К;
alpha - температурный коэффициент сопротивления меди, alpha=0.0039 К^-1;
T - температура, пока не рассчитано тепловое поле, пусть T=293 К.

В диалоговом окне Add/Edit Constant в строку редактирования Name of constant впишем имя константы r0, в строку редактирования Expression впишем определяющее выражение. Определяющим выражением может быть любое выражение MATLAB, операндами которого могут быть числовые константы MATLAB и константы FEMLAB. В данном случае туда нужно вписывать числовую константу 1.754e-8. Нажатие кнопки Set приведёт к включению определяемой константы в список констант FEMLAB данной модели. Теперь в строку редактирования Name of constant впишем имя T0 (старое имя стираем), в строку редактирования Expression впишем число 293 (старое стираем). Опять нажатием кнопки Set включаем новую константу в список. Аналогичные операции проделаем для оставшихся имён констант. Нажатие кнопки Apply приведёт к вычислению всех определяющих выражений. Значения этих выражений будут показаны в списке констант (рис. 1.71). Нажатие кнопки OK приведёт к включению списка констант в модель и закрытию диалогового окна Add/Edit Constant.

Рис. 1.71. Диалоговое окно создания и редактирования констант

Режим прорисовки геометрии

Рис. 1.72. Фигура femlab с прорисованной расчётной областью

Для создания квадрата нажмём кнопку инструментальной панели прорисовки геометрических объектов. Указатель мыши поместим в точку (0,1), нажмём левую кнопку и переместим указатель в точку (1,0). Отпускание кнопки мыши приведёт к созданию геометрического объекта. Чтобы создать круг в центре квадрата, нажмём кнопку инструментальной панели прорисовки геометрических объектов. Указатель мыши поместим в точку (0.5,0.5), нажмём левую кнопку, удерживая клавишу Ctrl, и переместим указатель в точку (0.6,0.5). Отпускание кнопки мыши приведёт к созданию геометрического объекта - круга. Нажав клавишу Ctrl-A, выделим оба геометрических объекта. Теперь нажатием кнопки инструментальной панели прорисовки геометрических объектов создадим композиционный объект, представляющий расчётную область модели. На рис. 1.72 показан вид фигуры femlab с созданным композиционным геометрическим объектом.

Режим задания граничных условий

Рис. 1.73. Задание ненулевых граничных условий Дирихле

Рис. 1.74. Задание нулевых граничных условий Неймана

Чтобы перейти в этот режим, надо нажать кнопку главной инструментальной панели. Щелчком мыши выделим левый вертикальный граничный сегмент. Выполним команду меню Boundary/ Boundary Settings, в результате чего раскроется диалоговое окно задания граничных условий. В списке граничных сегментов выделенным окажется сегмент № 1. В группе элементов управления Boundary Coefficients включим радиокнопку V и в соответствующую строку редактирования впишем значение скалярного электрического потенциала, равное 0.1 (в вольтах), см. рис. 1.73. Далее, не закрывая диалогового окна, щелчком мыши в фигуре femlab выделим правый вертикальный граничный сегмент. В списке граничных сегментов выделенным окажется сегмент № 4. В группе элементов управления Boundary Coefficients включим радиокнопку V=0. Теперь выделим все сегменты, кроме первого и четвёртого, и в группе элементов управления Boundary Coefficients включим радиокнопку g=q=0 (рис. 1.74). Нажатие кнопки OK приведёт к принятию граничных условий и закрытию диалогового окна.

Режим задания коэффициентов PDE и материальных свойств

Чтобы перейти в этот режим, надо нажать кнопку главной инструментальной панели. Выполнив команду меню Subdomain/ Subdomain Settings, раскрываем диалоговое окно ввода и редактирования коэффициентов PDE (параметров материальных свойств). В списке Domain selection выделим одну единственную зону и в строки редактирования, соответствующие параметрам материальных свойств, впишем нужные значения или выражения. В строку редактирования Q впишем 0.0, т.к. в системе нет сторонних источников тока. В строку редактирования впишем 1/(r0*(1+alpha*(T-T0))), т.к. это выражение приближённо отражает температурную зависимость удельной электрической проводимости металла, обладающего проводниковыми свойствами. На рис. 1.75 показан вид диалогового окна Subdomain Settings с заданными значениями материальных свойств. Нажатие кнопки OK приведёт к принятию коэффициентов материальных свойств и закрытию диалогового окна.

Рис. 1.75. Диалоговое окно Subdomain Settings с заданными коэффициентами PDE

Режим генерации конечноэлементной сетки

Для перехода в этот режим достаточно нажать кнопку главной инструментальной панели. Если сетка не построена, то при переходе в Mesh Mode происходит инициализация сетки. Можно также принудительно инициализировать сетку, нажав кнопку главной инструментальной панели. Если сетка слишком грубая, то можно её переопределить, нажав кнопку той же панели. На рис. 1.76 показана фигура femlab с изображением конечноэлементной сетки, которая после инициализации была переопределена один раз.

Рис. 1.76. Фигура femlab с изображением конечноэлементной сетки

Решение

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

Визуализация решения и постпроцессорная обработка

Рис. 1.77. Картина распределения потенциала, получаемая сразу после решения

По умолчанию в GUI-приложении femlab установлен режим автоматической визуализации решения. По окончании решения в поле axes будет показано в цвете распределение скалярного электрического потенциала в исследуемой пластине (режим Surface Plot). На рис. 1.77 показана картина распределения потенциала, получаемая сразу после решения.

Часто бывает удобно отображать распределение решения в виде изолиний. Добавим к режиму визуализации Surface Plot режим Contour Plot. Для этого выполним команду меню Post/ Plot Parameters. Эта команда разворачивает диалоговое окно Plot Parameters, имеющее несколько закладок. Откроем закладку Contour. Установим в ней флажок Contour plot. В группе Contour data установим флажок Max/ Min marker. В группе Contour parameters в ниспадающем меню Contour color выберем пункт inverted background, что соответствует чёрному цвету изолиний. В строку редактирования Contour levels впишем число 19, что соответствует шагу по потенциалу, равному 0.005 В. Вид закладки Contour с установленными параметрами показан на рис. 1.78. Нажатие кнопки OK приведёт к построению картины, изображённой на рис. 1.79.

Рис. 1.78. Диалоговое окно Plot Parameters (закладка Contour)

Рис. 1.79. Картина распределения потенциала с изолиниями

Чтобы построить картину распределения модуля напряжённости электрического поля, нужно в ниспадающих меню диалогового окна Plot Parameters: Surface expression и Contour expression выбрать пункт electric field (E). Чтобы построить картину распределения модуля плотности тока нужно в этих меню выбрать пункт current density (J). Поскольку удельная электрическая проводимость распределена в системе однородно, напряжённость электрического поля и плотность тока пропорциональны друг другу. Поэтому для краткости покажем только напряжённость электрического поля (рис. 1.80), В/м.

Рис. 1.80. Распределение модуля напряжённости электрического поля

Чтобы построить картину распределения объёмной плотности мощности тепловых потерь при протекании электрического тока через пластину, нужно в строки редактирования диалогового окна Plot Parameters: Surface expression и Contour expression вписать выражение (Vx^2+Vy^2)/(r0*(1+alpha*(T-T0))). Получаемая при этом картина показана на рис. 1.5.2.1. Единица измерения - Вт/м³. Тот же эффект будет, если в указанные строки редактирования вписать выражение Vx*cux+Vy*cuy.

Рис. 1.5.2.1. Картина распределения объёмной плотности мощности тепловых потерь

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

Как известно, электрический ток производит теплоту в медной пластине, так что попробуем смоделировать этот эффект в прикладном режиме Heat transfer, используя возможности мультифизики FEMLAB.

Пусть верхняя и нижняя грани медной пластины абсолютно теплоизолированы от окружающей среды. Пусть все остальные грани охлаждаются воздушным потоком при температуре 300 К.

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

Параметры и настройки

Прежде чем добавлять в модель второй прикладной субрежим, удалим переменную T из списка констант, определённых пользователем. Для этого выполним команду меню Options/ Add/Edit Constants, в развёрнутом диалоговом окне выделим константу T и нажатием кнопки Delete удалим эту константу.

Навигатор моделей

Для добавления второго прикладного субрежима выполним команду меню Multiphysics/ Add/Edit Modes. Эта команда раскрывает закладку Multiphysics Навигатора моделей. В левом списке прикладных режимов выбираем Heat transfer, нажимаем кнопку >> и в ниспадающем меню Solver type выбираем Time dependent. Вид закладки Multiphysics Навигатора моделей с выбором всех необходимых настроек показан на рис. 1.5.2.2. Нажатие кнопки OK добавляет новый прикладной субрежим в модель и закрывает окно Навигатора моделей.

Рис. 1.5.2.2. Закладка Multiphysics Навигатора моделей

Режим задания граничных условий

Командой меню Boundary/ Boundary Settings раскрываем диалоговое окно Boundary Settings. На границах № 2 и 3 задаём граничное условие (полная теплоизоляция). На всех остальных границах задаём граничное условие Дирихле T=300 (см. рис. 1.5.2.3).

Рис. 1.5.2.3. Диалоговое окно Boundary Settings с граничными условиями Дирихле для теплового поля

Режим задания коэффициентов PDE и материальных свойств

Вид диалогового окна Subdomain Settings с введёнными значениями параметров материальных свойств (закладка Coefficients) показан на рис. 1.5.2.4. Начальные условия задаются в этом же диалоговом окне (закладка Init, см. рис. 1.5.2.5).

Рис. 1.5.2.4. Диалоговое окно Subdomain Settings (закладка Coefficients)

Рис. 1.5.2.5. Диалоговое окно Subdomain Settings (закладка Init)

Коэффициенты материальных свойств, показанные на рис. 1.5.2.4, заданы в единицах СИ: - плотность вещества, кг/м³; С - удельная теплоёмкость, Дж/(кг К), k - коэффициент теплопроводности, Вт/(м² К); Q - объёмная плотность мощности тепловых источников, Вт/м³.

Прикладной режим теплопередачи - нестационарный. Прикладной режим электрического поля также приводится к нестационарному, поэтому для него тоже нужно задавать начальные условия. Командой меню Multiphysics/ Conductive Media DC переключимся в прикладной режим электрического поля. В диалоговом окне Subdomain Settings (закладка Init) зададим начальное условие для скалярного электрического потенциала V=0.1*(1-x). После этого переключимся в субрежим Heat Transfer. Теперь можно настраивать параметры решателя.

Решение

Выполним команду меню Solve/ Parameters. Раскроется диалоговое окно Solver Parameters. Раскрываем закладку Timestepping. В строку редактирования Output times введём выражение MATLAB: linspace(0,3000,41). Нажатие кнопки OK приведёт к принятию изменённых настроек и закрытию диалогового окна Solver Parameters. Теперь можно решать задачу. Для этого нажмём кнопку главной инструментальной панели. Во время процесса решения будет раскрыто окно ожидания завершения вычислений (рис. 1.5.2.6).

Рис. 1.5.2.6. Окно ожидания завершения вычислений

Визуализация решения и постпроцессорная обработка

Если в строки редактирования Surface expression и Contour expression вписать имя переменной T, то получится картина распределения температуры при t=3000 с (рис. 1.5.2.7).

Рис. 1.5.2.7. Картина распределения температуры при t=3000 с

Если в строки редактирования Surface expression и Contour expression вписать имя переменной Q_ht, то получится картина распределения объёмной плотности мощности тепловых потерь при t=3000 с (рис. 1.5.2.8). Как видно, это распределение отличается от распределения, изображённого на рис. 1.5.2.1.

Рис. 1.5.2.8. Картина распределения объёмной плотности мощности тепловых потерь при t=3000 с

Рис. 1.5.2.9. Картина распределения плотности потока тепловой мощности

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

Есть возможность одновременного построения трёхмерного графика градиента температуры и цветового поверхностного графика плотности тока. Для этого в диалоговом окне Plot Parameters нужно открыть закладку Surface, в строку редактирования Surface expression вписать имя переменной J, а в строку редактирования Height expression вписать имя переменной flux_ht. Фигура femlab с графиком показана на рис. 1.5.2.10.

Рис. 1.5.2.10. Трёхмерный график, у которого высота и цвет соответствуют разным физическим величинам

Выбор формы решения и решателя

Прежде чем решать нелинейную задачу, обычно бывает нужно скорректировать параметры решателя. Выбор решателя и настройка его параметров осуществляется командой меню Solve/ Parameters, которая раскрывает диалоговое окно Solver Parameters.

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

При решении линейной задачи предпочтительной является коэффициентная форма, т.к. она даёт более быстрое решение.

Закладка Timestepping диалогового окна Solver Parameters позволяет настраивать параметры алгоритма продвижения по времени при решении нестационарных линейных и нелинейных задач.

Стационарные нелинейные задачи могут решаться как в генеральной, так и в коэффициентной форме, однако предпочтительной является генеральная форма. Выбор коэффициентной формы может нарушать сходимость итерационных алгоритмов. Если всё же выбрана коэффициентная форма, то в закладке General в ниспадающем меню Jacobian нужно выбирать пункт Numeric. Коэффициентную форму в нелинейных задачах приходится задавать, например, если модель представлена m-файлом коэффициентов PDE.

Подходы к мультифизическому моделированию

Модель термоэлектрических эффектов поясняет один из путей мультифизического моделирования в GUI-приложении femlab. Здесь сначала в закладке New Навигатора моделей выбирался один из субрежимов, рисовалась геометрия, устанавливались граничные условия и параметры материальных свойств, генерировалась сетка и получалось частное решение. Затем добавлялись другие прикладные субрежимы, решение которых также доводилось до конца.

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

В диалоговом окне Solver Parameters тоже есть закладка Multiphysics. В этой закладке есть группа элементов управления Solve for variables. GUI-приложение femlab рашает только те субрежимы, которые в этой группе являются выделенными. Значения невыделенных переменных принимаются равными либо начальным условиям, либо значениям, полученным в ходе предыдущего решения. Таким образом, у пользователя есть возможность совершенно произвольно разрывать связи между отдельными субрежимами мультифизической модели.

Иногда бывает необходимо дезактивировать прикладные режимы в подобластях (зонах) расчётной области. Для этого в закладке Coefficients диалогового окна Subdomain Settings имеется флажок Active in this domain. По умолчанию обычно этот флажок бывает установлен. Если он сброшен, то данный прикладной субрежим не действует в выделенных зонах.

В оглавление книги \ К следующему разделу \ К предыдущему разделу