О восстановлении упругих модулей по экспериментальным значениям фазовых скоростей упругих $p$- и $s$-волн Текст научной статьи по специальности «Физика»

Известия Тульского государственного университета Естественные науки. 2015. Вып. 1. С. 82-90

ФизикА =

УДК 534.2, 534.16

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

И.Ю. Зель, Т.И. Иванкина, Д.М. Левин

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

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

1. Введение

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

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

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

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

2. Методы вычисления

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

8г = / (^Р^ав} У(Я1)тва81 У(Я2)тва81 пг)--/ (^Р^аС У(8!)са1су У(32)са1су пг) >

1 Г

х = 82 ^ ™п'

г=!

где г — количество независимых направлений измерения скоростей, /(Ур,УдьУд2) некоторая функция скоростей продольных (Ур) и поперечных (V?!, У^г) волн, которая будет определена ниже, п - волновая нормаль.

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

Методы вычисления упругих модулей могут различаться как выбором функции минимизации /, так и алгоритмом, обеспечивающим пошаговое решение обратной задачи. В качестве функции / в различных работах [6-8] использовались инварианты тензора Кристоффеля Г ¿г, которые являются коэффициентами кубического уравнения для определения фазовых скоростей:

ёе! (Г - А«I) = - Г!М + ГгХа - Гз = 0,

(1)

где Г1, Г2, Г3 — инварианты тензора Кристоффеля Г^ = С^иги^ии соответственно, Аа = рУ&, а — обозначение Р, 51, 52-волн, I — единичная матрица, р — плотность.

Главным преимуществом в выборе данных инвариантов является простая аналитическая связь скоростей и упругих модулей, вытекающая непосредственно из уравнения (1):

¡1 (С, п)=Г1 = Г„ = Ар + А51 + Аs2, (2)

/2 (С, п) = Г2 = - (ГнГц - ГЙГН) = АрА$1 + АрАв2 + Аз1Аз2, (3)

¡3 (С, п) =Гз = ёе1 (Г^) = АрА^А^• (4)

В то время как в Г2, Г3 включены все компоненты тензора упругих модулей, след Г1 содержит только диагональные элементы Г^, что соответствует 15 неизвестным компонентам тензора С . Поэтому дополнительно к уравнению (2) в работе [6] было предложено использовать полученное методом теории возмущения выражение для скоростей продольных волн в первом приближении:

Ар (С, п)=ГйигЩ. (5)

Это аналогично разложению в ряд Тейлора

