В.Е.Шмелев "Заметки по использованию системы FEMLAB"

В оглавление

Формы представления дифференциальных уравнений в частных производных в системе FEMLAB. Поддержка решения краевых задач, поставленных в проекционной формулировке

Системы дифференциальных уравнений в частных производных (PDE) составляют математическую основу системы FEMLAB. Иначе такие уравнения в математике называют уравнениями математической физики. Система дифференциальных уравнений, дополненная граничными и начальными условиями, называется краевой задачей. Такое представление краевой задачи называется дифференциальной формулировкой. Начиная с самых первых версий, система FEMLAB была ориентирована на решение стационарных и нестационарных краевых задач, описываемых в виде дифференциальных формулировок. Все PDE в системе FEMLAB решаются конечноэлементными методами (FEM). Существует вариационная и проекционная формулировка этих методов. В FEMLAB поддерживаются проекционные методы конечных элементов, которые являются подмножеством методов взвешенных невязок. Существование методов взвешенных невязок позволяет описывать краевые задачи непосредственно в проекционной формулировке без записи PDE. Такая формулировка краевых задач поддерживается в FEMLAB, начиная с версии 2.3. Цель данной заметки – кратко описать формы представления краевых задач, поддерживаемые в системе FEMLAB 2.3.

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

По числу измерений PDE бывают одномерные (1D), двумерные (2D), трёхмерные (3D) и в общем случае многомерные. В системе FEMLAB поддерживаются 1D, 2D, 3D уравнения.

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

По максимальному порядку временных дифференциальных операторов различают стационарные (“эллиптические”) (без временных дифференциальных операторов), “параболические” (с временными дифференциальными операторами первого порядка), “гиперболические” (с временными дифференциальными операторами второго порядка), а также PDE более высокого порядка. Первые три типа PDE поддерживаются в системе FEMLAB. Они имеют следующие названия в этой системе: стационарные (Stationary), нестационарные (Time dependent), волновое расширение нестационарных PDE (Wave extension).

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

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

В настоящее время для решения PDE широко применяются три пакета расширения MATLAB. Одномерные матричные параболические уравнения первого порядка решаются в ядре MATLAB (функции pdepe, pdeval). Ограниченный класс двумерных матричных эллиптических, параболических и гиперболических PDE решается в PDE Toolbox. Система FEMLAB поддерживает решение всех названных выше типов PDE и краевых задач на их основе. Рассмотрим формы представления краевых задач в системе FEMLAB.

Коэффициентная форма представления краевой задачи

Сначала эту форму рассмотрим для скалярного поля. Пусть в одномерном, двумерном или трёхмерном пространстве имеется некоторая область , в которой нужно рассчитать скалярное поле. Обозначим границу этой области . Система FEMLAB позволяет моделировать скалярные поля, описываемые уравнениями вида:
, (1)
где u – искомое скалярное поле; d_a, a, f – заданные скалярные поля; , , – заданные векторные поля; c – заданное скалярное или тензорное (второй валентности) поле. Поля d_a, c, , , , a, f называют коэффициентами PDE. Иначе их называют параметрами материальных свойств среды. Скалярные переменные h, r, q, g называют коэффициентами граничных условий.

Система уравнений (1) называется коэффициентной формой представления краевой задачи (Coefficient form). Если d_a=0, то краевая задача называется стационарной (Stationary), иначе – нестационарной (Time dependent). Если первую производную по времени заменить на вторую, то (1) будет называться волновым расширением нестационарной краевой задачи (Wave extension). Второе уравнение (1) называется граничным условием первого рода (граничным условием Дирихле). Третье уравнение (1) называется граничным условием второго рода (граничным условием Неймана). Если коэффициенты PDE и граничных условий не зависят от искомого поля u и его дифференциальных операторов, то краевая задача называется линейной, в противном случае – нелинейной. Коэффициентная форма пригодна для решения линейных и несущественно-нелинейных краевых задач.

