Замечание о вычислении скорости волны Рэлея и производной определителя Рэлея в упругих средах

Поверхностные волны или волны Рэлея [1] давно и широко применяются в науке и технике, например, для определения характеристик их излучателей и изучения физических процессов, происходящих при возбуждении и распространении колебаний [2, 3], фазовых переходов [4], изучения свойств веществ и состояния их поверхностей [5], в дефектоскопии и оценке остаточного ресурса [6–8], для передачи и обработки информации [9, 10], в геофизике и сейсмологии [11].

Задача вычисления скорости рэлеевской волны, сводящаяся к решению уравнения третьей степени, раньше решалась численно [2] или с помощью простых приближенных формул [6, 12– 15]. Поскольку к настоящему времени математики уже предложили аналитические методы решения уравнений невысоких степеней, например метод Лагранжа [16] или формулу Кардано [1719], то эти методы учитываются при исследовании выражений в программах компьютерной алгебры. Это позволяет записать и упростить аналитические выражения для скорости поверхностной волны и для корней соответствующего уравнения [20, 21] (см. также источники 4, 6–13 в [20]).

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

Согласно [22, с. 136], уравнение для определения скорости волны Рэлея

г

^ 6

- 8 ^ 4 + 8 ^ ²

V

2 А

³ - ² г cl 7

-

(

16 1

V

-

„ 2 А cl 7

= 0,

где § = to/ctk = er|et (0 < р < 1), to - циклическая частота колебаний, k = to/c - волновое число, cr,ct,"i - скорость поверхностных, поперечных и продольных волн соответственно.

Сделаем замену x = с2 = (k/kt )2, введем обозначение и2 = (ct I ci )2 и получим уравнение x3 - 8 x2 + 8 x (3 - 2 и 2)-16 (1 - и 2 ) = 0,                            (2)

имеющее в интервале x е [0,1) единственный действительный корень a , который и определяет скорость волны Рэлея:

c_r = c_t4a .                                           (3)

Введем обозначения

Физика

T = 4 (( - 64 u ² + 107) u ² - 62) u ² + 11,

U = 3/17 - 45 u ² - 3^ T                                (5)

и запишем решение уравнения (2), полученное в программах Maxima [23] и WolframAlfa [24] в виде

a = -

2       2(6 u ² - 1)

3        U

Это выражение вместе с (4) и (5) является решением поставленной задачи. Значение корня определяется только отношением квадратов скоростей объемных волн или коэффициентом Пуассона а , связанным с u ² известным соотношением u ² = c 2 / c ² = (1 - 2 а )/2(1 - а ) (табл. 1).

Таблица 1

Зависимость а ( а )

а	u 2	a , точное решение (6)	a , численный метод с )
–1,0	0,7500000000000000	0,4745724391564827	0,474572439156483
–0,9	0,7368421052631579	0,4960417626756930	0,4960417626756933
–0,8	0,7222222222222222	0,5191753282850295	0,5191753282850295
–0,7	0,7058823529411765	0,5440779615104171	0,5440779615104173
–0,6	0,6875000000000000	0,5708262701628667	0,570826270162867
–0,5	0,6666666666666667	0,5994463782101697	0,5994463782101699
–0,4	0,6428571428571429	0,6298836881418670	0,6298836881418671
–0,3	0,6153846153846159	0,6619660579475539	0,661966057947554
–0,2	0,5833333333333333	0,6953666629760182	0,6953666629760183
–0,1	0,5454545454545455	0,7295801516555792	0,7295801516555793
0,0	0,5000000000000000	0,7639320225002102 a )	0,7639320225002103
0,1	0,4444444444444444	0,7976383362116029	0,797638336211603
0,2	0,3750000000000000	0,8299135133739662	0,8299135133739663
0,3	0,2857142857142857	0,8600943341185433	0,8600943341185434
0,4	0,1666666666666667	0,8877322341853701 b )	0,88773223418537
0,5	0,0000000000000000	0,9126219746158490	0,9126219746158474

a ) 3 - 75, b ) а = 2(4 - 3^9)/3 [21], с ) – значения получены с помощью [23] при использовании функции find_root , пример кода: s : - 1.0; u 2:0.5 *(1 - 2* s )/(1 - s ); find _ root ((( x - 8)* x - 16* u 2 + 24)* x + 16* u 2 - 16, x ,0.47,0.92);

