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

ВВЕДЕНИЕ

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

При произвольном расположении частиц в пространстве анализ физического механизма изменения скорости электромагнитного поля (света) затруднен из-за сложного поведения дифракционного поля на совокупности частиц, поэтому механизм изменения скорости света за счет многократного рассеяния будем исследовать в средах, в которых частицы расположены регулярно. В зависимости от взаимного расположения частиц в среде с включениями в процессе распространения электромагнитных волн могут образовываться направляющие системы различного типа (см. [6]). В средах с низкой концентрацией частиц в качестве направляющей системы удобно использовать бесконечный волновод прямоугольного поперечного сечения с абсолютно проводящими стенками. На продольной оси z волновода в интервале –L ≤ z ≤ QL располагаются центры Q частиц сферической формы. Используя метод отражений от стенок волновода и условие периодичности по поперечным координатам, задачу распространения электромагнитных волн в волноводе можно свести к задаче дифракции волн на бесконечной Q-слойной двоякопериодической в поперечном направлении решетке частиц. Решение задачи дифракции поля на решетке будет моделировать процесс распространения электромагнитного поля внутри достаточно протяженного облака разреженных частиц в области вне окрестности края облака.

1. ЭЛЕКТРОМАГНИТНОЕ ПОЛЕ В ВОЛНОВОДЕ

Введем  в  трехмерном  пространстве  декартову  систему  координат  и  выберем в качестве модельной направляющей системы бесконечный вдоль оси z волновод W = $ = \;\{ (x,y,z): - {{l}_{1}} < x,y < {{l}_{1}},\; - {\kern 1pt} \infty < z < \infty \} $ квадратного поперечного сечения с идеально проводящими стенками. Исследуем метод многократного рассеяния на примере распространения электрических или магнитных волн в волноводе, на оси которого расположено конечное множество Q частиц сферической формы. Центры частиц расположены в интервале –L ≤ z ≤ QL. Если сравнить процесс распространения поля в волноводе, содержащем Q частиц, с подробно исследованным процессом распространения поля в волноводе, содержащем диэлектрический слой S = {(x, y, z): –l₁ ≤ x, y ≤ l₁, –L ≤ z ≤ QL} толщиной (Q + 1)L со специальными свойствами среды в слое (см. [7]), можно получить соотношения, определяющие скорость распространения поля в среде с частицами. Основное отличие процесса распространения электромагнитных волн от акустического случая [8] заключается в том, что при распространении электромагнитных волн не существует однородной плоской волны, распространяющейся как в свободном пространстве, так и в волноводе (см. [6]). Для волновода с частицами рассмотрим спектр излучения, для которого волновой размер элементов решетки (частиц), то есть произведение характерного размера частицы на волновое число электромагнитного поля, достаточно мал.

Пространство внутри волновода W будем характеризовать параметрами ε и μ, ε-диэлектрическая, μ-магнитная проницаемости. Для вакуума параметры ε = μ = 1. На оси волновода размещены частицы сферической формы D_q , q = 1, …, Q. Радиус частицы D_q равен а < l₁, L, ее центр расположен в точке с координатами (0, 0, (q – 1)L), так что начало декартовой системы совпадает с центром первой частицы. С целью упрощения решения считаем поверхность частицы D_q идеально проводящей. В гауссовой системе измерения будем решать задачу дифракции электромагнитного поля, распространяющегося в волноводе, на частицах D_q. Поскольку центры частиц расположены регулярно, можно найти не конфигурационное среднее, а определяемое точностью метода вычисления значение скорости распространения электромагнитного поля в зависимости от концентрации, размера частиц и толщины слоя, содержащего частицы.

Обозначим через

