Peculiarities of simulation of magnetic field in electromechanical nodes of magnetic-levitation transport system by the method of secondary sources

Full Text

Введение

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

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

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

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

Обратимся к конструкции индуктора линейного двигателя.

На рис. 1 изображено сечение индуктора типичного линейного двигателя, используемого в качестве тягового в высокоскоростном транспорте (рисунок построен на основании примеров из [1]). Чаще всего в таких двигателях применяется двухслойная обмотка. Как было сказано выше, ввиду наличия пазов индуктор создаёт пространственно-неоднородное магнитное поле, причём величина этого поля значительно возрастает при приближении к вершинам зубцов (рис. 2, 3). Данное обстоятельство приводит к тому, что в окрестности указанных точек при дискретизации расчётной области необходимо значительно сгущать сетку.

Рис. 1. Сечение индуктора линейного двигателя

Рис. 2. Распределение x-компоненты индукции магнитного поля в рабочем зазоре индуктора линейного двигателя

Рис. 3. Распределение y-компоненты индукции магнитного поля в рабочем зазоре индуктора линейного двигателя

В работах [2—4] для расчёта магнитного поля линейных двигателей применяется метод конечных элементов. Особенностью данного метода является необходимость построения дискретной сетки во всей расчётной области, включая окружающее электромагнитную систему пространство (вообще говоря, неограниченное) и воздушные зазоры. В линейных двигателях, применяемых в высокоскоростном транспорте, зазоры занимают области пространства, сопоставимые с объёмами магнитопроводов. Кроме того, ввиду разомкнутости магнитной цепи в пространстве создаются значительные поля рассеяния. Как следствие, при попытке построения сетки в таких областях возникает потребность в весьма значительных объёмах физической памяти компьютера. В связи с этим в данной работе предлагается математическая модель, основанная на концепции вторичных источников [5], согласно которой задача расчёта электромагнитного поля сводится к отысканию плотностей источников, индуцируемых в объёмах и на границах кусочно-однородных сред. Это приводит к тому, что расчётная область не выходит за пределы магнитопровода исследуемой электромагнитной системы.

Метод вторичных источников, в свою очередь, также обладает недостатками, которые проявляются при его численной реализации. В [6] показано, что при численном решении интегрального уравнения относительно плотности вторичных источников, даже при весьма точном решении аппроксимирующей его линейных алгебраических уравнений (СЛАУ), не выполняются условия, вытекающие из законов электромагнитного поля, записанных в интегральной форме (принципа непрерывности магнитного потока и закона полного тока). С целью повышения эффективности метода вторичных источников в [6] предложено применять дополнительные интегральные соотношения, в результате чего получаются переопределённые СЛАУ, псевдорешения которых близки к точному решению интегрального уравнения. Идея, реализованная в [6], получила дальнейшее развитие в работах [7, 8], в которых было предложено применять дополнительные интегральные соотношения, введённые в [6], для каждого элемента разбиения границы ферромагнетика, при этом интегральное уравнение не записывается вовсе. В полученных таким образом дискретных моделях для вторичных источников не задаётся система точек коллокаций, а граничные условия учитываются не в отдельных точках границы, а во всех её точках. Это приводит к повышению точности моделирования магнитного поля по сравнению с остальными методами.

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

Границу магнитопровода i разобьём на $N^{\left(L\right)}$ малых отрезков $\Delta l_k$ ($i=1,2,...,N^{\left(L\right)}$). Каждый такой отрезок охватим замкнутым контуром и запишем для него закон полного тока в интегральной форме. Будем полагать, что искомая функция — плотность токов намагниченности $j(Q,t)$ ― в пределах каждого отрезка разбиения для заданного момента времени t принимает постоянное значение:$j(Q,t)=j_k\left(t\right)$ , $Q\in \Delta l_k$, $k=1,2,...,N^{\left(L\right)}$.

Рассмотрим вывод уравнений относительно величин $j_k\left(t\right)$ на примере некоторого ферромагнитного тела с магнитной проницаемостью ${\mu }^{+}$, расположенного в пространстве с магнитной проницаемостью ${\mu }^{-}$ и заданным переменным магнитным полем ${\overrightarrow{B}}_0(Q,t)$ (рис. 4).

