Журнал вычислительной математики и математической физики, 2019, T. 59, № 8, стр. 1410-1418

1. ВВЕДЕНИЕ

Для решения многих практических задач из области магнетизма, например, задач неразрушающего магнитного контроля, актуальной является проблема создания обоснованных алгоритмов аналитического или численного решения задач магнитостатики применительно к магнитным телам различной формы с возможными инородными включениями. Обоснованные эффективные методы решения такого рода задач имеются по большей части для безграничных модельных осесимметричных тел относительно простой геометрической формы с бесконечно протяженными дефектными включениями в виде цилиндра с эллипсоидальным сечением, помещенных в постоянное внешнее магнитное поле определенного (удобного для исследования) направления, что позволяло во многих случаях свести задачу к двумерной (см., например, [1]–[4]). Такие ограничения в основном связаны с достаточно серьезными математическими проблемами, возникающими при исследовании и выработке оптимального аналитического (или численно-аналитического) алгоритма решения получающихся из системы уравнений Максвелла дифференциальных или интегральных уравнений на многообразиях сложной формы. Однако большой теоретический и практический интерес представляют задачи для тел с включениями конечных размеров в произвольном внешнем поле. Например, в [5] в этом направлении исследуется шаровой дефект в магнетике шаровой формы. В настоящей работе получены формулы для расчета напряженности результирующего магнитного поля для однородной магнитной пластины конечной толщины с произвольным дефектом (типа полости или инородного вкрапления) конечных размеров, помещенной в произвольное внешнее магнитное поле. Подобного рода задачи решались, например, в работах [6]–[8]. В работе [6] получены аналитические формулы для оценивания магнитостатических полей при намагничивании тонкой однородной ферромагнитной пластины, но при отсутствии дефектов в ней. В работе [7] получены формулы для расчета напряженности результирующего магнитного поля для однородного магнитного полупространства с конечным дефектом в виде полости, помещенного в произвольное внешнее магнитное поле. В работе [8] тоже исследуется поле от пластины с конечным дефектом, однако рассмотрен дефект только в виде воздушной полости и при специальной конфигурации внешнего поля.

2. ОБЩАЯ СХЕМА РАСЧЕТА НАПРЯЖЕННОСТИ РЕЗУЛЬТИРУЮЩЕГО МАГНИТНОГО ПОЛЯ

Для решения указанных выше задач мы исходим из так называемого основного уравнения магнитостатики (см., например, [9, с. 16])

Это уравнение эквивалентно системе уравнений Максвелла для случая магнитостатики и связывает искомую напряженность результирующего магнитного поля ${\mathbf{H}}({\mathbf{r}}) = \{ {{H}_{x}}({\mathbf{r}}),{{H}_{y}}({\mathbf{r}}),{{H}_{z}}({\mathbf{r}})\} $ в произвольной точке пространства ${\mathbf{r}} = (x,y,z)$ (не лежащей на границе магнетика) с заданным полем ${{{\mathbf{H}}}^{0}}({\mathbf{r}}) = \{ H_{x}^{0}({\mathbf{r}}),H_{y}^{0}({\mathbf{r}}),H_{z}^{0}({\mathbf{r}})\} $ внешнего источника. В данном уравнении $\Omega $ есть область в пространстве ${{R}^{3}}$, ограниченная поверхностью $S$ и занятая исследуемым магнетиком с заданной магнитной проницаемостью $\mu $, постоянной (однородный магнетик) или зависящей от координат (неоднородный магнетик).

Напряженность результирующего поля ${\mathbf{H}}({\mathbf{r}})$ внутри и вне магнетика определяется из (2.1) по следующей схеме. Полагая ${\mathbf{r}} \in \Omega $ в (2.1), получают интегро-дифференциальное уравнение относительно напряженности поля ${\mathbf{H}}({\mathbf{r}})$ внутри магнетика. Доказано (см., например, [10]–[13]), что это уравнение, рассматриваемое в вещественном гильбертовом пространстве ${{{\mathbf{L}}}_{2}}(\Omega )$ вектор-функций ${\mathbf{M}}({\mathbf{r}})$ с (конечной) нормой