Искомое скалярное поле u в терминологии FEMLAB называется зависимой переменной. Пространственные координаты x, y, z называются независимыми переменными. Есть ещё одна особая независимая переменная для нестационарных задач – это время t. Все эти имена переменных могут записываться в выражениях FEMLAB. В этих выражениях могут встречаться также и другие переменные. Перечень их сведём в таблицу.

Интерпретацию перечисленных переменных сведём в следующую таблицу.

Общий вид диалогового окна ввода коэффициентов PDE, разворачиваемого по команде меню Subdomain/ Subdomain Settings при работе GUI-приложения femlab, показан на рис. 1.

Рис. 1. Диалоговое окно ввода коэффициентов PDE в трёхмерном режиме моделирования

Порядок ввода скалярных и векторных коэффициентов PDE отчётливо виден на рис. 1. В трёхмерном режиме тензорный коэффициент c имеет 9 компонентов. В строке редактирования c компоненты этого тензора записываются через пробелы в столбцовом порядке: c_xx c_yx c_zx c_xy c_yy c_zy c_xz c_yz c_zz . Если тензор c симметричный, то можно в строку редактирования вписывать только его верхнюю треугольную часть: c_xx c_xy c_yy c_xz c_yz c_zz . Если тензор c ортотропный (диагональный), то можно вписывать только его диагональные компоненты: c_xx c_yy c_zz . Если коэффициент c скалярный (шаровой тензор), то можно вписывать только один компонент c.

FEMLAB поддерживает также решение систем произвольного числа скалярных PDE. Рассмотрим имена переменных в случае системы двух скалярных PDE.

Коэффициентная форма в случае системы двух скалярных PDE

Матричная коэффициентная форма краевой задачи, основанной на системе скалярных PDE, имеет следующий вид:

, (2)

где в случае системы двух скалярных PDE , , – матрица тензорных коэффициентов, , (данное матричное произведение представляет собой матрицу-столбец векторных полей), , – матрицы векторных полей, – матрица-столбец векторных полей, – матрица скалярных полей, – матрица-столбец скалярных полей, , – матрицы скалярных полей на границах, , – столбцовые матрицы скалярных полей на границе. Столбцовая матрица вычисляется автоматически.

Имена и интерпретацию переменных, поддерживаемых в выражениях FEMLAB, сведём в таблицу.

Форма (2) записана для произвольного числа зависимых переменных u₁, u₂, … u_n. Имена переменных FEMLAB в этом общем случае составляются аналогично тому, как это делается для двух зависимых переменных.

Генеральная форма представления краевой задачи

В случае одного скалярного поля генеральная форма (General form) краевой задачи имеет следующий вид:

, (3)

где вектор Г и скаляры F, R, G могут выражаться через независимые и зависимые переменные и их пространственные дифференциальные операторы.

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

Интерпретацию перечисленных переменных сведём в следующую таблицу.

Общий вид диалогового окна ввода выражений, задающих параметры PDE (3), показан на рис. 2. Выражения для компонентов вектора Г записываются в строку через пробелы.

Рис. 2. Диалоговое окно ввода параметров PDE в генеральной форме (General form)

FEMLAB поддерживает также решение систем произвольного числа скалярных PDE. Рассмотрим имена переменных в случае системы двух скалярных PDE.

Генеральная форма в случае системы двух скалярных PDE

Система скалярных PDE может быть записана в матричной форме
, (4)
где , , , , т.е. [Г] – столбцовая матрица векторов, остальные из перечисленных – столбцовые матрицы скаляров.

Имена переменных FEMLAB для формы (4) и определяющие их математические выражения сведём в таблицу.

Генеральная форма пригодна для решения линейных и существенно-нелинейных краевых задач.

Представление краевой задачи в виде ослабленной
проекционной формулировки

В системе FEMLAB дифференциальные уравнения в частных производных решаются проекционным методом Галёркина с конечными элементами. Этот метод является одним из важнейших частных случаев методов взвешенных невязок. Если говорить совсем коротко, то суть этого метода заключается в том, что весовые функции здесь равны функциям формы, с помощью которых осуществляется конечноэлементная интерполяция. Пусть имеется стационарный вариант скалярной краевой задачи (1). Если к дифференциальному уравнению применить метод Галёркина с конечными элементами, то получим

