Excitation of a rectangular resonator through communication windows in the conveyor installation of microwave heating

Full Text

Введение

Теория возбуждения резонаторов в форме Вайнштейна [1; 2] или Маркова [3] широко используется для стационарных и нестационарных задач электродинамики и оптики (см., например, [4‒8]). По сути, эти теории используют построение резольвенты на основе базиса самосопряженного оператора [9]. В случае наличия диссипации собственные значения (частоты) перестают быть действительными, базисы – полными, а оператор Максвелла – самосопряженным [10]. При этом становится возможным появление непрерывного спектра (например, для открытых граничных задач). В случае учета реальных металлических стенок резонатора в [1] предложено продлить оболочку вглубь металла на глубину нескольких скин-слоев, где можно поставить нулевые граничные условия и использовать полный ортогональный базис нового резонатора с расширенной оболочкой. Для этого в [1] предложено использовать такие «расширенные» собственные функции, а собственные частоты считать комплексными и определять по методу возмущений, например, как в работе [11]. Это приближенный подход, а строгий подход требует учета слоя с комплексной диэлектрической проницаемостью металла, которая в модели Друде – Лоренца может быть весьма большой по модулю. В ряде задач необходимо учитывать диэлектрические тела внутри резонаторов, причем часто с существенными потерями, а иногда и с нелинейными свойствами. В этом случае строгая теория возбуждения резонатора формулируется в виде интегральных уравнений [4‒7].

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

В конвейерной установке СВЧ-нагрева НПП фирмы «Этна плюс» в качестве камеры использован отрезок многомодового желобкового волновода (ЖВ) [12; 13]. Такие секции ЖВ соединены в периодическую структуру посредством щели подачи продукта на транспортере и расположены горизонтально. В конструкции используются две секции фильтров на прямоугольных волноводах (ПВ), подсоединенные на входном и выходном окне для ввода транспортера с продуктом, а также конические трансформаторы типов мод с моды H10 ПВ к моде H11 ЖВ. Трансформаторы нужны для подключения магнетронных источников излучения с коаксиально-волноводными переходами (КВП) и с выходом на ПВ. Трансформаторы мод подключены к торцам секций ЖФ по два на каждую [14; 15]. Для увеличения вводимой мощности и равномерности нагрева рассмотрен вариант дополнительного подключения источников к верхнему желобу. Рассмотренная конструкция достаточно крупногабаритная и металлоемкая.

Поскольку габариты указанной установки весьма велики, в данной работе рассматривается возможность снижения в несколько раз металлоемкости и габаритов резонаторной камеры с несколькими источниками с использованием теории возбуждения резонатора. Для этого камера предполагается выполненной как прямоугольный резонатор с размерами a, b, c, имеющий вид длинного отрезка прямоугольного волновода (рис. 1). В двух противоположных малых стенках $(a\times b)$ резонатора при z = 0 и z = c прорезаны щели, через которые входит транспортер с продуктом. Длина транспортера в камере и максимальный ее размер c порядка нескольких метров, как и в описанной выше установке. Щели переходят в два заграждающих фильтра, препятствующих излучению из камеры. Фильтры выполнены в виде волновода с диафрагмами, имеющими поглощающие элементы. Магнетронные источники питания подключаются непосредственно на боковые стенки камеры рис. 1 через окна 1‒6 на стенке x = 0 и симметричные им окна на стенке x = a. Также такие же источники подключаются на верхней стенке камеры (окна 1‒4). Внутри камеры расположена лента транспортера с диэлектрическим продуктом, имеющим максимальную толщину t порядка 3 см. Высота узкой стенки b = 10‒15 см (в расчетах использовано значение b = 10 см), размер широкой стенки a = 90 см, размер верхней и нижней стенок $a\times c$, причем длина камеры c может меняться в широких пределах до нескольких метров (в зависимости от числа подключаемых источников). Для подключения используются стандартные ПВ $a_0=90$ мм $b_0=45$ мм с магнетронными источниками на частоте 2,45 ГГц, включенными через КВП.

 

Рис. 1. Схематический вид конвейерной камеры СВЧ-нагрева: 1 – отрезок волновода с продуктом; 2, 3 – окна ввода энергии через боковые стенки (номера окон: 1, 2, 3, 4, 5, 6); 4 – устройства ввода энергии через верхнюю стенку (КВП) (1, 2, 3, 4 – окна связи на верхней стенке)

Fig. 1. Schematic view of the microwave heating conveyor chamber: 1 – waveguide section with the product; 2, 3 – windows for energy input through the side walls (window numbers: 1, 2, 3, 4, 5, 6); 4 – devices for introducing energy through the upper wall (TWT) (1, 2, 3, 4 – communication windows on the upper wall)

 

Целью работы является получение простой и строгой электродинамический модели, пригодной для инженерного расчета, включая расчет КПД нагрева, температуры нагрева, нелинейных свойств, интегральных параметров. Строгая электродинамическая модель анализа полей в рабочей камере с прямоугольными границами может быть реализована на основе метода интегральных уравнений [4]. Абсолютно строгий подход требует как расчета полей в многомодовом режиме в самой камере (резонаторе), так и в питающих его отрезках волноводов с учетом отражений от магнетронов и потерь в стенках. Описывая магнетрон и КВП эквивалентной схемой активного двухполюсника с неким резонансным контуром, определяющим добротность, и заданным активным нелинейным элементом (например, источником тока, зависящим от напряжения на контуре), такую модель формально можно записать, но реализовать в виде программы достаточно сложно. В коммерческих программных пакетах типа HFSS, CST и в ряде других реализуется линейный режим задания амплитуды падающей основной волны полубесконечного ПВ без отражений от источника, но с учетом отражения от резонатора. Более сложная модель – задание коаксиальных портов на КВП и рассмотрение полей внутри КВП. Даже для прямоугольной камеры СВЧ-нагрева с несколькими окнами такие программы считают распределение полей достаточно долго, тем более это касается оптимизации. Они не позволяют оценивать КПД нагрева, рост температуры со временем, интегральный баланс реактивной и активной мощностей, нелинейные эффекты и ряд других интегральных и локальных параметров. Расчет усложняется при наличии рупорных трансформаторов мод ПВ к ЖВ. Указанные трансформаторы важны при согласовании в режиме бегущей волны, однако реально имеет место резонаторный режим низкодобротного резонатора с движением энергии к продукту. В этом случае целесообразно оценить возбуждение без трансформаторов непосредственно через окна связи. При включении нескольких окон и большой длине (рис. 1) щели в виде окна подачи продукта влияют слабо, излучение через них в силу наличия продукта и фильтров мало, поэтому можно считать их отсутствующими. Излучение в фильтр на входе продукта (где больше влажность) существенно меньше, чем на выходе. Фильтры рассчитаны на подавление порядка 100 дБ, что также позволяет считать окна закороченными. Таким образом, конвейерная установка моделируется как резонатор с диэлектрическим слоем и идеально проводящими стенками с несколькими боковыми и верхними прямоугольными окнами, на которых задан поток мощности.

