четверг, 11 сентября 2014 г.

Скорости сейсмических волн из пробирки


Bridging the Gap
девиз DARPA

Самая безнадёжная, а потому и самая увлекательная задача при работе с какой-либо моделью (набором уравнений) – это её проверка всеми возможными методами. 
Я занимаюсь, как я писал в прошлых постах (крайний и далее по ссылкам), составлением компьютерной модели плавления вещества в условиях мантии Земли. На каждом этапе необходимо проверять, насколько сильно я вру. Один из неиспользованных мною до сих пор источников данных – это профили скоростей сейсмических волн. Поскольку для Земли есть подобные профили (та же модель PREM), и это единственный прямой источник данных о её структуре, его использовать необходимо. 
Что может дать данное сравнение? Первое – это дополнительная проверка, что вся модель запрограммирована верно. Второе – это выход на какие-то экспериментально измеряемые величины.

Итак, скорости объемных сейсмических волн. Их две - продольные \(V_p\)  и поперечные \(V_s\) (обозначения от слов primary и secondary - в соответствии с порядком, в котором они приходят на сейсмодатчик). Скорости волн определяются как:
\[V_p = \sqrt {\frac {K_s + \frac{4}{3}G}{\rho}} \]
\[V_s = \sqrt {\frac {G}{\rho}} \]
\[V_p > V_s \]

Как видно из формул, они зависят от модуля упругости \(K_S\) и модуля сдвига \(G\) и плотности \(\rho\).
Кстати, про плотность. Тут забавная история. Как нас учили на курсе общей геологии на первом курсе, скорость сейсмических волн зависит от плотности и увеличивается с её ростом. Нас так же продолжили учить и на кратком курсе геофизики на четвёртом курсе, но тут кто-то спросил “а как же тогда скорость сейсмических волн может увеличиваться с глубиной, если с глубиной плотность пород растёт, а она в знаменателе?”. И что я помню совершенно отчётливо – это как зам.декана Степанов П.Ю. был озадачен этим вопросом и отделался фразой “ну модули растут быстрее” без каких-либо пояснений. Искать ни на первом, ни на четвёртом курсе ответ времени у меня не было, но проблема в мозгу осталась. Теперь я могу дать ответ сам: потому что скорость сейсмических волн от плотности, как следует из формул, не зависит!
Применим Истинное Знание (термодинамику) к модулю упругости. Уравнение для него распишется как функция теплоёмкостей (\(C_p\), \(C_v\)), плотности \(\rho\), параметра Грюнайзена \(\gamma\) и коэффициента термического расширения \(\alpha\):
\[ K_S = \gamma \frac{C_p \rho}{\alpha} = \rho \frac{C_p - C_v}{\alpha^2 T} \frac{C_p}{C_v} \]
Остаётся вторая величина - модуль сдвига. Вот она уже не выводится из таких простых формул. Но. Можно сделать фокус. Есть немало работ (напр. [Ledbetter, 1977]), которые показывают, что для изотропных и около-изотропных сред (включая кристаллы кубической сингонии) есть определённое соотношение модуля упругости и модуля сдвига:
\[ G = K_S \left( 9\left(\frac{G}{E}\right) -3\right) = k K_S\]
где \( \left(\frac{G}{E}\right) \), отношение модуля сдвига (shear modulus) к модулю Юнга, является величиной почти постоянной. Теория предсказывает значение 0.375, из экспериментов следует 0.391.  Таки образом,
\[ G = 0.375 \div 0.519 K_S \]
Хорошая это оценка или плохая? Расчёт для мантии Земли по модели PREM даёт значения в пределах этого диапазона. Более того, мы знаем, что мантия Земли состоит из минералов с кубической или около того структурой. Мы можем использовать такие модели, если нас устраивает погрешность в первые проценты - десять процентов. В качестве первого приближения это прекрасный результат.

Итак, обоснование, почему скорость сейсмических волн напрямую не зависит от плотности (на самом деле, немножко сложнее, поскольку \(\alpha\) есть производная логарифма удельного объёма, но это уже высшие сферы):
\[V_p = \sqrt {\frac {1}{\rho} \left( K_S + \frac{4}{3} G \right)} = \sqrt {\frac {1}{\rho} \left( K_S + \frac{4}{3} kK_S \right)} = \sqrt {\frac {1}{\rho} \gamma \frac{C_p \rho}{\alpha} \left( 1 + \frac{4}{3} k \right)} \]
Или
\[V_p = \sqrt {\gamma \frac{C_p }{\alpha} \left( 1 + \frac{4}{3} k \right)}= \sqrt {\frac{C_p - C_v}{\alpha^2 T} \frac{C_p}{C_v}\left( 1 + \frac{4}{3} k \right)} \]
В итоговое выражение, как мы видим, плотность не входит. Вообще, забавно, что скорость сейсмических волн определяется теплоёмкостью, коэффициентом термического расширения и температурой.
Стоит отметить, что модуль объёмного сжатия убывает с ростом количества пор (в т.ч. заполненных жидкостью), например, [Russell, Smith, 2007]. Модуль сдвига падает по степенной функции до нуля примерно при 40% жидкости (напр. [Kováčik, 2008]). 
Модуль объёмного сжатия зависит от текстуры и структуры агрегата зёрен. По этому вопросу есть замечательная статья [Man, Huang, 2011], называющаяся "простая явная формула ...", от которой у меня начинала болеть голова примерно на третьей странице с суммами тензоров.
Теперь же о том, как всё это внедрялось. На этом графике рыжая линия - это то, что должно быть, синяя - это то, что получилось по старой методике расчёта, зелёная - по новой методике расчёта.