В научных работах, например в [20], есть подобные результаты, но они записаны в другой форме – в виде подробного алгоритма действий, которые к нему приводят, а не в виде конечной формулы (6). В [20] приведена таблица корней уравнения (2) для различных значений коэффициента Пуассона ( - 1 <  а < 0,5) и корней, полученных численным методом (см. табл. 1 из [20]). Из рассмотрения этой таблицы следует, что необходимо порядка 10 итераций для вычисления корня с абсолютной погрешностью, меньшей 10 ^{- 9}. Очевидно, что точное выражение существенно проще в использовании, чем применение численных методов, которые требуют многократного вычисления исходной функции (2).

Согласно [22, с. 137], u ² изменяется для различных веществ от 0 до 12 ( u от 0 до 1/V2). Решение (6) допускает подстановку любого значения u ² из указанного интервала и остается действительным. Вычисление не вызывает никаких трудностей в интервале 1/6 <  u ²< 0,3215 с использованием действительных чисел. При u ²< 1/6 ( а >  0,4), подкоренное выражение в (5) отрицательно и это требует аккуратного вычисление корня U = - 3 - 17 + 45 u ² + 3*j3T при работе с такими веществами, например, как свинец и золото. В самой точке u ² = 1/6 ( а = 0,4) и ее малой окрестности при вычислении функции может накопиться значительная вычислительная погрешность, так как для получения значения знаменателя в формуле (6) требуется существенно больший объем вычислений, чем для числителя. Неопределенность типа 0/0 легко устраняется разложением функции U ( u ²) в ряд Тейлора в окрестности этой точки

Гуревич С.Ю., Голубев Е.В.

U ( и2) «

£ ( u 2

1)   360 ( 21

- +--.,,— u —+

6 J  19319 I что дает из (6) значение корня a = 2(4 - 319)/3 ~ 0,887732234185370088.... Автором [21] получен аналогичный результат из (2) при u2 = 1/6. При u2 > 0,3215 (u = cd cl > 0 ,567) T (4) становится мнимым и в расчетах необходимо использовать переменные комплексного типа (например, для таких веществ как цинк, германий, бериллий), для получения действительного значения a необходимо выделить реальную часть, мнимая часть сравнима с вычислительной погрешностью.

В работе [21] с помощью формулы Кардано и программы MAPLE получено выражение (см. формулу после 2.14), определяющее действительный корень через коэффициент Пуассона, которое в заключительной части работы автор привел к (6), записанному не в упрощенном виде. Автором отмечено, что функция имеет разрыв в точке и ² = 1/6 ( ст = 2/5), что, как показано, не имеет места. В [21] также получено, что a = 2(4 - 3 19)/3, при ст = 2/5 = 0,4. В [25] приведены точные

15 2 5 3 77 20 55 114 значения друг их корней для ^{ст равн}ы1х у , ,у , 123,28,365,69,136,235 ^.

Результаты вычисления на основе выражения (6) можно сравнить с экспериментальными данными для величины скорости рэлеевской волны в различных материалах и признать удовлетворительным соответствие расчетных и табличных данных. Так, в работе [26] приведены данные измерений скорости: в алюминии (А-1) – 2990 м/с, в железе (АРМКО) – 2912 м/с. По данным [7] (см. приложение, табл. П2, 10³ м/с): свинец – 0,63; золото – 1,12; платина – 1,57; серебро – 1,48; висмут – 1,03; латунь – 1,95; вольфрам – 2,65; медь – 3,52 (указана c_t = 3,72); алюминий – 2,80; олово – 1,56; никель – 2,64; кадмий – 1,4; железо – 3,0; цинк – 2,22; бериллий – 7,87.

При решении задачи возбуждения и распространения акустических волн в сплошных упругих средах при использовании модели полупространства и методов интегральных преобразований для нахождения решения в аналитическом виде конечные квадратурные формулы для акустического поля (поля векторов деформаций и напряжений) содержат в знаменателе выражение [12, стр. 7] (I.6) (также см. [2]):

R ( к ) = ( к ² + 5 ² )² - 4 к ² qs = (2 к ² - к ₂² )² - 4 к ² 7 к ² - к 2 7 к ² - к 2 , (7) где к = < у/c - волновое число; q = 7 к ² - к2 , к = ^c; s = 7 к ² - к 2 , k_t = to!c_t . Согласно [27],

R(к) определяет четыре точки ветвления подынтегральной функции к = ±к1, ±к2 и три полюса к = 0, ±кг, где кг = <у/cr - волновое число волны Рэлея. Вычет в к = кг определяет вклад волны

Рэлея в акустическое поле. Для его определения необходимо вычислить производную

R'(кг) = dR dk

к = к_г

22   22

3 к          I О О 6 О О       1 XI к     к*        к XI - к]

