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

УДК 533.6.011+551.2+519.6

В. В. Шайдуров, Г. И. Щепановская ГАЗОДИНАМИЧЕСКАЯ МОДЕЛЬ ВНУТРЕННЕГО СТРОЕНИЯ ЗЕМЛИ1

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

Внутреннее строение Земли оценивается по известной массе, моменту инерции земного шара и на основе изучения упругих волн от землетрясений. Известно, что плотность вещества в центре Земли р> 12,2 г/см3 и ядро Земли на глубине 2 900 км отделено от лежащих выше слоев резким скачком плотности - порядка 4 г/см3. Таким образом, упругие свойства внутри Земли изменяются скачкообразно на некоторых определенных значениях глубин и плавно - в пределах слоев, разделенных этими границами. Важнейшими границами являются поверхность Мохоровичича, залегающая на глубине 10...70 км, и поверхность Вихерта-Гутенберга на глубине 2 900 км, резко преломляющая продольные упругие волны и не пропускающая поперечных волн. Эти границы разделяют земной шар на три главные зоны: кору, мантию и ядро. Земля образована твердым слоем толщиной примерно равной половине радиуса, и ядром, имеющим два агрегатных состояния - жидкое и твердое. Кора обладает наибольшей жесткостью, мантия характеризуется высокой вязкостью, а ядро находится в состоянии, близком к жидкому, и реагирует лишь на продольные волны изменением объема. Внутри трех главных зон земного шара имеются и менее четко выраженные границы. В последнее время в некоторых работах вводится предположение, что внутреннее вещество Земли ведет себя по отношению к длительно действующим силам как жидкость.

Авторами предлагается математическая модель, позволяющая адаптировать структуру Земли газодинамическим состоянием плотности, динамической вязкости, внутренней энергии и высоких температур. Динамика геосферы Земли исследуется в рамках модели вязкого сжимаемого газа, когда плотность среды меняется во времени и пространстве. Сферически-симметричное течение вязкого теплопроводного газа рассматривается с учетом гравитационных сил. Уравнения решаются в безразмерной постановке [2], когда линейные размеры отнесены к радиусу Земли, а все газодинамические величины -к соответствующим характерным значениям на поверхности земного шара.

Для описания процессов сферически-симметричного расширения, сжатия, разогревания и охлаждения Земли используются уравнения Навье-Стокса. Введем сферическую систему координат (r, 0, ф), связанную с декартовыми координатами следующим образом: x = r cos ф sin 0, y = r sin ф sin 0, z = r cos 0, 0 < r < + ^, 0 < 0 < n, 0 < ф < 2 п. После записи уравнений Навье-Стокса в сферической системе координат [2] учет сферической симметрии приводит к равенству нулю производных по ф и 0 и компо-

нент скоростей по ф и 0. В итоге уравнения неразрывности, импульса энергии записываются в следующем безразмерном виде:

до 1

д . 2 ч 1 до

+-------------------(r о u) + — u — = 0.

dt 2r dr 2 dr

o(t, r) = Vp(t, r),

д / ч 1 d . 2 2.

— (pu) + ^ — (r pu ) = dt r dr

= -дР + 1_ d (r ^ ) + Fm (r, p),

dr

d

dr

+P

(pe) + ^ (r 2upe) +

dt r dr

1 d 24 1 d 2 4

— — (r u ) = —(r qr) + x r dr r dr

du

dr

(1)

(2)

(3)

Искомыми функциями являются плотность р = о2, скорость и и внутренняя энергия единицы массы в, через которые выражаются давление Р и динамический коэффициент вязкости ц из уравнений состояния, и Еет (г, р) - потенциал гравитационных сил. Уравнения состояния представляют собой сложную алгебраическую зависимость плотности от давления и температуры с учетом фазовых переходов и метаморфизма основных составляющих веществ. Отметим, что запись уравнения (1) для о, а не для р переводит закон сохранения массы из нормы пространства Ь1 в норму гильбертового пространства Ь2. Впоследствии это значительно упростит обоснование устойчивости и сходимости.

Потоковые соотношения, связывающие тензор напряжения тгг с тензором скоростей деформации и определяющие тепловой поток qr через безразмерные параметры - число Рейнольдса Яе и число Прандтля Рг, выражаются зависимостями

4 ди 4

тгг =-----ц------------ци,

ЗЯе дг 3г Яе

. (4)

У

qr =