В чём разница?
В "голубой" модели я использовал простую старую простую параметризацию:
\[ V = (a_1T+a_2)P^2 + (b_1T + b_2)P + (c_1T+c_2) \]
\[ E = (d_1T+d_2)V^2 + (e_1T + e_2)V + (f_1T+f_2) \]
Как я уже писал в прошлых постах, эта схема неплохо описывала многие параметры для мантии типа её плотности и вполне вписывалась в допустимые диапазоны по теплоёмкости и коэффициенту термического расширения.
Почему она "сломалась" и стала рисовать эту странную голубую параболу?
Ответ довольно смешной. Потому что из того, что значения полинома описывают набор значений некоторого параметра, не следует того, что производная этого полинома будет описывать производную этого параметра.
Собственно, дифференцирование полинома дало параболу, зависящую от давления. Расчёт модуля упругости "усилил" эту зависимость, а окончательно она нарисовалась при расчёте скоростей сейсмических волн.
Когда я это понял, я поплакал в ванной, после чего потратил день и заменил параболы для объёма и энергии на уравнения Birch-Murnaghan (о них я немного писал, подробнее есть в той же википедии). И мы сразу получили "зелёный" вариант, который я считаю невероятно крутым результатом.  Особенно если учесть, что исходные коэффициенты мы берём из экспериментов над миллиметровыми количествами вещества или вообще в численных опытах на компьютере, а в итоге можем нарисовать реалистичную картину тех недр планет, которых мы никогда в жизни не видели и в руках не держали. Погрешность менее 5% при переходе через 9 порядков масштаба - это ничто.
Почему я сразу не сделал хорошо?
Интересный вопрос. Во-первых, в корректном решении новых уравнений есть много своих подводных камней, начиная с корректной итерационной схемы и кончая подходом Mie-Gruneisen (см. прошлый пост). Когда я начинал работать над моделью, я про них не знал. Это сейчас я потратил день. Тогда я бы потратил бы неделю-другую-третью, и не факт, что успешно. И это было правильное решение - начать с простого, перейти к сложному. Сам себя не похвалишь - весь день как оплёванный.

Мораль сей басен такова:
  1. Всё ложь. Если в уравнении [сейсмических волн] нарисована буковка [плотность], это ещё не значит, что она имеет значение.
  2. Всё ложь. Из того, что значения полинома описывают набор значений некоторого параметра, не следует то, что производная этого полинома будет описывать производную этого параметра.
    К слову, частично эта проблема решается сплайнами (см. напр. тут).
  3. Мы смело можем шагать из пробирки в планеты через 9 порядков масштаба, получая ошибки менее 5%. И это можно считать действительно потрясающим результатом.

2 комментария:

  1. Что тут сказать. Я, если ты помнишь, сейсмик, поэтому мне эти вещи довольно интересны (хотя у тебя сугубо теоретическая часть). Скажу лишь, что на мой (прикладной, без формул) взгляд скорость зависит от плотности лишь до определённой степени (так сказать, в реальных средах). В мантии зависимость больше выражается величиной давления и пластичности пород. На этом этапе плотность уже не играет роли, поскольку физически она уже больше и не повышается.
    Надеюсь, мой ответ не звучит по-дилетантски, повторюсь, это сугубо практический подход.
    Спасибо за пост!
    З.Ы. А почему не учитывается, что поперечных волн два типа?

    ОтветитьУдалить
    Ответы
    1. Ну да, в конечном счёте именно давление и оказывается переменной. От температуры мало что зависит, фазовый состав в первом приближении тоже постоянный.

      Два типа волн - хе-хе, это уже совсем высший пилотаж. Чтобы это рассчитать, надо знать, как минимум, размер зёрен и текстуру породы.
      Сейчас, кстати, есть статьи, где показана анизотропия как минимум для верхней мантии, и даже соотнесена с текстурами пород - то есть вдоль какой кристаллографической оси соориентированы зёрна минералов. Но всё это, в той или иной мере, пока +/- валенок.

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

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

      Удалить