${{\bar {D}}_{q}}$ – замыкание области D_q. В области D из полупространства z < 0 на конечное множество частиц D_q падает электромагнитное поле. Представим полное поле в виде суммы падающего поля и рассеянной части. Векторы напряженности полного электрического E и магнитного H поля в области D удовлетворяют уравнениям Максвелла где ω – круговая частота, с – скорость распространения света в вакууме. Из уравнений (1.1) следует, что в кусочно-однородной изотропной среде без локальных источников векторы E и H удовлетворяют векторным уравнениям Гельмгольца (см. [6]) здесь k – волновое число поля в области D. На стенке волновода и на поверхностях S_q шаров D_q, q = 1, …, Q, должно выполняться условие обращения в ноль касательной составляющей векторного поля E. Рассеянная часть векторов E и H должна быть ограничена по модулю в области D при Jmk > 0. Особенность распространения электромагнитных волн в волноводах заключается в том, что существует критическая длина волны, больше которой волна в волноводе не распространяется. Действительно, отыскивая частное решение уравнений (1.2) в виде где h – неопределенная константа, легко обнаружить, что $h = {{h}_{{nm}}} = \sqrt {k_{{}}^{2} - g_{{nm}}^{2}} ,$ g_nm – собственные значения задачи Дирихле или Неймана для уравнения Лапласа в поперечном сечении волновода, целые числа n,m или больше нуля (задача Дирихле) или одновременно в нуль не обращаются (задача Неймана) в случае электрических или магнитных полей соответственно. Для распространяющихся в волноводе волн $h > 0$. Если μ*, ε* – произвольные константы, характеризующие среду с включениями, то для вычисления этих констант среды нужно исследовать задачу распространения как электрической волны в волноводе (H_z-компонента равна нулю), так и задачу распространения магнитной волны (Е_z-компонента равна нулю). Рассмотрим далее случай, когда материал среды и частиц не ферромагнетик, т.е. с высокой точностью можно положить µ* = 1.0. Поскольку скорость света в среде с частицами определяется только диэлектрической проницаемостью, достаточно исследовать задачу дифракции в волноводе или электрических, или магнитных волн. Далее исследуем магнитные волны, для которых E_z-компонента электромагнитного поля равна нулю. Пусть распространяющаяся из бесконечности при z < 0 в области D магнитная волна имеет вид Множитель exp(–iωt) далее опускается. Если выписать представление для поперечных составляющих электромагнитного поля при произвольных ε и µ через продольные составляющие с учетом (1.3) то из формул (1.4) следует, что падающее поле в волноводе содержит составляющие $H_{x}^{0},{\text{ }}E_{y}^{0},{\text{ }}H_{z}^{0}$. Остальные компоненты поля равны нулю. Потребуем, чтобы поперечный волновой размер волновода лежал в диапазоне ${{g}_{{10}}}{{l}_{1}} < k{{l}_{1}} < {{g}_{{20}}}{{l}_{1}}$, g₁₀, g₂₀ – первое и второе в порядке возрастания, отличные от нуля собственные значения задачи Неймана уравнения Лапласа, определенные в поперечном сечении волновода. В этом случае тип волны, распространяющейся в волноводе при z < –L , сохранится при z > QL, т.е. новых нормальных волн в волноводе в процессе дифракции на частицах не возникнет. Поскольку волновой размер частицы ka достаточно мал по постановке задачи, параметр kl₁ лежит в диапазоне ${{g}_{{10}}}{{l}_{1}} < k{{l}_{1}} < {{g}_{{20}}}{{l}_{1}}$, то параметр a/l₁ одного порядка малости с параметром ka. Поскольку величина (a/l₁)³ пропорциональна объемной концентрации среды, то процесс распространения электромагнитного поля в волноводе в указанном диапазоне частот действительно отвечает распространению поля в среде с малым содержанием частиц. Через U обозначим H_z-компоненту полного магнитного поля в волноводе и положим U = U₀+ U₁, U₀ = $H_{z}^{0}$, U₁ – рассеянная на частицах часть магнитного поля. Поле U₁ внутри волновода удовлетворяет скалярному однородному уравнению Гельмгольца (см. [9]) где k – волновое число в области D, ${{\left. {\partial {{U}_{1}}{\text{/}}\partial n} \right|}_{\Gamma }} = 0,$ Г – боковая поверхность волновода, n – нормаль к боковой поверхности, направленная внутрь области определения поля. Поле U₁ ограничено в области D при Jmk > 0. На поверхности сфер D_q должны выполняться условия Здесь ∂/∂n – производная по нормали к сфере S_q, направленная вне шара D_q.

