Одномерная задача о нестационарной связанной упругой диффузии для слоя

using the finite element method Am. J. Orthod. Dentofac. Orthop., 1989, vol. 95, pp. 200-207.

29. Tindlund R. S., Rygh P., Boe O. E. Intercanine widening and sagittal effect of maxillary transverse expansion in patients with cleft lip and palate during the deciduous and mixed dentitions. Cleft Palate Craniofac. J., 1993, vol. 30, pp. 195-207.

30. Yu H. S., Baik H. S., Sung S. J., Kim K. D., Cho Y. S. Three-dimensional finite-element analysis of maxillary protraction with and without rapid palatal expansion. Eur. J. Orthod., 2007, vol. 29, pp. 118-125.

31. Landes C. A., Laudermann K., Petruchin O., Mack M. G., Kopp S., Ludwig B., Sader R. A., Seitz O.

Comparison of bipartite versus tripartite osteotomy for maxillary transversal expansion using 3-dimensional preoperative and postexpansion computed tomography data. J. Oral. Maxillofac. Surg., 2009, vol. 67, pp. 22872301.

32. Zimring J. F., Isaacson R. J. Forces produced by rapid maxillary expansion III. Forces present during retention. Angle Orthod., 1965, vol. 35, pp. 178-186.

33. Wood S. A., Strait D. S., Dumont E. R., Ross C. F., Grosse I. R. The effects of modeling simplifications on craniofacial finite element models: The alveoli (tooth sockets) and periodontal ligaments. J. Biomech., 2011, vol. 44, pp. 1831-1838.

УДК 539.3

ОДНОМЕРНАЯ ЗАДАЧА О НЕСТАЦИОНАРНОЙ СВЯЗАННОЙ УПРУГОЙ ДИФФУЗИИ ДЛЯ СЛОЯ

А. Р. Гачкевич1, А. В. Земсков2, Д. В. Тарлаковский3

1 Доктор физико-математических наук, профессор, Институт прикладных проблем механики и математики им. Я. С. Под-стригача НАН Украины, Львов, dept13@iapmm.lviv.ua

2Кандидат физико-математических наук, доцент кафедры математического моделирования, Московский авиационный институт (национальный исследовательский университет), azemskov1975@mail.ru

3Доктор физико-математических наук, заведующий кафедрой сопротивления материалов и деталей машин, профессор, Московский авиационный институт (национальный исследовательский университет), tdvhome@mail.ru

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

Ключевые слова: упругая диффузия, нестационарные задачи, ряды Фурье, преобразование Лапласа.

ВВЕДЕНИЕ

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

Соотношения между параметрами, характеризующими такую систему, строятся с помощью метода потенциалов, разработанного Дж. Гиббсом. При этом в качестве потенциалов выступают внутренняя энергия, свободная энергия и химический потенциал (или удельный потенциал Гиббса). Основополагающими законами, позволяющими выразить термодинамические потенциалы через термодинамические параметры, являются первый и второй законы термодинамики. Для вывода уравнений, описывающих физико-механические процессы в деформируемой среде, применяются уравнения движения упругого тела и законы массопереноса [1-4]. При этом используются следующие предположения:

1) изменения концентрации, сопутствующие деформированию, малы;

2) возрастание концентрации не приводит к существенным изменениям упругих и диффузионных характеристик материала, поэтому эти константы рассматриваются как величины, не зависящие от концентраций;

3) деформации малы.

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

1. ПОСТАНОВКА ЗАДАЧИ

Рассматривается однородный плоский слой, ограниченный плоскостями х = 0 и х = X. Предполагается, что физико-механические процессы в среде являются одномерными и описываются геометрически линейной моделью упругой диффузии без учёта температурных эффектов [1,4]. Соответствующие уравнения имеют следующий вид (компоненты векторов и тензоров указаны в прямоугольной декартовой системе координат 0X1 х2х3, для которой х3 = х):

„ д2 и д 2 и дп д2п дп . д3и

03333^-2 = Р^ГГ + а33о-' = + '

