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

В оглавление

Интерполяция измеренных материальных свойств и ввод их в Femlab

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

Интерполяция относительно пространственных координат параметров свойств неоднородных материалов

Пусть распределение какого-либо параметра материальных свойств в расчётной области получено экспериментальным путём и задано таблично. Пусть эта таблица хранится в базовой рабочей области MATLAB. Пусть x1, y1, z1 – массивы пространственных координат точек измерения; f1 – массив измеренных значений некоторого параметра материальных свойств. Пусть таблично заданное распределение этого параметра нужно ввести в модель Femlab. Для простоты будем полагать, что массивы x1, y1, z1, f1 описывают распределение коэффициента f (правую часть) PDE в скалярной “уравнение-основанной” краевой задаче (PDE Mode, Coefficient form, Linear stationary, 3D).

Первое, что нужно сделать – это внедрить массивы x1, y1, z1, f1 в модель Femlab. Для этого выполним команду меню Options/ Add/Edit Constants. По данной команде развернётся диалоговое окно редактирования имён предопределённых констант Femlab и соответствующих определяющих выражений. Далее выполняем следующую последовательность действий:

1. В строку редактирования “Name of constant” вписываем, например, имя xx1, в строку редактирования “Expression” вписываем выражение evalin(‘base’,’x1’), нажимаем кнопки “Set” и “Apply”.
2. В строку редактирования “Name of constant” вписываем, например, имя yy1, в строку редактирования “Expression” вписываем выражение evalin(‘base’,’y1’), нажимаем кнопки “Set” и “Apply”.
3. В строку редактирования “Name of constant” вписываем, например, имя zz1, в строку редактирования “Expression” вписываем выражение evalin(‘base’,’z1’), нажимаем кнопки “Set” и “Apply”.
4. В строку редактирования “Name of constant” вписываем, например, имя ff1, в строку редактирования “Expression” вписываем выражение evalin(‘base’,’f1’), нажимаем кнопки “Set” и “Apply”.
5. Нажатием кнопки “OK” закрываем диалоговое окно Add/Edit Constants.

На этом внедрение массивов x1, y1, z1, f1 в модель Femlab закончено.

Теперь предположим, что геометрия прорисована, граничные условия заданы. Нужно только задать коэффициенты PDE в расчётной области. Для этого выполняем следующую последовательность действий:

1. Командой меню Subdomain/ Subdomain Mode переходим в режим ввода и редактирования коэффициентов PDE.
2. Командой меню Subdomain/ Subdomain Settings разворачиваем диалоговое окно редактирования коэффициентов PDE и выделяем нужные зоны расчётной области.
3. В строку редактирования коэффициента f вписываем выражение interp3(xx1,yy1,zz1,ff1,x,y,z,'spline').
4. В строки редактирования остальных коэффициентов вписываем необходимые выражения.
5. Нажатием кнопки “OK” закрываем диалоговое окно Subdomain Settings.

Все представленные здесь выражения являются выражениями MATLAB. Функции evalin и interp3 являются функциями ядра MATLAB. В режиме 2D функция interp3 заменяется на interp2. После выполнения указанных действий можно переходить к следующим этапам моделирования: генерации сетки, решению PDE, визуализации и постпроцессорной обработке.

Пример решения двумерной задачи в PDE Mode

Пусть требуется рассчитать распределение скалярного электрического потенциала в квадратной расчётной области –1<=x<=1, –1<=y<=1. Пусть это скалярное поле обусловлено неоднородным распределением объёмной плотности электрического заряда при нулевых граничных условиях Дирихле. Скалярную краевую задачу электростатики в PDE Mode решим при аналитическом и при табличном задании объёмной плотности заряда, а затем сравним полученные решения. Более того, для сравнительного анализа можно задать разные методы интерполяции табличных данных.

Уравнение электростатики относительно скалярного электрического потенциала имеет вид:

–div(grad–P_r/) = /, (1)

где – скалярный электрический потенциал;
– относительная диэлектрическая проницаемость вещества;
P_r – векторное поле остаточной электрической поляризованности вещества;
– абсолютная диэлектрическая проницаемость вакуума (=8,85419*10^–12Ф/м);
– объёмная плотность электрического заряда. Пусть в рассматриваемой системе =1, P_r=0, /= exp(x+y).

Для прорисовки геометрии расчётной области выполним следующие действия. Создадим объект “прямоугольник” в рабочей области MATLAB:

R1=rect2(-1,1,-1,1);

Выполним команду меню File/ Insert from Workspace/ Geometry Object(s), в строку редактирования в развернутом диалоговом окне впишем имя вставляемого геометрического объекта R1 (см. рис.1).

Рис. 1. Вид диалогового окна вставки геометрических объектов

Вид поля axes GUI-приложения Femlab после вставки геометрического объекта показан на рис. 2.

Рис. 2. Вид геометрии расчётной области в режиме Draw Mode

Поскольку нулевые граничные условия Дирихле задаются по умолчанию, из Draw Mode перейдём сразу в Subdomain Mode и аналитически зададим неоднородное распределение правой части уравнения (1) с помощью команды меню Subdomain/ Subdomain Settings (вид соответствующего диалогового окна см. на рис. 3).

Рис. 3. Аналитическое задание неоднородного распределения коэффициентов PDE

Теперь можно генерировать сетку и решать PDE. По умолчанию сетка будет состоять из Лагранжевых элементов второго порядка. Решение выполним на сетке, параметры которой показаны на рис. 4.

Рис. 4. Параметры конечноэлементной сетки.

Цветовая визуализация с линиями равного уровня скалярного поля при аналитьическом задании правой части PDE показана на рис. 5.

Рис. 5. Вид решения при аналитическом задании правой части

Массивы таблично заданного распределения правой части (1) получим, выполняя в командном окне MATLAB следующую последовательность операторов:

[x1,y1]=meshgrid(-1:0.2:1,-1:0.2:1);

f1=exp(x1+y1);

Теперь массивы таблично заданного распределения правой части PDE нужно внедрить в модель Femlab, выполнив команду меню Options/ Add/Edit Constants (см. рис. 6).

Рис. 6. Вид диалогового окна определения констант модели Femlab

Таблично заданное распределение объёмной плотности электрического заряда введём в модель Femlab командой меню Subdomain/ Subdomain Settings (см. рис. 7).

Рис. 7. Кусочно-линейная интерполяция таблично заданных коэффициентов PDE

Вид решения PDE при кусочно-линейной интерполяции таблично заданных коэффициентов показан на рис.8.

Рис. 8. Вид решения при кусочно-линейной интерполяции табличных данных

На рис. 9 показано пространственное распределение ошибки решения, представленного на рис. 8.

Рис. 9. Ошибка решения при кусочно-линейной интерполяции табличных данных

Вид решения PDE при сплайновой интерполяции таблично заданных коэффициентов показан на рис.10.

Рис. 10. Вид решения при сплайновой интерполяции табличных данных

На рис. 11 показано пространственное распределение ошибки решения, представленного на рис. 10.

Рис. 11. Ошибка решения при сплайновой интерполяции табличных данных

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

В оглавление