Перевод статей(Сницерук Павел). Расчет дифракции шести пучков в слоистых средах с использованием полиномов основных миноров
![]()
|
Расчет дифракции шести пучков в слоистых средах с использованием полиномов основных миноров Юрий Н. Беляев Кафедра математического моделирования и кибернетики Сыктывкарский государственный университет, Сыктывкар 167001 Октябрьский пр. 55, Russia ybelyayev@mail.ru Получено 5 января 2017 г. Принято 5 декабря 2017 г. Опубликовано 10 мая 2018 г. Рассмотрено рассеяние плоских упругих волн на анизотропной слоистой среде в случае шестилучевой дифракции. Разработан матричный метод решения волновых уравнений. Определены коэффициенты преобразования Cij для трех типов падающих волн (горизонтально-поляризованная поперечная волна, вертикально-поляризованная поперечная волна и продольная волна). Найдены представления коэффициентов Cij через элементы матрицы переноса. Представлен метод для расчета коэффициентов. Это не требует решения алгебраической задачи о собственных значениях для волн в анизотропном слое. Некоторые особенности функциональных зависимостей Cij по углам падения, частоте волн и толщине слоя демонстрируются на нескольких примерах кристаллов в трехслойной модели. Очевидно, что преобразования SH волны в SV волны и SV волны в SH волны эквивалентны. Ключевые слова: упругие волны; матричная экспонента; многочлены; рекуррентные соотношения; ошибка усечения. Введение Когда продольная (P-тип) или сдвиговая упругая волна падает на твердотельный слой, некоторая часть энергии отражается и передается слоем в форме волн того же типа, что и падающая. В дополнение к этим волнам изотропный слой может преобразовывать P-волну в вертикальную волну сдвига (SV-типа) и наоборот. Анизотропия упругих свойств в слое приводит к возможности взаимодействия (преобразования) упомянутых волн с горизонтально горизонтальной волной ( SH-типа). В результате этих преобразований каждый упомянутый тип падающей плоской волны может создавать шесть упругих волн в анизотропном слое. В какой пропорции энергия падающей волны распределяется между рассеянными волнами? Ответ на этот вопрос применительно к многослойным средам включает в себя решения волновых уравнений для каждого из слоев и использование граничных условий для полученных полевых переменных ψ, j = 1, ..., 6 (компоненты вектора смещения и тензора напряжений ). Это может быть эффективно достигнуто с использованием метода матрицы переноса. Согласно этому методу, так называемая передаточная матрица связывает значения вектор-функции Ψ = Ψ (x3) при толщинах di и di + 1: Ψ (di + 1) = Ti + 1Ψ (di). Таким образом, Ψ (d) = TΨ (0), где T = i1i = ∏NTi является глобальной передаточная матрица и d = d1 + ··· + dN - толщина пластины N-слоя. Формирование матрицы переноса из решения волнового уравнения второго порядка для однородного слоя является относительно простой процедурой. Вот почему такой матричный подход был лучше всего развит в приложениях двухволновой дифракции света, двухмикровых волн и сдвиговых волн в слоистых средах 4,5, а также в исследованиях электронных состояний полупроводниковых сверхрешеток6,7. Аналитические решения были получены и для некоторых четвертого порядка. матрицы переноса (например, для волн P − SV в изотропном слое, 8 для электромагнитных волн в магнитооптическом слоистом мультимедиа9). В общем случае матрица переноса порядка > 4 определяется и вычисляется как матрица экспонент T = exp (Wd) где матрица n-го порядка состоит из коэффициентов соответствующих дифференциальных уравнений. В обычном подходе матрица передачи 4-го порядка для однородного слоя представляется в виде функции от числа волновых чисел kj (собственных значений матрицы W): T = T (exp (k1d), ..., exp (knd)). Такой подход предполагает решение уравнения дисперсии для каждого слоя рассматриваемой структуры. Волновые числа определяются путем решения алгебраической задачи о собственных значениях. Эти аспекты рассматриваются во многих работах, ссылки на которые можно найти в обзоре10 и в монографии11. Результат расчета матрицы переноса зависит, прежде всего, от точности определения волновых чисел и, во-вторых, от достоверности расчетов экспоненциальное выражение (kjd). Влияние ошибок в решениях этих двух задач на точность расчета матрицы переноса возрастает с ростом параметров kjd. Таким образом, расчет матрицы переноса для неоднородных (или затухающих) волн при больших значениях всей толщины среды или частоты падающей волны может привести к неудовлетворительным результатам. Указанные трудности при расчете экспоненциальной матрицы exp (Wd) - для интерполяционной формулы Лагранжа – Сильвестра и других методов, основанных на использовании собственных значений матрицы W. Обсуждение специфики применения этих методов и сравнение различных подходов для вычисления экспоненты матрицы было проведено в обзоре. Показано, что наиболее конкурентным подходом является масштабирование метода и возведение в квадрат в сочетании с аппроксимацией Pade. Для вычисления матрицы переноса в настоящем документе используется более высокоэластичная дифракция волн в анизотропном слое. Метод расчетов. Основные уравнения. Напряженно-деформированное состояние упругих твердых сред описывается уравнениями движения и законом Гука, которые в декартовой системе координат могут быть представлены соответственно как ρ∂2ttu1 = ∂1P1 + ∂2P6 + ∂3P5, ρ∂2ttu2 = ∂1P6 + ∂2P2 + ∂3P4, (1) ρ∂2ttu3 = ∂1P5 + ∂2P4 + ∂3P5, и ∂1u1 = S11P1 + S12P2 + S13P3 + S14P4 + S15P5 + S16P6, (2a) ∂1u1 = S122u + S22P2 + S23P3 + S24P4 + S25P5 + S26P6, (2b) ∂3u3 = S13P1 + S23P2 + S33P3 + S34P4 + S35P5 + S36P6, (2c) ∂3u2 + ∂2u3 = S14P1 + S24P2 + S34P3 + S44P4 + S45P5 + S46P6, (2d) u1u3 + ∂3u1 = S15P1 + S25P2 + S35P3 + S45P4 + S55P5 + S56P6, (2e) u2u1 + ∂1u2 = S16P1 + S26P2 + S36P3 + S46P4 + S56P5 + S66P6 ,(2f) Где p плотность по времени;t это время; x1, x2, x3 являются декартовыми координатами; ∂j и ∂t -дифференцирование по времени и дифференцирование по пространственной координате xj; u1, u2, u3 являются компонентами вектора смещения; Pl и Sgl соответственно компоненты тензора напряжений и тензора соответствия в матричной форме. что плотность и упругие параметры среды могут зависеть только от одной координаты x3 вдоль оси, перпендикулярной поверхности анизотропного слоя, ρ = ρ (x3), Sgl = Sgl (x3). (3) Деформации в анизотропном слое 0 x3≤d генерируется плоскость исходной волны u0 = A0exp [i(k0 · r-ωt)], (4) где u0, A0, k0 и ω вектор смещения, амплитуды вектора, волновой вектор и угловая частота падающей волны, соответственно, и i это мнимая единица. Матрица переноса упругих волн. Нетрудно показать, что решения задач (1) - (4) имеют вид ‖u2P4u1u3P5P3P1P2P6‖ = ‖ψ1 (x3) ψ2 (x3) ... ψ9 (x3) ‖exp [i ( k01x1 + k02x2 − ωt)], (5) где ψj (x3), j = 1, ..., 9, выражают неизвестные x3-зависимости функций в левой части равенства (5), а k01 = k0sinθ0cosα, k02 = k0sinθ0sinα - проекции волнового вектора на оси x1 и x2 (см. рис. 1 (б). Отметим, что такой выбор последовательностей неизвестной функции - как в формуле (5) позволяет легко идентифицировать случаи точных аналитических решений для матрицы переноса. Возможные направления для распространения волн в слоистой среде определяется соотношениями kjsinθj = k0sinθ0, j = 1, ..., 6, (6) где θj - углы между волновыми векторами kj и осью x3. Подстановка выражений (5) в уравнения. (2a), (2b), (2f) и решение этих уравнений относительно функций ψ7, ψ8, ψ9 дают следующие равенства: ψi = ai1ψ1 + ai2ψ2 + ai3ψ3 + ai5ψ5 + ai6ψ6, i = 7,8,9, (7) где a71 = ik01γ3-ik02γ2, a72 = S24γ2-S14γ1-S46γ3, a73 = ik01γ1 + ik02γ3, a75 = S25γ2-S15γ1-S56γ3, a76 = S23γ2-S13γ1-S36γ3, a81 = ik02γ4-ik01γ5, a82 = S14γ2-S24γ4 + S46γ5 , a83 = -ik01γ2-ik02γ5, a85 = S15γ2-S25γ4 + S56γ5, a86 = S13γ2-S23γ4 + S36γ5, a91 = ik01γ6-ik02γ5, a92 = S24γ5-S14γ3-S46γ6, a93 = ik01γ3 + ik02γ6, a95 = S25γ5-S15γ3- S56γ6, a96 = S23γ5-S13γ3-S36γ6. Используем обозначения: γj = ∆j / (S11∆1-S12∆2 + S16∆3), j = 1, ..., 6; ∆1 = S22S66-S226, ∆2 = S12S66-S16S26, ∆3 = S12S26-S16S22, ∆4 = S11S66-S216, ∆5 = S11S26-S16S12, ∆6 = S11S22-S212. Если принять во внимание зависимости ( 5) и заменить функции ψ7, ψ8, ψ9 в уравнениях. (1), (2c), (2d) и (2e) с помощью (7), то в результате элементарных преобразований получим дифференциальное уравнение для вектора-столбца Ψ (x3) = ‖ψj (x3) ‖ (j = 1, ..., 6), ∂3Ψ = WΨ. (8) Здесь ненулевые элементы матрицы W≡‖wgh‖ имеют значения w11 = w22 = −i (k01a92 + k02a82), w12 = S14a72 + S24a82 + S46a92 + S44, w13 = w52 = -i (k01a72 + k02a92), w14 = w62 = -ik02, w15 = w32 = S14a75 + S24a85 + S46a95 + S45, w16 = w42 = S14a76 + S24a86 + S46a96 + S34, w21 = k201γ6 + k202γ4 -2k01k02γ5-ρω2, w23 = W51 = k201γ3 + k01k02 (γ6-γ1) -k202γ5, w25 = w31 = -i (k01a95 + k02a85), w26 = w41 = -i (k01a96 + k02a86), w55 = w33 = -i (k01a75 + k02a95), w34 = w65 = -ik01, w35 = S15a75 + S25a85 + S56a95 + S55, w36 = w45 = S15a76 + S25a86 + S56a96 + S35, w43 = w56 = -i (k01a76 + k02a96), w46 = S13a76 + S23a86 + S36a96 + S33, w53 = k201γ1 + 2k01k02γ3 + k202γ6 − ρω2, w64 = −ρω2. . Интегрирование системы (8) методом решения Пеано – Бейкера17 дает в форме Ψ (d2) = TΨ (d1), d2 > d1, (9) где T = I + ∫d2d1Wdx3 + ∫d2d1W (x3) ∫x3d1W (ξ1) ∫d1 d3 d1dx2x (d1) dx1x2x1x1x2x1x1x2x1x2x1x2x1x2x1x2x1x1x2x2x1x1x1x2x1x1x1x2x1x1x2x1x1x2x1x1x1xxxxx) (ξ1) ∫ξ1d1W (ξ2) dξ2dξ1dx3 + ... (10) называется передаточной матрицей и является единичной матрицей. Очевидно, что если слоистая структура состоит из набора слоев 0≤x3 Передача матрицы в равномерный слой Передача матрицы в однородный слой 0≤x3≤d является результатом интегрирования в (10) из 0tod при условии, что W = ‖wgh‖ = const. T = ∞∑j = 0Wjdjj! ≡exp (Wd). ( 11) Некоторые случаи, для которых сумма в формуле (11) может быть вычислена аналитически, рассматриваются в статье16. В общем случае анизотропной среды матрица (11) может быть найдена численно. Для этого мы используем представление с помощью Bg (n) -полиномов. T = exp (Wd) = 5∑h = 0 Th (Wd) h, (12) где Th = 1h! + ∞∑j = 61j! H∑g = 0p6-Н + GBJ-1-г (6), ч = 0,1, ..., 5, J = 0,1, ..., (13) pm(m = 1, ..., 6 ) являются коэффициентами характеристического уравнения матрицы Wd, а полиномы Bg (6) определяются рекуррентными соотношениями B0 (6) = ... = B4 (6) = 0; B5 (6) = 1; Bg (6) = 6 ∑h = 1phBg − h (6), g≥ 6. (14) Следует отметить, что при расчете по формулам (12) - (14) не используются собственные значения матрицы Wd, в отличие от методов Бейкера и Лагранжа – Сильвестра, который был применен в методе матрицы пропагатора. Мы используем здесь хорошо известные представления коэффициентов pm через основные миноры матрицы Wd= ‖wghd‖. p1 = 6∑g = 1wggd, p2 = −∑j> i∣∣∣∣∣wiidwijdwjidwjjd∣∣∣∣∣, ..., p6 = −det (Wd). (15) Поэтому в настоящей статье функции Bg (6), определяемые формулой (14), называются многочленами главных миноров. Расчет коэффициентов pm удобно проводить по рекуррентным формулам (см. Метод Le Verrier), mpm = sm − p1sm − 1− ··· −pm − 1s1, m = 1,2, ..., 6, где sm = tr (Wd), следы матриц (Wd) m, m = 1, .. ., 6. Для численных расчетов матрицы переноса по формуле (12) необходимо ввести в формулах (13) аппроксимацию рядов ∞j = 6 для сумм конечного числа термов, скажем J. Th≈1h! + J − 1∑j = 61j! H∑g = 0p6 − h + gBj − 1 − g (6), J = 6 + N, N = 0,1, .... (16) Относительная ошибка усечения. Очевидно, что значения (16) зависят от J: Th = Th (J). Равенство в формуле (16) становится точным, когда J → ∞, т.е. Th (∞), h = 0, ..., 5 соответствуют точным значениям координат матрицы T = exp (Wd) в базисе I, Wd, ..., (Wd) . Давайте оценим значение относительной ошибки при вычислении координаты Thh = ∣∣∣∣Th (∞) −Th (J) Th (∞) ∣∣∣∣. ( 17) Очевидно, что ch = h!/ | 1 + G | ∣∣∣∣∣∣∞∑j = J1j! Cjh∣∣∣∣∣∣ |