Решение задачи в Octave. Влияние грунта

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

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

Параметры грунта
Пример расчёта влияния грунта
Предварительные выводы
Упрощённый расчёт влияния грунта
Больше примеров
Поиск в морской воде

Решение задачи в Octave. Приложение: справочник по функциям

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

Параметры грунта

Грунт способен оказать существенное влияние на работу чувствительного металлоискателя. Расчёт влияния сопряжён с определёнными трудностями, обусловленными тем, что параметры грунта могут изменяться в широком диапазоне значений, в зависимости от типа грунта в данной местности и условий. Что ещё хуже, грунт может быть крайне неоднородным даже в пределах одной площадки.

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

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

Чтобы выполнить какие-либо оценочные расчёты, необходимо задаться конкретными параметрами среды. Свойствам грунтов разных типов посвящено множество работ, так что не составляет особого труда найти массу интересной информации среди имеющихся в свободном доступе публикаций. Например, для первичных оценок, можно воспользоваться ГОСТ Р 50571.5.54-2013/МЭК 60364-5-54:2011^*. В таблице %tbl:gnd приведены некоторые данные из указанного ГОСТа и ряда других источников.

^* ГОСТ Р 50571.5.54-2013/МЭК 60364-5-54:2011. Электроустановки низковольтные. Часть 5-54. Заземляющие устройства, защитные проводники и защитные проводники уравнивания потенциалов.

В таблице даётся удельное сопротивление, удельная проводимость - обратная величина, \(\sigma=1/\lambda\).

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

Чем меньше удельная проводимость (чем больше удельное сопротивление), тем слабее влияние грунта. Исходя из этого, лучше всего искать в вечной мерзлоте или гранитных скалах (это, конечно, пока дело не дойдёт до копки). Неблагоприятные среды для поиска - болота, сильно увлажённые глинистые и минерализованные грунты, морская вода. Следует отметить, что проводимость одного и того же грунта сильно зависит от условий. Так, при высыхании грунта, его удельная проводимость снижается. Ещё сильнее она снижается при промерзании грунта (впрочем, при этом возникают гораздо более серьёзные проблемы - с извлечением найденных объектов). К сожалению, когда наступает наиболее удобное время для "полевых работ", т.е. весной и осенью (когда минимум растительности, а почва хорошо копается), грунт очень увлажнён, а потому имеет высокую удельную проводимость.

Пример расчёта влияния грунта

Воспользуемся методикой, описанной в предыдущей статье "Решение задачи в Octave. Примеры расчётов". При этом пойдём на ряд упрощений. Вместо "бесконечной" области (полупространства), заполненной грунтом, будем рассматривать влияние достаточно большого цилиндрического блока, находящегося непосредственно под катушкой металлоискателя (рис. %img:gnd). Грунт, заполняющий этот блок, будем считать однородным и немагнитным.

Рис. %img:gnd

Пример. Пусть заданы следующие параметры:
a₁ = 0.1 м (радиус катушки);
R = 1.0 м (радиус цилиндрического блока грунта);
H = 1.0 м (высота цилиндрического блока);
h = 0.01 м (расстояние от катушки металлоискателя до поверхности грунта);
f = 20 кГц (рабочая частота, w = 2 * pi * f);
sigma = 0.02 (Ом*м)^-1 (удельная проводимость грунта исходя из типичного значения удельного сопротивления 50 Ом*м для не слишком жёстких, но и не самых благоприятных условий поиска).

Катушку считаем расположенной достаточно близко к поверхности земли (в данном примере на расстоянии 1 см), что вполне соответствует действительности: во время поиска катушку подносят как можно ближе к земле, чтобы максимально использовать возможности металлоискателя по глубине поиска.

Выполняем вычисления и в зависимости от выбранной "густоты" сетки, с большей или меньшей точностью получаем результат (в расчёте на одновитковую катушку металлоискателя). Например, для
nrows = 50;
ncols = 50;
выполняя
dZ1_g=solve_md(a1, R, H, h, w, sigma, mu0, nrows, ncols),