de и—. Re Pr dr

Потенциал гравитационных сил в данном случае представлен в следующем виде:

Р г

=-4nFr ч

l 2 r2

Jpx2 dx,

(5)

где Гг - число Фруда; х - безразмерная величина, которая определяется гравитационной постоянной, ускорением свободного падения и плотностью на поверхности Земли.

1 Работа выполнена при финансовой поддержке программы Президиума РАН № 16 (проект № 9).

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

Г 2 до 7

\г -*■

1

2 Ф/ ом

1і

Г1 2 дф/

. - I — г ом—'-ёг г = Ь дг

ти (см. таблицу) [1].

На расчетной области радиус Земли занимает одну четвертую часть, остальные три четвертых приходятся на атмосферу. Введем равномерную сетку гі = (і +1/2)И, і = -1, п, с шагом И = 2/(2п +1).Обозначим ОИ = {гі, і = -1, ..., п} и введем множество внутренних узлов і ОИ = {гі, і = 0,1, ..., п -1}. В результате Подставим приближенное решение (7) в равенство (11) расчетная область □ будет разбита на п + 1 интервал вместо °

11 2

+| 2 г Ф/м

до

дг

ёг = 0, / = 0, ..., п.

(11)

(Г, Г+1), г = -1, 0, ..., п-1. Для каждого узла г е Ол определим базисную функцию ф., которая равна единице в г., нулю во всех остальных узлах Оа и линейна на каждом отрезке (г, г+1) (рис. 1):

Фі(г)=

г - г

И

ес^и г є [г,.-, г ],

ес^и г є (г, г+1 ],

0, ес^и г й [г,.-, г+1 ].

I- 0 4

до, (і) 2 12 дф‘ 12 дфі

—д^~г Ф/Фі- ^г ом ^ Ф'+2г оімФ‘

+2(м°і- Ф/ Фі)

= 0, / = 0, ..., п.

ёг +

(12)

(6)

С учетом выбора базисных функций (6) отличными от нуля в (12) остаются следующие слагаемые:

I рО; (/) 2 12 дФ‘ 12 дФ(

II д_г Ф/Ф --г ои-ТФ +2г ОмФ/ "дг"

ёг -

п

Ои(і,г) = £о, (і)Фі (г)

Я =

до

1 д

ді 2г2 дг

(г оим) + — м

1 до,

2 дг

+ 2(мО-Ф/Ф-) I=1

= 0, / = 0,..., п.

(13)

Для вычисления интегралов в уравнении (13) применим квадратурную формулу трапеций и равенства 0-1 = О0, и_1 = -и0, вытекающие из симметрии задачи. Аппроксимируя производную по времени правой разностной производной, в итоге получим разностные уравнения:

/1 1 \

И

Рис. 1. Вид базисной функции ф.

С помощью введенных обозначений сформулируем метод Бубнова-Галеркина. Будем искать приближение Ои (/, г) к функции о(/, г) в виде

т

г 2 о1+1 + '0 и0 ^

1 9 1 1 9 1

4 г^! + 4 г0м0

V У

= — ок0, I = 0,

4т 0

+1

°1 =

(7)

с некоторыми коэффициентами о( (і), зависящими от времени. Подставим пробное решение оИ в уравнение (1) вместо точного решения и умножим скалярно на Ф^ (Я, Ф‘) = 0, / = 0, ..., п (8)

где

4

1

'/-1 +- г,Ио,

— г}+, и1+| +— гім]

о1+1 = — г Ио

и/+1 Ч пи/ -

I = 1, ..., п -1,

-1 г2 мк 1 -1 г2мк

4 п-1 п-1 4 п п

V /

+1 , оп-1 +

(9)

И 2 1 2 1

^ гп +0 гпМп ч2т 2 У

о1+1 =

И

2т

/ = п.

(14)

Здесь в нелинейных членах заморожены коэффициенты,

Используя редукцию скалярного произведения из К3 в известные с предыдущего шага по времени.

.К1 при условии сферической симметрии, получим следующее выражение:

Аналогично (7) будем искать функцию ик (/, г), аппроксимирующую функцию и(/, г) в виде

Iг 2 Ф/

до 1 д

1 до

----1--- —(г ом) + — м —

ёг = 0, / = 0,..., п. (10)

і (і, г) = Х м/ (і) Фі (г)

(15)

Значение плотности, коэффициента динамической вязкости, давления и температуры внутреннего вещества Земли в размерных величинах