и скалярным произведением имеет единственное решение как для ограниченных, так и неограниченных областей $\Omega $. После решения указанного уравнения напряженность поля ${\mathbf{H}}({\mathbf{r}})$ вне магнетика вычисляется непосредственно из соотношения (2.1), в котором точка наблюдения ${\mathbf{r}}$ берется вне магнетика:

Понятно, что основная сложность при реализации описанной схемы заключается в ее первой части, т.е. получении решения уравнения (2.1) для ${\mathbf{r}} \in \Omega $.

3. УРАВНЕНИЕ ДЛЯ ПОЛЯ МАГНИТНОЙ ПЛАСТИНЫ С ВНУТРЕННИМ ДЕФЕКТОМ

Рассмотрим ситуацию, когда магнетик, занимающий область $\Omega $ (конечную или нет) является неоднородным в том смысле, что внутри $\Omega $ есть конечная подобласть ${{\Omega }_{d}} \subset \Omega $, в которой постоянная магнитная проницаемость ${{\mu }_{d}}$ отлична от постоянной магнитной проницаемости $\mu $ остальной части магнетика ${{\Omega }_{1}}: = \Omega {\backslash }{{\Omega }_{d}}$ (см. фиг. 1). Такую область ${{\Omega }_{d}}$ можно интерпретировать как дефект (инородное включение или воздушная полость), обнаружение и описание которого являются основной задачей неразрушающего контроля изделий методами магнитной дефектоскопии. Для такой конфигурации магнетиков соотношение (2.1) может быть приведено к виду (см., например, [7], [13, с. 257]):

где ${{S}_{1}}$ – граница области $\Omega $, $S$ – граница области ${{\Omega }_{d}}$, ${{H}_{n}}$– предельные значения изнутри области ${{\Omega }_{1}}$ нормальных составляющих на ${{S}_{1}}$ и $S$ вектора напряженности ${\mathbf{H}}$, а векторы единичной нормали ${\mathbf{n}}$ на ${{S}_{1}}$ и $S$ для удобства выбраны внешними по отношению к областям $\Omega $ и ${{\Omega }_{d}}$ соответственно. Как видно из (3.1), для вычисления напряженности ${\mathbf{H}}({\mathbf{r}})$ в любой точке вне границы раздела сред ${{S}_{1}}$ и $S$ необходимо найти значения нормальных составляющих искомой напряженности на этих границах. Для такого поиска можно поочередно скалярно умножить обе части (3.1) на вектор единичной нормали ${\mathbf{n}}$ на ${{S}_{1}}$ и $S$ соответственно, совершить последовательно предельный переход на ${{S}_{1}}$ и $S$ изнутри области ${{\Omega }_{1}}$ с учетом формулы предельного значения на поверхности нормальной производной потенциала простого слоя. В результате получается система интегродифференциальных уравнений относительно искомых нормальных составляющих ${{H}_{n}}$ на ${{S}_{1}}$ и $S$. Однако эта система оказывается в общем случае слишком сложной для построения эффективного алгоритма ее решения. Далее рассматривается один из возможных случаев, когда упомянутую систему можно заменить одним уравнением.

Пусть область $\Omega $, занятая магнетиком с постоянной магнитной проницаемостью μ, заключена между двумя параллельными плоскостями ${{S}_{1}}$ и ${{S}_{2}}$. Внутри имеется произвольная “дефектная” область ${{\Omega }_{d}} \subset \Omega $, ограниченная поверхностью $S$ и занятая магнетиком с постоянной проницаемостью ${{\mu }_{d}}$ (случай ${{\mu }_{d}} = 1$ соответствует воздушной полости, а ${{\mu }_{d}}$>$1$ соответствует инородному вкраплению). Введем пространственную декартову систему координат, начало которой выберем внутри ${{\Omega }_{d}}$, а ось $z$ направим перпендикулярно параллельным плоскостям ${{S}_{1}}$ и ${{S}_{2}}$ (см. фиг. 2). Пусть параметры ${{d}_{1}}$ и ${{d}_{2}}$ имеют смысл расстояния от начала координат до плоскостей ${{S}_{1}}$ и ${{S}_{2}}$ соответственно, а потому соотношения $z = - {{d}_{1}}$ и $z = {{d}_{2}}$ суть уравнения этих плоскостей. Для такой конфигурации магнетиков уравнение (3.1) имеет вид