1. Постановка задачи и формулы возбуждения

Рассмотрим прямоугольную камеру с металлическими стенками (прямоугольный параллелепипед) с размерами a, b, c. Пусть в его двух боковых и в верхней стенке прорезаны прямоугольные окна, запитываемые ПВ одинакового поперечного сечения. Конструктивно удобно подключать КВП. Можно поворачивать ПВ на одной стенке на 90 градусов, изменяя поляризацию вводимого излучения, но это усложняет модель, не принося в данном случае реального выигрыша. Магнетроны имеют полосу генерации (нестабильность частоты) порядка 1 % и не фазированы, что существенно улучшает равномерность поля. На двух гранях $a\times b$ реализован ввод продукта. Считаем ширину продукта равной $a$ (реально лента и продукт несколько меньше), а его толщину – равной t. На двух гранях $b\times c,$ а также на верхней грани $a\times с$ имеются периодически расположенные окна. Можно запитывать окна синфазно, противофазно и раздельно во времени относительно каких-то других окон. Однако для расчета мощностных характеристик удобно считать ввод энергии одновременным с одинаковыми фазами. Считаем, что на каждом из окон задан бегущий режим волновода с одной основной модой. Окна расположены так, что электрическое поле в них имеет компоненту $E_z$ (вдоль движения продукта). Касательное электрическое поле запитываемого от магнетрона окна есть заданный поверхностный магнитный ток на электрической стенке камеры. Этот ток возбуждает резонатор, при этом мощность, втекающая в него, есть интеграл от вектора Пойнтинга $S=\mathrm{Re}(E\times H^{*})/2$ по поверхности окна. Для окна на боковой поверхности это втекающая плотность мощности $S_x=Z_{10}{\left|E_z\right|}^2/2.$ Здесь $Z_{10}=\sqrt{{\mu }_0/{\epsilon }_0}/\sqrt{1-{\left(\pi /k_0{a}_0\right)}^2}$ ‒ волновой импеданс основной моды ПВ. В локальной системе координат окна имеем $E_z\left(y\right)=E_0\sin \left(\pi y/a_0\right).$ Задавая вводимую мощность, мы определяем $E_0$ и тем самым компоненту поверхностного магнитного тока $j_y^h\left(y,z\right)=E_z\left(y\right)=E_0\sin \left(\pi y/a_0\right).$ Формально можно ввести объемный магнитный ток с плотностью $J_y^h\left(x,y,z\right)=j_y^h\left(y\right)\delta \left(x\right).$ Аналогично на верхней стенке вводятся магнитные поверхностные токи $j_x^h\left(x,z\right)=-E_z\left(x\right)=-E_0\sin \left(\pi x/a_0\right).$ Модель можно модифицировать, когда часть окон неактивны. Окно с подключенным, не работающим в данный момент магнетроном есть импедансная поверхность с поверхностным импедансом, удовлетворяющим соотношению $E_{\mathbf{\tau }}\mathbf{=}{\mathbf{Z}}_{\mathbf{s}}\mathbf{\nu }\mathbf{\times }\mathbf{H}\mathbf{.}$ Энергия через нее как бы вытекает из резонатора и поглощается в резонаторе магнетрона. Однако поскольку добротность последнего существенно выше, можно считать этот импеданс реактивным и равным импедансу короткозамкнутого отрезка ПВ в КВП. Его влияние мало, и им можно пренебречь. Мы не рассматриваем режим, когда какие-либо из окон не работают на излучение.

Для решения задачи применим метод возбуждения резонатора в сочетании с теоремой Пойнтинга [2]. Теория возбуждения резонатора основана на представлении полей в виде

$E\left(r,\omega \right)=\sum _{k=1}^{\infty }{\alpha }_k\left(\omega \right)E_k\left(r,{\omega }_k\right),$ (1)

$H\left(r,\omega \right)=\sum _{k=1}^{\infty }{\beta }_k{H}_k\left(r,{\omega }_k\right).$

Удобно вводить только соленоидальные поля полого резонатора, считая, что диэлектрик эквивалентен некому заданному току поляризации $J^{\mathbf{p}}\left(r,\omega \right)\mathbf{=}{\mathbf{i\omega \epsilon }}_{\mathbf{0}}\left(\epsilon -1\right)\mathbf{E}\left(r,\omega \right)\mathbf{,}$ также возбуждающему поле. Потери в стенках считаем отсутствующими. Тогда поля (1) полны, ортогональны и удовлетворяют однородным уравнениям Максвелла

$\nabla \times E_k\left(r,{\omega }_k\right)=-i{\mu }_0{\omega }_k{H}_k\left(r,{\omega }_k\right),$ (2)

$\nabla \times H_k\left(r,{\omega }_k\right)=i{\epsilon }_0{\omega }_k{E}_k\left(r,{\omega }_k\right).$

Введем полную плотность магнитного тока

$J^{\mathbf{h}}\left(r,\omega \right)\mathbf{=}{\underset{j}{\sum J}}_j^h\left(r,\omega \right)\mathbf{,}$

где в сумме учтены все плотности, создаваемые поверхностными источниками в виде касательных электрических полей на окнах связи рис. 1. Уравнения Максвелла приобретают вид

$\nabla \times E\left(r,\omega \right)=-i{\mu }_0\omega H\left(r,\omega \right)-J^h\left(r,\omega \right),$ (3)

$\nabla \times H\left(r,\omega \right)=i{\epsilon }_0\omega E\left(r,\omega \right)+J^p\left(r,\omega \right).$

Подставляя в них разложения (2) и пользуясь ортогональностью функций по объему резонатора, выразим коэффициенты разложения. В результате получим представления полей через токи в виде функций Грина