і=0

И

і=0

м

і=0

к, км 15 60 100 150 200 300 400 600 1 100 1 600 2 700 2 870 2 900 4 700 6 371

Р, кбар 3,5 14 29 46,5 64,9 99 130 190 30 530 1 240 1 340 1 350 3 040 3 632

Т, К 673 873 1 103 1 343 1 473 1 623 1 673 1 723 1 873 2 700 3 500 4 000

с 2,85 3,34 3,37 3,37 3,36 3,48 3,54 4,13 4,74 5,03 5,55 5,68 9,89 12,3 13,01

м 0,35 0,72 0,72 0,72 0,7 0,75 0,82 1,3 1,7 2,2 2,9 2,95 0 0

с некоторыми коэффициентами и1 (/), зависящими от времени. Вновь, используя редукцию скалярного произведения из К3 в К1, в случае сферической симметрии получим следующее выражение:

+

1 4 дФі 4

- Рммі-Ф/ Фі - Ф/ — ци, — + Ф/ — ц мі Фі

2 3 Яе дг 3 Яе

V 1

дм м Эр

) — + -

ді 2 ді

1 д 2 2ч м Э 2 Эр + —^г(г рм2)--—^—(г2рм)+ -— г Эг 2г Эг Эг

= 1 Р

2гф/ + г

дФ/

дг

ёг - (Ф/Р) |г=1 + |г 'фЛ,^,

/ = 0, ..., п.

(19)

_1 _д_

г2 дг

дм 4

-г* ц--------------------------гц м

-

ёг = 0,

(16)

3Яе дг 3Яе

\

/ = 0, ..., и.

Отметим, что в (16) использовано исходное уравне- образом ние (2), из которого вычтено уравнение неразрывности, умноженное на и / 2. Это приведет и в непрерывном, и в дискретном виде к устойчивости в норме пространства Ь2, а не в Ь1, как в уравнении (2). Применяя к (16) формулу интегрирования по частям, получим

Для вычисления интегралов в уравнении (19) применим квадратурную формулу трапеций и равенства м-1 = -м0, ц-1 = ц0, вытекающие из симметрии задачи. Аппроксимируем производную по времени следующим

<м + м =

ді 2 ді

/~\1+1*» 1+1 I +1 /

Р; мі -\]Рі \Р/м

т

|г2 Ф/ |р^ +

дм м Эр

ёг + (ф/ р м2)

В итоге получим

ді 2 ді

У

= 1 -|г2рм2 Эф^ёг - 1(Ф/рм2) 0 дг 2

1 1 д

+1—г 2р м — (ф/м )ёг + (г2 ф/Р)

2 дг

^И 2 2 л

— г 2р1+1 +--------г 2ц1 +------------г 2ц*

0^0 , Т! 0^0 Т Т! 1^1

т И Яе 3И Яе

м„+1 +

2 2 1 2 1 12 1+1 1 2

-----г12ц1 + — ^ +Тг02р0 +1мо -ТГ^-

4 3 И Яе

-г12р11 X 1

4 1 1 3И Яе 1 3Яе

г 2ц1

'0

-I Р

2г ф/ + г2 Эф^ Т/ дг

/р° >/р0 м±Иг2 +1 _2р* + 2ИгР -1 Р^г2 + Иг„2Я*„, / = 0,

ёг -

4 дм

ф/^~ ц^~-ф/^~ цм

3Яе дг 3Яе

1 1г ^ 9 ІГ ^ ІГ 7 ІГ 4-1 1г ^

—г 1 р1 їм і--------------------г\ц? 1-----------------------г Х1 — г рГХ----------

4 /-1Км /-1 3ИЯе м ^ 3Яе /-1 м 4 11 1 3Ие

г 2ц1 —1+ +