получаем
dZ1_g = 1.2235e-007 - 7.8815e-012i.

Сразу замечаем, что изменение реактивного сопротивления (мнимая часть результата) ничтожно по сравнению с изменением активного (действительная часть), т.е. грунт практически не влияет на эквивалентную индуктивность катушки (по крайней мере, при умеренной удельной проводимости), но вносит заметное активное сопротивление.

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

Для сравнения, сопоставим влияние рассматриваемого "усреднённого" грунта с откликом на полезный объект. Допустим, мы ищем свою тестовую монету 5 копеек, которая находится на небольшой глубине 10 см. Рассчитаем отклик на монету без учёта влияния грунта для той же катушки на той же частоте:

a1=0.1;
R=12.5e-3;
H=1.5e-3;
h=0.1;
w=2*pi*20000;
sigma=sigma_brass;
nrows=20;
ncols=30;
dZ1_c=solve_md(a1, R, H, h, w, sigma, mu0, nrows, ncols),

получаем
dZ1_c = 6.5650e-007 - 2.5035e-006i.

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

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

Предварительные выводы

В сравнении с удельной проводимостью металлов, для поиска которых и предназначен металлоискатель, удельная проводимость даже самого высокопроводного грунта крайне низка. Вихревые токи оказываются крайне малы по величине, слабым оказывается создаваемое ими магнитное поле. Поэтому, те эффекты, которые характерны для металлов, такие например, как вытеснение поля на поверхность проводника, здесь практически не проявляются. С этим же связано то, что изменение индуктивности катушки под влиянием грунта, как правило, оказывается незначительным, меньшим чувствительности прибора и лишь в худших условиях достигает предела чувствительности или превышает его. Зато вносимое активное сопротивление оказывается достаточно большим. Если рассматривать полное сопротивление катушки переменному току, то его изменение под воздействием грунта вполне может значительно превзойти изменения, получаемые в ответ даже на достаточно "мощные" полезные цели. Потому важно контролировать не просто изменение сопротивления, но и по отдельности отслеживать его реактивную и активную составляющие, т.е. с высокой точностью следить за фазовыми соотношениями сигнала-отклика. При этом погрешности в определении фазы в условиях значительного влияния грунта, вполне способны привести к ложному определению целей. Поэтому, грунт с высокой удельной проводимостью способен снизить чувствительность металлоискателя (ну или мы сами вынуждены снижать чувствительность для устранения ложных срабатываний).

Упрощённый расчёт влияния грунта

Удельная проводимость грунта мала, а возникающие в нём вихревые токи слабы. Поэтому, в первом приближении, при расчёте величины плотности вихревых токов, можем пренебречь магнитным полем этих вихревых токов и считать, что магнитное поле в грунте определяется только внешним магнитным полем, $$ \dot{\vec A} \approx \dot{\vec A}_0. $$ Внешнее поле A₀ известно нам в любой точке пространства. Зная векторный потенциал, вычисляем плотность вихревых токов $$ \dot{\vec j} = - i \omega \sigma \dot{\vec A} \approx - i \omega \sigma \dot{\vec A}_0, $$ где
\(\omega\) - циклическая частота внешнего поля;
\(\sigma\) - удельная проводимость среды в рассматриваемой точке.
Теперь, зная распределение плотности вихревых токов в грунте, сможем вычислить создаваемое ими магнитное поле и наводимую этим (переменным) полем ЭДС индукции в катушке. А значит, и изменение эквивалентного активного сопротивления катушки под влиянием грунта.

В случае, когда источником внешнего магнитного поля является круговой контур с током, как нам известно, линии вектора A₀ имеют вид концентрических окружностей. Такого же вида будут линии вектора плотности вихревых токов. Тогда, для вычисления магнитного поля, создаваемого вихревыми токами, мы можем разделить их на множество круговых токов, вычислить поле каждого такого тока в интересующей нас точке и просуммировать полученные значения.

