CC BY
38
Методика расчета упругих волновых полей для мезообъема, содержащего несколько разноориентированных кристаллов с различными свойствами
Н.А. Мельникова, М.М. Немирович-Данченко1,2
Северская государственная технологическая академия, Северск, 636070, Россия
1 Томский политехнический университет, Томск, 634050, Россия 2 Институт геофизики СО РАН, Новосибирск, 630090, Россия
В работе рассмотрены некоторые приемы компьютерного моделирования деформирования анизотропных упругих твердых тел. Приведено численное решение задачи о деформировании кристалла, главная ось которого не совпадает с направлением нагружения. Кроме того, в результате решения модельной задачи для случая плоской деформации выявлено на мезоуровне существенное искажение воздействующего поля, деформирующего анизотропную поликристаллическую среду.
Calculation procedure of elastic wave fields for the mesovolume with several misoriented crystals of different properties
N.A. Melnikova and M.M. Nemirovich-Danchenko12
Seversk State Technological Academy, Seversk, 636070, Russia 1 Tomsk Polytechnic University, Tomsk, 634034, Russia
2 Institute of Geophysics SB RAS, Novosibirsk, 630090, Russia
The paper considers methods for computer simulation of deformation of anisotropic elastic solids. It gives the solution of a problem on deformation of a crystal whose main axis does not coincide with the loading direction. Besides, by solving a model problem for the case of plane deformation we have revealed considerable mesoscopic distortion of the affecting field that deforms the anisotropic polycrystalline medium.
1. Введение
Материалы, прочность которых обычно составляет предмет исследований, являются по сути поликристал-лическими. Как правило, при физико-математической постановке задач для таких материалов свойства каждого кристалла не учитываются, а проводится осреднение этих свойств с построением эффективной однородной модели. Однако целый ряд проблем прочности решается только в масштабах нескольких зерен, блоков. Для такого масштаба характерно изучение свойств каждой неоднородности и их взаимодействия между собой [1].
Многие задачи в механике сплошных сред рассматриваются для поликристаллических тел с хаотичной ориентацией. Такие тела изотропны, матрица упругих модулей, как известно, имеет вид:
А + 2ц А А 0 0 0
А А + 2ц А 0 0 0
А А А + 2ц 0 0 0
0 0 0 ц 0 0
0 0 0 0 ц 0
0 0 0 0 0 ц
причем упругие модули получаются осреднением свойств большого числа разноориентированных кристаллов [2]. Этот подход типичен для задач макроуровня.
Но, особенно в последнее время, ряд проблем деформирования твердых тел формулируется на уровне нескольких кристаллов, зерен, блоков (так называемый мезоуровень). В этом случае кристаллы будут ориенти-
© Мельникова Н.А., Немирович-Данченко М.М., 2GG6
рованы, скорее всего, различным образом, но в силу их малого количества осреднение свойств проводить некорректно.
Для расчета упругого деформирования мезообъема, содержащего несколько разноориентированных кристаллов с различными свойствами, была разработана специальная методика расчетов.
Очевидно, что вид матрицы упругих модулей для кристаллов зависит от ориентации кристалла в системе координат. Итак, для численного решения задач деформирования твердых тел в этом случае необходимо иметь процедуры преобразования компонент тензоров 4-го ранга при повороте системы координат на произвольный угол.
Преобразования такого рода довольно сложны, поэтому остановимся на этом вопросе более подробно.
2. Формулы преобразования матрицы упругих модулей в анизотропной среде при повороте системы координат и их использование в численном методе
Рассмотрим случай поворота одного кристалла с вмороженной, подвижной системой координат х1, х2 относительно неподвижной к, I. Ось х2 направлена вдоль оси симметрии кристалла.
При повороте кристалла изменятся направления осей х1, х2 на угол ю.
Формулы, выражающие старые компоненты тензора напряжений (без штриха) через новые (со штрихом), будут иметь вид:
п gh = Т%Т^, (1)
где
Т = С°8(х', X] ). (2)
Формула (2) определяет матрицу направляющих косинусов размерностью (3x3).
Формулы, выражающие старые компоненты тензора деформаций через новые:
ртп = ТктТ1прк1. (3)
Подставим (1) и (3) в обобщенный закон Гука:
^ = С р
gh сghmn тп>
получим:
т Т п' = с Т Т р'
1ig1 ]к у с ghmn1 кт11п к1 ’
что, вследствие ортогональности преобразования, можно записать в виде:
Приведем начало выражения для элемента с1/111 по формуле (4):
п' = Т Т с Т Т р'
І 1 і§Т jhcg^mn1km1 Ы Ш ’
где
ТщТ]}гсфшпТЫТ1п суЫ ‘
(4)
Здесь g, h, т, п — индексы суммирования.
Выражение (4) — это и есть формула преобразования компонент тензора упругих модулей 4-го ранга при повороте системы координат (см., например, [3]).
с1111
Т11Т11с1 11Т11Т11 + Т11Т11с1112Т11Т12 +
+ 11 11Т Ті 113Т11Т13 + Т11Т11с1121Т12Т11 +
+ Т11Т11с 122Т12Т12 + Т11Т11с1123Т12Т13 +
+ Т11Т11с 131Т13Т11 + Т11Т11с1132Т13Т12 +
+ Т11Т11с 133Т13Т13 + Т11Т12 с1211Т11Т11 +
+ Т11Т12 с 212Т11Т12 + Т11Т12 с1213Т11Т13 +
+ Т11Т12 с 221Т12Т11 + Т11Т12с1222Т12Т12 +
+ <0 Т112 Т11 223Т12Т13 + Т11Т12с1231Т13Т11 +
+ Т11Т12 с 232Т13Т12 + Т11Т12 с1233Т13Т13 +
+ с1 ГГ) Т11 Т11 311Т11Т11 + Т11Т13с1312Т11Т12 +
+ с1 ГГ) Т11 Т11 313Т11Т13 + Т11Т13с1321Т12Т11 +
+ с1 ГГ) Т11 Т11 322Т12Т12 + Т11Т13с1323Т12Т13 +
+ Т11Т13с 331Т13Т11 + Т11Т13с1332Т13Т12 +
+ Т11Т13с 333Т13Т13 + ...
(5)
В выражении (5) приведены первые 27 слагаемых, полученных при пробегании индексами h, т и п значений от 1 до 3, при этом первый индекс g фиксирован. Следующие 27 слагаемых отличаются от приведенных тем, что индекс g фиксирован и равен 2, последние 27 слагаемых образуются при g = 3. Таким образом, в формуле (4) в вычислении каждого элемента матрицы с'н участвует 81 слагаемое.
Рассмотрим пример поворота гексагонального монокристалла вокруг оси у.
Ось симметрии 6-го порядка совместим с осью г, т.е. с 3-й осью декартовой ортогональной системы координат. Тогда полное представление (9x9) тензора модулей упругости для кристалла с гексагональной решеткой будет следующим:
с11 с12 с13 0 0 0 0 0 0
с12 с22 с13 0 0 0 0 0 0
с13 с13 с33 0 0 0 0 0 0
0 0 0 с44 0 0 с44 0 0
0 0 0 0 с55 0 0 с55 0
0 0 0 0 0 с66 0 0 6 6 с
0 0 0 с44 0 0 с44 0 0
0 0 0 0 с55 0 0 с 4 4 0
0 0 0 0 0 с66 0 0 6 6 с
упругий модуль с55 = с 44, а с12 = с11 - 2с66
(6)
Таким образом тензор упругих модулей для кристалла с гексагональной осью симметрии содержит 21 ненулевую компоненту, из них 5 являются независимыми.
Направляющие косинусы вычисляются по формуле
1 = С°8(X', X] ).
Например, при повороте на 10° матрица Ту имеет вид:
" ^10° 0 sin 10° "
0 1 0 - sin 10° 0 cos10°
Т =
Новые упругие константы вычисляются по формуле (4), причем эти вычисления очень трудоемкие, как это видно из выражения (5). Приведем сразу же окончательный вид преобразованной матрицы:
С11 с12 С13 0 с15 0 0 с15 0
с12 С22 С23 0 с25 0 0 с25 0
с13 С13 С зз 0 с35 0 0 с35 0
0 0 0 С44 0 с46 С47 0 с46
с15 с25 с35 0 С55 0 0 С55 0
0 0 0 с46 0 С66 с4б 0 С66
0 0 0 С44 0 с46 С44 0 с46
с15 с25 с35 0 С55 0 0 С55 0
0 0 0 с4б 0 С66 с4б 0 С66
(7)
Сравнение (6) и (7) наглядно показывает, что при повороте системы координат в матрице упругих модулей в общем случае появляются дополнительные члены, которые до этого были нулями. При этом симметричность матрицы относительно главной диагонали сохранилась [4].
Если осуществлять поворот на любой угол вокруг оси г, которую мы совместили с осью симметрии 6-го порядка, то матрица упругих модулей не изменится, так как ось симметрии 6-го порядка является осью симметрии бесконечного порядка.
В целом определяющие соотношения для численного моделирования волновых полей в неоднородной среде с анизотропными включениями произвольной ориентации примут вид:
6 XX = С'и(х’ 2) & XX + С13 (X 2) & 22 + С15 (X 2)& Х2 ,
622 = с1з(X, 2) &XX + сзз(X, 2) &+ СЗ5(х, 2) &х2, (8)
6XX = С15 (X, 2) &XX + С35 (X, 2) &22 + С55 (X 2) &XX •
Для сравнения приведем формулы, использованные в расчетах для среды с анизотропными включениями без учета поворота:
6 XX = 2) & XX + С13 (X 2) & 22 ,
622 = С13 (x, 2) &XX + С33 (X, 2) &22,
62X = С55(X, 2) &2X •
3. Результаты расчета волновых полей
для монокристалла с поворотом оси симметрии
Использование формул (8) позволяет сделать алгоритм расчета волновых полей универсальным. В тех расчетных ячейках (или областях), в которых угол поворота равен нулю, матрица с = с, для других расчетных ячеек матрица упругих модулей должна вычисляться
по формуле (4). Таким образом, при расчете волновых полей входными параметрами будут не только значения упругих модулей и плотности в каждой точке, но и значения угла поворота оси симметрии кристалла. Чтобы не перегружать основной алгоритм большим объемом вычислений, связанным с формулой (4), вычисления были разбиты на два этапа. На первом этапе по заданным упругим модулям и углам поворота подготавливается матрица с для отдельных неоднородностей среды. На втором этапе входным параметром является уже матрица с •
Задачи, обсуждаемые в данном разделе, приводятся в порядке увеличения сложности получаемого волнового поля. Конечной целью моделирования здесь было создание и отладка пакета программ для расчета волновых полей в произвольно неоднородных средах с различной ориентацией кристаллических зерен. В таком расчете волновая картина осложнена волнами различных типов.
Была решена задача о распространении продольной волны в монокристалле цинка с повернутой осью симметрии.
Сначала осуществляли поворот оси симметрии кристалла относительно лабораторной оси симметрии (оси дальнейшего расчета) на 30°. Система координат (х, у) является неподвижной фиксированной (лабораторной), а система (X!, у1) связана с кристаллом. Поворот осей подразумевает применение формулы преобразования (4) к системе (x1, у1). В результате получается матрица (7), компоненты которой подставляются в определяющие соотношения. После всех преобразований волновой процесс рассчитывают в системе (х, у).
Так как пока мы решаем задачу о распространении только продольной волны, картина волнового поля для кристалла с осью симметрии совпадающей с осью симметрии расчета не будет меняться на протяжении всего процесса распространения. После поворота происходят значительные изменения.
На рис. 1 показано векторное поле скоростей в один из начальных моментов времени. Цифрой 1 обозначена распространяющаяся продольная волна, цифрой 2 — поперечная волна, цифрами 3, 4, 5 — головные волны.
По результатам этого расчета можно сделать вывод о том, что в кристалле с повернутой осью симметрии наряду с продольной волной возникает поперечная волна. Убедиться в том, что полученная волна является поперечной, нам позволяют численные снимки горизонтальной и вертикальной компонент скоростей смещения в различные моменты времени (рис. 2).
Этот рисунок представляет собой процесс распространения волны во времени. На численных снимках вертикальных компонент более сильной является продольная волна 1, а на горизонтальных компонентах — поперечная волна 2. Небольшое загибание фронта продольной волны 1 обусловлено влиянием границ вследствие поворота кристалла. Волна, отмеченная цифрой 3, яв-
Рис. 1. Векторное поле скоростей смещения частиц среды в повернутом кристалле при распространении плоской волны. Обозначения даны в тексте
ляется головной продольной волной, а отмеченная цифрой 4—головной поперечной. Природа происхождения этих волн объясняется влиянием углов модели в повернутом кристалле. В неповернутом кристалле головных волн не наблюдается и фронт не загибается.
Более сложная волновая картина получается при действии сосредоточенного источника. Была решена задача о распространении упругих волн от источника вертикальной силы в анизотропном монокристалле, у которого ось симметрии не совпадает с осью симметрии задачи. Это, вообще говоря, аналог известной в эласто-
динамике задачи Лэмба. Результаты этого расчета представлены на рис. 3.
Полученная волновая картина для повернутого кристалла в точности повторяет волновую картину для непо-вернутого кристалла, но с поворотом на угол 30°. Характерные области рефракции поперечной волны, которые у неповернутого кристалла расположены вдоль оси х3, смещены относительно нее также на 30°. Смещение (поворот) фронта продольной волны происходит на такой же угол. Иными словами, по волновой картине можно определить угол поворота кристалла.
Рис. 2. Численные снимки вертикальной (а-в) и горизонтальной (г-е) компонент скорости смещения повернутого монокристалла в различные моменты времени
Рис. 3. Задача Лэмба для кристалла цинка с повернутой осью симметрии
4. Исследование волновых процессов в среде с анизотропными включениями
Модель среды с включением с круто наклоненными границами представляет особый интерес. При рассмотрении волновых процессов в таких средах, в первую очередь, изучают распространение волн вблизи границ раздела. Так как при распространении волн в среде с включением угол раздела также является дифрагирующим объектом, наблюдается волновая картина, типичная для моделей с включениями.
Результаты, ставшие классическими, были получены Б.Г. Михайленко [5] для целого ряда моделей (задачи о падении плоской волны на двухкомпонентные среды). Авторы сочли необходимым для тестирования программы решить сходную задачу и сравнить с результатами Б.Г. Михайленко. Серия исследований была проведена на моделях, состоящих из основной среды с включением, границы которых образуют между собой угол 90°.
В качестве модели была взята среда, представленная на рис. 4. На рисунке структура АВС является средой 2 с задаваемыми различными свойствами. Грани этой среды образуют между собой угол 90°. На верхней части модели задается плоская продольная волна с заданным импульсом.
В среде 1 скорость распространения продольной волны обозначим УР1, скорость распространения поперечной волны У81, в среде 2 — скорость распространения продольной волны УР2, скорость распространения поперечной волны У82. В зависимости от изменений скоростей в средах 1 и 2 изменялся материал среды.
Среди всевозможного разнообразия моделей были выбраны несколько, при это модули упругости с11, с33 и с55 взяты из различных источников (см., например, [6, 7]).
Модель 1. Среда 1 имеет скорости УР1 = 6300 м/с, уи = 3 240 м/с, среда 2 имеет скорость продольной волны вдоль направления оси х УР2 х = 4750 м/с, а скорость поперечной волны У82 = 2320 м/с. Среды 1 и 2 представляют собой модельный материал, изотропный и по продольной, и по поперечной волне.
Модель 2. Среда 1 имеет скорости: УР1 = 6300 м/с, У81 = 3 240 м/с, среда 2 имеет скорости продольной волны: вдоль направления оси х УР2х = 4750 м/с, вдоль направления оси z УР22 = 2920 м/с, а скорость распространения поперечной волны У82 = 2320 м/с. Среда 1 — изотропный материал, а среда 2 — гексагональный кристалл (цинк), ось симметрии которого параллельна оси симметрии задачи (константы цинка в ГПа с11 = 161, с33 = 61, с44 = 38.3, с13 = 50.1).
1
V Р1’''81 90^
А О О Г
2
Vp2, Vs2
О О Т
С
Рис. 4. Геометрия среды
Модель 3. Среда 1 имеет скорости УР1 = 6300 м/с, У81 = 3 240 м/с, среда 2 — анизотропный материал (гексагональный кристалл цинка) с осью симметрии, повернутой на угол 30° относительно оси симметрии всей модели. Константы среды 2 в системе координат, связанной с задачей, таковы: а11 = 141.9, а13 = 44.2, а15 = = -18.2, а33 = 92, а35 = -25, а55 = 32.4.
Следует подробнее остановиться на описании полученных волновых картин для каждой модели.
Для полученных численных снимков были введены следующие обозначения: 1 — падающая продольная волна; 2 — отраженная продольная волна от границы АВ; 3 — преломленная продольная волна через границу АВ; 4 — дифрагированная продольная волна от угла А; 5 — дифрагированная поперечная волна от угла А; 6— головная продольная волна; 7 — головная поперечная волна; 8 — коническая волна в среде 1; 9 — дифрагированная продольная волна от угла А в среде 2; 10 — дифрагированная поперечная волна от угла А в среде 2.
Стрелками показаны направления движения падающей, отраженной и преломленной продольных волн.
На полученных волновых картинах сплошной линией показано реальное расположение границы раздела двух сред.
На рис. 5 представлена полученная волновая картина для модели 1. Видно, что скорость продольной волны 1 в среде 1 больше скорости преломленной волны 2 в среде 2. При распространении продольной волны 1 проис-
ходит отражение от границы АВ, и в противоположную сторону распространяется отраженная продольная волна 3. Как было сказано выше, угол включения является дифрагирующим объектом и от него начинается излучение сферических фронтов волн. Сферическая продольная волна 4 и сферическая поперечная волна 5 в среде имеют вид окружности. По исследованиям с монокристаллом в изотропной среде мы и должны наблюдать именно такую картину. Фронт плоской преломленной волны 2 отстает от фронта падающей плоской волны 1, так как скорость распространения продольной волны в среде 1 больше скорости распространения продольной волны в среде 2. Волна 1, распространяясь в среде 1 вдоль границы АС, генерирует в среде две интенсивные головные волны 6 и 7. Их возникновение объясняется поперечной диффузией энергии вдоль фронта падающей волны 1. Взаимодействие фронта волны 3 с границей АС также вызывает образование конической волны 8 в среде 1. Дифрагированная продольная волна 9 и дифрагированная поперечная волна 10, возникающие в результате взаимодействия падающей волны с неоднородностью (в данном случае угол раздела двух сред), распространяются с меньшими скоростями, чем волны 4 и 5. Вследствие изотропности материала в среде 2 фронты волн 9 и 10 распространяются также по окружности, но меньшего радиуса. Приведенный анализ волновых полей показывает хорошее соответствие полученных данных с результатами Б.Г. Михайленко [5].
Рис. 5. Волновое поле вертикальной компоненты смещения для модели 1 в фиксированный момент времени
Обсуждение результатов расчетов для моделей 2 и 3 лучше проводить, сравнивая полученные волновые картины. Поэтому рассмотрим несколько волновых полей в последовательные моменты времени. В модели 2 включением является анизотропный материал (цинк) с осью симметрии, совпадающей с ось симметрии всей модели, а модель 3 содержит включение, у которого ось симметрии повернута относительно оси симметрии расчета на 30° против часовой стрелки. Механизм осуществления поворота оси симметрии и связанные с этим преобразования были описаны выше.
По приведенным волновым картинам можно сделать вывод о том, что значительные отличия между волнами в модели 2 и волнами в модели 3 происходят в среде 2. Наличие поворота оси симметрии среды 2 в модели 3 приводит к возникновению волн, которые не наблюдались в модели 2. В результате поворота оси симметрии
кристалла (среда 2) возникает поперечная волна, которая распространяется на некотором расстоянии от преломленной продольной волны 2. Так как имеется угол В, то от него возникает волна дифракции.
На рис. 6 показан процесс развития волновых картин во времени для горизонтальных и вертикальных компонент смещения для модели 2. Рисунки 6, а, в, д, ж являются волновыми картинами вертикальных компонент смещения, а рис. 6, б, г, е, з — волновые картины горизонтальных компонент смещения. На рис. 6, а, б показан момент времени, в который плоская волна достигла границы АВ включения и находится прямо на ней. На рис. 6, в, г плоская волна миновала границу, и от границы в разные стороны начали распространяться преломленная и отраженная волны. В основной среде от угла зарождается фронт сферической продольной волны. По мере распространения плоских волн в среде
Рис. 6. Сравнение волновых полей вертикальной (а, в, д, ж) и гори- Рис. 7. Сравнение волновых полей вертикальной (а, в, д, ж) и горизонтальной (б, г, е, з) компонент смещения для модели 2 в различные зонтальной (б, г, е, з) компонент смещения для модели 3 в различные
моменты времени моменты времени
0201024801230002310002020100000200010102000248
0100005353020100020201010202023102020101000202
Рис. 8. Модель кристаллической среды с наложенным волновым полем вертикальной компоненты
и во включении зарождаются головные и конические волны (рис. 6, д, е). От преломленной плоской волны во включении в связи с влиянием границы включения зарождается коническая волна (рис. 6, ж, з).
На рис. 7 показан процесс развития волновых картин во времени для горизонтальных и вертикальных компонент смещения для модели 3. В этой модели, как описывалось выше, включение представляет собой кристалл цинка, у которого ось симметрии повернута относительно оси симметрии основной среды на 30°. Вследствие этого наблюдаются волны, которых не было на волновых картинах для модели 2. На рис. 7, б видим, что присутствует поперечная волна, которая возникает в этой модели.
Итак, исследования волновых полей позволяют понять особенности распространения упругих волн в реальных кристаллитах, оси симметрии которых в общем случае могут не совпадать с нормалью нагружения. Для демонстрации того, как может повести себя при плоском деформировании некоторый мезообъем, состоящий из нескольких зерен, было рассчитано волновое поле от падения плоской волны на модельную среду, состоящую из 12 кристаллов, повернутых на различные углы (см., например, [8]). На рис. 8 приводится исходная структура среды с наложенным на нее численным снимком вертикальной компоненты скорости смещения. Хорошо видно, что первоначально плоский фронт в такой среде преобразуется в искривленный. Подчеркнем, что здесь все 12 кристаллов — одно вещество, исходные модули упругости во всех 12 кристаллах идентичны.
Но различные повороты порождают в матрице (7) различные компоненты. По существу, исходными данными для расчета явились 12 матриц упругих модулей.
Волновая картина на рис. 8 существенно осложнена волнами дифракции, которые выше наблюдались в единичном кристалле.
Таким образом, разработанный алгоритм может быть использован для расчетов распространения упругих волн в поликристаллических средах на мезоуровне.
5. Выводы
Выполнена модификация линеаризованного метода Уилкинса для применения его в решении задач эласто-динамики для анизотропных неоднородных сред. В модели среды возможно наличие произвольного числа разноориентированных кристаллов, зерен, блоков зерен. Конечно-разностная схема упрощена за счет предположения о малости деформаций, естественного в задачах о распространении упругих волн. При этом возможно полагать расчетные ячейки недеформируемыми и упростить выражения для пространственных производных.
Для тестирования алгоритма, модифицированного для расчетов в анизотропных средах, решен аналог задачи Лэмба для однородного анизотропного монокристалла. Результат расчета (фронты продольной и поперечной волн) совпадает с аналитическими поверхностями лучевых скоростей, что доказывает применимость алгоритма.
Изучены процессы распространения упругих волн в кристаллах в случае, когда главная ось анизотропии повернута относительно оси симметрии задачи. Показано, что дифракция плоской волны на таком кристалле (с несовпадением осей симметрии кристалла и задачи) имеет существенные отличия от дифракции плоской волны на кристалле в случае совпадения осей. Этот вывод важен при изучении деформирования поликристаллических тел на мезоуровне, а также при интерпретации волновых полей в задачах физической акустики.
Работа выполнена при поддержке гранта РФФИ № 04-05-64547-а.
Литература
1. Физическая мезомеханика и компьютерное конструирование материалов / Под ред. В.Е. Панина. - Новосибирск: Наука, 1995. -Т. 1. - 298 с., Т. 2. - 320 с.
2. Шермергор Т.Д. Теория упругости микронеоднородных сред. -М.: Наука, 1977. - 400 с.
3. Хирт Дж., Лоте И. Теория дислокаций. - М.: Атомиздат, 1972. -
600 с.
4. Крючкова В.В., Немирович-Данченко М.М. Преобразование матри-
цы упругих модулей в анизотропной среде при повороте системы координат: Методическое пособие. - Томкс: Изд-во ТГУ, 1999. -12 с.
5. Михайленко Б.Г. Сейсмические поля в сложнопостроенных средах. - Новосибирск: Наука, 1988. - 310 с.
6. Сиротин Ю.И., Шаскольская М.П. Основы кристаллофизики. -М.: Наука, 1975. - 680 с.
7. Федоров Ф.И. Теория упругих волн в кристаллах. - М.: Наука, 1965. - 388 с.
8. Смолин И.Ю., Соппа Э., Шмаудер 3., Макаров П.В. Двумерное моделирование пластической деформации в матрице металлокерамического композита на мезоуровне: оценка напряженных состояний и численных методов // Физ. мезомех. - 2000. - Т. 3. -№ 1. - С. 17-22.