Введём единичные нормальный ${\overrightarrow{n}}_k$ и касательный ${\overrightarrow{\tau }}_k$ векторы к отрезку $\Delta l_k$. Вектор ${\overrightarrow{n}}_k$ направлен во внешность магнетика с магнитной проницаемостью ${\mu }^{+}$. Ось z направлена перпендикулярно плоскости рисунка. Векторы ${\overrightarrow{\tau }}_k$, ${\overrightarrow{n}}_k$ и ${\overrightarrow{e}}_z$ образуют правую тройку векторов (рис. 4).

Рис. 4. Ферромагнитное тело, расположенное в плоскопараллельном магнитном поле

В соответствии с законом полного тока в интегральной форме для замкнутого контура, охватывающего отрезок $\Delta l_k$, можно записать

${\int }_{\Delta l_k}\overrightarrow{H^{+}}(Q,t){\overrightarrow{\tau }}_{k}dl-{\int }_{\Delta l_k}\stackrel{}{{\overrightarrow{H}}^{-}(Q,t)}{\overrightarrow{\tau }}_{k}dl=0$,  $k=1,2,...,N^{\left(L\right)}$. (1)

Векторы ${\overrightarrow{H}}^{+}(Q,t)$ и ${\overrightarrow{H}}^{-}(Q,t)$ представим в виде:

${\overrightarrow{H}}^{+}(Q,t) =\left(\overrightarrow{B_0}\right(Q,t)+{\overrightarrow{B}}_{jk}^{+}(Q,t)+\overrightarrow{B_j}(Q,t\left)\right)/\mu$,  $k=1,2,...,N^{\left(L\right)}$, (2)

${\overrightarrow{H}}^{-}(Q,t) =\left({\overrightarrow{B}}_0\right(Q,t)+{\overrightarrow{B}}_{jk}^{-}(Q,t)+{\overrightarrow{B}}_j(Q,t\left)\right)/{\mu }_0$,  $k=1,2,...,N^{\left(L\right)}$, (3)
где ${\overrightarrow{B}}_{jk}^{+}(Q,t)$ — индукция, обусловленная токами полоски $\Delta l_k$; ${\overrightarrow{B}}_j(Q,t)$ — индукция, обусловленная всеми остальными токами на боковой поверхности ферромагнетика.

Далее, соотношения (2) и (3) подставляются в (1).

При выбранной кусочно-постоянной аппроксимации плотности токов намагниченности, ввиду симметрии, для каждого отрезка $\Delta l_k$ можно получить следующие выражения:

${\int }_{\Delta l_k^{+}}{\overrightarrow{B}}_{jk}^{+}(Q,t){\overrightarrow{\tau }}_{k}dl=-{\int }_{\Delta l_k^{-}}{\overrightarrow{B}}_{jk}^{-}(Q,t){\overrightarrow{\tau }}_{k}dl=\frac{1}{2}\underset{\Delta l_k^{+}+\Delta l_k^{-}}{\oint }{\overrightarrow{B}}_{jk}(Q,t)\overrightarrow{dl}=\frac{{\mu }_0{j}_k\left(t\right)\Delta l_k}{2}$. (4)

Учитывая последнее соотношение, для мгновенных значений плотности токов намагниченности получаем следующую СЛАУ:

$j_k\left(t\right)-\frac{\lambda }{\pi \Delta l_k}\sum _{i=1}^{N^{\left(L\right)}}{j}_i\left(t\right){\int }_{\Delta l_k}{\int }_{\Delta l_i}\frac{({\overrightarrow{r}}_{PQ},{\overrightarrow{n}}_k)}{r_{PQ}^2}d{l}_{p}d{l}_Q=\frac{\lambda }{\Delta l_k{\mu }_0}{\int }_{\Delta l_k}{\overrightarrow{B}}_{0}9Q,t){\overrightarrow{\tau }}_{k}d{l}_Q$,  $k=1,2,...,N^{\left(L\right)}$, (5)
где $\lambda =(\mu -{\mu }_0)/(\mu +{\mu }_0)$.

После решения системы уравнений (5) может быть рассчитана индукция магнитного поля в любой точке  по формуле:

$\overrightarrow{B}(M,t)=\frac{{\mu }_0}{2\pi }\sum _{i=1}^{N^{\left(L\right)}}{j}_i\left(t\right){\int }_{\Delta l_i}\frac{\left[{\overrightarrow{e}}_z,{\overrightarrow{r}}_{PQ}\right]}{r_{PQ}^2}d{l}_p+{\overrightarrow{B}}_0(Q,t)$ . (6)