где ${{H}_{z}}({\mathbf{r}})$ в поверхностных интегралах суть предельные (изнутри области ${{\Omega }_{1}}$) значения на ${{S}_{1}}$ и ${{S}_{2}}$ z-компоненты напряженности ${\mathbf{H}}({\mathbf{r}})$. Переходя в первых двух интегралах в (3.2) от поверхностных к двойным, получаем где для краткости и удобства введены обозначения

Наша дальнейшая цель – избавиться в (3.3) от функций ${{\psi }_{1}}(x,y)$, ${{\psi }_{2}}(x,y)$ (т.е. от ${{H}_{z}}(x,y, - {{d}_{1}})$, ${{H}_{z}}(x,y,{{d}_{2}})$) и в результате получить выражение для вектора ${\mathbf{H}}({\mathbf{r}})$ только через его нормальную составляющую ${{H}_{n}}({\mathbf{r}})$ на поверхности $S$ дефекта ${{\Omega }_{d}}$. В следующем разделе получим систему уравнений относительно ${{\psi }_{1}}(x,y)$ и ${{\psi }_{2}}(x,y)$.

4. СИСТЕМА УРАВНЕНИЙ ДЛЯ НОРМАЛЬНЫХ СОСТАВЛЯЮЩИХ ПОЛЯ НА ГРАНИЦАХ ПЛАСТИНЫ

Перейдем к z-компоненте векторов в правой и левой части (3.3). При этом градиент перейдет в производную ${\partial \mathord{\left/ {\vphantom {\partial {\partial z}}} \right. \kern-0em} {\partial z}}$, после внесения которой под знак интеграла получим

где согласно (3.6) z-компонента вектора ${\mathbf{G}}({\mathbf{r}})$ имеет вид Перейдем в (4.1) к пределу по ${\mathbf{r}}$ изнутри ${{\Omega }_{1}}$ последовательно на ${{S}_{2}}$ (т.е. $z \to {{d}_{2}} - 0$), а затем на ${{S}_{1}}$ (т.е. $z \to - {{d}_{1}} + 0$), учитывая следующую общую формулу (см., например, [13, с. 265]): В результате получаем следующую систему уравнений относительно искомых функций ${{\psi }_{1}}(x,y)$ и ${{\psi }_{2}}(x,y)$: где введены обозначения С учетом (4.2) выражения для функций ${{G}_{1}}(x,y)$ и ${{G}_{2}}(x,y)$ имеют вид

5. ВЫВОД ФОРМУЛ ДЛЯ ФУРЬЕ-ОБРАЗОВ ФУНКЦИЙ ${{\psi }_{1}}(x,y)$ И ${{\psi }_{2}}(x,y)$

Поскольку система уравнений (4.3), из которой ищутся функции ${{\psi }_{1}}(x,y)$ и ${{\psi }_{2}}(x,y)$, содержит свертки искомых функций с функцией вида

