Метод конечных элементов в моделировании сферически-симметричных пульсаций Земли Текст научной статьи по специальности «Физика»

УДК 519.6:533.6:551.2

Метод конечных элементов в моделировании сферически-симметричных пульсаций Земли

Александр В.Вяткин*

Институт фундаментальной подготовки, Сибирский федеральный университет,

пр. Свободный 79, Красноярск, 660041, Россия

Владимир В.Шайдуров^

Институт математики, Сибирский федеральный университет,

пр. Свободный 79, Красноярск, 660041, Россия

Галина И.ЩЩепановская*

Институт вычислительного моделирования СО РАН, Академгородок, 660036, Красноярск, Россия

Получена 15.03.2008, окончательный вариант 15.05.2008, принята 05.06.2008 Рассматриваются вопросы математического и численного моделирования геодинамических процессов 'расширения, сжатия, разогревания и охлаждения Земли.

Ключевые слова: математическое и численное моделирование, геодинамические процессы

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

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

* e-mail: gi@icm.krasn.ru t e-mail: shidurov@icm.krasn.ru te-mail: gi@icm.krasn.ru

© Siberian Federal University. All rights reserved

Предложенная в данной работе модель позволяет рассмотреть не только кору и мантию Земли, но её полную внутреннюю структуру, включая ядро. Модель является сферически-симметричной, поэтому при задании свойств соответствующих слоёв Земли и атмосферы строение Земли и атмосферы в различных направлениях радиус-вектора, исходящего из центра Земли, считается одинаковым. Задача рассматривается в нестационарной постановке. Уравнения решаются в безразмерной постановке [3], [4], когда линейные размеры отнесены к радиусу Земли, а все газодинамические величины — к соответствующим характерным значениям на поверхности земного шара.

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

да 1 д , 2 „ 1 да . . /—---. .

ж + 2Т2дТ + 2udr = °, a(tr) = ^^ (1)

д 1 д дР 1 д

+ — тг(r2pu2) = + (r2 тгг)+ Fem(r,p), (2)

дг r2 дr дr r2 дr

д , 1 д 2 n1 д 2 n 1 д 2 дм

дг (pe) + r2 зт(r uPe)+дг(r u) = - r2 дг(r qr)+Trr дГ- (3)

Искомыми функциями являются плотность р, скорость м и внутренняя энергия единицы массы — e, через которые выражаются давление Р и динамический коэффициент вязкости у из уравнений состояния, Fem(r, р) — потенциал гравитационных сил. Уравнения состояния представляют собой сложную алгебраическую зависимость плотности от давления и температуры с учётом фазовых переходов и метаморфизма основных составляющих веществ.

Потоковые соотношения, определяющие компоненту тензора напряжений тгг и тепловой поток qr через безразмерные параметры — число Рейнольдса Re и число Прандтля Pr, — выражаются зависимостями [5], [6]:

4 du 4 y de

Trr = 3Re ^dr- sTrM u qr = - Repr ^ (4)

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

r

Fem = - 4п Fr XI • р J (5)

0

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

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

На расчётной области радиус Земли занимает одну четвёртую часть. Введём равномерную сетку ri = (i + 1/2)h, i = —l,...,n, с шагом h = 2/(2n + 1). Обозначим

11 н = г = — 1, ..., п} и введём множество внутренних узлов Пн = г = 0,1,... ,п — 1}. В результате расчётная область О. разбилась на п +1 интервал (г^ ^+1), г = —1,0,... ,п — 1. Для каждого узла г^ € 11 н определим базисную функцию которая равна единице в г^, нулю во всех остальных узлах 11 н и линейна на каждом отрезке (г^ г^+х) (рис. 1).

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

г — г 1-1

если г € [гi-l, г];

^®(г) 1 —-, если г € (гi, гi+l]; Н

0,

если г € К-ь гi+l].

(6)

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

= ^2 а^)уч(г)

(7)

i=0

