Расчёт вихревых токов в общем случае

Перейдём к выводу уравнений, описывающих вихревые токи в проводящих телах под действием сторонних (внешних) магнитных полей.

Оглавление
Влияние грунта на металлоискатель. Отклик на предметы, не являющиеся малыми (основной документ)

Задача о расчёте вихревых токов
Расчёт вихревых токов в общем случае
Решение задачи при наличии осевой симметрии
Численное решение задачи о вихревых токах при наличии осевой симметрии
Решение задачи в Octave. Примеры расчётов
Решение задачи в Octave. Влияние грунта
Решение задачи в Octave. Приложение: справочник по функциям

Смотрите также
Металлоискатели (общие вопросы)

Предположим, что имеется некоторый источник переменного магнитного поля, например катушка L₁, через которую протекает известный переменный ток. Магнитное поле катушки также считаем известным: в любой точке пространства X имеем возможность вычислить векторный потенциал $ \vec A_0 $ поля, создаваемого катушкой, рис. %img:f.

Катушка L₁ находится в среде 1 (воздух), параметры которой для решения нашей задачи можно считать такими же, как у вакуума $ \mu_1=1, \sigma_1=0, \text{ где } \sigma_1 \ -$ удельная проводимость первой среды.

В магнитном поле катушки находится тело 2 из проводящего материала (метал, грунт, морская вода или др.) с отличной от нуля удельной проводимостью $ \sigma_2 $. Поскольку в нашей задаче только проводимость второй среды ненулевая, далее её будем обозначать просто как $\sigma_2=\sigma $.

Будем считать, что также как для первой, для второй среды $ \mu_2=1$, т.е. рассматриваем случай, когда материал проводящего тела не обладает выраженными магнитными свойствами. Тогда получаем $ \mu_1=\mu_2=\mu=1$, следовательно, на границе раздела двух сред не происходит преломления линий векторов магнитного поля $ \vec H, \vec B$. Этот факт несколько упрощает нашу задачу.

Под действием переменного магнитного поля в проводящем теле возникают вихревые токи. Вихревые токи создают своё магнитное поле, которое обозначим как $\vec A_e$ (индекс "e" - сокращение от "eddy currents"). Считая каждую из сред 1, 2 линейной, однородной, изотропной и с учётом того, что, как было сказано ранее, на границе сред в данном случае не происходит преломления линий магнитного поля, магнитное поле в любой точке пространства будет складываться из поля катушки и поля, создаваемого вихревыми токами $$ \vec A = \vec A_0 + \vec A_e. $$ В том числе, это справедливо и для области 2. То есть, в вихревые токи вносит вклад не только поле катушки-индуктора, но и магнитное поле самих вихревых токов.

Для решения задачи воспользуемся уравнениями Максвелла. $$ \def\op{\operatorname} \begin{equation} \op{rot} \vec H = \vec j + \frac {\partial \vec D} {\partial t}, \label{mx1} \end{equation} $$ $$ \begin{equation} \op{div} \vec B = 0, \label{mx2} \end{equation} $$ $$ \begin{equation} \op{rot} \vec E = - \frac {\partial \vec B} {\partial t}, \label{mx3} \end{equation} $$ $$ \begin{equation} \op{div} \vec D = \rho. \label{mx4} \end{equation} $$

Дополним эту систему материальными уравнениями $$ \begin{equation} \vec B = \mu \mu_0 \vec H, \label{mx5} \end{equation} $$ $$ \begin{equation} \vec D = \varepsilon \varepsilon_0 \vec E, \label{mx6} \end{equation} $$ $$ \begin{equation} \vec j = \sigma \vec E. \label{mx7} \end{equation} $$

В уравнениях \eqref{mx1}..\eqref{mx7} используются обозначения:
$\vec H$ - вектор напряжённости магнитного поля;
$\vec B$ - вектор индукции магнитного поля;
$\vec E$ - вектор напряжённости электрического поля;
$\vec D$ - вектор электрической индукции;
$\varepsilon$ - относительная диэлектрическая проницаемость;
$\varepsilon_0$ - электрическая постоянная, $$ \varepsilon_0 = 8.8541878128(13) \cdot 10^{−12} \text{ Ф/м} $$ (смотрите изменения СИ 2018..2019 годов);
$\mu$ - относительная магнитная проницаемость;
$\mu_0$ - магнитная постоянная, $\mu_0=1.25663706212(19) \cdot 10^{-6} \text{ Гн/м}$ (смотрите изменения СИ 2018..2019 годов);
$\vec j$ - вектор плотности тока;
$\sigma$ - удельная проводимость среды.