На основании предложенной математической модели был произведён расчёт магнитного поля индуктора, изображённого на рис. 1. Результаты этих расчётов приведены на рис. 2, 3. Поскольку в [1] геометрические размеры индуктора указаны не были, то на осях абсцисс рисунков отмечены безразмерные величины (где L — длина индуктора). В расчётах были условно заданы следующие параметры индуктора:L=1000 мм, полюсное деление $\tau =98$ мм, ширина зубца — 22 мм, ширина паза — 27 мм, плотность тока в обмотке — 37 МА/м2. Была задана двухслойная обмотка по схеме, приведённой в [1]. Магнитная проницаемость стали индуктора полагалась постоянной и равной $\mu =3500{\mu }_0$.

Перейдём к рассмотрению вопроса построения математической модели магнитного поля электродинамического подвеса.

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

Исследования устройств подвеса, применяемых в высокоскоростном транспорте, показали, что наибольшей электродинамической стабильностью обладают так называемые системы отталкивающего типа [9]. В устройстве подобных систем могут применяться как управляемые электромагниты (в том числе со сверхпроводящими катушками), так и постоянные магниты.

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

Как показано, например, в [10], на границах ферромагнитных проводников, расположенных в плоскопараллельном переменном магнитном поле, отсутствуют электрические заряды. Поэтому плоскопараллельная постановка позволяет уменьшить не только геометрическую размерность задачи, но и количество неизвестных величин, описывающих исследуемый процесс.

Большинство современных научных исследований и разработок систем подвеса наземного транспорта, основанных на использовании постоянных магнитов, предполагает применение массивов Хальбаха. Особенностью подобных систем является несимметричность магнитного потока в области над массивом и под ним, что позволяет сконцентрировать магнитное поле в рабочей области и уменьшить его в зоне пребывания пассажиров. Существуют различные модификации систем электродинамического подвеса, основанных на применении массивов Хальбаха. В [11] предложена система, которая представляет собой массив из постоянных магнитов, в промежутках между которыми расположены полюсные магнитомягкие вставки (ярмо). Массивы с постоянными магнитами устанавливаются на подвижном составе. С целью повышения левитационного качества подвеса массивы часто располагают по обе стороны относительно путевого полотна [12].

Получим уравнения относительно плотностей вторичных источников, индуцируемых в системе подвеса, предложенной в [11].

Будем полагать, что все материалы системы подвеса являются однородными и изотропными. Магнитные свойства стального ярма характеризуются магнитной проницаемостью $\mu =const$. Постоянные магниты изготовлены из высококоэрцитивного магнитного материала с однородно распределённой по всему объёму намагниченностью. Путевое полотно представляет собой массивную проводящую плиту из материала с удельной электрической проводимостью $\gamma =const$.

Неподвижную систему координат (СК) свяжем с путевым полотном. Тогда относительно рассматриваемой СК поле, создаваемое перемещающимся массивом с постоянными магнитами, будет зависеть от времени. Вследствие плоскопараллельности поля плотности вихревых токов и токов намагниченности будут иметь только одну компоненту. Мгновенное значение векторного потенциала магнитного поля в сечении путевого полотна будет иметь вид:

$A(M,t)=A_0(M,t)+\frac{{\mu }_0}{2\pi }{\int }_{S_3}\delta (N,t)\ln \left(1/r_{MN}\right)d{S}_N+\frac{{\mu }_0}{2\pi }\underset{L_2}{\oint }j(P,t)\ln \left(1/r_{MP}\right)d{L}_p+C\left(t\right)$, (7)
где $M\in S_3$ — точка сечения путевого полотна; $A_0(M,t)$ — мгновенное значение векторного потенциала перемещающихся магнитов; $\delta (N,t)$ — мгновенное значение плотности вихревых токов в сечении путевого полотна; $j(P,t)$ — мгновенное значение плотности токов намагниченности на границе ярма; $L_2$ — замкнутый контур, ограничивающий сечение ярма; $r_{MN}$ и $r_{MP}$ — расстояния между точками M и N и точками M и P соответственно; $C\left(t\right)$ — функция, не зависящая от пространственных координат.

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

