CC BY
56
Вычислительные технологии
Том 16, № 5, 2011
Нестационарная модель конвективных мантийных течений в приближении слабосжимаемой жидкости*
Г. Г. Лазарева1 , О. П. Полянский2, М. П. Федорук3, A.B. Бабичев2 , В. А. Вшивков1, В. В. Ревердатто2 1 Институт, вычислительной математики и математической геофизики СО РАН, Новосибирск, Россия
2
3Институт вычислительных технологии, СО РАН, Новосибирск, Россия e-mail: [email protected], [email protected], [email protected], [email protected], [email protected], [email protected]
Представлена новая нестационарная модель течений в мантии Земли в приближении слабосжимаемой жидкости. Цель работы — моделирование конвективных течений с сильно изменяющимися реологическими и транспортными свойствами, такими как вязкость, плотность, теплопроводность и т. п. В отличие от традиционного подхода, основанного на приближении Буссинеска, модель основана на решении системы полных классических уравнений Навье — Стокса, описывающих динамику слабосжимаемой жидкости с переменными плотностью и вязкостью. Для улучшения сходимости построенной численной модели при малых числах Маха использован метод предобуславливания. Приведены обоснование и верификация модели. Модель разработана в качестве основы для параллельного алгоритма.
Ключевые слова: полные уравнения Навье — Стокса, метод предобуславливания, гравитационная термоконвекция.
Введение
Исследование эволюции гравитационно-неустойчивых систем в мантии Земли является одной из актуальных задач геодинамики. Решение этой задачи связано с чрезвычайно важной проблемой эндогенной геологии — анализом процессов тепломассопереноса в земной коре. Проблема устойчивости или перераспределения вещества и энергии в недрах Земли в настоящее время изучается в рамках равновесной термодинамики или неизотермической механики сплошной среды. Один из возможных плодотворных подходов состоит в том, что если рассматривать вещество земной коры как сплошную среду, то при изучении процессов деформации в геологическом времени можно с некоторыми ограничениями использовать законы сохранения для упругоплаетичного/вязкого материала. Среди эндогенных механизмов переноса вещества в земной коре и в верхней мантии главное место занимают движения, причинами которых являются гравитационная неустойчивость и тепловые потоки.
Формирование крупных изверженных провинций на континентальной и океанической плитах связывают с нижнемантийными суперплюмами, поднимающимися с глубин
* Работа выполнена при финансовой поддержке Интеграционного проекта СО РАН № 2в рамках реализации ФЦП "Научные и научно-педагогические кадры инновационной России" на 2009^2013 годы.
границы мантии и ядра |1|, Физические аспекты процесса венльшапия илюма термической или термохимической природы достаточно хорошо изучены с помощью физического и математического моделирования |2, 3|. Эта концепция включает ьспльшапие в локальной области легкого высокотемпературного и маловязкого мантийного материала — илюма на фоне крупномасштабных конвективных течений. Неясными остаются вопросы о конечном этане эволюции дианиров: как высоко они могут подниматься и каково соотношение подъемной силы и вязкого сопротивления вещества при подъеме на верхние уровни литосферы. Таким образом, дианировый механизм транспорта магмы в наиболее вязкой и холодной части мантийной литосферы требует тщательного изучения. Важными вопросами являются следующие: 1 — возможно .ни проникновение частично расплавленного вещества сквозь жесткую литосферу?, если да, то 2 — каковы механизм подъема и форма всплывающих тел? и 3 — какова роль реологии литосферы, испытывающей фазовый переход при плавлении, в процессе дианиризма?
В статье представлена новая нестационарная модель геодинамических течений в приближении слабосжимаемой жидкости. Модель строилась дня описания процесса плавления и вызванного им венльшапия легкого вещества в результате андерн.нейтинга базитовой высокотемпературной магмы иод основанием коры. Цель моделирования — нахождение параметров процесса плавления и дианиризма в нижней коре, определение структуры течения всплывающей гранитной магмы и предсказание возможной формы гранитогнейсовых дианировых тол.
Традиционно прямое численное моделирование процессов в недрах Земли основано па решении уравнений Навье-Стокса в приближении Обербека-Буссипоска |4-7|. К решению уравнений, записанных в переменных вихрь — функция тока, применяется как конечно-разностный метод, так и метод конечных элементов 18 131. В геодинамике приближение Обербека — Буссипоска является достаточно общепринятым и его применение обоснованным |14, 151 дня случая малых изменений плотности. Модели, основанные на
100 200
Н, км
300 400 500
4.1
3.5-
и
и. з.з -
о.
2.9-
3.7"
3.1-
/
2.7-
2.5
0 2 4 6 8 10 12 14 16 18 20 Р, ГПа
Рис. 1. Плотность водосодержащих и "сухих" иеридотитовых расплавов в Р— Т-условиях литосферы и верхней мантии в сравнении с профилем плотности средней мантии
приближении Обербека — Буссинеска, используются при решении задач о конвекции в мантии, где рассматривается только изменение плотности вследствие теплового расширения (около 30 кг/м3), В рассматриваемой задаче учитываются скачки плотности из-за фазовых переходов при плавлении, что составляет до 350 кг/м3, т.е. на порядок больше. На рис, 1 представлен график плотности водоеодержащих и "сухих" перидо-титовых расплавов в Р — Т-условиях литосферы и верхней мантии [16] в сравнении с профилем плотности средней мантии АК-135 по данным [17], принятым в модели. Серым показана увеличивающаяся разность плотностей при подъеме к поверхности. В работе предпринята попытка создания модели в приближении елабоежимаемой жидкости. Модель основана на решении полной системы уравнений Навье —Стокса с учетом переменной вязкости, дополненной уравнением состояния. Для улучшения сходимости построенной численной модели при малых числах Маха использован метод предобу-елавливания.
1. Математическая модель
Рассмотрим систему уравнений, описывающую динамику елабоежимаемой жидкости, замкнутую уравнением состояния
| + У(„и)=0,
<9и 3 / дщ ^ дик 2 ^ <9ит\
дЬ р ^ рдхк^\дхк дхг 3 гкдхт) ^ у'
^ + (иУ)Т = ЪкЪ ■ Т,
которое является прямым следствием выражения для плотности р = ру=0(1 — аТ+ в(Р — Ру=о))-
Р = РУ=о + о ( + аТ-1
в \Ру=о
где р — плотность, и = (пх,пу) — вектор скорости, Т — температура, р — давление, П — вязкость, д — ускорение свободного падения, к — температуропроводность, ру=0 и ру=0 — плотность и давление на поверхности.
Характерные значения переменных задачи следующие: характерная длина Ь0 = 3 ■ 104 м, ро = 2.8 ■ 103 кг/м3, То = 550 °С, ро = Поко/Ц Па, к = 10-6 м2/с, ¿о = ¿о/к по = £0Д0. В уравнении состояния приняты параметры а = 3-10-5 °С-1, в = 10-11 Па-1. Для
определения вязкости в модели использовано уравнение Аррениуса п = А- ехр (- ),
\ ЯпТ )
где А =1.2 ■ 1017, Е = 2.16 ■ 104, п = 2.6 — экспериментальные данные, Я — универсальная газовая постоянная. В представленной постановке задачи значение коэффициента вязкости находится в диапазоне п = 1018 — 1020 Па/с, следовательно, число Прандт-
ля Рг= -Щ- = 3.6 • Ю20 4- 3.6 • 1022, число Рэлея Еа= адр0^1° = 2.7 • 102 4- 2.7 • 104,
роко Поко
где 9 — характерный перепад температуры. В модели рассматривается нормальная кора
с экспоненциальным распределением радиоактивных источников тепла (см. 1 на рис. 2).
Рис. 2. Типы стационарных геотерм, характерные для разных термических режимов в земной коре: 1 нормальная кора с экспоненциальным распределением радиоактивных источников тепла, 2 утолщенная кора с гранитоидными интрузиями, 3 нормальная кора с повышенным тепловым потоком
Выбор модели слабосжимаемой жидкости определяется как желанием использовать более полную модель процесса с учетом скачков плотности, вызванных фазовыми переходами при плавлении, так и возможностью создания численной технологии решения с привлечением хорошо апробированных конечно-разностных схем. Известен |18| критерий применимости классической модели Обербека —Буссинеска дня описания тепловой
т? t gL0po ^
гравитационном конвекции. Ьс.ни параметр £ = —-— имеет порядок < 1. то модель
Voko
Обербека —Буссинеска не применима. Величина £ характеризует относительный вклад факторов плавучести и объемного расширения жидкости в формирование поля скоростей. В рассматриваемом классе задач значение коэффициента вязкости может варьироваться в диапазоне n = 1014 ^ 1028 Па/с, тараметр £ принимает значения от 105 до 10-9.
Известно, что численное интегрирование полной системы уравнений Навье — Стокса является весьма сложной и трудоемкой вычислительной задачей, требующей дня своего решения разработки специальных конечно-разностных схем и численных алгоритмов, причем проблема разработки эффективных методов численного моделирования дозвуковых течений вязких газов |19, 20| на основе полной системы уравнений Навье — Стокса связана именно с решением внутренних задач. В связи с этим широко используются приближенные математические модели, не обладающие общностью, характерной дня полной системы уравнений Навье —Стокса 118 211. Однако в тех случаях, когда допущения, .нежащие в их основе, справедливы, применение таких моделей в методическом отношении в сравнении с системой уравнений Навье — Стокса является более оправданным.
Исторически первой и наиболее простой является классическая модель вязкой несжимаемой жидкости, основанная па системе уравнений Стокса |6|, В дальнейшем с цолыо расширения границ применимости этой модели была сформулирована приближенная система уравнений дня описания конвекции жидкостей и газов, получившая название приближения Буссинеска. Последняя модель широко используется при численном моделировании геодипамических процессов, определяющихся конвективными течепи-
ями среды [22-25], Всесторонний анализ методов численного интегрирования систем уравнений Стокса и Буссинеска содержится в монографиях [26, 27], Для расчета стационарных течений более эффективно использование методов, основанных на введении искусственной сжимаемости [28], Примером успешного применения такого подхода для моделирования существенно трехмерных процессов в мантии Земли является работа [29], В рамках приближения Буссинеска изменения плотности в потоке учитываются лишь частично. Плотность считается зависящей от температуры только при определении массовой силы, К числу допущений, отказ от которых не затрагивает основ модели вязкой несжимаемой жидкости, следует отнести также постоянство теплофизичееких и переносных свойств газа или пренебрежительно малую роль процессов диссипации механической энергии потока и работы сил давления. Очевидно, что модификации модели несжимаемой жидкости не могут применяться для течений с существенно изменяющейся плотностью из-за использования в них уравнений несжимаемости, В связи с этим начиная с 1980-х гг. для описания рассматриваемого класса течений развиваются новые приближенные подходы, свободные от недостатков модели несжимаемой жидкости и приближения Буссинеска,
Наибольшее распространение при численном моделировании внутренних течений получили параболизованные модели, в том числе приближение пограничного слоя [30], приближение узкого канала [31] и модель вязкого ударного слоя, основанная на параболической аппроксимации полных уравнений Навье-Стокса [32], Последняя модель особенно эффективна в случае параболизации в криволинейной системе координат, продольное направление которой совпадает с локальным направлением линий тока [33], Для описания нестационарных течений многокомпонентных газовых смесей был разработан подход [34], позволяющий, оставаясь в рамках приближенных моделей, учитывать динамику удельного объема среды за исключением изменений, вызванных воздействием неоднородности поля давления.
Нестационарные математические модели конвекции в большинстве случаев основаны на различных приближенных моделях, хотя нестационарная двумерная численная модель сжимаемого вязкого теплопроводного газа реализована в 1967 г, [27, 35], Наряду с этим были широко распространены модели гипозвуковых течений сжимаемого газа в задачах обтекания и разработан ряд моделей для описания химических процессов в сжимаемых средах с учетом температуры. Отличительными чертами задач мантийной конвекции являются значение числа Прадтля Рг « 1020, нелинейное уравнение состояния, разномасштабность различных процессов, сильно (на десятки порядков) изменяющиеся значения вязкости. Эти особенности создают дополнительные трудности при численной реализации модели с учетом сжимаемости среды. Тем не менее в настоящее время возможности вычислительной техники позволяют проводить численное моделирование самых подробных моделей конвективных течений с учетом сильно изменяющихся реологических и транспортных свойств, таких как вязкость, плотность и теплопроводность,
2. Алгоритм решения
При небольших скоростях геодинамических процессов (несколько сантиметров за миллион лет) для рассматриваемого типа задач характерны высокая скорость звука и малое
число Маха: с
вро ' с
^ = 104 м/с, М - 10"13.
Геодинамика рассматривает очень медленные течения, поэтому здесь в отличие от сейсмологии и других разделов геофизики в теории ранее не использовалось число Маха, Формально вычислив число Маха, можно обратиться к опыту расчетов в области существенно дозвуковых течений,
К настоящему времени разработан широкий класс моделей, основанных на численном решении полной системы уравнений Навье — Стокса, Известно [36], что при уменьшении характерного числа Маха потока до значений ниже 0.1 происходит замедление сходимости итераций по времени в методе установления и ухудшение точности получаемого решения. Замедление сходимости в этом случае объясняется возрастающей при М! —у 0 жесткостью уравнений динамики сжимаемого газа, определяемой как отношение максимального и минимального собственных значений матрицы Якоби векторов конвективных потоков. Детальное исследование этого эффекта с учетом используемых разностных схем, сетки и параметров течения проведено в [37],
Для избежания проблем со сходимостью при М < 1 был предложен метод предо-
—0
дификации члена с производной по времени в исходных уравнениях движения путем домножения его на матрицу, которая подбирается так, чтобы скорости распространения —0
новлении система описывает стационарное решение исходных уравнений Эйлера, Метод предобуелавливания применительно к уравнениям газовой динамики при малых числах Маха впервые использован в работе [36], В 1987 г, была опубликована теоретическая работа [38], посвященная методам предобуелавливания для уравнений несжимаемой и сжимаемой жидкости, В ней автор распространил метод искусственной сжимаемости [39] на случай сжимаемого газа, получив двухпараметричеекое семейство матриц предобуелавливания, отличных от найденной в [36], но также обеспечивающих одинаковый порядок собственных значений модифицированной матрицы Якоби при всех М < 1, В это же время в работе [40] был успешно применен "метод масштабирования сжимаемости" (по сути — метод предобуелавливания) для расчета вязких течений сжимаемого <1
етве одной из зависимых переменных относительного давления p = p — pxap, где pxap —
<1
также гибридного способа аппроксимации производных конвективных членов и членов с давлением. Близкий подход предложен для численного моделирования конвективных течений и теплопереноеа в жидкостях с параметрами вблизи термодинамической критической точки [41], где выполняется двухмаештабное расщепление давления: полное давление заменяется суммой среднего по области и приращения, преобразуется член с давлением в уравнении импульсов. Это сохраняет систему полной, но весьма полезно при создании эффективного численного кода. Таким образом, широкое применение находят подходы, близкие к методу предобуелавливания, но адаптированные для рассматриваемого класса задач.
Для масштабирования элементов решения системы уравнений введена операция приведения уравнений к безразмерному виду с привлечением безразмерных параметров
течения К = , ¡3 = ; ¿х = aT0, р = ру=о/ро, р = ру=о/ро [42], Все параметры ар06 L0
течения представлены в виде f = f0fгде f — физическая характеристика, f0, f' — ее характерное и безразмерное значения. Для решения нестационарной задачи введем итерации по фиктивному времени [43] для уравнений движения. Следуя идеологии рас-
смотренных подходов, умножим производную но фиктивному времени на коэффициент крг, имеющий значение порядка величины числа Прандтля, На каждом шаге по реальному времени выполняются итерации дня уравнений движения по фиктивному времени. Модифицированные уравнения не будут жесткими, что позволит обеспечить сходимость к искомому решению па каждом шаге по реальному времени. Итоговая система уравнений в безразмерных переменных имеет вид
др др др (ди,х ди% д1 11 х дх 11 у ду ^ \ дх ду
дих дих дих дих 1 / др дохх даху
КРгТГ" + -77Г + Их^— + = РГ~ -Т- + +
дт дЬ дх ду р \ дх дх ду
диу дщ диу диу 1 ( др даху , д&уу\ -
дт дЬ дх ду р \ ду дх ду
^ и ^ и ^ - А —
д1 х дх у ду дх \ дх) ду \ ду /
р = р + — (- + аТ — 1 в \Р
где тензор вязких напряжении
@хх П
4 дих 2 диу,\ /4 диу 2 дих\ (дих диу"
3 дх 3 ду ) \3 ду 3 дх ) \ ду дх
Численное интегрирование полной системы уравнений Навье — Стокса представляет собой чрезвычайно сложную и трудоемкую задачу. Поэтому эта система стана полигоном, па котором проходят испытания многие современные численные методы решения |44|, Особенность численной реализации полной системы уравнений Навье —Стокса дня течений с малыми числами Маха состоит в значительном различии двух характер-пых дня существенно дозвуковых течений масштабов времени: характерного времени
Рис. 3. Зависимость невязки е от числа итераций N по фиктивному времени Дт для трех шагов по времени Дt
конвективных процессов £ = Ь/ю н характерного времени распространения акустических возмущений т = Ь/(ю + с). При использовании явных конечно-разностных схем в соответствии с условием устойчивости Куранта шаг по времени не может превышать характерного времени наиболее быстрого процесса передачи акустических возмущений. Поэтому применение явных конечно-разностных схем за исключением ряда специфических задач оказывается совершенно неоправданным. В настоящей работе система уравнений реализована неявным конечно-разностным методом стабилизирующей поправки первого порядка по времени и пространству на регулярной прямоугольной сетке в декартовой системе координат. Сходимость невязки для вектора скорости иллюстрирует рис. 3.
3. Решение модельной задачи
Рассматривается прямоугольная область земной коры глубиной утах = 30 км и шириной хтах = 60 км (рис. 4). На верхней границе области задана свободная поверхность с постоянным нулевым значением температуры, приняты плотность 2.8 • 103 кг/м3 и давление 105 Па. Боковые границы области изолированы для передачи тепла и выхо-
20
до температуры 1200 °С, т.е. рассматривается процесс "андерплейтинга". На остальной части нижней границы принята температура 550 °С. Коэффициент температуропроводности линейно зависит от начального распределения температуры и не зависит от времени. Для согласования граничных условий температура вдоль нижней границы изменяется по экспоненциальному закону. В ходе расчета скорости в зоне прогрева задаются согласованно в соответствии с исходными уравнениями, что позволяет получать гладкие значения скоростей. Давление на нижней границе в начальный момент времени принято равным 108 Па. Расчеты проводились при Н = 0.5, Д£ = 10-3, Дт = 10-4.
В начальный момент времени среда неподвижна, распределения плотности и давления задаются согласно уравнению состояния с учетом замены переменных для давления. Начальное распределение температуры принимается с учетом потока:
Т\г=ъ = Т0 + (Т550 - Т0)у/ут^ + +(Т12оо - Т550) ехр ^ ^ ехр У) ^ ?
где Яу — ширина распространения температуры вдоль х и у, х0 — -^тах, Т0 = 0 °С,
2
Т550 = 550 °С, Т1200 = 1200 °С. В течение всего расчета происходит разогрев нижней
10-
7
Рис. 4. Начальное распределение температуры в расчетной области (зона прогрева выделена черным цветом)
Рис. 5. Результаты расчета с постоянной (а) и переменной (б) вязкостью; изотермы в моменты времени £ = 0, £ =1.5 • 1013 с и £ = 2.3 • 1013 с (сверху вниз)
границы. В результате конвективных процессов и переноса разогретого вещества предполагаются плавление и всплывание легкой гранитной магмы.
На рис. 5 показаны результаты расчетов в форме меняющегося температурного поля в разные моменты геологического времени (1 млн лет = 3.15 • 1013 с). Из рис. 5, а, где представлен процесс распространения начального температурного возмущения при постоянной не зависящей от температуры вязкости (1019 П а •с), следует, что процесс теплопереноса характеризуется теплопроводностью без заметного влияния конвективных течений. На рис. 5, б изображена эволюция температурного поля при переменной зависящей от температуры вязкости, изменяющейся в диапазоне 1018 ^ 1020 Па • с. Видно, что снижение вязкости в теле частично расплавленного материала (диапира) влияет на режим теплопереноса в коре: вклад конвективных течений становится существенным. Подобные режимы конвекции могут наблюдаться в зонах срединно-океанических хребтов, где на незначительных глубинах в коре (5-7 км) сейсмическими методами фиксируются области плавления (магматические камеры).
Заключение
В статье предложена новая нестационарная модель мантийных течений в приближении слабосжимаемой жидкости. Необходимость создания такой модели основана на специфике задачи, предполагающей значительные изменения плотности, подтверждена критерием Пухначева. Показан вектор развития математических моделей мантийных течений от модели Стокса к модели Обербека —Буссинеска и далее к более сложным
моделям, адаптированным к особенностям решаемых задач. Представлен краткий обзор работ по численному моделированию дозвуковых течений вязких газов. Рассмотрены различные подходы к численному решению внутренних задач, описывающих существенно дозвуковые течения. Приведен алгоритм решения задачи с использованием интегрирования уравнений движения по фиктивному времени. Для улучшения сходимости построенной численной модели при малых числах Маха использован метод про. i-обуелавливания, адаптированный к специфике задачи с большим числом Прандтля,
Выполнено численное моделирование процесса плавления и вызванного им веплы-вания легкого вещества в результате андерплейтинга базитовой высокотемпературной магмы под основанием коры. Анализ эволюции течения позволил определить характеристики процесса гравитационной неустойчивости в коре над локальным тепловым источником, Результаты моделирования показали, что скорость подъема и форма диапира определяются зависимостью вязкости от температуры. При заданном законе вязкости скорость веплывания составляет 0,02 м/год, что хорошо согласуется с геологическими наблюдениями.
Авторы выражают благодарность д-ру физ.-мат, наук проф. Г.Г. Черных и канд. физ.-мат, наук Д.В, Чиркову за участие в обсуждении представленной численной модели.
Список литературы
[1] Доврецов Н.Л. Глобальная геодинамическая эволюция Земли и глобальные геодинамические модели // Геология и геофизика. 2010. Т. 51(6). С. 761-787.
[2] Кирдяшкин A.A., Доврецов Н.Л., Кирдяшкин А.Г. и др. Гидродинамические процессы при подъеме мантийного плюма и условия формирования канала излияния // Там же. 2005. Т. 46(9). С. 891-907.
[3] Farnetani C.G., Richards М.А. Numerical investigations of the mantle plume initiation model for flood basalt events //J. Geophvs. Res. 1994. Vol. 99. F. 13813-13833.
[4] Андреев B.K., Капцов O.B., Пухначев B.B., Родионов A.A. Применение теоретико-групповых методов в гидродинамике. Новосибирск: Наука, 1994. 319 с.
[5] Ландау Л.Д., Лифшиц Е.М. Гидродинамика. М.: Наука, 1986. 736 с.
[6] Лойцянский Л.Г. Механика жидкости и газа. М.: Наука, 1978. 736 с.
[7] Монин A.C., Яглом A.M. Статистическая гидромеханика. Т. 1. СПб.: Гидрометеоиз-дат, 1992. 264 с.
[8] Полежаев В.П., Бунэ A.B., Верезув H.A. и др. Математическое моделирование тепломассообмена на основе уравнений Навье — Стокса. М.: Наука, 1987. 271 с.
[9] Берковский Б.М., Ноготов Е.Ф. Газностные методы исследования задач теплообмена. Минск: Наука и техника, 1976. 144 с.
[10] Берковский Б.М., Полевиков В.К. Вычислительный эксперимент в конвекции. Минск: Изд-во Университетское, 1988.
[11] Тарунин Е.Л. Вычислительный эксперимент в задачах свободной конвекции. Иркутск: Изд-во Иркутского гос. ун-та, 1990. 228 с.
[12] Воеводин А.Ф., Остапенко В.В., Пивоваров Ю.В., Шугрин С.М. Проблемы вычислительной математики. Новосибирск: Изд-во СО ГАИ, 1995. 154 с.
[13] Марчук Г.И. Методы расщепления. М.: Наука, 1988. 264 с.
[14] McKenzie D.P., Weiss N.O. Speculations on the thermal and tectonics history of the Earth // Geophvs. J. Roy. Astr. Soc. 1975. Vol. 48. P. 131 17!.
[15] Schubert G., Turcotte D.L. Phase changes and mantle convection // J. Geophvs. Res. 1971. Vol. 76. P. 1424.
[16] Отани Э., Дапэн Чжао. Роль воды в глубинных процессах в верхней мантии и переходном слое: Дегидратация стагнирующих субдукциоппых плит и ее значение для "большого мантийного клина" // Геология и геофизика. 2009. Т. 50(12). С. 1385-1392.
[17] Kennet B.L.N., Engdahl E.R., Buland R. Constraints on seismic velocities in the Earth from travel times // Geophvs. J. Intern. 1995. Vol. 122. P. 108-124.
[18] Андреев В.К., Гапоненко Ю.А., Гончарова О.Н., Пухначев В.В. Современные математические модели конвекции. М.: Физматлит, 2008. 368 с.
[19] Лапин Ю.В., Нехамкина О.А., Поспелов В.А и др. Численное моделирование внутренних течений вязких химически реагирующих газовых смесей // Итоги науки и техники. МЖГ. Т. 19. М.: ВИНИТИ, 1985. С. 86485.
[20] Нехамкина О.А., Никулин Д.А., Стрелец М.Х. Об иерархии моделей тепловой естественной конвекции совершенного газа // Теплофизика высоких температур. 1989. Т. 27, № 6. С. 1115-1125.
[21] Березин Ю.А., Жуков В.П. О влиянии вращения на конвективную устойчивость крупномасштабных возмущений в турбулентной жидкости // Изв. АН СССР. МЖГ. 1989. № 4. С. 3-9.
[22] van К к к ex P. Evolution of starting mantle plumes: A comparison between numerical and laboratory models // Earth and Planetary Sci. Lett. 1997. Vol. 148. P. 1-11.
[23] Fullsack P. An arbitrary Lagrangian-Eulerian formulation for creeping flows and its application in tectonic models // Geophvs. J. Intern. 1995. Vol. 120. P. 1-23.
[24] Gerya Т., Yuen D.A. Characteristics-based marker-in-cell method with conservative finite-differences schemes for modeling geological flows with strongly variable transport properties // Phvs. Earth Planet. Intern. 2003. Vol. 140. P. 293-318.
[25] blttner D., Schmeling H. Numerical modelling of melting processes and induced diapirism in the lower crust // Geophvs. J. Intern. 1995. Vol. 123. P. 59-70.
[26] Роуч П. Вычислительная гидродинамика. М.: Мир, 1980. 612 с.
[27] Пасконов В.М., Полежаев В.П., Чудов Л.А. Численное моделирование процессов тепло- и массообмена. М.: Наука, 1984. 285 с.
[28] Владимирова Н.Н., Кузнецов Б.Г., Яненко Н.Н. Численный расчет симметричного обтекания пластинки плоским потоком вязкой несжимаемой жидкости // Некоторые вопросы прикладной и вычислительной математики. Новосибирск, 1966. С. 186-192.
[29] Червов В.В. Моделирование трехмерной конвекции в мантии Земли с применением неявного метода слабой сжимаемости // Вычисл. технологии. 2009. Т. 14, № 3. С. 86-92.
[30] Лапин Ю.В. Турбулентный пограничный слой в сверхзвуковых потоках газа. М.: Наука, 1982. 312 с.
[31] Williams J.С. Viscous compressible and incompressible flow in slender channels // AIAA J. 1963. Vol. 1, No 1. P. 186-195.
[32] Davis R.T. Numerical solution of the hypersonic viscous shock-layer equation // Ibid. 1970. Vol. 8, No. 5. P. 843-851.
[33] Ковеня В.М., Черный С.Г. Метод решения стационарных упрощенных уравнений вязкого газа. Препр. ИТПМ СО АН СССР. 1981. № 42. 51 с.
[34] Стрелец М.Х. О численном моделировании существенно дозвуковых течений газов и газовых смесей при наличии значительных изменений плотности // Динамика неоднородных и сжимаемых сред. Л.: Изд-во ЛГУ, 1984. С. 70-83.
[35] Полежаев В.И. Численное решение системы двумерных нестационарных уравнений На-вье^Стокса для сжимаемого газа в замкнутой области // Изв. АН СССР. МЖГ. 1967. № 2. С. 103-111.
[36] Чой Д., Меркл Ч.Л. Применение метода установления для расчета низкоскоростных течений // Аэрокосмическая техника. 1986. № 7. С. 29-37.
[37] Чирков Д.В. Моделирование гипозвуковых течений с использованием предобусловлен-ных уравнений Эйлера и Навье — Стокса // Материалы конф. "Вычислительные и информационные технологии в науке, технике и образовании". Ч. 3. Новосибирск; Алматы; Усть-Каменогорск, 2003. С. 262-271.
[38] Turkel Е. Preconditioned methods for solving the incompressible and low speed compressible equations //J. Сотр. Phvs. 1987. Vol. 72. P. 277-298.
[39] Chorin A.J. A numerical method for solving incompressible viscous flow problems // Ibid. 1967. Vol. 2. P. 12-26.
[40] Стрелец M.X., Шур M.Л. Метод масштабирования сжимаемости для расчета стационарных течений вязкого газа при произвольных числах Маха // Журнал вычисл. математики и матем. физики. 1988. Т. 28. С. 254-266.
[41] Полежаев В.П., Соболева Е.Б. Нестационарные эффекты тепловой гравитационной конвекции околокритической жидкости при боковом нагреве и охлаждении // Изв. РАН. МЖГ. 2002. № 1. С. 81-93.
[42] Turcotte D.L., Torrance К.Е., Hsui А.Т. Convection in the earth's mantle in methods // Comput. Phvs. 1973. Vol. 13. P. 431-454.
[43] Weiss J.M., Smith W.A. Precondition applied to variable and constant density flows // AIAA J. 1995. Vol. 33. No. 11. P. 2050-2057.
[44] Ковеня В.M., Яненко H.H. Метод расщепления в задачах газовой динамики. Новосибирск: Наука, 1981. 304 с.
Поступила в редакцию 15 апреля 2011 г., с доработки — 15 июня 2011 г.