Г дф/ ( 4 2 дм 4

+ 1 I г2ц~-:^^гцм

дг I 3Яе дг 3Яе

ёг -

І I И 2^ 1+1 . 2 ? 1 1 1

+ 1 —г, р, +--------------------г 1ц/ 1 +----------------г ці +------------------г,+1ц,+1 —

1 т 11 3И Яе / -1 /-1 3И Яе 11 3И Яе / + 1+1 I '

1 2-1+1 1 2 ^, 1 2 ..1 . 1 2.1+1 1 2 1 I 1+1

4гіег+1ц+1 + ііег+1ц-1 + 4гр “/-ІИ| =

1

г2ф1ретёг = 0, / ^ ..., п.

(17)

л/рЛ/р[

Подставим приближенное решение (15) в равенство (17) вместо м:

1

~Иг? +1 г/2-1Р/-1 + 2Иг1Р/1 -

г Ф/ Фі

рдмі + м^ Эр ді 2 ді

-1 г 2р мм. ф. +

2 * дг

— г2 Р+ + Иг-2., / = 1, ..., п -1,

2 / +1 /+1 і ет,і ’ ’ ’

-1 г21рп+1і —— гп іцп і —~ гп іцп і -1 г2рп+ім; —— г,2цп 1+ +

4 ^п-1 1 3И Яе п-і 3Яе п-1 4 п 3АЯе I -1

1 2 дфі. 4 2 дф/ дфі

+ — г омм —^ф/ +----------г цм. ————

2 р і дг/ 3Яе ^ дг дг

+1 _г2р*+1 г21ц* 1-----— г2ц*

1 2т 3И Яе ”- ”-і 3И Яе ” 3Яе

гТцп +^г,цп ^г,2рпX - + =

4 дФ/

------гц м. т,-?1-

3Яе ™ дг

ёг +

УР^л/рп Иг 2 + 1 г 2 р1 +

2т Иг” + 2 п-іРп-1

1 4 дфі 4

-рммі Ф/Фі-Ф/ — цм — + Ф/ — цмі Фі

2 3Яе дг 3Яе

1

(20)

2гф/ + г2 — Т/ дг

ёг - (ф/Р)

1

=1 + |г ^Ф/^ет ё^, (18)

Будем искать приближение еи (і, г) к функции е(і, г) в виде

і (і, г) = Х е (і)фі (г)

_ . (21)

/ = 0, ..., и.

С учетом выбора базисных функций (6) отличными от с некоторыми коэффициентами е.(), зависящими от

времени.

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

нуля в (18) остаются следующие слагаемые:

Эм1 мі Эр 1 1 2 Эф; 1 2 дфі

ґф,Ф I р—L + — — \-~г рмм ф —— + — г рмм —— Ф + ді 2 ді 2 Н ^ дг 2 У 1 дг ^

дФ/ дФі 4

дФ/

г цму ————-гц мі Фі. ——

3Яе ' дг дг 3Яе

дг

ёг +

і ( д д д

|1 г >/— (ре) + Ф Р эг (г2м) + Ф/ дг (г 2мре)-

Яе Рг дг I дг

+

і=0

3Яе

г Ф/ ц

ґди * Чдгу

4 дм

— г Ф/ ц м э-3Яе дг

ёг = 0,

(22) / = 0, ..., п.

Применяя формулу интегрирования по частям к (22), получим

ГІ 2 д , ^ „ 2 „дм 2 дф/ 7 2

11 Г Ф/ ді (ре) + 2гф/мР + г Ф/Р "д--г "ре-г + ЯР г ц

,дФ/ , У 2М де дФ/

д/

3Яе

г Ф/ ц

дм

Чдгу

дг

4 дм

+—гФ/ цмт-3Яе дг

дг ЯеРг дг дг

ёг + (мреф/)

' у деЛ ЯеРг дг

V У

= 0, / = 0, ..., п.

(23)

Подставим в равенство (23) вместо е приближенное решение (21):

м р Ф/Є Фі

Эф

дг

-п ГІ 2 ^ Ч 2 Эф/ У 2 дфі дф/ 1

£ IIг ффіді(реі)-г м р ^ фі +ярцг еі -дТ:-д: \ёг+

1. Яе Рг'

дфі Эф дг дг

У дфі

---!—цф/е,——

ЯеРг дг

3Яе

г°ц°м°1+1(м11+1 - м1+1), / = 0,

--г2 мk+1рk+1____________г-____г2 и1 -

Г/-и/-^ /-і Г/-1^/-1

2 2И Яе Рг

у г2 ц1

-г 2р1 +1 + т / / 2И Яе Рг

7 г- цГ-і +т^р- г2 ц1 +

2И Яе Рг

У

е1+1 +

И ЯеРг

2И Яе Рг

г+1 ц;+і

1

У

--------------г2 ц, +1 г+1мlk++11рk++11 +-------------------

2И ЯеРг / п 2 т тт 2И Яе Рг

г+1 ц1+1

е1+1 = е!+1

= Иг2р1е1 - 2Иг1м1+1 Р -1 г/Р1 (м^1 - мЦ) +

1

3И Яе 2

г2ц 1 (мк1+1 - м*+1)2 -

3Яе

г ц 1 м1 +1(м^1 - мЦ),

/ = 1, ..., п -1, У

1 гі-мЖ1-2, Я Р г„2-іцп-і-И Я Р гХ

2 2И ЯеРг 2И ЯеРг

е1+1 +

еи-1 т

I ( 2гФ/мР+г 2ф,р *Эг-

дм

3Яег “ф/ ц1дг

2

4 дм

+—г Ф/ цм1Г 3Яе дг

ёг,

м р Ф/ЄФі

У дфі

------цфе —-

ЯеРг Т/ ‘ дг

=- (2гФмР+г 2фр 'дг ~

3Яе

г ф/ ц

дм дг

4 дм

- — г ф/ ц м — 3Яе дг

(25)

(24) / = 0,..., п.

С учетом выбора базисных функций (6) отличными от нуля в (24) остаются следующие слагаемые:

Г Г Г 2 ^ ч 2 дф/ У 2 дф. дф/ V

5, IIг фф т>,,Ре 1 -г и^-э;геф +яергц е ЭГ 1Г Г+

— гУ+1 +---------¥--г21ц1 1 +1 ггм1+1р1+І---------У-----г2ц1

2т п 2И Яе Рг "- 2 п п Гп 2И ЯеРг п

= И гії/, - Иг,м‘,+Р1 - 2 г; Р1 (м;+1 - м“)+ 2 гХ (м;•' - м;+1)- -

3Яе

3И Яе

г„ цХ+1(м:+1 - мП), / = Й.

(26)

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

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

Для вычисления интегралов в уравнении (25) приме- намический коэффициент вязкости, а также температура ним квадратурную формулу трапеций и равенства отнесены к с°°тветствующим характерным величинам е-1 = е0, ц-1 = ц0, вытекающие из симметрии задачи. на поверхности земного шара.

Аппроксимируя производную по времени правой раз-

Представим результаты вычислительного эксперимен-

ностной производной, в итоге получим разностные урав- та в виде графиков (рис. 2, 3). Эти графики показывают, нения: что со временем наша планета сжимается под действием

гравитационных сил и уменьшается в радиусе, плотность

И

1 + 31 г 2ц1 +_________1 г 2ц1

2И ЯеРг 0 ц0 2И Яе Рг 1 ці

Є1+1 +

1

------------г02ц0 +- г2мk+1р1+1 +

2И Яе Рг 0 2 1 1 Р1 2И Яе Рг

Є1+1 =

= И г02р0Є01 - 2Иг,м\+1Р01 -1 г2Р1 (мі1+1 -и°k+1) + т 2

+ Го2Ц°(Mlk+1 - и°°+1)2 -

3И Яе

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

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

Рис. 2. Распределение плотности по радиусу Земли для различных моментов времени под действием гравитационных сил (выделена линия начального распределения плотности, соответствующая важнейшим граничным поверхностям Мохоровичича и Вихерта-Гутенберга)

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

1. Добрецов, Н. Л. Глубинная геодинамика / Н. Л. Доб-рецов, А. Г. Кирдяшкин, А. А. Кирдяшкин. Новосибирск : Гео. Филиал изд-ва Сиб. отд-ния Рос. акад. наук, 2001.

2. Ушакова, О. А. Метод конечных элементов для уравнений Навье-Стокса в сферической системе координат / О. А. Ушакова, В. В. Шайдуров, Г. И. Щепановская // Вестн. Краснояр. гос. ун-та. 2006. № 4. С. 151-156.

Р

0,05 0,09 о,13 0,17 0,21 0,25 0,29 г

Рис. 3. Распределение безразмерной температуры по радиусу Земли для тех же моментов времени, что и на рис. 2, под действием гравитационных сил (выделена линия начального распределения температур)

V. V. Shaidurov, G. I. Shchepanovskaya GASEOUS MODEL OF THE INTERIOR STRUKTURE OF EARTH

In present paper computing model allowing considering geodynamic process of expansion, compression, heating and cooling of the Earth is suggested. Dynamic of geosphere is investigated in the context of viscous heat-conducting coercible gas when density and viscosity of medium depend on time and coordinates. Suggested model allows considering not only crust and mantle of Earth but also internal structure including Earth core.