[ «ДИНАМИКА РОСТА КРИСТАЛЛА В ОЧАГАХ И КАНАЛАХ ВУЛКАНА ...»
МОСКОВСКИЙ ГОСУДАРСТВЕННЫЙ УНИВЕРСИТЕТ
имени М.В. Ломоносова
МЕХАНИКО-МАТЕМАТИЧЕСКИЙ ФАКУЛЬТЕТ
Кафедра гидромеханики
На правах рукописи
УДК 532.2:536.421.4
Горохова Наталья Владимировна
ДИНАМИКА РОСТА КРИСТАЛЛА В ОЧАГАХ И КАНАЛАХ ВУЛКАНА
Специальность 01.02.05 – Механика жидкости, газа и плазмы Диссертация на соискание учной степени кандидата физико-математических наук
Научный руководитель:
доктор физико-математических наук, член корреспондент РАН О.Э. Мельник
Научный консультант:
доктор геолого-минералогических наук, П.Ю. Плечов Москва, Оглавление Введение Глава 1. Математическая модель роста кристалла и е реализация Основные предположения................................. 1. Равновесная модель роста кристалла........................ 1. Неравновесная модель диффузионного роста кристалла........ 1. Обезразмеривание системы уравнений диффузии............. 1. Переход в подвижную систему координат................... 1. Аппроксимация уравнений................................ 1. Линеаризация уравнений................................. 1. Программная реализация.................................. 1. Сеточная независимость.................................. 1. 1.10 Выводы к главе 1......................................... Глава 2. Качественный анализ кристаллизации плагиоклаза Постановка задачи о росте кристалла плагиоклаза............. 2. Параметры модели........................................ 2. Анализ размерностей 2.3.................................... Равновесная кристаллизация плагиоклаза.
2. Аналитическое решение.................................. Задача Стефана (диффузионный рост кристалла).
2. Общие закономерности процесса кристаллизации............ Закономерности процесса кристаллизации в реакционноконтролируемых процессах................................ Выводы к главе 2......................................... 2. Глава 3. Рост кристалла плагиоклаза в остывающем магматическом расплаве Постановка задачи и основные предположения............... 3. Результаты в случае отсутствия растворения кристалла......... 3. Результаты в случае с учетом растворения кристалла........... 3. Выводы к главе 3......................................... 3. Глава 4. Рост кристалла плагиоклаза при подъме магмы по Постановка задачи и основные предположения............... Результаты расчетов по неравновесной модели............... Введение В недрах Земли при определнных условиях может образоваться силикатный расплав – магма. За счет сил плавучести магма начинает подниматься к земной поверхности и обычно скапливается перед извержением в приповерхностных коровых очагах, хотя возможен и непосредственный подъем из зоны магмогенерации. В коровых очагах часто происходит тепловое и механическое взаимодействие магм различного состава. В дальнейшем, в процессе извержения, магма поднимается из очага к поверхности. Различный состав первоначальных пород, процессы, протекающие внутри самой магмы, и внешние условия (например, изменение давления) могут значительно изменить химический состав магмы и е физические свойства.
Процессы смешения магм, кристаллизацию минералов из расплава и их растворение непосредственно в магматическом очаге или канале вулкана наблюдать невозможно. Это во многом понижает достоверность теорий, построенных на изучении природного и экспериментального материала, и делает их в первую очередь качественными, а не количественными. Поэтому становится необходимым использование математических методов моделирования как инструмента для изучения процессов, происходящих в вулканических системах.
Магма – сложная смесь большого числа компонент (отдельных атомов и стабильных комплексов), которые перераспределяются между различными минералами, образующимися при кристаллизации, или остаются в расплавной фазе. В работе изучаются процессы кристаллизации, происходящие в магме в условиях падения температуры (в очаге вулкана) или давления (при подъеме магмы по каналу). Динамика роста кристалла во многом определяется диффузией соответствующих компонент из расплава к его границе: кристалл не может расти, если в пограничном слое отсутствуют в достаточном количестве составляющие его компоненты.
А. Фиком в 1855 году [16]. Первый закон Фика устанавливает связь между потоком массы диффундирующего вещества и градиентом его концентрации в расплаве J D grad C, где D – коэффициент диффузии компоненты. Второй закон Фика есть следствие закона сохранения массы. Он определяет изменение концентрации во времени:
В случае постоянного коэффициента диффузии уравнение имеет вид:
Это уравнение имеет аналитические решения при нескольких типах граничных условий. В многокомпонентном расплаве диффузия одних компонент зависит от градиентов концентраций других, поэтому для описания многокомпонентной диффузии вводится диффузионная матрица D [42]:
где N – число компонент (концентрация N-ой компоненты является зависимой переменной). Подтверждением этого является существование возвратной диффузии или диффузии вверх по потоку (градиенту концентрации, uphill diffusion), наблюдаемой как в экспериментах [8,14,15,52], так и при исследовании петрологических данных [24].
Определение полной диффузионной матрицы для конкретной среды, содержащей большое число компонент, представляется весьма сложным. Для многокомпонентной диффузии ионов в кристаллах выведены формулы, определяющие коэффициенты диффузионной матрицы в зависимости от концентраций и зарядов ионов [29]. В [53] коэффициенты диффузионной матрицы записаны с использованием атомной мобильности компонент и скорости решетки («lattice velocity»). В магматических расплавах диффундирующими единицами являются оксиды. Для некоторых многокомпонентных расплавов коэффициенты диффузионной матрицы измерены экспериментально [6,8,9,50,51]. В [33] модели многокомпонентной диффузии, использующие полную диффузионную матрицу (как, например, выведенную в [29]), протестированы с использованием экспериментальных данных. Экспериментально показано, что они могут зависеть от концентраций компонент, температуры и давления системы, а также от вязкости расплава [3,9,50,51]. Такие эксперименты весьма сложны и позволяют определить диффузионную матрицу только для конкретного состава смеси компонент. В связи с этим предпринимаются попытки определять диффузионные потоки более простым образом.
Наиболее простая модель многокомпонентной диффузии основана на предположении двухкомпонентности расплава. В этом случае диффузионная матрица имеет единственный коэффициент, называемый эффективным бинарным коэффициентом диффузии. Это предположение используют многие авторы как при моделировании диффузии (например [30]) так и в экспериментах.
Эффективный коэффициент диффузии каждой компоненты в предположении о том, что диффузия этой компоненты определяется е собственным потоком массы, экспериментально может быть определен и для систем с бльшим количеством компонент [32,51]. Однако этот подход невозможно использовать для компонент, которые диффундируют вверх по потоку.
Для систем, состоящих из большего числа компонент, многие авторы предполагают диагонализируемость диффузионной матрицы или изначально используют диагональную матрицу [8,18]. Собственные векторы диффузионной матрицы определяют собственные компоненты, объединяющие исходные компоненты системы в соответствующих пропорциях. Собственные компоненты диффундируют независимо, их коэффициенты диффузии – собственные значения диффузионной матрицы. Это позволяет рассматривать диффузионный процесс в многокомпонентной системе как объединение независимых обменных процессов, описываемых уравнением Фика. Недостатком такого рассмотрения является то, что описываемое объединение компонент в группы – лишь способ интерпретации и визуализации измерений, но не имеет под собой физической основы [8].
Для описания многокомпонентной диффузии с возможной возвратной диффузией Y.Oishi для трехкомпонентной смеси предложена система уравнений, учитывающая перекрестную диффузию компонент [41]. Перекрестные члены в выражениях потоков массы возникают в случае различных скоростей диффузии компонент. В [41] проводится анализ решений системы в зависимости от соотношения коэффициентов диффузии компонент. В [40] проводится анализ модели Y.Oishi, получены критерии возникновения возвратной диффузии на примере конкретного расплава.
Проблема выбора диффузионной матрицы в настоящее время не решена. В диссертации модель [41] адаптирована для моделирования многокомпонентной диффузии в системе, состоящей из произвольного числа компонент.
В работе исследуется рост одиночного кристалла плагиоклаза – одного из основных породообразующих минералов. Плагиоклаз состоит из молекул двух типов: анортита (An - CaAl2Si2O8) и альбита (Ab - NaAlSi3O8). При понижении температуры или давления кристалл начинает расти за счет диффузии своих компонент (альбита и анортита) из расплава к границе кристалла. При этом другие компоненты расплава, наоборот, диффундируют от границы кристалла.
Этот минерал был выбран по той причине, что диффузия компонент внутри кристалла на несколько порядков медленнее, чем в расплаве. Поэтому распределение компонент внутри кристалла сохраняется на значительных временах после окончания кристаллизации, а анортитовый номер плагиоклаза – величина An=Ca/(Ca+Na) (содержание анортита в кристалле) – важный индикатор условий кристаллизации магмы. Для многих изверженных пород существуют детальные измерения этого параметра.
Моделирование неравновесной кристаллизации магматических расплавов (когда скорость диффузии компонент в расплаве оказывает существенное влияние на процесс роста кристалла) началось в 70-х годах прошлого века. Наиболее двухкомпонентной системе, предложена в [30]. Предполагается простая линейная зависимость концентрации компоненты в кристалле от ее концентрации в расплаве (их отношение называется коэффициентом распределения). Важность этой работы в том, что в ней предложены основанные на экспериментальных данных формулы для вычисления температуры равновесия кристалл-расплав и скорости роста кристалла. Температура равновесия и скорость роста определены в зависимости от состава расплава в пограничном слое и температуры системы. В работах [37,38] проводится исследование применимости предложенной в [30] формулы для вычисления скорости роста в сравнении с экспериментами.
Большой интерес представляет наблюдаемая в изверженных породах и при проведении экспериментов осцилляторная зональность (немонотонное изменение состава концентрации анортита в кристалле) химического состава ([27], рис.1), в кристаллах плагиоклаза. Существуют различные точки зрения на происхождение этой зональности. Классическим является представление, связанное с колебаниями температуры, вызванными магматической конвекцией в очаге извержения [22]. В [46] в двухкомпонентной системе моделируется кристаллизация, вызванная скачком температуры на границе кристалл-расплав.
Решаются одновременно уравнения теплопроводности и диффузии, из решения которых при этом находится состав кристалла. Процесс контролируется диффузией компонент и выделением скрытой теплоты кристаллизации. Показано, что в такой системе возможно возникновение периодически-зональных кристаллов. Проведено исследование в зависимости от безразмерных параметров (определяемых коэффициентами диффузии компонент и скрытой теплотой кристаллизации). В результате определены области параметров, в которых образуются кристаллы разных типов зональности. Недостатком работы является то, что при моделировании не учитывается зависимость состава кристалла и скорости роста от состава расплава и других параметров. Кроме того, как будет показано ниже, для характерных значений размеров ячеек магматического расплава, из которых растут кристаллы плагиоклаза, коэффициент температуропроводности на шесть порядков больше коэффициента диффузии, поэтому выравнивание температуры внутри ячейки на временах роста кристалла происходит мгновенно. Предложенный в [46] механизм не может объяснить зональности природных плагиоклазов.
Ниже, как и в ряде других работ, показано, что колебания состава кристалла могут быть вызваны конкуренцией процессов роста и диффузии компонент, тогда как температура системы меняется во времени монотонно. При моделировании зональности авторами рассматривается только двухкомпонентная система.
Рис.1: зональность в кристалле плагиоклаза [45], а – срез зонального кристалла плагиоклаза, б – зависимость состава кристалла (концентрации анортита, C1S ) от расстояния от границы кристалла.
возникновения периодической зональности в кристалле плагиоклаза. В [1] предполагается, что такая зональность в кристалле возникает в результате того, что зависимость скорости роста от концентрации альбита в приграничном слое расплава имеет вид трехзначной функции. В результате, может происходить экспериментального подтверждения. Эксперименты Kirkpatrick et al. [23], Hammer and Rutherford [19], Couch et al. [12] и других исследователей показывают, что зависимость скорости роста кристалла от переохлаждения расплава имеет куполообразную форму.
В [30] показано, что если скорость роста кристалла – однозначная функция состава расплава, то концентрации компонент в кристалле изменяются монотонно. Предполагается, что осциллирующее зонирование в кристалле может возникнуть, если для одного и того же состава расплава при постоянной температуре скорость роста кристалла может иметь несколько различных значений, в результате чего возможны переходы между режимами быстрого и медленного роста кристалла. Также было получено немонотонное изменение состава кристалла в случае, когда скорость роста – монотонно убывающая ступенчатая функция времени. Такой закон роста в статье предполагается более вероятной причиной возникновения сложной зональности, чем многозначность функции скорости роста, однако экспериментального подтверждения данный закон роста также не нашел.
В [48] для вычисления скорости роста предлагается формула, учитывающая влияние шероховатости растущей грани кристалла на скорость его роста. При росте кристалла грань заполняется, и скорость роста скачком изменяется от режима, соответствующего шероховатой стенке, к более медленному режиму, соответствующему гладкой стенке. Далее в режиме гладкого роста из-за обогащения расплава скорость роста начинает увеличиваться и скачком переходит на режим шероховатого роста. Модель объясняет колебания состава кристалла, однако период колебаний определяется положением точек перехода между режимами, физическая обоснованность которых не ясна.
В [25,26] на основе модели роста кристалла [30] проводится исследование устойчивости решения линеаризованного в окрестности стационарного решения уравнения диффузии. При анализе получено, что периодическая зональность возможна только при значениях коэффициента распределения компоненты между кристаллом и расплавом (отношения концентраций компоненты в кристалле и расплаве) меньших единицы. Такие значения коэффициента распределения не являются типичными для роста кристалла плагиоклаза из магматического зональности.
Считается, что коэффициенты диффузии анортита и альбита определяются коэффициентами диффузии Ca и Na в расплаве. В [69] по результатам экспериментов с образцами изверженных пород были оценены коэффициенты диффузии основных элементов, составляющих магматический расплав. На основе этих оценок определены элементы, определяющие диффузию анортита и альбита.
Подобный подход существенно упрощает реальную диффузию в магматическом расплаве, где основные диффундирующие единицы – оксиды соответствующих элементов. В [31,32] предлагается последовательное термодинамическое описание многокомпонентной диффузии в модельных магматических расплавах, однако, в настоящее время отсутствует необходимый экспериментальный материал для описания системы плагиоклаз-расплав.
Из анализа опубликованных моделей можно заключить, что моделирование многокомпонентной диффузии (в системах, состоящих больше, чем из двух компонент) проводилось многими авторами для задач, не связанных с исследованием роста кристаллов из многокомпонентных расплавов. В тех работах, где проводится моделирование диффузионного роста кристалла, рассматривается только двухкомпонентная система, хотя в реальности, магматический расплав состоит из большого числа компонент. Авторами переменном коэффициенте распределения для определения состав кристалла в зависимости от состава расплава. В работе Wang Y. и Merino E. [49] при моделировании осцилляторной зональности в кристалле хотя и записываются уравнения диффузии для многокомпонентной смеси (больше двух компонент), но используется диагональная диффузионная матрица (следовательно, возвратная диффузия не учитывается). Предложенная в [49] модель не учитывает сложную зависимость состава кристалла от состава расплава и влияние на процесс кристаллизации условий, в которых растет кристалл (изменение температуры расплава, давления и т.д.).
Целью данной работы является построение модели кристаллизации многокомпонентного расплава и исследование процесса кристаллизации на примере плагиоклаза. Одной из важных задач является объяснение, с помощью предложенной в работе модели кристаллизации, возникающей в кристаллах зональности. Научная новизна состоит в том, что, в отличие от предыдущих работ, при исследовании кристаллизации плагиоклаза используется модель многокомпонентной диффузии, учитывающая влияние перекрестных членов (модель [41]). Состав кристалла плагиоклаза впервые определяется с помощью экспериментально-полученных зависимостей. Скорость роста кристалла при этом однозначно определяется составом расплава.
Научная значимость работы состоит в том, что построены модели равновесного и неравновесного (диффузионного) роста кристалла из произвольного расплава с учетом перекрестных членов. Впервые проанализировано влияние граничного условия в неравновесной модели на процесс кристаллизации на примере кристаллизации плагиоклаза. Результаты численного моделирования показали, что обычно используемое предположение о существовании на границе локального равновесия и по составу, и по температуре (условие Стефана) не дает объяснения явлениям, наблюдаемым при экспериментах и в реальных кристаллах. Тем не менее, это условие может применяться при моделировании процессов, близких к равновесным. Показано, что зональность в кристалле может возникать и при однозначно-определяемой скорости роста.
Практическая значимость работы в том, что с помощью предложенных в ней моделей при использовании реальных зависимостей состава кристалла от состава расплава (и других параметров) можно реконструировать условия роста кристаллов в магматическом расплаве, тем самым получая уникальную информацию об условиях внедрения и подъема магмы.
Методом исследования в неравновесном случае является построение системы уравнений в частных производных, описывающей многокомпонентную диффузию в расплаве. Производится построение конечно-разностной краевой задачи, решаемой численно методом матричной прогонки. В равновесном случае система уравнений может быть численно решена явным образом. При использовании модельных зависимостей, задающих состав кристалла и температуру равновесия, получено аналитическое решение задачи о равновесном росте кристалла.
Основные результаты диссертации отражены в публикациях [60-70].
Апробация результатов была проведена на следующих конференциях:
1. Моделирование диффузионного роста кристалла плагиоклаза (устный доклад) // Конференция-конкурс молодых ученых НИИ механики МГУ, 14октября 2009 г., Москва.
2. Моделирование диффузионного роста кристалла плагиоклаза (устный доклад) // Ломоносовские чтения, 16-23 апреля 2010 г., Москва.
3. Модель диффузионного роста кристалла: механизм образования зональности в кристалле плагиоклаза (устный доклад) // 16 школа-семинар «Современные проблемы аэрогидродинамики», 6-16 сентября 2010 г., Сочи.
4. Модель диффузионного роста кристалла (устный доклад) // Конференцияконкурс молодых ученых НИИ механики МГУ, 13-15 октября 2010 г., 5. Модель диффузионного роста кристалла плагиоклаза (устный доклад) // Всероссийская научная школа молодых ученых «Механика неоднородных жидкостей в полях внешних сил», 30 ноября - 2 декабря 2010 г., Москва.
6. Модель диффузионного роста кристалла плагиоклаза (устный доклад) // X Всероссийский съезд по фундаментальным проблемам теоретической и прикладной механики «Актуальные проблемы механики», 24-30 августа 2011 г., Нижний Новгород.
7. Модель роста кристалла плагиоклаза при его подъеме по каналу вулкана (устный доклад) // Конференция-конкурс молодых ученых НИИ механики МГУ, 12-14 октября 2011 г., Москва.
8. Модель роста кристалла плагиоклаза при подъеме магмы по каналу вулкана (устный доклад) // XIX Международная конференция студентов, аспирантов и молодых ученых «Ломоносов», 9-13 апреля 2012 г., Москва.
9. Numerical simulation of plagioclase rim growth during magma ascent at Bezymianny volcano, Kamchatka (poster; the second author: O. Melnik) // European Geosciences Union, General Assembly 2012. 22-27 April 2012, Austria, Vienna.
10.Горохова Н.В. Исследование влияния параметров модели диффузионного роста кристалла на процесс кристаллизации плагиоклаза (устный доклад) // Конференция-конкурс молодых ученых НИИ механики МГУ, 8-9 октября 2013 г., Москва.
неравновесного роста кристалла в общем случае. Предполагается, что неравновесный рост кристалла контролируется диффузией его компонент к границе кристалл-расплав. Модель диффузионного роста кристалла основана на системе уравнений многокомпонентной диффузии, предложенных в [41].
Описаны методы решения полученных уравнений и алгоритм программной реализации.
В главе 2 построенные в главе 1 модели применены для исследования кристаллизации плагиоклаза, вызванной линейным падением температуры расплава, как в равновесном, так и в неравновесном случае. При определении коэффициентов диффузии, состава расплава и температуры равновесия использованы модельные зависимости. Исследовано влияние параметров задачи на процесс роста кристалла. Показано, что возможны режимы, когда температура расплава становится ваше равновесной температуры и, несмотря на падение температуры, происходит растворение кристалла. Объяснены причины возникновения зональности в кристалле плагиоклаза при линейном падении температуры. Получено, что существуют области значений параметров, когда образуются зональные кристаллы.
В главе 3 также рассматривается рост кристалла плагиоклаза в остывающем расплаве, но используются реальные зависимости параметров модели от процесса кристаллизации. Коэффициенты диффузии – экспоненциальные функции температуры, а состав кристалла и температура равновесия кристалл-расплав определяются эмпирическими функциями. Исследуется влияние скорости остывания и состава расплава, а также коэффициентов диффузии на процесс кристаллизации. Получен критерий возникновения зональности в кристаллах и найдена связь между шириной зон разного состава и режимом охлаждения.
В главе 4 моделируется рост кристалла плагиоклаза при подъеме магмы по каналу вулкана. Состав кристалла и равновесная температура определяются с помощью пакета ПЕТРОЛОГ [13], разработанного на геологическом факультете МГУ, в зависимости от состава расплава и давления. Исследуется влияние условий подъема на процесс кристаллизации. Результаты расчетов сравниваются с реальными кристаллами плагиоклаза из извержений вулкана Безымянный (Камчатка). В результате исследования реконструированы условия подъема магмы для отдельных извержений этого вулкана.
Глава Математическая модель роста кристалла и е реализация Основные предположения.
1. Для описания процесса кристаллизации рассматривается одномерная система, состоящая из кристалла и расплава, с подвижной границей между ними – фронтом кристаллизации (рис.1.1). Система координат вводится так, что центру кристалла соответствует плоскость (предполагается симметрия относительно этой плоскости). На подвижной границе x s между расплавом и кристаллом выполнены законы сохранения массы. На внешней границе предположение справедливо, когда с правой стороны существует аналогичная ячейка с растущим кристаллом, так что на границе ячеек отсутствует переток массы. Конечно, такая геометрия системы не отражает в полной мере сложности исследуемого явления, но позволяет качественно изучить влияние определяющих параметров на процесс роста кристалла.
Рассматривается рост кристалла, состоящего из K S компонент, в расплаве, состоящем из K m компонент ( K S K m ). Начальное распределение концентрации в расплаве для каждой компоненты задается ее общим содержанием в системе, т.е.
Здесь Ci – концентрация i-й компоненты в расплаве. При этом, кристалл состоит из компонент с номерами от 1 до K S, а компоненты с номерами от ( K S 1) до K m входят только в расплав.
Рост кристалла начинается, когда температура расплава опускается ниже температуры равновесия. Когда температура системы выше равновесной температуры, то существует только расплав, когда она ниже – то могут одновременно существовать и расплав, и кристаллическая фаза. Температура равновесия и состав кристалла могут зависеть от состава расплава, давления и других параметров. Способы их определения описаны в следующих главах.
Предполагается, что время диффузии в кристалле на несколько порядков превышает время диффузии в расплаве, поэтому диффузией компонент в кристалле можно пренебречь. Это предположение верно, например, для кристалла плагиоклаза [11]. В [17] коэффициенты диффузии отдельных элементов в кристалле и коэффициент совместной диффузии ионов CaAl-NaSi, определяющий совместную диффузию анортита и альбита в кристалле плагиоклаза, оценены по результатам экспериментов. Плотность магматического расплава принимает значение коэффициент теплопроводности 1 Дж мKс (средние значения) [10]. Откуда коэффициент температуропроводности равен коэффициенты диффузии имеют порядок D0 ~1013 м2 с (характерное значение коэффициентов диффузии компонент). Поэтому выравнивание температуры в пределах ячейки можно считать мгновенным. В модели не учитываются возникающие в результате роста (или растворения) кристалла напряжения и изменение плотности, поскольку они мало влияют на рассматриваемый процесс кристаллизации. Предполагается, что плотности всех компонент одинаковы, а макроскопическое движение магмы вокруг растущего кристалла отсутствует.
Задача о кристаллизации будет рассматриваться в двух постановках:
равновесной и неравновесной. Критерием для выбора постановки может являться диффузионное число Пекле Ped, характеризующее соотношение между характерное время диффузии, tV L V0 – характерное время роста кристалла, V0 – характерное значение скорости роста кристалла. Если число Пекле мало ( Ped 1 ), т.е. t D tV, при изменении размера кристалла компоненты в расплаве мгновенно (по сравнению с характерным временем роста) перераспределяются за счет диффузии. В этом случае можно использовать равновесную модель, не учитывающую диффузионные процессы. Зависимости состава кристалла от состава расплава и других параметров получены в результате равновесных экспериментов. Если времена роста и диффузии сравнимы ( Ped ~ 1 ), то диффузию компонент необходимо учитывать, и задача о кристаллизации решается в неравновесной постановке. В неравновесной модели предполагается локальное равновесие по составу на границе кристалл-расплав, так что зависимость состава кристалла от состава расплава остается той же, что и в равновесном приближении.
Равновесная модель роста кристалла.
1. В равновесной постановке предполагается, что характерное время диффузии много меньше характерного времени роста (диффузия происходит мгновенно). Эта модель описывается следующей системой уравнений:
Здесь C iS – концентрация i-ой компоненты в кристалле, определяется как функция состава расплава и давления; M iS ( s ) - удельная масса i-й компоненты в слое кристалла толщиной s; Ci - концентрация i-й компоненты в расплаве; Ci0 общее содержание (концентрация) i-й компоненты в системе, равное начальной концентрации соответствующей компоненты в расплаве; s – положение границы кристалла; f i ( x ) – профиль концентрации i-ой компоненты в кристалле – (рис.
1.2).
Уравнения (1.2.1, 1.2.2) описывают зависимость состава кристалла и температуры равновесия от состава расплава и давления, уравнение (1.2.3) определяет функции f i ( x ), уравнение (1.2.4) вычисляет удельную массу i-ой компоненты в слое кристалла как интеграл от функции f i ( x ), уравнения (1.2.5, 1.2.6) описывают законы сохранения массы каждой компоненты в системе, уравнения (1.2.7, 1.2.8) есть стехиометрические уравнения для кристалла и расплава соответственно. В общем виде эти уравнения были сформулированы М.Я.Френкелем [59]. Они широко используются в петрологических программах моделирования, таких как ПЕТРОЛОГ [13] или COMAGMAT [2].
В процессе роста кристалла за счет выделения скрытой теплоты кристаллизации температура изолированной системы кристалл-расплав растет:
где Lh 3.5 105 J kg – скрытая теплота кристаллизации, cV 1200 J ( kg K ) – удельная теплоемкость (взяты характерные значения для кристалла плагиоклаза в андезитовом расплаве). В случае учета теплообмена, вязкой диссипации и других неизотермических процессов закон изменения температуры может быть более сложным.
Используя уравнения (1.2.3 – 1.2.6), можно найти C1,, CKm как функции s и, подставив их в уравнения (1.2.1, 1.2.2) получить систему уравнений (1.2.10):
В силу (1.2.7) из K S уравнений для концентраций компонент в кристалле будут независимы ( K S 1) уравнений. Таким образом, система (1.2.10) содержит K S независимых уравнений. Если задать давление и равновесную температуру, то можно вычислить профили концентрации всех компонент в кристалле.
Чтобы решить обратную задачу, необходимо задать, например, профиль первой компоненты в кристалле. Тогда, зная C1S при всех значениях s, из уравнения системы для этой компоненты, решая обратную задачу, можно найти давление как функцию s, а затем, из остальных уравнений – температуру равновесия и концентрации других компонент. При этом решение обратной задачи может быть не единственно (в качестве первой компоненты можно выбрать любую из компонент), возможен поиск решения методом минимизации отклонения измеренного состава кристалла от вычисленного.
Неравновесная модель диффузионного роста кристалла.
1. кристалла. Одномерная многокомпонентная диффузия описывается следующими уравнениями:
где Ji – поток массы i-й компоненты.
Потоки Ji находятся из соотношения Фика-Онзагера [42]:
или в матричной форме:
где J - вектор-столбец потоков массы компонент ( J1,, J Km 1 ), D - диффузионная матрица. При этом компонента K m является зависимой и вычисляется из второго уравнения системы (1.3.1).
В простейшем случае диффузионная матрица имеет диагональный вид, т.е.
Dii Di, Dij 0 при i j. При этом коэффициенты Di ( i 1,..., K m 1 ) являются коэффициентами диффузии компонент, а сами компоненты диффундируют независимо друг от друга. Тогда выражения для потоков массы примут следующий простой вид (первый закон Фика):
и при подстановке их в уравнения (1.3.1) будут получены уравнения диффузии Однако, если диффузионная матрица диагональна, то не учитывается перекрестная диффузия, возникающая при различных скоростях диффузии используется следующее выражение для потоков массы [41]:
То есть, диффузионная матрица примет вид:
Таким образом, в неравновесной постановке диффузия компонент в расплаве описывается системой уравнений, состоящей из (1.3.1) и (1.3.3). Эта система уравнений решается только в расплаве, т.е. в области с подвижной границей s x L. Начальные условия заданы уравнениями (1.1.1). Граничным условием на подвижной границе кристалла x s (рис. 1.3) являются условия сохранения массы на разрыве V Ci J i 0 : к границе, движущейся со скоростью V, притекает со стороны расплава поток массы i-й компоненты V Ci 0 J i 0, и в кристалл уходит масса V Ci S ( Ci и J i – скачки концентрации и потока i-й компоненты на границе между расплавом и кристаллом; в кристалле поток массы за счет диффузии отсутствует). Здесь C iS – концентрация в кристалле i-ой компоненты (кристалл состоит только из компонент 1,…, K S, т.е.
CiS 0 при i K S 1, K m ), V – скорость роста кристалла.
При x L выполнено условие равенства нулю потоков массы J i. Таким образом, граничные условия для системы ((1.3.1),(1.3.3)) имеют вид:
Обезразмеривание системы уравнений диффузии.
1. Для изучения общих закономерностей процесса кристаллизации магмы в результате моделирования запишем систему уравнений и граничные условия в безразмерном виде. Вводятся следующие параметры: характерное значение D (порядок) коэффициентов диффузии и характерное время диффузии tD L2 D0.
V ' Vt D L – безразмерная скорость роста кристалла, s ' s L – размер кристалла.
Уравнения записываются в безразмерных переменных: x x ' x L, t t ' t t D.
где Из того, что получим:
Таким образом, получена система уравнений в безразмерном виде:
и граничные условия:
В дальнейших преобразованиях рассматриваются уравнения, записанные только в безразмерном виде, поэтому штрихи опущены.
Переход в подвижную систему координат.
1. Так как граница кристалла движется с переменной скоростью V(t), то перейдем в подвижную систему координат: (t, x) (t1 t, x1 ( x s) (1 s)) (здесь t и x – безразмерные координаты, а индексом 1 обозначены новые координаты). Границей кристалла в новой системе координат становится x1 = 0, а внешняя неподвижная граница x=L переводится в дифференцирования:
Так как скорость роста кристалла зависит только от времени и не зависит от пространственной координаты, и ds Vdt, то Тогда в подвижной системе координат (t1, x1 ) получим систему уравнений Граничные условия в подвижной системе координат имеют вид (1.4.2), но используются выражения для потоков массы, записанные в подвижной системе координат:
Начальные условия в подвижной безразмерной системе координат запишутся в том же виде, как и в (1.1.1).
Аппроксимация уравнений.
1. Построенная в п.1.5 система уравнений решается численно средствами языка программирования FORTRAN. В этом разделе строится разностная схема математической модели. Для этого аппроксимируются записанные в подвижной системе координат безразмерные уравнения системы (1.5.1) и граничные условия (1.5.2). В дальнейших выкладках индекс 1, обозначающий подвижную систему, опущен.
В рассматриваемой области вводится разностная сетка [57]: по времени используется равномерное разбиение с шагом dt, по пространству – сильно сгущающееся к границе кристалла. Такое разбиение выбрано для того, чтобы диффузионный пограничный слой, в котором особенно сильно проявляется нелинейность, был покрыт достаточным количеством точек. Отрезок x [0,1] разбивается на N отрезков [ xi 1, xi ] так, что точка разбиения x0 соответствует границе x 0 (границе кристалл-расплав), точка разбиения x N соответствует границе x 1 (границе ячейки). Первый шаг по пространству h1 – самый маленький. Каждый последующий шаг в q раз больше предыдущего, где q – близкий к единице параметр (в расчетах принято, что q = 1.01): hi 1 q hi. Т.е.
Отсюда разбиение отрезка [0,1]:
При аппроксимации системы используется полностью неявная схема для обеспечения устойчивости. Так как распределение концентрации CKm находится из последнего уравнения системы (стехиометрического уравнения), если известны распределения концентраций других компонент ( Ci, i 1,, K m 1 ), то для построения разностно-краевой задачи достаточно аппроксимировать уравнения системы, описывающие диффузию первых K m 1 компонент, подставив в них выражение для концентрации CKm :
производные по t и по x:
Здесь верхний индекс 0 означает значение концентрации в предыдущий момент времени.
Поскольку выражения для потоков массы Ji нелинейны, то используется аппроксимация их производных вверх по потоку Где Здесь узлы сетки, в которых взяты значения концентраций в выражениях для потоков, зависят от направления потоков массы соответствующих компонент:
Таким образом, с помощью выкладок (1.6.3) – (1.6.8), аппроксимированные уравнения системы (1.5.1) для каждого внутреннего узла сетки i 1,, N можно записать в виде:
Для замыкания системы аппроксимируются граничные условия (1.5.2).
Выражение для потоков массы компонент в расплаве при х = 0 (i = 0):
для определения C j 0 используются формулы (1.6.8). Таким образом, граничное условие при i = 0 имеет вид:
При х = 1 (i = N) с помощью (1.6.4) получено граничное условие Таким образом, в результате аппроксимации системы (1.5.1) и граничных условий (1.5.2) получена нелинейная система разностных уравнений в следующем виде:
соответствуют уравнениям системы первых K m компонент; Ci C1 i,, CK m 1 i – векторы (столбцы): C j i, j 1,, K m 1, соответствуют значениям концентраций компонент в i-ом узле сетки. Концентрации при i 1 и i N 1 используются только для однообразия формул, их можно считать равными нулю. Система соответствующим условию (1.1.1). Во всех выкладках концентрация CKm определяется по формуле (1.6.2).
Линеаризация уравнений.
1. Нелинейная система уравнений (1.6.13) далее линеаризуется путем разложения в ряд Тейлора и отбрасывания квадратичных членов:
В итоге получена линейная система уравнений относительно неизвестных векторов (столбцов) вариаций Ci C1 i,, CKm 1 i, в качестве граничных условий выступают соответствующие уравнения при i = 0 и i = N. Матрицы использованием малого приращения концентрации dx (в расчетах принято, что dx 105 ).
Скорость роста и состав кристалла в общем случае переменные величины, которые могут зависеть как от условий кристаллизации (например, скорости определения состава и скорости роста кристалла будут описаны в следующих главах.
Программная реализация.
1. Полученная в п.1.7 краевая задача решается итерационным методом. В начале каждого просчета, соответствующего новому шагу по времени, имеется распределение C i, полученное на предыдущем шаге по времени и полагается Ci Ci. При каждой итерации решается система (1.7.1) методом матричной прогонки, описанным в [57]. В результате, на каждой итерации получается новое распределение Ci. Итерации проводятся до тех пор, пока i сумма станет меньше заданного, предполагается, что найдены новые распределения концентраций компонент на этом шаге по времени, и далее аналогично просчитывается следующий шаг по времени. Программа, решающая эту задачу, написана на языке FORTRAN, на рисунке 1.4 изображена ее блоксхема.
и другие (задающие усл-я крист-ии) Состав кристалла на границе C jS, скорость роста V Здесь, критерий прекращения расчетов зависит от задачи. В задаче, описанной в главе 2, это остывание расплава на определенное количество градусов (когда процесс роста кристалла приближается к равновесному и прекращаются возможные колебания состава кристалла); в главе 3 –достижение температурой расплава определенного значения, когда рост кристалла практически завершается. В задачах из главы 4 критерий остановки – достижение определенного давления.
Сеточная независимость.
1. Была исследована зависимость полученных результатов от количества точек в разбиении области вычисления (отрезка [0,1] ). На рисунке 1.5 построены зависимости скорости роста кристалла от времени при трех значениях N (числа узлов сетки): 150, 200, 250. Данный график приведен лишь в качестве иллюстрации сеточной независимости численного решения. Причины возникновения колебаний будут исследованы в следующих главах. Было получено, что при дальнейшем увеличении N результаты перестают зависеть от числа узлов (происходит полное наложение соответствующих линий). Поэтому, для уменьшения объема вычислений в расчетах принято, что N=250.
1.10 Выводы к главе 1.
Построена модель роста кристалла из многокомпонентного расплава.
Рассмотрен случай равновесного роста кристалла, когда состав кристалла полностью определяется содержанием компонент в расплаве и внешними условиями. Основным законом, определяющим равновесную кристаллизацию, является условие сохранения массы компонент в системе. Можно считать, что кристалл растет равновесно, когда характерное время диффузии много меньше характерного времени роста ( Ped 1 ). Именно равновесная постановка обычно используется в петрологии при моделировании кристаллизации минералов в магматическом расплаве. В случае неравновесного роста основную роль играют диффузия компонент из расплава к границе кристалла. Характер протекания процесса кристаллизации определяется соотношением скорости роста кристалла и скорости диффузии компонент. В результате, состав кристалла в каждый момент времени зависит от самого процесса.
Модель построена в общем случае и может быть применена для исследования кристаллизации различных минералов произвольного состава в условиях, например, остывания расплава или падения давления. Изменение внешних условий (температуры и давления) может быть задано произвольно.
Чтобы исследовать кристаллизацию конкретного минерала, необходимо задать соотношения, определяющие состав кристалла в зависимости от состава расплава (и, возможно, других параметров), а также скорость роста кристалла.
кристаллизации плагиоклаза (одного из основных породообразующих минералов) в магматическом расплаве, вызванной остыванием расплава или падением давления.
Глава Качественный анализ кристаллизации плагиоклаза.
В главе 1 построена математическая модель диффузионного роста произвольного кристалла. Кристаллизация может происходить как по причине охлаждения расплава, так и по причине падения давления в водонасыщенных магматических расплавах, поскольку снижение давление приводит к повышению равновесной температуры. При этом параметры модели – коэффициенты диффузии, состав кристалла, температура равновесия и скорость роста кристалла могут определяться самыми различными способами. Это делает возможным использовать предложенную модель для исследования роста различных кристаллов в произвольных условиях. Здесь возникает другая сложность:
определение перечисленных выше параметров. Эксперименты по определению состава кристалла и равновесной температуры в зависимости от состава расплава, давления очень сложны и с технической стороны (в силу высоких температур и давлений) и с точки зрения интерпретации результатов (поскольку процессы, происходящие в магме, зависят от очень многих параметров, в том числе неучтенных в предложенной модели).
В этой и последующих главах будет рассматриваться кристаллизация плагиоклаза. Во второй главе предполагаются упрощенные зависимости параметров модели от процесса кристаллизации, которые, тем не менее, могут использоваться в качестве «локальных» зависимостей (например, при малом изменений температуры). Проводится исследование зависимости процесса кристаллизации от вида граничного условия на границе кристалл-расплав, коэффициентов диффузии, и параметров, характеризующих законы изменения состава кристалла и температуры плавления.
Постановка задачи о росте кристалла плагиоклаза.
2. Исследуется рост одиночного кристалла плагиоклаза. Плагиоклаз состоит из молекул двух типов: анортита (An - CaAl2Si2O8) и альбита (Ab - NaAlSi3O8).
Поэтому предполагается, что кристалл состоит из двух компонент ( K S 2 ), а расплав – из трех компонент ( K m 3 ). Первой и второй компонентой являются анортит и альбит соответственно, третьей компонентой – остаточный расплав.
Система уравнений {(1.3.1),(1.3.3)} в неподвижной системе координат примет вид Уравнения диффузии решаются только для анортита и альбита в силу третьего уравнения системы (см. п. 1.3). Граничные условия (1.3.4) примут вид рассматривать только концентрацию анортита в кристалле.
расплава. Давление считается постоянным. Предполагается, что температура расплава линейно убывает с течением времени:
Здесь T0 [ K ] – начальная температура, VT [ K ч] – скорость остывания расплава.
Система уравнений решается методами, описанными в главе 1. Для е решения необходимо доопределить неизвестные величины и функции:
коэффициенты диффузии, состав кристалла (концентрацию анортита в кристалле) и скорость роста кристалла в зависимости от состава расплава и, возможно, других параметров.
Параметры модели.
2. Коэффициенты диффузии В этой главе предполагается, что коэффициенты диффузии компонент постоянны; в главе 3 будет рассмотрен случай, когда коэффициенты диффузии зависят от температуры. При этом D1, D2 ~ 1013 м2 с. Основным параметром, определяющим взаимную диффузию анортита и альбита, является отношение коэффициентов диффузии CD D2 D1, которое может принимать значения как большие 1, так и меньшие. Коэффициент диффузии анортита D1 и отношение CD могут принимать различные значения, коэффициент диффузии остаточного расплава считается малым по сравнению с коэффициентами диффузии анортита и альбита ( D3 1017 м2 с ) и в данном исследовании не варьировался.
Отметим, что диффузия в магматическом расплаве определяется движением катионов, а не отдельных молекул анортита или альбита, которые, однако, образуют устойчивые образования в кристалле, поэтому предложенная схема во многом является упрощенной, однако недостаток экспериментальных данных не позволяет заполнить диффузионную матрицу для системы катионов, составляющих магматический расплав. Это также является причиной, почему в расчетах коэффициенты диффузии варьируются в широких пределах.
Граничное условие на границе кристалл-расплав.
Классическим считается представление граничных условий в форме Стефана: вместо одного из K m 1 граничных условий (1.3.2) на границе температуре расплава, то есть предположение об отсутствии локального переохлаждения на границе растущего кристалла:
Тогда остальными K m 2 условиями на границе кристалл-расплав будут условия, выражающие законы сохранения массы на разрыве для любых K m компонент. Таким образом, граничные условия на границе кристалл-расплав имеют вид (при x=s):
Условие для ( K m 1) -ой компоненты позволит вычислить скорость роста кристалла:
Для задачи о кристаллизации плагиоклаза граничные условия примут вид:
и скорость роста кристалла определяется по формуле кристалла встраиваются в его решетку существенно быстрее, чем подходят в диффузионный пограничный слой, так что любое локальное переохлаждение мгновенно компенсируется переходом вещества в кристалл с изменением состава пограничного слоя и, соответственно, температуры равновесия. Таким образом, фактором, лимитирующим кристаллизацию, является диффузия компонент как более медленный процесс [54, 56].
Если же скорость встраивания частиц в растущий кристалл является лимитирующим фактором (осуществляется реакционный рост кристалла), то состав приграничного слоя может отличаться от равновесного с кристаллом, а равновесная температура не будет равна температуре расплава. Тогда в качестве K m 1 граничных условий на границе кристалл-расплав берутся условия сохранения массы компонент, а скорость роста кристалла определяется из дополнительного соотношения:
Здесь T T Teq – переохлаждение расплава, а Tu const – переохлаждение, при котором скорость роста достигает максимума V0 (фиксированная величина).
Teq – температура равновесия. Если температура выше температуры равновесия, кристалл растворяется, если ниже – растет. Зависимость скорости роста от переохлаждения представлена на рис. 2.1. При этом на величину скорости роста существенно влияет лишь переохлаждение расплава, а влияние температуры расплава пренебрежимо мало. Формула (2.2.6) для определения скорости выведена путем адаптации выражения, полученного из термодинамических соображений для магматических расплавов [20,35] ( Tu 62 K для плагиоклаза).
При этом при положительных переохлаждениях, когда температура расплава больше температуры равновесия, скорость роста кристалла равна нулю.
Во многих задачах необходимо учитывать растворение кристалла при температурах, больших равновесной температуры [47]. Тогда в вычислениях скорость роста полагается отрицательной при положительном переохлаждении. В этом случае считается, что при 0 T Tu скорость растворения по модулю равна скорости роста при таком же по модулю переохлаждении расплава. При T Tu скорость растворения максимальна и равна по модулю V0. Когда кристалл растворяется, состав расплава на границе с кристаллом определяется содержанием компонент в растворившемся слое кристалла.
Состав кристалла В простейшем случае состав кристалла плагиоклаза определяется долей анортита в расплаве на границе с кристаллом – величиной Xan C1 (C1 C2 ). В этой главе не учитывается зависимость состава кристалла от других параметров:
температуры, давления и доли остаточного расплава. Зависимость концентрации анортита в кристалле от величины Xan можно описать квадратичной функцией (рис. 2.2):
При малом содержании анортита в расплаве по сравнению с альбитом ( Xan мал) кристалл растет преимущественно из альбита и наоборот. Такое приближенное описание хорошо аппроксимирует зависимости, представленные в гл. 3, при незначительном изменении валового состава расплава (степени его кристаллизации).
Температура равновесия Равновесная температура также определяется долей анортита в расплаве на границе с кристаллом. Предполагается, что температура равновесия линейно зависит от Xan и определяется формулой Здесь Teq 0 – начальная температура равновесия, Xan0 – доля анортита в начальный момент времени, kT – параметр, задающий наклон прямой. При любом kT начальная температура расплава равна начальной температуре равновесия (расплав начинает кристаллизоваться из состояния равновесия).
В случае, когда зависимость равновесной температуры от Xan имеет более сложный вид, всегда можно считать зависимость (2.2.8) линейным приближением для малого изменения состава расплава (за небольшой промежуток времени).
Поэтому для исследования поведения процесса кристаллизации в зависимости от функции, определяющей температуру равновесия, на качественном уровне можно использовать линейную формулу.
Анализ размерностей 2. При решении задачи о кристаллизации находятся концентрации анортита и альбита в расплаве как функции расстояния от границы кристалла и времени, а также температура равновесия, концентрация анортита в кристалле на границе и размер кристалла s как функции времени. Далее, с точки зрения теории размерностей [58] будут определены независимые безразмерные параметры, характеризующие процесс кристаллизации.
В силу уравнений (2.2.1) процесс кристаллизации определяется коэффициентами диффузии D1, D2, D3 всех компонент и размером ячейки L (фиксированная величина в данном исследовании). Из начальных условий добавляются C10, C20 (начальный состав расплава), а из граничных условий (2.2.2) – скорость роста кристалла V, которая задатся своим максимумом V0 и переохлаждением Tu, при котором этот максимум достигается (ф-ла (2.2.6)).
Если используются граничные условия в форме Стефана, то эти два параметра (V и Tu ) не входят в систему определяющих. В формулу (2.2.6) входит также равновесная температура, которая определяется коэффициентом kT, начальными температурой T0 и содержанием анортита Xan0. Поскольку рост кристалла происходит в результате охлаждения расплава, то дополнительный параметр – скорость остывания VT.
При постоянных коэффициентах диффузии от температуры расплава явно процесс кристаллизации не зависит, а влияние оказывает только переохлаждение расплава. Поэтому начальную температуру T0 можно исключить из набора определяющих параметров. Таким образом, при моделировании кристаллизации плагиоклаза определяющими параметрами для всех искомых величин являются:
Размерности определяющих параметров:
Размерности искомых величин:
Из них размерно-независимыми являются L, D1 и kT. Если выразить характерное время t0, то с помощью них можно составить следующие безразмерные комбинации определяющих параметров:
Следовательно, в силу -теоремы [58] все безразмерные искомые величины есть функции этих комбинаций и безразмерных параметров:
(безразмерных), своя для каждой искомой величины. В этом исследовании начальный состав расплава фиксирован (в формуле соответствующие величины выделены), поэтому при записи формулы параметры C10, C20 можно опустить:
ищутся функции g для конкретного начального состава.
Таким образом, в случае реакционного роста кристалла (скорость роста кристалла определяется по формуле (2.2.6)) на характер протекания процесса кристаллизации влияют пять параметров: CD, C D 3, Ped, Tu ' и VT ' :
Параметры CD, C D 3 и Ped определяются коэффициентами диффузии анортита и альбита и отвечают за соотношение скоростей диффузии компонент и остаточного расплава, и скоростью роста кристалла. Если зафиксировать коэффициенты диффузии компонент и максимум скорости роста кристалла V0, то искомые функции определяются двумя параметрами (формула (2.3.2)). Скорость остывания VT и коэффициент kT меняются при этом независимо друг от друга.
При использовании граничных условий в форме Стефана (осуществляется диффузионный рост кристалла) скорость роста вычисляется из законов сохранения массы компонент на границе кристалл-расплав. В этом случае параметры V0 и Tu не являются определяющими и в формуле (2.3.1) не будет соответствующих им безразмерных комбинаций Ped и Tu '. Искомые величины будут тогда определяться тремя безразмерными параметрами:
При фиксированных коэффициентах диффузии компонент характер протекания процесса определяется единственным параметром VT ', т.е. зависит только от отношения скорости остывания VT и углового коэффициента kT :
Равновесная кристаллизация плагиоклаза.
2. Аналитическое решение.
При использовании описанных выше упрощенных зависимостей состава кристалла и температуры равновесия от состава расплава задачу о равновесной кристаллизации плагиоклаза (п.1.2) можно решить аналитически. Предполагаем, что кристалл растет равновесно из расплава с содержанием анортита и альбита C10 и C 20 соответственно и имеет начальный размер s 0 (положение границы кристалла). Система уравнений (1.2.1) – (1.2.8) для плагиоклаза примет вид:
Из системы, состоящей из уравнений (2.4.1), (2.4.2), (2.4.3)-(2.4.7), можно найти составы кристалла и расплава в зависимости от положения границы s.
Температура равновесия как функция s находится из уравнения (2.4.3).
Таким образом, из системы {(2.4.1), (2.4.2), (2.4.4)-(2.4.7)} ищем функции C1 ( s ), C2 ( s ), C1S ( s ) и Xan( s ), для чего продифференцируем уравнения (2.4.6) по переменной s (подставив M 2 S ( s ) s M 1S ( s ) ):
Из (2.4.5) следовательно, получим:
Сложив эти два уравнения и упростив второе, получим в итоге систему уравнений на функции C1 ( s ), C2 ( s ) :
В первом уравнении системы (2.4.8) можно сделать замену переменных X C1 C2, тогда оно примет вид:
Следовательно, Или, окончательно Константа А находится из условия, что C1 s C10 и C2 s C20 :
Таким образом, Подставляя выражение (2.4.9) во второе уравнение системы (2.4.8), получим дифференциальное уравнение на C2 ( s ) :
Преобразовав, получим:
Для решения этого уравнения (уравнения Бернулли) необходимо сделать вид Однородное решение уравнения (2.4.10):
подстановке в (2.4.10) получим уравнение на функцию B ( s ) :
Следовательно, И тогда Поскольку C2 s C20, следовательно, То есть окончательную формулу:
Из (2.4.9) имеем откуда Откуда, используя (2.4.1), найдем состав кристалла:
Таким образом, получено полное решение равновесной задачи:
Или в безразмерном виде ( s ' s L, s0 ' s0 L ):
Температура равновесия в зависимости от размера кристалла находится из уравнения (2.4.3): Teq Teq 0 kT ( Xan Xan0 ).
На рисунке 2.3 построены графики зависимости найденных функций от размера кристалла (положения границы кристалл-расплав) для C10 C20 0. (т.е. Xan0 0.5 ), s0 ' 0.1 :
Вне зависимости от закона изменения равновесной температуры в равновесном случае концентрации анортита и альбита в расплаве убывают в процессе роста кристалла. Относительное содержание анортита в расплаве (Xan) и, соответственно, концентрация анортита в кристалле также убывают. Кристалл растет до тех пор, пока анортит и альбит полностью не вычерпываются из расплава. Температура плавления при этом может как возрастать (при kT 0 ), так и убывать (при kT 0 ).
Задача Стефана (диффузионный рост кристалла).
2. Общие закономерности процесса кристаллизации.
В неравновесном случае сначала рассмотрена задача Стефана: задание равенства равновесной температуры температуре расплава в качестве граничного условия. Такое условие, в принципе, может иметь место при близкой к равновесной кристаллизации (п. 1.1).
Далее перечислены характерные значения величин, используемых в этом разделе. Кристаллизация вызвана остыванием расплава со скоростью VT 5 K ч 0.0014 K c на 200 К с начальной температуры T0 1100 K (для определенности). Поскольку коэффициенты диффузии постоянны, а температура равновесия и состав кристалла не зависят явно от температуры расплава (зависят только от его состава), то начальная температура может быть выбрана произвольно. Начальный состав расплава C10 C20 0.2, т.е. Xan0 0.5. Размер ячейки L 104 м, начальный размер кристалла s0 105 м (для возможности растворения кристалла). Коэффициент диффузии анортита C D 0.5, что соответствует реальным значениям коэффициентов диффузии анортита и альбита при этой температуре. Коэффициент диффузии остаточного расплава D3 1017 м2 с. Коэффициент kT варьируется в достаточно широком диапазоне – от -100 К до 1000 К. Таким образом, при постоянных коэффициентах диффузии протекание процесса кристаллизации определяется единственным безразмерным параметром VT ' (п. 2.3), принимающим в расчетах значения от 0.14 до 14.
На рисунке 2.4 показаны характерные профили концентраций анортита С1 и альбита С2 в расплаве. В силу того, что кристалл при данных условиях растет преимущественно из анортита, концентрация анортита в расплаве убывает сильнее, чем альбита. Поскольку альбит диффундирует медленнее, градиенты его концентрации больше. При этом наблюдается диффузия альбита вверх по потоку – увеличение концентрации альбита в расплаве с течением времени за счет влияния перекрестной диффузии.
Рис.2.4: профили концентраций анортита (а) и альбита (б) в расплаве На рисунке 2.5 представлены результаты расчетов при разных значениях VT '. Зависимость относительного содержания анортита в расплаве от времени (рис. 2.5-г) будет линейной, поскольку температура равновесия равна температуре расплава, убывающей с постоянной скоростью.
При анализе графиков на рис. 2.5-а,в можно заключить, что при увеличении безразмерной скорости остывания VT ' (когда она положительна) концентрация анортита в расплаве на границе с кристаллом быстрее убывает со временем, что следует из меньшего угла наклона соответствующей прямой на рис. 2.5-в. При больших значениях VT ' график Anm (t ) принимает положение, близкое к вертикальному. В этом случае компоненты кристалла практически мгновенно полностью вычерпываются из расплава, и в дальнейшем решение описывает нефизичную ситуацию, когда концентрации компонент в расплаве отрицательны.
Полное вычерпывание анортита и альбита из расплава происходит при значениях параметра, бльших 0.3. При уменьшении VT ', наоборот, график зависимости концентрации анортита в расплаве от времени становится более пологим. При отрицательных VT ' происходит смена режима и кристалл растворяется.
В силу (2.3.4) на процесс кристаллизации влияет только отношение скорости остывания VT и углового коэффициента kT, решение при отрицательном угловом коэффициенте и положительной скоростью остывания совпадет с решением для таких же по модулю положительного kT и отрицательной скорости остывания (расплав нагревается). Поскольку в природных условиях кристаллизация обычно вызвана остыванием магматического расплава (при постоянном давлении), в этом исследовании скорость остывания предполагается положительной, т.е. отрицательные значения VT ' не рассматриваются. Таким образом, для моделирования можно зафиксировать некоторое положительное значение углового коэффициента и варьировать только скорость остывания (при неизменных коэффициентах диффузии).
На рисунке 2.6 показано сравнение расчетов при граничных условиях, заданных в форме Стефана, с аналитическим решением равновесной задачи (2.4.11). Здесь оба параметра, отвечающих за диффузию компонент ( CD D2 D1 и CD 3 D1 D3 ) варьируются. Поскольку решение зависит от отношения kT и VT, зафиксируем kT, а VT будем варьировать.
В результате расчетов получено, что при любых значениях параметров кристалл растет монотонно. Параметр C D 3 (и, соответственно, коэффициент диффузии анортита D1 ), определяющий диффузионное время, не влияет на решение. На процесс кристаллизации оказывают влияние только соотношение коэффициентов диффузии компонент и безразмерная скорость остывания (параметр приближаются к аналитическому решению с увеличением отношения CD (рис.
2.6-б). При (т.е. когда коэффициент диффузии альбита больше коэффициента диффузии анортита), от величины CD результат практически не аналитическому равновесному решению при уменьшении скорости остывания.
Таким образом, расчеты показывают возможность использования граничных условий в форме Стефана в процессах близких к равновесным (п. 1.1).
Закономерности процесса кристаллизации в реакционноконтролируемых процессах.
В реакционно-контролируемых процессах кристаллизации (п.2.1) скорость роста кристалла определяется из дополнительного соотношения (2.2.6) в зависимости от переохлаждения расплава. В качестве граничных условий на границе кристалл-расплав используются законы сохранения массы каждой компоненты.
Используемые в расчетах значения величин:
постоянные величины:
T0 1100 K (для определенности), C10 C20 0.2.
варьируемые величины:
D1 1014,1013,1012 м2 с, CD 0.01, 0.1, 0.5, 10, VT 1, 2, 5, 10, 20 K ч, Этому соответствуют следующие значения безразмерных параметров:
Расчеты проводились как для случая без учета растворения кристалла, так и для случая с растворением. При этом при разных значениях варьируемых параметров могут быть получены различные сценарии протекания процессов кристаллизации. При одних значениях кристалл монотонно растет, при других – наблюдается периодический рост кристалла со сменой периодов роста и растворения (отсутствия роста). При этом профиль анортита в кристалле имеет ступенчатый вид. Также существуют значения параметров, при которых кристалл сначала растет, а затем начинает растворяться, но при этом профиль анортита в кристалле – монотонная функция.
положительном переохлаждении скорость роста равна нулю.
периодического роста представлен на рис. 2.7. При численном моделировании получено, что периодический рост кристалла наблюдается только при C D 1.
Здесь показаны результаты расчетов для C D 0.01, CD3 105, Ped 0.1, Tu ' 3.2 и VT ' 0.014 (что соответствует, например, D1 1012 м2 с, VT 1 K ч и kT 200 K ). Рост кристалла начинается, когда температура расплава становится ниже равновесной температуры. В процессе роста концентрации анортита и альбита в расплаве в пограничном слое убывают, причем концентрация анортита уменьшается быстрее, поскольку для данных условий в кристалл уходит преимущественно анортит. Изменение анортитового номера Xan на границе соответственно, переохлаждения. В результате, наступает момент, когда температура равновесия становится ниже температуры расплава и рост кристалла прекращается.
прекращение роста Рис.2.7: общий характер протекания процесса кристаллизации (а), зависимость концентраций анортита C1m и альбита C2m на границе кристалла (б, кривые 1,2), безразмерного переохлаждения T / TU (в) и безразмерной скорости роста кристалла VL / D1 (г) от времени.
Пока кристалл не растет, происходит перераспределение анортита и альбита в диффундирует быстрее, поэтому, в то время как его концентрация в пограничном слое достигает практически постоянного значения, концентрация альбита продолжает расти. Это вызывает дальнейшее понижение равновесной температуры, причем с уменьшающейся скоростью, так как рост концентрации альбита замедляется. При дальнейшем остывании расплава наступает момент, когда температура равновесия становится больше температуры расплава и снова начинается рост кристалла, соответствующего новому составу расплава в приграничном слое. Период роста кристалла значительно короче периода отсутствия роста, так как процесс диффузии протекает медленнее процесса кристаллизации. Общее содержание анортита в расплаве убывает быстрее, чем общее содержание альбита, так как кристалл при данных параметрах растет преимущественно из анортита.
На рисунке 2.8 приведены характерные профили концентраций анортита и альбита в неподвижной системе координат для одного из периодов роста. Линии проведены через равные промежутки времени, начало каждой линии соответствует текущему положению границы кристалла. По профилям хорошо видно быстрое убывание концентраций компонент в приграничном слое во время роста кристалла (расстояние между линиями относительно велико) и затем более длительное перераспределение концентраций в расплаве при отсутствии роста кристалла. Причем перераспределение анортита происходит значительно быстрее, по сравнению с альбитом, за счет более быстрой диффузии.
В случае с учетом растворения также существует диапазон значений параметров, при которых происходит периодическая смена режимов роста и растворения. На рисунке 2.9 представлен общий характер протекания процесса кристаллизации с возможным растворением кристалла.
На рисунках 2.10 и 2.11 показано характерное поведение системы кристаллрасплав для одного из периодов роста (первого).
Рис.2.10: профили концентраций анортита (а) и альбита (б) в расплаве На рис. 2.11-а,б показано изменение в процессе кристаллизации основных величин на границе кристалл-расплав (концентраций компонент в расплаве и кристалле, положения границы). Этот расчет проводился при тех же значениях параметров, что и в случае без учета растворения.
Рис.2.11: а – концентрация анортита в расплаве на границе от положения границы кристалла, б – концентрация анортита в кристалле (показано начало второго периода); в – изменение размера кристалла, концентрации анортита в кристалле и концентраций в расплаве анортита (1) и альбита (2) на границе По графикам видно, что в процессе кристаллизации на одном периоде роста условно можно выделить несколько этапов. Сначала преобладает процесс роста кристалла, и концентрации анортита и альбита на границе убывают. При достижении достаточно больших градиентов концентраций в пограничном слое начинает преобладать процесс диффузии и концентрации анортита и альбита в расплаве на границе начинают расти. При этом концентрация альбита растет быстрее и начинает расти раньше, поскольку за счет его более медленной диффузии градиент концентрации альбита значительно больше (для данного расчета). За счет этого температура равновесия убывает с увеличивающейся скоростью. Т.к. переохлаждение было отрицательным, то кристалл еще какое-то время продолжает расти. Но, поскольку скорость остывания меньше скорости убывания равновесной температуры, наступает момент, когда переохлаждение становится положительным, рост кристалла прекращается и начинается его растворение.
В процессе растворения граница кристалла сдвигается обратно внутрь кристалла, и в расплав идет поток компонент из растворившихся слоев кристалла.
При этом продолжается рост концентраций компонент в расплаве на границе кристалла. Концентрация альбита растет быстрее за счет того, что перепад его концентрации в ячейке еще достаточно большой, тогда как анортит, за счет более быстрой диффузии, достаточно быстро уравновешивается. Постепенно процесс диффузии замедляется, и в какой-то момент диффузия компонент в пограничном слое начинает идти в обратную сторону за счет поступления вещества из растворившихся слоев в пограничный слой, и концентрации компонент в расплаве на границе уменьшаются. За счет медленного изменения концентраций, скорость изменения температуры равновесия мала, и при постоянной скорости остывания наступает момент, когда переохлаждение становится отрицательным и снова начинается рост кристалла уже с другим составом, соответствующим новому составу расплава в пограничном слое (переход к следующему периоду).
Для сравнения двух случаев: без учета и с учетом растворения – были проведены расчеты при одних и тех же начальных условиях и одинаковых значениях определяющих параметров. Сравнительные графики зависимостей концентрации анортита и альбита в расплаве и в кристалле на границе от положения границы и зависимости размера кристалла от времени приведены на рисунке 2.12.
Рис.2.12: а – зависимость от положения границы концентрации анортита в кристалле; б – зависимость размера кристалла (положения границы) от времени. Сплошная линия – учет растворения, пунктирная линия – отсутствие растворения Эти графики отражают уже описанные выше процессы, происходящие во время роста и растворения кристалла. График на рисунке 2.11-а также показывает распределение концентраций в образующемся кристалле в обоих случаях. Для случая с учетом растворения на графике присутствуют возвратные линии:
участки, соответствующие смещению границы внутрь кристалла при его растворении. Если их не учитывать, то получим окончательное распределение концентрации анортита в кристалле. В реальных условиях и при проведении экспериментов наблюдается растворение кристалла, когда температура выше температуры равновесия. Если растворение кристалла не учитывается при моделировании, то получено, что периоды отсутствия роста достаточно продолжительны, поэтому предположение об отсутствии растворения приводит к значительному искажению результатов и противоречит наблюдениям.
На рисунках 2.13-15 показано влияние различных параметров на процесс кристаллизации (для случая, когда скорость роста определяется в зависимости от переохлаждения расплава, растворение кристалла учитывается). Результаты расчетов сравниваются с аналитическим решением для равновесного роста кристалла (п. 2.4). Расчеты приближаются к равновесному решению при уменьшении VT ' и увеличении параметров C D 3, CD (рис. 2.13, 2.15). При этом, в отличие от предыдущего случая (п. 2.7), существенное влияние на процесс кристаллизации оказывает параметр C D 3, при неизменном максимуме скорости роста отвечающий за величину числа Пекле (в этом случае число Пекле можно диффузионное число Пекле убывает. Для CD3 105 число Пекле равно Ped 0.1, и решение близко к равновесному, приближаясь к нему с убыванием VT ' (рис.
2.13-б) и ростом CD (рис. 2.15-б).
В случае сильно неравновесного решения при некоторых значениях C D (например, CD3 104 ) и CD возможен рост кристалла со ступенчатым профилем (сложнозонального). Кристалл растет немонотонно при значениях безразмерной скорости остывания VT ', меньших некоторого максимального значения. При уменьшении VT ' количество полос в кристалле увеличивается, а профиль анортита, в среднем, приближается к равновесному решению (рис. 2.13-а). При отрицательных VT ' кристалл также монотонно растет, в отличие от результатов, полученных для задачи Стефана. Влияние параметра Tu ', отвечающего за величину углового коэффициента kT, менее существенно: при его росте незначительно увеличивается число полос в кристалле с сохранением профиля анортита в среднем (рис. 2.14). Параметр не влияет на процесс кристаллизации при C D 1 (рис. 2.15), когда C D 1 решение приближается к равновесному с ростом этого параметра. При значениях CD, близких к 1, кристалл растет монотонно; при меньших значениях CD возможен рост кристалла со сложной периодической зональностью.
Рис.2.13: влияние VT ' и C D 3.
Получено, что при некоторых значениях коэффициентов диффузии компонент (параметров C D 3 и CD соответственно) в плоскости ( Tu ',VT ') можно выделить область значений этих параметров, в которой наблюдается периодический рост кристалла (немонотонное изменение концентрации анортита в кристалле) – рис. 2.16. Для фиксированных коэффициентов диффузии только эти два безразмерных параметра определяют характер протекания процесса кристаллизации (п. 2.3). По графикам видно, что между областью, в которой растет сложнозональный кристалл, и областью, где кристалл растет монотонно, можно провести четкую границу. Хотя эти области получены лишь для некоторых значений коэффициентов диффузии, можно увидеть некоторые закономерности явления. Существует максимальное значение безразмерной скорости остывания VT ' (при фиксированных коэффициентах диффузии), выше которого кристалл всегда растет монотонно. При увеличении C D 3, т.е. при приближении к равновесному процессу, минимальное значение VT ' возрастает. Параметр CD существенно не влияет на величину минимального VT ' (в случае, если существует область периодического роста), хотя при малых значениях CD с уменьшением VT ' возможно полное растворение кристалла (соответствующая область отмечена треугольниками на графиках). В этом случае кристалл полностью растворяется после первого периода роста. И хотя, в дальнейшем, по результатам расчетов, он растет со сменой режимов роста и растворения, предполагается, что такая ситуация не описывается предложенной моделью, поскольку в ней не учитывается нуклеация кристаллов.
2. Построенная в главе 1 модель диффузионного роста кристалла применена для качественного исследования кристаллизации плагиоклаза в остывающем магматическом расплаве. Использованы упрощенные зависимости температуры равновесия и состава кристалла от состава расплава. Предполагалось, что коэффициенты диффузии компонент постоянны, а температура расплава убывает линейно.
Были рассмотрены различные способы определения скорости роста: из условия равенства нулю локального переохлаждения на границе кристалл-расплав (условие переохлаждения расплава на границе с кристаллом. Показана необходимость учета растворения кристалла при втором способе задания граничных условий.
В обоих случаях исследовался процесс кристаллизации при различных значениях определяющих параметров. Показано, что граничные условия в форме Стефана могут использоваться при моделировании процессов, близких к равновесным. Однако они не позволяют объяснить возникновение зональности (периодического изменения содержания анортита в кристалле). В случае, когда на границе допускается ненулевое переохлаждение, получены области значений параметров, при которых происходит смена режимов роста и растворения кристалла (отсутствия роста) и образуются сложнозональные кристаллы.
Глава Влияние реальных зависимостей состава кристалла и температуры равновесия на кристаллизацию плагиоклаза Основные предположения и постановка задачи.
3. В этой главе, как и в предыдущей, моделируется рост кристалла плагиоклаза в остывающем магматическом расплаве. Предполагается, что давление постоянно, а температура расплава линейно убывает со временем, скорость ее убывания – один из свободных параметров задачи. Вместо модельных зависимостей температуры равновесия и состава кристалла от состава расплава используются реальные зависимости для плагиоклазов в водонасыщенных магматических расплавах.
Состав кристалла, температура равновесия и скорость роста кристалла Как показано в предыдущей главе, кинетику кристаллизации определяет, в основном, переохлаждение расплава – разность температуры расплава и равновесной температуры. Температура равновесия Тeq зависит от состава расплава (определяется величиной Xan ), эта зависимость получена экспериментально [30] (рис. 3.1-а). Состав кристалла плагиоклаза помимо Xan зависит также и от температуры расплава [43] (рис. 3.1-б). При высоких температурах зависимость концентрации анортита в кристалле от Xan хорошо описывается формулой (2.2.7).
Рис.3.1: зависимость температуры ликвидуса (а) и состава кристалла (б) от состава расплава на границе кристалла.
В этой главе скорость роста кристалла определяется по формуле (2.2.6) с заданным значением Tu 62 K. Рассмотрены оба случая: с учетом и без учета растворения кристалла (п. 2.2).
Определение коэффициентов диффузии экспоненциально убывают с уменьшением температуры (рис. 3.2), а коэффициент коэффициентами диффузии анортита и альбита.
Далее будет исследовано влияние различных параметров на процесс кристаллизации. В расчетах задан начальный состав расплава – общее содержание анортита и альбита в расплаве – как C10 C20 0.2 (т.е. Xan0 0.5 ). Температура равновесия при таком составе расплава равна 1322.4 K. Температура расплава убывает линейно с начального значения, равного температуре равновесия до значения в 1200 K со скоростью VT.
В расчетах были приняты следующие значения величин: V0 108 м / с – максимум скорости роста, L 104 м – размер ячейки. Характерное значение коэффициента диффузии D0 0.1 10 1014 м2 с (предполагается фиксированное отношение коэффициентов диффузии анортита и альбита и D0 ) и скорость остывания VT 0.1 100 K час варьировались в широких пределах.
В случае без учета растворения начальное положение границы кристалла находилось в нуле: s0 0. Когда учитывалось растворение кристалла, во всех расчетах принято начальное положение границы в середине ячейки: s0 0.5.
Результаты в случае отсутствия растворения кристалла.
3. В широкой области значений параметров происходит периодическая смена режимов роста и отсутствия роста. Общий характер протекания процесса диффузионного роста кристалла описан в предыдущей главе.
При увеличении скорости остывания количество периодов роста кристалла уменьшается (притом, что размер образовавшегося в итоге кристалла изменяется мало), а перепад концентраций компонент в кристалле между двумя периодами роста становится больше (рис. 3.3,а).
Рис.3.3: зависимость концентрации анортита в кристалле C1S от расстояния s от центра кристалла (а) и размера кристалла от времени Это связано с тем, что переохлаждение увеличивается по абсолютному значению и дольше остается отрицательным (рис. 3.4-а). В результате, скорость достигает бльших значений, а ее величина дольше остается сравнимой с максимумом (рис.
3.4-б). За счет более быстрого и продолжительного роста кристалла, во-первых, в пограничном слое расплав более обеднен, поэтому для перераспределения концентраций в расплаве требуется больше времени (рис. 3.3-б). Во-вторых, большая доля анортита и альбита из расплава переходит в кристалл, и их общее содержание в расплаве уменьшается сильнее, что приводит к большему скачку состава кристалла при переходе к следующему периоду роста кристалла.
Рис.3.4: зависимость переохлаждения (а) и скорости роста кристалла (б) от расстояния от центра кристалла для D0 1014 м2 / с и Xan0 0.5.
При уменьшении начального значения Xan и неизменной скорости остывания количество периодов роста практически не изменяется, но снижаются скачки состава расплава (рис. 3.5). Это связано с тем, что при меньшем Xan равновесная температура меньше, а скорость ее изменения больше (рис. 3.1,а).
Результаты в случае с учетом растворения кристалла 3. Как и в случае без учета растворения, получено, что, как следствие периодического роста, в кристалле образуются полосы разных составов. В пределах одной полосы состав меняется незначительно по сравнению со скачком состава при начале следующего периода роста. В обоих случаях образующиеся в кристалле полосы имеют примерно одинаковую ширину (за исключением первой), а периоды роста имеют одинаковую продолжительность. Можно заметить, что ширина полос (как и длительность одного периода роста) для случая с учетом растворения значительно меньше, чем для случая без учета растворения (примерно в 2 раза на основном участке – рис. 3.6-а) Рис.3.6: а – профиль анортита в кристалле, б – зависимость размера кристалла (положения границы) от времени для D0 1014 м2 / с и VT 5 K / ч. Сплошная линия – учет растворения, пунктирная линия – отсутствие растворения.
Это объясняется тем, что растворение ускоряет процесс диффузионного перераспределения расплава. На графике (рис. 3.7) точками показаны результаты вычисления средней ширины колец для разных скоростей остывания (график построен в безразмерных величинах). Здесь Tm – характерный диапазон изменения температур (разница между начальной и конечной температурами расплава), все расчеты проводились для Tm 120 K. Зависимость средней ширины полос от скорости остывания хорошо аппроксимируется полиномом 4-й степени (сплошная линия).
При уменьшении скорости остывания до очень небольших значений ширина полос уменьшается, хотя перепады состава между полосами настолько малы, что становятся практически неразличимы, а продолжительность периодов роста изменяется мало.
При увеличении скорости остывания средняя ширина полос практически не зависит от скорости остывания. Когда скорость остывания достигает некоторого критического значения, происходит качественное изменение характера протекания процесса кристаллизации: кристалл растет монотонно, отсутствуют смена режимов роста и растворения и, соответственно, полосы разных составов в кристалле. Были проведены расчеты для разных значений безразмерных параметров (определяемых характерным коэффициентом диффузии и скоростью остывания), которые являются независимыми параметрами, определяющими процесс кристаллизации. Получена (рис. 3.8) область значений параметров, в которой кристалл растет со сменой режимов роста и растворения (образуются зональные кристаллы). При переходе границы области происходит качественное изменение характера протекания процесса, особенно существенно проявляющееся при малых D0: кристалл начинает расти монотонно. При быстрой диффузии (характерное значение D0 велико) процесс кристаллизации близок к равновесному и кристалл растет монотонно. При существенной неравновесности процесса (значения D0 малы, на графике соответствует области с образуются кристаллы с выраженной зональностью: количество полос разного состава в кристалле велико. Зависимость средней ширины полос от скорости остывания аналогична показанной на рис. 3.7. При приближении значения безразмерного D0 к верхней границе существования зональности (около 0.065) количество полос в кристалле уменьшается и при значениях, близких к верхней границе, происходит единственная смена режима роста на режим растворения кристалла (и обратно). Профиль анортита в кристалле при этом мало отличается от профиля при монотонном росте кристалла. В этом диапазоне значений D0 (от 0.025 до 0.065) при малых скоростях остывания происходят периодические колебания скорости роста кристалла, но при этом она остается положительной и кристалл растет монотонно. Помимо верхней границы значений безразмерной скорости диффузии ( D0 V0 L ), существует также максимальная безразмерная скорость остывания, выше которой не могут образовываться зональные кристаллы. Этому значению безразмерного параметра, равному примерно 3.1, соответствует скорость остывания 130 К/ч (при выбранных значениях параметров задачи). Существование максимальной скорости остывания, при которой может существовать зональность в кристаллах, было получено и при качественном исследовании процесса кристаллизации, проведенном в главе 2.
3. Исследован процесс роста кристалла плагиоклаза в остывающем магматическом расплаве для реальных зависимостей температуры равновесия и состава кристалла от состава и температуры расплава. При этом зависимость концентрации анортита в кристалле от анортитотового номера (величины Xan) расплава на границе с кристаллом при высоких температурах хорошо описывается формулой (2.2.6), используемой в главе 2 (рис. 3.1,б). При более низких температурах эта зависимость сохраняется на качественном уровне.
Зависимость температуры равновесия от состава расплава – нелинейная функция Xan. Для значений Xan, достигаемых в процессе кристаллизации угловой коэффициент kT изменяется от 200 до 1000 K (рис. 3.1,а) В обоих случаях (без учета растворения и в случае с растворением кристалла) в широком диапазоне параметров при линейном падении температуры получен ступенчатый рост кристалла: режимы роста сменяются режимами растворения (отсутствия роста в первом случае). В результате в кристалле образуются полосы разного состава (зональность). Этот результат соотносится с результатами, полученными в главе 2 для таких же значений параметров.
Для случая с учетом растворения кристалла получена зависимость средней ширины полос в кристалле от скорости остывания. При очень медленном остывании полосы также образуются, но они практически неразличимы. При быстром остывании расплава зависимость ширины полос от скорости остывания слабая.
Проведено исследование зависимости характера протекания процесса кристаллизации (смена режимов роста и растворения или непрерывный рост кристалла) от скорости остывания и скорости диффузии компонент (характерного значения коэффициента диффузии). Получен критерий появления зональности:
когда скорость остывания достигает определенного значения, полосы в кристалле исчезают, и кристалл непрерывно растет. При скорости остывания ниже критической образуются зональные кристаллы.
Глава Рост кристалла плагиоклаза при подъме магмы по каналу вулкана Постановка задачи и основные предположения.
4. В процессе роста состав кристалла существенно зависит от условий, в которых происходит кристаллизация. В кристаллах плагиоклаза выделяют 2 зоны:
ядро (основная часть кристалла, вырастающая в очаге вулкана) и кайма, которая нарастает в процессе подъема по каналу вулкана. Рост кристалла, вызванный падением давления (при подъеме магмы по каналу) – широко распространенное явление в водонасыщенных магмах [4,5,7,19]. Экспериментальные исследования [7,34,39] показывают, что с ростом скорости падения давления растет содержание анортита в кристалле плагиоклаза. Это позволяет оценить скорость падения давления по результатам измерения состава кристалла. В этой главе по изменению концентрации анортита в каймах кристаллов плагиоклаза реконструируются условия подъема магмы при извержениях вулкана Безымянный (Камчатка).
Вулкан Безымянный поставляет материал на поверхность примерно каждые полгода. Зональность плагиоклаза была исследована для извержений из диапазона 2000-2007 гг. [45]. Были выбраны извержения 2000, 2006, 2007 годов как наиболее типичные для последнего цикла активности (рис. 4.1). Различия между профилями может быть вызвано различными условиями после извержения.
Процесс кристаллизации моделируется в системе предположений, описанной в п. 1.1. При этом, плоскость x 0 (рис. 1.1) соответствует плоскости раздела между каймой и ядром кристалла и ось направлена от этой плоскости к границе кристалла. Проводится сравнение результатов расчетов с профилями анортита в каймах выбранных кристаллов плагиоклаза. Размер ячейки определяется средним расстоянием между фенокристами и в расчетах берется равным 0.35 мм.
Исследуется рост кристалла плагиоклаза при подъеме магмы по каналу вулкана. Предполагается, что в процессе подъема давление падает линейно, но скорость падения давления может меняться в процессе подъема. Давление изменяется от 100 МПа в начале подъема до атмосферного (0.1 МПа). Далее будет исследоваться влияние условий подъема (изменение скорости падения давление, возможное остывание расплава) на процесс кристаллизации. Начальное давление было оценено по содержанию воды в расплавных включениях [45].
Состав кристалла и температура равновесия.
Температура равновесия плагиоклаза описывается функцией состава расплава, давления и содержания воды. Состав плагиоклаза в равновесии с расплавом определяется теми же параметрами и равновесной температурой.
Моделирование неравновесной кристаллизации требует дополнительного предположения, что в процессе кристаллизации существует локальное равновесие по составу на границе кристалл-расплав. Другими словами, температура и состав плагиоклаза могут быть определены как равновесные с расплавом, находящимся вблизи с границей растущего кристалла.
Общее содержание компонент в системе оценивается из стекловатых расплавных включений в фенокристах ортопироксена и плагиоклаза с извержения вулкана Безымянный в 2006 г. [45]. Рассчитывался средний состав расплавных включений и разлагался на 3 компоненты: альбит, анортит и некий гипотетический остаточный расплав (в котором не содержатся компоненты плагиоклаза). Затем была вычислена 5% трехмерная сетка с различными пропорциями альбита, анортита и остаточного расплава. Для каждой точки этой сетки было оценено содержание воды по модели [36] в диапазоне давлений от 0. до 200 МПа. Температуры кристаллизации и составы плагиоклаза определялись для водонасыщенных условий по модели равновесия плагиоклаз-расплав [43] для каждой точки трехмерной сетки гипотетических расплавов. Эта модель плагиоклаз-расплав была выбрана потому, что среди всех известных моделей только модель [43] учитывает влияние активности воды в расплаве на состав плагиоклаза, и это влияние может быть основной причиной зональности плагиоклаза в островодужных магмах [21,28,44]. Все вычисления были произведены в пакете программ PETROLOG-3 [13].
В результате расчетов, была получена дискретная функция, ставящая в соответствие каждому набору значений концентраций анортита, альбита и давления соответствующие температуру равновесия и состав кристалла (содержание анортита в кристалле). Для определения равновесной температуры и состава кристалла для произвольных состава расплава и давления были получены гладкие функции 3-х переменных с помощью трикубической интерполяции сплайнами соответствующих дискретных функций.
На рис. 4.2. приведены графики зависимости состава кристалла и температуры равновесия от содержания анортита в расплаве при разных давлениях и фиксированном содержании альбита. Давление меняется в диапазоне от атмосферного до 100 МПа.
Рис.4.2: графики зависимости состава кристалла и температуры На рис. 4.3 показано сравнение зависимости концентрации анортита в кристалле от Xan C1 (C1 C2 ) при давлении 100 MPa с квадратичной зависимостью, используемой в главе 2. Как видно, квадратичная функция хорошо описывает зависимость концентрации анортита в кристалле от Xan при концентрациях M остаточного расплава близких к 0.5.
Определение коэффициентов диффузии.
Коэффициенты диффузии зависят от температуры и содержания воды в магме. Эта зависимость выражается уравнением Аррениуса D D0e. При повышении температуры за счет выделения скрытой теплоты кристаллизации скорость диффузии растет. При понижении давления водонасыщенность магмы уменьшается, что ведет к уменьшению скорости диффузии. Как будет показано ниже, изменение температуры при подъеме магмы мало – порядка 20 градусов, поэтому в первом приближении можно считать коэффициенты диффузии анортита и альбита постоянными. Далее анализируется влияние отношения CD D1 D2 коэффициентов диффузии компонент на характер протекания процесса кристаллизации. За характерное значение коэффициентов диффузии было принято значение коэффициента диффузии анортита D1, которое также варьируется для определения влияния скорости диффузии в целом на процесс кристаллизации. Коэффициент диффузии остаточного расплава, как и ранее, мал по сравнению с коэффициентами диффузии анортита и альбита.
переохлаждения и вычисляется по формуле (3.1.1).
Результаты расчетов по равновесной модели.
4. (бесконечная скорость диффузии). При применении модели, описанной в п.1.2., для задачи о кристаллизации плагиоклаза получена следующая система уравнений:
В силу уравнений (4.2.6) и (4.2.7) достаточно использовать только одно условие для определения состава кристалла (например, для анортита) и можно считать, что функции для определения состава кристалла и температуры равновесия зависят от 3 переменных: давления и концентраций анортита и альбита в расплаве. Уравнение (1.1.6) выполнено автоматически. Эти уравнения можно переписать в виде (см. (1.2.10)):
В результате, имеется 2 уравнения на 3 неизвестные величины ( C1S, s, p).
Для замыкания системы будет задан профиль анортита в кристалле. Для этого построим отобранные в п. 4.1 профили на одном графике (рис. 4.5), начиная с максимума концентрации анортита в кайме ( s 0 соответствует максимуму концентрации, направление оси – от центра к периферии кристалла). На начальном этапе роста кайм профили анортита практически совпадают (до s 0.04 мм ), а затем начинаются значительные расхождения.
«