дх2 д£2 дх дх2 д£ дх3

а диффузионный поток в направлении оси 0х задаётся соотношением

д2

3 = Лдх2 — ^ах' Л = —таТ"'

где £ — время, и — перемещение в направлении оси 0х; п = п — п0 — приращение концентрации; п0 — начальная концентрация вещества, п — текущая концентрация вещества; 03333 — соответствующая компонента тензора упругих постоянных; р — плотность среды; а33 — коэффициент, определяемый типом кристаллической решётки таким образом, что величина а33е33п есть относительное объёмное изменение, выражающее зависимость объёма вещества от состава (концентрации исходного вещества); е^ — компоненты тензора деформации; Я33 — коэффициент самодиффузии; д — термодинамический множитель Даркена; Я — универсальная газовая постоянная; Т0 — температура.

Полагается, что на поверхностях слоя заданы перемещения и диффузионные потоки:

и|х=о = /1

J |х=0 = /3

а начальные условия однородные:

4=0 =

ди

ж

= пи = 0.

¿=0

Используются следующие безразмерные величины (при одинаковом начертании они обозначены символом «~», который в дальнейшем изложении опускается):

х7 = х и 7 = X'

Л. П0' Я = Я33д сХ

с£ 2 03333 П0 «33 т = —, с =-, а = ^

3333

X' Р ' 0

'/ п -оо^ т Л 7 3 г, ч /к

п = —, Я = -~г' Л = ^—' з7=-' 7 (т ) =-•

п0 сх схп0 сп0 сп0

Тогда начально-краевая задача (1)-(5) записывается так (здесь и далее штрихи обозначают производные по координате, а точки — производные по безразмерному времени т):

и'' = и + ап', пп" = п + Ли''',

3 = Ли'' — Яп',

и|х=0 = /1

2. ИНТЕГРАЛЬНОЕ ПРЕДСТАВЛЕНИЕ РЕШЕНИЯ

Пусть 0ик = и, 0пк = п — функции Грина задачи (6)-(9), а именно решения четырех задач (к = 1,4 — их номер), включающих в себя уравнения (6), соотношение (7), начальные условия (9) и следующие граничные условия:

0ик х=0 = &1к6(т), 0„к |х=1 = $2к6(т),

где 6(т) — дельта-функция Дирака, — символ Кронекера.

Изв. Сарат. ун-та. Нов. сер. Сер. Математика. Механика. Информатика. 2013. Т. 13, вып. 4, ч. 1 Тогда решение задачи (6)-(9) имеет вид (звездочка обозначает свертку по времени):

4 4

и = ^) * /

k=1 k=1

Отметим, что из (6), (7) и (8) с использованием замены пространственной переменной y = 1 — x следует следующая связь функций Грина:

Gu2 (x,T) = Gui (1 — х,т), Gu4(x,T) = Gu3 (1 — Х,т), Gn2(x,T ) = — Gni (1 — x,T ), G^4 (x,T ) = — Gn3(1 — х,т ).

Отсюда вытекает, что формулы (11) можно записать так:

2

u = [Gu,2k-i(x,T) * /2fc-l(T) + Gu,2k-1 (1 — x,T) * /2k (т)] ,

— x,T ) *

1

п = X] [^п,2А-1 (х, т) * /2Й-1 (т) - Сп,2Л-1 (1 - X, т) * /2^ (т)] .

Таким образом, для определения перемещения и приращения концентрации достаточно найти только функции Guk и Gnk при к = 1,3.

3. ИЗОБРАЖЕНИЯ ФУНКЦИЙ ГРИНА

Рассмотрим задачу (6)-(9) при /2 = /4 = 0, и применим к ней преобразование Лапласа по времени [4, 5] (^ — параметр преобразования, индекс Ь обозначает трансформанту):

52- '' + ап^' = 0, - -п^" + ''' = 0,

иЬ|х=0 = #, иь|х=1 = 0,

Для решения краевой задачи (13), (14) будем использовать разложения искомых функций в ряды собственным функциям. Для этого редуцируем краевую задачу к однородным краевым условиям. С этой целью положим