$E\left(r,\omega \right)=\underset{V}{\int }\left[{\widehat{\Gamma }}_{\omega }^{ee}{\left(r,r^{\text{'}}\right)}^p\left(r^{\text{'}}\right)+{\widehat{\Gamma }}_{\omega }^{eh}\left(r,r^{\text{'}}\right)J^h\left(r^{\text{'}}\right)\right]d{r}^{\text{'}},$ (4)

$H\left(r,\omega \right)=\underset{V}{\int }\left[{\widehat{\Gamma }}_{\omega }^{he}{\left(r,r^{\text{'}}\right)}^p\left(r^{\text{'}}\right)+{\widehat{\Gamma }}_{\omega }^{hh}\left(r,r^{\text{'}}\right)J^h\left(r^{\text{'}}\right)\right]d{r}^{\text{'}}.$ (5)

Четыре тензорные функции Грина приведены, например, в работах [4‒8; 16]. Они построены из собственных функций резонатора (содержат диады типа $E_{\mathbf{k}}\left(r,{\omega }_k\right)\mathbf{\otimes }{\mathbf{E}}_{\mathbf{k}}^{\mathbf{*}}\left(r,{\omega }_k\right)\mathbf{,}$ $H_{\mathbf{k}}\left(r,{\omega }_k\right)\mathbf{\otimes }{\mathbf{H}}_{\mathbf{k}}^{\mathbf{*}}\left(r,{\omega }_k\right)\mathbf{,}\mathbf{ }{\mathbf{E}}_{\mathbf{k}}\left(r,{\omega }_k\right)\mathbf{\otimes }{\mathbf{H}}_{\mathbf{k}}^{\mathbf{*}}\left(r^{\text{'}},{\omega }_k\right)\mathbf{,}\mathbf{ }{\mathbf{H}}_{\mathbf{k}}\left(r,{\omega }_k\right)\mathbf{\otimes }{\mathbf{E}}_{\mathbf{k}}^{\mathbf{*}}\left(r^{\text{'}},{\omega }_k\right)\mathbf{)}\mathbf{.}$

Собственные функции прямоугольного резонатора известны и могут быть определены как E-моды и H-моды ПВ относительно оси z [2; 3; 11]. Собственные частоты даются формулой ${\omega }_{lmn}={\omega }_k={\left({\epsilon }_0{\mu }_0\right)}^{-1/2}\sqrt{{\left(l\pi /a\right)}^2+{\left(m\pi /b\right)}^2+\left(n\pi /c\right)},$ с учетом диапазона изменения индексов. Уравнения (4) и (5) перепишем в виде

$E\left(r,\omega \right)=E_0\left(r,\omega \right)+\underset{V_d}{\int }{\widehat{\Gamma }}_{\omega }^{ee}\left(r,r^{\text{'}}\right)J^p\left(r^{\text{'}},\omega \right)d{r}^{\text{'}},$ (6)

$H\left(r,\omega \right)=H_0\left(r,\omega \right)+\underset{V_d}{\int }{\widehat{\Gamma }}_{\omega }^{he}\left(r,r^{\text{'}}\right)J^p\left(r^{\text{'}},\omega \right)d{r}^{\text{'}}.$ (7)

Здесь выделенные индексом ноль поля созданы поверхностными источниками, т. е. они известны. Эти поля требуют однократного определения. Интегралы дают вклад от тока поляризации. Эти уравнения являются интегральными, поскольку неизвестное электрическое поле входит под интегралы. Решение (6), (7) следует искать только в области диэлектрика V_d. Поскольку она мала по сравнению с объемом резонатора, весьма точным является приближенное решение с подстановкой вместо поля E под интегралы выделенного поля:

$\begin{array}{l}E\left(r,\omega \right)=E_0\left(r,\omega \right)+\\ +i\omega \left(\epsilon -1\right)\underset{V_d}{\int }{\widehat{\Gamma }}_{\omega }^{ee}\left(r,r^{\text{'}}\right)E_0\left(r^{\text{'}},\omega \right)d{r}^{\text{'}},\end{array}$ (8)

$\begin{array}{l}H\left(r,\omega \right)=H_0\left(r,\omega \right)+\\ +i\omega \left(\epsilon -1\right)\underset{V_d}{\int }{\widehat{\Gamma }}_{\omega }^{he}\left(r,r^{\text{'}}\right)E_0\left(r^{\text{'}},\omega \right)d{r}^{\text{'}}.\end{array}$ (9)

Уравнения (8), (9) для расчета полей явные. Мощность потерь определяется выражением $P_{\epsilon }\left(r,\omega \right)=\mathrm{Re}\left(E\left(r,\omega \right)J^{P*}\left(r,\omega \right)\right)/2$

или $P_{\epsilon }\left(r,\omega \right)=\omega {\epsilon }_0{\epsilon }^{\text{'}\text{'}}\left(r,\omega \right){\left|E\left(r,\omega \right)\right|}^2/2.$

Заданное поле E⁰ и H⁰ в (8) и (9) можно записать как поверхностные интегралы от известных поверхностных источников и легко вычислить, используя при этом разложения (1) и следующие формулы для коэффициентов:

${\alpha }_k^{\left(0\right)}=\frac{-i}{N_k\left({\omega }^2-{\omega }_k^2\right)}{\omega }_k\underset{S}{\oint }{j}^h{H}_{k}dS,$ (10)

Диэлектрический продукт описываем поляризацией $P_{\epsilon }={\epsilon }_0\left(\epsilon \left(\omega \right)-1\right)E,$ которая при возбуждении резонатора функционирует как поглощающий энергию ток поляризации с плотностью $J^{\mathbf{p}}\left(r,\omega \right)\mathbf{=}{\mathbf{i\omega \epsilon }}_{\mathbf{0}}\left(\epsilon \left(\omega \right)-1\right)\mathbf{E}\mathbf{.}$

Рассмотрим интегральный баланс мощности в виде комплексной теоремы Пойнтинга:

$\begin{array}{l}\mathrm{Re}\left(\sum _n\underset{S_n}{\int }{S}_n\nu dS-\sum _m{\underset{S_m}{\int }}_m\nu dS-\right.\\ \left.-\sum _{j=\mathrm{1,2}}\underset{S_d}{\int }{\text{S}}_{d_j}\nu dS-\underset{V}{\int }E{J}^{p^{*}}dV\right)/2=0.\end{array}$ (11)

Здесь индексом n обозначены все прямоугольные окна (и соответствующие векторы Пойнтинга), через которые мощность втекает, индексом m обозначены оставшиеся окна, поглощающие мощность (к ним относятся окна с неработающими магнетронами), а индексом j – окна подачи продукта). Нормаль $\nu$ считаем направленной вдоль потока мощности (т. е. внутренней при втекании и внешней при вытекании). В первом приближении поглощающие окна можно не учитывать, поскольку поглощение в них пренебрежимо мало по сравнению с поглощением продукта. Поток мощности $S_{{\mathbf{d}}_{\mathbf{j}}}$ поглощается диэлектриком в окне его подачи. Считаем, что канал подачи между двумя секциями достаточно протяженный, и попавшая в него мощность поглощается, не доходя до следующей секции. В этом случае можно использовать квазипериодическую модель и рассчитывать секцию из одного периода, что упрощает модель в случае длинной камеры. В (11) взята внутренняя нормаль и рассматривается втекающая мощность. Моды прямоугольного резонатора приведены в работе [11]. Будем их рассматривать как E_mnl- и H_mnl-моды (ТМ и ТE) относительно оси z. Индексы l_mn обозначим мультииндексом k. Они пробегают значения $l,m=\mathrm{1,2,3...,}n=\mathrm{0,1,2...}$ для первых и $l,m=\mathrm{0,1,2...}$ (за исключением $l=m=0)$ и $n=\mathrm{1,2,3,...}$ для вторых [11]. Соответствующие комплексные поля будем писать как $E_{\mathbf{k}}^{\mathbf{e}}\mathbf{,}\mathbf{ }{\mathbf{H}}_{\mathbf{k}}^{\mathbf{e}}$ и $E_{\mathbf{k}}^{\mathbf{h}}\mathbf{,}\mathbf{ }{\mathbf{H}}_{\mathbf{k}}^{\mathbf{h}}\mathbf{.}$ Считая толщину продукта малой по сравнению с размером b, можно пренебречь потенциальными полями и оставить только соленоидальные. Тогда не нужно решать уравнение Пуассона для наведенных на поверхности диэлектрика зарядов и связанных объемных зарядов (в случае неоднородного диэлектрика) и сразу записать полные поля в виде разложений по соленоидальным базисам

$E=\sum _k\left(A_k^e{E}_k^e+A_k^h{E}_k^h\right), H=\sum _k\left(B_k^e{H}_k^e+B_k^h{H}_k^h\right).$   (10)

Отметим, что представления (1) и (10) для соленоидальных полей (полого резонатора) эквивалентны, поскольку колебания резонатора вырождены, а в (1) они предполагаются невырожденными и ортонормированными. Подход к задаче возбуждения с использованием потенциальных подсистем сложно реализуем [2; 4; 10]. Далее использован метод, не требующий нахождения потенциальных полей. В нем влияние диэлектрика рассмотрено как ток поляризации, что является строгим подходом и не требует учета потенциальных подсистем базисных функций.

ДП влажного продукта запишем в виде $\epsilon \left(\omega \right)=\tilde{\epsilon }\left(1-c_w\right)+c_w{\epsilon }_w\left(\omega \right),$

где $\tilde{\epsilon }$ ‒ ДП сухого продукта, которую считаем положительной константой больше единицы, $c_w$ ‒ положительный коэффициент влагосодержания $(0\le c_w<1), {\epsilon }_w\left(\omega \right)$ ‒ ДП воды. Пусть $k_S$ ‒ соленость воды (концентрация солей в граммах на литр). Для морской воды соленость 3,47 % (в промилях 34,7 ‰), т. е. в литре морской воды содержится примерно 35 граммов солей. Считая, что вся соль – это NaCl, получаем 0,6 моль/литр, т. е. при диссоциации в литре содержится $\mathrm{3,6}\cdot {10}^{23}$ зарядов. Поэтому для квадрата плазменной частоты получим ${\omega }_p^2=\mathrm{2,5}\cdot {10}^{23}$ Гц². В обрабатываемых продуктах соленость на 1‒3 порядка меньше. Поскольку вода – весьма плотно упакованная жидкость, частота столкновений для водных растворов обычно порядка плазменной. Используя формулу Дебая и формулу Друде, для воды получим

${\epsilon }_w\left(\omega \right)=1+\frac{\kappa }{1+{\left(\omega \tau \right)}^2}-i\frac{\kappa \omega \tau }{1+{\left(\omega \tau \right)}^2}+k_s\frac{{\omega }_p^2}{{\omega }^2-i\omega {\omega }_c}.$ (11)

Здесь $k_s$ ‒ величина безразмерная и равная отношению веса солей к весу воды. Формула (11) удобна для моделирования ДП при сушке, если ее использовать с коэффициентом влагосодержания $k_w$ и ДП ${\epsilon }_d$ сухого продукта. Тогда $\epsilon \left(\omega \right)=1+\left({\epsilon }_d-1\right)\left(1-k_w\right)+k_w\left({\epsilon }_w\left(\omega \right)-1\right).$

В процессе движения продукта он высыхает, т. е. влагосодержание $k_w\left(z\right)$ уменьшается и становится зависимым от продольной координаты. Определить указанную зависимость можно, зная скорость движения, выделяемую плотность мощности, удельную теплоемкость воды и теплоту параобразования. Для начального определения плотности мощности в первом приближении берем влагосодержание постоянным и равным входному $k_w\left(0\right).$ Далее можно использовать одну-две итерации. Таким способом рассмотренная модель была адаптирована к процессу сушки, для которого при заданной мощности для получения требуемого влагосодержания на выходе весьма важна скорость движения.

Итак, есть три пути расчета полей. Первый – пренебречь влиянием диэлектрика, т. е. взять поля с индексом нуль (нулевое приближение). Второй – применить метод возмущений в первом порядке (первое приближение) и рассчитать поля (8), (9). Третий – строго численно решать систему интегральных уравнений (6), (7). Поскольку коэффициенты разложений полей зависят от источников и стоков, которые сами зависят от полей, возникает задача типа интегральных уравнений или интегральных соотношений. Решить такую задачу можно на введенном базисе. Проекционный или вариационный подходы приводят к бесконечной системе линейных алгебраических уравнений (СЛАУ), разрешенной относительно неизвестных в форме $X=\widehat{M}X+Y.$ Здесь вектор X составлен из коэффициентов $A_k^{\left(e,h\right)}, B_k^{\left(e,h\right)}$ или ${\alpha }_k$ и ${\beta }_k.$ Источники определяются вектором Y. Решать такую систему удобно путем редукции итерационным методом, например методом минимальных невязок [17]. Для этого следует умножить систему уравнений на эрмитово сопряженную матрицу ${\left(\widehat{I}-\widehat{M}\right)}^H.$

В работе выбран третий подход, но с использованием формул возбуждения, что приводит к более простому алгоритму, чем матричный. Примененим разложение (1). Для окон сечения $a_0\times b_0$ в локальной координатной системе бокового окна имеем $E_z=E_0\sin \left(\pi y/a_0\right), H_y=Z_{10}{E}_0\sin \left(\pi y/a_0\right).$

Тогда нормальная компонента вектора Пойнтинга $S_x=Z_{10}{E}_0^2{\sin }^2\left(\pi y/a_0\right).$

Интегрируя по окну, получаем $P_n=Z_{10}{E}_0^2{b}_0{a}_0/4.$

Это мощность магнетрона, которая в среднем поступает в секцию. Здесь $Z_{10}=\sqrt{{\mu }_0/\left[{\epsilon }_0\left(1-{\pi }^2/{\left(k_0{a}_0\right)}^2\right)\right]},$ поэтому известна амплитуда $E_0.$ Собственные моды E_k и H_k имеют компоненты:

$E_{xk}=A_{lmn}\cos \left(k_{l}x\right)\sin \left(k_{m}y\right)\sin \left(k_{n}z\right),$ (12)

$E_{yk}=B_{lmn}\sin \left(k_{l}x\right)\cos \left(k_{m}y\right)\sin \left(k_{n}z\right),$

$E_{zk}=C_{lmn}\sin \left(k_{l}x\right)\sin \left(k_{m}y\right)\cos \left(k_{n}z\right),$

$H_{xk}={\tilde{A}}_{lmn}\sin \left(k_{l}x\right)\cos \left(k_{m}y\right)\cos \left(k_{n}z\right),$

$H_{yk}={\tilde{B}}_{lmn}\cos \left(k_{l}x\right)\sin \left(k_{m}y\right)\cos \left(k_{n}z\right),$

$H_{zk}={\tilde{C}}_{lmn}\cos \left(k_{l}x\right)\cos \left(k_{m}y\right)\sin \left(k_{n}z\right).$

В силу уравнений Максвелла $H_{\mathbf{k}}\mathbf{=}{\left(-{\mathrm{i\omega }}_k{\mu }_0\right)}^{\mathbf{-}\mathbf{1}}\mathbf{\nabla }\mathbf{\times }{\mathbf{E}}_{\mathbf{k}}$ коэффициенты подчинены соотношениям:

${\tilde{A}}_{lmn}=\frac{\left(k_m{C}_{lmn}-k_n{B}_{lmn}\right)}{-i{\omega }_k{\mu }_0},$ (13)

${\tilde{B}}_{lmn}=-\frac{\left(k_l{C}_{lmn}-k_n{A}_{lmn}\right)}{-i{\omega }_k{\mu }_0},$

${\tilde{C}}_{lmn}=\frac{\left(k_l{B}_{lmn}-k_m{A}_{lmn}\right)}{-i{\omega }_k{\mu }_0}.$

Поскольку имеет место вырождение колебаний, вводим комбинацию E-мод и H-мод, обозначая

$A_{lmn}=-\left(k_l{k}_n{C}_{lmn}-k_0{k}_m{C^{\text{'}}}_{lmn}\right)/\left(k_l^2+k_m^2\right)$ и $B_{lmn}=-\left(k_m{k}_n{C}_{lmn}+k_0{k}_l{C^{\text{'}}}_{lmn}\right)/\left(k_l^2+k_m^2\right).$

Как нетрудно видеть, такое представление удовлетворяет тождеству $\nabla \cdot E_k=\mathrm{0,}$ т. е. поля соленоидальные. Оставшиеся два произвольных коэффициента подчиняем условию ${C^{\text{'}}}_{lmn}=C_{lmn}$ и нормируем собственные колебания. Для нормировки берем условие

 $\begin{array}{l}{N}_k^e=1={\epsilon }_0{\underset{V}{\int }}_k^{2}dV=-N_k^h=-{\mu }_0{\underset{V}{\int }}_k^{2}dV=\\ =\frac{{\epsilon }_{0}abc}{8}\left[A_{lmn}^2\left(1+{\delta }_{l0}\right)+\right.\\ \left.+B_{lmn}^2\left(1+{\delta }_{m0}\right)+C_{lmn}^2\left(1+{\delta }_{n0}\right)\right],\end{array}$

откуда имеем

$\begin{array}{l}{C}_{lmn}^2=\frac{8}{{\epsilon }_{0}abc}\left[\left(1+{\delta }_{l0}\right){\left(k_l{k}_n-k_0{k}_m\right)}^2/{\left(k_l^2+k_m^2\right)}^2+\right.\\ {\left.+\left(1+{\delta }_{m0}\right){\left(k_m{k}_n+k_0{k}_l\right)}^2/{\left(k_l^2+k_m^2\right)}^2+\left(1+{\delta }_{n0}\right)\right]}^{-1}.\end{array}$

Норма пропорциональна энергии колебания. Теперь магнитные функции нормированы условием $N_k^h=-1$ (см. [2], стр. 357). Отрицательная норма соответствует наличию мнимой единицы в знаменателях соотношений (13) и означает переход в разные полупериоды электрической энергии колебаний в энергию магнитную и наоборот, т. е. сдвиг по фазе $\pi ).$ Тем самым мы получили невырожденные ортогональные системы соленоидальных векторных функций.

Рассмотрим два симметричных излучателя на боковых стенках (при x = 0 и x = a) с магнитной поверхностной плотностью тока $j_{sy}^h\left(y,z\right)=E_0\sin \left(\pi \left(y-y_0\right)/a_0\right)\Pi \left(z_s,z\right),$ т. е. с объемной плотностью $J_{sy}^h=j_{sy}^h\left(y,z\right)\left(\delta \left(x\right)+\delta \left(x-a\right)\right).$

Здесь $\Pi \left(z_s,z\right)$ ‒ прямоугольная функция (формфактор), равная единице в области $z_s\le z\le z_s+b_0$ и нулю в оставшейся. Индекс $s=\mathrm{1,2,...,}{N}_s$ выполняет нумерацию боковых диафрагм (для рис. 1 $N_s=6).$ Выполняя интегрирование с $H_{yk}\left(\mathrm{0,}y,z\right)$ и $H_{yk}\left(a,y,z\right),$ имеем вклад в коэффициенты

$\begin{array}{l}{\alpha }_{ks}^{\left(0\right)}=-i{E}_0{\tilde{B}}_{lmn}\left(1+{\left(-1\right)}^l\right){\omega }_k\times \\ \times \frac{\left[\sin \left(k_n\left(z_s+b_0\right)\right)-\sin \left(k_n{z}_s\right)\right]I_m}{2N_k^e{k}_n\left({\omega }^2-{\omega }_k^2\right)},\end{array}$ (14)

$\begin{array}{l}{\beta }_{ks}^{\left(0\right)}=i{E}_0{\tilde{B}}_{lmn}\left(1+{\left(-1\right)}^l\right)\omega \times \\ \times \frac{\left[\sin \left(k_n\left(z_s+b_0\right)\right)-\sin \left(k_n{z}_s\right)\right]I_m}{2N_k^e{k}_n\left({\omega }^2-{\omega }_k^2\right)},\end{array}$ (15)

где обозначены интегралы, приведенные в Приложении (П. 1). Рассмотрим теперь излучатель наверху. Для него $j_{px}^h\left(y,z\right)=E_0\sin \left(\pi \left(x-x_0\right)/a_0\right)\Pi \left(z_p,z\right),$ и вклад в коэффициенты разложения полей вычисляется как интеграл с $H_{xk}\left(x,b,z\right)$ в виде

$\begin{array}{l}{\alpha }_{kp}^{\left(0\right)}=i{E}_0{\tilde{A}}_{lmn}{\left(-1\right)}^m{\omega }_k\times \\ \times \frac{\left[\sin \left(k_n\left(z_p+b_0\right)\right)-\sin \left(k_n{z}_p\right)\right]I_l}{2N_k^e{k}_n\left({\omega }^2-{\omega }_k^2\right)},\end{array}$

$\begin{array}{l}{\beta }_{kp}^{\left(0\right)}=i{E}_0{\tilde{A}}_{lmn}{\left(-1\right)}^m\omega \times \\ \times \frac{\left[\sin \left(k_n\left(z_p+b_0\right)\right)-\sin \left(k_n{z}_p\right)\right]I_l}{2N_k^e{k}_n\left({\omega }^2-{\omega }_k^2\right)},\end{array}$

$I_l=2\underset{x_0}{\overset{x_0+a_0}{\int }}\sin \left(k_{l}x\right)\sin \left(\pi \left(x-x_0\right)/a_0\right)dx.$

Здесь $p=\mathrm{1,2,...,}{N}_p$ соответствует верхним окнам (в нашем случае $N_p=4),$ а интеграл дается той же формулой (П.1) $I_l=I\left(k_l,a_0,x_0\right).$ Теперь имеем коэффициенты разложений

${\alpha }_k^{\left(0\right)}=\sum _{s=1}^{N_s}{\alpha }_{ks}^{\left(0\right)}+\sum _{p=1}^{N_p}{\alpha }_{kp}^{\left(0\right)},{\beta }_k^{\left(0\right)}=\sum _{s=1}^{N_s}{\beta }_{ks}^{\left(0\right)}+\sum _{p=1}^{N_p}{\beta }_{kp}^{\left(0\right)},$

учитывающие вклад от всех окон ввода энергии. Нулевой порядок теории возмущений означает ${\alpha }_k={\alpha }_k^{\left(0\right)}, {\beta }_k={\beta }_k^{\left(0\right)}.$ Компоненты плотности тока поляризации имеют вид

$\begin{array}{l}{J}_x^p=i\omega {\epsilon }_0\left(\epsilon -1\right)\sum _k{\alpha }_k{A}_{lmn}\times \\ \times \cos \left(k_{l}x\right)\sin \left(k_{m}y\right)\sin \left(k_{n}z\right),\end{array}$ (16)

$\begin{array}{l}{J}_y^p=i\omega {\epsilon }_0\left(\epsilon -1\right)\sum _k{\alpha }_k{B}_{lmn}\times \\ \times \sin \left(k_{l}x\right)\cos \left(k_{m}y\right)\sin \left(k_{n}z\right),\end{array}$

$\begin{array}{l}{J}_z^p=i\omega {\epsilon }_0\left(\epsilon -1\right)\sum _k{\alpha }_k{C}_{lmn}\times \\ \times \sin \left(k_{l}x\right)\sin \left(k_{m}y\right)\cos \left(k_{n}z\right).\end{array}$

Множитель $\left(\epsilon -1\right)$ отличен от нуля только в области диэлектрика $\left(b-t\right)/2<y<\left(b+t\right)/2.$ Интегрируя по этой области, получаем поправки к полям, обусловленные заполнением. Указанные поправки определяются коэффициентами ${\alpha }_{kd}$ (П. 2) и ${\beta }_{kd}$ (П. 3). Уточненные поля определяются коэффициентами ${\alpha }_k^{\left(1\right)}={\alpha }_{kd}+{\alpha }_k^{\left(0\right)}, {\beta }_k^{\left(1\right)}={\beta }_{kd}+{\beta }_k^{\left(0\right)}.$ Полагая ${\alpha }_k={\alpha }_k^{\left(1\right)}, {\beta }_k={\beta }_k^{\left(1\right)},$ имеем первый порядок теории возмущений. Его можно использовать, чтобы снова пересчитать плотности тока поляризации (16) и коэффициенты (П. 2), (П. 3). При этом пересчитывать коэффициенты ${\alpha }_k^{\left(0\right)}, {\beta }_k^{\left(0\right)}$ не нужно, они вычисляются один раз, что определяет простоту метода. В случае не слишком большой величины $\epsilon -1$ и малого коэффициента заполнения $t/b\ll 1$ возмущение поля диэлектриком мало, и итерационный процесс сходится весьма быстро за 1‒2 итерации. Интегралы в (16) вычисляются весьма просто в виде однократных сумм:

$\begin{array}{l}{I}_{xk}=\underset{V_d}{\int }{J}_x^p{E}_{xk}dxdydz=\frac{i\omega {\epsilon }_0\left(\epsilon -1\right)ac}{8}\times \\ \times \left(1+{\delta }_{l0}\right)\sum _{m^{\text{'}}}{\alpha }_k{A}_{lmn}{A}_{l{m}^{\text{'}}n}{I}_{xm{m}^{\text{'}}},\end{array}$ (17)

$\begin{array}{l}{I}_{yk}=\underset{V_d}{\int }{J}_y^p{E}_{yk}dxdydz=\\ =\frac{i\omega {\epsilon }_0\left(\epsilon -1\right)ac}{8}\sum _{m^{\text{'}}}{\alpha }_k{B}_{lmn}{B}_{l{m}^{\text{'}}n}{I}_{ym{m}^{\text{'}}},\end{array}$

$\begin{array}{l}{I}_{zk}=\underset{V_d}{\int }{J}_z^p{E}_{zk}dxdydz=\frac{i\omega {\epsilon }_0\left(\epsilon -1\right)ac}{8}\times \\ \times \left(1+{\delta }_{n0}\right)\sum _{m^{\text{'}}}{\alpha }_k{C}_{lmn}{C}_{l{m}^{\text{'}}n}{I}_{zm{m}^{\text{'}}},\end{array}$

где обозначенные интегралы приведены в (П. 4), (П. 5).

2. Результаты и выводы

На рис. 2‒4 приведены результаты расчета квадрата поля |E|² в плоскости $y=b/2$ от x и z в камере длиной 2 м, шириной 0,9 м и высотой 0,1 м с продуктом толщины 3 см и $\epsilon =3-\mathrm{0,3}i$ при 6 окнах ввода энергии на каждой из боковых стенок и 4 окнах ввода сверху для линейной частоты 2,45 ГГц. Вводимая мощность через каждое окно полагалась равной 500 Вт, т. е. полная мощность равнялась 8 кВт. Результат сходится для трех итераций. В приведенном расчете использовано 11 вариаций по каждой из осей, т. е. более 1000 мод, что привело ко времени счета порядка нескольких секунд на персональном компьютере малой производительности. Моды высоких порядков вносят малый вклад, поэтому в расчетах можно в несколько раз снизить число используемых мод и существенно понизить время счета при сохранении порядка точности. Поле симметрично относительно плоскостей $x=a/2$ и $z=c/2.$ На рис. 4 также приведены результаты моделирования той же структуры в CST Microwave Studio с заданием портов в глубине отрезков волноводов. Как видно, имеют место хорошее совпадение и относительно равномерный нагрев. Также имеет место хорошее соответствие распределения ${\left|E\right|}^2$ в сечениях для рис. 2 и 3. Моделирование в CST занимает примерно час времени на компьютере высокой производительности. Предложенный метод удобен тем, что позволяет легко оптимизировать конструкцию, изменяя положение окон, поскольку не требует большого времени счета.

 

Рис. 2. Зависимость ${\left|E\right|}^{\mathbf{2}}$ (В² / м²) от координаты z (м) при разных значениях координаты x (м): 0,0 (кривая 1), 0,1 (2), 0,2 (3) и y = b / 2

Fig. 2. Dependence of ${\mathbf{\left|}\mathbf{E}\mathbf{\right|}}^{\mathbf{2}}$ (V² / m²) on the z coordinate (m) for different values of the x (m) coordinate: 0,0 (curve 1), 0,1 (2), 0,2 (3) and y = b / 2

 

Рис. 3. Зависимость ${\left|E\right|}^{\mathbf{2}}$ (В² / м²) от координаты x (м) при разных значениях координаты z (м): 0,25 (кривая 1), 0,3 (2), 1,0 (3), 1,75 (4), 1,7 (5) и y = b / 2

Fig. 3. Dependence of ${\mathbf{\left|}\mathbf{E}\mathbf{\right|}}^{\mathbf{2}}$ (V² / m²) on the x coordinate (m) for different values of the z coordinate (m): 0,25 (curve 1), 0,3 (2), 1,0 (3), 1,75 (4 ), 1,7 (5) and y = b / 2

 

Рис. 4. Распределение ${\left|E\right|}^{\mathbf{2}}$ (В² / м²) в плоскости x, z при y = b / 2 на основе теории возбуждения (сверху) и моделирование в пакете CST (снизу)

Fig. 4. Distribution of ${\mathbf{\left|}\mathbf{E}\mathbf{\right|}}^{\mathbf{2}}$ (V² / m²) in the x, z plane at y = b / 2 based on excitation theory (top) and CST modeling (bottom)

 

Результаты моделирования приводят к следующим выводам. При движении транспортера вдоль оси z каждая точка продукта проходит через локальные минимумы и максимумы (рис. 2), при этом среднее (усредненное) поле почти не меняется, что обеспечивает в среднем однородный нагрев. При некоторых координатах z вдоль оси x реализуется достаточно однородное распределение выделяемой мощности, а неравномерность нагрева на других линиях x компенсируется за счет движения, поскольку каждая точка проходит области минимумов и максимумов (рис. 3). В целом движение для распределения поля (рис. 4) обеспечивает достаточно равномерный нагрев. Реально температура при близко расположенных чередующихся многочисленных локальных минимумах и максимумах выравнивается за счет теплопроводности. Следует учесть разброс (нестабильность) порядка 1 % частот генерации магнетронов и случайность их фаз, что приводит к перестройке во времени распределений поля. Улучшить равномерность нагрева можно, механически меняя граничные условия, например вводя вентилятор с металлическими лопастями, улучшающий также конвекцию и вывод влажного воздуха.

Заключение

В работе получена модель возбуждения резонатора, имеющая ряд достоинств. Кроме весьма малого объема памяти и малого объема вычислений отметим следующие. Учет диэлектрика не требует пересчета коэффициентов ${\alpha }_k^{\left(0\right)}$ и ${\beta }_k^{\left(0\right)},$ а пересчет ${\alpha }_{kd}$ и ${\beta }_{kd}$ занимает существенно меньше ресурсов. Неравномерный нагрев может приводить к локальному изменению ДП в зависимости от мощности, например за счет расширения диэлектрика или уменьшения влажности. Эти явления легко учесть в алгоритме как неоднородную ДП, при этом интегралы по области диэлектрика следует определять численно. Это же относится и к использованию нелинейного диэлектрика в виде $\epsilon \left(r\right)={\epsilon }_l+\tilde{\epsilon }\left({\left|\mathit{E}\mathbf{\left(}\mathit{r}\mathbf{\right)}\right|}^2\right),$ где ${\epsilon }_l$ ‒ линейная часть ДП, $\tilde{\epsilon }\left(0\right)=0.$ Изначально однородный диэлектрик при введении в поле становится неоднородным и влияет на само поле. Если функция $\tilde{\epsilon }$ возрастающая, то в силу выталкивания поля из диэлектрика и ограниченной мощности возникает насыщение поля. Примером сильно нелинейного диэлектрика на СВЧ может служить слабо ионизированная газовая плазма низкого давления в стеклянной колбе. Большая длина свободного пробега электронов позволяет им приобретать за полупериод поля энергию порядка 15‒20 эВ и более, достаточную для ударной ионизации нейтральных атомов. Это приводит к «поджигу» плазмы СВЧ-полем даже в относительно слабых полях, что может быть применено для магнетронного плазменного напыления из газовой фазы [18]. Предложенная модель может быть использована для моделирования магнетронного напыления на достаточно большую поверхность в прямоугольном СВЧ плазменном реакторе.

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

Таким образом, в работе предложены модель и итерационный алгоритм для решения задачи о возбуждении объема в виде прямоугольного параллелепипеда с диэлектриком через окна ввода энергии, которую можно использовать как в линейном, так и в нелинейном случаях. Алгоритм реализован для моделирования установки СВЧ-нагрева. Учет потерь в стенках возможен как введение дополнительных поверхностных стоков с импедансными условиями Лентовича – Щукина. Более точный учет неидеальных стенок требует использования базиса расширенной области и включения в рассмотрение поверхностного слоя диэлектрика с толщиной 2–3 скин-слоя. Полученная модель требует минимальных вычислительных ресурсов и удобна для численной оптимизации (см. приложение).

Финансирование

Работа выполнена при финансовой поддержке Министерства образования и науки Российской Федерации в рамках выполнения государственного задания (проект № FSRR-2020-0004).

Приложение

Интегралы в (14) и (15) имеют вид

$I_m=2\underset{y_0}{\overset{y_0+a_0}{\int }}\sin \left(k_{m}y\right)\sin \left(\pi \left(y-y_0\right)/a_0\right)dy$

и значение

$\begin{array}{c}{I}_m=I\left(k_m,a_0,y_0\right)=\\ =\frac{\sin \left(\left(k_m-\pi /a_0\right)\left(y_0+a_0\right)+\pi y_0/a_0\right)-\sin \left(\left(k_m-\pi /a_0\right)y_0+\pi y_0/a_0\right)}{k_m-\pi /a_0}-\\ -\frac{\sin \left(\left(k_m+\pi /a_0\right)\left(y_0+a_0\right)-\pi y_0/a_0\right)-\sin \left(\left(k_m+\pi /a_0\right)y_0-\pi y_0/a_0\right)}{k_m+\pi /a_0}.\end{array}$(П.1)

Коэффициенты, учитывающие влияние диэлектрика для разложений (16), определены как

${\alpha }_{kd}=i\omega \frac{{\displaystyle \underset{V_d}{\int }\left(J_x^p{E}_{xk}+J_y^p{E}_{yk}+J_z^p{E}_{zk}\right)dxdydz}}{N_k^e\left({\omega }^2-{\omega }_k^2\right)}=i\omega \frac{I_{xk}+I_{yk}+I_{zk}}{N_k^e\left({\omega }^2-{\omega }_k^2\right)},$ (П.2)

${\beta }_{kd}=i{\omega }_k\frac{{\displaystyle \underset{V_d}{\int }\left(J_x^p{E}_{xk}+J_y^p{E}_{yk}+J_z^p{E}_{zk}\right)dxdydz}}{N_k^e\left({\omega }^2-{\omega }_k^2\right)}=i{\omega }_k\frac{I_{xk}+I_{yk}+I_{zk}}{N_k^e\left({\omega }^2-{\omega }_k^2\right)}.$ (П.3)

Интегралы в (17) имеют вид

$I_{xm{m}^{\text{'}}}=2\underset{\left(b-t\right)/2}{\overset{\left(b-t\right)/2}{\int }}\sin \left(k_{m^{\text{'}}}y\right)\sin \left(k_{m}y\right)dy=I_x\left(k_m,k_{m^{\text{'}}}\right),$ (П.4)

$I_{ym{m}^{\text{'}}}=2\underset{\left(b-t\right)/2}{\overset{\left(b-t\right)/2}{\int }}\cos \left(k_{m^{\text{'}}}y\right)\cos \left(k_{m}y\right)dy=I_y\left(k_m,k_{m^{\text{'}}}\right),$ (П.5)

при этом $I_{zm{m}^{\text{'}}}=I_{xm{m}^{\text{'}}}.$ Здесь

$I_{xm{m}^{\text{'}}}=2\left[\frac{\sin \left(\left(k_m-k_{m^{\text{'}}}\right)t\right)\cos \left(\left(k_m-k_{m^{\text{'}}}\right)b\right)}{\left(k_m-k_{m^{\text{'}}}\right)}-\frac{\sin \left(\left(k_m+k_{m^{\text{'}}}\right)t\right)\cos \left(\left(k_m+k_{m^{\text{'}}}\right)b\right)}{\left(k_m+k_{m^{\text{'}}}\right)}\right],$ (П.6)

$I_{ym{m}^{\text{'}}}=2\left[\frac{\sin \left(\left(k_m-k_{m^{\text{'}}}\right)t\right)\cos \left(\left(k_m-k_{m^{\text{'}}}\right)b\right)}{\left(k_m-k_{m^{\text{'}}}\right)}+\frac{\sin \left(\left(k_m+k_{m^{\text{'}}}\right)t\right)\cos \left(\left(k_m+k_{m^{\text{'}}}\right)b\right)}{\left(k_m+k_{m^{\text{'}}}\right)}\right],$ (П.7).

About the authors

Mikhail V. Davidovich

Saratov State University; LLC «NPF “Etna plus”»

Email: davidovichmv@info.sgu.ru

Doctor of Physical and Mathematical Sciences (2000), professor of the Department of Radio Engineering and Electrodynamics (since 2000)

Russian Federation, 83, Astrakhanskaya Street, Saratov, 410012; 11, Vishnevaya Street, Saratov, 410040

Alexander K. Kobetz

Saratov State University; LLC «NPF “Etna plus”»

Email: kobetzak@info.sgu.ru

quality director

Russian Federation, 83, Astrakhanskaya Street, Saratov, 410012, Russia; 11, Vishnevaya Street, Saratov, 410040, Russia

Kirill A. Sayapin

Saratov State University

Author for correspondence.
Email: sayapin_kirill@mail.ru

postgraduate student of the Department of Radio Engineering and Electrodynamics

Russian Federation, 83, Astrakhanskaya Street, Saratov, 410012

References

Supplementary files