Метод конечных элементов и COMSOL Multiphysics · ЭМЭЭА

1. Зачем полевой расчёт в ВКР

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

распределение магнитной индукции B в воздушном зазоре и стали;
токи Фуко и потери на вихревые токи в массивных проводящих частях;
силу электродинамического взаимодействия токопроводов при коротком замыкании;
распределение напряжённости электрического поля E и локальное усиление поля в изоляции; тангенс угла диэлектрических потерь задают свойствами материала, а частичные разряды требуют отдельной модели или экспериментальной диагностики;
температурное поле в обмотках и сердечнике;
электромагнитные силы, которые после передачи в механическую модель используют для расчёта деформаций, собственных частот и вибрации.

Что приветствует ГЭК

Аналитическая оценка + численная модель + проверка результата. Допустимое расхождение задают для конкретной величины с учётом допущений, неопределённости исходных данных и требуемой точности. МКЭ не освобождает от баланса размерностей.

2. Уравнения Максвелла

Всё электромагнитное моделирование сводится к решению системы Максвелла. В дифференциальной форме:

\nabla \cdot \mathbf{D} = \rho

(1)

\nabla \cdot \mathbf{B} = 0

(2)

\nabla \times \mathbf{E} = -\frac{\partial \mathbf{B}}{\partial t}

(3)

\nabla \times \mathbf{H} = \mathbf{J} + \frac{\partial \mathbf{D}}{\partial t}

(4)

где $\mathbf{E}$ — напряжённость электрического поля, В/м;
$\mathbf{D} = \varepsilon\mathbf{E}$ — электрическая индукция, Кл/м²;
$\mathbf{B}$ — магнитная индукция, $\text{Тл}$;
$\mathbf{H} = \mathbf{B}/\mu$ — напряжённость магнитного поля, А/м;
$\mathbf{J} = \sigma\mathbf{E}$ — плотность тока проводимости, А/м²;
$\rho$ — плотность заряда, Кл/м³.

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

3. Потенциалы и формулировки

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

3.1. Магнитный векторный потенциал

Из $\nabla\cdot\mathbf{B}=0$ следует, что существует векторный потенциал $\mathbf{A}$:

\mathbf{B} = \nabla \times \mathbf{A}

(5)

Подставляя в уравнение Ампера, получаем основное уравнение магнитостатики:

\nabla \times \left( \frac{1}{\mu}\, \nabla \times \mathbf{A} \right) = \mathbf{J}

(6)

В плоско-параллельной двумерной задаче, инвариантной вдоль оси $z$, принимают $\mathbf{A}=(0,0,A_z)$ и решают скалярное уравнение для $A_z$. Осесимметричная постановка в координатах $r$–$z$ является отдельной формулировкой.

3.2. Электрический скалярный потенциал

В электростатике из $\nabla\times\mathbf{E}=0$ следует $\mathbf{E}=-\nabla\varphi$, а из закона Гаусса получаем уравнение Пуассона:

\nabla \cdot ( \varepsilon\, \nabla \varphi ) = -\rho

(7)

При $\rho = 0$ — уравнение Лапласа $\Delta\varphi = 0$. В COMSOL это формулировка Electrostatics (es).

3.3. Гармонический режим (AC)

Для синусоидального тока частоты $f$ переходим в комплексную область: $\partial/\partial t \to j\omega$, где $\omega = 2\pi f$.

\nabla \times \left( \frac{1}{\mu}\, \nabla \times \tilde{\mathbf{A}} \right) + j\omega\sigma\,\tilde{\mathbf{A}} = \tilde{\mathbf{J}}_s

(8)

Это упрощённая $\mathbf A$-формулировка для гармонического квазистатического режима при выбранной калибровке и без явно записанного скалярного потенциала. В COMSOL состав уравнений интерфейса Magnetic Fields (mf) зависит от областей, источников и настроек проводимости.

4. МКЭ: дискретизация

Метод конечных элементов превращает дифференциальное уравнение в систему алгебраических уравнений. Кратко:

Слабая формулировка. Умножаем уравнение на тестовую функцию v и интегрируем по объёму, перенося часть производных на тестовую функцию (формула Грина).
Дискретизация. Расчётная область разбивается на конечные элементы — треугольники / тетраэдры (P1, P2) или четырёхугольники / гексаэдры. Внутри каждого элемента решение аппроксимируется полиномом малой степени.
Сборка. Локальные матрицы собираются в глобальную систему $\mathbf{K}\mathbf{u}=\mathbf{f}$, где $\mathbf{u}$ — вектор узловых значений.
Граничные условия. Дирихле (значение на границе) — корректировка строк матрицы. Нейман (поток на границе) — добавляются в правую часть.
Решение. Прямой или итерационный решатель выбирают по размеру, обусловленности, нелинейности и доступной памяти; пригодность подтверждают сходимостью и устойчивостью интересующей величины.

