Рис. 2. Результаты моделирования всплывания диапира с постоянным пределом текучести материала коры. Приведены картины поля температуры в теле всплывающего диапира (в цвете) и вне его (в изолиниях), от начала всплывания (а) до конечного момента времени 2.162 млн. лет (д). Шкала в диапазоне 650-1200°С показана слева, изотермы в °С.

Рис. 3. Результаты моделирования всплывания диапира при температурно-зависимом пределе текучести. Показано поле температуры в теле всплывающего диапира (в цвете) и вне его (в изолиниях), а также форма диапирового тела в модели при переменном пределе текучести материала коры для этапа от начала всплывания (7.36 млн. лет) (а) и до финальной формы остывшего диапира (7.95 млн. лет) (д).


Вторая серия тестов проводилась с целью исследования влияния температурно-зависимого предела пластичности на динамику всплывания диапира. Для этого в выражении температурной зависимости предела текучести σy = 106 (Па)+A exp(-E/RT) варьировались параметры А и Е и определялась высота подъема и форма диапира. В этой серии экспериментов расчет велся только в верхней 30-км области моделирования, а жесткое основание (слой 8 км) заменялось соответствующими граничными условиями. Граничные условия на нижней границе области моделирования переносятся на границу y=0 без жесткого основания. Расчеты показывают, что это упрощение не сильно сказывается на динамике всплывания и форме диапира. Варьируя параметры А и Е, предел текучести изменяется от 1 МПа в частично расплавленной породе до максимального значения (10 МПа), соответствующего приповерхностной температуре. В этих расчетах оказалось, что высота всплывания ограничена некоторым уровнем глубины, выше которого деформации становятся малыми ввиду низкой температуры.

НЕ нашли? Не то? Что вы ищете?

На рис. 3 показаны результаты расчетов по модели температурно-зависимой пластичности, когда предел текучести изменяется в диапазоне от 1 до 10 МПа для интервала от температуры плавления (650ºС) до поверхностной. В сравнении с предыдущим вариантом подъем диапира ограничен уровнем средней коры 14-16 км. Его двухмерная форма при всплывании меняется от арочной до грибовидной. Частичного расплавленный материал растекается в горизонтальном направлении, формируя тело шириной до 31 км; канал диапира (“ножка”) существенно короче, чем в модели с постоянной пластичностью, его высота – 3.5 км, ширина – 2.3 км (рис. 3 г-д). Форма диапира похожа на вертикальное сечение лополита, что может вызвать сложности в идентификации природы магматического тела.

Численные эксперименты позволяют сделать следующие основные выводы. 1) Для того чтобы в гравитационном поле началось всплывание, должен сформироваться критический объем частично-расплавленного вещества. По модельным оценкам высота области плавления в гранитной коре должна быть не менее 6-7 км. 2) Независимо от размера теплового источника (фиксированной или переменной ширины) во всех моделях наблюдалась грибовидная форма всплывающего тела: формируется канал высокотемпературной магмы (магмопроводник) и головное тело диапира. 3) Высота всплывания диапира зависит от реологических свойств окружающей коры: увеличение предела текучести на порядок (от 1 до 10 МПа) при снижении температуры ограничивает возможный уровень подъема до глубины 15-16 км. 4) Над осевой частью диапира в рельефе дневной поверхности формируется поднятие высотой около 750 м.

IIб. Результаты моделирования гравитационной неустойчивости в предположении о вязкой реологии вещества. Рассмотрена задача о тепловой конвекции в земной коре. Расчеты выполнены с использованием двумерной математической модели тепловой конвекции, в которой движение породы описывалось уравнениями Стокса. Распределение температуры находилось из уравнения переноса тепла. Ввиду малости изменений плотности от температуры в уравнениях Стокса использовалось приближение Буссинеска. Было выполнено преобразование исходной задачи в естественных переменных «скорость – давление» в переменные «функция тока - вихрь».

Рис.3. Геометрия области и краевые условия.

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

(1)

где (2)

Подстановка (2) в (1) дает уравнение четвертого порядка для функции тока:

(3)

Уравнение переноса тепла имеет вид:

(4)

Здесь x, z – оси координат, Vx, Vz – компоненты вектора скорости, Y - функция тока, W - вихрь, n - коэффициент кинематической вязкости, t – время, q – температура. Координаты масштабированы на ширину L и глубину H области соответственно. Скоростные переменные обезразмерены на характерную скорость c/H, где c - коэффициент температуропроводности. Остальные величины нормировались следующим образом: время - на H2/c, температура - на максимальное значение T1, коэффициент кинематической вязкости - на значение вязкости n0 при минимальной температуре T0. Число Рэлея, характеризующее взаимодействие подъемных сил и сил вязкости, вычислялось как:

,

