Радиотехника и электроника, 2022, T. 67, № 9, стр. 908-915

ВВЕДЕНИЕ

Задачи о распространении волн различной природы через слоистые и неоднородные среды издавна привлекают внимание исследователей [1–6]. В случае плавной неоднородности традиционным приемом расчета является ее приближение в виде ступенчатой слоистой структуры, каждый слой в которой является однородным, а параметры среды от слоя к слою изменяются. При этом для достижения желаемой точности приходится число слоев выбирать достаточно большим, хотя бы в пределах нескольких десятков. Типичным способом решения таких задач является прямой метод, состоящий в нахождении волновой функции в каждом слое по отдельности с последующим их сшиванием на границах между слоями. С увеличением числа слоев громоздкость решения резко возрастает, поэтому прямой метод при числе слоев более пяти для практики становится непригодным [7]. В этом случае можно использовать методы последовательного перехода от слоя к слою, из которых наиболее распространены метод матрицы [8–11] и метод импеданса [1, 2, 12]. В цитированной литературе относительно обоих этих методов, как правило, даются только общие рекомендации или производится расчет для небольшого числа слоев, не более двух–пяти, после чего расчет усложняется. При этом рассматривается падение на структуру единственной волны одного направления.

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

1. ПОСТАНОВКА ЗАДАЧИ

Будем рассматривать волноведущую структуру, содержащую N плоскопараллельных слоев, разделенных $N - 1$ границами. Ось координат вдоль структуры, перпендикулярную слоям и их границам, обозначим через Ox.

Решаем задачу в приближении одномерных волн типа волны в струне [7, 11]. Координатная часть волновой функции, являющейся решением волнового уравнения второго порядка, имеет вид

где ${{k}_{n}}$ – волновое число волны в среде (слое) с номером $n$, ${{A}_{n}}$ и ${{B}_{n}}$ – соответствующие произвольные постоянные, которые должны быть определены из граничных условий. Первое слагаемое в этом выражении описывает волну, распространяющуюся в положительном направлении оси $Ox$, второе – в отрицательном.

Общая геометрия задачи представлена на рис. 1. Здесь в верхнем ряду показана волноведущая структура, слои которой представляют собой среды с волновыми числами ${{k}_{{1,2,...,N}}}$. Координата начала структуры – ${{x}_{0}}$, конца – ${{x}_{N}}$. Координаты границ между средами – от ${{x}_{1}}$ до ${{x}_{{N - 1}}}$. Длины участков между границами, т.е. толщины слоев, обозначим ${{d}_{n}}$, где ${{d}_{n}} = {{x}_{n}} - {{x}_{{n - 1}}}$, причем $n$ меняется от 1 до $N$.

На структуру в положительном направлении оси $Ox$ на слой 1 падает волна $A$, выходит из этого слоя в отрицательном направлении волна B. На ту же структуру со стороны слоя $N$, т.е. в отрицательном направлении оси $Ox$, падает волна $D$, выходит их этого слоя в положительном направлении волна $C$.

Следуя [11], сопоставим структуре эквивалентную цепочку из четырехполюсников, показанную в нижнем ряду рис. 1. То есть каждой среде или каждой границе с номером n соответствует четырехполюсник с тем же номером, где номера n для сред меняются от 1 до $N$, а для границ – от 1 до $N - 1$.

2. СТРУКТУРА ЦЕПОЧКИ ИЗ ЧЕТЫРЕХПОЛЮСНИКОВ

На рис. 1 четырехполюсники, соответствующие средам, обозначены буквами ${{C}_{n}}$. Эти четырехполюсники описываются матрицами распространения, имеющими вид

Четырехполюсники, соответствующие границам, обозначены буквами ${{A}_{n}}$. Эти четырехполюсники описываются матрицами связи, имеющими вид

Полная цепочка четырехполюсников описывается результирующей матрицей ${{\hat {G}}_{N}}$, представляющей собой последовательное соединение чередующихся матриц распространения и связи:

В записи через элементы результирующая матрица имеет вид

Матрица ${{\hat {G}}^{{\left( N \right)}}}$ связывает волны A и B, входящие в начале структуры, с волнами C и D, выходящими в ее конце, в соответствии с соотношением

В алгебраической записи эта связь имеет вид

Для реализации такой связи в полном виде надо найти элементы матрицы $g_{{i\,k}}^{{\left( N \right)}}$ в выражении через параметры структуры, такие как волновые числа ${{k}_{n}}$ и протяженности участков структуры ${{d}_{n}}$.