Пусть функция F(a, r, z) вычисляет в точке с координатами (r, z) азимутальную проекцию векторного потенциала, создаваемого круговым током с радиусом a в расчёте на силу тока 1 А в круговом контуре. Координаты точки задаются в цилиндрической системе, привязанной к контуру с током. Если магнитное поле создаётся круговым контуром радиуса a₁ с током I₁, то проекция векторного потенциала в некоторой точке $$ A_0 = I_1 F(a_1, r, z), $$ а проекция плотности вихревого тока в этой точке $$ j = - i \omega \sigma I_1 F(a_1, r, z). $$ Для вычисления ЭДС индукции, наводимой в контуре-индукторе магнитным полем вихревых токов, нам нужно определить проекцию векторного потенциала поля, создаваемого вихревыми токами в любой точке контура индуктора (в данном случае во всех точках контура эта величина одна и та же). Проекция для поля, создаваемого кольцевым вихревым током с малым сечением dr dz, который проходит через рассматриваемую точку (r, z), с учётом перехода к новой локальной системе координат, привязанной к круговому току, может быть рассчитана как $$ dA_1 = F(r, a_1, z) j dr dz, $$ где учтено, что $$ F(r, a_1, -z) = F(r, a_1, z). $$ Подставляя значение для j, получаем $$ dA_1 = - i \omega \sigma I_1 F(r, a_1, z) F(a_1, r, z) dr dz, $$ тогда для поля, создаваемого вихревыми токами в цилиндрической области (соосной с контуром-индуктором), заполненной грунтом, радиус которой R, высота H, а расстояние от плоскости контура-индуктора до поверхности грунта h, рассчитывается как $$ A_1 = - i \omega \sigma I_1 \int_h^{H+h}\int_0^R F(r, a_1, z) F(a_1, r, z) dr dz. $$ Теперь без труда вычисляем наводимую магнитным полем вихревых токов ЭДС в контуре индуктора и вычисляем вносимое активное сопротивление: $$ \Delta E_1 = - 2 \pi i \omega a_1 A_1, \\ \Delta Z_1 = - \Delta E_1 / I_1, \\ \Delta Z_1 = 2 \pi \omega^2 \sigma a_1 \int_h^{H+h}\int_0^R F(r, a_1, z) F(a_1, r, z) dr dz. $$

Как видим, результат - изменение эквивалентного сопротивления контура - содержит только вещественную составляющую, т.е. изменяется только активная составляющая сопротивления. Это подтверждает сделанный ранее вывод о том, что если удельная проводимость среды настолько мала, что магнитным полем возникающих в ней вихревых токов можно пренебречь по сравнению с внешним магнитным полем, то воздействие такой среды на катушку будет проявляться только в изменении активного сопротивления катушки. Реактивное сопротивление остаётся неизменным (в реальности - практически неизменным).

Подынтегральное выражение F(r, a, z) F(a, r, z) можно несколько упростить, если учесть, что для кругового тока данная функция имеет вид $$ F(a, r, z) = \frac {\mu \mu_0} {2 \pi b} \sqrt {\frac {a} r} \left( (2-b^2) K(b)-2E(b)\right), \\ b^2 = \frac {4 a r} {z^2 + r^2 + a^2 + 2 a r}, $$ где K(b), E(b) - полные эллиптические интегралы первого и второго рода соответственно.

Нетрудно заметить, что параметр b не изменяется при взаимной перестановке переменных a и r. Тогда $$ F(r, a, z) F(a, r, z) = \frac {\mu^2 \mu_0^2} {4 \pi^2} \frac 1 {b^2} \left( (2-b^2) K(b)-2E(b)\right)^2. $$

В Octave описанный алгоритм можно записать приблизительно следующим образом.

Здесь __A0q(a, r, hz) - вспомогательная функция, используемая для задания подынтегрального выражения.

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

Больше примеров

Ранее мы предположили, что влияние грунта на катушку можно оценить по влиянию достаточно большого цилиндрического блока из этого грунта. Хотелось бы знать, насколько большим должен быть блок, чтобы точность результата оказалась приемлемой. Чтобы немного прояснить этот вопрос, рассмотрим следующую задачу.