(при $b{\text{ }} = {\text{ }}d{\text{ }}$ и $a{\text{ }} = $ 3/2), то, естественно, можно перейти к поиску фурье-образов этих функций ${{\hat {\psi }}_{1}}({{k}_{1}},{{k}_{2}})$ и ${{\hat {\psi }}_{2}}({{k}_{1}},{{k}_{2}})$, для которых в этом случае будут получаться линейные алгебраические уравнения. Отметим, что преобразование Фурье функции (5.1) имеет следующий вид (см., например, [14]): Мы используем следующую форму преобразования Фурье функции $f(x,y)$: Обратное преобразование Фурье, примененное к функции $\hat {f}({{k}_{1}},{{k}_{2}})$, приводит, как известно, к исходной функции $f(x,y)$: Применяя преобразование Фурье к обеим частям уравнений в (4.3), с учетом (5.2) и формул преобразования свертки, получаем систему уравнений для ${{\hat {\psi }}_{1}}({{k}_{1}},{{k}_{2}})$ и ${{\hat {\psi }}_{2}}({{k}_{1}},{{k}_{2}})$Решая эту систему, получаем Если ввести для выражения в квадратной скобке (3.3) обозначение то соотношение (3.3) представляется в виде Таким образом, поле ${\mathbf{H}}({\mathbf{r}})$ удовлетворяет соотношению (5.8), а в выражении (5.7) для $D({\mathbf{r}})$ у функций ${{\psi }_{1}}(x,y)$ и ${{\psi }_{2}}(x,y)$ известны лишь их преобразования Фурье (5.5) и (5.6). Поиск функций ${{\psi }_{1}}(x,y)$ и ${{\psi }_{2}}(x,y)$ через обратное преобразование Фурье (5.4) с последующей их подстановкой в (5.7) показался не оптимальным. Поэтому нашей дальнейшей целью является выражение $D({\mathbf{r}})$ в (5.7) именно через ${{\hat {\psi }}_{1}}({{k}_{1}},{{k}_{2}})$ и ${{\hat {\psi }}_{2}}({{k}_{1}},{{k}_{2}})$ с последующей подстановкой вместо них найденных выражений (5.5), (5.6).

6. ВЫРАЖЕНИЕ НАПРЯЖЕННОСТИ РЕЗУЛЬТИРУЮЩЕГО ПОЛЯ ЧЕРЕЗ ЕЕ НОРМАЛЬНОЕ ЗНАЧЕНИЕ НА ПОВЕРХНОСТИ ДЕФЕКТА

Для достижения поставленной цели используем тот факт, что для любой квадратично суммируемой в ${{R}^{2}}$ функции $f(x,y)$ и произвольных чисел p, q и b (b > 0) выполнено