с некоторыми коэффициентами оч(£), зависящими от времени. Подставим пробное решение ан в уравнение (1) вместо точного решения и умножим скалярно на у;:

где

(Д, у;) = 0, 1 = 0, ..., п, дан 1 д 2 N . 1 дан

Д = ~дГ + 2Т2дг(г ани) + 2

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

1

2

0

да

1д

1 да

У ' + 2г2 дг(г а и) + 2идг ) йг = 0 1 = 0 ..., п.

2 дг

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

1 1 /г2у; ¿>?г+(2аЧ=1—/

ду;

1

1 2 да

и — —г а и—— йг + —г уци^—йг = 0, г=1 .] 2 дг ] 2 дг

00

(8) (9)

(10) (11)

Н

0

l = 0, ..., n.

Подставляя приближённое решение (7) в равенство (11) вместо а, получим

Е

/da¿(t) 2 1 2 1 2 ,

(v_^ - 2r + 2Г ' +

1

r = 1J

0,

l = 0, ..., n.

(12)

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

l+1 Е

¿=г-1

2 1 2 dtPl . 1 2 д^М J ,

(v_^ - 2r ai« + 2r ' +

r=U

0, l = 0, ...,

(13)

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

hr2afc+1 + f1 r2ufc + 1 Л,« |аЛ Tro ао + r1 «1 + 4 ro «о J а1

h

2„ k r 0 а0 ,

l = 0,

-1 r2 «k _k+1 + ^ r2hak+1 + /1rl_1Ml_1 лГ1 «l )al_1 + rl hal +

+ I4 r'2+1«k+1

+ 1 r2«k ak+11 = irfhaf, l = 1, ..., n — 1,

1

2k

4 l al+1 , 1 2 A k+1 / h 2 1 9 k\ fc+1 _ h

— тгп«n ak_1 + ( — rn + ¿rn< )аП+1 = ^rnak, l = n.

_ _r ^ « _

4 n—1 n—1 4 n^ny ^n_ 1 1 ^2t ' n 1 2' n^n n 2t ' n^n'

2k

(14)

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

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

«h(t,r) = ^ «¿(t) ^i(r)

(15)

i=0

с некоторыми коэффициентами «¿(¿), зависящими от времени.

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

r2 <£>l

дм « др \ 1 д 2 ^ « д 2 дР

+ 2 ж) + r2дГ(r р« ) — дТ(r р«) + аТ —

о

1 д / 4 2 д« 4 F

—r2 дД3Rer ^— 3Re—

dr = 0,

l = 0, ...,

(16)

Отметим, что здесь использовано не исходное уравнение (2), а из него вычтено уравнение неразрывности, умноженное на и/2. Это приведёт и в непрерывном и в дискретном виде к

1

о

о

n.

n

1

n

устойчивости в норме пространства ¿2, а не в как в исходном уравнении (2). Применяя к нему формулу интегрирования по частям, получим

1 1

/2 ( ди и др\ 2\ ¡ 2 2 дуг , 1, 2\

Г Уг\Р~д1 + 2~дг) + (УгРи ) Г=1 - г Ри ~дгс1т - 2(УгРи )

о

+

+ У 1 ^ри^т (уг и)Л'Т + (т2угР) 1 - J Р (2т уг + г2 ^ут1) ¿т-

4 ди 4

уг зке - уг зкери