где h0, r0 - коэффициент динамической вязкости и плотность при T0, Q = T1-T0 - характерная разность температур, g – ускорение свободного падения, a - температурный коэффициент объемного расширения. Использовалось линейное уравнение состояния: , где p, p0 - давление на глубине и на поверхности соответственно, b - коэффициент изотермической сжимаемости. Давление, соответствующее механическому равновесию при постоянных температуре T и плотности r, меняется с глубиной по гидростатическому закону: p = p0-rgz.

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

где l - коэффициент теплопроводности, q – тепловой поток на нижней грани земной коры. Краевое условие на температуру на границе g1 записывалось в виде баланса тепловых потоков, в котором конвективная составляющая моделировалась как теплоотдача горизонтальной плиты при естественной конвекции. Для числа Нуссельта, характеризующего интенсивность теплообмена, использовалась зависимость следующего вида: .

В задачах геофизики вязкость претерпевает существенные изменения в полосе малой толщины земной коры. Большие градиенты вязкости и большие числа Рэлея являются характерными особенностями данной задачи. Для твердых пород вязкость описывается экспоненциальной зависимостью от температуры вида: , что требует специальных методов решения при рассмотрении задачи тепловой конвекции. В расчетах применялась аппроксимационная формула Франк-Каменецкого в виде степенной зависимости вязкости от температуры

Последняя сохраняет высокий градиент вязкости, но несколько снижает требования к численной схеме. Наличие вторых производных от вязкости в уравнении переноса вихря (1) обычно приводит к расходимости итерационных конечно-разностных алгоритмов. В настоящей работе при больших числах Рэлея (Ra > 106) и больших градиентах вязкости (n1/n0 < 10-3) вместо уравнений (1), (2), решалось уравнение четвертого порядка с переменными коэффициентами (3). Решение сформулированной задачи осуществлялось численным методом с использованием неявной итерационной конечно-разностной схемы со стабилизирующей поправкой. Недостатком такого подхода является сравнительно большое время расчета (особенно на неравномерных сетках), поскольку требует больших объемов вычислений. При рассмотрении задачи на неравномерной сетке, и в частности, в двух компонентной постановке «гранит – расплавленный базальт» представляется эффективным использование численных методов «частиц - в - ячейках».

Рассматривалось влияние числа Рэлея и начальных условий на эволюцию движения в земной коре. На рис. 4-6 представлен пример результатов численного моделирования конвективной неустойчивости с высоким градиентом вязкости n1/n0 = 10-3 и числом Ra = 5×105 при t = 0.2 M лет. В расчетах использовались следующие значения параметров: L = 60 км, H = 30 км, Li = 25 км, Lo = 35 км, T0 = 293 K, T1 = 823 K, r0 = 2800 кг/м3, p0 = 105 Па, l = 3.5 Вт/(м×K), c = 10-6 м2/сек, a = 3×10-5 1/K, b = 10-11 1/Па, q = 0.03 Вт/м2. Вычисления проводились на прямоугольной сетке с числом узлов 101x51. В качестве начальных данных бралось линейное распределение температуры по глубине z. Структура течения имеет характерный вид восходящего конвективного потока, по обеим сторонам которого образуются интенсивные вихревые течения, закрученные в направлении потока. С удалением от центра потока расположены вихри меньшей интенсивности, которые с увеличением числа Рэлея вытесняют основные вихри вверх и способствуют формированию системы более мелких вихрей. Распределение температуры имеет вид, характерный для эффекта диапиризма, при котором горячий температурный фронт проникает в более холодные слои земной коры.

Рис. 4. Распределение температуры.

Рис. 5. Структура течения: изолинии функции тока и поле скоростей.

Рис. 6. Распределение плотности.

Интенсивность проникновения и форма диапира определяются числом Рэлея, и в меньшей степени отношением вязкостей n1/n0. Распределение плотности повторяет форму диапира и показывает ее понижение на фронте тепловой волны. Решение задачи о конвективной неустойчивости в земной коре показало, что при больших числах Рэлея и градиентах вязкости возникает восходящий конвективный поток, формирующий тепловую волну к поверхности Земли, характерную для эффекта диапиризма. Для детального решения данной задачи, в частности, в двух компонентой постановке «гранит – расплавленный базальт» представляется перспективным использование метода «частиц - в - ячейках» на многопроцессорных вычислительных системах.

IIв. Создана новая 2D нестационарная модель геологических течений в приближении слабосжимаемой жидкости.

Рассмотрим систему уравнений, описывающую динамику слабосжимаемой жидкости, замкнутую уравнением состояния:

где r – плотность, – скорость, T – температура, p – давление, – вязкость, g – ускорение свободного падения, k – температуропроводность.

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

,

,

,

,

,

,

,

,

,

,

Для улучшения счетных характеристик численной модели введена замена переменных:

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

В итоге для безразмерных величин получена следующая система уравнений:

Из за большого объема этот материал размещен на нескольких страницах:
1 2 3