Математическое моделирование относительной динамики твердого и жидкого ядер Земли Текст научной статьи по специальности «Математика»

УДК 517.958,550.31,551.16

Ю. М. Григорьев, О. Е. Скрябина

МАТЕМАТИЧЕСКОЕ МОДЕЛИРОВАНИЕ ОТНОСИТЕЛЬНОЙ ДИНАМИКИ ТВЕРДОГО И ЖИДКОГО ЯДЕР ЗЕМЛИ

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

Ключевые слова: математическая модель, приливная деформация, ядро Земли, дифференциальное вращение.

Введение. Изучение динамики внутренних масс Земли может осуществляться несколькими способами. Первый способ основан на обработке сейсмических данных. Современные представления о структуре Земли (кора, мантия, жидкое ядро, твердое ядро) получены именно таким образом [1]. Этим способом получено и наиболее значимое в динамике геосфер открытие последнего времени, заключающегося в том, что внутреннее ядро имеет скорость вращения несколько большую, чем Земля в целом, т. е. имеет место так называемый восточный дрейф ядра [2]. Второй способ изучения осуществляется с помощью лабораторных экспериментов. Так, в Институте горного дела СО РАН (Новосибирск) на плоской и трехмерной математических моделях показана справедливость гипотезы существования переноса внутренних масс для случая однородной Земли. Имеются экспериментальные данные о том, что в зависимости от параметров материала и области возможен как западный (скорость вращения ядра немного меньше скорости вращения Земли), так и восточный (скорость ядра немного больше скорости Земли) дрейф ядра [3].

Третий способ основан на построении математических моделей [4; 5]. Представляется, что все три указанных способа изучения динамики внутренних масс Земли взаимно дополняют друг друга.

В литературе наблюдается широкий разброс значений скорости дифференциального вращения ядра Земли, полученных из сейсмологических данных [6]: от 0,2 до 3 градусов в год с восточным направлением [7; 8], причем в последних работах наблюдается тенденция уменьшения этой величины [9; 10; 11]. Разработана математическая модель планеты, подверженной приливных деформациям, состоящей из двух твердых тел - ядра и оболочки, между которыми находится сферический слой вязкой несжимаемой жидкости [5]. Получено значение средней угловой скорости вращения ядра относительно оболочки в восточном направлении 0,38 минуты в год. Проведена оценка эффектов приливного переноса внутренних масс Земли без учета наличия внутреннего твердого ядра [4].

Постановка задачи. Ставилась задача оценить эффекты переноса масс в ядре Земли приливными волнами. Остальные существенные факторы процесса деформирования Земли (температура, неоднородность мантии и др.) не учитывались.

При построении математической модели представим Землю в виде тонкой оболочки, заполненной вязкой не-

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

|j.Av — Vp = 0, Vv = 0 , (1)

где v(r) - вектор скорости. р(г) - давление, М- - динамический коэффициент вязкости.

Применим следующий кинематический подход [12]: заранее зададим форму деформируемого тела и граничные скорости, имитирующие движение приливной волны. Для простоты будем считать, что возмущающая масса только одна (в данном случае ею является Луна), и она лежит в экваториальной плоскости Земли. При таких предположениях для упрощения используемых методов можно применять плоскую модель. Правомерность применения плоской модели объясняется тем, что из-за физических условий и симметрии задачи не должно происходить переноса вещества через экваториальную плоскость. Следовательно, явление можно промоделировать только на экваториальной плоскости, причем в условиях плоской деформации. Хотя плоские модели дают только качественные результаты, но для случая однородной Земли они совпадают с полученными экспериментальными лабораторными данными [4].

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

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

|j.Vv - Vp = 0, Vv = 0, vr = const, vn = 0.

На внутренней границе (окружности) задается постоянная касательная скорость, которая является заранее неизвестной угловой скоростью вращения внутреннего ядра Земли. Метод определения этой скорости будет рассмотрен ниже.

Вестник Сибирского государственного аэрокосмического университета имени академика М. Ф. Решетнева