Ар (а п) » Ар (Со.п) + (^ (С - ед = (Щ & С

где Со - упругие модули изотропной среды (нулевое приближение).

Для скоростей поперечных 5-волн также можно найти аналитическое представление, основываясь на том, что скорости продольных волн определены (например, формулой (5)). Тогда, решая уравнение (1) относительно скоростей 51 и 52 волн, получим

А^1,52 = 2 (Г1 - Ар ± сА), (6)

где сА = ^(Г - Ар)2 + 4Ар (Г1 - Ар) - 4Г2.

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

1 Р,Я!,Я2 г

X = Г X] (А(а)тваз - А(а)са1с) ^ ™1п • (7)

а г=!

Во всех случаях пошаговое приближение С[к+!] = С[к + ДС[к строится на основе известных алгоритмов, таких как методы Гаусса-Ньютона, Левенберга-Марквардта или различные градиентные методы [9,10].

Вычисление поправок к компонентам тензора упругости на каждом этапе итерационного процесса осуществляется посредством выражения ДС = = [3Т3 + тI) Зт ■ Д/п, где 3 = — матрица Якоби, /п обозначает вектор-столбец, элементы которого составляют значения функции / (Ур, Ув 1, У82) для направления п, Д/ - разница между экспериментальными и вычисленными значениями /, т - скалярный параметр, который задается на каждой итерации. Способ задания параметра т определяется выбранным алгоритмом. Например, для метода Гаусса-Ньютона т = 0, а в методе Левенберга-Марквардта т выбирается так, чтобы на к-й итерации выполнялось условие хк1 < х[и_!].

В работе [10] параметр т предлагается определять как тк+!] = тк1 /р, где р — некая постоянная.

Для метода наискорейшего спуска ДС = -ц ■ дс, где Ц принимает такие

значения, чтобы минимизировать функцию х (ц) = X (С - ц ■ дс

Различия в подходах к вычислению поправок на каждой итерации для указанных алгоритмов может сильно влиять на сходимость и точность вычислений. В методе Гаусса-Ньютона симметричная матрица 3Т3 может быть плохо обусловленной, а для метода Левенберга-Марквардта задание начального значения параметра и постоянной р является в достаточной степени произвольным и зависит, в частности, от качества входных данных.

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

1. ДС[к < ДСт1п или ||ДС[и]| < е,

2. и1! < £*,

3. х[к] < Хт1п,

4. N ^ ^тах,

где ||А|| = ^¿г(ЛТА), е, £*, хтт> Nmax — пороговые значения параметров.

Условия 2, 3 и 4 неоднозначны, поскольку выбор ограничений для этих неравенств достаточно произволен и зависит от качества экспериментальных данных. Критерий 1 также может привести к остановке программы до фактического достижения минимума х. Тем не менее, значение е можно оценить исходя из точности ультразвуковых экспериментов и положить равной 10-4 ГПа.

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

3. Численный анализ

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

Для проведения расчетов выберем 15 ° сетку направлений, реализуемую в экспериментах на сферических образцах [3,4]. Далее для численного анализа требуются табличные или другие хорошо известные данные упругих модулей анизотропных сред. В нашем случае проведем моделирование на примере образца, представляющего собой монокристалл кварца, упругие модули которого приведены в [11]. Используя уравнение (1), вычислим скорости упругих P- и 5-волн для выбранной сетки направлений. Эти значения скоростей заменяют в данном моделировании значения скоростей, полученные в эксперименте, и являются исходными данными (input). Используя алгоритм наискорейшего спуска, рассчитаем компоненты тензора упругих модулей (output) по различным методикам (calculus) (рис.1):

1) метод Xi: решение системы из 2r уравнений (2) и (5) для скоростей поперечных и продольных волн;

2) метод X2: решение системы r уравнений (3);

3) метод Хз: решение системы r уравнений (4);

4) метод Х4: решение 3-х систем r уравнений:

{\р = \р (Cijkt, n), Asi = Asi {Cijkh n)

AS2 = AS2 {Cijkh n)

на основе аналитических зависимостей фазовых скоростей от упругих модулей (5), (6).

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

Степень отклонения значений восстановленных скоростей V* от исходных табличных данных V определяется относительными коэффициентами

* \V* - Va| |AV|

en = max-—- = max ———

V V

для входных данных V и

\V[k]- Va\

max

для вычисленных на к-итерации Уа 1 после завершения работы программы.

На рис. 2 показаны распределения значений е[к для скоростей Р и 5-волн, полученные в результате восстановления тензора упругих модулей кварца методикой (4), в зависимости от погрешности входных данных для скоростей поперечных волн. При этом для продольных волн максимальная погрешность для входных данных составила 60 м/с.

Рис. 1. Схема численного анализа методов восстановления упругих

модулей

Для сравнения различных методик вычисления при одинаковых отклонениях в начальных данных на рис. 3 представлены пространственные распределения скоростей упругих волн до (input) и после (output) восстановления упругого тензора.

Рис. 2. Изолинии значений е[к] при различных относительных отклонениях во входных данных для скоростей поперечных волн

4. Обсуждение

Рассматривая задачу восстановления упругих модулей по значениям скоростей упругих Р- и 5-волн, следует отметить зависимость метода вычисления от самой экспериментальной методики ультразвуковых измерений. Для методов, использующих х1, х2, хз, требуется определение скоростей трех типов волн в каждом направлении, но в случае х1 количество измерений

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

Р 51 52

п 11 шшш пшииии пинав

5.30 7.33 3.57 5.43 2.88 4.58

Рис. 3. Пространственные распределения скоростей Р и ¿"-волн, рассчитанные на основе компонент тензора упругости, восстановленного с помощью различных методик, км/с

Влияние погрешности входных (экспериментальных) данных на результаты вычислений для различных методов (рис. 3) обусловлено непосредственно видом функции f. Высокие значения относительных отклонений скоростей Р и 5-волн для методик Х2, Хз (рис. 3) связаны с оптимизацией разницы произведений квадратов скоростей, что, во-первых, повышает степень нелинейной зависимости f от упругих модулей, а, во-вторых, ведет к увеличению не только для 5, но и для Р-волн. По сравнению с остальными методами использование методики Х4 позволяет получить наименьшие значения е[к^.

Выражение (7) можно представить в виде х4 = Хр + Хяг + Хя2.

Приведение Х4 к минимуму с помощью последовательных итераций возможно только при достижении минимальных значений Х для всех трех типов волн, что, в свою очередь, ведет к оптимальным значениям упругих модулей,

ответственных за скорости продольных и поперечных волн. Данный метод также устойчив относительно увеличения погрешностей входных данных (рис. 2), а результирующие значения e[k] для P и 5-волн оказываются в 2-3 раза меньше, чем исходные отклонения e*.

5. Заключение

Проведенный численный анализ позволил определить оптимальный метод восстановления компонент тензора упругости из значений скоростей P- и 5-волн при решении обратной задачи уравнения Кристоффеля: решение 3-х систем уравнений для скоростей продольных и поперечных волн, измеренных на 15° сетке направлений, с использованием алгоритма наискорейшего спуска.

Список литературы

1. Анизотропия и текстура оливиносодержащих мантийных пород при высо-ких давлениях / А.Н. Никитин, Т.Н. Иванкина, Д.Е. Буриличев, К. Клима, Т. Ло-каичек, З. Прос // Физика Земли. 2001. №1. С. 64-78.

2. King M. S. Recent developments in seismic rock physics // Int. J. Rock Mech. Min. Sci. 2009. V.46. P. 1341-1348.

3. Pros Z., Lokajicek T, Klima K. Laboratory approach to the study of elastic anisotropy on rocks samples // Pure appl. geophys. 2008. № 151. P. 619-629.

4. Lokajicek T, Svitek T. Laboratory measurement of elastic anisotropy on spherical rock samples by lon-gitudinal and transverse sounding under confining pressure // Ultrasonics. 2015. V. 56. P. 294-302.

5. Александров К.С., Продайвода Г.Т. Анизотропия упругих свойств минералов и горных пород. Новосибирск: Изд-во СО РАН, 2000. 354 с.

6. Klima K. The computation of the elastic constants of an anisotropic medium from the velocities of body waves // Studia geoph. et geod. 1973. V. 17. P. 115-122.

7. Every A. G, Sachse W. Sensitivity of inversion algorithms for recovering elastic constants of anisotropic solids from longitudinal wavespeed data // Ultrasonics. 1992. V. 30. №. 1. P. 43-48.

8. Every A. G. Determination of the elastic constants of anisotropic solids // NDT&E International. 1994. V. 27. №. 1. P. 3-10.

9. Мину М. Математическое программирование. Теория и алгоритмы. М.: Наука, 1990. 488с.

10. Pujol J. The solution of nonlinear inverse problems and the Levenberg-Marquardt method // Geophysics. 2007. V. 72. №. 4. P. 1-16.

11. Беликов Б. П., Александров К. С., Рыжова Т. В. Упругие свойства породообразующих минералов и горных пород. М.: Наука, 1970. 276 с.

Зель Иван Юрьевич ([email protected]), аспирант, кафедра физики, Тульский государственный университет.

Иванкина Татьяна Ивановна ([email protected]), к.ф.-м.н., старший научный сотрудник, лаборатория нейтронной физики им. И. М. Франка, Объединенный институт ядерных исследования, Дубна.

Левин Даниил Михайлович ([email protected]), д.ф.-м.н., профессор, кафедра физики, Тульский государственный университет..

On the recovery of elastic constants from phase velocity data of elastic P and S-waves

I.Yu. Zel, T.I. Ivankina, D.M. Levin

Abstract. The analysis of different inversion techniques of elastic moduli of anisotropic materials from phase velocities is presented. In the described calculation procedure different methods and minimization algorithms were applied. Based on the numerical analysis the optimal recovery method for components of the tensor of elasticity has been installed.

Keywords: elastic anisotropy, inversion of elastic moduli, phase velocities.

Zel Ivan ([email protected]), postgraduate student, department of physics, Tula State University.

Ivankina Tatiana ([email protected]), candidate of physical and mathematical sciences, senior researcher, Frank laboratory of neutron physics, Joint Institute for Nuclear Research, Dubna.

Levin Daniil ([email protected]), doctor of physical and mathematical sciences, professor, department of physics, Tula State University.

Поступила 10.11.2014