Считая процессы в рассматриваемой системе квазистационарными (достаточно медленно изменяющиеся токи и поля), будем пренебрегать токами смещения (считаем, что $\partial \vec D / \partial t \approx 0$). Это допустимо, поскольку магнитное поле, создаваемое токами смещения, в данном случае крайне слабое по сравнению с магнитным полем катушки L₁ или магнитным полем вихревых токов. Тогда уравнение \eqref{mx1} примет вид $$ \begin{equation} \op{rot} \vec H \approx \vec j, \label{mx1a} \end{equation} $$ а с учётом \eqref{mx7} можем записать $$ \begin{equation} \op{rot} \vec H \approx \sigma \vec E. \label{mx1b} \end{equation} $$

Нашу задачу будет проще решить, если перейти к описанию магнитного поля с помощью векторного потенциала. По определению, векторный потенциал - это такой вектор $\vec A$, что $$ \vec B = \op{rot} \vec A, $$ откуда $$ \begin{equation} \vec H = \frac 1 {\mu \mu_0} \op{rot} \vec A. \label{h_from_a} \end{equation} $$

Подставим B, выраженное через A, в уравнение \eqref{mx3}: $$ \op{rot} \vec E = - \frac {\partial \vec B} {\partial t}, \\ \vec B = \op{rot} \vec A, \\ \op{rot} \vec E = - \frac {\partial} {\partial t} \op{rot} \vec A. $$ В правой части выражения можем поменять местами операторы $\partial / \partial t \text{ и } \op{rot}$ как операторы дифференцирования по независимым переменным (по времени и по координатам). Получим $$ \op{rot} \vec E = - \op{rot} \frac {\partial \vec A} {\partial t}, \\ \op{rot} \vec E + \op{rot} \frac {\partial \vec A} {\partial t} = 0, \\ \begin{equation} \op{rot} \left( \vec E + \frac {\partial \vec A} {\partial t} \right) = 0. \label{rtsz} \end{equation} $$

Ротор градиента любой скалярной функции равен 0, поэтому можем считать, что выражение в скобках в \eqref{rtsz} есть градиент некоторой функции, что запишем в следующем виде $$ \vec E + \frac {\partial \vec A} {\partial t} + \op{grad} \phi =0. $$ В нашем случае поле E создаётся исключительно переменным магнитным полем и поэтому имеет только вихревую составляющую; потенциальная составляющая отсутствует. Так что $$ \op{grad} \phi=0 $$ и $$ \begin{equation} \vec E = - \frac {\partial \vec A} {\partial t}. \label{e_from_a} \end{equation} $$

Подставим в \eqref{mx1b} выражения для H и E из \eqref{h_from_a} и \eqref{e_from_a} соответственно $$ \op{rot} \left( \frac 1 {\mu \mu_0} \op{rot} \vec A \right ) \approx - \sigma \frac {\partial \vec A} {\partial t}, \\ \begin{equation} \op{rot} \op{rot} \vec A \approx - \sigma \mu \mu_0 \frac {\partial \vec A} {\partial t}. \label{rra1} \end{equation} $$ Как мы помним, $\vec A = \vec A_0 + \vec A_e$, поэтому левую часть \eqref{rra1} можно преобразовать следующим образом $$ \op{rot} \op{rot} \vec A = \op{rot} \op{rot} \vec A_0 + \op{rot} \op{rot} \vec A_e, $$ но $$ \op{rot} \op{rot} \vec A_0 = \op{rot} \vec B_0 = \frac 1 {\mu \mu_0} \op{rot} \vec H_0 \approx 0, $$ всюду, кроме точек внутри провода катушки L₁, которые мы не рассматриваем (применять уравнения Максвелла и следствия из них к точкам катушки не требуется, так как поле катушки нам и так уже известно). Последнее равенство следует из \eqref{mx1a}, в соответствии с которым в любой точке $$ \op{rot} \vec H_0 \approx \vec j_0, $$ т.е. ротор напряжённости магнитного поля, создаваемого катушкой, в любой точке пространства равен плотности тока в проводе катушки в данной точке, которая, очевидно, равна нулю всюду вне провода катушки. Тогда получаем, что $\op{rot}\op{rot}\vec A_0 \approx 0 $ и $$ \op{rot} \op{rot} \vec A \approx \op{rot} \op{rot} \vec A_e, $$ и уравнение \eqref{rra1} принимает вид $$ \begin{equation} \op{rot} \op{rot} \vec A_e \approx - \sigma \mu \mu_0 \frac {\partial (\vec A_0 + \vec A_e)} {\partial t}. \label{rra2} \end{equation} $$

