Кинетика движения остывающей магмы по вертикальной дайке постоянной толщины

1. Введение

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

В настоящей работе анализируется нижний, “однофазный” участок движения магмы. При этом считается, что канал, по которому магма поднимается к поверхности, представляет собой узкую вертикальную щель постоянного сечения (дайку), причем толщина w щели много меньше высоты Н дайки, т.е. расстояния от питающего очага до дневной поверхности. Таким образом, дайку условно можно представить в виде вертикальной длинной плоской ленты, соединяющей магматический очаг с поверхностью.

2. Вывод уравнений тепломассопереноса

Рассмотрим ситуацию в тот момент времени, когда вся дайка заполнена движущейся магмой, причем “однофазный” участок занимает значительную часть высоты дайки (обычно глубины питающих очагов трещинных извержений существенно больше 3 км). Характер движения однофазной магмы будет в большинстве случаев ламинарным (действительно, даже для w=10 м, скорости потока v=10 м/с и вязкости m= 10 Па*с получим Re=100< <Reкр=1300-20000, при менее напряженных режимах течения условие ламинарности выполняется с еще большим запасом). Это позволяет сравнительно просто вывести основные уравнения тепло- и массопереноса для дайки в любой ее внутренней точке.

Рис.1. К выводу уравнения теплового баланса

Очевидно, задачу можно рассматривать в двумерном варианте (т.е. пренебречь краевыми эффектами на кромках дайки-ленты). Поместив начало координат на осевой линии дайки в самом ее начале (т.е. в вершине питающего очага) и направив ось Z вертикально вверх (вдоль оси дайки), а ось Х - по направлению к ближайшей боковой стенке дайки, можно теперь рассмотреть произвольную точку с координатами (X,Z) (0£X£W/2, 0£Z£ L) и описать около нее пробный прямоугольник со сторонами D X и DZ так, чтобы точка (X,Z) была в центре этого прямоугольника (см. Рис.1).

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

(1)

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

(2)

Теперь проанализируем двумерные динамические уравнения Навье-Стокса при условии ламинарности течения несжимаемой магмы:

 | (3) (4)
где |

p- избыточное давление, m- динамическая вязкость магмы.

Уравнения (3) и (4) можно развернуть в виде:

(5)

(6)

Условие ламинарности потока и несжимаемости магмы позволяет довольно точно “рассортировать” члены уравнений (2), (5) и (6) по порядкам малости. Для этого нужно всего лишь учесть, что:

для любой из переменных
(уравнение непрерывности)

Результат такой “сортировки” позволяет переписать (2), (5) и (6) в виде:

нулевой порядок малости

первый порядок малости

второй и четвертый порядок малости

нулевой порядок малости

первый порядок малости

второй и третий порядок малости

нулевой порядок малости

первый порядок малости

второй порядок малости

(2а)

(5а)

(6а)

“Порядок малости” в этих уравнениях измеряется величиной W/2L и (поскольку W обычно составляет метры или доли метра, а L измеряется километрами) оценочно составляет 10-3-10-6, , т.е. различие во вкладе разных групп членов уравнений (2а), (5а) и (6а) очень велико. Это позволяет без существенной ошибки отбросить все члены первого и более высоких порядков малости и окончательно получить уравнения гидродинамики нестационарного течения магмы:

(5б)

(6б)

Оценить величину [image] без дополнительных предположений довольно трудно, поскольку для высокочастотных колебаний производная, т.е. скорость изменения даже очень маленькой величины vx может оказаться достаточно большой. Ограничимся для начала случаем квазистационарного режима течения, когда колебания отсутствуют и режим течения меняется достаточно медленно; это позволяет приравнять нулю первые слагаемые в обоих уравнениях. Уравнение (5б) при таком предположении свидетельствует о том, что изобары в потоке представляют собой горизонтальные линии, т.е. давление p является функцией только высоты z. Кроме того, горизонтальная составляющая скорости vx из обоих уравнений исчезла, т.е. скорость течения v (в первом приближении!) имеет только z-компоненту и также зависит только от одной переменной x.

Таким образом, получаем:

(7)

где А не зависит от х (т.е. является функцией только от z).

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

(8)

(Q - объемный расход магмы на единицу ширины дайки).

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

(А)

(Б)

(В)

(Г)

Условие (А) соответствует нормальному поведению пристеночного слоя жидкости, условия (Б) и (В) следуют из симметричной геометрии потока, условие (Г) для температуры задано через тепловой поток на стенку. Величина градиента температуры на границе зависит от времени t существования заполненной магмой дайки [image]. Это значит, что стационарным течение магмы в принципе быть не может, т.к. граничное условие само меняется со временем. Однако начиная с некоторого t изменение [image] становится малым, что позволяет основным переменным T и v “отслеживать” это изменение в квазистационарном режиме.

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

(9)

(10)

(11)

Эта система записана для отрезка [image] при граничных условиях (А) - (Г).

Уравнение (10) интегрируется; с учетом граничного условия (Б) получается:

или

(12)

(здесь А не зависит от х).

В системе есть масштабная величина [image], что позволяет перейти к системе с безразмерной переменной [image]:

(13)

(14)

справедливой на отрезке [image] при условиях:

(А’)

(Б’)

(Г’)

(Д’)

(уравнение (11), записанное в форме (Д’), перенесено к граничным условиям).

Система в общем виде аналитически не решается, поскольку в ней присутствуют довольно сложные функции температуры (в работе [1] принято

. (*)

Положив [image], где [image] - не известная нам функция (о ее виде пока можно лишь сказать, что она - четная функция от y), запишем ( с учетом граничного условия (А’):

(15)

где

Нормировочное условие (Д’) дает:

(16)

где

Н(1) - число, поэтому (16) позволяет получить связь между неизвестными m 0 и А:

(17)

Подставив (15) и (17) в (13), получим:

(18)

при граничных условиях (В’) и (Г’).

Дальнейший анализ этого уравнения проведем для нескольких частных случаев. Дело в том, что реально весь интервал подъема 0<z<h’ (h’ - высота, на которой магма начинает выделять летучие и уравнение (17) становится недействительным) можно разделить на три участка: 1) начальный участок движения горячей и сравнительно маловязкой магмы, где последним членом в (17) можно пренебречь по сравнению с предыдущим; 2) конечный участок, где основной механизм теплопритока - вязкий разогрев, описываемый последним членом (17), а дальнейшее охлаждение, описываемое вторым членом уравнения, практически отсутствует; 3) небольшой по протяженности промежуточный участок, где роль обоих членов приблизительно одинакова.

3. Анализ уравнения теплопереноса в дайке при условии полной компенсации теплопотерь на стенках за счет вязкого разогрева

(18а)

(0<y<1), граничные условия: Ty(0)=0; Ty(1)=-Fw/2.

Интегрируем (18а) с учетом (В’):

(19)

Граничное условие (Г’) позволяет определить m 0:

(20)

Учитывая (20), перепишем (19) в виде:

(19а)

Интегрируя (19а) еще раз, получим:

(21)

где [image] - величина, определяемая реальной зависимостью m (Т) для данного вида магм.

Итак, полное решение системы для любого сечения z=const на последнем участке дайки таково:

(21)

(22)

(23)

(24)

Обращает на себя внимание обратно пропорциональная зависимость перепада от объемного расхода. Это значит, что дифференциальное гидродинамическое сопротивление участка дайки, где доминирующим механизмом теплопритока является вязкий разогрев, всегда отрицательно (dA/dQ<0), т.е. такое движение принципиально неустойчиво.

Еще одна особенность решения - то, что температура Т0 на оси дайки “следит” за расходом Q и интенсивностью теплоотвода l F, но никак не связана с температурой в очаге или на предшествующем участке дайки. Это именно самонастраивающийся режим.

4. Анализ теплопереноса при простом охлаждении жидкой магмы

(25)

(0<y<1, граничные условия Ty(0)=0, Ty(1)=-Fw/2).

Участок движения без вязкого разогрева имеет некоторую протяженность по высоте, а значит, уравнение (25) должно выполняться не только в рассматриваемом сечении z=z0, но и в соседних более или менее удаленных от него сечениях. Поскольку граничное условие Ty(1)=-Fw/2 изменяется по координате z довольно слабо (производная dF/dz по модулю меньше геотермического градиента), то естественно предположить, что картина распределения температуры по х на близких высотах меняется либо параллельным переносом, либо пропорциональным умножением всего профиля на величину, близкую к единице. А поскольку относительное изменение температуры [T(1)-T(0)}/T(0) сравнительно невелико, то оба крайних варианта изменения профиля (и все промежуточные) практически равноценны. Поэтому для облегчения математических выкладок положим Tz=-a = const(x). Тогда уравнение (25) интегрируется и с учетом граничного условия (В’) дает:

(26)

Граничное условие (Г’) позволяет определить величину a:

(27)

и (26) можно переписать в виде:

(26а)

Нам осталось еще одно интегрирование (26а), а граничные условия исчерпаны, константу интегрирования определять не из чего. Однако есть простой выход. Речь ведь идет о начальном, нижнем участке дайки, а температура втекающей в него магмы равна температуре питающего очага Тп0. В то же время единственной причиной изменения температуры на этом участке является теплоотток в стенки, а этот отток на всем пути от z=0 до рассматриваемой высоты z можно подсчитать, зная F(z), и, значит, зная

(28)

(здесь [image]). Это позволяет проинтегрировать (26а) и получить:

(29)

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

(29)

(30)

(24)

Здесь, в отличие от предыдущего случая, температурный профиль связан с начальной температурой магмы в очаге и со всей предысторией ее остывания. Это принципиально влияет на устойчивость течения. Вычислим dA/dQ:

(31)

(здесь [image]). В режиме, далеком от застывания, определяющим в формуле вязкости (* ) является первое слагаемое, т.е. m (Т)» m Аexp(TA/T), тогда dm 0/dT0=-m 0TA/T02. С другой стороны, [image]. Подставляя все это в (31), получим:

(32)

В нижних участках дайки, где R(z) мало, dA/dQ>0, гидродинамическое сопротивление положительно, течение устойчиво. Граница устойчивости простирается от корня дайки до высоты zy, определяемой из уравнения:

(33)

(его нетрудно получить, приравняв нулю выражение в фигурных скобках (32)). Если левая часть (33) даже при zy=L продолжает оставаться меньше правой части, то движение будет устойчиво вплоть до самой верхней границы рассматриваемого отрезка дайки. Если zy<L, неустойчивость режима в верхней части дайки обеспечена. До поры до времени эта неустойчивость компенсируется перераспределением давления между нижней частью дайки и верхней. Поскольку это перераспределение происходит с конечной скоростью (равной скорости звука vз в магме), следует ожидать, что в дайке начнут генерироваться колебания (основная частота - не ниже fн» vз/L).

Недостаточность компенсации возникает, когда
, т.е.

(34)

В линейном (достаточно грубом) приближении F(z)=F0=const и T(L)»T0 можно получить:

(35)

При Q£ Qмин и фиксированном РS = [image] гидродинамическая неустойчивость потока в целом неизбежно приводит к срыву квазистационарного режима и появлению значительных ускорений.

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

1. Режим спокойного истечения. Его условие:

(33а)

2. Стабильный режим с колебаниями магматического столба в вертикальном направлении:

(36)

3. Срыв стабильного режима (истечение со значительным ускорением того или иного знака):

(35а)

5. Анализ уравнения теплопереноса при сравнимой роли вязкого разогрева и охлаждения.

Переходный участок между режимами, рассмотренными выше, очевидно, сравнительно невелик по координате z. Прежде чем анализировать полное уравнение (18), полезно отметить, что распределение теплопритока в двух рассмотренных выше “крайних” режимах различается довольно сильно, а именно: при охлаждении теплоприток пропорционален скорости, т.е. максимален при у=0 и равен нулю у стенки у=1; при вязком разогреве теплоприток (точнее, тепловыделение) пропорционален vvy, т.е. равен нулю и на оси, и у стенки, а максимума достигает на некоторой промежуточной координате. Это должно вызвать соответствующее изменение в распределении температуры Т(у), причем при переходе от режима охлаждения к режиму вязкого разогрева профиль распределения станет более плоским в центральной части и более крутым - в области максимального градиента скорости. Это значит, что, во-первых, считать [image] постоянным вдоль всего поперечного сечения мы уже не имеем права, а во-вторых, относительный вклад охлаждения и вязкого разогрева здесь меняется не только по z, но и по х.

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

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

6. О приближенных формулах вязкости.

Попробуем исследовать качественную картину распределения скорости потока по сечению дайки. Для заданного расхода согласно (15) и (16) можно записать:

(15а)

(16а)

где [image] не известная нам функция, входящая в соотношение m =m 0f(y), m 0 - вязкость магмы в центре сечения дайки.

Рассмотрим три варианта аналитического выражения для f(y).

1. f(y)=1. Такая функция соответствует постоянству (или очень малому изменению) температуры магмы по сечению.

из (16а) получаем:

(37)

и из (15):

(38)

(39)

Диссипация энергии происходит также неравномерно:

(40)

2. f(y)=exp(y2/y02). Такая функция соответствует вязкости, плавно увеличивающейся по направлению от центральной плоскости дайки к ее боковой стенке. Уменьшая у0, можно постепенно перейти ко все более резкому изменению температуры по поперечной координате у.

(41)

(42)

(здесь [image] - интеграл вероятности)

(43)

(44)

3. f(y)=(1-y2)-n. Функция приближенно соответствует течению магмы при температурах, обеспечивающих “намерзание” магмы на стенках. Показатель степени n характеризует толщину почти неподвижного “намерзшего” на стенках слоя магмы.

(45)

(46)

n | 1 | 2 | 3
H(1) | 8/15 | 16/35 | 128/315

(47)

(48)

Рис.2.

На рисунке 2 показано распределение диссипации энергии для этих трех функций f(y). Как видно из рисунка, кривые тепловыделения за счет вязкого разогрева качественно довольно похожи. Определяющим фактором в реальном процессе является не качественное, а количественное значение Nдис, т.е. в сущности величина вязкости m 0 в центре сечения дайки. Если исключить экстремальный случай намерзания магмы на стенках, когда реальное “живое сечение” дайки значительно уменьшается (и это можно учесть соответствующим изменением w), то с точки зрения точности расчетов не будет большой ошибкой использование простейшей функции f=const (кривая 1 на рис.2). Таким образом, методика машинного счета, принятая в [1] и основанная именно на этом предположении, может быть признана достаточно корректной. Тем не менее, в свете результатов настоящей работы выявляется еще один “порог” (возникновение колебаний столба магмы), который в [1] не обнаружен (да и не мог быть обнаружен, т.к. там анализировалась внешняя, а не внутренняя задача тепломассопереноса). Использование критерия (36) без труда позволяет обозначить на графиках рис.6 в работе [1] режимы колебательного истечения.

Литература:

1. В.В.Титов, С.А.Федотов. Математическое моделирование подъема вязкой магмы по вертикальной дайке. Препринт ИАЭ-3470/16, 1981.