В процессе рассеяния поля на частицах образуется не только дифракционное поле U₁, но и так называемая эффективная скорость распространения волн в среде с частицами. Эту скорость определим следующим образом. Пусть внутри волновода с параметрами среды ε, μ размещен слой S ={(x, y, z): –l₁ ≤ x, y ≤ l₁, –L ≤ z ≤ QL} диэлектрика с параметрами среды $\varepsilon _{s}^{'} = {{\varepsilon }_{s}} + {\text{4}}\pi {\text{i}}{{\sigma }_{{\text{s}}}}{\text{/}}\omega ,{\text{ }}{{\mu }_{{\text{s}}}},$ ${{\sigma }_{{\text{s}}}}$ – проводимость среды слоя. Материал среды слоя будем считать не ферромагнетиком, т.е. µ_s = 1.0. Из полупространства z < 0 на слой S падает стороннее поле в виде магнитной волны $H_{z}^{0} = H{\text{cos}}(\pi (x - {{l}_{{\text{1}}}}){\text{/2}}{{l}_{{\text{1}}}}){\text{exp}}(i{{h}_{{{\text{1}}0}}}z)$. Требуется найти амплитуду прошедшей через слой волны при z > QL. Будем считать, что поперечный волновой размер волновода, содержащего диэлектрический слой, лежит в диапазоне ${{g}_{{10}}}{{l}_{1}} < k{{l}_{1}} < {{g}_{{20}}}{{l}_{1}}$. Магнитное поле удовлетворяет вне S однородному уравнению Гельмгольца (1.5), внутри S уравнению (1.5), если ε и µ заменить на $\varepsilon _{s}^{'},\;{{\mu }_{{\text{s}}}}$, и однородным условиям Неймана на границе поперечного сечения. В процессе распространения волны $H_{z}^{0}$ в волноводе со слоем диэлектрика при z < 0 образуется отраженная волна, внутри слоя – стоячая волна с волновым числом $k_{s}^{2} = k_{0}^{2}\varepsilon _{s}^{'}{{\mu }_{s}}$, а вне слоя при z > QL возбуждается прошедшая волна. Амплитуды этих волн вычисляются из условий равенства касательных составляющих магнитного поля H_z на границах слоя при z = –L и z = QL. Решение этой задачи известно [6]. Прошедшая через слой волна имеет вид

Если потребовать, чтобы в волноводе при z > QL амплитуда H_zs прошедшего через слой S поля равнялась амплитуде H_zd падающего поля плюс поля дифракции на конечном множестве частиц D_q,

то при известных параметрах задачи ω, ε, μ = μ_s = 1.0, соотношения (1.7), (1.8) есть нелинейное уравнение для вычисления электрической проницаемости $\varepsilon _{s}^{'} = \varepsilon {\text{*}}$ или скорости распространения электромагнитного поля ${{c}_{s}} = {\text{c/}}\sqrt {\varepsilon {\kern 1pt} *} $, которая, вообще говоря, комплексна. Скорость распространения поля c_s в слое S назовем эффективной скоростью распространения электромагнитного поля в среде волновода, содержащей включения.

2. РЕШЕНИЕ ЗАДАЧИ ДИФРАКЦИИ

Чтобы вычислить эффективную скорость света, необходимо знать амплитуду магнитной волны H_z = U₀ + U₁ в волноводе при z > QL, когда на множество частиц D_q падает волна $H_{z}^{0} = {{U}_{0}}$. Функция U₁ должна быть ограничена по модулю в области D при Jmk > 0 и должна быть дважды непрерывно дифференцируемым внутри области D и непрерывным вплоть до границы области решением задачи (1.5), (1.6). Поскольку функция U₁ есть решение задачи (1.5), (1.6), ее можно представить в виде суммы потенциалов простого слоя для уравнения Гельмгольца (колебательных потенциалов по терминалогии Купрадзе) [10] со специальным ядром в виде функции Грина для бесконечного волновода квадратного поперечного сечения с идеально проводящими стенками:

где ${{\nu }_{{\text{q}}}}$ – неизвестная плотность потенциала, заданная на поверхности S_q. Символ ∑∑* означает, что целые числа n, m одновременно в нуль не обращаются. Функция Грина $G(\bar {x},{{\bar {\xi }}_{q}})$ ищется методом разделения переменных и обладает всеми свойствами колебательного потенциала (см. [10]), поскольку представляет собой сумму фундаментального решения уравнения Гельмгольца плюс непрерывное слагаемое. Подставим представление (2.1) в краевые условия (1.6). Устремим точку наблюдения $\bar {x}$ на поверхность S_p и воспользуемся свойствами колебательного потенциала простого слоя. В результате получим систему интегральных уравнений для вычисления плотностей ${{\nu }_{q}}$ вида где ∂/∂n₁ – производная по внешней к поверхности S_p нормали в точке ${{\bar {\xi }}_{{1p}}}$. Введем малый параметр $\delta = {\text{max}}(a{\text{/}}{{l}_{{\text{1}}}},\left| {ka} \right|,\left| {{{h}_{{{\text{1}}0}}}a} \right|)$. Исследования, проведенные в работе [8] показали, что главный член низкочастотной асимптотики достаточно хорошо описывает дифракционное поле в волноводе. Будем далее искать асимптотику плотностей ${{\nu }_{q}}$ по малому параметру δ, сохраняя за асимптотиками те же обозначения. Из теории дифракции на сфере известно, что для вычисления главного члена асимптотики при ka → 0 достаточно взять не более двух первых членов разложения плотностей по полной системе сферических функций, записанных в системе сферических координат, связанной с центром сферы S_q . Представление (2.1) удобно при исследовании поля вне области –L₁ < z < QL. Для вычисления асимптотики решения системы уравнений (2.2) при достаточно малом δ удобно другое представление решения задачи. Построим функцию Грина G бесконечного вдоль оси z прямоугольного волновода с условием на стенке ∂U/∂n = 0 методом отражений источника, расположенного в точке $(\xi ,\eta ,\varsigma ) \in W$, от стенок волновода. Эта функция имеет следующий вид: В результате для функции U₁ получим представление Правую часть формулы (2.3) можно интерпретировать как поле дифракции падающего поля U₀ на двоякопериодической по осям x,y многослойной по оси z решетке, элементами которой являются тела D_qnm с поверхностями S_qnm. Из условия периодичности следует, что в представлении (2.3) ${{\nu }_{{qnm}}} = {{\nu }_{{q00}}} = {{\nu }_{q}},$ если равенство понимать в том смысле, что в системе координат, связанной с центром сферы S_qnm, коэффициенты Фуре функции ${{\nu }_{{qnm}}}$ совпадают с коэффициентами Фурье функции ${{\nu }_{q}}_{{00}}$ в системе координат, связанной с центром сферы S_q00. Для выбранного представления функции Грина система интегральных уравнений (2.2) перепишется в форме где $({{\theta }_{{1p}}},{{\varphi }_{{1p}}})$ – сферические координаты с полюсом в центре сферы S_p. Удобство представления (2.3) при решении системы интегральных уравнений (2.4), (2.5) заключается в том, что можно применить известную процедуру переразложения сферических функций в разных системах координат (см. [13]). Разложение фундаментального решения в ряд по сферическим функциям определяется выражением $R = \sqrt {{{r}^{2}} + {{\rho }^{2}} - 2\rho r(\cos \theta \cos \theta {\kern 1pt} {\text{'}} + \sin \theta \sin \theta {\kern 1pt} {\text{'}}\cos (\varphi - \varphi {\kern 1pt} ')} ,$ ${{j}_{s}}(x),\;h_{s}^{{(1)}}(x)$ – сферические функции Бесселя и Ханкеля, $P_{s}^{t}(\cos \theta )$ – присоединенные полиномы Лежандра, которые определены в [14]. У некоторых авторов функции $P_{s}^{t}$ с таким же обозначением отличаются знаком. При вычислении интегралов в системе (2.4), (2.5) используем представление сферических функций в разных системах координат здесь $({{r}_{i}},{{\theta }_{i}},{{\varphi }_{i}})$ – сферические координаты в i-й системе координат, $({{r}_{j}},{{\theta }_{j}},{{\varphi }_{j}})$ – сферические координаты в j-й системе координат, $(r_{j}^{i},\theta _{j}^{i},\varphi _{j}^{i})$ – координаты j-го центра в i-й системе координат, Коэффициенты Клебша–Гордана $b_{\sigma }^{{(qpnm)}}$ определяются громоздкими, но простыми формулами, приведенными в [13].

Перепишем правую часть системы (2.4), (2.5) в виде

Для исследования влияния на поле объемной концентрации частиц в среде в формуле (2.8) сохранен член первого порядка малости по параметру a/l₁. Из формулы (2.8) следует, что асимптотику плотности ${{\nu }_{p}}$ следует искать в виде где a_ip, i = 1, 2, p = 1, …, Q – неопределенные коэффициенты, $({{\theta }_{p}},{{\varphi }_{p}})$ – сферические координаты с полюсом в центре сферы S_p. Подставим выражение (2.9) в систему интегральных уравнений (2.4), умножим левую и правую части (2.8) последовательно на $P_{i}^{1}(\cos {{\theta }_{{1p}}})\sin ({{\theta }_{{ip}}})\exp ( - i{{\varphi }_{{1p}}}),$ $i = 1,2$, и проинтегрируем по $d{{\varphi }_{{1p}}}d{{\vartheta }_{{1p}}}$, используя представления (2.7) сферических функций в разных системах координат. Интегрирование удобно проводить в системе координат, связанной с центром сферы S_p. Обозначая через получим алгебраическую систему для вычисления амплитуд асимптотик неизвестных плотностей ${{\nu }_{q}}$ и вида Правая часть алгебраической системы (2.10) задается выражением Матричные коэффициенты алгебраической системы (2.10) вычисляются по формулам где $j_{i}^{'}(ka),\;h_{i}^{{(1)}}{\kern 1pt} '(ka),\;i = 1,2$ – производная сферических функций Бесселя и Ханкеля по аргументу, $r_{{qnm}}^{j} = \sqrt {{{{(j - q)}}^{2}}L_{{}}^{2} + ({{n}^{2}} + {{m}^{2}})4l_{1}^{2}} $, $\cos \theta _{{qnm}}^{j} = (q - j)L{\text{/}}r_{{qnm}}^{j}$, $\sin \varphi _{{qnm}}^{j} = m{\text{/}}\sqrt {{{n}^{2}} + {{m}^{2}}} $, сферические координаты сферы S_qnm в системе координат, связанной с центром сферы S_p, ${{\delta }_{{qp}}}$ – символ Кронекера,

Из представления функции Ханкеля по большому аргументу следует, что ряды в формулах (2.11) сходятся. Заменим в представлении для коэффициентов системы (2.11) сферические функции Бесселя и Ханкеля их асимптотиками при малых ka (см. [13], [14]). Анализ матричных коэффициентов показывает, что при ka → 0 и конечных kl₁, kL в выражении для коэффициентов C_qp слагаемое при δ_pp порядка О(1), слагаемое при (1 – δ_qp) порядка О(k³a³) и менее. Этим слагаемым можно пренебречь по сравнению со слагаемым при δ_pp и правой частью. Таким образом, с точностью до О(k³a³) алгебраическая система (2.10) имеет диагональную матрицу, решение системы (2.10) равно

а плотность потенциала ${{\nu }_{{q{\text{ }}}}}$ вычисляется из соотношения

3. ВЫЧИСЛЕНИЕ СКОРОСТИ РАСПРОСТРАНЕНИЯ ЭЛЕКТРОМАГНИТНОГО ПОЛЯ

Поскольку падающая в волноводе нормальная волна и дифрагированная в области z > (2Q – 1)L волна, а также прошедшая через слой волна представляют собой магнитную волну типа H₁₀, то для сравнения поля U = U₀ + U₁, прошедшего в волноводе через множество частиц, и поля H_zs, прошедшего через слой S, нужно приравнять амплитуды этих нормальных волн. Для этого нужно выписать поле дифракции U₁ при z > (2Q – 1)L в основной системе координат, связанной с центром сферы S₁. Воспользуемся для поля U₁ представлением (2.1), в котором вместо плотности ν_j следует подставить ее асимптотику (2.9), вычисленную при решении системы (2.10). Пренебрегая экспоненциально затухающими членами, получим следующее выражение для амплитуды магнитной волны полного поля, прошедшего через область, содержащую множество Q частиц, расположенных на оси волновода

Условие $\sqrt {{{k}^{2}} - {{{(\pi {\text{/}}2{{l}_{1}})}}^{2}}} \ne 0$ означает, что электромагнитная волна с волновым числом k = π/2l₁ не может возбудиться в волноводе. В формуле (3.1) величины x₁ и x_{Q + 1} определены по формулам (2.12). Опустим проблемы, связанные со специальным случаем распространения волн через прозрачный слой и слой четверти длины волны [7], поскольку всегда можно выбрать другую частоту волны из рассматриваемого диапазона. Из представления (3.1) следует, что с точностью до О((a/l₁)³) амплитуда электромагнитного поля H_z, прошедшего через множество частиц в волноводе, равна амплитуде падающей волны $H_{z}^{0}$. В этом случае уравнение (1.7) удовлетворяется при ${{h}_{{{\text{1}}0}}} = h_{{10}}^{1}$, т.е. скорость света в среде с частицами равна скорости света в вакууме. С точки зрения явления дифракции полученный результат означает, что в исследуемой задаче при заданной геометрии в указанном диапазоне частот электромагнитного поля эффект многократного рассеяния на частицах пренебрежимо мал, поскольку поле прямой дифракции на сферах порядка О(δ³), а поле вторичной дифракции порядка О(δ⁴). Эффект изменения скорости света в среде с включениями следует ожидать в средах с высокой концентрацией частиц.