CC BY
17
Беляев Ю.Н. Метод расчета акустических напряжений при шестилучевой дифракции в слоистых средах // Вестник Пермского национального исследовательского политехнического университета. Механика. - 2018. - № 4. - С. 82-92. Б01: 10.155 93/регш.шесЬ/2018.4.07
Belyayev Yu.N. Method for calculating acoustic stresses in six-beam diffraction in layered media. PNRPU Mechanics Bulletin, 2018, no. 4, pp. 82-92. DOI: 10.15593/perm.mech/2018.4.07
ВЕСТНИК ПНИПУ. МЕХАНИКА
№ 4, 2018 PNRPU MECHANICS BULLETIN
http://vestnik.pstu.ru/mechanics/about/inf/
Б01: 10.15593/регш.шесЬ/2018.4.07 УДК 534.2+512.64+519.622
МЕТОД РАСЧЕТА АКУСТИЧЕСКИХ НАПРЯЖЕНИЙ ПРИ ШЕСТИЛУЧЕВОЙ ДИФРАКЦИИ В СЛОИСТЫХ СРЕДАХ
Ю.Н. Беляев
Сыктывкарский государственный университет, Сыктывкар, Россия
О СТАТЬЕ
АННОТАЦИЯ
Получена: 03 июля 2018 г. Принята: 09 сентября 2018 г. Опубликована: 28 декабря 2018 г.
Ключевые слова:
упругие волны, матричная экспонента, метод масштабирования и квадрирования, многочлены главных миноров, рекуррентные соотношения, относительная погрешность усечения, резонансы напряжений.
Теоретически исследуются напряжения, возникающие в слоистой среде в результате воздействия акустической волны. В общем случае под действием падающей упругой волны в анизотропном слое формируются шесть волн, три из которых направлены в область отражения и три - в область прохождения. Напряженно-деформированное состояние слоя является результатом суммарного воздействия этих волн и описывается уравнениями движений сплошной среды и обобщенным законом Гука. Эта система дифференциальных уравнений решается относительно компонент вектора смещения и тензора напряжений в декартовой системе координат в матричной форме. Компоненты вектора смещений и тензора напряжений на двух противоположных границах слоя толщиной С выражаются друг через друга с помощью матрицы переноса шестого порядка Т = ехр(М СС). Вычисление этой экспоненты проводится с помощью многочленов главных миноров матрицы М и не требует нахождения собственных значений матрицы М . Этот метод обеспечивает более точное и надежное вычисление матрицы переноса М-слойной среды Т = ТлТлм-.-Т в сравнении с другими известными алгоритмами.
Амплитуды волн, рассеянных анизотропным слоем, выражены через элементы матрицы переноса. Распределение акустических напряжений по толщине анизотропного слоя определяется амплитудами рассеянных волн и элементами соответствующих матриц переноса. Этот способ расчета акустических напряжений продемонстрирован для падающих волн SH-, вУ- и Р-типа на трехслойной модели: изотропный слой - кристаллический слой -изотропный слой.
Приведено сравнение спектров рассеяния упругих волн и зависимостей напряжений от углов рассеяния для кристаллических слоев кремния и молибдата свинца. Дана интерпретация резонансов акустических напряжений, возникающих в кристаллическом слое под действием сдвиговых волн.
©ПНИПУ
© Беляев Юрий Николаевич - к.ф.-м.н., доц., зав. каф., e-mail: [email protected], 0000-0002-9182-8721
Yuriy N. Belyayev - CSc in Physical and Mathematical Sciences, Associate Professor, Head of Department, e-mail: [email protected], 0000-0002-9182-8721
-„ „—I Эта статья доступна в соответствии с условиями лицензии Creative Commons Attribution-NonCommercial 4.0 International
(j)^ I License (CC BY-NC 4.0) к-^^Н^НЗ^Н This work is licensed under a Creative Commons Attribution-NonCommercial 4.0 International License (CC BY-NC 4.0)
METHOD FOR CALCULATING ACOUSTIC STRESSES IN SIX-BEAM DIFFRACTION IN LAYERED MEDIA
Yu.N. Belyayev
Syktyvkar State University, Syktyvkar, Russian Federation
ABSTRACT
The stresses arising in a layered medium as a result of an acoustic wave are investigated theoretically. In the general case, under the action of an incident elastic wave in an anisotropic layer, six waves are formed, three of which are directed to the reflection region and three of them are directed to the region of transmission. The stress-strain state of the layer is caused by the combined effect of these waves and is described by the equations of motion of a continuous medium and the generalized Hooke's law. This system of differential equations is solved with respect to the components of the displacement vector and the stress tensor in the Cartesian coordinate system in the matrix form. The components of the displacement vector and the stress tensor at two opposite boundaries of a layer of thickness d, are expressed through each other by means of a sixth-order transfer matrix T= exp(W d,). The calculation of this exponential is carried out using polynomials of the principal minors of the matrix W, and does not require finding the eigenvalues of the matrix W. This method provides a more accurate and reliable calculation of the transfer matrix of the W-layer medium T = TNTN-1^T1 in comparison with other known algorithms.
The amplitudes of the waves scattered by the anisotropic layer are expressed in terms of the elements of the transfer matrix. The distribution of acoustic stresses along the thickness of an anisotropic layer is determined by the amplitudes of the scattered waves and the elements of the corresponding transfer matrices. This method of calculating acoustic stresses is demonstrated for the incident SH-, SV- and P-type waves on the three-layer model: isotropic layer-crystal layer-isotropic layer.
We present the comparison of the scattering spectra of elastic waves and the dependence of the stresses on the scattering angles for the crystalline layers of silicon and lead molybdate. The interpretation of the resonances of acoustic stresses arising in the crystalline layer due to the action of shear waves is given.
©PNRPU
ARTICLE INFO
Received: 03 July 2018 Accepted: 09 September 2018 Published: 28 December 2018
Keywords:
elastic waves, matrix exponential, scaling and squaring method, polynomials of principal minors, recurrence relations, relative truncation error, stress resonances.
Введение
Определение упругих постоянных тонких пленок [1, 2], исследование критических напряжений в пленках [3, 4], разработка новых электроакустических и акусто-оптических преобразователей [5, 6] и некоторые другие проблемы тесно связаны с расчетом напряженно-деформированного состояния пленки, вызванного упругими волнами. Распределение напряжений в твердом слое зависит от вида возможных упругих волн. Как известно [7, 8, 9], в любом направлении изотропной среды могут распространяться упругие волны трех типов, различающихся направлением колебаний частиц среды. В продольной волне (обозначаемой как волна Р-типа) колебания происходят вдоль направления волнового вектора. В волнах сдвига амплитудный вектор колеблется перпендикулярно волновому вектору. Если в сдвиговой волне колебания направлены перпендикулярно или параллельно плоскости падения, то волна называется соответственно горизонтально поляризованной ^Н-тип) и вертикально поляризованной (SV-тип) волной сдвига.
При падении продольной или сдвиговой волны на слой твердого тела, часть энергии отражается и проходит через слой в виде волн того же типа, что и падающая волна. В дополнение к этим волнам изотропный слой может преобразовать Р-волну в вертикальную волну сдвига и наоборот [8, 10]. Анизотропия упругих
свойств среды распространения приводит к возможности взаимного преобразования волн всех трех типов. Соответствующие исследования для сред, обладающих орторомбической и тригональной симметриями проводились, например, в работах [11, 12]. В любом направлении анизотропной среды скорости распространения трех волн могут различаться друг от друга, и их поляризации не являются ни чисто продольными, ни чисто поперечными. Волна, у которой направления колебаний наиболее близки к направлению волнового вектора, называется квазипродольной. Две другие волны называются квазипоперечными [5]. В анизотропном слое конечной толщины ё в результате воздействия плоской упругой волны образуются шесть волн (по три в направлении отражения и пропускания), суммарное воздействие которых формирует напряженное состояние слоя.
В расчетах распространения волн в слоистых средах получил развитие подход, основанный на матричном методе решения системы дифференциальных уравнений [13, 14]. Описание этого метода применительно к упругим волнам имеется в обзорах [11, 15] и монографиях [16-18]. Согласно этому методу значения компонент вектора смещения и тензора напряжений на одной границе слоя выражаются через аналогичные компоненты на противоположной границе с помощью так называемой матрицы переноса Т = ехр(^У), где Ж - матрица
п-го порядка, составленная из коэффициентов соответствующей системы дифференциальных уравнений. Относительно просто решается задача рассеяния горизонтальной волны сдвига изотропным слоем [18], которая описывается матрицей переноса второго порядка. Для матрицы переноса четвертого порядка, определяющей рассеяние Р-8У-волн изотропным слоем, также известна аналитическая формула [19]. Эти матрицы переноса, как и их аналоги из оптики (см. например [20, 21]), находятся вычислением матричной экспоненты ехр(Жё) по формуле Лагранжа-Сильвестра [13] с помощью собственных значений /.,, ] = \,...,п, матрицы й7 в виде функции Т = Т(ехр(к1ё),...,ехр(кпс1)). В случае матрицы Ж общего вида, порядок которой п > 4, решение алгебраической задачи на собственные значения возможно только численными методами. Поэтому результат вычислений матрицы переноса по формуле Лагранжа-Сильвестра отягощен погрешностями вычислений собственных значений и экспонент ехр(к).
Похожие проблемы возникают, если для вычислений матричной экспоненты воспользоваться методом Бэке-ра-Вандермонда [22], интерполяционной формулой Ньютона [23], методом Лапласа [24], канонической формулой Жордана [13]. Указанные погрешности накапливаются при вычислении матрицы переноса всей структуры Т = ТХТХ ■■■Т1 с ростом числа слоев N и величин параметров распространения к^ = к, где й1, - это толщины слоев. В рассматриваемой задаче к] имеют смысл волновых чисел, и, следовательно, размерность модуля параметра распространения пропорциональна произведению частоты излучения и толщины. Возрастание именно этих величин может приводить к неудовлетворительным результатам при вычислении матрицы переноса и зависящих от нее величин (см. статью [25] и библиографические указания в ней). Это одна из причин, по которым перечисленные выше именные методы, известные из учебной и справочной литературы, практически не применяются для вычислений матричных экспонент большого порядка.
Проблема вычисления матричной экспоненты является центральной при решении больших систем дифференциальных уравнений в различных приложениях. Это стимулирует развитие известных методов расчетов, сравнительный анализ которых имеется в работах [26, 27], и разработку новых подходов, примеры которых рассмотрены в статьях [28-35]. Из набора известных методов наилучшим численным подходом к нахождению матричной экспоненты ехр(Жё), по мнению авторов обзора [28], является алгоритм масштабирования и кратного квадрирования (МКК): 1) масштабирование матрицы Жё ^ Жё / (от2) с целью понижения ее нормы, 2) вычисление экспоненты масштабированной матрицы ехр(Жё/(от2)), 3) вычисление Т = [ехр(Жё/(от2))]"2
с помощью кратного квадрирования. В методе МКК при выполнении второго этапа подразумевается использование аппроксимации Паде или аппроксимации разложения в ряд Тейлора конечной суммой. Примеры реализаций этих способов расчетов представлены в работах [36, 37]. Метод расчета [36] реализован в математическом пакете МАТЬАВ и используется некоторыми исследователями [35] для сравнительного тестирования своих алгоритмов.
В данной работе для вычисления матрицы переноса упругих волн применяется метод многочленов главных миноров [38], который, по крайней мере для матриц порядка п < 9 [39], обеспечивает более надежное и точное вычисление матричной экспоненты в сравнении с методами [36, 37].
1. Метод вычислений
1.1. Основные уравнения
Динамика напряженно-деформированного состояния упругой среды описывается уравнениями движения и законом Гука. Эти уравнения в декартовой системе координат выражаются формулами
^ Ч = £ да,
„=1
8
3 3
дг
ды± + д^ | = уу дХ дХ ) 8=1 у=1 к,1, у = 1,2,3,
(1)
где р - плотность; г - время; х, X, X - декартовы координаты; и, а, и к,1, у = 1,2,3) - компоненты вектора смещений, тензора напряжений и тензора упругой податливости соответственно.
Уравнения (1) решаются применительно к плоскослоистой среде, у которой плотность и упругие параметры среды зависят только от одной координаты х3 вдоль оси, перпендикулярной поверхности анизотропного слоя (рис. 1):
р = р( Х\ = (х3). (2)
Деформации в анизотропном слое 0 < х < ё порождаются плоской волной:
Ыо = Ао ехр[/'(к0 • г -ю/)], (3)
падающей из области х < 0 . Здесь ц0, А, к = = к01 е + к02е2 + к03е3 и ю являются соответственно вектором смещений, амплитудным вектором, волновым вектором и циклической частотой падающей волны; - мнимая единица.
Решение задачи (1)-(3) методом разделения переменных показывает, что все величины и и а имеют
одинаковую зависимость от переменных х, х2 и ' в виде ехр[г'(к01х + к02х2 - юг)], как в падающей волне. Поэтому если обозначить через (х,) (¿ = 1,...,9) неизвестные зависимости компонент вектора смещений и тензора напряжений от координаты х, то указанные компоненты принимают вид
= Vj (X3)eXP[i(k01 X1 + k02x2 -<B0L
(4)
где индекс} принимает значения от 1 до 9 в том же порядке, в каком перечислены компоненты вектора смещений и тензора напряжений в левой части равенства (4).
С учетом равенств (4) уравнения (1) редуцируются в систему уравнений
dТ(x) _
dx-
Т (x3) =
= WТ(x),
V1( X3)
V2( X3)
V3( X3) V4( X3) V5( X3) V6( X3)
W=
W11 W12 W13 W14 W15 W16
W21 Wii W23 0 W25 W26
W25 W15 W33 W34 W35 W36
W26 W16 W43 0 W36 W46
W23 W13 W53 0 W33 W43
0 W14 0 W64 W34 0
V (X3) = aa V (X3) + ai2y2 (X3) + any3(X3) +
+ ai5y5(X3) + a,6Уб(X3), i = 7,8,9,
(5)
(6)
где м>п = -¿(к01а92 + к02а82), м>12 = Б14а72 + Б24а82 + Б46а92 +
+ Б44 , ^13 =-'(к01а72 + к02а92\ ™14 = -iк02, ™15 = Б14а75 + + ^24а85 + ^4ба95 + Б45, ^16 = Б14 а76 + Б24 а86 + Б46 а96 + Б34 ,
™21 = к021 Уб + к022У4 - 2к01к02Уз - РЮ > ™23 = 4Уз - ^2У 5 + + к01к02(Уб У1 ™25 =-'(к01а95 + к02а85 X ^2б =-'(к01а9б + + к02а8б X ™33 = -'(к01а75 + к02а95 ), ™34 = -''к01 , ™35 = = Б15а75 + Б25а85 + Б56а95 + Б55, ™36 = Б15а76 + Б25а86 + Б56а96 + + Б35 , ^43 = -'(к01а7б + к02а9б ), ^46 = ^13Я7б + Б23а86 + + ^6 + ^33, ^53 = к021У1 + к022У6 + 2к01к02У3 -Р®2= ™64 = -Р®2
и использованы обозначения: - компоненты тензора упругой податливости в матричной форме [40];
а71 = '(к01У3 - к02У2 а72 = Б24У 2 - Б14У1 - »Я^ «73 =
= '(к01У1 + к02У3Х а75 = Б25У2 - Б15У1 - S5бУз, а76 = Б23У2 -- Б13У1 - Б36 У3> а81 = ' (к02У4 - к01У5), 082 = Б25У2 - Б24У4 + + Б46У5 , = 'к01У2 - 'к02У5> °85 = Б15У2 - Б25У4 + Б56У5 >
а86 = Б13У2 - Б23У4 + Б36У5' а91 = - 'к02У5 > а92 = Б24 У5 -
- 514У3 - 546Уб, а93 = 'к01У3 + 'к02Уб, а95 = Б25У5 - Б15У3 - Б56Уб, а96 = Б23У5 - Б13У3 - Б36У6; У7- =А7- Д^А - Б12Д2 + Б16Aз),
7=1,.. ,6; А1=522566-5226, А2 = 512566-5165,26,
^3 = Б12Б26 - Б16Б22, = Б11Б66 - Б16, ^5 = Б11Б26 - Б16Б12,
^6 = Б11Б22 - Б12 '
1.2. Матрица переноса
Решение системы уравнений (5) в матричной форме имеет вид
T(d2) = T T(dj),
d > d,
(7)
где
T = / +
d 2 d 2
J W (X )dX + J W (X) J W (^ )d ^ dX3
} Г(х3) J ) J )d^dxз +...
(8)
в теории дифференциальных уравнений называется матрицантом, а в теории волн в слоистых средах - матрицей переноса. Здесь и далее I обозначает единичную матрицу шестого порядка.
Если структура толщиной d состоит из N однородных слоев, т.е. W = W = const, dw < x3 < dt, то = ^(0), где T = TNTN_, ■■■Tl,
2]=exp[]Vi(di -d^)], i = \,...,N и матричная экспонента определяется формулой
i ша Wldl ^WJdJ exp(Wd) = I + Wd н---1---h... = > -. (9)
2! 3! £ /'
Соотношения (7)-(9) позволяют вьфазить значения функций \|/1(х3),...,\|/6(х3) через значения этих функций на глубине х3 = 0 . После этого нахождение значений функций у7(х3), у8(х3), ^ (х) по формулам (6) не представляет труда.
1.3. Вычисление матрицы переноса однородного слоя
Т(х) = ехр(^х3) в данной работе проводится по формуле
T(X3) =
(WX3 )h f 1 6^1 Л в
SI W"1 I + S -Л X P6-h+ sBj-1-g (6)
,=6 j!g=0
(10)
Здесь рг] = 1.....(■>. являются коэффициентами характеристического уравнения матрицы Жх3 / т. Как известно, коэффициенты характеристического уравнения матрицы равны с точностью до знака суммам соответствующих главных миноров этой матрицы, в частности
U2, ^23 , М1,М3,СТ13,СТ33,СТ11,СТ22,СТ12 =
+
d
d
d
d
d
m
Р равен следу матрицы, а р6 =- йе^Жх3 / т). Вследствие этого функции В. (6), определяемые рекуррентными уравнениями [38]
Вх (6) = В2(6) = ■■■ В4 (6), Я5(6) = 1,
6 (11) В у (6) = £ Р,Ву(6), ] > 6, ( )
,=1
называются многочленами главных миноров.
Для выполнения расчетов по формулам (10)—(11) требуется предварительное нахождение матриц (Жх3 / т)2, (Жх / т)3, (Жх / т)4, (Жх3 / т)5 и коэффициентов р.. Эти коэффициенты вычисляются рекуррентно по методу Леверье [13]. Для каждой из найденных матриц (Жх3 / т)у определяется след и применяются формулы
jPj = ^ ~ А V!-----р^Л, 7 = 1, • • •, 6. (12)
Выражение (10) получено из точной формулы в результате усечения ряда ДО конечной суммы
^6 " (•••). Порождаемая таким усечением относительная погрешность б в вычислениях элементов матрицы переноса зависит от числа N1 и удовлетворяет неравенству [41]
s<-
N +
шц 1
(13)
(6 + ¿(6 + h)
где ^ = 8,5dю/ (mmin v), min v ~1/^pmax - минимальная скорость распространения упругих волн в рассматриваемом слое. Параметр масштабирования m яв-
ляется минимальным целым, которое выбирается в соответствии с условием ^ < 1. Пример оценки параметра масштабирования для кристаллического слоя кремния представлен в статье [41]. Выбор М = 14 обеспечивает вычисление матрицы переноса с помощью алгоритма (10)—(12) с двойной арифметической точностью 10-16, а Ni = 25 соответствует четверной арифметической точности 2,5 х10-32. Для реализации этого метода требуется выполнение четырех матричных перемножений. Для сравнения Алгоритм 2.3 [36], использующий для вычисления матричной экспоненты аппроксимацию Паде, требует выполнения шести матричных перемножений плюс решение одного матричного уравнения. При этом нормализованная относительная погрешность превышает 1016.
1.4. Геометрия рассеяния
Из равенств (4) следует, что лучи всех волн, возникающих в слоистой среде, лежат в одной плоскости (x, 4) (рис. 1, а), и проекции волновых векторов на оси x и x определяются соответственно формулами k.¡ = k01 = k0 sin 0О cos а, & = k02 = k0 sin 0О sin а. Если обозначить 0j угол между волновым вектором kj и осью х, то kj3 = kj cos 0j, где kj =| k j | и углы 0 j определяются законом преломления (отражения):
kj sin 0j = k sin 0. (14)
Рассмотрим случай, когда области х3 < 0 их > d, ограничивающие анизотропную слоистую структуру, являются изотропными средами.
а б в г
Рис. 1. К определению направляющих косинусов падающей волны: a - углы падения 0О и а; б - падающая волна SH-типа; в - падающая волна SV-типа; г - падающая волна Р-типа Fig. 1. To defining the direction cosines of an incident wave: a - angles of incidence 0O and а; b - SH-type of the incident wave; c - SV-type of the incident wave; d - P-type of the incident wave
Как известно, в любом направлении изотропной среды могут распространяться волны Р-, 8И- и SV-типа, скорости распространения которых ур = ю^(2ц + к)/р,
уБН = у^г = ю^/ц/р, определяются упругими постоянными Ламе: к и ц. В результате воздействия падающей волны на анизотропный слой в изотропных областях х < 0 и х > ё могут возникнуть от двух до шести волн (рис. 2):
U = A ■ exp[/(k • r -ю/)],
kj = кjft + kJ2e2 + kJ3e3, j = 1,... 6,
где А. и к. - соответственно амплитудные и волновые
векторы; индексами 1 и 4 обозначены горизонтальные волны сдвига, индексы 2 и 5 соответствуют вертикальным волнам сдвига, наконец, индексы 3 и 6 отмечают продольные волны. Волновые числа этих волн имеют
следующие значения: | к1 |=| к2 |= ю^р0 / ц0 ,
I к3 |= ^Р0/(2И0 + ^0)> | к4 |=| к5 |=
| к6 |= рЛ / (2цЛ + ) . На рисунках направления колебаний амплитудных векторов А] волн 8И-, 8У-
и Р-типа показаны стрелками зеленого, синего и красного цветов соответственно.
u = u4 + u5 + u6, y.d * 0;
Рис. 2. К определению направляющих косинусов рассеянных волн: Ai и A4 - амплитудные векторы SH-волн, А2 и А5 -амплитудные векторы SV-волн, Аз и Аб - амплитудные векторы Р-волн Fig. 2. To defining the direction cosines of the scattered waves, Ai and A4 are the amplitude vectors of SH waves, A2 and А5 are the amplitude vectors of SV waves, Аз and Аб are the amplitude vectors of P waves
1.5. Граничные условия
Упругие свойства изотропных сред (компоненты тензора упругой податливости) выражаются через постоянные Ламе. С помощью последних уравнения закона Гука (1) существенно упрощаются, и напряжения агз
в изотропных областях x < 0 и x3 > d выражаются через компоненты вектора смещения формулами
(
стй = ц
du,. du
кдх3
_-
dxi у
d u3
i = 1,2,
ст33 = (2ц+ Х)—1+ Х
dx
^ du du, ^ —1 + —-
vdxi
(15)
dx
где ^ = , ц = ц0 ' если х3 < 0 и Х = ^, Ц = Ц^,
если х > d.
Смещения частиц среды в этих изотропных средах обусловлены волнами, которые в них распространяются. Таким образом, в области х ^ 0
цо *
Цо = о,
(16)
в области х ^ d
u = Ш
Vj = 0
(17)
Подстановка векторов (16) и (17) в уравнения (15) позволяет найти компоненты тензора напряжений на границах анизотропного слоя. Следовательно, используя соотношения (3), (4), (15), (16) и (17), несложно выразить значения неизвестных функций у ] (х3) на границах X = 0 и х = Л анизотропного слоя через амплитуды рассеянных волн. Так, для случая, когда ц0 ^ 0, ^ 0, такой подход дает следующий результат:
(0) = С02 А0 + С12 А1 + С22 А2 + С32 А3 >
3
У2(0) = М-0 А Е (СИ2 кИ3 + СИ3кИ2)СИИ
И=0
У3 (0) = ^А + СпА + С2^2 + С31А >
У 4 (0) = С03 А + СпА + С23 А + С33 А3 ,
3
(0) = Ц А Е (СИ1кИ3 + СИ3кИ1 )СИ >
V6 (0) = A0 Z 2ЦоCh3kh3 + ^0 Z chgkhg
h=0 g=1
Vi(d ) = C42 л + c52 A + c62 A,
6
V2(d) = M0 Z (Ch2kh3 + Ch3kh2)Ch
C
V-(d ) = C41A4 + C51A5 + C6lA,
V4(j ) = C 43A4 + C53A5 + C63 A6,
6
V5(j ) = VdA0 Z (Ch1kh3 + Ch3kh1)Cih
V6(j ) = A0 Z
2VdCh3kh3 d Zchgkhg
C
(18)
где cgh - направляющие косинусы амплитудных векторов Ag волн. А, =|Л/;|, h = 1,...,6, Сш = 4,'Цр В последней формуле индексы i = 1, 2, 3 указывают три возможных типа падающей волны. Значение i = 1 соответствует SH-волне (см. рис. 1, б), в которой направляющие косинусы амплитудного вектора с01 = sin а, с02 =- cos а, с03 = 0; i = 2 соответствует SP-волне (см. рис. 1, в) с направляющими косинусами с01 = cos 90 cos а, с02 = cos 90 sin а, с03 = - sin 90; и i = 3 соответствует Р-волне (см. рис. 1, г) с направляющими косинусами с01 = sin 90 cos а, с02 = sin 90 sin а, с03 = cos 90. Рис. 2 определяет направляющие косинусы амплитудных векторов рассеянных волн: сп = sin а, с12 =- cos а, с13 = 0,
с21 = - cos 9 cos а, с22 = - cos 9 sin а, с23 = - sin 9, с31 = sin 9 cos а, с32 = sin 9 sin а, с33 =- cos 9, с41 = sin а, с42 =- cos а, с43 = 0, с51 = cos 9 cos а,
с52 = cos 9 sin а, с53 =- sin 9, с62 = sin 9 sin а, с63 = cos 9 •
с61 = sin 9 cos а,
h=0
h=4
h=4
u = u + u, + u + u3
u = u + u3
1.6. Вычисление коэффициентов преобразований волн
В результате подстановки выражений функций ^ (0), ^ (ё), из формул (18) в уравнение ) = = Т(ё (0) получается система алгебраических уравнений относительно амплитуд волн А. Решение этой системы, например, по методу Гаусса дает значения коэффициентов преобразований Сй падающей волны в рассеянные волны. Энергия волны пропорциональна квадрату модуля амплитуды. Поэтому величина Са2 характеризует долю энергии падающей волны, передаваемой к-й волне, и называется интенсивностью к-й волны.
1.7. Вычисление акустических напряжений
Вычиление акустических напряжений на глубине х3 производится по формуле х) = Т(х )Т(0),
где значения компонент матрицы Т (0) определяются из (18) с учетом найденных коэффициентов преобразований С .
2. Результаты расчетов
Представленный алгоритм расчетов был использован в исследовании акустических напряжений в кристаллических слоях. Некоторые результаты вычислений для слоев молибдата свинца и кремния показаны на рис. 3-6. Все значения амплитуд | а.3 | колебаний компонент тензора напряжений даны в единицах А1015 Па. Циклическая частота ю = 2л:107 Гц. Плоскость х3 = 0 совпадает с гранью кристалла (001). Параметры кристаллов, использованные в расчетах, представлены в таблице, и области х < 0, х3 < ё предполагались твердыми с пара-
метрами =Xd = 1,67-10 Н/м , х1010 Н/м2, p0=prf = 2,65-103 кг/м3.
Ц = Hd = 3,27 х
Параметры кристаллов Parameters of crystals
Кристалл P Sii Sl2 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Па
На рис. 3-5 представлены результаты расчетов напряжений на границах слоя РЬМо04 толщиной ё = 0,000924 м, возникающих при падении на этот слой горизонтальной волны сдвига (см. рис. 3), вертикальной волны сдвига (см. рис. 4) и продольной волны (см. рис. 5). Зависимости |агз (0) =|агз (0)|(0О) отмечены
на рисунках цифрами 1, 2, 3, соответствующими значению индекса / у компонент тензора напряжений. На рис. 3-5, а показаны зависимости энергетических характеристик С?(0О) волн, рассеянных в область х3 < 0 , а на рис. 3-5, б - соответствующие характеристики волн, рассеянных в область х > ё . На этих рисунках
кривые С2(0О) отмечены цифрами 1, 2, 3, 4, 5, 6, соответствующими значениям индекса ]. Кроме этого, функциональные зависимости Сд(00) и С24(00) интен-сивностей горизонтальных волн сдвига изображены зеленым цветом, интенсивности вертикальных волн сдвига С2(00) и С2(00) нарисованы синим цветом, а для показа интенсивностей С23(00) и С26(00) продольных волн использован красный цвет. Такие же цвета используются в рис. 6.
В результате перерассеянии волн на границах кристаллического слоя рассеянные волны могут усиливаться или, наоборот, ослабляться за счет интерференции. Число таких экстремумов возрастает с увеличением
толщины и уменьшением скорости распространения волн. Наряду с системой интерференционных максимумов и минимумов на рис. 3 и 4 видны экстремумы резонансного характера. Причиной их возникновения является следующее.
Рис. 3. Амплитуды напряжений, возникающих на границах слоя PbMoO4 под воздействием горизонтальной волны сдвига: а = 38° , d = 0,000924 м Fig. 3. Amplitudes of stresses arising at the boundaries of PbMoO4 layer under the influence of a horizontal shear wave: а = 38° , d = 0,000924 m
Рис. 4. Амплитуды напряжений, возникающих на границах слоя PbMoO4 под воздействием вертикальной волны
сдвига: а = 38° , d = 0,000924 м Fig. 4. Amplitudes of stresses arising at the boundaries of PbMoO4 layer under the action of a vertical shear wave: а = 38° , d = 0.000924 m
Рис. 5. Амплитуды напряжений, возникающих на границах слоя PbMoO4 под воздействием продольной волны:
а = 38° , d = 0,000924 м Fig. 5. Amplitudes of stresses arising at the boundaries of PbMoO4 under the action of a longitudinal wave: а = 38° , d = 0.000924 m
Если сдвиговая волна падает на кристалл под углом 9 > 90C (критический угол 0OC определяется в соответствии с законом (14) равенством sin 90С = k3 / k0), то
продольные волны, генерируемые кристаллом, становятся неоднородными, распространяющимися вдоль оси £ . Для рассматриваемой структуры 0ОС = 39,15°. При определенной толщине кристаллического слоя и углах падения суммарное смещение вертикальных волн сдвига внутри кристалла создает поперек слоя мо-
дулированную структуру. Дифракция продольных волн на этой структуре порождает резко выраженные экстремумы в спектрах рассеянных волн при некоторых углах падения 0О > 0ОС. Именно этим объясняются резко выраженные экстремумы напряжений в спектрах рассеянных волн. Корреляции между экстремумами кривых С = С2(0О), с одной стороны, и экстремумами
функций |стй(0)| = | ^/3 (0)| (00 ) и ^ =К-3(Л)(00) -с другой, определяются формулами (18) и хорошо видны при сравнении рис. 3, a и рис. 3, в, рис. 3, б и рис. 3, г, рис. 4, a и рис. 4, в, рис. 4, б и рис. 4, г.
На распределение энергии падающей волны между волнами, рассеянными кристаллом, большое влияние оказывает также вращение плоскости поляризации сдвиговой волны по мере прохождения кристаллической среды. Одним из проявлений этого эффекта являются равенства С\2 = С21 и С = С24, которые выполнялись с высокой точностью во всех численных экспериментах, проведенных в рамках данной работы и [41]. Иллюстрациями сказанного являются синяя кривая на рис. 3, a и зеленая кривая на рис. 4, a, синяя кривая на рис. 3, б и зеленая кривая на рис. 4, б. В результате вращения плоскости поляризации сдвиговой волны падающая волна 5^-типа ^К-типа) при отражении или пропускании полностью может превращаться в волну £К-типа (SH-типа). Пример такого преобразования горизонтальной волны сдвига кристаллическим слоем кремния показан на рис. 6. При углах падения а = 44°
и 0= 42,7232° коэффициент С2 = 0, т.е. отраженная от кристалла сдвиговая волна является вертикально поляризованной.
Рис. 6. Интенсивности волн, излучаемых слоем Si под воздействием горизонтальной волны сдвига а = 44°, d = 0,000845 м Fig. 6. The intensities of waves emitted by a Si layer under the influence of a horizontal shear wave а = 44°, d = 0.000845 m
Амплитуды всех других волн резонансно увеличены. Продольные волны являются неоднородными, их волновые векторы направлены вдоль поверхностей кристалла x3 = 0 и x3 = d. Эти волны энергию от кристаллического слоя не переносят. При указанных условиях на границах кристаллического слоя кремния возникают
резонансы напряжений. Амплитуды колебаний этих напряжений имеют следующие значения: |а13 (0) = 16,3,
|а23(0) = 15,8, |а33(0) = 4,8, |а13(ё)| = 16,4 , |а23(ё)| = 16,2, |а13(0) = 5,4. Ширина резонансов спектров рассеяния, показанных на рис. 6, и соответствующих им резонансов напряжений равна примерно пяти угловым секундам.
Заключение
Такие эффекты, как полное преобразование поляризации сдвиговой волны, дифракционные резонансы, которые возможны при шестилучевой дифракции упругих волн в анизотропном слое, происходят в очень узком диапазоне Д0О углов падения 0. Положение такого диапазона Д0О может быть изменено за счет не-
больших изменений значений второго угла падения а и параметра распространения к = ёю / у. Аналитический расчет таких эффектов для среды с тензором упругой податливости общего вида, по-видимому, невозможен. Результаты численного моделирования спектров рассеяния упругих волн и создаваемых ими напряжений зависят от точности расчетов матрицы переноса. Метод многочленов главных миноров позволяет находить эту матрицу с заданной точностью.
Важная особенность метода расчета, представленного в данной работе, состоит в том, что при исследовании распространения волн в многослойной среде требуется знание волновых чисел только для двух крайних, ограничивающих структуру слоев (на рис. 2 это слои х < 0 и х > ё). Эти волновые числа учитываются в граничных условиях, но не влияют на результат вычислений матриц переноса внутренних слоев Т1,...,ТЫ.
Библиографический список
1. Thickness Dependence of the Properties of Epitaxial Barium Strontium Titanate Thin Films / V.B. shirokov, Yu.I. Golovko, V.M. Mukhortov, Yu.I. Yuzyuk, P.E. Janolin, B. Dkhil // Physics of the Solid State. - 2015. - Vol. 57. - No. 8. - Р. 1529-1534. DOI: 10.1134/S1063783415080314
2. Material Constants of Barium Titanate Thin Films / V.B. Shirokov, V.V. Kalinchuk, R.A. Shakhovoi, Yu.I. Yuzyuk // Physics of the Solid State. - 2015. - Vol. 57. - No. 8. - Р. 15351540. DOI: 10.1134/S1063783415080302
3. Physics of thin films. Advances in research and development. Vol. 1. Ed. G. Hass. - New York: Academic Press, 1963. - 350 p.
4. Егоров Г.П., Волков А.А. Определение критического уровня внутренних напряжений в тонких пленках // Композиты и наноструктуры. - 2016. - Т. 8, № 3. - С. 187-203.
5. Dieulesaint E., Royer D. Ondes élastiques dans les solides. Application au traitment du signal. - Paris: Masson, 1974. - 424 c.
6. Блистанов А.А. Кристаллы квантовой и нелинейной оптики. - М.: Изд-во МИСИС, 2000. - 432 с.
7. Lamb G. The dynamical theory of sound. - London: Edward Arnold & Co, 1931. - 307 p.
8. Ewing W.M., Jaidetzky W.S., Press F. Elastic waves in layered media. - New York: McGraw-Hill Book Company, 1957. - 380 p.
9. Виноградова М.Б., Руденко О.В., Сухоруков А.П. Теория волн. - М.: Наука, 1979. - 384 с.
10. Zhag P., Wei P., Li Y. Wave propagation through a micropolar slab sandwiched by two elastic half-spaces // Journal of Vibration and Acoustics. - 2016. - Vol. 138. - Р. 04100804100817. D0I:10.1115/1.4033197
11. Crampin S. A review of wave motion in anisotropic and cracked elastic media // Wave Motion. - 1981. - Vol. 3 - Р. 343391. DOI: 10.1016/0165-2125(81)90026-3
12. Polikarpova N.V., Mal'neva P.V., Voloshinov V.B. The ani-sotropy of elastic waves in a tellurium crystal // Acoustical Physics. -2013. - Vol. 59. - Р. 291-296. DOI: 10.1134/S1063771013010144
13. Gantmacher F.R The Theory of Matrices. Vol. 1. - New York: Chelsea, 1959. - 276 p.
14. Michal A.D. Matrix and tensor calculus with application to mechanics? Elasticity and aeronautics. - New York: Dover Publication, 2008. - 132 p.
15. Беляев Ю.Н. Методы вычислений матриц переноса упругих деформаций // Вестник Пермского национального
исследовательского политехнического университета. Механика. - 2013. - No. 3. - С. 63-109.
16. Молотков Л.А. Матричный метод в теории распространения волн в слоистых упругих и жидких средах. - Л.: Наука, 1984. - 201 с.
17. Brekhovskikh L.M., Godin O.A. Acoustics of layered media. - Berlin: Springer-Verlag, 1990. - 416 с.
18. Aki K., Richards P.G. Quantitative seismology, Sausali-to. - CA: University Science Books, 2002. - 700 p.
19. Thomson W.T. Transmission of elastic wave through a stratified solid material // J. Appl. Phys. - 1950. - Vol. 21. -Р. 89-93. DOI: 10.1063/1.1699629
20. Knittl Z. Optics of thin films. - London: J. Wiley, 1975. -548 p.
21. Abdulhalim I. Analytic propagation matrix method for anisotropic magnetooptic layered media // J. Opt. A: Pure Appl. Opt. - 2000. - Vol. 2. - Р. 557-564.
22. Angot A. Compléments de mathématiques a l'usage des ingénieurs de l'élektrotechnique et des telecommunications. - Paris: Masson, 1982. - 868 p.
23. Dehghan M., Hajarian M. Determination of a matrix function using the divided difference method of Newton and the interpolation technique of Hermite // Journal of Computational and Applied Mathematics. - 2009. - Vol. 231. - Р. 67-81. DOI: 10.1016/j.cam.2009.01.021
24. Bellman R. Introduction to matrix analysis. - New York: McGaw-Hill Book Company, 1960. - 348 p.
25. Potel C., Gatignol P., de Belleval J.-F. Energetic criterion for the radiation of floquet waves in infinite anisotropic periodically multilayered media // Acta Acustica. Acustica. - 2001. -Vol. 87. - Р. 340-351.
26. Moler C., Van Loan C. Nineteen dubious ways to compute the exponential of a matrix, Twenty-five years later // SI-AM Review. - 2003. - Vol. 45. - Р. 3-49. DOI: 10.1137/S00361445024180
27. Higham N.J. Functions of matrices. Theory and computations. - Philadelphia: SIAM, 2008. - 425 p.
28. Popolizio M., Simoncini V. Acceleration techniques for approximating the matrix exponential operator // SIAM J. Matrix Anal. Appl. - 2008. - Vol. 30. - Р. 657-683. DOI:10.1137/060672856
29. Hochbruck M., Ostermann A. Exponential integrators // Acta Numer. - 2010. - Vol. 19. - P. 209-286. DOI: 10.1017/S0962492910000048
30. Al-Mohy A., Higham N. Computing the action of the matrix exponential, with an application to exponential integrators // SIAM J. Sci. Comput. - 2011. - Vol. 33. - Iss. 2. - P. 488-511 DOI: 10.1137/100788860}
31. Eiermann M., Ernst O., Guttel S. Deflated restarting for matrix functions // SIAM J. Matrix Anal. Appl. - 2011. - Vol. 32. -Iss. 2. - P. 621-641. DOI: 10.1137/090774665
32. Soodhalter K., Szyld D., Xue F. Krylov subspace recycling for sequences of shifted linear systems // Appl. Numer. Math. -
2014. - Vol. 81. - P. 105-118. DOI: 10.1016/j.apnum.2014.02.006
33. Archid A., Bentbib A. Approximation of the matrix exponential operator by a structure-preserving block Arnoldi-type method // Appl. Numer. Math. - 2014. - Vol. 75. - P. 37-47. DOI: 10.1016/j.apnum.2012.11.008
34. Soodhalter K. Block Krylov subspace recycling for shifted systems with unrelated right-hand sides // SIAM J. Sci. Comput. -2016 - Vol. 38. - P. A302-A324. DOI: 10.1137/140998214
35. Wu G., Pang H.-K., Sun J.-L. A shifted block FOM algorithm with deflated restarting for matrix exponential computations //
References
1. Shirokov V.B., Golovko Yu.I., Mukhortov V.M., Yuzyuk Yu. I., Janolin P.E., Dkhil B. Thickness dependence of the properties of epitaxial barium strontium titanate thin films. Physics of the .solid state, 2015, vol. 57, no. 8, pp. 1529-1534. DOI: 10.1134/S1063783415080314
2. Shirokov V. B., Kalinchuk V.V., Shakhovoi RA., Yuzyuk Yu.I. Material constants of barium titanate thin films . Physics of the solid state,
2015, vol. 57, no. 8, pp. 1535-1540. DOI: 10.1134/S1063783415080302
3. Physics of thin films. Advances in research and development. vol. 1. Ed. G. Hass. New York, Academic Press, 1963, 350 p.
4. Egorov G.P., volkov A.A. Opredeleniye kriticheskogo urovnya vnutrennikh napryazheniy v tonkikh plonkakh [Determination of the critical level of internal stresses in thin films]. Composites and nanostructures, 2016, vol. 8, no. 3, 187-203.
5. Dieulesaint E., Royer D. Ondes élastiques dans les solides. Application au traitment du signal. Paris, Masson, 1974, 424 p.
6. Blistanov A.A. Kristally kvantovoy i nelineynoy optiki [Crystals of quantum and nonlinear optics]. Moscow, MISIS, 2000, 432 p.
7. Lamb G. The dynamical theory of sound. London. Edward Arnold & Co, 1931, 307 p.
8. Ewing W.M., Jardetzky W.S., Press F. Elastic waves in layered media. New York, McGraw-Hill Book Company, 1957, 380 p.
9. Vinogradova M.B., Rudenko O.V., Sukhorukov A.P. Teor-iya voln [Theory of waves]. Moscow, Nauka, 1979, 384 p.
10. Zhag P., Wei P., Li Y. Wave propagation through a micropolar slab sandwiched by two elastic half-spaces. Journal of vibration and acoustics, 2016, vol. 138, pp. 041008-041008-17. DOI:10.1115/1.4033197
11. Crampin S. A review of wave motion in anisotropic and cracked elastic media. Wave Motion, 1981, vol. 3, pp. 343-391. DOI: 10.1016/0165-2125(81)90026-3
12. Polikarpova N.V., Mal'neva P.V., voloshinov V.B. The ani-sotropy of elastic waves in a tellurium crystal. Acoustical Physics, 2013, vol. 59, pp. 291-296. DOI: 10.1134/S1063771013010144
13. Gantmacher F.R The Theory of Matrices. New-York, Chelsea, 1959, vol. 1, 276 p.
14. Michal A.D. Matrix and tensor calculus with application to mechanics? Elasticity and aeronautics. New-York, Dover Publication, 2008, 132 p.
15. Belyayev Yu.N. Methods for computing transfer matrices of elastic deformations. PNRPU Mechanics Bulletin, 2013, no. 3, pp. 63-109.
16. Molotkov L.A. Matrichnyy metod v teorii rasprostraneniya voln v sloistykh uprugikh i zhidkikh sredakh [The matrix method in the theory of
Applied Numerical Mathematics. - 2018. - Vol. 127. - P. 306323. DOI: 10.1016/j.apnum.2018.01.015
36. Higham N.J. The scaling and squaring method for the matrix exponential revisited // SIAM Review. - 2009. - Vol. 51. -Р. 747-764. DOI: 10.1137/090768539
37. New Scaling-Squaring Taylor Algorithms for Computing the Matrix Exponential / J. Sastre, J. Ibanez, E. Defez, P. Ruiz // SIAM J. Sci. Comput. - 2015. - Vol. 37. - No. 1. - Р. A439-A455. DOI: 10.1137/090763202
38. Belyayev Yu.N. On the calculation of functions of matrices // Mathematical Notes. - 2013. - Vol. 94. - No. 2. - Р. 177184. DOI: 10.1134/S0001434613070171
39. Belyayev Yu.N. On the calculation of matrix exponential of a large order // Proceedings of the International Conference DAYS on DIFFRACTION, St. Petersburg. - 2017. - Р. 55-59. DOI: 10.1109/DD.2017.8167995
40. Сиротин Ю.И., Шаскольская М.П. Основы кристаллофизики. - М.: Наука, 1979. - 640 с.
41. Belyayev Y.N. Conversion of elastic waves as a result of diffraction in anisotropic layer // IOP Conf. Series: Materials Science and Engineering. - 2017. - Vol. 208. - Р. 012003-012993-7. DOI:10.1088/1757-899X/208/1/012003
wave propagation in layered elastic and liquid media]. Leningrad, Nauka, 1984, 201 p.
17. Brekhovskikh L.M., Godin O.A. Acoustics of layered media. Berlin, Springer-Verlag, 1990, 416 p.
18. Aki K., Richards P.G. Quantitative seismology. Sausali-to, CA, University Science Books, 2002, 700 p.
19. Thomson W.T. Transmission of elastic wave through a stratified solid material. J. Appl. Phys, 1950, vol. 21, pp. 89-93. DOI: 10.1063/1.1699629
20. Knittl Z. Optics of thin films. London, J. Wiley, 1975, 548 p.
21. Abdulhalim I. Analytic propagation matrix method for anisotropic magnetooptic layered media. J. Opt. A: Pure Appl. Opt, 2000, vol. 2, pp. 557-564.
22. Angot A. Compléments de mathématiques a l'usage des ingénieurs de l'élektrotechnique et des telecommunications. Paris, Masson, 1982, 868 p.
23. Dehghan M., Hajarian M. Determination of a matrix function using the divided difference method of Newton and the interpolation technique of Hermite. Journal of computational and applied mathematics, 2009, vol. 231, pp.67-81. DOI:10.1016/j.cam.2009.01.021
24. Bellman R. Introduction to matrix analysis, New York, McGaw-Hill Book Company, 1960, 348 p.
25. Potel C., Gatignol P., , de Belleval J.-F. Energetic criterion for the radiation of floquet waves in infinite anisotropic periodically multi-layered media. Acta acustica. Acustica. 2001, vol. 87, pp. 340-351.
26. Moler C., Van Loan C. Nineteen dubious ways to compute the exponential of a matrix, Twenty-five years later. SIAM Review, 2003, vol. 45, pp. 3-49. DOI: 10.1137/S00361445024180
27. Higham N.J. Functions of matrices. Theory and computations. Philadelphia, SIAM, 2008, 425 p.
28. Popolizio M., Simoncini V. Acceleration techniques for approximating the matrix exponential operator. SIAM J. Matrix Anal. Appl, 2008, vol. 30, pp. 657-683. DOI:10.1137/060672856
29. Hochbruck MM., Ostermann A. Exponential integrators. Acta Numer, 2010, vol. 19, pp. 209-286. DOI:10.1017/S0962492910000048
30. Al-Mohy A., Higham N. Computing the action of the matrix exponential, with an application to exponential integrators. SIAM J. Sci. Comput, 2011, vol. 33, iss 2, pp. 488-511 DOI: 10.1137/100788860}
31. Eiermann M., Ernst O., Guttel S. Deflated restarting for matrix functions. SIAM J. Matrix Anal. Appl, 2011, vol. 32, iss. 2, pp. 621-641. DOI: 10.1137/090774665
32. Soodhalter K., Szyld D., Xue F. Krylov subspace recycling for sequences of shifted linear systems. Appl. Numer. Math, 2014, vol. 81, pp. 105-118. DOI: 10.1016/j.apnum.2014.02.006
33. Archid A., Bentbib A. Approximation of the matrix exponential operator by a structure-preserving block Arnoldi-type method. Appl. Numer. Math, 2014, vol. 75, pp. 37-47. DOI: 10.1016/j.apnum.2012.11.008
34. Soodhalter K. Block Krylov subspace recycling for shifted systems with unrelated right-hand sides. SIAMJ. Sci. Com-put. 2016, vol. 38, pp. A302-A324. D0I:10.1137/140998214
35. Wu G, Pang H.-K., Sun J.-L. A shifted block FOM algorithm with deflated restarting for matrix exponential computations. Applied Numerical Mathematics, 2018, vol. 127, pp. 306-323. DOI: 10.1016/j.apnum.2018.01.015
36. Higham N.J. The scaling and squaring method for the matrix exponential revisited. SIAM Review, 2009, vol. 51, pp. 747764. DOI: 10.1137/090768539
37. Sastre J., Ibanez J., Defez E., Ruiz P. New Scaling-Squaring Taylor Algorithms for Computing the Matrix Exponential. SIAM J. Sci. Comput, 2015, vol. 37, no. 1, pp. A439-A455. DOI: 10.1137/090763202
38. Belyayev Yu.N. On the calculation of functions of matrices. Mathematical Notes, 2013, vol. 94, no. 2, pp. 177-184. DOI: 10.1134/S0001434613070171
39. Belyayev Yu.N. On the calculation of matrix exponential of a large order. Saint Petersburg, Proceedings of the International Conference DAYS on DIFFRACTION, 2017, pp. 55-59. DOI: 10.1109/DD.2017.8167995
40. Sirotin YU.I., Shaskol'skaya M.P. Osnovy kristallofiziki [Fundamentals of Crystallophysics]. Moscow, Nauka, 1979,640 p.
41. Belyayev Y.N. Conversion of elastic waves as a result of diffraction in anisotropic layer. IOP Conf. Series: Materials Science and Engineering, 2017, vol. 208, pp. 012003-012993-7. DOI:10.1088/1757-899X/208/1/012003