где функции $g_{\alpha }^{b}$ и $\hat {g}_{a}^{b}$ представлены в (5.1) и (5.2). Формула (6.1) следует из свойства преобразования Фурье сохранять скалярное произведение функций (формула Планшереля) и известной формулы преобразования Фурье при сдвиге аргументов преобразуемой функции. Применяя формулу (6.1) для каждого из интегралов в (5.7), с учетом (5.2) получаем Подставляя сюда выражения для ${{\hat {\psi }}_{1}}({{k}_{1}},{{k}_{2}})$ и ${{\hat {\psi }}_{2}}({{k}_{1}},{{k}_{2}})$ из (5.5), (5.6), получаем где для компактности введены обозначения Учитывая разложение из (6.3) получаем где обозначено Применяя к правой части (6.5) формулу (6.1) (читая ее справа налево) при $a = 1{\text{/}}2$ с учетом (5.1) и (5.2) будем иметь Преобразуем выражение для $P_{1}^{{(m)}}(a)$. Подставив в (6.6) (при j = 1) выражение для ${{G}_{1}}(x,y)$ из (4.6), получим где Значение такого вида двойного интеграла приведено в [15, с. 487], а потому имеем Подставляя это значение в (6.7) с учетом a > 0, $z{\kern 1pt} {\text{'}} + {{d}_{1}} > 0$ (поверхность $S$ выше плоскости ${{S}_{1}}$), получаем Аналогично для $P_{2}^{{(m)}}(a)$ получается формула Подставляя эти выражения для $P_{1}^{{(m)}}(a)$ и $P_{2}^{{(m)}}(a)$ в (6.4), а затем полученные выражения для ${{P}_{j}}(a)$ в (6.2), получаем где введены следующие обозначения: Подставляя выражение (6.8) для $D({\mathbf{r}})$ в (5.8) с учетом (4.2), получаем для ${\mathbf{r}} \in {{R}^{3}}{\backslash }({{S}_{1}} \cup {{S}_{2}} \cup S)$ соотношение где Отметим, что правая часть в (6.18) не содержит искомой напряженности магнитного поля, а определяется только конфигурацией внешнего поля ${{{\mathbf{H}}}^{0}}({\mathbf{r}})$ и параметрами самой пластины (при отсутствии в нем дефекта). Из (6.18) видно, что для вычисления из этого соотношения напряженности результирующего поля ${\mathbf{H}}({\mathbf{r}})$ необходимо иметь выражение для нормальной составляющей ${{H}_{n}}({\mathbf{r}})$ этой напряженности на поверхности дефекта $S$. Далее получим уравнение для определения этой составляющей.

7. ВЫВОД УРАВНЕНИЯ ДЛЯ НОРМАЛЬНОЙ СОСТАВЛЯЮЩЕЙ НАПРЯЖЕННОСТИ НА ПОВЕРХНОСТИ ДЕФЕКТА

Для получения такого уравнения зафиксируем единичный вектор внешней нормали ${\mathbf{n}}$ в произвольной точке на поверхности $S$ и умножим (скалярно) на этот вектор обе части (6.18). Затем в полученном соотношении устремим ${\mathbf{r}} \in {{\Omega }_{{{\text{ }}1}}}$ к той точке поверхности $S$, в которой был взят вектор ${\mathbf{n}}$. Тогда, учитывая формулу предельного значения нормальной производной потенциала простого слоя [16, с. 267], получаем следующее уравнение для определения ${{H}_{n}}({\mathbf{r}})$:

Первый интеграл в левой части суть прямое значение нормальной производной потенциала простого слоя на поверхности $S$, ядро в этом интеграле имеет вид [16, с. 266] (в числителе косинус угла между указанными векторами).

8. ОКОНЧАТЕЛЬНЫЕ ФОРМУЛЫ. СРАВНЕНИЕ С ЧАСТНЫМ СЛУЧАЕМ

Таким образом, из (6.18) следует, что вне границ раздела сред ${\mathbf{r}} \in {{R}^{3}}{\backslash }({{S}_{1}} \cup {{S}_{2}} \cup S)$ напряженность результирующего поля ${\mathbf{H}}({\mathbf{r}})$ вычисляется по формуле

В этой формуле $F({\mathbf{r}})$ определяется в (6.13) через значения z-компонент $H_{z}^{0}(x,y, - {{d}_{1}})$ и $H_{z}^{0}(x,y,{{d}_{2}})$ напряженности внешнего поля на границах ${{S}_{1}}$ и ${{S}_{2}}$ пластины; параметры $\lambda $ и ${{\lambda }_{d}}$ определены в (4.4) и (6.19), а выражения для $R_{m}^{{(1)}}$, $R_{m}^{{(2)}}$, $R_{m}^{{(3)}}$, $R_{m}^{{(4)}}$ приведены в (6.9)–(6.12). Входящая в формулу (8.1) нормальная составляющая ${{H}_{n}}({\mathbf{r}})$ на поверхности дефекта $S$ предварительно определяется из уравнения (7.1).

Отметим, что аналогичные формулы были получены в [8] для частного случая, когда дефект является воздушной полостью (т.е. ${{\mu }_{d}} = 1$), а внешнее поле ${{{\mathbf{H}}}^{0}}({\mathbf{r}})$ является плоским и направленным параллельно границам ${{S}_{1}}$ и ${{S}_{2}}$ пластины (т.е. $H_{z}^{0}({\mathbf{r}}) \equiv 0$). Если в формулах (8.1) и (7.1) положить ${{\mu }_{d}} = 1$ и $H_{z}^{0}({\mathbf{r}}) \equiv 0$, то они в точности переходят в формулы, полученные в указанной выше работе, что является дополнительной проверкой выведенных соотношений.

Для получения из (8.1) и (7.1) расчетных формул для составления соответствующих компьютерных программ необходимо конкретизировать форму поверхности дефекта $S$. В упомянутой работе [8] такие формулы были выведены для сферической формы поверхности $S$ в указанном выше случае дефекта в виде полости и внешнего поля, параллельного плоским границам пластины. В этой же работе проведены соответствующие расчеты и построены графики зависимости компонент напряженности результирующего поля от геометрических параметров пластины и дефекта в ней.