Точность МКЭ зависит от аппроксимации геометрии и материалов, размера и порядка элементов, корректности граничных условий, нелинейного и линейного решателей. Контролируют не только норму невязки, но прежде всего сходимость целевой величины: потока, потерь, силы, момента или температуры.

5. Типы полевых задач, актуальных для ВКР ЭМЭЭА

Типовые полевые задачи для ВКР ЭМЭЭА
Задача	Что решаем	Модуль COMSOL
Электростатика	Распределение E в изоляции аппаратов, ёмкости	AC/DC · Electrostatics (es)
Магнитостатика	B в зазоре машины, реакторе, экранах	AC/DC · Magnetic Fields (mf, stationary)
Гармоническое поле	Вихревые токи, индуктивность, потери в стали	AC/DC · Magnetic Fields (mf, frequency domain)
Переходные процессы	Включение трансформатора, КЗ, удар тока	AC/DC · Magnetic Fields (mf, time-dependent)
Электромеханика	Расчёт сил, моментов, движущийся ротор	AC/DC · Rotating Machinery, Magnetic (rmm)
Тепловой режим	Нагрев обмоток, сердечника	Heat Transfer (ht) + Joule heating
Связанная задача	Электромагнетизм + тепло + механика	Multiphysics: Electromagnetic Heating

6. Порядок работы в COMSOL Multiphysics

Типовая последовательность для задачи магнитостатики (например, индуктор / реактор):

Шаг 1. Выбор пространства и физики

File → New → Model Wizard
  Space Dimension: 2D Axisymmetric  (для осесимметричной катушки)
                   или 2D            (плоско-параллельная задача)
                   или 3D            (несимметричная геометрия)
  Physics → AC/DC → Magnetic Fields (mf)
  Study → Stationary  (или Frequency Domain / Time Dependent)

Шаг 2. Геометрия

Геометрию создают в узле Geometry или импортируют из CAD. Внешний воздушный домен выбирают по исследованию чувствительности к положению границы либо заменяют бесконечными элементами. Универсального коэффициента увеличения области нет: он зависит от геометрии, источников и граничного условия.

Шаг 3. Материалы

Из библиотеки COMSOL: Materials → Add Material from Library → AC/DC → Soft Iron (with HB-curve). Для собственных материалов:

Materials → Add Material → Blank Material
  Basic property → Relative permeability:        μ_r = ...
                   Electrical conductivity:       σ = ... [S/m]
                   Density:                       ρ = ... [kg/m³]
Magnetic field → B-H curve from file (interpolation) для нелинейной стали

Числа из примера интерфейса нельзя использовать как свойства реального материала. Для итогового расчёта вводят измеренную характеристику $B(H)$ либо данные изготовителя с указанием марки, температуры, частоты и направления прокатки, если они значимы.

Шаг 4. Граничные условия и источники

Magnetic Fields (mf):
  Ampère's Law (default) — основное уравнение во всём объёме
  Magnetic Insulation (default) — на внешней границе (B·n = 0)
  + Coil — задаёт ток в катушке
    Conductor model: Single conductor / Homogenized multi-turn
    Coil excitation: Current  I = 100 [A]
    Number of turns: 250

В гармонической задаче проводимость, нелинейная характеристика $B(H)$, шихтование и модель магнитных потерь задаются раздельно. Наличие кривой $B(H)$ само по себе не включает вихревые токи и не определяет потери в стали.

Шаг 5. Сетка (Mesh)

Mesh:
  Element size: Normal / Fine / Extra Fine — общее качество
  + Free Triangular (для 2D)
  + Boundary Layers — на поверхности проводников при скин-эффекте
  + Size — локальное измельчение в зазоре и углах

Для проводника полезно сначала оценить характерную глубину проникновения:

\delta = \sqrt{\frac{2}{\omega \mu \sigma}}

(9)

где $\delta$ — глубина скин-слоя при заданной угловой частоте $\omega$. Необходимое число элементов на глубине $\delta$ определяют исследованием сходимости. Для ферромагнитной стали эта оценка является ориентировочной, поскольку $\mu$ зависит от поля и частоты.