Поставленная задача решается методом малого параметра. Естественным малым параметром является эксцентриситет эллипса. После решения краевой задачи можно оценить реально наблюдаемые эффекты.

Геометрия метода малого параметра. Для системы Стокса рассмотрим плоскую задачу в двусвязной области 5, ограниченной изнутри окружностью радиуса , снаружи - эллипсом Ь2. Уравнение внешней границы Ь2 запишем в полярных координатах г и & :

Ь * , & ,

™ ~ = 1 —-, а > Ъ.

s COS &

V,. = со^ = const, vn = О, Г Е Lx, V,. = v0 = const, vn = О, Г Е L2,

r = 2>V„O) = r0+^r, г = £Г

n=0 n=0

Если взять r0 = R2, to

r0 =^2,ri(9) = -r0sm2 a,

,2 on ,

l(3).

1

r2(&) = r0 sin &(1 + 3C0S &)...

Для угла справедливы соотношения

п dy . dx

COSÍk =—, sin =--------------.

1 ds ds

(8)

где (к - дифференциал дуги. Если уравнение границы записать в виде х = х(&) и у - у(&), то

У' • ~ х'

cos^ =

х'2 + у'2

sin 9, = -

х'2 + у’2

(9)

(3)

Граничное поле скоростей в принятой кинематической постановке имеет следующий вид:

(4)

где уг , уп - соответственно касательная и нормальная компоненты скорости, со - заранее неизвестная угловая скорость вращения внутреннего твердого ядра.

Рассмотрим часть области £ (рис. 1), в которой а -большая полуось эллипса, Ъ ~ малая полуось эллипса, - радиус внутренней границы, - внешний радиус,

равный большой оси; угол 0 является полярным углом; угол образован нормалью к внешнему эллипсу; &* = 3'1 - & - угол между нормалью и полярным радиусом г;т,п- орты касательной и нормали; пг, пд - орты полярной системы координат.

где штрих означает дифференцирование по &. Согласно (6) можно записать

х = (г0 + Xr) cos &, у = (r^+X^sinO, х' = ~(r0 + Xr)sinQ + Xr 'cos&, у' = (r0 + Xr) cos & + Xr ’ sin &.

Учитывая, что

cos&* = cos(^ -&) =

= cos cos & + sin sin &, sin&* = sin(^ -S) = s = in ^ sin & - eos ^ eos &,

из (11), (9), (10) получим

(10)

(11)

cos 9- =

sin & =

r0 +Xr

r0+Xr)2+(Xrf

Xr'

r0+Xr)2+(kr'f

Теперь можно выписать разложения

cos^ — 1 + ХС\ + X С2 +..., sin =1 + XS1+X2S2+...

где

с, = о,С2= --sin2 29.5, = sin 29, S2 = 4cos3 &sin&,....

(12)

(13)

(14)

Метод малого параметра. Суть метода малого параметра заключается в том, что решение задачи (2) ищем в следующем виде:

или

Проектируя граничную скорость v0 на внешней границе Ь2, которая направлена вдоль т, на направления пг и пд, напишем граничные условия (2.2) в полярных координатах:

Lr2 = ~v0 sin 9’, V9 ЦЙ2 = v0 cos 9’. (5)

Для решения первой краевой задачи для системы Стокса в области S с граничным полем скоростей (5) применим метод малого параметра. Естественным малым параметром X поставленной задачи является эксцентриситет эллипса s:

x=¿.

2

Следуя [13], разложим уравнение внешней границы в ряд по степеням X :

v(r) = ¿rv("4r),

v,(r, 9) = £rví"4r,9),

п=0

00

у9(г,9) = 2>Чй)(,-,Э).

(15)

Подставим (15) в граничное условие на внешней границе (эллипсе) (5) с учетом (6):

£ X"v(rn> (r0 + Хг, 9) = -v0 sin 9*,

п=О

00 __

X h"ví"} (го + ^, 9) = —v0 cos 9*.

(16)

В формуле (16) разложение правых частей получается с помощью (13). Теперь найдем разложения левых частей (16) в ряды по степеням X:

у*и)(г0+АЯ») =

