Численное моделирование распространения акустических волн в анизотропных средах Текст научной статьи по специальности «Физика»

Численное моделирование распространения акустических волн

в анизотропных средах

В.В. Крючкова, М.М. Немирович-Данченко

Томский филиал Института геологии нефти и газа СО РАН, Томск, 634055, Россия

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

1. Введение

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

о представительности расчетной области. Целый ряд проблем корректно рассматриваются только в масштабах нескольких зерен, нескольких блоков зерен [1]. Это тот масштабный уровень, на котором говорить о «поликристалличности» (т.е. о макроуровне) нельзя. В то же время, это, с точки зрения кристаллографических особенностей, и не микроуровень.

Хорошо известно, что кристаллы титана, никеля, кадмия, цинка, магния, различные виды керамики обладают значительной анизотропией по упругим свойствам. Так, упругие модули с11 и с33 (это модули сжатия в направлении 1 и 3 осей соответственно) для Zn отличаются в 2.56 раза, в изотропных же средах эти константы одинаковы. Вообще, для кристаллов характерной является анизотропия их свойств, лучше всего это проявляется при распространении упругих (акустических) волн, анизотропия оптических свойств выражена слабее

[2].

Фундаментальные основы для решения задач крис-таллоакустики заложены во второй половине нашего столетия Федоровым Ф.И. (см. монографию [3]). Удается получать и исследовать общее решение для широкого класса практических случаев. Эти результаты находят свое применение в физической акустике, дефектоскопии, пьезотехнике, сейсморазведке.

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

торые аспекты численного моделирования динамических задач для анизотропных сред.

2. Метод расчета

Здесь необходимо для связности изложения привести некоторые положения и обозначения теории волн в анизотропных средах.

Рассмотрим обычную систему уравнений движения:

д2щ д°у

Р—Г = ^— , X = X У, г. дt дху

Здесь щ — вектор смещений; остальные обозначения очевидны.

Обобщенный закон Гука имеет вид о у = с щ є ы , причем модули упругости образуют тензор четвертого ранга, в общем случае в трехмерном пространстве обладающий 81 ненулевой компонентой. Симметрия тензоров деформации и напряжений (т.е. инвариантность по перестановкам у ^ ]і и к1 ^ 1к) позволяет записать закон Гука в виде

1 дик 1 ди1 ди1

°іу = 2 Сі]к1 ~дХ~ + 2 Сі]к1 дХ~ = Сі]к1 дХ~ ■ (1)

Если в общем случае одна компонента напряжений выражается через 9 компонент деформаций (и поэтому 9 х 9 = 81 компонента тензора упругих модулей), то при условии симметричности обоих тензоров различных компонент тензоров напряжений и деформаций будет по 6 и ненулевых упругих модулей — не более 36. В этом случае вводят обозначения

°11 = °1> 022 = 02> 033 = 03>

023 = 032 = 04> °13 = 031 = 05, °12 = 021 = 06>

є11 = Є1> Є22 = Є2> Є33 = Є3>

2є23 = 2є32 = Є4 , 2є13 = 2є31 = Є5, 2є 12 = 2є21 = Є6 .

© Крючкова В.В., Немирович-Данченко М.М., 1999

Двойки связаны с суммированием равных (вследствие симметричности) компонент тензора деформаций. Закон Гука теперь запишется в виде аа = сарер .

Мы будем рассматривать лишь упругие волны, поэтому процесс изоэнтропийный (адиабатический) и тензор упругих модулей инвариантен по перестановкам у ^ И (соотношение Максвелла) и, стало быть, а ^ в. Поэтому в достаточно общем случае анизотропии ненулевых компонент тензора сар — 36, а независимых — 21 (кристаллы триклинной сингонии, например дистен Al2OSiO4, медный купорос Си^04).

При повышении степени симметрии число ненулевых и независимых упругих модулей уменьшается. Так, для тетрагональных кристаллов (олово, индий) ненулевых компонент остается 12, а независимых — 6. Для сред с гексагональной симметрией (кристаллы титана, цинка, кадмия, трансверсально-изотропные материалы — текстолит, керамики) ненулевых компонент 12, различных — 6, независимых — 5. Кубические кристаллы — самые «правильные», у них 3 независимых упругих модуля, ненулевых по-прежнему 12 (алюминий, молибден, свинец, никель, литий). Высшая симметрия — у поликристаллических изотропных сред: ненулевых компонент — 12, различных — 3, независимых — 2. Отметим, что к низшей категории относят триклинную, моноклинную и ромбическую сингонии, к средней категории — тригональную, тетрагональную и гексагональную, и к высшей категории относят кубическую сингонию.