3. ПОСЛЕДОВАТЕЛЬНОЕ НАХОЖДЕНИЕ ЭЛЕМЕНТОВ РЕЗУЛЬТИРУЮЩЕЙ МАТРИЦЫ

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

Так, для структуры, состоящей из одного участка среды без границ, имеем только матрицу распространения

то есть

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

Для структуры, состоящей из двух сред, разделенных одной границей, результирующая матрица содержит одну матрицу распространения, затем одну матрицу связи и еще одну матрицу распространения:

Элементы такой результирующей матрицы вычисляем по формулам

Действуя далее таким же порядком, на основе элементов матрицы ${{\hat {G}}^{{\left( n \right)}}}$ получаем элементы матрицы ${{G}^{{\left( {n + 1} \right)}}}$:

И, наконец, на основе элементов матрицы ${{G}^{{\left( {N - 1} \right)}}}$, находим элементы матрицы ${{G}^{{\left( N \right)}}}$:

Таким образом, матрицу ${{\hat {G}}^{{\left( N \right)}}}$ можно считать найденной. Эта матрица связывает волны A и B, входящие в структуру, с волнами $C$ и $D$, выходящими из нее в соответствии с соотношением (6).

4. ПРЕОБРАЗОВАНИЕ МАТРИЦЫ ДЛЯ ЧЕТЫРЕХПОЛЮСНИКОВ В МАТРИЦУ ДЛЯ ВОЛНОВЕДУЩЕЙ СТРУКТУРЫ

Теперь перейдем от схемы четырехполюсников (см. рис. 1, нижняя часть) к схеме многослойной волноведущей структуры (см. рис. 1, верхняя часть). Здесь входящими в структуру являются волны A и D, а выходящими – волны B и C. Чтобы найти выходящие волны через входящие надо систему (7), (8) разрешить относительно B и C, т.е. преобразовать ее к виду

В матричной записи эта система имеет вид

где матрица ${{\hat {H}}^{{\left( N \right)}}}$ имеет вид

При этом выходящие из структуры волны $B$ и $C$ при входящих волнах $A$ и $D$ определяются соотношениями

где элементы $g_{{ik}}^{{\left( N \right)}}$ определяются формулами (23)–(26).

5. КОЭФФИЦИЕНТЫ ОТРАЖЕНИЯ И ПРОХОЖДЕНИЯ

Рассмотрим сначала случай падения на структуру только одной волны с амплитудой $A$, распространяющейся в положительном направлении оси $Ox$. От структуры отражается волна $B$ и проходит волна $C$. Для расчета коэффициентов отражения и прохождения по энергии воспользуемся следующими формулами [7, 11]:

Так, при $D\,\, = \,\,0$ в соответствии с (31), (32) получаем коэффициент отражения

и прохождения

Из структуры формул (35) и (36) с учетом (23)–(26) видим, что при $N\,\, = \,\,2$ и действительных ${{k}_{{1,\,2}}}$ формулы (35), (36) принимают вид

что соответствует классическому виду коэффициентов отражения и прохождения по энергиям для случая двух сред, разделенных одной границей [7, 11].

В случае падения одной волны встречного направления, т.е. при A = 0, аналогичным образом в соответствии с (31), (32) получаем

Для двух сред в этом случае, подобно (37), (38), получаем

6. АЛГОРИТМ НАХОЖДЕНИЯ КОЭФФИЦИЕНТОВ ОТРАЖЕНИЯ И ПРОХОЖДЕНИЯ

Из приведенной выше последовательности действий видно, что вычисления по формулам (19)–(22) выполняются многократно циклическим образом с увеличением верхнего индекса элементов на единицу. Для удобства машинного программирования представим полный расчет амплитуд выходящих волн $B$ и $C$ в виде алгоритма с последовательными шагами.

1. В качестве подготовительного шага следует задать начальные параметры структуры – волновые числа ${{k}_{n}}$ и толщины участков ${{d}_{n}}$ для всех $n$ от 1 до $N$, а также амплитуды входящих волн $A$ и $D$.

2. Находим элементы $g_{{i\,k}}^{{\left( 1 \right)}}$ в соответствии с формулами (10)–(13).