$\delta (M,t)-\frac{\gamma {\mu }_0}{2\pi }\frac{\partial }{\partial t}\underset{S_3}{\int }\delta (N,t)F(M,N)d{S}_N-\frac{\gamma {\mu }_0}{2\pi }\frac{\partial }{\partial t}\underset{L_2}{\oint j(p,t)F9M,P,t)d{L}_p=W_0(M,t)}$. (8)

В (8) обозначено: $W_0(M,t)=K{w}_0$, $F(M,N)=Kf(M,N)$, $F(M,P,t)=Kf(M,P,t)$; $w_0=-\gamma \partial A_0(M,t)/\partial t$, $f(M,N)=-\ln \left((1/r_{MN}\right)$; K — оператор, действие которого на некоторую скалярную функцию $\psi$ приводит к выражению: $K\psi =\psi -(1/S_3)\underset{S_3}{\int }\psi dS$.

При численном решении уравнения (8) сечение проводящего путевого полотна разбивается на элементы, в пределах которых мгновенные значения плотности вихревых токов считаются неизменными: $\delta (N,t)={\delta }_k\left(t\right)$, $k=1,2,...,N^{\left(S\right)}$ ($N^{\left(S\right)}$ ― число элементов разбиения сечения). Уравнение (8) необходимо дополнить уравнениями относительно плотности токов намагниченности, которые могут быть получены аналогично уравнению (5).

Рассматриваемый период времени разбивается на , в общем случае, неравных интервалов. Мгновенные значения искомых плотностей токов для каждого интервала времени обозначим следующим образом: ${\delta }_k\left(t_i\right)={\delta }_{k,i}$ ($k=1,2,...,N^{\left(S\right)}$; $i=0,1,...,N_T$), $j_k\left(t_i\right)=j_{k,j}$ ($k=1,2,...,N^{\left(L\right)}$; $i=0,1,///,N_T$). Соответствующая итерационная схема примет вид:                

$j_{k,i}^{\left(s\right)}-\frac{\lambda }{\pi \Delta l_k}\sum _{m=1}^{N^{\left(L\right)}}{j}_{m,i}^{\left(s\right)}\underset{\Delta l_k}{\int }\underset{\Delta l_m}{\int }\frac{({\overrightarrow{r}}_{PQ},{\overrightarrow{n}}_k)}{r_{PQ}^2}d{l}_{p}d{l}_{Q=}$ $\frac{\lambda }{\pi {\mu }_0\Delta l_k}\sum _{m=1}^{N^{\left(S\right)}}{\delta }_{m,i}^{(s-1)}\underset{\Delta l_k}{\int }\underset{\Delta S_m}{\int }\frac{({\overrightarrow{r}}_{NQ},{\overrightarrow{n}}_k)}{r_{NQ}^2}d{S}_N+\frac{\lambda }{\Delta l_k{\mu }_0}\underset{\Delta l_k}{\int }{\overrightarrow{B}}_0(Q,t){\overrightarrow{\tau }}_{k}d{l}_Q$,  $k=1,2,...,N^{\left(L\right)}$;  (9) ${\tilde{\delta }}_{k,i}^{\left(S\right)}=\frac{\gamma {\mu }_0}{2\pi }\sum _{m=1}^{N^{\left(S\right)}}\frac{\partial {\delta }_m^{(s-1)}\left(t_i\right)}{\partial t}\underset{\Delta S_m}{\int }F(M_k,N)d{S}_N$$+\frac{\gamma {\mu }_0}{2\pi }\frac{\partial }{\partial t}\sum _{m=1}^{N^{\left(L\right)}}{j}_m^{\left(s\right)}\left(t_i\right)\underset{L_2}{\oint F(M,P,t)d{L}_p+W_0\left(M_{k,t_i}\right)}$ . (10)

В (9), (10) индекс s в скобках означает номер итерации.

Порядок выполнения итерационной процедуры следующий:

1. Задаётся начальное приближение для плотности вихревых токов${\delta }_{k,i}^{\left(0\right)}$ , равное ${\delta }_{k,i}^{\left(0\right)}=W_0(M_k,t_i)$, где $k=1,2,...,N^{\left(S\right)}$; $0,1,...,N_T$.
2. Для каждого момента времени решается СЛАУ (9), в результате чего определяется первое приближение для плотности токов намагниченности $j_{k,i}^{\left(1\right)}$($k=1,2,...,N^{\left(L\right)}$, $i=0,1,...,N_T$).
3. Имеющиеся значения ${\delta }_{k,i}^{\left(1\right)}$ и $j_{k,i}^{\left(1\right)}$ подставляются в правую часть формулы (10). После осуществления всех вычислений определяется функция ${\delta }_{k,i}^{\left(1\right)}$, $k=1,2,...,N^{\left(S\right)}$, $i=0,1,...,N_T$.
4. По формуле ${\delta }_{k,i}^{\left(1\right)}={\delta }_{k,i}^{\left(0\right)}+\beta ({\delta }_{k,i}^{\left(1\right)}-{\delta }_{k,i}^{\left(0\right)})$ вычисляется новое приближение для ${\delta }_{k,i}^{\left(1\right)}$ ($k=1,2,...,N^{\left(S\right)}$, $i=0,1,...,N_T$), где $\beta$ выбирается в интервале (0;1) исходя из численного эксперимента. Подбор оптимального значения коэффициента $\beta$ обеспечивает монотонную сходимость процесса, даже если при $\beta =1$ этот процесс расходился.

На этом заканчивается первая итерация. Далее процесс переходит к пункту 2 и совершается вторая итерация и т.д.

5. Итерационный процесс прерывается после выполнения условия:$100·||{\delta }^{\left(s\right)}(M,t)-{\delta }^{(s-1)}(M,t)||/||{\tilde{\delta }}^{\left(s\right)}(M,t)||\le \epsilon$ , где $\epsilon$ — заданное малое число; $||{\delta }^{\left(s\right)}(M,t)||=\sqrt{\underset{S_r}{\int \underset{0}{\overset{T}{\int }}}\left({\delta }^{\left(s\right)}\right(M,t){)}^{2}dtdS}$.

Производные, входящие в (10), вычисляются численно по следующим формулам:

$\partial {\delta }_m^{(s-1)}\left(t_0\right)/\partial t=(-3{\delta }_{m,0}^{(s-1)}+4{\delta }_{m,1}^{(s-1)}-{\delta }_{m,2}^{(s-1)})/\left(2\Delta T\right)$ (при $i=0$);

$\partial {\delta }_m^{(s-1)}\left(t_i\right)/\partial t=({\delta }_{m,i+1}^{(s-1)}-{\delta }_{m,i-1}^{(s-1)})/\left(2\Delta T\right)$ (при $0<i<N_T$);

$\partial {\delta }_m^{(s-1)}\left(t_0\right)/\partial t=({\delta }_{m,N_T-2}^{(s-1)}-4{\delta }_{m,N_T-1}^{(s-1)}+3{\delta }_{m,N_T}^{(s-1)})/\left(2\Delta T\right)$ (при $i=N_T$).

Необходимо также учитывать, что во втором интегральном члене (10) от времени зависит не только искомая функция , но и ядро соответствующего интеграла, а также пределы интегрирования ввиду того, что ярмо перемещается вместе с массивом магнитов относительно путевого полотна.

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

Силовые характеристики электродинамического подвеса вычисляются по формуле ампера.

Заключение

1. Ввиду особенности конструкции индуктора тягового линейного двигателя и характера распределения магнитного поля в его рабочем зазоре с целью снижения размерности расчётной области целесообразно применять метод вторичных источников. Дискретные модели для вторичных источников предлагается строить на основании записи закона полного тока для контуров, охватывающих каждый элемент разбиения границы магнитопровода. Данный метод также эффективен при расчёте магнитного поля электродинамического подвеса.
2. Разработанная дискретная модель может быть применена к моделированию магнитного поля линейного двигателя при учёте насыщения участков магнитопровода. Для этого достаточно применить методику расчёта нелинейности, разработанную в [13].
3. Разработанная математическая модель магнитного поля электродинамического подвеса позволяет производить расчёты как установившегося режима работы подвеса (при постоянной скорости движения транспортного средства), так и динамического (при ускоренном движении транспортного средства). Соответствующий режим определяется заданным законом механического движения массива постоянных магнитов и ярма.

About the authors

Dmitry M. Filippov

Author for correspondence.
Email: grantarc@mail.ru

Ph.D., a senior pre-feeder of the Department of Computer Engineering and Modeling

References

Supplementary files