close

Вход

Забыли?

вход по аккаунту

?

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

код для вставкиСкачать
Известия Тульского государственного университета
Естественные науки. 2015. Вып. 1. С. 82–90
Физика
УДК 534.2, 534.16
О восстановлении упругих модулей
по экспериментальным значениям
фазовых скоростей упругих P - и S-волн
И.Ю. Зель, Т.И. Иванкина, Д.М. Левин
Аннотация. Проведен анализ методов восстановления компонент тензора упругих модулей анизотропных материалов по экспериментальным результатам измерений фазовых скоростей в ультразвуковых экспериментах. Представлена схема расчета упругих
модулей с помощью различных методов и алгоритмов. На основе численного анализа влияния погрешностей на результаты вычислений
установлен оптимальный метод восстановления компонент тензора
упругости.
Ключевые слова: анизотропия упругих свойств, упругие модули,
скорости упругих волн.
1. Введение
Для решения задач геофизики и физического материаловедения в изучении явления упругой анизотропии неоднородных поликристаллических
материалов (металлов, сплавов, керамик, горных пород и др.) наибольшее
распространение получил метод импульсного прозвучивания [1,2]. При этом
с использованием пьезоэлектрических преобразователей производится генерация и регистрация упругих импульсов, прошедших через образец, и
определяются различные параметры зарегистрированных упругих волн, такие как времена вступления и значения скорости, а также спектральные
характеристики акустических сигналов. В работах [3,4] для изучения анизотропии упругих свойств геологических материалов был предложен метод
измерения скоростей продольных P и поперечных S волн на сферических
образцах. Данная методика дополнительно включает систему пошагового
двухкоординатного вращения образца и систему окружения образца (камера
высокого давления, термостат и др.).
Фазовая скорость упругих волн, распространяющихся в анизотропной
среде, определяется тензором модулей упругости этой среды. Поэтому экспериментальные значения скоростей P и S волн, измеренные в различных
направлениях, могут быть использованы для определения компонент тензо-
О восстановлении упругих модулей
83
ра упругих модулей. Эта задача является обратной решению классического
уравнения Кристоффеля, связывающего скорости упругих волн и компоненты тензора упругости [5]. В произвольно выбранном направлении зависимость скоростей от упругих модулей является нелинейной, поэтому решение
обратной задачи основывается на применении численных алгоритмов.
На решение и на вычислительный процесс могут влиять несколько факторов: непосредственно выбранный алгоритм, погрешности значений экспериментальных данных, особенности экспериментальной методики (регулярность сетки направлений и ее частота, типы измеренных волн) и степень
анизотропии образца. В данной работе проводится анализ влияния перечисленных выше факторов на процедуру и погрешность определения упругих
модулей анизотропных сред на сферических образцах, во-первых, для установления необходимых требований к методике эксперимента и соответствующему методу вычисления упругих модулей и, во-вторых, для определения
оптимального алгоритма вычисления.
2. Методы вычисления
Расчет модулей упругого тензора Cijkl по данным экспериментально
измеренных скоростей упругих волн представляет собой обратную задачу,
общее решение которой заключается в минимизации разницы между экспериментальными и расчетными данными:
¢
¢
¡
¡
δi = f V(P )meas , V(S1)meas , V(S2)meas , ni − −f V(P )calc , V(S1)calc , V(S2)calc , ni ,
χ=
r
1X 2
δi → min,
r
i=1
где r — количество независимых направлений измерения скоростей,
f (VP ,VS1 ,VS2 ) некоторая функция скоростей продольных (VP ) и поперечных
(VS1 , VS2 )волн, которая будет определена ниже, n – волновая нормаль.
Процедура расчета строится на итерационном приближении к оптимальному значению тензора упругих модулей C с использованием первых и
вторых производных функции f по упругим модулям (Якобиан и Гессиан).
Поскольку в общем случае неизвестен 21 упругий модуль, необходимо решить определенную (21 уравнение) или переопределенную (>21 уравнений)
систему уравнений.
Методы вычисления упругих модулей могут различаться как выбором
функции минимизации f , так и алгоритмом, обеспечивающим пошаговое
решение обратной задачи. В качестве функции f в различных работах [6-8]
использовались инварианты тензора Кристоффеля Γil , которые являются
коэффициентами кубического уравнения для определения фазовых скоростей:
det (Γ − λα I) = λ3α − Γ1 λ2α + Γ2 λα − Γ3 = 0,
(1)
84
И.Ю. Зель, Т.И. Иванкина, Д.М. Левин
где Γ1 , Γ2 , Γ3 — инварианты тензора Кристоффеля Γil = Cijkl nj nk соответственно, λα = ρVα2 , α — обозначение P , S1, S2–волн, I — единичная матрица,
ρ — плотность.
Главным преимуществом в выборе данных инвариантов является простая
аналитическая связь скоростей и упругих модулей, вытекающая непосредственно из уравнения (1):
f1 (C, n) = Γ1 = Γll = λP + λS1 + λS2 ,
f2 (C, n) = Γ2 =
1
(Γll Γii − Γil Γli ) = λP λS1 + λP λS2 + λS1 λS2 ,
2
f3 (C, n) = Γ3 = det (Γil ) = λP λS1 λS2 .
(2)
(3)
(4)
В то время как в Γ2 , Γ3 включены все компоненты тензора упругих
модулей, след Γ1 содержит только диагональные элементы Γil , что соответствует 15 неизвестным компонентам тензора C . Поэтому дополнительно
к уравнению (2) в работе [6] было предложено использовать полученное
методом теории возмущения выражение для скоростей продольных волн в
первом приближении:
λP (C, n) = Γil ni nl .
(5)
Это аналогично разложению в ряд Тейлора
µ
¶
µ
¶
∂λP
∂λP
λP (C, n) ≈ λP (C0 , n) +
(C − C0 ) =
C,
∂C C0
∂C C0
где С0 – упругие модули изотропной среды (нулевое приближение).
Для скоростей поперечных S−волн также можно найти аналитическое
представление, основываясь на том, что скорости продольных волн определены (например, формулой (5)). Тогда, решая уравнение (1) относительно
скоростей S1 и S2 волн, получим
λS1,S2 =
1
(Γ1 − λP ± dλS ) ,
2
(6)
q
где dλS = (Γ1 − λP )2 + 4λP (Γ1 − λP ) − 4Γ2 .
Выражение (6) является общим решением уравнения (1), если известны
скорости продольных P -волн. Возможность аналитического представления
производных скоростей по упругим модулям позволяет использовать вид
суммы квадратов невязок χ, традиционный для методов, использующих
численное дифференцирование:
О восстановлении упругих модулей
χ=
P,S1,S2 r
¢2
1 X X¡
λ(α)meas − λ(α)calc → min .
r α
85
(7)
i=1
Во всех случаях пошаговое приближение C[k+1] = C[k] + ∆C[k] строится на основе известных алгоритмов, таких как методы Гаусса–Ньютона,
Левенберга–Марквардта или различные градиентные методы [9,10].
Вычисление поправок к компонентам тензора упругости на каждом этапе итерационного процесса осуществляется посредством выражения ∆C =
¡
¢−1 T
n
= JT J + τ I
J · ∆fn , где J = ∂f
∂C — матрица Якоби, fn обозначает
вектор-столбец, элементы которого составляют значения функции f (Vp,
Vs1, V s2) для направления n, ∆f – разница между экспериментальными и
вычисленными значениями f , τ – скалярный параметр, который задается
на каждой итерации. Способ задания параметра τ определяется выбранным алгоритмом. Например, для метода Гаусса–Ньютона τ = 0, а в методе
Левенберга–Марквардта τ выбирается так, чтобы на k-й итерации выполнялось условие χ[k] < χ[k−1] .
В работе [10] параметр τ предлагается определять как τ [k+1] = τ [k] /p, где
p — некая постоянная.
∂χ
Для метода наискорейшего спуска ∆C = −µ · ∂C
, где
такие
³ µ принимает
´
∂χ
значения, чтобы минимизировать функцию χ (µ) = χ C − µ · ∂C
.
Различия в подходах к вычислению поправок на каждой итерации для
указанных алгоритмов может сильно влиять на сходимость и точность вычислений. В методе Гаусса–Ньютона симметричная матрица JT J может быть
плохо обусловленной, а для метода Левенберга–Марквардта задание начального значения параметра и постоянной p является в достаточной степени
произвольным и зависит, в частности, от качества входных данных.
Поскольку сходимость, вообще говоря, не будет конечной (из-за погрешностей в экспериментальных данных и ограничений в машинной точности),
необходимо установить критерии остановки программы. Наиболее употребляемыми критериями являются
следующие:
°
°
1. ∆C[k] 6 ∆Cmin или °∆C[k] ° 6 ε,
° °
2. °J[k] ° 6 ε∗ ,
3. χ[k] 6 χmin ,
4. N > N
qmax ,
где kAk = tr(AT A), ε, ε∗, χmin , Nmax — пороговые значения параметров.
Условия 2, 3 и 4 неоднозначны, поскольку выбор ограничений для этих
неравенств достаточно произволен и зависит от качества экспериментальных
данных. Критерий 1 также может привести к остановке программы до фактического достижения минимума χ. Тем не менее, значение ε можно оценить
исходя из точности ультразвуковых экспериментов и положить равной 10−4
ГПа.
86
И.Ю. Зель, Т.И. Иванкина, Д.М. Левин
В качестве алгоритма вычисления выберем метод наискорейшего спуска
по двум причинам: во-первых, этот алгоритм обеспечивает устойчивую сходимость, во-вторых, вычисление параметра µ производится аналитически.
3. Численный анализ
Для проведения анализа влияния различных факторов на результаты
вычисления упругих модулей из данных, полученных в ультразвуковых экспериментах, необходимо рассмотреть различные виды анизотропных сред,
погрешности экспериментальных данных и установить сетку направлений, в
которых определяются скорости в эксперименте.
Для проведения расчетов выберем 15˚ сетку направлений, реализуемую
в экспериментах на сферических образцах [3,4]. Далее для численного анализа требуются табличные или другие хорошо известные данные упругих
модулей анизотропных сред. В нашем случае проведем моделирование на
примере образца, представляющего собой монокристалл кварца, упругие модули которого приведены в [11]. Используя уравнение (1), вычислим скорости
упругих P - и S-волн для выбранной сетки направлений. Эти значения скоростей заменяют в данном моделировании значения скоростей, полученные в
эксперименте, и являются исходными данными (input). Используя алгоритм
наискорейшего спуска, рассчитаем компоненты тензора упругих модулей
(output) по различным методикам (calculus) (рис.1):
1) метод χ1 : решение системы из 2r уравнений (2) и (5) для скоростей
поперечных и продольных волн;
2) метод χ2 : решение системы r уравнений (3);
3) метод χ3 : решение системы r уравнений (4);
4) метод χ4 : решение 3-х систем r уравнений:
(
λP = λP (Cijkl , n),
λS1 = λS1 (Cijkl , n),
λS2 = λS2 (Cijkl , n),
на основе аналитических зависимостей фазовых скоростей от упругих модулей (5), (6).
Для учета влияния погрешностей расчета к входным данным значений
скоростей добавляются случайные отклонения ∆V (errors).
Степень отклонения значений восстановленных скоростей Vα∗ от исходных табличных данных Vα определяется относительными коэффициентами
e∗α = max
|Vα∗ − Vα |
|∆V |
= max
Vα
Vα
для входных данных Vα и
[k]
e[k]
α = max
|Vα − Vα |
Vα
О восстановлении упругих модулей
87
[k]
для вычисленных на k-итерации Vα после завершения работы программы.
На рис. 2 показаны распределения значений e[k] для скоростей P и S-волн,
полученные в результате восстановления тензора упругих модулей кварца
методикой (4), в зависимости от погрешности входных данных для скоростей
поперечных волн. При этом для продольных волн максимальная погрешность для входных данных составила 60 м/с.
Рис. 1. Схема численного анализа методов восстановления упругих
модулей
Для сравнения различных методик вычисления при одинаковых отклонениях в начальных данных на рис. 3 представлены пространственные распределения скоростей упругих волн до (input) и после (output) восстановления
упругого тензора.
Рис. 2. Изолинии значений e[ k] при различных относительных
отклонениях во входных данных для скоростей поперечных волн
4. Обсуждение
Рассматривая задачу восстановления упругих модулей по значениям скоростей упругих P - и S-волн, следует отметить зависимость метода вычисления от самой экспериментальной методики ультразвуковых измерений. Для
методов, использующих χ1 , χ2 , χ3 , требуется определение скоростей трех
типов волн в каждом направлении, но в случае χ1 количество измерений
88
И.Ю. Зель, Т.И. Иванкина, Д.М. Левин
скоростей поперечных волн может быть меньше, чем для продольных волн.
Эти методы также ограничены необходимостью измерения скоростей двух
поперечных S-волн (быстрая S1 и медленная S2). Преимуществом использования методики χ4 в этом случае является независимость от количества и
типа волн, измеренных в одном направлении.
Рис. 3. Пространственные распределения скоростей P и S-волн,
рассчитанные на основе компонент тензора упругости, восстановленного
с помощью различных методик, км/с
Влияние погрешности входных (экспериментальных) данных на результаты вычислений для различных методов (рис. 3) обусловлено непосредственно видом функции f . Высокие значения относительных отклонений
скоростей P и S-волн для методик χ2 , χ3 (рис. 3) связаны с оптимизацией
разницы произведений квадратов скоростей, что, во-первых, повышает степень нелинейной зависимости f от упругих модулей, а, во-вторых, ведет к
увеличению e[k] не только для S, но и для P -волн. По сравнению с остальными методами использование методики χ4 позволяет получить наименьшие
значения e[k] .
Выражение (7) можно представить в виде χ4 = χP + χS1 + χS2 .
Приведение χ4 к минимуму с помощью последовательных итераций возможно только при достижении минимальных значений χ для всех трех типов
волн, что, в свою очередь, ведет к оптимальным значениям упругих модулей,
О восстановлении упругих модулей
89
ответственных за скорости продольных и поперечных волн. Данный метод
также устойчив относительно увеличения погрешностей входных данных
(рис. 2), а результирующие значения e[k] для P и S-волн оказываются в 2-3
раза меньше, чем исходные отклонения e∗ .
5. Заключение
Проведенный численный анализ позволил определить оптимальный метод восстановления компонент тензора упругости из значений скоростей P - и
S-волн при решении обратной задачи уравнения Кристоффеля: решение 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 с.
Зель Иван Юрьевич (ivangreat2009@gmail.com), аспирант, кафедра физики, Тульский государственный университет.
90
И.Ю. Зель, Т.И. Иванкина, Д.М. Левин
Иванкина Татьяна Ивановна (iti@jinr.ru), к.ф.-м.н., старший научный
сотрудник, лаборатория нейтронной физики им. И. М. Франка, Объединенный институт ядерных исследования, Дубна.
Левин Даниил Михайлович (levin@physics.tsu.tula.ru), д.ф.-м.н., профессор, кафедра физики, Тульский государственный университет..
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 (ivangreat2009@gmail.com), postgraduate student, department of
physics, Tula State University.
Ivankina Tatiana (iti@jinr.ru), candidate of physical and mathematical
sciences, senior researcher, Frank laboratory of neutron physics, Joint Institute
for Nuclear Research, Dubna.
Levin Daniil (levin@physics.tsu.tula.ru), doctor of physical and mathematical
sciences, professor, department of physics, Tula State University.
Поступила 10.11.2014
Документ
Категория
Без категории
Просмотров
6
Размер файла
1 134 Кб
Теги
экспериментальной, значения, скоростей, восстановлен, волна, модулем, упругие, фазовые
1/--страниц
Пожаловаться на содержимое документа