Full Text

Введение Определение упругих постоянных тонких пленок [1, 2], исследование критических напряжений в пленках [3, 4], разработка новых электроакустических и акустооптических преобразователей [5, 6] и некоторые другие проблемы тесно связаны с расчетом напряженно-деформированного состояния пленки, вызванного упругими волнами. Распределение напряжений в твердом слое зависит от вида возможных упругих волн. Как известно [7, 8, 9], в любом направлении изотропной среды могут распространяться упругие волны трех типов, различающихся направлением колебаний частиц среды. В продольной волне (обозначаемой как волна P-типа) колебания происходят вдоль направления волнового вектора. В волнах сдвига амплитудный вектор колеблется перпендикулярно волновому вектору. Если в сдвиговой волне колебания направлены перпендикулярно или параллельно плоскости падения, то волна называется соответственно горизонтально поляризованной (SH-тип) и вертикально поляризованной (SV-тип) волной сдвига. При падении продольной или сдвиговой волны на слой твердого тела, часть энергии отражается и проходит через слой в виде волн того же типа, что и падающая волна. В дополнение к этим волнам изотропный слой может преобразовать P-волну в вертикальную волну сдвига и наоборот [8, 10]. Анизотропия упругих свойств среды распространения приводит к возможности взаимного преобразования волн всех трех типов. Соответствующие исследования для сред, обладающих орторомбической и тригональной симметриями проводились, например, в работах [11, 12]. В любом направлении анизотропной среды скорости распространения трех волн могут различаться друг от друга, и их поляризации не являются ни чисто продольными, ни чисто поперечными. Волна, у которой направления колебаний наиболее близки к направлению волнового вектора, называется квазипродольной. Две другие волны называются квазипоперечными [5]. В анизотропном слое конечной толщины d в результате воздействия плоской упругой волны образуются шесть волн (по три в направлении отражения и пропускания), суммарное воздействие которых формирует напряженное состояние слоя. В расчетах распространения волн в слоистых средах получил развитие подход, основанный на матричном методе решения системы дифференциальных уравнений [13, 14]. Описание этого метода применительно к упругим волнам имеется в обзорах [11, 15] и монографиях [16-18]. Согласно этому методу значения компонент вектора смещения и тензора напряжений на одной границе слоя выражаются через аналогичные компоненты на противоположной границе с помощью так называемой матрицы переноса где W - матрица n-го порядка, составленная из коэффициентов соответствующей системы дифференциальных уравнений. Относительно просто решается задача рассеяния горизонтальной волны сдвига изотропным слоем [18], которая описывается матрицей переноса второго порядка. Для матрицы переноса четвертого порядка, определяющей рассеяние P-SV-волн изотропным слоем, также известна аналитическая формула [19]. Эти матрицы переноса, как и их аналоги из оптики (см. например [20, 21]), находятся вычислением матричной экспоненты по формуле Лагранжа-Сильвестра [13] с помощью собственных значений матрицы в виде функции В случае матрицы W общего вида, порядок которой решение алгебраической задачи на собственные значения возможно только численными методами. Поэтому результат вычислений матрицы переноса по формуле Лагранжа-Сильвестра отягощен погрешностями вычислений собственных значений и экспонент Похожие проблемы возникают, если для вычислений матричной экспоненты воспользоваться методом Бэкера-Вандермонда [22], интерполяционной формулой Ньютона [23], методом Лапласа [24], канонической формулой Жордана [13]. Указанные погрешности накапливаются при вычислении матрицы переноса всей структуры с ростом числа слоев N и величин параметров распространения где - это толщины слоев. В рассматриваемой задаче имеют смысл волновых чисел, и, следовательно, размерность модуля параметра распространения пропорциональна произведению частоты излучения и толщины. Возрастание именно этих величин может приводить к неудовлетворительным результатам при вычислении матрицы переноса и зависящих от нее величин (см. статью [25] и библиографические указания в ней). Это одна из причин, по которым перечисленные выше именные методы, известные из учебной и справочной литературы, практически не применяются для вычислений матричных экспонент большого порядка. Проблема вычисления матричной экспоненты является центральной при решении больших систем дифференциальных уравнений в различных приложениях. Это стимулирует развитие известных методов расчетов, сравнительный анализ которых имеется в работах [26, 27], и разработку новых подходов, примеры которых рассмотрены в статьях [28-35]. Из набора известных методов наилучшим численным подходом к нахождению матричной экспоненты , по мнению авторов обзора [28], является алгоритм масштабирования и кратного квадрирования (МКК): 1) масштабирование матрицы с целью понижения ее нормы, 2) вычисление экспоненты масштабированной матрицы 3) вычисление с помощью кратного квадрирования. В методе МКК при выполнении второго этапа подразумевается использование аппроксимации Паде или аппроксимации разложения в ряд Тейлора конечной суммой. Примеры реализаций этих способов расчетов представлены в работах [36, 37]. Метод расчета [36] реализован в математическом пакете MATLAB и используется некоторыми исследователями [35] для сравнительного тестирования своих алгоритмов. В данной работе для вычисления матрицы переноса упругих волн применяется метод многочленов главных миноров [38], который, по крайней мере для матриц порядка [39], обеспечивает более надежное и точное вычисление матричной экспоненты в сравнении с методами [36, 37]. 1. Метод вычислений 1.1. Основные уравнения Динамика напряженно-деформированного состояния упругой среды описывается уравнениями движения и законом Гука. Эти уравнения в декартовой системе координат выражаются формулами (1) где - плотность; - время; - декартовы координаты; и - компоненты вектора смещений, тензора напряжений и тензора упругой податливости соответственно. Уравнения (1) решаются применительно к плоскослоистой среде, у которой плотность и упругие параметры среды зависят только от одной координаты вдоль оси, перпендикулярной поверхности анизотропного слоя (рис. 1): (2) Деформации в анизотропном слое порождаются плоской волной: (3) падающей из области . Здесь и являются соответственно вектором смещений, амплитудным вектором, волновым вектором и циклической частотой падающей волны; - мнимая единица. Решение задачи (1)-(3) методом разделения переменных показывает, что все величины и имеют одинаковую зависимость от переменныхи в виде , как в падающей волне. Поэтому если обозначить через неизвестные зависимости компонент вектора смещений и тензора напряжений от координаты то указанные компоненты принимают вид (4) где индекс j принимает значения от 1 до 9 в том же порядке, в каком перечислены компоненты вектора смещений и тензора напряжений в левой части равенства (4). С учетом равенств (4) уравнения (1) редуцируются в систему уравнений (5) (6) где и использованы обозначения: - компоненты тензора упругой податливости в матричной форме [40]; 1.2. Матрица переноса Решение системы уравнений (5) в матричной форме имеет вид (7) где (8) в теории дифференциальных уравнений называется матрицантом, а в теории волн в слоистых средах - матрицей переноса. Здесь и далее I обозначает единичную матрицу шестого порядка. Если структура толщиной d состоит из N однородных слоев, т.е. то где и матричная экспонента определяется формулой (9) Соотношения (7)-(9) позволяют выразить значения функций через значения этих функций на глубине . После этого нахождение значений функций по формулам (6) не представляет труда. 1.3. Вычисление матрицы переноса однородного слоя в данной работе проводится по формуле (10) Здесь являются коэффициентами характеристического уравнения матрицы . Как известно, коэффициенты характеристического уравнения матрицы равны с точностью до знака суммам соответствующих главных миноров этой матрицы, в частности равен следу матрицы, а . Вследствие этого функции определяемые рекуррентными уравнениями [38] (11) называются многочленами главных миноров. Для выполнения расчетов по формулам (10)-(11) требуется предварительное нахождение матриц и коэффициентов . Эти коэффициенты вычисляются рекуррентно по методу Леверье [13]. Для каждой из найденных матриц определяется след и применяются формулы (12) Выражение (10) получено из точной формулы в результате усечения ряда до конечной суммы . Порождаемая таким усечением относительная погрешность в вычислениях элементов матрицы переноса зависит от числа N1 и удовлетворяет неравенству [41] (13) где - минимальная скорость распространения упругих волн в рассматриваемом слое. Параметр масштабирования m является минимальным целым, которое выбирается в соответствии с условием Пример оценки параметра масштабирования для кристаллического слоя кремния представлен в статье [41]. Выбор N1 = 14 обеспечивает вычисление матрицы переноса с помощью алгоритма (10)-(12) с двойной арифметической точностью 10-16, а N1 = 25 соответствует четверной арифметической точности 2,5×10-32. Для реализации этого метода требуется выполнение четырех матричных перемножений. Для сравнения Алгоритм 2.3 [36], использующий для вычисления матричной экспоненты аппроксимацию Паде, требует выполнения шести матричных перемножений плюс решение одного матричного уравнения. При этом нормализованная относительная погрешность превышает 10-16. 1.4. Геометрия рассеяния Из равенств (4) следует, что лучи всех волн, возникающих в слоистой среде, лежат в одной плоскости (рис. 1, а), и проекции волновых векторов на оси и определяются соответственно формулами Если обозначить угол между волновым вектором и осью , то где и углы определяются законом преломления (отражения): (14) Рассмотрим случай, когда области и ограничивающие анизотропную слоистую структуру, являются изотропными средами. а б в г Рис. 1. К определению направляющих косинусов падающей волны: a - углы падения и б - падающая волна SH-типа; в - падающая волна SV-типа; г - падающая волна P-типа Fig. 1. To defining the direction cosines of an incident wave: a - angles of incidence and b - SH-type of the incident wave; c - SV-type of the incident wave; d - P-type of the incident wave Как известно, в любом направлении изотропной среды могут распространяться волны P-, SH- и SV-типа, скорости распространения которых определяются упругими постоянными Ламе: и В результате воздействия падающей волны на анизотропный слой в изотропных областях и могут возникнуть от двух до шести волн (рис. 2): где и - соответственно амплитудные и волновые векторы; индексами 1 и 4 обозначены горизонтальные волны сдвига, индексы 2 и 5 соответствуют вертикальным волнам сдвига, наконец, индексы 3 и 6 отмечают продольные волны. Волновые числа этих волн имеют следующие значения: , . На рисунках направления колебаний амплитудных векторов волн SH-, SV- и P-типа показаны стрелками зеленого, синего и красного цветов соответственно. Рис. 2. К определению направляющих косинусов рассеянных волн: А1 и А4 - амплитудные векторы SH-волн, А2 и А5 - амплитудные векторы SV-волн, А3 и А6 - амплитудные векторы P-волн Fig. 2. To defining the direction cosines of the scattered waves, А1 and А4 are the amplitude vectors of SH waves, А2 and А5 are the amplitude vectors of SV waves, А3 and А6 are the amplitude vectors of P waves 1.5. Граничные условия Упругие свойства изотропных сред (компоненты тензора упругой податливости) выражаются через постоянные Ламе. С помощью последних уравнения закона Гука (1) существенно упрощаются, и напряжения в изотропных областях и выражаются через компоненты вектора смещения формулами (15) где если и если Смещения частиц среды в этих изотропных средах обусловлены волнами, которые в них распространяются. Таким образом, в области (16) в области (17) Подстановка векторов (16) и (17) в уравнения (15) позволяет найти компоненты тензора напряжений на границах анизотропного слоя. Следовательно, используя соотношения (3), (4), (15), (16) и (17), несложно выразить значения неизвестных функций на границах и анизотропного слоя через амплитуды рассеянных волн. Так, для случая, когда такой подход дает следующий результат: (18) где - направляющие косинусы амплитудных векторов Аg волн, В последней формуле индексы указывают три возможных типа падающей волны. Значение соответствует SH-волне (см. рис. 1, б), в которой направляющие косинусы амплитудного вектора соответствует SV-волне (см. рис. 1, в) с направляющими косинусами и соответствует P-волне (см. рис. 1, г) с направляющими косинусами Рис. 2 определяет направляющие косинусы амплитудных векторов рассеянных волн: 1.6. Вычисление коэффициентов преобразований волн В результате подстановки выражений функций из формул (18) в уравнение получается система алгебраических уравнений относительно амплитуд волн Решение этой системы, например, по методу Гаусса дает значения коэффициентов преобразований падающей волны в рассеянные волны. Энергия волны пропорциональна квадрату модуля амплитуды. Поэтому величина характеризует долю энергии падающей волны, передаваемой h-й волне, и называется интенсивностью h-й волны. 1.7. Вычисление акустических напряжений Вычиление акустических напряжений на глубине производится по формуле где значения компонент матрицы определяются из (18) с учетом найденных коэффициентов преобразований . 2. Результаты расчетов Представленный алгоритм расчетов был использован в исследовании акустических напряжений в кристаллических слоях. Некоторые результаты вычислений для слоев молибдата свинца и кремния показаны на рис. 3-6. Все значения амплитуд колебаний компонент тензора напряжений даны в единицах Па. Циклическая частота Гц. Плоскость совпадает с гранью кристалла (001). Параметры кристаллов, использованные в расчетах, представлены в таблице, и области предполагались твердыми с параметрами Параметры кристаллов Parameters of crystals Кристалл ρ S11 S12 S13 S16 S33 S44 S66 PbMoO4 6950 21 -12,4 -4,93 -15,5 16,6 37,5 40,6 Si 2329 7,69 -2,14 -2,14 0 7,69 12,58 12,58 Примечание: плотность кристалла ρ, кг/м3; компоненты тензора Sgh, 10-12Па-1. На рис. 3-5 представлены результаты расчетов напряжений на границах слоя PbMoO4 толщиной , возникающих при падении на этот слой горизонтальной волны сдвига (см. рис. 3), вертикальной волны сдвига (см. рис. 4) и продольной волны (см. рис. 5). Зависимости отмечены на рисунках цифрами 1, 2, 3, соответствующими значению индекса i у компонент тензора напряжений. На рис. 3-5, a показаны зависимости энергетических характеристик волн, рассеянных в область , а на рис. 3-5, б - соответствующие характеристики волн, рассеянных в область . На этих рисунках кривые отмечены цифрами 1, 2, 3, 4, 5, 6, соответствующими значениям индекса j. Кроме этого, функциональные зависимости и интенсивностей горизонтальных волн сдвига изображены зеленым цветом, интенсивности вертикальных волн сдвига и нарисованы синим цветом, а для показа интенсивностей и продольных волн использован красный цвет. Такие же цвета используются в рис. 6. В результате перерассеянии волн на границах кристаллического слоя рассеянные волны могут усиливаться или, наоборот, ослабляться за счет интерференции. Число таких экстремумов возрастает с увеличением толщины и уменьшением скорости распространения волн. Наряду с системой интерференционных максимумов и минимумов на рис. 3 и 4 видны экстремумы резонансного характера. Причиной их возникновения является следующее. Рис. 3. Амплитуды напряжений, возникающих на границах слоя PbMoO4 под воздействием горизонтальной волны сдвига: , м Fig. 3. Amplitudes of stresses arising at the boundaries of PbMoO4 layer under the influence of a horizontal shear wave: , m Рис. 4. Амплитуды напряжений, возникающих на границах слоя PbMoO4 под воздействием вертикальной волны сдвига: , м Fig. 4. Amplitudes of stresses arising at the boundaries of PbMoO4 layer under the action of a vertical shear wave: , m Рис. 5. Амплитуды напряжений, возникающих на границах слоя PbMoO4 под воздействием продольной волны: , м Fig. 5. Amplitudes of stresses arising at the boundaries of PbMoO4 under the action of a longitudinal wave: , m Если сдвиговая волна падает на кристалл под углом (критический угол определяется в соответствии с законом (14) равенством ), то продольные волны, генерируемые кристаллом, становятся неоднородными, распространяющимися вдоль оси . Для рассматриваемой структуры . При определенной толщине кристаллического слоя и углах падения суммарное смещение вертикальных волн сдвига внутри кристалла создает поперек слоя модулированную структуру. Дифракция продольных волн на этой структуре порождает резко выраженные экстремумы в спектрах рассеянных волн при некоторых углах падения . Именно этим объясняются резко выраженные экстремумы напряжений в спектрах рассеянных волн. Корреляции между экстремумами кривых , с одной стороны, и экстремумами функций и - с другой, определяются формулами (18) и хорошо видны при сравнении рис. 3, a и рис. 3, в, рис. 3, б и рис. 3, г, рис. 4, a и рис. 4, в, рис. 4, б и рис. 4, г. На распределение энергии падающей волны между волнами, рассеянными кристаллом, большое влияние оказывает также вращение плоскости поляризации сдвиговой волны по мере прохождения кристаллической среды. Одним из проявлений этого эффекта являются равенства и , которые выполнялись с высокой точностью во всех численных экспериментах, проведенных в рамках данной работы и [41]. Иллюстрациями сказанного являются синяя кривая на рис. 3, a и зеленая кривая на рис. 4, a, синяя кривая на рис. 3, б и зеленая кривая на рис. 4, б. В результате вращения плоскости поляризации сдвиговой волны падающая волна SH-типа (SV-типа) при отражении или пропускании полностью может превращаться в волну SV-типа (SH-типа). Пример такого преобразования горизонтальной волны сдвига кристаллическим слоем кремния показан на рис. 6. При углах падения и коэффициент , т.е. отраженная от кристалла сдвиговая волна является вертикально поляризованной. Рис. 6. Интенсивности волн, излучаемых слоем Si под воздействием горизонтальной волны сдвига , м Fig. 6. The intensities of waves emitted by a Si layer under the influence of a horizontal shear wave , m Амплитуды всех других волн резонансно увеличены. Продольные волны являются неоднородными, их волновые векторы направлены вдоль поверхностей кристалла и . Эти волны энергию от кристаллического слоя не переносят. При указанных условиях на границах кристаллического слоя кремния возникают резонансы напряжений. Амплитуды колебаний этих напряжений имеют следующие значения: , , , , , . Ширина резонансов спектров рассеяния, показанных на рис. 6, и соответствующих им резонансов напряжений равна примерно пяти угловым секундам. Заключение Такие эффекты, как полное преобразование поляризации сдвиговой волны, дифракционные резонансы, которые возможны при шестилучевой дифракции упругих волн в анизотропном слое, происходят в очень узком диапазоне углов падения . Положение такого диапазона может быть изменено за счет небольших изменений значений второго угла падения и параметра распространения . Аналитический расчет таких эффектов для среды с тензором упругой податливости общего вида, по-видимому, невозможен. Результаты численного моделирования спектров рассеяния упругих волн и создаваемых ими напряжений зависят от точности расчетов матрицы переноса. Метод многочленов главных миноров позволяет находить эту матрицу с заданной точностью. Важная особенность метода расчета, представленного в данной работе, состоит в том, что при исследовании распространения волн в многослойной среде требуется знание волновых чисел только для двух крайних, ограничивающих структуру слоев (на рис. 2 это слои и ). Эти волновые числа учитываются в граничных условиях, но не влияют на результат вычислений матриц переноса внутренних слоев .

About the authors

Yu N Belyayev

Syktyvkar State University

References

Statistics

Views

Abstract - 149

PDF (Russian) - 57