Имеется катушка радиуса a₁, на которой проверяется воздействие цилиндрических блоков разного размера из грунта с проводимостью sigma. Требуется определить величину этого воздействия, если размеры блоков равны R = H = k * a₁, где k - заданный коэффициент (или несколько заданных значений коэффициента).

Решаем, используя функцию calc_dZ1_gnd(a1, H, R, h, w, sigma, mu_abs, tol).

a1=0.1;
h=0.01;
w=2*pi*20000;
sigma=0.02;
tol=[]; % точность - по умолчанию

k=[0.5, 1, 2, 4, 8, 12, 25, 50, 100];
kmax=200;

dZ1max=calc_dZ1_gnd(a1, kmax*a1, kmax*a1, h, w, sigma, mu0, tol);

for p=1:length(k)
    dZ1(p)=calc_dZ1_gnd(a1, k(p)*a1, k(p)*a1, h, w, sigma, mu0, tol);
endfor

plot(k, dZ1/dZ1max);

По результатам вычислений стоим график (рис. %img:gnd_v) зависимости вносимого активного сопротивления от относительного размера блока (относительно радиуса катушки). Вносимое сопротивление также измеряется в относительных единицах (относительно значения, соответствующего выбранному достаточно "большому" блоку). Как видим, в нашем примере, влияние всего грунта можно оценить по влиянию цилиндрического блока из грунта, линейные размеры которого в 5..10 раз превышают радиус катушки. Можно предположить, что и в других ситуациях это останется справедливым, по крайней мере, когда катушка расположена достаточно близко к поверхности земли. При желании можно продолжить изучение этого вопроса - теоретически или путём вычислительных экспериментов.

Рис. %img:gnd_v

Также интересно было бы определить, что важнее для получения большей точности - взять блок "потолще" или "пошире". Для ответа на вопрос, сначала возьмём цилиндр с заведомо большим радиусом и рассмотрим зависимость влияния на катушку от высоты цилиндра (толщины слоя грунта). Затем, возьмём цилиндр с заведомо большой высотой и посмотрим, какова зависимость влияния на катушку от радиуса цилиндра.

a1=0.1;
h=0.01;
w=2*pi*20000;
sigma=0.02;
tol=[]; % точность - по умолчанию

k=[0.5, 1, 2, 4, 8, 12, 25, 50];
kmax=200;

dZ1max=calc_dZ1_gnd(a1, kmax*a1, kmax*a1, h, w, sigma, mu0, tol);

for p=1:length(k)
    dZ1(p)=calc_dZ1_gnd(a1, k(p)*a1, kmax*a1, h, w, sigma, mu0, tol);
endfor

plot(k, dZ1/dZ1max);

В случае блока из грунта большого радиуса, значения отклика, близкие к предельным, достигаются достаточно быстро, уже при толщине слоя грунта (высоте цилиндра) в 5..10 радиусов катушки (рис. %img:gnd_h).

Рис. %img:gnd_h

a1=0.1;
h=0.01;
w=2*pi*20000;
sigma=0.02;
tol=[]; % точность - по умолчанию

k=[0.5, 1, 2, 4, 8, 12, 25, 50];
kmax=200;

dZ1max=calc_dZ1_gnd(a1, kmax*a1, kmax*a1, h, w, sigma, mu0, tol);

for p=1:length(k)
    dZ1(p)=calc_dZ1_gnd(a1, kmax*a1, k(p)*a1, h, w, sigma, mu0, tol);
endfor

plot(k, dZ1/dZ1max);

Если взят блок большой толщины (высота цилиндра велика), то значения отклика, близкие к предельным, достигаются при радиусе блока в 10..20 радиусов катушки (рис. %img:gnd_r).

Рис. %img:gnd_r

Итак, анализ отклика катушки на грунт можно заменить анализом отклика на достаточно большой цилиндр из грунта. Высота цилиндра должна быть 5..10 раз больше радиуса катушки, радиус цилиндра можно выбрать немного побольше высоты. Так, в нашем примере, если взять H=5*a₁, R=15*a₁, получим точность оценки влияния грунта порядка 10%, что можно считать весьма неплохим результатом.

