CC BY
17Вестник Сыктывкарского университета. Сер. 1. Вып. 18.2013
УДК 517.926+538.951
МАТРИЦА ПЕРЕНОСА УПРУГИХ ДЕФОРМАЦИЙ В
КРИСТАЛЛАХ
Ю.Н. Беляев, С.А. Попов
Дифференциальные уравнения упругих волн в кристаллах решаются с помощью симметрических многочленов шестого порядка и метода масштабирования. Исследовано влияние толщины слоя и частоты волны на масштабирующий фактор. Получено аналитическое решение, описывающее перенос упругих напряжений в кристаллическом слое кубической сингонии. Ключевые слова: слоистые среды, волны, матрица, симметрические многочлены, погрешность усечения, масштабирование.
1. Введение
Теоретические исследования распространения волн упругих напряжений в твёрдом теле основаны на использовании тензора Коши, закона Гука и уравнения движений в сплошной среде, которые могут быть записаны, соответственно, следующим образом:
1 / ¿и „Л
£=г( Тг + (1>
3 3 3 3
Рц = ^ ^ или е^ = ^ ^ Б^ы'Ры, (2)
к=1 1=1 к=1 1=1
д2и
= ^Р, (3)
где = 1, 2, 3, ё = (е^) — линейный тензор деформации, и — вектор смещения, г - радиус-вектор точки, Р = (рц) — тензор напряжений, Сук1 — компоненты тензора коэффициентов упругой жёсткости, Б^ы —
© Беляев Ю.Н., Попов С.А., 2013.
компоненты тензора коэффициентов упругой податливости, р — плотность тела, Ь — время.
Вычисления матрицы переноса упругих волн в изотропном слое [1] основано на представлении смещений с помощью скалярного и векторного потенциалов. Другой ПОДХОД состоит в преобразовании правой части (3) с помощью (1) и первого из соотношений (2) и получении, таким образом, дифференциального уравнения второго порядка относительно смещений. Волновые числа плоских волн, удовлетворяющих последнему уравнению, ищутся из решения уравнения Кристоффеля. Этот метод был использован применительно к трансверсально-изотропным средам [2, стр. 149]. Расчёт этим методом волновых характеристик в кристаллах [3] в общем случае достаточно громоздок и возможен только численными методами.
В данной работе развивается метод вычисления матрицы переноса упругих напряжений в многослойных кристаллических плёнках, основанный на решении уравнений упругих волн методом симметрических многочленов (МСМ) [4]. Система (1) - (3) решается в декартовых координатах х\ = х, х2 = у, х3 = г для кристаллического слоя, ограниченном плоскостями г = 0 и г = а, при условии, что в начальные моменты времени Ь < Ь0
и0 ехр[г(к0 •г — шЬ)], для г< 0,
и = 1 п \ л ^
[ 0, для г > а.
Здесь и0 и ш — соответственно амплитуда и циклическая частота колебаний плоской волны (4), падающей из однородного полупространства г < 0 та кристаллический слой; I — мнимая единица; волновой вектор к0 лежит в плоскости хг и направлен под углом 90 к ос и 0г.
2. Определяющее уравнение матрицы переноса
Для упрощения записи соотношений упругости часто используется так называемые матричные обозначения [5, стр. 618], согласно которым пары индексов ] заменяются одним индексом, скажем д, по схеме:
11 22 33 23 31 12
i = 1 , 1 , 1 , 1 , 1 , 1 •
q 1 2 3 4 5 6
^ Sn, S1232 ^ S64 S3313 ^ S35 и т.д. При такой
компонент тензора Sjki на двухиндексные обозначения Sqr тензорIdI дО"* формации и напряжений преобразуются следующим образом: е11 ^ е1}
е12 = е21 ^ 1 е6, е13 = е31 ^ 2е5) е22 ^ е2, е23 = е32 ^ 2е4> е33 ^ е3';
Pli ^ Pl, Р12 = P21 ^ Рб, Р13 = P31 ^ P5, P22 ^ P2, P23 = P32 ^ P4, p33 ^ P3. В результате вторая из формул (2), принимает вид
iq =
6
"Y^SqrPr , ^ = 1
r=1
или, с учётом определения (1),
dux
dx duy,
$11 рхх + S12Pyy + S13 Pzz + S14 Pyz + S15 Pzx + S16Pxy, (5)
S21 Рхх + S22Pyy + S23 Pzz + S24 Pyz + S25 Pzx + S26Pxy, (6)
S31 Рхх + S32Pyy + S33 'Pzz + S34 Pyz + S35 Pzx + S36Pxy, (7)
S41 рхх + S42Pyy + S43 Pzz + S44 Pyz + S45 pzx + S46Pxy, (8)
S51 рхх + S52Pyy + S53 Pzz + S54 Pyz + S55 Pzx + S56Pxy, (9)
S61 рхх + S62Pyy + S63 Pzz + S64 Pyz + S65 Pzx + S66Pxy ■ (10)
дУ дщ
дг
диу диг дг ду диг дих дх дг
дих диу дУ дх
Система из девяти уравнений (3) и (5)^(10) содержит 9 неизвестных функций: иу,Руг,их,иг,РхХ,Рхх,рхх,руу,рху■ Будем полагать, что плотность тела р и элементы матрицы коэффициентов упругой податливости Бд] зависят только от координаты г\ р = р(г), Бд] = Бд](г). Решая систему уравнений (3) и (5)^(10) методом разделения переменных, несложно показать, что искомые величины имеют такую же функциональную зависимость от х, у и Ь, что и падающая плоская волна, т.е. в рассматриваемом случае,
иу, Руг, их, и2, Ргх, Ргг, Рхх, Руу, Рху = Ф] (г) ехр[1(к0Х вШ 0О - шЬ)], (11)
где номер ] соответствует последовательности функций в (11).
Подставляя выражения (11) в уравнения (3), (5)—(10) получаем уравнения, определяющие функции ф1, ■ ■ ■ , ф9:
аФ
аФ = ^ Ф, (12)
аг
6 6 6 ф7 = ^2 а7]ф], ф8 = ^2 а8]ф], фд = ^2 Фз, (13)
3 = 1 3 = 1 3 = 1
6
где Ф — вектор-столбец, элементами которого являются функции 01(г),..., "06 (¿0; элемент ы ■ид. матрицы Ж и коэффициенты а7з- ,а8., а9., д] = 1,..., 6, алгебраически выражаются через компоненты БдГ тензора упругой податливости.
Очевидно, что когда решение уравнения (12) найдено, определение функций 07,08,09 по формулам (13) не составляет труда.
3. Метод вычисления матрицы переноса
Матрица переноса Т связывает значения вектор-функции Ф на границе г = 0 и глуби не г, т.е. является решением задачи Коши для уравнения (12)
Ф(г) = ТФ(0).
Если слой однородный, то
Т = ехр(Жг). (14)
Если пластина толщиной й состоит из N слоёв, первый из которых лежит в области 0 > г > г1; второй — г1 > г > г2 и т.д. N-11 слой — г^-1 > г > г^ = й, то компоненты Ф(0) связаны с Ф(й) матрицей переноса
1
Т = Дехр^-(г. - г—)], (15)
з=м
где Ж. — определяющие матрицы отдельных слоёв.
Одним из способов расчёта матрицы переноса неоднородного слоя является аппроксимация последнего набором достаточно большого числа однородных слоёв и использование формулы (15).
Наиболее эффективным методом вычисления матричной экспоненты ехр(Жг), порядок которой больше четырёх, является метод симметрических многочленов в сочетании с масштабированием матрицы Ж г [4]. Согласно этому методу матрица переноса (14) упругих волн в кристалле, имеющая шестой порядок, может быть найдена по формулам:
Т = {ехр (^)Г (16)
ехр
(Щ=¿(
V т / ^—' \ т
4 / 1=0 4
ЖЛ11
' и
т
I те 1
1 + Рб-г+^ л&-1-д (6)
д=о 3=6 ]'
(17)
где т — целое число, называемое масштабирующим фактором, р. = (—1)3-1оз, ] = 1,..., 6, — элементарные симметрические многочлены
матрицы Шг/к, а
В (6) =
0,
1,
если ] если 3
0,1, 3, 4, 5,
(18)
рВ-1(6) + р2В-2(6) + ... + РбВ-б(6), еслиз > 6,
называются симметрическими многочленами шестого порядка.
Выбор значения масштабирующего фактора т определяется необходимой точностью вычислений. Рассмотрим алгоритм вычислений матрицы переноса по формулам (16) - (18) более подробно на примере кристалла кубической сингонии.
3.1. Определяющая матрица для кристалла кубической сингонии. Наличие у сплошной среды элементов симметрии приводит к уменьшению числа независимых компонент тензоров (¿дг). Так матрица коэффициентов упругой податливости Бп для кристалла кубической сингонии имеет лишь три различных ненулевых элемента [5]:
(19)
Соотношения (13) и уравнение (12) принимают в данном случае следу-
(¿11 ¿12 ¿12 0 0 0
¿12 ¿11 ¿12 0 0 0
¿12 ¿12 ¿11 0 0 0
0 0 0 ¿44 0 0
0 0 0 0 ¿44 0
0 0 0 0 Б44/
Фт =
-ъкохБифз ¿1206 , 08 = гкоа ¿1203 ¿1206
С2 С2 ¿12 ¿11 ¿12 + ¿11 =2 ¿12 С2 - ¿11 ¿12 + ¿11
Ф1 0 ^12 0 0 0 0 01
Ф2 ^21 0 0 0 0 0 02
й Фз 0 0 0 W34 W35 0 03
йг Ф4 0 0 W43 0 0 W46 04
05 0 0 W53 0 0 W5б 05
Фб 0 0 0 Wб4 Wб5 0 0б
09 = гкох ф1
¿44
(20)
Ненулевыми элементами матрицы Ш = II 1з здесь являются:
1ц — ^44,
к2
121 = -рш2 + -0х ,
Б44
134 = -гкох = 71,
135 = 112 = 72, Б12
146
Б11 —
2Б12 _
= 74
153 = -РШЛ
Б12 +2 Б; кПх Б
11
0х°И _
= 75
22 Б12 Б11
156 = 143,
2
164 = -рш2 = 76,
143
гкл
0х
Б12 + Б;
= 73 , 165 = 134.
11
(21)
Новые обозначения 71 ,^2, ■ ■ ■ введены в этих формулах для упрощения итоговой записи матрицы переноса.
3.2. Симметрические многочлены матрицы Шг. Элементарные симметрические многочлены а3 ,■■■, а6 матрицы Ш равны суммам главных миноров матрицы Шг соответствующего порядка. Один из способов вычисления коэффициентов Рз = (—1)3 1аз состоит в использовании формул Ньютона
Рз =-(«3 - Р18]-1 - Р28]-2 - ■ ■■ - Рз-в), 3 = 1, 2, ■ ■ ■, 6, (22)
где вз = Ьг(Шг)3 — след матрпцы (Шг)3. Так, для определяющей матрицы (21) кристалла кубической сингонии
Р1 Р2 Р4
0, Р3 = 0, Р5 = 0,
(121112 + 2134143 + 112153 + 146164^, [124146153 + 123112164 - 1^31^4 - 112153164146 - 112121(2143134 + 164146 + 112153)]г4,
(23)
Р6 = [(13421432 - 1342146153 - 1121432164 + 146153112164)1211ц]г6■
Эти равенства имеют простой физический смысл. Корнями характеристического уравнения матрицы Ш являются нормальные компоненты волновых векторов кь, к2г, ■ ■ ■, к6г шести волн упругих деформаций, которые могут распространяться в кристалле. Эти компоненты в рассматриваемом случае попарно равны с точностью до знака: к1г = -к4г, к2г = -к5г, к3г = -к6г (здесь волнам, имеющим положительные проекции волновых векторов присвоены номера 1,2,3, а отрицательные — 4,5,6). Коэффициенты характеристического уравнения аз являются элементарными симметрическими многочленами относительно корней это-
Рз 3
нулю, а с четными —
Р2 = -02 = -(к1хг)2 - (к2хг)2 - ^г)2,
Р4 = -04 = -(к1гг)2(к2хг)2 - (кьг)2^г)2 - (к^г)2^г)2, Рб = -06 = - (к^ г)2(к2г г)2(к3г г)2.
Будем считать, что волны пронумерованы в порядке убывания модуля проекции волнового вектора на ось г, т.е. \к1г\ > \к2х\ > \к3х\. Следова-
|Р2| < 3|(кьг)2\, \р4\ < 3\(к^г)4\, \рб\ < \(киг)б\, (24)
и, если \(к1гг)2\ < 1, то
\Рб\ < р4\ < \Р2\ < 3\(киг)2\. (25)
Соотношения (23) позволяют переписать рекуррентные формулы (18) применительно к симметрическим многочленам матрицы (21)
В (6) = Р2В-2 (6) + Р4В-4(6) + РбВ_б(6). (26)
Отсюда и из (18) следует, что все многочлены (6) с чётными индексами равны нулю. Оценим модули ненулевых многочленов (с индексами 3 = 2д + 1 д = 2, 3,...). Из формулы (26) и неравенств (24) следуют соотношения
Ь (6)\ < \Р2\|В^_2(6)\ + \Р4\\В_4(6)\ + \Рб\\В_б(6)\ (27)
< \Р2 \(Р\Щ_А(6)\ + |Р4|\В_б(6)| + \Рб\\В^_8(6)\)
+ \Р4\\В,-_4(6)\ + \Рб|\В_б(6)|
|)
■(\Р2\\В^_4(6)\ + \Р4\\В_б(6)\ + \Рб\\В^_8(6)^ . (28)
\Р2\
Используя рекуррентный характер неравенств (27) и (28) получаем следующую оценку величин модулей симметрических многочленов матрицы Шг:
|«) ^^^ \ . 3 > 3. (29)
3.3. Величина масштабирующего фактора m в формулах (16) и (17) выбирается из условия быстрой сходимости ряда Ej=6 1Bj-1-g(6), такой, чтобы аппроксимация этого ряда суммой первых N+1 слагаемых
j 6+N
£ Bj--g(6) ^ Bj--g(6), (30)
j=6 j' j=6 j'
обеспечивала достаточную точность вычислений.
Для матрицы Wz/m по аналогии с неравенствами (24) получаем
Ы< 3|(fcizz/m)X |Р4|< 3|(k\zz/m)41, Ы < \(klzz/m)%
и, согласно (30), |B2j+1(6)| < (4|k1z\z/m)j-2 , j > 3. Если в качестве m взять наименьшее целое, удовлетворяющее условию m > 40|k1z z,
Wz/m
неравенствам
|Br(6)| < 10-1, |Bg(6)| < 10-2, |Bii(6)| < 10-3,... (31) При выполнении этих условий суммирование в (30) до N = 3 обеспечи-
вает вычисление
в формуле (17) с ОТ НО-6
¿'Е P6-i+gE9=61 Bj-i-g (6)
g=0
10
Пример. Плотность кристалла GaAs р = 5.32 • 103кг/м2, а коэффициенты упругой податливости имеют следующие значения: Sn = 11.72 • 10-12м2/Н, 5*12 = -3.65 • 10-12м2/Н, S44 = 16.83 • 10-12м2/Н. При нормальном падении ультразвука частотой v =15 • 106Гц на пластинку этого материала толщиной z = 10-3м, |k1zz < 2nv^JpS44z ~ 28.3. Следовательно, для выполнения соотношении (31) необходимо выбрать m > 1131.
3.4. Вычисление экспоненты exp(Wz/m). После выбора мас-
m не составляет боль-
( Wz/m) 2
(Wz/m)3, (Wz/m)4 и (Wz/m)5 и величины pj для матрицы Wz/m. Последняя задача может быть решена параллельно первой, если использовать формулы (23).
3.5. Повторное квадрирование. Заключительным этапом нахождения матрицы переноса по методу симметрических многочленов с масштабированием является нахождение матрицы (16). Наименьшую погрешность вычислений, обусловленных округлением чисел при возведении матрицы в m-ю степень (если m > n + 1), в сравнении с другими
известными подходами, обеспечивает метод симметрических многочле-
нов [4]. Это справедливо и для случая т = 2-7', где 3 -примере, рассмотренном выше, целесообразно т
целое. Так, в 2048 = 211.
4. Заключение
Метод симметрических многочленов (МСМ), в отличие от методов вычисления матричной экспоненты, основанных на матричных преобразованиях подобия или формулах Лагранжа, Ньютона, Бэкера не требует н ахож д е н и я собственных значений матрицы Ш, которые для матрицы шестого порядка общего вида могут быть определены только численными методами. Это обеспечивает большую надежность и точность вычислений матрицы переноса по формулам (16)—(18), в сравнении с другими подходами. При этом для матриц переноса 2 - 4-го порядков МСМ, наравне, например, с формулой Лагранжа, позволяет находить аналитические формулы. Примером является матрица переноса упругих деформаций в кристаллах кубической сингонии. Для них, в рассматриваемой здесь геометрии рассеяния, определяющая матрица Ш имеет блочно-диагональную структуру: Ш = №зн 0
Шяи = \\° т12\\,Шр^у =
0 Wp-sv
0 ад34 адзб 0
•Ю43 0 0 -Ю46
-ш5з 0 0 т56 0 т64 т65 0
Такая структура определяющей матрицы, показывает, что для выбранной геометрии рассеяния плоской волны на пластине из кристалла кубической сингонии процессы распространения горизонтальной волны сдвига (волна б^-тииа) и волны Р- ¿У-типа друг с другом не связаны. Система уравнений (20) распадается на два матричных уравнения
dФ
SH
dz
Wsh Фsн,
dФр-sv dz
Wp -SV Фр-SV,
где компонентами вектор-столбцов Ф^ и Фр-SV соответственно являются ф1,ф2 И. фз,ф4,фб ,Фб-
Элементарные симметрические многочлены матрицы второго порядка WSH имеют следующие значения: a1 = 0, a2 = рш2S44 — = k2 cos2 0. Здесь k1 = ш^/pS44 и k1 sin9 = k0 sin90. Следуя [4], находим матрицу переноса упругих деформаций горизонтальной волной сдвига в кристаллическом слое:
Tsh = exp(WsH z)
cos(k1z cos 9)
S44 sin(k1z cos 9)
k1 cos 9
S
44
sin(k1z cos 9)
k1 cos 9 cos(k1z cos 9)
Элементарные симметрические многочлены матрицы WP —SV равны: Ох = 0, 02 = -2W43W34 - W35W53 - W46W64, О3 = 0 О4 = W342W432 -
w64w35w432 - w53w342w46 + w53w64w35w46. Полагая в формулах (16) и (17) m = 1 и выполняя суммирование в правой части формулы (16) по методу, представленному в статье [6], находим следующие значения элементов tij, i,j = 1, 2, 3, 4, матрицы переноса TP-SV = exp(WP-SVz) волны P-SV типа: tn = t33 = So + (7173+7275)#2, ti2 = 71S1 + [(Ti73 + 7275) 7i + (7174 + Y2Y3) Y6]S^ ti3 = 72S1 + [(7173 + Y2Y5) Y2 + (YiY4 + Y2Y3) Yi]S3, ti4 = t23 = (7174+7273) S2, t21 = 73Si + [(YI73 + 7476 ) 73 + (7174 + 7273 ) 75]S3, t22 = t44 = So + (7i73 +7476)S2, t24 = Y4S1 + [(7173 + 7476)Y4 + (YIY4 + 7273)73] S^ t3i = 75S1 + [(7571 + 7376)73 + (Yi73 + 7275)75]S3, t32 =
t41 = (7571+7376)S^ t34 = 73S1 + [(7571 + 7376)74 + (7173 + 7275)73]S3, t42 = 76Si + [(7571 + 7376)71 + (7173 + 7476)Y6]S3, t43 = 7iSi + [(7571 + 7376)72 + (7173 + 7476 )Yi]S3-
Здесь 71 , . . . , 76 (21);
S = klz cos (zk3z) - k1z cos (zk2z) S = kfz sin (z^3z) - sin (zk2z)
k2z - k3z k2z k3z (k2z - k3z)
S = cos (zk3z) - cos (zk2z) S = k2z sin (zk3z) - k'3z sin (zk2z) ;
k2z - k3z k2z k3z (k2z - k3z)
1/2
>'2,3z
7173+7275+7476 ^ у/ (7275 - 7476)
± V 2 \ 4 6--(7175 +7376)(Y273+7471)
2 V 4
это проекции волновых векторов Р — БУ волн на ось 0г. Литература
1. М о л о т к о в Л. А. Матричный метод в теории распространения волн в слоистых упругих и жидких средах. Л.: Наука. 1984. 201 с.
2. БреховскихЛ.М., ГодинО.А. Акустика слоистых сред. М.: Наука. 1989. 416 с.
3. Красильников В. А., К р ы лов В. В. Введение в физическую акустику. М.: Наука. 1984. 400 с.
'-А. £) 6 Л Я 6 В ю.н. к вычислению функций матриц // Математические заметки, 2013. Т. 94, Вып. 2, С. 175-182.
2
5. Сиротин Ю.И., Шаскольская М.П. Основы кристаллофизики. М.: Наука. 1979. 640 с.
6. Беляев Ю. Н. Симметрические многочлены в расчётах матричной экспоненты // Вестник СыктГУ, ('ср. 1 Математика, механика, информатика, 2012. Вып. 16, С. 28-41.
Summary
Belyayev Yu.N., Popov S.A. Transfer matrix of elastic deformations in crystals
Differential equations of elastic waves in crystals are solved using sixth-order symmetric polynomials and scaling method. Influence of layer thickness and frequency of the wave by the scaling factor is investigated. Analytic solution describing the transfer of elastic stresses in the crystalline layer of the cubic system is obtained.
Keywords: layered media, waves, matrix, symmetric polynomials, truncation error, scaling.
Сыктывкарский университет Поступила 12.10.2013
