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

Femlab

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

В оглавление

Трёхмерная скалярная краевая задача магнитостатики: цилиндрическая катушка с ферромагнитным сердечником

Пусть имеется цилиндрическая катушка, имеющая следующие геометрические параметры: внутренний радиус слоя обмотки Rвн = 25 мм; наружный радиус слоя обмотки Rн = 30 мм; радиус сердечника Rс = 20 мм; ось сердечника и ось катушки совпадают с осью z системы координат; высота сердечника и катушки hс = 100 мм. Катушка с сердечником расположены в пространстве симметрично относительно плоскости x0y. Расчётная область – прямоугольный параллелепипед -60<=x<=60 мм, -60<=y<=60 мм, 0<=z<=150 мм. Плоскость симметрии x0y имеет постоянное значение скалярного магнитного потенциала (пусть оно будет равно нулю), значит, весь расчёт можно выполнять только для верхнего полупространства. Плотность тока в слое обмотки Jо = 1 А/мм2.

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

Известно, что магнитостатическое поле описывается следующими уравнениями:

rot H = J, (1)
div B = 0, (2)
B = mua H + Br, (3)
где H – вектор напряжённости магнитного поля, А/мм; J – вектор плотности тока, А/мм2; B – вектор магнитной индукции, Тл; mua – абсолютная магнитная проницаемость среды, мкГн/мм; mua = mu0 mu; mu0 = 4E-4*pi мкГн/мм – абсолютная магнитная проницаемость вакуума; mu – относительная магнитная проницаемость среды; Br = mu0 Mr – вектор остаточной магнитной индукции, Тл (для областей, занятых постоянными магнитами); Mr – вектор остаточной намагниченности вещества, А/мм.

Уравнение (1) называется законом полного тока, (2) – уравнением непрерывности линий магнитной индукции, (3) – уравнением материальной связи. Для расчёта магнитостатического поля уравнения (1), (2), (3) нужно решать совместно. Конечноэлементная технология плохо пригодна для непосредственного решения этих уравнений, поэтому для выполнения практических расчётов вводят понятия потенциалов поля. Из уравнений (1), (2), (3) можно получить уравнение магнитостатики относительно векторного магнитного потенциала. Такой метод моделирования поддерживается как в ядре FEMLAB, так и в модуле Электромагнетизма. Недостаток этого метода моделирования заключается в том, что в итоге для трёхмерной задачи приходится решать систему трёх скалярных PDE относительно трёх компонентов векторного магнитного потенциала.

Выведем уравнение магнитостатики относительно скалярного магнитного потенциала.

Из (1) следует уравнение непрерывности линий плотности полного тока div J = 0, значит, плотность тока можно представить как ротор некоторого векторного поля:

J = rot Mф , (4)
где Mф – фиктивная намагниченность вещества, этим векторным полем можно эквивалентировать электрические токи. Выражение (4) подставим в (1), отбросив значок ротора из обеих частей, получим:
H = Mф – grad vH, (5)
где vH – скалярный магнитный потенциал, измеряемый в единицах тока (амперы).
Подставив (5) в (3), а затем (3) в (2), получим скалярное уравнение магнитостатики:
–div(mua grad vH – (mua Mф + Br)) = 0 (6)

Для удобства расчётов в GUI-приложении femlab обе части уравнения (6) можно разделить на mu0:

–div(mu grad vH – (mu Mф + Mr)) = 0 (7)

Немагнитные материалы обладают линейными магнитными свойствами, их магнитная проницаемость не зависит от напряжённости магнитного поля. Пусть для простоты магнитные свойства ферромагнитного сердечника описываются уравнением mu = 1 + (mu1 – 1)*H0/(H0+H), где mu1 = 500; H0 = 0,5 А/мм.

Моделирование

После запуска GUI-приложения femlab в Навигаторе моделей выбираем режим 3D, PDE Modes/ Coefficient/ Nonlinear Stationary и нажимаем кнопку OK. После этого в командном окне MATLAB выполним следующую последовательность команд:

>> R0=rect2(-60,60,-60,60); % Горизонтальное сечение всей расчётной области
>> Cs=circ2(0,0,20); % Горизонтальное сечение сердечника
>> C1=circ2(0,0,25); % Внутренний радиус слоя обмотки
>> C2=circ2(0,0,30); % Наружный радиус слоя обмотки

В GUI-приложении femlab выполним команду меню Draw/ Add/Edit/Delete Work Plane, в результате развернётся диалоговое окно, изображённое на рис. 1. Поскольку нам нужна рабочая плоскость z=0, можно сразу нажимать кнопку OK в этом диалоговом окне. femlab перейдёт в режим показа этой рабочей плоскости.


Рис. 1. Вид диалогового окна задания рабочей плоскости

Теперь выполним команду File/ Insert from Workspace/ Geometry Object(s), в результате развернётся диалоговое окно, показанное на рис. 2. Ввод в строку редактирования имён R0 C1 C2 Cs и нажатие кнопки OK приведёт к вставке геометрических объектов из рабочей области MATLAB в рабочую плоскость.


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

Вид рабочей плоскости после вставки геометрических объектов показан на рис. 3.


Рис. 3. Геометрические объекты в рабочей плоскости

Выделим объект R0 и экструдируем его на высоту 150 мм (команда Draw/ Extrude). Объекты C1, C2, Cs экструдируем на высоту 50 мм. В результате получим трёхмерное изображение, показанное на рис. 4.


Рис. 4. Прорисованная трёхмерная геометрия

Выполним команду Boundary/ Boundary Mode и перейдём в режим ввода граничных условий. На плоскости z=0 зададим граничное условие Дирихле vH = 0 (u=0). На остальных граничных плоскостях зададим условие антисимметрии (условие Неймана B n = 0). Далее командой Subdomain/ Subdomain Mode переключаемся в режим ввода коэффициентов PDE.

В воздухе, окружающем катушку с сердечником, (зона №1) задаём c=1, f=0, gamma = 0 0 0. В слое обмотки (зона №2) задаём c=1, f=0, gamma = 0 0 30-sqrt(x^2+y^2). В зазоре между сердечником и слоем обмотки (зона №3) задаём c=1, f=0, gamma = 0 0 5. В сердечнике (зона №4) задаём c= 1+(500-1)*0.5/(0.5+sqrt(ux^2+uy^2+(5-uz)^2)), f=0, gamma = 0 0 5*(1+(500-1)*0.5/(0.5+sqrt(ux^2+uy^2+(5-uz)^2))). Теперь можно генерировать конечноэлементную сетку командой Mesh/ Mesh Mode (рис. 5).


Рис. 5. Трёхмерная конечноэлементная сетка

По команде Solve/ Solve Problem начинается процесс решения PDE (7). По умолчанию запускается итерационный нелинейный решатель (в Навигаторе моделей выбран нелинейный режим). Перед запуском решения нужно настроить некоторые параметры итерационного решателя, т.к. по умолчанию итерационный процесс прекращается после 25 нелинейных итераций и получается неправильный результат из-за медленной сходимости вычислений. Командой Solve/ Parameters раскроем диалоговое окно настройки параметров решателя и в закладке Nonlinear установим предельное число итераций 150. При таких настройках время вычисления на компьютере с процессором Athlon 950 составит примерно 2 часа. Распределение скалярного магнитного потенциала в плоскости x=0 показано на рис. 6.


Рис. 6. Распределение скалярного магнитного потенциала, полученное в трёхмерной модели


Рис. 7. Распределение аксиальной составляющей напряжённости магнитного поля, рассчитанное в трёхмерной модели

Распределение осевой составляющей скалярного напряженности магнитного поля показано на рис. 8.

Для сравнения представим распределение исследуемого поля в осевом сечении магнитной системы, полученное в двумерной осесимметричной модели (рис. 8, 9). Известно, что двумерные модели при прочих равных условиях рассчитываются точнее. Сравнение показывает незначительное расхождение результатов.


Рис. 8. Распределение скалярного магнитного потенциала, рассчитанное в двумерной осесимметричной модели


Рис. 9. Распределение аксиальной составляющей напряжённости магнитного поля, рассчитанное в двумерной осесимметричной модели

В оглавление


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

Система Orphus

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