3. Используя элементы $g_{{i\,k}}^{{\left( 1 \right)}}$, находим элементы $g_{{i\,k}}^{{\left( 2 \right)}}$ с помощью формул (15)–(18).

4. Следующие действия выполняются в цикле из таких же последовательных шагов с увеличением текущего индекса элемента на единицу по формулам (19)–(22). Цикл продолжается до получения элементов $g_{{i\,k}}^{{\left( N \right)}}$, соответствующих формулам (23)–(26). (При машинном счете эти действия удобно выполнять с помощью единой подпрограммы.)

5. Теперь, когда элементы $g_{{i\,k}}^{{\left( N \right)}}$ найдены, их следует использовать для получения амплитуд выходящих волн $B$ и $C$ в соответствии с формулами (31)–(32). На этом шаге заканчивается получение амплитуд выходящих из структуры волн в случае падения двух волн встречных направлений.

6. В качестве заключительного шага может быть определение коэффициентов отражения и прохождения по формулам (35), (36) и (39), (40).

7. ПРИМЕНЕНИЕ АЛГОРИТМА ДЛЯ РЕШЕНИЯ ЗАДАЧИ О СТРУКТУРЕ СО СТУПЕНЧАТОЙ НЕОДНОРОДНОСТЬЮ

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

На рис. 2 представлено распределение волнового числа по ступенькам при различной степени крутизны барьера. Рассмотрим распространение одномерной волны через барьерную ступенчатую структуру, разбитую вдоль оси $Ox$ на 16 однородных участков (в общем случае это число может быть произвольным). Полная длина структуры $L$ равна 16 условным единицам, начало структуры соответствует координате х₀, конец – координате ${{x}_{{16}}}$. Координаты концов участков обозначены через ${{x}_{n}}$, где $n$ – номер участка. Длины участков ${{d}_{n}}$ определяют по формуле

Волновое число$,$ будучи действительным, от участка к участку линейно возрастает по закону

где $\Delta \,k$ – величина постоянная, отсчет начинается с $n\,\, = \,\,2$ и продолжается до $n = 16$.

Длины первого и последнего участков равны друг другу, т.е. ${{d}_{1}}$ = ${{d}_{{16}}}$. Длины всех остальных участков равны между собой, но от ${{d}_{{1,\,16}}}$ могут отличаться. При этом часть структуры от участка 2 до участка 15, т.е. между координатами ${{x}_{1}}$ и ${{x}_{{15}}}$ разделена на 14 равных по длине участков длиной $\Delta \,d$:

Часть структуры, где волновое число увеличивается, т.е. длина барьера от участка 2 до участка 15, обозначена через ${{L}_{B}}$:

Из рис. 2 видно, что по мере уменьшения $\Delta \,d$ крутизна барьера возрастает, а его длина ${{L}_{B}}$ соответственно уменьшается.

Ограничимся случаем падения на структуру только одной волны с амплитудой $A$, распространяющейся в положительном направлении оси $Ox$. От структуры отражается волна $B$ и проходит волна $C$. Для расчета коэффициентов отражения и прохождения по энергии воспользуемся формулами (33), (34) при $N = 16$:

Из рис. 2 видно, что крутизна барьера при уменьшении его длины возрастает. При максимальной крутизне барьера, т.е. при ${{L}_{B}} = 0$, коэффициенты отражения и прохождения определяются формулами, соответствующими задаче для двух сред, разделенных одной границей [7, 11]:

С целью нормировки в данной работе значение волнового числа на первом участке ${{k}_{1}}$ было принято равным $\Delta k$, так что при изменении $\Delta \,k$ весь интервал значений волнового числа пропорционально изменялся. Как видим, с учетом ${{k}_{1}} = \Delta k$ и ${{k}_{{16}}} = \,\,16\,\Delta \,k\,$ значения $R_{e}^{{\left( 2 \right)}}$ и $T_{e}^{{\left( 2 \right)}}$ остаются неизменными при любой величине $\Delta \,k,$ т.е. $R_{e}^{{\left( 2 \right)}}$ = 0.7785 и $T_{e}^{{\left( 2 \right)}}$ = 0.2215. Сумма этих значений равна единице.

Результаты расчета представлены на рис. 3. Показаны зависимости коэффициентов отражения ${{R}_{e}}$ и прохождения ${{T}_{e}}$ от длины барьера ${{L}_{B}}$ при различных значениях перепада волнового числа $\Delta k$. Значения длины барьера на горизонтальной оси отсчитываются в обратном порядке, т.е. длина барьера слева направо уменьшается, так что его крутизна увеличивается.

