Решение задачи в Octave. Примеры расчётов
В предыдущей статье была подробно рассмотрена идея численного решения задачи о расчёте вихревых токов с теоретической точки зрения. Теперь перейдём к практической реализации этой идеи. В качестве инструмента, как можно догадаться из заголовка, воспользуемся системой Octave.
metal_detector.zip
m-файл Octave для решения задачи расчёта вихревых токов.
Предлагаемое решение оформлено в виде m-файла для системы Octave. Файл содержит необходимые функции для решения задачи расчёта вихревых токов для случая, когда рассматриваются вихревые токи в проводящем цилиндре, который находится во внешнем переменном магнитном поле, обладающем осевой симметрией относительно оси цилиндра. Файл также содержит некоторые полезные для работы константы (магнитная постоянная, удельные проводимости некоторых материалов). Файл достаточно подробно комментируется и включает всю необходимую информацию о назначении функций. Наиболее важными функциями файла являются A0field, solve_pde, Aeddy, solve_md.
Рассмотрим применение функций на нескольких примерах. Во всех предлагаемых задачах будем считать, что источником внешнего магнитного поля является круговой контур с переменным током амплитудой 1 A и заданной частотой (если ток иной или количество витков в катушке более одного, можно просто выполнить масштабирование результата). Случай, когда источником внешнего магнитного поля является круговой контур, наиболее интересен с точки зрения анализа металлоискателей. Для расчёта такого поля в документе имеется функция A0field. Источник внешнего магнитного поля может быть и иным, тогда просто нужно будет задавать внешнее поле самостоятельно. Важно, чтобы оно обладало осевой симметрией.
В заданном внешнем переменном магнитном поле находится цилиндр из проводящего ток материала. Известны радиус, высота цилиндра, удельная проводимостью материала.
Задача состоит в определении распределения вихревых токов в цилиндре и в расчёте изменения индуктивности и активного сопротивления контура под влиянием этого проводящего цилиндра.
Рис. %img:grid
Как мы уже выяснили, распределение плотности вихревых токов, распределение магнитных полей (полного и его составляющих - внешнего поля и поля, создаваемого вихревыми токами), одинаковое во всех сечениях цилиндра, проходящих через его ось. Поэтому рассматриваемая задача является двумерной. Все величины задаются и вычисляются для одного из сечений цилиндра (половины сечения, если быть точным); вычисления производятся для точек - узлов сетки, которая вводится для разделения сечения на маленькие прямоугольные ячейки (рис. %img:grid).
Все рассматриваемые векторы (плотность вихревого тока, векторные потенциалы) здесь имеют только одну проекцию, которая может быть отлична от нуля. Это азимутальная проекция. Всюду будем оперировать с этими проекциями, т.е. со скалярными величинами. Кроме того, величины считаем переменными, изменяющимися по синусоидальному закону с заданной частотой, что даёт возможность применить метод комплексных амплитуд. Итак, при решении задачи используются комплексные скалярные величины, по которым без труда могут быть найдены амплитуды и фазы величин (относительно фазы внешнего магнитного поля).
Вначале рассмотрим характер распределения вихревых токов в объёме массивных тел. Допустим, у нас есть медный диск диаметром 30 см и толщиной 1 см, который мы помещаем в магнитное поле кругового контура с током. Требуется установить распределение плотности вихревых токов по сечению цилиндра. Для решения задачи, запустим предложенный выше m-файл, после чего выполним следующие команды.
% Радиус контура с током.
a1=0.1;
% Расположение цилиндра относительно контура (z-координата верхнего
% торца цилиндра; здесь цилиндр расположен ниже плоскости контура на
% расстоянии |z0|).
z0=-0.01;
% Высота и радиус цилиндра.
H=0.01;
R=0.15;
% Удельная проводимость материала цилиндра (медь).
sigma=sigma_cu; % значение sigma_cu определено в основном m-документе.
% Задаём количество узлов в сетке.
nrows=40;
nclos=50;
% Задаём частоту переменного тока в контуре; начинаем с заведомо низкой
% частоты, чтобы наблюдать проникновение поля на значительную глубину.
f=50; w=2*pi*f;
% Рассчитываем магнитное поле индуктора в узлах сетки (поле кругового
% контура радиуса a1 с переменным током амплитудой 1А и циклической
% частотой w).
A0=A0field(a1, z0, H, R, nrows, ncols, mu0);
% Рассчитываем магнитное поле вихревых токов в узлах сетки.
Ae=solve_pde(A0, w, sigma, H, R, mu0);
% Полное магнитное поле в узлах сетки.
A=A0+Ae;
% Плотность вихревых токов в узлах сетки.
j=-i*w*sigma*A;
Теперь отобразим полученные данные в виде карты распределения плотности тока (амплитудные значения) по сечению цилиндра командой pcolor(abs(j)); или, для большей наглядности получаемого результата, выполним дополнительные действия: допишем слева к матрице j симметричную ей матрицу, чтобы получить матрицу, описывающую полное сечение цилиндра; укажем диапазоны изменения величин r, z, которым соответствуют значения в матрице j; построим дополнительные графические элементы - изображение контура-источника внешнего магнитного поля для цилиндра, нулевой уровень.
% Дополняем матрицу, чтобы она описывала полное сечение цилиндра.
j=cat(2, fliplr(j), j);
% Указываем область построения - формируем векторы с координатами
% точек, соответствующими элементам дополненной матрицы j.
rrange=linspace(-R, R, 2*ncols);
hrange=linspace(z0, z0-H, nrows);
[X Y]=meshgrid(rrange, hrange);
% Выводим карту распределения абсолютных значений j (плотность вихревых токов
% пропорциональна A и распределение A качественно будет таким же).
pcolor(X, Y, abs(j));
% Задаём цвет фона - серый.
set(gca(), "color", [0.65 0.65 0.65]);
% Сглаживаем изображение.
shading("interp");
% Для наглядности "пририсовываем" схематическое изображение контура-индуктора.
ld=0.005; % диаметр поперечного сечения контура на изображении
rectangle("Position", [-a1-ld/2, -ld/2, ld, ld], "Curvature", [1, 1], "FaceColor", "red");
rectangle("Position", [+a1-ld/2, -ld/2, ld, ld], "Curvature", [1, 1], "FaceColor", "red");
line([-a1, a1], [-ld/2, -ld/2], "linestyle", "--", "color", "black");
line([-a1, a1], [+ld/2, +ld/2], "linestyle", "--", "color", "black");
% Прорисовываем 0-уровень.
ax=axis;
line(ax(1:2), [0, 0], "linestyle", "--", "color", "blue");
Подобрав размеры окна с графиком так, чтобы масштабы по горизонтали и по вертикали примерно совпадали, получим изображение как на рис. %img:cu_50a.
Рис. %img:cu_50a (50 Гц)
Из-за многократной разницы между диаметром диска и его толщиной, характер распределения токов на этом изображении виден не очень хорошо. Поэтому на рис. %img:cu_50 предложен другой вариант, с сильным растяжением по вертикали.
Рис. %img:cu_50 (50 Гц)
Теперь повторим те же самые действия, но увеличим частоту с 50 Гц до 500 Гц. Распределение вихревых токов окажется уже несколько иным (рис. %img:cu_500). Как и следовало ожидать, с ростом частоты, глубина проникновения электромагнитного поля в проводник уменьшается; вихревые токи сосредотачиваются вблизи поверхности проводника.
Рис. %img:cu_500 (500 Гц)
Теперь увеличим частоту до 5000 Гц, т.е. до значений из диапазона рабочих частот типичного металлоискателя. Нетрудно заметить (рис. %img:cu_5000), что на таких частотах глубина проникновения переменного магнитного поля в проводник становится совсем малой; вихревые токи сосредоточены в тонком слое вблизи поверхности проводника.
Рис. %img:cu_5000 (5000 Гц)
Мы рассмотрели характер распределения вихревых токов в массивных металлических телах. Теперь перейдём к изучению вихревых токов в телах, малых по сравнению с катушкой. Таковыми являются, например, монеты. В качестве объекта для исследования выберем диск диаметром 25 мм и толщиной 1.5 мм из латуни (что примерно соответствует монете 5 копеек СССР). Выполняемые действия совершенно аналогичны тем, что были рассмотрены в предыдущем пункте, поэтому не будем повторяться и сразу перейдём к результатам (рис. %img:bc_5000e).
Рис. %img:bc_5000e (5000 Гц)
Изображение монеты получилось слишком маленьким. На следующем рисунке (рис. %img:bc_5000) монета даётся крупным планом, без дополнительных графических элементов.
Рис. %img:bc_5000 (5000 Гц)
Для поиска таких мелких объектов лучше подходят более высокие частоты; для данной монеты оптимальны частоты около 20 кГц или выше. На частоте 20 кГц получаем следующую картину распределения вихревых токов: рис. %img:bc_20k.
Рис. %img:bc_20k (20 кГц)
Как мы и предполагали, вихревые токи сосредоточены, преимущественно, вблизи ободка диска.
Перейдём к самому важному вопросу - вычислению отклика катушки на объект, т.е. вычислению относительного изменения индуктивности катушки под воздействием находящегося рядом проводящего предмета. Решать задачу будем следующим образом.
Сначала задаём параметры катушки, объекта; указываем прочие условия.
% Задаём параметры катушки и объекта:
a1=0.1; % радиус катушки
R=12.5e-3; % радиус цилиндра-объекта (радиус монеты, например)
H=1.5e-3; % высота цилиндра (толщина монеты, например)
h=0.05; % глубина (расстояние от плоскости катушки до
% ближайшей поверхности монеты)
w=2*pi*20000; % циклическая частота
sigma=sigma_brass; % удельная проводимость (материала монеты)
% Задаём параметры сетки (больше узлов, мельче ячейка - выше точность
% и больше время счёта).
nrows=20; ncols=25;
Затем приступаем собственно к решению задачи.
% Рассчитываем внешнее поле в узлах сетки.
A0=A0field(a1, -h, H, R, nrows, ncols, mu0);
% Решаем задачу по вычислению магнитного поля вихревых токов.
Ae=solve_pde(A0, w, sigma, H, R, mu0);
% Полное магнитное поле в узлах сетки.
A=A0+Ae;
% Вычисляем магнитное поле, создаваемое вихревыми токами проводящего
% предмета в одной из точек на контуре индуктора (для вычисления
% магнитного потока через поверхность индуктора, создаваемого
% вихревыми токами).
Ap1=Aeddy(A, w, sigma, H, R, h, a1, mu0);
% ЭДС индукции, наводимая в контуре-индукторе под влиянием магнитного
% поля вихревых токов.
dE1=-2*i*w*pi*a1*Ap1;
% Изменение комплексного сопротивления контура (катушки), с учётом
% того, что амплитуда тока в контуре считается равной 1 А.
dZ1=-dE1/1.0
% Абсолютное изменение индуктивности контура (разделив на начальную
% индуктивность контура, получим относительное изменение).
dL1=imag(dZ1)/w
% Абсолютное изменение активного сопротивления контура
dR1=real(dZ1)
В этом примере получаем результат:
dZ1 = 2.7366e-006 - 1.0370e-005i,
dL1 = -8.2521e-011,
dR1 = 2.7366e-006.
Следует иметь в виду, что результаты (dZ1, dL1, dR1) получены для контура, т.е. в расчёте на один виток. В случае многовитковой катушки, dZ1 должно быть умножено на квадрат количества витков в катушке, соответственно, во столько же раз изменятся значения dL1 и dR1.
Ввиду важности вычисления влияния объектов на катушку, перечисленные здесь действия объединены в одной функции solve_md(a1, a2, d2, h, w, sigma, mu_abs, nrows, ncols). Поэтому рассмотренная выше задача может быть решена за один вызов этой функции:
dZ1=solve_md(a1, R, H, h, w, sigma, mu0, nrows, ncols)
dL1=imag(dZ1)/w
dR1=real(dZ1)
Для приблизительной оценки относительного изменения индуктивности катушки можно воспользоваться функцией calc_el1(a1, a2, d2, h, w, sigma, mu_abs, nrows, ncols). Что очень удобно, относительное изменение уже не зависит от количества витков в катушке. Для нашего примера получаем результат -1.5480e-004.
Воспользовавшись функцией calc_el1(a1, a2, d2, h, w, sigma, mu_abs, nrows, ncols), мы можем решить множество важных задач. Например, исследовать зависимость отклика от размеров катушки, размеров объекта, проводимости его материала, используемой рабочей частоты.
Рассмотрим зависимость отклика катушки от рабочей частоты. В качестве объекта для исследования возьмём всё ту же монету 5 копеек.
a1=0.1; R=12.5e-3; H=1.5e-3; h=0.15;
sigma=sigma_brass;
nrows=20; ncols=25;
% Частоты, для которых вычисляем отклик.
f=[1, 2, 4, 8, 12, 16, 20, 25, 40, 60]*1000;
w=2*pi*f;
% Вычисляем отклик на заданных частотах.
for k=1:length(w) resp(k)=abs(calc_el1(a1, R, H, h, w(k), sigma, mu0, nrows, ncols)); endfor
% Приблизительное значение предельного отклика вычисляем как отклик
% на заведомо высокой частоте (возможно, потребуется увеличить
% количество элементов сетки, чтобы размер ячейки не оказался больше
% глубины проникновения поля в объект).
respmax=abs(calc_el1(a1, R, H, h, 2*pi*4e6, sigma, mu0, 30, 50));
% Строим график - зависимость отклика в относительных величинах
% (относительно предельного значения) от частоты.
plot(f, resp/respmax);
Рис. %img:fr
Получаем зависимость, изображённую на рис. %img:fr. График построен для относительных значений отклика (значений отклика относительно предельного значения). Что удобно, относительные значения практически не зависят от радиуса катушки, расстояния между катушкой и объектом, поэтому значения этих параметров могут быть заданы произвольно. Можете проверить это, например, с помощью вычислительного эксперимента.
При низких частотах имеем довольно резкую зависимость отклика от частоты, работа на этом участке неэффективна из-за недоиспользования чувствительности металлоискателя. С увеличением частоты, рост отклика замедляется. На частоте около 20 кГц отклик на рассматриваемую монету достигает уровня 0.8 от предельного значения, т.е. с дальнейшим увеличением частоты он уже не увеличивается существенным образом. Таким образом, для объектов данного типа, рабочую частоту в районе 20 кГц можно считать оптимальной. Допустимы отклонения от этой частоты, как в сторону увеличения, так и в сторону снижения: увеличение частоты только увеличивает отклик, а умеренное снижение не приводит к слишком сильному падению отклика - в данной области зависимость отклика от частоты сравнительно слабая.
Смотрите далее: "Решение задачи в Octave. Влияние грунта".
author: hamper; date: 2020-03-23