8 к (2 к ² - к /) - 8 Цк ² - к^ 7 к ² - к 2 - 4 к ³        t - 4 к ^{3        1}

4 к ² - к2      к к ² - к 2

к = к_г

2 к 8 - 8 к ² к 6 + 16 к ⁶( к 2 - к2) к (2 к ² - к /)²

2[ a ³ - 4 a ² + 8(1 - и ²)] ₃

7 a ( a - 2)²        t

/rs 7 2 7 2 \ 2 А 7 2 / 7 2 7 2 1 7 2 7 2 7 7 / Г" 7 7 ГХ . .

где учтено, что (2 к_г - к₂) = 4 к_г к к- - к ₁ ^к- - к ₂ , к_г = k_t|^a , к_г = uk_t . Это выражение, записанное в аналитическом виде, пропорционально частоте в третьей степени. Безразмерный коэффициент перед k ³ (приведен для некоторых веществ в табл. 2), от которого зависит амплитуда волны Рэлея, определяется только отношением скоростей распространения объемных упругих волн или коэффициентом Пуассона.

Аналогичные результаты представлены в [2], где приведены графики расчетных зависимостей величин скорости волны Рэлея и отношения производной определителя и k ³ от коэффициента Пуассона (см. рис. 2 и 3 в [2]).

Таблица 2

Физика

Скорости звука для некоторых веществ ( c , c – данные [28], столбцы 4–7 – расчет по (6), (3) и (8))

	cl	ct	2 = 2/2 u = ctcl	a	cr	1 dR          1 k 3 dk k = kr
1	2	3	4	5	6	7
Свинец	2160	700	0,105024	0,898355	663,47	–8,09648
Золото	3240	1200	0,137174	0,893063	1134,02	–7,64224
Платина	3960	1670	0,177845	0,885579	1571,56	–7,07566
Серебро	3600	1590	0,195069	0,882107	1493,34	–6,83877
Висмут (кристалл)	2140	960	0,20124	0,880815	900,98	–6,75439
Нейзильбер	4760	2160	0,205918	0,879817	2026,05	–6,69059
Латунь	4430	2123	0,229664	0,874507	1985,33	–6,36919
Вольфрам	5460	2620	0,230259	0,874369	2449,9	–6,36118
Медь	4700	2260	0,231218	0,874145	2113	–6,3483
Алюминий	6260	3080	0,242077	0,871556	2875,4	–6,20289
Олово	3320	1670	0,253021	0,868846	1556,64	–6,05731
Константен	5240	2640	0,253831	0,868641	2460,5	–6,04658
Висмут	2180	1100	0,254608	0,868444	1025,09	–6,03629
Никель	5630	2960	0,276418	0,862689	2749,28	–5,7496
Чугун	4500	2400	0,284444	0,860454	2226,26	–5,6452
Свинец (кристалл)	2350	1266	0,290223	0,858804	1173,22	–5,57043
Кадмий	2780	1500	0,291134	0,858541	1389,86	–5,55867
Олово (кристалл)	3480	1900	0,298091	0,8565	1758,4	–5,46916
Железо	5850	3230	0,304855	0,854465	2985,72	–5,3826
Цинк	4170	2410	0,334012	0,845069	2215,46	–5,01527
Германий (кристалл)	5390	3540	0,431349	0,80447	3175,11	–3,87086
Бериллий	12660	8900	0,494211	0,767838	7798,7	–3,2182

Выводы

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