Заметим еще, что практически все металлы и их сплавы кристаллизуются в различных классах кубической или гексагональной сингоний [2]. Примерами сред, обладающих симметрией как гексагональный кристалл, служат однородные изотропные тела с плоскопараллельной системой трещин и поликристаллические образцы с коллинеарной ориентировкой одной оси кристаллов и хаотической ориентировкой других осей.

В дальнейшем будем рассматривать среды гексагональной системы (если говорить о кристаллах) или трансверсально-изотропные среды (если говорить о поликристаллических телах) с вертикальной осью симметрии 6-го порядка.

Матрица упругих модулей сар имеет вид

(

С11 С12 С13 0 0 0 Л

С12 С11 С13 0 0 0

С13 С13 С33 0 0 0

0 0 0 С44 0 0

0 0 0 0 С44 0

0 0 0 0 0 С11 - с12

(2)

Вернемся к уравнениям движения. В качестве оси симметрии 6-го порядка выберем вертикальную ось г и будем рассматривать плоскую задачу (плоская деформация). Уравнения движения примут вид

да хх да р их = —^ +

дх

дz

да да х

Р% =-г^ + -

(3)

дz дх

Закон Гука будем записывать так, как он обычно используется для построения дискретного аналога:

а хх _ С1І&хх + С13Є zz ’

аzz = С13Єхх + С33Єzz > (

аzx = с44^zx .

Оговоримся, что мы отдаем себе отчет, насколько неполной является плоская постановка для кристаллических сред. Нами был разработан алгоритм и отлажена программа для решения такой задачи в трехмерной постановке. Удалось провести несколько пробных расчетов на сетке 50x50x50. Но на таких сетках волны разных скоростей не успевают разделиться, и такие расчеты нельзя признать законченными в смысле качественного описания явлений. Увеличение числа расчетных ячеек невозможно из-за нехватки оперативной памяти. В тоже время некоторые характерные особенности можно изучать и на двумерных моделях (которые, кстати, используются и в физическом моделировании). С другой стороны, часто поликристаллическое тело с повторяющейся структурой заменяют эффективной моделью трансверсально-изотропной среды, в этом случае двумерная постановка корректна.

Расчетная область разбивается на прямоугольные в общем случае ячейки. Скорости смещений частиц будем определять в узлах (], k) сетки, один такой узел помечен квадратиком на рис. 1. По значениям скоростей

к+1

к-1

О о 1 V

Л 1 о 1 л о

І+1

Иногда в качестве независимого элемента выбирают с66, тогда, очевидно, с12 = с11 - 2с66.

Рис. 1. Шаблон для аппроксимации пространственных производных в уравнениях движения

смещений в узлах компоненты тензора скоростей деформаций єхх, є^, єбудут определяться в центрах ячеек (кружочки на рис. 1) и, стало быть, напряжения также подсчитываются в центрах.

да хх

Естественно при этом производную -д определить через значения в точках, отмеченных на рис. 1 крестиком:

> j+1/2, к

-la г

'j-1/2, к

Ах

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

(axx) j+1/2, к 2 {xx)-

/ у+1/2, к+1/2 V ^ хх ! у+1/2, к-1/2 ] •

Окончательно правая часть первого уравнения движения выглядит так:

1 ГГ

2AxAz I