Теперь рассмотрим зависимость влияния грунта на катушку от частоты.

a1=0.1;
H=5*a1;
R=15*a1;
h=0.01;
sigma=0.02;
tol=[]; % точность - по умолчанию

f=[5, 10, 15, 20, 25, 30]*1000;
w=2*pi*f;

for p=1:length(w)
    dZ1(p)=calc_dZ1_gnd(a1, H, R, h, w(p), sigma, mu0, tol);
endfor

plot(f, dZ1);

Рис. %img:gnd_f

Как и следовало ожидать, влияние усиливается с ростом частоты, зависимость является квадратичной (это легко проверить, вычислив отношение \(dZ1/f^2\), на разных частотах оно оказывается одним и тем же).

Поиск в морской воде

Одной из наиболее неблагоприятных сред, если оценивать по удельной проводимости, где иногда приходится осуществлять поиск, является морская вода. Рассмотрим, насколько сильное она оказывает влияние. Будем считать, что катушка погружена в воду (рис. %img:wtr).

Рис. %img:wtr

a1=0.1;
H=10*a1;
R=15*a1;
h=-H/2;  % катушка находится в толще воды
w=2*pi*20000;
sigma=5; % исходя из оценки удельного сопротивления среды 0.2 Ом*м
tol=[];  % точность - по умолчанию

dZ1=calc_dZ1_gnd(a1, H, R, h, w, sigma, mu0, tol)

Тот же результат можно получить с учётом того, что воздействие среды сверху и снизу вдвое больше воздействия только с одной стороны:

a1=0.1;
H=10*a1;
R=15*a1;
w=2*pi*20000;
sigma=5; % исходя из оценки удельного сопротивления среды 0.2 Ом*м
tol=[];  % точность - по умолчанию

dZ1=2*calc_dZ1_gnd(a1, H/2, R, 0, w, sigma, mu0, tol)

И в том, и в другом случае получаем:
dZ1 = 7.7216e-005.

Это очень большое значение. Вначале документа мы вычисляли влияние монеты 5 копеек, находящейся на расстоянии 10 см от катушки и получили значение dZ1_c = 6.5650e-007 - 2.5035e-006i. Сравнивая, видим, что изменение активного сопротивления под влиянием среды (морской воды) здесь на 2 порядка больше, чем изменение под влиянием полезного объекта. Более того, изменение активного сопротивления катушки за счёт воздействия среды в данном случае оказалось в 30 раз больше, чем изменение реактивного сопротивления под влиянием монеты. Что ещё раз подчёркивает важность не просто измерения изменения полного сопротивления, но и учёта его активной и реактивной составляющих.

Впрочем, в столь тяжёлых условиях поиска, влияние среды на индуктивное сопротивление катушки также становится ощутимым. Рассчитаем изменение комплексного сопротивления катушки с использованием функции solve_md(a1, a2, d2, h, w, sigma, mu_abs, nrows, ncols)

a1=0.1;
H=10*a1;
R=15*a1;
w=2*pi*20000;
sigma=5; % исходя из оценки удельного сопротивления среды 0.2 Ом*м

nrows=50;
ncols=150;

dZ1=2*solve_md(a1, R, H/2, 0, w, sigma, mu0, nrows, ncols)

Результат:
dZ1 = 7.7232e-05 - 1.0899e-06i
Активная составляющая с хорошей точностью согласуется с результатом упрощённого вычисления с помощью функции calc_dZ1_gnd. Кроме того, видим появление заметной реактивной составляющей в изменении сопротивления катушки. По сравнению с величиной, полученной выше при рассмотрении примера с монетой (dZ1_c = 6.5650e-007 - 2.5035e-006i), изменение реактивной составляющей сопротивления под влиянием среды всего чуть более чем в 2 раза меньше, чем отклик на монету. Что очень влияет на поиск - при увеличении расстояния до объекта или с уменьшением размеров объекта, полезный отклик становится много меньше отклика на среду.