(6)

= ví’V»,S) + *.

dv{:\r0+Xr,$)l

_дХ

V d2v(;\r0+Xr, 9)

(17)

+ ...

2 дХ2

Используя (6), найдем коэффициенты этого ряда:

ду(гп\г0+Хг,§)

дХ

( dv™(r,3) дгЛ

дг дк

(18)

dvfjr^)

дг

d2v{;\r0+'kr,&)

Л (9),

52v^(r,Q)f дг_

Ж

d2v

д2г (и)(г,9)

д2Х 2

^(3) +

dv(rn)(r,&) д2г

дг ¥х

dv(f(r,§)

2 г2(3).

= 0,vr=v0;v® =

а>/0)

= -v0Sl -r^S)

vr=->i(S)

dr dr

Отсюда получаем следующее:

„(я)

= 0, п = 0, оо,

І =С0І?1, « = °,

„(и)

= 0, п > 1.

ц(и +/у) = ср(г) -2у'(2) + \4/(г), (21)

где и и V - декартовы компоненты скорости. Такое представление аналогично формулам Колосова-Мусхелиш-вили из теории упругости.

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

Уу(0) =0, Vр- цДу(0) = О,

= о, 1 у

= 0,

1,(0)

(22)

1,(0)

д1г г=г° 1 дг

Подставляя в граничные условия (5) разложения (13) и (17) с учетом (18) и приравнивая выражения при одинаковых степенях А,, получаем систему «зацепляющихся» граничных условий для поправок . Таким образом, краевые условия на эллипсе сводятся к краевым условиям на границе г0 =К2. Ограничимся первым порядком приближения и выпишем эти условия:

■ \r=R2 ’ ■» \r=R2~

где v0 = const - граничная скорость, имитирующая движение приливной волны.

Решение задачи (22) известно [15] и имеет следующий вид:

vr =0,

-r2R2v0 + R2 Г(г2 -Rl)co + R2v0 ] (23)

VS“ r{Rl-R22) '

Задача (22) решена также с помощью комплексного представления (21), при этом комплексные потенциалы имеют вид:

(19)

cp(z) = iz

R2со -R2v0

В силу линейности задачи каждая поправка \ также удовлетворяет системе Стокса, а граничное условие для него вычисляется из предыдущих поправок и задается на окружности радиуса г = Л2. Таким образом, последовательно решая ряд краевых задач для системы Стокса внутри кругового кольца радиусом /?2 и внутренним мы можем найти решение в виде рядов (3.1). Когда параметр X достаточно мал, следует ожидать достаточно быструю сходимость этих рядов.

Запишем теперь краевые условия на внутренней границе (окружности):

=1>УИ)(ЗД = 0,

=Х^Чи)(^^) = со^.

ад2 -r22Y

^z) = -l3RL (R2®~Vo)-

Z J\ j — J\ 2

(24)

Если потенциалы (24) подставить в формулу (21), то после перехода к полярным координатам, получим поле скоростей вида (23).

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

Уу(1) = 0,

Ур - |аДу(1) = О,

V® U = o,

,0)

U = o,

U* = ~v0 sin 23,

(25)

2 n v0(R2 +R2)-2R2R2g)

К -Щ

(20)

Таким образом, для решения исходной задачи с граничным полем скоростей вида (5) получен ряд краевых задач для кругового кольца с внутренним радиусом Rl , внешним - Р2 и с краевыми условиями (19) и (20).

Решение краевых задач методом комплексных функций. Для решения краевых задач для системы Стокса внутри кругового кольца, возникающих в методе малого параметра, применим метод комплексных функций. В [14] показано, что в плоском случае кватернионное представление общего решения системы Стокса переходит в комплексное, выраженное через две аналитические функции <р(г) и \|/(г) комплексного переменного г = х + ¡у:

Введенные в (21) аналитические функции <р(г) и \|/(г) ищутся в виде разложений Лорана с учетом их многозначности в двусвязной области. Коэффициенты ряда находят из систем линейных алгебраических уравнений, получаемых из краевых условий. В первом порядке приближения комплексные потенциалы получены в следующем виде:

1 з

(р(г) = а 1— + а1г + а3г ,

| ! (26)

у (г) = Ь_3— + Ь_1— + Ь1г, г г

где ах, а , аъ, Ь_ъ, /> , />, - выражены через коэффициенты разложения в комплексный ряд Фурье краевых условий из (25).

Полученные комплексные потенциалы <р(г) и \|/(г) подставляем в выражение (21), затем, отделяя реальную и мнимую части полученного выражения, получаем декартовы компоненты искомого поля скоростей V®, В ко-

рисунков нету-файл испорчен

Вестник Сибирского государственного аэрокосмического университета имени академика М. Ф. Решетнева

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

Метод определения угловой скорости движения внутреннего ядра. Угловая скорость движения твердого внутреннего ядра определяется из условия стационарности движения нашей системы, которая заключается в следующем: суммарный момент всех сил М, приложенных к ядру, должен быть равен нулю. Это условие выразим через тензор вязких напряжений следующим образом:

М = (^ г х (стгапа)<й = Щог$с!1 = 0 . (27)

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

<7,!» = 2цб,.9 = 1 dvr(<D)+_dv9(o>)

(со)

= 0,

да волна подходит к нему, он совершает вертикальное колебание. Один оборот внешней точки в кинематической постановке соответствует одному обороту приливной волны. Картина, видимая наблюдателю, будет более наглядной, если следить за концами траекторий частиц, которые в начальный момент времени (в кинематической постановке) находились на одном радиусе. Реальный наблюдатель, который в начальный момент времени находился на вершине приливной волны, увидеть эффект переноса внутренних масс под действием приливной волны после ее оборотов (рис. 3).

(28)

(29)

д& дг

Подставляя уравнение (28) в (27) и интегрируя по полярной координате d&, мы получаем уравнение с одной неизвестной щ:

2nRR?

---9--

f [со0 (^ - 2) - 2а>(^ -1>] ^ (

|^+Т?2—2со + соq (X- + 2) j где со0 -угловая скорость вращения внешней границы.

Из (29), как и следовало ожидать, в нулевом приближении угловая скорость ядра равна внешней угловой скорости со0, т. к. имеем задачу о вращении кругового кольца как жесткого целого, т. е.

— ^о/ ^2’

где v0 = const - скорость на внешней границе (эллипсе), R\ - внешний радиус кольца, в котором решаются краевые задачи (большая полуось эллипса).

Решая алгебраическое уравнение (29), находим угловую скорость внутреннего ядра в виде:

Д2а-2) + Д?(А, + 2)

со = coQ lV 7 , -----(30)

2(Д22+Д,2(?1-1))

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

Приведем траектории частиц жидкости после одного полного оборота внешней частицы X = 0,08; а = 3, Ъ = 2,598 (рис. 2); покажем также положения частиц, находившихся на большой полуоси у = 0 в момент времени t = 0, после 30 оборотов внешней частицы Х = 0,08;а = 3, Ь = 2,598 (рис. 3).

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

Рис. 3

Оценим скорость вращения внутреннего твердого ядра Земли. Для этого разложим выражение (30) по малому параметру X (до первого порядка):

со = coJ 1 +

(31)

Высота приливной волны для реальной Земли составляет 0,5 м, но так как областью нашего исследования является жидкое ядро с твердым внутренним ядром, берем высоту порядка 0,2...0,25 м. Тогда малый параметр X получается порядка Ю“7 • Согласно формуле (31) получаем, что опережение внутреннего твердого ядра составляет 0,39 минуты в год. Эта цифра хорошо согласуется с оценкой из [5], полученной другим, более сложным методом.

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

Библиографический список

1. Ботт, М. Внутреннее строение Земли / М. Ботт. М.: Мир, 1974.

2. Song,X. Seismological evidence for differential rotation of the Earth’s inner core / X Song, P. G. Richards // Nature. 1996. Vol. 382. P. 221-224.

3. Ревуженко, А. Ф. О течении сыпучей среды с возможными неограниченными скольжениями по поверхности локализации / А. Ф. Ревуженко, А. П. Бобряков,

B. П. Косых//Физико-техн. пробл. разраб. полез, ископаемых. 1997. № 3. С. 37-12.

4. Григорьев, Ю. М. Плоская задача о переносе внутренних масс приливными волнами / Ю. М. Поляков // Мат. заметки ЯГУ 1999. Т. 6. Вып. 2. С. 9-20.

5. В ильке, В. Г. Об относительном движении ядра и оболочки планеты в гравитационном поле точечной массы / В. Г. Вильке // Прикл. математика и механика. 2006. Т. 70. Вып. 4. С. 617-630.

6. Кузнецов, В. В. Анизотропия и супервращение внутреннего ядра Земли / В. В. Кузнецов // Науки о Земле. Современ. пробл. сейсмологии : сб.ст. М., 2001.

C. 68-83.

7. Овчинников, В. М. О скорости относительного вращения внутреннего ядра Земли / В. М. Овчинников,

B. В. Адушкин, В. А. Ап // Докл. АН. 1998. Т. 362.

C. 1000-1003.

8. Laske, G. Limits on differential rotation of the inner core from an analysis of the Earth’s free oscillations / G. Laske, G. Masters//Nature. 1999.402. №. 6757. P. 66-69.

9. Souriau, A. Inner core rotation: a test of the wordwide scale / A. Souriau, G. Poupinet // Physics of the Earth and Planetary Interiors. 2000. Vol. 118. P. 13-27.

10. Song, X. Support for differential inner core superrotation from earth-quakes in Alaska recorded at South Pole station / X. Song, A. Li // Journal of geophysical research. 2000. Vol. 105.B1.P. 623-630.

11. Xu, X. Evidence for inner core super-rotation from time-dependent differential PKP traveltimes observed at Beijing Seismic Network / X. Xu, X. Song // Geophys. J. Int. 2003.152. №3. P. 509-514.

12. Ревуженко, А. Ф. Механика упруго-пластических сред и нестандартный анализ / А. Ф. Руженко ; Новосиб. гос. ун-т. Новосибирск, 2000.

13. Григорьев, Ю. М. Пространственная задача о переносе масс приливными волнами / Ю. М. Григорьев, А. Ф. Ревуженко //Вычисл. технологии. 2000. Т. 4.5. С. 40-54.

14. Валландер, С. В. Лекции по гидроаэромеханики /

С. В. Валландер//Л.: Изд-во Ленинград, гос. ун-т, 1978.

15. Григорьев, Ю. М. Математические проблемы моделирования направленного переноса внутренних масс Земли приливными деформациями / Ю. М. Григорьев, О. Е. Скрябина // Динамика сплошной среды. Вып. 122. Новосибирск, 2004. С. 57-62.

Y. M. Grigoryev, O. E. Skryabina

MATHEMATICAL MODELING OF RELATIVE DYNAMICS OF THE EARTH LIQUID AND SOLID INNER CORES

On the base of the kinematical approach two-dimensional mathematical model of the Earth tidal deformations taking into account inner solid core is presented. A method of the least parameter for solving boundary value problems in this model is developed. Solving boundary value problems for the Stokes system in a biconnected domain a complex function method is used. The model describes the effect of the differential rotation of the Earth inner core.

Keywords: mathematical modeling, tidal deformations, inner core, differential rotation.

УДК004.052

А. В. Котенок

К ПРОБЛЕМЕ ВЫБОРА МУЛЬТИВЕРСНОННОЙ МОДЕЛИ ПОВЫШЕНИЯ НАДЕЖНОСТИ ПРОГРАММНОГО ОБЕСПЕЧЕНИЯ

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

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

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