Шаг 6. Решение (Study)

Study 1 → Step 1: Stationary
  Solver Configurations → Stationary Solver 1
    Relative tolerance: обосновать по сходимости целевой величины
    Nonlinear method: подобрать для заданной B-H характеристики
  Compute

Время и устойчивость расчёта зависят не только от числа степеней свободы, но и от обусловленности, нелинейности, типа решателя и вычислительных ресурсов. Предварительную сетку используют для отладки постановки, а итоговую выбирают по исследованию сеточной сходимости.

Шаг 7. Постобработка (Results)

В разделе Results доступны:

2D / 3D Plot Group — Surface plot (магнитная индукция mf.normB), Contour, Arrow, Streamline;
1D Plot Group — графики по линии или дуге (для построения распределения B вдоль зазора);
Derived Values — интегральные величины:
- Surface Integration → магнитный поток через сечение;
- Volume Integration → интегрирование плотности потерь или энергии; имя и размерность переменной проверяют в текущем физическом интерфейсе;
- Force Calculation → электродинамическая сила.
Export — графики в PNG/EPS/PDF, таблицы в CSV/Excel для ВКР.

7. Практические примеры с формулами

7.1. Индуктивность катушки через энергию

Для линейной магнитной системы с одной токовой координатой индуктивность можно вычислить через магнитную энергию $W_m$:

L = \dfrac{2\,W_{\text{m}}}{I^2}, \qquad W_{\text{m}}=\frac{1}{2}\int_V \mathbf{B} \cdot \mathbf{H}\,\mathrm{d}V

(10)

Равенство справедливо для линейной среды. В нелинейной системе различают потокосцепление, статическую и дифференциальную индуктивности; способ расчёта выбирают по задаче. Перед интегрированием проверьте, является ли выбранная переменная плотностью энергии или уже интегральной величиной.

7.2. Сила через тензор натяжений Максвелла

Электродинамическую силу можно вычислить интегрированием тензора натяжений Максвелла по замкнутой поверхности в однородной немагнитной области вокруг объекта:

\mathbf{F} = \oint_S \left[ \mathbf{H}(\mathbf{B}\cdot\mathbf{n}) - \tfrac{1}{2}\,\mathbf{n}(\mathbf{H}\cdot\mathbf{B}) \right]\,\mathrm{d}S

(11)

В COMSOL: Definitions → Component Couplings → Integration on a surface, либо через готовый Force Calculation feature.

7.3. Потери в стали (модель Бертотти)

Расширенная формула удельных потерь в шихтованной стали включает три составляющие: гистерезисные, вихревые и аномальные:

p_{\text{ст}} = k_{\text{h}}\, f\, B_{\text{m}}^{\alpha} + k_{\text{c}}\, (f B_{\text{m}})^2 + k_{\text{e}}\, (f B_{\text{m}})^{1{,}5}

(12)

где $k_{h}$, $k_{c}$, $k_{e}$ — коэффициенты, идентифицированные по данным измерений для применимого диапазона частоты и индукции;
$\alpha$ — показатель Штейнмеца, обычно $1{,}6 \ldots 2{,}2$;
$B_{m}$ — амплитуда индукции, $\text{Тл}$;
$f$ — частота, $\text{Гц}$.

7.4. Скин-эффект в массивном проводнике

Глубина проникновения переменного тока в проводящую среду:

\delta = \sqrt{\frac{2}{\omega \mu \sigma}} = \sqrt{\frac{1}{\pi f \mu \sigma}}

(13)

Для плоской поверхности однородного проводящего полупространства сравните результат МКЭ-расчёта с оценкой $J(x)\approx J_s e^{-(1+j)x/\delta}$. Для проводника конечного сечения используют соответствующее аналитическое решение или исследование сходимости.

7.5. Индукция в зазоре электрической машины

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

B_{\delta} \approx \dfrac{N I}{A_{\delta}\,\mathcal{R}_{\Sigma}}, \qquad \mathcal{R}_{\Sigma}=\sum_i \dfrac{l_i}{\mu_i A_i}

(14)

где $\mathcal{R}_{\Sigma}$ — суммарное магнитное сопротивление принятой эквивалентной цепи; $l_i$, $\mu_i$, $A_i$ — длина, магнитная проницаемость и площадь соответствующего участка. Коэффициент Картера, рассеяние и насыщение учитывают отдельно.