, (5)
где [N] – матрица-строка функций формы конечных элементов, [N]^T – матрица-столбец весовых (пробных) функций. Такой выбор пробных функций гарантирует их линейную независимость.

В левой части матричного уравнения (5) внесём весовые функции под знак дифференциального оператора и применим теорему о дивергенции:

.

Перенесём все интегралы полученного уравнения в правую часть, учитывая обобщённый вид граничных условий , в результате получим:
. (6)

Интегральное матричное уравнение (6) называется ослабленной проекционной формулировкой (Weak form) стационарной скалярной краевой задачи (1). В данном случае (6) представляет собой преобразование коэффициентной формы в Weak form. Для поддержки этой формы GUI-приложение femlab имеет специальные прикладные режимы Weak Modes. При использовании коэффициентной и генеральной форм в диалоговых окнах ввода коэффициентов и параметров PDE и граничных условий имеются также закладки Weak. Это видно даже на рис. 1,2. Для заполнения соответствующих строк редактирования нужны специальные обозначения переменных FEMLAB. Имена этих переменных мы сейчас рассмотрим.

Имена переменных в Weak Mode для случая системы двух скалярных PDE сведём в таблицу.

u1_test, u2_test – пробные функции для зависимых переменных u1, u2;
u1x_test, u1y_test, u1z_test, u2x_test, u2y_test, u2z_test – производные пробных функций по пространственным координатам;
u1_time, u2_time – производные зависимых переменных по времени;
u1x_time, u1y_time, u1z_time, u2x_time, u2y_time, u2z_time – смешанные производные зависимых переменных по времени и по пространственным координатам.

Если краевая задача двумерная, коэффициент c=c – скаляр и зависимая переменная только одна, то для решения интегрального матричного уравнения (6) в строку редактирования weak диалогового окна Subdomain Settings нужно вписать

ux_test*(-c*ux-alx*u+gax)+uy_test*(-c*uy-aly*u+gay)+u_test*(f-bex*ux-bey*uy-a*u)

в строку редактирования weak диалогового окна Boundary Settings нужно вписать

u_test*(-q*u+g)

Здесь вместо c, alx, aly, gax, gay, f, bex, bey, a, q, g нужно писать конкретные значения или выражения для выбранных зон расчётной области или границ. Матрицы связи h и вводятся в граничные условия автоматически.

Представление краевой задачи в виде ослабленной проекционной формулировки существенно расширяет возможности моделирования нестационарных процессов. Пусть динамические члены нестационарного PDE содержат смешанные производные искомого поля по времени и по пространственным координатам, например, . Запишем выражения FEMLAB, которые нужно вписать в строку редактирования dweak диалогового окна Subdomain Settings.

В случае одномерного PDE:

ux_test*(d_c*ux_time+d_al*u_time)+u_test*d_be*ux_time

В случае двумерного PDE:

ux_test*(d_c*ux_time+d_al1*u_time)+uy_test*(d_c*uy_time+d_al2*u_time) +u_test*(d_be1*ux_time+d_be2*uy_time)

Аналогично в случае трёхмерного PDE.

Пусть граничное условие Неймана также содержит производную по времени:
.
Тогда в строку редактирования dweak диалогового окна Boundary Settings нужно вписать
-u_test*d_q*u_time

В приведённом примере вместо переменных d_c, d_al, d_be, d_al1, d_al2, d_be1, d_be2, d_q нужно вписывать конкретные значения или выражения для выбранных зон и участков границ. Указанные коэффициенты фирма Comsol планирует внедрить в коэффициентную и генеральную формы представления краевой задачи в будущих версиях FEMLAB.

Как видно, Weak form существенно расширяет возможности моделирования стационарных и нестационарных физических полей.

В оглавление