= иь + р, пЬ = Нь +

где функции р и ф выбираются так, чтобы они удовлетворяли граничным условиям (14):

р

Тогда для функций иь и получаем следующую задачу:

52иь - иь" + аНь' = ^х, 5Нь - -Н+ Ли= F3,

Fi (x, s) = p*(x) [—s2/f (s) + / (s)] , F3(x, s) = [1 — s^* (x)] /f(s),

D'

lx=0 lx=1

AUL" — DHL M =0, (AUL// — DHL '

= 0.

x=0 V / x=1

Ее решение представляем в виде рядов по собственным функциям:

те те

UL

n=1 n=0

При этом записываем в виде аналогичных рядов функции р*, ф* и правые части уравнений в (17):

те

_ \ Л , .L \ _ ,/,*

р* (x) = ^^ pf sin Anx, (x) = ^2 cos Anx,

n=1 n=0

те

F.(x, s) = F1n(s) sin Anx, F3(x, s) = F3n(s) cos Anx,

n=1 n=0

где

=/ Я^*(х) ¿х = —1, ^30(5)=/ ^3(х,5) ¿х = (1 — /Ь(*), и0 3 ./0

1 о Г1

2

= 2 I Р^*(х)сов Лпх^х = ^, = 2 I ^*(Х)Б1П Лпх^х =

0

0

2

Лп

Рщ(5) = 2 Р^х^ВШ Лпх ¿х = [ — 52/Ь(5) + /(«)] ^п,

Я'

^3п(й) = 2 / Р3(х, 5) СОБ Апх ¿х = —п ^ 1.

0

Подставляя ряды (18) и (19) в (17), получаем следующие алгебраические уравнения для коэффициентов и ЯЬ: - при п = 0

^ = ^30,

- при п ^ 1

(52 + ЛП)ЦЬ — аЛп ЯЬ = Рщ, —ЛАП ЦЬ + (5 + -ЛП)ЯЬ = Р„3.

Решение последней системы имеет вид

5 + -Лп

аЛП

ГтЬ _ " 1 "га -цр | тр ТгЬ _ ^"П ТР I

= ---Г1п + —-Гп3, Яп = —— +

Ь ЛЛП тп , 52 + Лп

р

3п

р

3п

п р - 1п ' р Р3п р

Рп3,

3п

Р3п (5) = (в2 + ЛУ (5 + — аЛА:

Учитывая равенства (20), для изображений и Я^ получаем следующий результат:

ЦЬ = —^ /Ь +

Ь , 2Лп(5 + (- — аЛ)Лп) „, 2аЛ Р + Р

р3п р3п

тть , 2ЛЛп52 2(52 + лп) ЯЬ = —^п /Ь--^Т*— /1Ь + -— /Ь.

/Ь + -^/Ь,

р3п

р

3п

р

3п

После подстановки (16) и (18) в (15) с учетом формул (19), (21) и (22) получаем решение задачи (13), (14):

иЬ =

2 Е Ап [(5 + (- — аЛ)Лу /Ь + а/3

б1п Лп х

п=1

р3

3п

пЬ = 1 — ^ [лап 52 /Ь — (52 + лп)/Ь

соб Лпх

п) /

п=1

р

3п

Сравнивая его с (11) с учетом равенства (20) и свойств преобразования Лапласа, находим изображения функций и при к = 1,3:

где

(х,5) = (5)й1п Лпx' 0^3(х,5) = 0Ь3п(Ф1П Лпх

п=1 п=1

те те

(х, 5) = ^2 (5) СОБ Лпх, 0^3 (х,5) = С^в) + ^ (в) СОБ Л„х,

п=1

, 1 п (5) -

п=1

П30

(5) =

и1п( 1

5

2

Р3п

0п3п(5) =

2ЛЛП 52

р3

3п

Р3п

2аЛп Р

п > 1.

3п

Ь

4. ОРИГИНАЛЫ ФУНКЦИЙ ГРИНА

Наиболее просто вычисляется оригинал функции GL30 [5]:

Gnзо

где Н(т) — функция Хевисайда.

Остальные функции в (23) являются рациональными функциями параметра 5. Их оригиналы вычисляются достаточно просто с использованием разложений на элементарные дроби. Для этого необходим анализ нулей их знаменателя, т. е. корней кубического уравнения:

Рзп

Поскольку для большинства реальных сред имеют место неравенство [1-3]

0 < аЛ < - < 1,

то анализ проведем асимптотически, раскладывая корни »1п, »2п, »3п этого уравнения в степенные ряды по малому параметру е = аЛ (ограничиваемся первыми двумя членами):

»¡¡п = 4п + е4П + к = 1, 2, 3.

Подставляя эти равенства в (25) и приравнивая коэффициенты при нулевой и первой степенях

малого параметра, получаем следующие уравнения относительно »¡п и »¡п:

(»¡Г + А?,) (»¡2 + -АП) = 0, (3.<Г + 20А» .!0) + ап) »¡1 - А

Их решения имеют вид

Л0) = а »(0) = 7А »(°) = ПА2 »(1) = АП ро (ч) = 3»2 + 2ПА2 » + А 2

'1п = гА™, »2п = -гАп, »3п = - — Ап, »¡п = ~-/ (0)\ , Р2п(5) = 3» + 2—Ап» + Ап

Р2п (»¡п

Отсюда следует, что уравнение (25) имеет два комплексно сопряженных корня »1п и »2п, один действительный корень »3п < 0. Поскольку

Р2п (»10Л = -3АП + 2г-АП + АП = -2АП + 2ШАп = 2АП(-1 + ШАп),

то Яе »1п < 0, что в совокупности с отрицательностью »3п соответствует затуханию решений.

Теперь можем построить разложения рациональных функций, входящих в формулы (23), на элементарные дроби:

» + (— - аЛ)АП АП„в + А.12п А13п »2 А21„» + А^п , А>3п

+-, =—^^— +

Р3п »2 + ^П » »3п Р3п »2 + ^П » »3п

1 = А31п8 + А32п + А33п »2 + АП = А41„з+А42„ + А43п

Р3п »2 + 8 - »3п , Р3п *2+*П ^

где

»3п + (- - аЛ)АП . 4 - »3п(- - аЛ)АП . »3п + (- - аЛ)А2

Л _ "3п ' V" ""У п /I _ п ^- /"п Л __' у-1-^ ^ г'п

А11п = 2 I 2 , А12п = 2 I 2 , А13п = 2 I 2

2 2 2 2 2 2

+ »3п + »Зп + »3п

2 2 2

л _ л _ »3п^п л _ »3п

А21п = 2 +—~, А22п = 2 +—~, А23п = 2 +—~,

+ »3п + »3п + »3п

1 »3п 1

А31п =--^-2Г~, А32п =--^-^, А33п = _^-,

2 2 2 2 2 2

+ »3п + »2п + »3п

г2 - А 2 -2 - а 2 »2 + А 2

л __А п л _ „ — п А п л _ »3п + А п „ ^ л

А41п = 2 + 2 , А42п = »3 п 2 + 2 , А43п = 2 + 2 , п ^ 1

+ »3 п + »3п + »3 п

Подставляя теперь разложения (26) в (23) и учитывая (24), с использованием таблиц обратного преобразования Лапласа (см., например, [5]) получаем оригиналы функций влияния:

œ

Gui = 2 An [Aiin cos znT + Ai2nz-1 sin znT + Ai3nes3nr] ,

n=1

œ

Gni = -2Л Y^ ЛП [A2in cos znT + A22nz-i sin znT + A23neS3nr] ,

n=i

œ

Gu3 = 2a ^ An [A3in cos znt + А32Пz-1 sin znt + Аззпes3nr] ,

n=i

œ

Gn3 = H(t) + 2 [A4in cos znt + А42Пz-1 sin znt + ^43™es3nr] .

n=i

5. ПРИМЕРЫ

Положим, что материал слоя имеет следующие безразмерные характеристики: D = 0.01, a = 0.1, Л = 0.001. Рассмотрим два варианта граничных условий (8):

А. Обе поверхности слоя зафиксированы, на границе x = 1 диффузионный поток отсутствует, а при x = 0 он задан в виде функции Хевисайда, что соответствует следующим правым частям равенств (8):

fi (T ) = /2 (T ) = /4 (T ) = 0, /3 (T ) = H (T ).

Тогда, вычисляя интегралы в (12), получаем, что перемещения и концентрация определяются равенствами

u(x, T) = Gu3(x, T) * /3(t) =

œ

= 2a ^ An [A3inz-1 sin znt + А32Пz-2(1 - cos znt) - ^33ns—1 (1 - es3nr)] sin Anx,

n=1

n(x, t) = Gn3(x, t) * /3(t) =

œ

= T + 2 ^2 [A41n z-1 sin zn T + A42n z-2 (1 - cos zn T ) - A43n s-1 (1 - es3n r )] cos An x.

n=1

Результаты расчётов по этим формулам продемонстрированы на рис. 1, а, б, где представлены графики зависимостей перемещений и концентрации от времени в различных точках (сплошная линия соответствует x = 0.1, линия с точками — x = 0.5, пунктирная — x = 0.9).

20

15

0.2

0.1

0

10

10 20 30 40

а

4 6

б

Рис. 1. Зависимость перемещений (а) и приращения концентрации (б) от времени. Вариант А

0

2

8

0

Б. Положим, что /1 (т) = /2(т) = /4(т) = 0, /3(т) = в_т. Тогда для искомых функций получаем следующие формулы:

ГО . Л

и(х, т) = -2а £ АпА33п (в_т - в53«т) +

п=1

+ »3п

п=1

+2а Ап [(^31п - А32п) сов + (А31п-п + ^п^ эт -пт + (А32п - А31п) в т]

ГО

П (х,т ) = 1 - в_т - 2^2 А43п (<

_т1 бШ А пХ 1 + - 2'

п=1

(е-т _ в^з„т} СОВ Апх +

+ »3п

+2 [(^41п - А42п) сов -пт + (А41п-п + ^42п^ ЭШ ^т + (^42п - ^41п) в т] —

т1 совАпх

п=1

+

Соответствующие графики приведены на рис. 2, а, б (сплошная линия соответствует х = 0.1, линия с точками — х = 0.5, пунктирная — х = 0.9).

Все указанные выше результаты получены при учете 50 членов рядов. При учёте большего количества членов ряда максимальная относительная погрешность составляет менее 1%.

и-10

о Ш

10 15 20 25

л

3

2.5 2 1.5 1

0.5 0

0

10 15

б

20 25

Рис. 2. Зависимость перемещений (а) и приращения концентрации (б) от времени. Вариант Б

6

4

2

2

0

5

5

а

Работа выполнена при финансовой поддержке РФФИ (проекты 11-08-00064-а и 11-08-90453-Укр_ф_а).

Библиографический список

1. Еремеев В. С. Диффузия и напряжения. М. : Изд-во Моск. ун-та, 1978. 287 с.

2. Князева А. Г. Введение в локально-равновесную термодинамику физико-химических превращений в деформируемых средах. Томск : Томский гос. ун-т, 1996. 146 с.

3. Подстригая Я. С., Павлина В. С. Дифференциальные уравнения термодинамических процессов в п-ком-понентном твёрдом растворе // Физико-химическая механика материалов. 1965. № 4. С. 383-389.

4. Бугаев Н. М., Гачкевич А. Р., Земсков А. В., Тар-лаковский Д. В. Приближённое решение одномерной задачи связанной термоупругой диффузии для полупространства // Проблеми обчислювально1 мехашки 1 мицност1 конструкц1й: зб1рник наукових праць / Дн1-пропетровський нац1ональний ун1верситет. Дн1пропет-ровськ : 1МА-прес, 2011. Вип. 16. С. 60-68.

5. Горшков А. Г., Медведский А. Л., Рабинский Л. Н., Тарлаковский Д. В. Волны в сплошных средах. М. : Физматлит, 2004. 472 с.

The One-dimensional Problem of Unsteady-related Elastic Diffusion Layer A. R. Gachkevich1, A. V. Zemskov2, D. V. Tarlakovsky2

1 Ya. S. Pidstryhach Institute for Applied Problems of Mechanics and Mathematics, Ukraine, 79060, L'vov, Naukova st., 3 b, dept13@iapmm.lviv.ua

2Moscow Aviation Institute (State University of Aerospace Technologies), Russia, 125993, Moscow, GSP-3, A-80, Volokolamskoe Shosse, 4, azemskov1975@mail.ru, tdvhome@mail.ru

The problem of determining the stress strain state of an elastic medium, taking into account the structural changes caused by the presence of diffusion fluxes. The influence of diffusion processes on the stress-strain state of the environment is taken into account by using the locally equilibrium model of thermoelastic diffusion, which includes the coupled system of equations of motion of an elastic body and the equations of heat and mass transfer. For solutions used decompositions of the unknown functions in Fourier series and then applying the integral Laplace transform with respect to time. We construct a fundamental solution of the problem. For examples the cases where the diffusion flux at the boundary is constant, or decays exponentially are considered.

Key words: elastic diffusion, time-dependent problems, Fourier series, Laplace transform. References

1. Eremeev V. S. Diffuziia i napriazheniia [Diffusion and stresses]. Moscow, Moscow Univ. Press, 1978, 287 p. (in Russian).

2. Knjazeva A. G. Vvedenie v lokal'no-ravnovesnuiu termodinamiku fiziko-khimicheskikh prevrashchenii v deformiruemykh sredakh [Introduction to the locally-equilibrium thermodynamics of physical and chemical transformations in deformable environments]. Tomsk, 1996, 146 p. (in Russian).

3. Podstrigach Ja. S., Pavlina V. S. Differential equation of thermodynamic processes in the n-component solid solution. Fiziko-himicheskaja mehanika materialov, 1965, no. 4, pp. 383-389. (in Russian).

УДК 51-73

МЕТОДИКА ОПРЕДЕЛЕНИЯ ОБЛАСТЕЙ, ТРЕБУЮЩИХ КВАНТОВОГО ОПИСАНИЯ В РАМКАХ ГИБРИДНОГО МЕТОДА

О. Е. Глухова1, А. С. Колесникова2, М. М. Слепченков3, Г. В. Савостьянов4

1 Доктор физико-математических наук, заведующий кафедры радиотехники и электродинамики, Саратовский государственный университет им. Н. Г. Чернышевского, oeglukhova@yandex.ru

2 Программист отдела математического моделирования, Саратовский государственный университет им. Н. Г. Чернышевского), kolesnikova.88@mail.ru

3Ассистент кафедры радиотехники и электродинамики, Саратовский государственный университет им. Н. Г. Чернышевского), slepchenkovm@mail.ru

4Аспирант кафедры радиотехники и электродинамики, Саратовский государственный университет им. Н. Г. Чернышевского), follow.a.white.rabbbitt@gmail.com

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

Ключевые слова: квантово-химические методы, молекулярно-механические методы, поле локальных напряжений, силы, углеродные наноструктуры.

4. Bugaev N. M., Gachkevych A. R., Zemskov A. V., Tarlakovsky D. V. Approximated solution of one-dimensional bound thermo elastic diffusion problem for half space. Problemi obchisljuval'noi mehaniki i micnosti konstrukcij: zbirnik naukovih prac', Dnipropetrovs'kij nacional'nij universitet, Dnipropetrovs'k, IMA-pres, 2011, iss. 16, p. 60-68. (in Russian).

5. Gorshkov A. G., Medvedsky A. L., Rabinsky L. N., Tarlakovsky D. V. Volny v sploshnykh sredakh Waves in Continuous Media. Moscow, Fizmatlit, 2004, 472 p. (in Russian).