8. Качество сетки

Качество сетки оценивают не по одному универсальному порогу, а по совокупности признаков:

геометрическое качество элементов — контролируют по метрике, определённой в используемой программе, особенно в узких зазорах и у острых кромок;
разрешение характерных масштабов — сетка должна описывать зазор, скин-слой, кривизну и зоны большого градиента;
сеточная сходимость — последовательно уточняют сетку и фиксируют изменение целевой величины;
баланс точности и стоимости — критерий прекращения уточнения связывают с требуемой неопределённостью результата.

В отчёте приводят определение использованной метрики качества, число степеней свободы, способ уточнения и таблицу сходимости для интересующей величины. Большая сетка — ещё не большая точность: зачёт получает сходящийся результат, а не число элементов.

9. Граничные условия — что выбрать

Граничные условия в COMSOL
Условие	Когда применять	Что означает физически
Magnetic Insulation	Плоскость симметрии или удалённая внешняя граница после проверки чувствительности	$\mathbf{B}\cdot\mathbf{n}=0$: нормальная составляющая магнитной индукции равна нулю.
Perfect Magnetic Conductor	Плоскость симметрии нормальная к $\mathbf{B}$	$\mathbf{n}\times\mathbf{H}=0$: зеркальное отражение поля.
Periodic Condition	Периодический сектор машины	Поля на парных границах связаны периодическим или антипериодическим условием в соответствии с геометрической и электромагнитной симметрией.
Surface Current Density	Тонкая катушка / лента	$\mathbf{n}\times\mathbf{H}=\mathbf{J}_s$: ток задаётся без модели проводника.

10. Расчёт потерь в COMSOL

В частотной задаче набор доступных величин зависит от физического интерфейса, материала и включённых моделей потерь. Перед использованием проверьте определение и размерность переменной. Обычно анализируют:

mf.Qrh — удельные потери на джоулев нагрев, Вт/м³;
mf.Wm_avg — средняя плотность магнитной энергии, Дж/м³;
интеграл по объёму даёт полные потери (Вт) и полную энергию (Дж).

Для потерь в шихтованной стали:

Magnetic Fields → материал и закон намагничивания
  B-H characteristic: данные для выбранной марки материала
  Magnetic losses: модель, согласованная с видом возбуждения
  Loss data: измерения или технические данные изготовителя
             для применимых частот, индукций и температуры

Альтернативно — постпроцессорный расчёт по формуле Бертотти на основе амплитуды $B_{\text{m}}$.

11. Силы и моменты

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

Тензор натяжений Максвелла — интегрирование по замкнутой поверхности; результат чувствителен к положению поверхности и сетке, поэтому требует проверки сходимости.
Метод виртуальных перемещений — производная энергии или коэнергии по обобщённой координате; численная реализация и требования к сетке зависят от постановки.

В COMSOL предпочтительно использовать предусмотренный выбранным интерфейсом узел Force Calculation. Имена переменных зависят от версии и постановки; их физический смысл, систему координат и единицы необходимо проверить в документации модели.

12. Связанные задачи — электромагнетизм + тепло

Самая частая связь для ВКР — нагрев обмоток джоулевыми потерями.

Add Physics → Heat Transfer in Solids (ht)
Multiphysics → Electromagnetic Heating (emh1)
  Источник тепла: mf.Qrh (автоматически)
Study → Stationary (если установившийся режим) или Time Dependent

Граничные условия теплопередачи:

Convective Cooling на поверхности корпуса: коэффициент теплоотдачи задают по измерениям, критериальной зависимости или проверенному источнику для конкретной поверхности и режима течения;
Surface-to-Ambient Radiation, когда лучистый теплообмен значим для температур и степени черноты рассматриваемых поверхностей;
Temperature — заданная температура (например, обмотка статора как граница).

13. Альтернативы COMSOL

Если в вузе нет лицензии COMSOL, для ВКР подойдут:

ANSYS Maxwell (Electronics Desktop) — коммерческая среда для электромагнитного моделирования; доступность учебной лицензии проверяйте на момент выполнения работы;
FEMM (David Meeker) — программа для двумерных и осесимметричных задач; область применимости определяется доступными физическими постановками;
JMAG-Designer / JMAG-Express — специализированные средства расчёта электрических машин; условия лицензирования необходимо проверить у разработчика;
Elmer FEM — программное обеспечение с открытым исходным кодом, использующее ElmerGUI или прямую настройку файла модели;
облачные расчётные среды — применимы при наличии разрешения на передачу модели и исходных данных внешнему сервису;
Python + scikit-fem / FEniCS / NGSolve — для ВКР с алгоритмической направленностью и воспроизводимым исходным кодом.