[ дуг ( 4 2 ди 4 , ,

г=1 + Сзйет - зке три^т-

1

-у Т уг^вт ¿Т = 0, 1 = 0, ..., П.

о

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

п 1

( 2 ( ди1 щ др\ 1 2 дуг

'Т У'УЛ ^ + Т^ - 2трии*у^ +

¿=о

(17)

1 2 дуг 4 2 дуг дуг 4 дуг ~

+2т рии^^уг + зкет - зке тригуг^) ¿т+

4

1 4 дуг

2 рииг угуг - уг — +уг з^еу

г = и

1 1 У Р (2туг + т2 £) ¿т - (угР) + ^ т2уг^

¿т 1 = 0, ..., п.

(18)

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

г+1

Е

г=г-1

2 / диг щ др\ 1 2 дуг 1 2 дуг

т угуЛ Р^--1-----т рии;у;—--1--т ЙИИ^— уг +

д£ 2 д^ 2 Г г^г дт 2 ' г дг

4 2 дуг дуг 4 дуг \ , .

+--т р<иг —------тр,игуг —— ¿т+

з Ке дт дт з Ке дт

/1 4 дуг +4

П 2 Р и иг уг уг - уг зКе ^'+ уг зКеМ иг ^

= IР (2туг + т2 ^дд^) ¿т - (угР) + ^ т2уг

о

¿т,

(19)

1 = 0, ....

1

Г

о

1

1

о

о

1

о

о

о

о

1

о

1

Г

1

1

о

п.

Для вычисления интегралов в уравнении (19) применим квадратурную формулу трапеций и равенства «_' = —«о, "_' = "о, вытекающие из симметрии задачи. Аппроксимируем производную по времени следующим образом:

д« « 3Р1 д, Р?+Ч+' — ^0+7 Р? «? Р11к + т "Ж = ^ д« "-Т-■

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

(НГ2Р0 + 1 | 2 Г2..о , 2 „2, А « 0+1+ 1гГо Ро + Н Ие Го "о +3 Н Ие Г1 "Ч «о +

Н]

— зтлйе + зие +1 г^'«« — знк -«"О«0+' =

а/р^1^-о 1 i

У о Нг2 + 2г2Ро0 + 2НгоРо? — 2Р0„2 + Н^о, 1 = 0,

т 2 2

1 „2Р0+1—о__'2_„2'М «0+1_

"„ГРГ 1 — Гг2'°) «°-11+

4 1 3НИе

+('т2^1 + эТни г2_1"-1 + зл^ве г'2'° + г2*'"0*') -о+1+

3 НИе з НИе ^ ' 3 Н]

У1 гг2+1рГ++11-Г+1 — ^ гг2+1"?+1 + г,+ 1"?+1+ + 1 Г^1«0 — -VГ2"0) и?*!1 =

4 1 3 Н Ие гРгУ 1+1

р?+\/р0 «0 1

1 У 1 1 ■ 2 , „2 о0 , ои„ г,0 ^ 2 о? , г__2 т?0

Нг°2 + 2 гг2_1Рг°_1 + 2НггРг0 — - гг°+1Р/+1 + Нгг2Ре0т,г

2

1 = 1, ...,п — 1, (20) 2 Р0 + 1 «0 _ 2 „2 '0__Г _'0

4

п_1рп_1-п_' 3 н ЕеГп_'"п_' 3 ИеГп_'"п_1 1 ,о__2 о А„0+'_

2 0+1 0 2 2 0 + 1 Гпрп «п — 3 н Ие Гп"п^ «п_1+

4' прп «п 3 н Ие «п_1

/ Н 2 0+1 2 2 0 2 2 о 2 О 1 2 0+1 О\ 0+1

^2ТГпРп +3НкеГп_'"п_'— Гп"п + 3ЁеГп"п + 2ГпРп «V-п = = У7рП+^У/рМНг2 + 1 г2 р0 + Нг р0 _ 1 г2р0 + Нг2р0 1 = п

= 2Т п +2 'п_1рп_1 + Н'прп 2 'прп + Н'пгет,п, 1 =

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

п

г) = £ е^Ыг) (21)

¿=о

с некоторыми коэффициентами еД«), зависящими от времени.

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

1

/V 2 д . . г, д , 2 \ д 2 \ 7 д / 2 де)

У 1г ^дй(ре) + ^рдГ(г «) + ^ дГ(г «Ре) — йерТдАГ

1

о

4 2 /дм\2 4 дм\

ЧдтУ + 31е „мдг>) = 0' 1 = 0'-'п- (22)

УдгУ ' 3]

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

1