( xx j j+12, k+1/2 + (a xx j j+1/2, k-1/2 Az

jy-1/2, k+1/2 + (axx )j-1/2, к-1j2 Az +

+ (aAx -

> j-1/2, к+1/2

Ax

хх > ]+1!1, к+1/2

Тхх)у+1/2, И-1/2 + (аж)у-1/2, к-1/2

Производные компонент тензора скоростей деформаций определяются аналогично, например

/. \ ди

> j+1/2, к+1/2 дz

j+1/2, к+1/2

( j j+1, k+1 + (Uz j j, k+1) ( j j+1, к + (Uz j j, k)

Az

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

Еще несколько замечаний об алгоритме. При проверке условия устойчивости число Куранта нужно считать по формуле Дt *ота^/Дх , где vmгx — максимальная продольная скорость во всей расчетной области с учетом анизотропии. Расчетная область может состоять из любого числа различным образом ориентированных блоков, из матрицы, содержащей анизотропные вклю-

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

Перейдем к обсуждению собственно расчетов и их результатов.

3. Примеры расчета волновых полей

Какие типы волн и с какими скоростями будут распространяться в соответствии с уравнениями движения

(3) в среде с матрицей упругих модулей (2)? Для ответа на этот вопрос подставим закон Гука в виде (1) в уравнения движения. Получим уравнения движения в перемещениях (смещениях)

д 2 щ

Pu i = cijkl

дxjдxk

()

Решение уравнений (5) будем искать в виде семейств плоских волн

u(r, t) = Ap exp i[ - rotj

(6)

Здесь r — радиус-вектор частицы среды; p — вектор поляризации; k — волновой вектор (вектор пространственных частот); ю — круговая частота; А— амплитуда волн. Стало быть, мы ищем решения в виде волн, распространяющихся в положительном направлении оси времени с частотой f (ю = 2nf). В пространстве эти волны распространяются с «частотами» кх, ку, к2, характеризующими степень наклона фронта волны к соответствующей оси. Например, волновой вектор (, 0, 0) соответствует плоскости с нормалью, колли-неарной оси х. Вводят единичный вектор волновой нормали m так, что к = (2л/А)т . В оптике этот вектор m называется вектором рефракции. Другие встречающиеся названия — медленность (slowness), обратная скорость. Здесь А — длина волны. Возмущение данной частоты в среде пробегает путь А за время T = 1/ f , поэтому величину v = А/ называют фазовой скоростью. Учитывая вышесказанное, формулу (6) запишем в виде

u(r, t) = Ap exp

2n .

А

i( mr - vt)

(7)

Подставив (7) в уравнения (5), получим уравнения Кристоффеля

Р^ 2Р; = СуиЩткР1 > где свертка Сук1ШуШк — тензор Кристоффеля Гй . Так

как р{ = 8яр1, то запишем окончательно

{ -Р^2§а}Р1 = 0 .

Это дисперсионное (секулярное) уравнение позволяет определить все возможные для данного тензора

Кристоффеля (и, стало быть, для среды с тензором Сщ ) фазовые скорости. Точнее говоря, квадраты фазовых скоростей являются собственными значениями приведенного тензора Кристоффеля Гй/р , а вектор поляризации — его собственным вектором.

При вычислении компонент тензора Кристоффеля нужно учитывать, что, например сш3 = сШ1 = с 3Ш = с 3113=

= С44-

В нашем случае матрица тензора Гй имеет размерность 2x2 и составлена из элементов с индексами 1 и 3

( С11т2 + С44т2 (С13 + С44 )тхт 1

(С13 + С44)тхт С44тХ + С33т1 J ’ (8)

а корни биквадратного уравнения

г„ =

Гіі - V2

А 13

выглядят так:

і3

Г33 - v

= 0

^,2 =

Г11 + Г33 і д/(Г11 Г33) 4Г13

Для кристалла цинка [3] с11= 161 ГПа, с33= 62.1 ГПа, с44= 38.3 ГПа, с13= 50.1 ГПа, плотность р = 7.14 г/см3. На рис. 2 приводятся кривые ^(9) и ^ (0) в полярной системе координат. Угол 9 в плоскости (х, z) отсчитывается от оси г, поэтому для вектора волновой нормали имеем тх = sin(9), mz = cos(9).

Внешняя кривая соответствует скорости квазипро-дольной волны, внутренняя — скорости квазипопереч-ной волны. По фазовым скоростям элементарно строятся кривые векторов рефракции (они же кривые нормальных скоростей, кривые обратных скоростей, кривые медленностей). Эти кривые приведены на рис. 3.

Подчеркнем, что в изотропной среде и кривые фазовых скоростей, и кривые медленностей — концентрические окружности.

В работе [5] рассмотрен пространственный случай для кристалла цинка и построены уже не кривые, а поверхности— фазовая, рефракции, а также лучевая (рис. 4). В трехмерном случае имеются три волны — одна квазипродольная (обозначена на рис. 4 буквой Ь) и две квазипоперечных (обозначены Т и Т2). Квази-поперечная волна Т поляризована в плоскости, нор-

Рис. 4. Поверхности фазовых скоростей (а), медленностей (б) и лучевых скоростей (в) кристалла цинка [5]

Р

Рис. 5. Волновая картина для кристалла цинка: R — волна Рэлея; С — «коническая» волна; Р — квазипродольная волна; SV—квазипоперечная волна. Пунктирными стрелками обозначены лучи, проведенные из точки излучения; сплошными линиями обозначены нормали, восстановленные к волновым фронтам в месте пересечения последних с лучами. Жирной стрелкой сверху обозначен источник

мальной к оси 6-го порядка, поэтому в нашей плоской постановке отсутствует.

Показанные на рис. 4 кривые есть сечения поверхностей — поверхности фазовых скоростей (а), поверхности медленностей (б) и поверхности лучевых скоростей (в). Последние поверхности есть не что иное, как фронты распространяющихся волн, ибо энергия в волне переносится по лучам (лучевые скорости иногда

называют групповыми скоростями). Участки вогнутости поверхностей рефракции определяют области рефракции на фронтах.

В изотропной среде лучевые и фазовые скорости совпадают — фронты (волновые поверхности) представляют собой сферы, и нормаль в любой точке фронта коллинеарна лучу, проведенному к этой точке из места излучения.

Рис. 6. Волновая картина для модельной изотропной среды: R — волна Рэлея; С — «коническая» волна; Р — продольная волна; SV— поперечная волна. Пунктирными стрелками обозначены лучи, проведенные из точки излучения. Жирной стрелкой сверху обозначен источник

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

Было рассчитано волновое поле от «точечного» импульсного источника, действующего на поверхности модельной среды с константами кристалла цинка. На рис. 5 приводится численный снимок, состоящий из векторов скоростей смещений и(х, I), изображенных в каждой расчетной точке. Чтобы лучше понять особенности волновой картины для анизотропной среды, необходимо видеть тот же результат для изотропной среды. Поэтому на рис. 6 приводится результат аналогичного расчета для модельной изотропной среды. Выделяются продольная волна, поперечная волна, волна Рэлея и волна, соответствующая в пространственном случае конической. Проведены лучи из точки приложения нагрузки. Видно, что в тех точках, в которых лучи пересекают фронты продольной и поперечной волн, нормаль, восстановленная к фронту, совпадет с направлением луча. Кроме того, можно видеть, что продольная волна всюду поляризована по лучу, а поперечная почти всюду — поперек луча. Фронты (волновые поверхности в пространственном случае и линии в нашем случае) — дуги окружностей.

Здесь необходимо сделать одно замечание по поводу формулы (7). Она задает плоские волны, нормали mj к этим плоскостям однозначно их определяют. В результате же расчета мы получили отнюдь не плоский фронт. Здесь нет противоречия, так как любую функцию, удовлетворяющую условиям Дирихле, можно разложить в ряд Фурье по функциям вида (7). Суперпозиция плоских волн может составить, например, сферический фронт.

Перейдем к расчету для кристалла цинка (рис. 5).

В результате расчета получены следующие волны: квазипродольная, квазипоперечная, коническая и волна Рэлея. Две последние волны появились, как и в случае с изотропной средой, из-за решения задачи в полупространстве. Дисперсионное уравнение не содержит корней для таких волн, так как составлено без учета граничных условий. Если уравнения Кристоффеля дополнить условиями на границах и задать функцию источника, в результирующие формулы для компонент смещений будет входить так называемый знаменатель Рэлея. При расчетах полного волнового поля по таким формулам получаются и волны Рэлея, и конические волны [6].

На рис. 5 указаны области рефракции— это области неоднозначности волновых поверхностей. Рассмотрим луч, опущенный из точки излучения вертикально вниз

(0 = 0). Двигаясь по лучу, мы встретим квазипопереч-ные колебания, распространяющиеся с разными скоростями, затем — продольные колебания (в этом направлении квазипродольная волна является чисто продольной). Видно, что, если восстановить нормали из точек пересечения вертикального луча с фронтами, то для продольной волны нормаль коллинеарна лучу, а для двух ветвей квазипоперечной волны проходит под значительным углом к вертикали (около 30°). Если обратиться к рис. 3, можно видеть, что углы наклона 60°, 120° и 240°, 300°, то есть отклоняющиеся от вертикали на 30°, как раз ограничивают вогнутости. Вогнутость поверхности обратных скоростей соответствует рефракции на фронтах. На рис. 3 видно так же, что у квази-продольной волны областей рефракции нет.

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

Подытоживая обсуждение результатов, можно сказать, что наблюдатель, фиксирующий приход колебаний на противоположной стороне от источника (нижняя сторона на рис. 5), будет регистрировать первое вступление почти продольных колебаний и несколько (до четырех) вступлений колебаний, связанных с квази-поперечными волнами. Эти особенности распространения упругих волн в анизотропных средах должны учитываться при проведении расчетов на мезоуровне и при обсуждении результатов таких расчетов.

Таким образом, построенный алгоритм позволяет точно рассчитывать распространение упругих (акустических) волн в анизотропных средах.

Литература

1. Физическая мезомеханика и компьютерное конструирование материалов: В 2-х т. / Под ред. В.Е. Панина. - Новосибирск: Наука, 1995.- 297 с. и 320 с.

2. Сиротин Ю.И., Шаскольская М.П. Основы кристаллофизики. -М.: Наука, 1975.- 680 с.

3. Федоров Ф.И. Теория упругих волн в кристаллах. - М.: Наука, 1965.- 388 с.

4. Немирович-Данченко М.М., Стефанов Ю.П. Применение конечно-разностного метода в переменных Лагранжа для численного расчета волновых полей в сложнопостроенных средах // Геология и геофизика. - 1995. - № 11. - С. 94-105.

5. Musgrave M.J.P On the propagation of elastic waves in aeolotropic media // Proc. Roy. Soc. - 1954. - V. 226. - No. 1166. - P. 339-366.

6. Немирович-Данченко М.М. Методика расчета P-, SV-, и SH-волн в дальней зоне при вертикальном и горизонтальном воздействиях на поверхности трансверсально-изотропного полупространства // Исследования распространения сейсмических волн в анизотропных средах. - Новосибирск: Наука, 1992. - C. 71-81.