Далее будем считать, что ток в катушке L₁ изменяется по синусоидальному закону. С одной стороны, именно этот случай сейчас нам наиболее интересен (чаще всего ток в катушке металлоискателя является синусоидальным либо равен сумме нескольких синусоидальных составляющих). С другой стороны, предположение о синусоидальности тока в катушке существенно упрощает задачу. В этом случае все поля в системе и плотности вихревых токов также будут изменяться по синусоидальному закону с той же частотой (так как порождающим их первоисточником является ток в контуре и его переменное магнитное поле). Тогда имеет смысл воспользоваться методом комплексных амплитуд. При переходе к комплексным амплитудам, из уравнений исключается время (что понижает размерность задачи), а оператор дифференцирования по времени преобразуется просто в умножение на константу $ i \omega$, где $ i=\sqrt{-1}$, а $ \omega = 2 \pi f$ - циклическая частота колебаний тока в контуре (f - частота).

Так, уравнение \eqref{rra2}, записанное для комплексных амплитуд, примет вид $$ \begin{equation} \op{rot} \op{rot} \dot{\vec A_e} \approx - i \omega \sigma \mu \mu_0 (\dot{\vec A_0} + \dot{\vec A_e}). \label{rra1i} \end{equation} $$

Далее, для преобразования левой части уравнения, воспользуемся известным векторным тождеством $$ \op{rot} \op{rot} \vec a = \op{grad} \op{div} \vec a - \nabla^2 \vec a, $$ получим $$ \begin{equation} \op{grad} \op{div} \dot{\vec A_e} - \nabla^2 \dot{\vec A_e} \approx - i \omega \sigma \mu \mu_0 (\dot{\vec A_0} + \dot{\vec A_e}). \label{base_eq1} \end{equation} $$

Возвращаясь к соотношению \eqref{e_from_a}, $\vec E = - \partial \vec A / \partial t$, для комплексных амплитуд можем записать $$ \op{div} \dot{\vec E} = - \op{div} (i \omega \dot{\vec A}) = - i \omega \op{div} \dot{\vec A}. $$ С другой стороны, в соответствии уравнениями Максвелла (а именно, используем \eqref{mx4}, \eqref{mx6}), имеем $$ \op{div} \dot{\vec E} = \op{div} \frac {\dot{\vec D}} {\varepsilon \varepsilon_0} = \frac 1 {\varepsilon \varepsilon_0} \op{div} \dot{\vec D} = \frac {\dot \rho} {\varepsilon \varepsilon_0} = 0, $$ так как в нашем случае свободные заряды отсутствуют и всюду $\rho=0$.

Итак, получили $$ - i \omega \op{div} \dot{\vec A} = 0, $$ значит, если $\omega \ne 0$, то можем утверждать, что $$ \op{div} \dot{\vec A} = 0. $$ Полученное равенство есть не что иное, как калибровка Кулона, которая, как оказалось, применима в данном случае. Если величина равна нулю, то и градиент от неё также будет нулевым $$ \op{grad} \op{div} \dot{\vec A} = 0. $$ То же самое справедливо для составляющей векторного потенциала $\dot{\vec A_e}$. С учётом чего, уравнение \eqref{base_eq1} примет вид $$ \begin{equation} \nabla^2 \dot{\vec A_e} \approx i \omega \sigma \mu \mu_0 \dot{\vec A_0} + i \omega \sigma \mu \mu_0 \dot{\vec A_e}. \label{base_eq2} \end{equation} $$

Уравнение, связывая известное нам магнитное поле катушки и неизвестное магнитное поле вихревых токов, позволяет это неизвестное поле рассчитать. После чего можем определить полное магнитное поле как $\dot{\vec A} = \dot{\vec A_0} + \dot{\vec A_e}$.

Если мы вычислили A в некоторой точке проводящего тела, сможем определить плотность вихревых токов в этой точке, воспользовавшись уравнениями \eqref{mx7} и \eqref{e_from_a}, которые для комплексных амплитуд примут соответственно вид $$ \dot{\vec j} = \sigma \dot{\vec E}, \\ \dot{\vec E} = - i \omega \dot{\vec A}, $$ откуда получаем $$ \dot{\vec j} = - i \omega \sigma \dot{\vec A} = - i \omega \sigma (\dot{\vec A_0} + \dot{\vec A_e}). $$