/( 2 д, 2 л дм 2 д^г 7 2 де

(г ^-(ре) + 2г^мР + г ^Р- - тре— + —г

4 2 /дм\2 4 дм\

ЗЁе г п \дг/ 3

7 де

Г=1

= 0, 1 = 0, ..., п. (23)

Г=1

^Ёе Рг дг

Подставляя в равенство (23) вместо е приближённое решение (21), получим ( П 2 д / л 2 д^г , 7 2 д<Л д^г\ , ,

¿-Д / Г ^т(рег) -г + Ё-РГ„г а7Гг+

1

/ 7 \ Г ( 2 „дм

+ (мр^г^ - ^-^Р-^¿^ ^ = - / (2г^гмР + г -

2 дм

1

о

4 2 /дм\2 4 дм

2 дм 2 4 дм

зие г^г Чд^; +зё- г^г рмдггг, 1=0,...,п. (24)

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

слагаемые:

1

V1 ( [ ( 2 д / \ 2 д^г , 7 2 д^г \ , ,

¿=г-1

о

2 дм

г=1У

о

4 2 дм 2 4 дм

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

(^г2р^+1 | г2 „0 , 7 г2 . 0\ е0 + 1 +

1тГ°Ро +2 Н Ёе РгГо „ + 2 Н Ёе РгГ1 N ео +

+ ( - 7 г2 „0 + 1 г2мй+1рй+1 + 7 г2 .Л е0 + 1 =

+ 1 2 Н Ёе РгГо „ +2Г1 м1 Р1 +2 Н Ёе Рг Г1 „V е1 =

= НгоРое0 - 2 Нгом^Р^ - 1 ^(м^1 - м0^+1) + ^г^(м^+1 - м0^+1)2-2

го „0 м0+1(м0+1 - м0+1), 1 = 0,

3 Ёе

1

Г 1 г2 м0 + 1р0 + 1 7 г2 0 7 г2 „ А е0+1 +

I - 2гг-1мг-1 Рг-1 - 2 Н Ёе Рггг-1 „г-1 - 2 Н Ёе Ргг N ег-1 +

о

1

Н

2 к+1 ,__I 2 к ,__I 2 к,__I

+ 1Тг; Рг +2 Н Не Ргг;-1 Мг-1 + НЁерТг; и +2 НИеРг

,2 ( к \ ек + 1 + г+1 иг+1) ег +

+ ( 7 г2 ,,к + 1 г2 „к+1 к+1 , 7 г2 к Л ек+1 + 1 2 Н ИеРгг; и +2 г1+1и;+1 Р;+1 +2 Н Ие Рг г;+1 иг+Ч 6;+1

Н „2„к_к ои „,ик + 1п 1 г2рк(и к+1 ик+1- 1

к+1

'"2 к к о и к + 1 о 2 г>к( к+1 к+1\ ,

= - гг Р г е г — 2 Нгг и г + рг — о гг рг (и г +1 — и г -1) +

2 к/ к + 1 к + Ь2

зННе' (иг+1 — и )

г-1)

—г, и к и к+1(и к+1 - и к+1 3Не 'гиг "г (иг+1 "г-1

+1

2 +1 +1 гп-1ип-1рп-1

7

2 Н Ие Рг

г-1

г2 и к гп-1 и п — 1

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

(26)

7

2 Н Ие Рг

гп ип I еп-1

2 ,, к \ „ к +1

1+

+ ( 2т г"Рп+1 +2 Н Ие Ргг"-1 иП-1 + 2 г>Г±РГ± 2 Н Ие е» ~ _

2

2 +1 +1

7

2 +1

_г2 Ркек _ Нг ик+1Рк _ 1 г2Рк(ик+1 _ ик+1) + 2 г2ик (и^ _ ик+1 )2_

2т

3 Н Ие

2

+1 +1 +1 гп ип ип+ (ип+ — ип—1) , 1 = п.

+1

3Ие

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

соУо — ЬоУ1 = /о, г = 0, — ОгУг-1 + СгУг — 6^+1 = /¿, г = 1, 2, ..., N — 1, —а« У« -1 + с« У« = /«, г = N,

определяется следующим образом. Сначала задаются начальные значения для С, А, Г и Ф:

С = со, А = аь Г = /о, Ф = /1

и параметра Ж1 = 0.

Далее последовательно для г = 0,1,..., N — 1 выполняется прямой ход. В зависимости от ситуации проводятся следующие действия: если |С| ^ 1|, то

^+1 = Ь^С, &+1 = Г/С, С = ^+1 — Ао^+ь Г = Ф + А&+1, ©i+l = ^+1 = г +1,

А = №¿+2, Ф = /i+2,

если |С| < 1|, то

Ог+1 = С/^, вг+1 = — Г/б^ С = ^+1^+1 — А, Г = Ф — Сг+2вг+1, ©¿+х = г +1, ^+1 =

А = ^+2 «¿+Ъ Ф = /а+2 + ai+2 вi +1.

1

2

1

Замечание. Для г = N — 1 переопределение А и Ф в обоих случаях не приводится.

После окончания прямого хода осуществляется обратный ход по определению неизвестных в обратном порядке. Сначала вычисляется неизвестное уп, где п = Жд по формуле

Уп = Р/С,

а затем последовательно для г = N — 1, N — 2,..., 0 вычисляются остальные неизвестные

Ут = «¿+1Уп + вг+1, где т = ©¿+1, п = жт, г = N — 1, N — 22,..., 0.

Таблица 1

Н, км 15 60 100 200 300 400 600 1100 1600 2700 2870 2900 4700 6371

Р, кбар 3,5 14 29 64,9 99 130 190 30 530 1240 1340 1350 3040 3632

Т, К 673 873 1103 1473 1623 1673 1723 1873 2700 3500 4000

Р, г/см3 2,85 3,34 3,37 3,36 3,48 3,54 4,13 4,74 5,03 5,55 5,68 9,89 12,3 13,0

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

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

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

Рис. 2. Распределение плотности по радиусу Земли для различных моментов времени под действием гравитационных сил

Рис. 3. Распределение температуры по радиусу Земли для различных моментов времени под действием гравитационных сил

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

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

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

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

[1] Г.Джеффрис, Земля, её происхождение, история и развитие, М., Изд-во иностр. лит., 1960.

[2] Н.Л.Добрецов, А.Г.Кирдяшкин, А.А.Кирдяшкин, Глубинная геодинамика, Новосибирск, Изд-во СО РАН, филиал "ГЕО", 2001.

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

[4] В.В.Шайдуров, Г.И. Щепановская, Математическое и численное моделирование нестационарного распространения импульса энергии большой мощности в вязком теплопроводном газе, Часть 1. Математическая постановка задачи, Красноярск, ИВМ СО РАН, 2003, 50 с. (Деп. в ВИНИТИ 24.10.03, № 1860-В2003).

[5] E.D.Karepova, A.V.Malyshev, V.V. Shaidurov, G.I. Shchepanovskaya, The parallel realization of the finite element method for the Navier-Stokes equations for a viscous heat conducting gas, Notes on Numerical Fluid Mechanics and Multidisciplinary Design, Springer, 91(2006), 41-54.

[6] В.А.Малышев, В.В.Шайдуров, Г.И.Щепановская, Математическое и численное моделирование на многопроцессорной вычислительной системе сверхзвукового взаимодействия тепловых импульсов, Вестник КрасГУ, 2006, № 4, 109-116.

[7] В.В.Шайдуров, Г.И.Щепановская, Газодинамическая модель внутреннего строения Земли, Вестник СибГАУ, 2008, № 1, 79-83.

The finite element method for modeling spherically symmetric pulsations of the earth

Alexander V.Vyatkin Vladimir V.Shaidurov Galina I.Shchepanovskaya

Some problems of mathematical and numerical modeling of geodynamic processes of expansion, compression, heating and cooling of the Earth are considered. Keywords: mathematical and numerical modeling, geodynamic processes