14. Оформление полевых расчётов в ВКР

Что должно быть в разделе «Полевой расчёт» ВКР:

Постановка задачи — какое поле решаем, в какой области, какие физические упрощения сделаны.
Расчётная область — размеры воздушного домена, плоскости симметрии. Рисунок 2D-сечения с размерами по ГОСТ 2.307-2011 → см. Рисунки и таблицы.
Уравнения и формулировка — приведите уравнение (5)–(8) с пояснением символов под формулой (см. Формулы). Не цитируйте всю систему Максвелла — только используемое.
Граничные условия — таблица: какие границы, какие условия, обоснование.
Свойства материалов — таблица с $\mu$, $\sigma$, $\rho$, $\lambda$ и ссылкой на технические данные изготовителя или результаты измерений.
Сетка — число степеней свободы, тип элементов, определение показателя качества и изображение сетки в зоне интереса.
Сходимость — таблица зависимости целевой величины от числа степеней свободы или характерного размера элемента.
Результаты — цветовые карты B/E/T, графики B(x) вдоль линий, численные значения индуктивности / потерь / температуры.
Сравнение с аналитикой — расхождение в процентах, объяснение.
Визуализации результатов — со шкалой, легендой, единицами и подписью.

Все формулы в ВКР должны быть набраны по правилам оформления формул: переменные курсивом, операторы прямым шрифтом, единицы измерения после числа с неразрывным пробелом.

15. Частые ошибки в полевых расчётах ВКР

Воздушный домен слишком мал. Внешняя граница влияет на целевую величину. Решение: изменять размер домена до достижения принятого критерия чувствительности либо применять бесконечные элементы.
Линейная сталь в нелинейной задаче. Если рабочая точка приближается к области насыщения выбранного материала, используйте его характеристику $B(H)$ и проверьте влияние нелинейности на результат.
Слишком грубая сетка в зазоре. Число элементов поперёк зазора выбирайте по сходимости распределения индукции, силы или момента.
Скин-слой не разрешён. Если рассчитывается распределение переменного тока в массивном проводнике, уточняйте сетку по глубине проникновения и подтверждайте сходимость потерь или полного сопротивления.
Двумерная постановка без обоснования. Оцените влияние концевых частей, скоса, лобовых соединений и трёхмерного рассеяния; при необходимости примените корректирующую модель или трёхмерный расчёт.
Отсутствие исследования сходимости. Одного запуска недостаточно. Число уровней уточнения должно позволять убедительно показать устойчивость целевой величины.
Нет независимой проверки. Без верификации модель остаётся «чёрным ящиком», даже если цветовая карта выглядит убедительно.
Изображения без шкалы. Цветовая карта без легенды, единиц и диапазона не позволяет интерпретировать результат.
Версия программы не указана. В тексте ВКР укажите версию COMSOL (например, 6.2), модуль (AC/DC) и тип лицензии (студенческая / коммерческая).

Связь с другими разделами

— Полевой расчёт обычно — это главы 3 и 4 ВКР (модель и эксперимент); см. Структура ВКР.
— Все формулы оформляйте по правилам типографики формул.
— Рисунки результатов — по ГОСТ Р 2.105-2019, см. Оформление рисунков и таблиц.
— Список ГОСТов на полевые расчёты → Все ГОСТы для ВКР.
— Расчёты как глава ВКР → Расчётная часть и Моделирование.

16. Источники и дополнительная литература

COMSOL AB. AC/DC Module User's Guide — официальная документация.
Демирчян К.С., Чечурин В.Л. Машинные расчёты электромагнитных полей. — Москва : Высшая школа, 1986.
Иванов-Смоленский А.В. Электромагнитные поля и процессы в электрических машинах. — Москва : Энергия, 1969.
Jin J.-M. The Finite Element Method in Electromagnetics. 3rd ed. — Hoboken : Wiley-IEEE Press, 2014. — 876 p.
Bianchi N. Electrical Machine Analysis Using Finite Elements. — Boca Raton : CRC Press, 2017. — 304 p.
Meeker D. Finite Element Method Magnetics (FEMM). — femm.info
ГОСТ Р 2.105-2019. Общие требования к текстовым документам — для оформления раздела «Полевой расчёт» в ВКР.