Из рисунка видно, что в любом случае при уменьшении ${{L}_{B}}$, т.е. при увеличении крутизны барьера, отражение увеличивается, а прохождение уменьшается. Кривые прохождения 4–6 являются симметричным отражением в уровне 0.5 отн. ед. кривых отражения 1–3, причем сумма тех и других значений равна единице, что соответствует закону сохранения энергии [7, 8, 11].

При обращении барьера в вертикальный перепад, т.е. при ${{L}_{B}} \to 0$, все зависимости для отражения 1–3 приходят в одну точку, соответствующую ${{R}_{e}}$ = 0.7785 отн. ед., а зависимости прохождения 4–6 – в точку ${{T}_{e}}$ = 0.2215. Эти значения соответствуют перепаду волновых чисел в 16 раз на структуре из двух сред, разделенных одной границей [7, 11].

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

При увеличении перепада отражение (кривая 2) при большой длине барьера (${{L}_{B}}$ = 14 отн. ед.) значительно меньше, чем при малом перепаде и составляет 0.2946. По мере уменьшения длины барьера отражение возрастает, стремясь к величине 0.7785. Прохождение (кривая 5) ведет себя зеркальным образом, составляя при большой длине барьера 0.7054, а при малой спадает до 0.2215. Обе кривые ведут себя полностью монотонно.

При значительном перепаде (кривые 3, 6) отражение (кривая 3) при большой длине барьера (${{L}_{B}}$ = 14 отн. ед.) еще меньше, чем в предыдущем случае, и составляет 0.0759 отн. ед. По мере уменьшения длины барьера отражение в целом возрастает, но вплоть до величины барьера 2 отн. ед. проходит через несколько (около трех) осцилляций, после чего резко устремляется вверх к величине 0.7785 отн. ед. Прохождение (кривая 6) ведет себя зеркальным образом, начинаясь при 0.9241 и, проходя через те же осцилляции, резко падает до величины 0.2215 отн. ед. При увеличении перепада волнового числа выше 1.0 отн. ед. осцилляции зависимостей для обоих коэффициентов становятся беспорядочными, приобретая близкий к хаотическому характер.

В ходе контрольной проверки было установлено, что наличие осцилляций обусловлено периодическим характером распространяющихся волн, определяемом формулой (1). При этом решающее значение приобретает средняя длина волны, соответствующая интервалу волнового числа в 8 раз больше единичного перепада волнового числа (т.е. 16 : 2). Осцилляции проявляются только в том случае, когда длина волны становится меньше полной длины барьера, что имеет место для кривых 3, 6 и отсутствует для кривых 1, 2 и 4, 5. При этом хаотизация осцилляций начинается тогда, когда длина волны становится меньше длины отдельного единичного участка, т.е. менее 1 отн. ед.

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

При практическом определении коэффициентов отражения и прохождения при заданной длине барьера следует выбирать условия, при которых кривые, типа приведенных на рис. 3, не имеют осцилляций, т.е. при принятых здесь числовых параметрах величина $\Delta k$ должна быть меньше 0.10 отн. ед. Поскольку эта величина определяется длиной барьера, то при значительной его высоте разбиение должно выполняться на число участков, равное или большее отношению высоты барьера к величине $\Delta k$.

ЗАКЛЮЧЕНИЕ

Предложена алгоритмическая реализация метода матрицы для задачи о падении на многослойную структуру двух встречных волн.

В приближении одномерных волн рассмотрена волноведущая структура, содержащая произвольное число однородных слоев, разделенных границами. Для структуры предложена эквивалентная схема цепочки из четырехполюсников, соответствующих элементам структуры – однородным участкам и границам между ними. Каждому четырехполюснику поставлена в соответствие квадратная матрица второго порядка. Получены матрицы двух видов: для участков – матрицы распространения, для границ – матрицы связи. Результирующая матрица представлена в виде последовательно произведения чередующихся матриц того и другого вида. Отмечено, что вся цепочка состоит из последовательно соединенных пар двух матриц, первая из которых является матрицей распространения, а вторая – матрицей связи. На основании анализа произведения двух матриц, составляющих одну пару, выявлено, что элементы матрицы такой пары для любого участка цепочки могут быть выражены через элементы матрицы такой же пары для участка, предшествующего данной паре.

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

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