MATLAB.Exponenta
MATLAB и Simulink на русском
Технологии разработки и отладки
		сложных технических систем

Femlab

В.Е.Шмелев "Заметки по использованию системы 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 – искомое скалярное поле; da, a, f – заданные скалярные поля; , , – заданные векторные поля; c – заданное скалярное или тензорное (второй валентности) поле. Поля da, c, , , , a, f называют коэффициентами PDE. Иначе их называют параметрами материальных свойств среды. Скалярные переменные h, r, q, g называют коэффициентами граничных условий.

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

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

Области

Размерность PDE

1D

2D

3D

Точка

 

u

u

Ребро

   

u

Граница

alux, au, beu, cux, f, g, gax, nalu, ncu, nga, qu, u, ux

alux, aluy, au, beu, cux, cuy, f, g, gax, gay, nalu, ncu, nga, qu, u, ux, uy

alux, aluy, aluz, au, beu, cux, cuy, cuz, f, g, gax, gay, gaz, nalu, ncu, nga, qu, u, ux, uy, uz

Зона

abscu1x, absux, alux, au, beu, cux, dau, f, gax, u, ux, uxx

abscu1x, absux, alux, aluy, au, beu, cux, cuy, dau, f, gax, gay, u, ux, uy, uxx, uxy, uyx, uyy

abscu1x, absux, alux, aluy, aluz, au, beu, cux, cuy, cuz, dau, f, gax, gay, gaz, u, ux, uy, uz, uxx, uxy, uxz, uyx, uyy, uyz, uzx, uzy, uzz

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

Имя переменной в системе FEMLAB

Определяющее математическое выражение

abscu1x

absux

alux, aluy, aluz

au

beu

cux

cuy

cuz

dau

da u

f

f

g

g

gax, gay, gaz

nalu

ncu

nga

qu

q u

u

u

ux, uy, uz

uxx, uxy, uxz, uyx, uyy, uyz, uzx, uzy, uzz

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


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

Порядок ввода скалярных и векторных коэффициентов PDE отчётливо виден на рис. 1. В трёхмерном режиме тензорный коэффициент c имеет 9 компонентов. В строке редактирования c компоненты этого тензора записываются через пробелы в столбцовом порядке: cxx cyx czx cxy cyy czy cxz cyz czz . Если тензор c симметричный, то можно в строку редактирования вписывать только его верхнюю треугольную часть: cxx cxy cyy cxz cyz czz . Если тензор c ортотропный (диагональный), то можно вписывать только его диагональные компоненты: cxx cyy czz . Если коэффициент c скалярный (шаровой тензор), то можно вписывать только один компонент c.

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

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

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

, (2)

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

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

Имя переменной в системе FEMLAB

Определяющее математическое выражение

abscu1x

abscu2x

absu1x

absu2x

alu1x, alu1y, alu1z

alu2x, alu2y, alu2z

au1

au2

beu1

beu2

cu1x

cu1y

cu1z

cu2x

cu2y

cu2z

dau1

da11 u1+da12 u2

dau2

da21 u1+da22 u2

f1

f1

f2

f2

g1

g1

g2

g2

ga1x, ga1y, ga1z

ga2x, ga2y, ga2z

nalu1

nalu2

ncu1

ncu2

nga1

nga2

qu1

q11 u1+ q12 u2

qu2

q21 u1+ q22 u2

u1

u1

u2

u2

u1x, u1y, u1z

u2x, u2y, u2z

u1xx, u1xy, u1xz, u1yx, u1yy, u1yz, u1zx, u1zy, u1zz

u2xx, u2xy, u2xz, u2yx, u2yy, u2yz, u2zx, u2zy, u2zz

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

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

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

, (3)

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

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

Области

Размерность PDE

1D

2D

3D

Точка

 

u

u

Ребро

   

u

Граница

f, g, gax, nga, u, ux

f, g, gax, gay, nga, u, ux, uy

f, g, gax, gay, gaz, nga, u, ux, uy, uz

Зона

absga1x, absux, dau, f, gax, u, ux, uxx

absga1x, absux, dau, f, gax, gay, u, ux, uy, uxx, uxy, uyx, uyy

absga1x, absux, dau, f, gax, gay, gaz, u, ux, uy, uz, uxx, uxy, uxz, uyx, uyy, uyz, uzx, uzy, uzz

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

Имя переменной в системе FEMLAB

Определяющее математическое выражение

absga1x

absux

gax, gay, gaz

dau

da u

f

F

g

G

nga

u

u

ux, uy, uz

uxx, uxy, uxz, uyx, uyy, uyz, uzx, uzy, uzz

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


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

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

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

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

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

Имя переменной в системе FEMLAB

Определяющее математическое выражение

absga1x

absga2x

absu1x

absu2x

ga1x, ga1y, ga1z

ga2x, ga2y, ga2z

dau1

da11 u1+da12 u2

dau2

da21 u1+da22 u2

f1

F1

f2

F2

g1

G1

g2

G2

nga1

nga2

u1

u1

u2

u2

u1x, u1y, u1z

u2x, u2y, u2z

u1xx, u1xy, u1xz, u1yx, u1yy, u1yz, u1zx, u1zy, u1zz

u2xx, u2xy, u2xz, u2yx, u2yy, u2yz, u2zx, u2zy, u2zz

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

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

В системе 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 сведём в таблицу.

Размерность PDE

1D

2D

3D

Имена переменных

u1_test, u2_test, u1x_test, u2x_test, u1, u2, u1x, u2x, u1_time, u2_time, u1x_time, u2x_time

u1_test, u2_test, u1x_test, u1y_test, u2x_test, u2y_test, u1, u2, u1x, u1y, u2x, u2y, u1_time, u2_time, u1x_time, u1y_time, u2x_time, u2y_time

u1_test, u2_test, u1x_test, u1y_test, u1z_test, u2x_test, u2y_test, u2z_test, u1, u2, u1x, u1y, u1z, u2x, u2y, u2z, u1_time, u2_time, u1x_time, u1y_time, u1z_time, u2x_time, u2y_time, u2z_time

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 существенно расширяет возможности моделирования стационарных и нестационарных физических полей.

В оглавление


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

Система Orphus

Яндекс.Метрика