Анализ режимов ритмической кристаллизации при направленном затвердевании магматических расплавов Текст научной статьи по специальности «Физика»

УДК 536, 551.2

Анализ режимов ритмической кристаллизации при направленном затвердевании магматических расплавов

В.К. Черепанова

Новосибирский государственный технический университет, Новосибирск, 630092, Россия

Проведен подробный качественный и количественный анализ условий развития колебательные режимов кристаллизации при направленном затвердевании квазибинарного расплава. Показана возможность установления автоколебаний в природной системе Di-An и оценена область изменения тепло- и массообменных параметров, при которык может реализовываться ритмический характер затвердевания.

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

Rhythmic crystallization modes in magmatic melts during directional solidification

V.K. Cherepanova Novosibirsk State Technical University, Novosibirsk, 630092, Russia

The paper presents the results of thorough qualitative and quantitative analysis of oscillatory crystallization modes of a quasi-binary melt during directional solidification. The possibility of self-oscillations in a natural Di-An system is demonstrated and the order of variation of heat and mass transfer parameters at which rhythmic solidification can be realized is estimated.

Keywords: crystallization, heat and mass transfer, mathematical modeling, magmatic melt

1. Введение

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

представления о механизмах, вызывающих появление в них слоистой структуры, пока нет [8-14].

Авторы многих работ [5, 10, 12 и др.] связывают образование слоистых структур с ритмическим характером кристаллизации, которая в техногенных условиях затвердевания сплавов обуславливается высоким положительным градиентом температуры, конвекцией жидкой фазы и достаточно узким температурным интервалом кристаллизации. Последовательный обзор результатов эмпирических наблюдений над процессом возникновения слоистых неоднородностей в слитках был сделан в [6]. Было показано, что важная роль при развитии такого рода явлений принадлежит динамике переохлаждения, являющегося движущей силой роста кристаллов. С помощью математической модели неравновесного процесса в работе [10] было получено замкнутое решение для задачи тепломассопереноса в активной двухфазной области, которое позволило показать условия су-

© Черепанова В.К., 2009

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

2. Постановка задачи и ее решение

Рассмотрим направленную кристаллизацию псевдо-бинарного магматического расплава, который заполняет полуограниченную интрузивную камеру: -х0 < х < х0, 0 < г < ^ -то < у < го. Считаем, что фронт двухфазной зоны движется вдоль оси г, а ее начало отсчета находится на поверхности, отделяющей вмещающую породу от расплава (рис. 1, а). Боковые поверхности интрузива, контактирующие с вмещающей породой, средняя температура которой Тс ниже температуры расплава, обмениваются с ней теплом (вдоль оси х) по закону Ньюто-на-Рихмана с коэффициентом теплообмена ас. Примем некоторые допущения, не меняющие существа рассматриваемого процесса. Будем считать, что 1) поперечный размер 2х0 интрузивного тела достаточно мал, так что его внутренним термическим сопротивлением можно пренебречь по сравнению с внешним (х0/X << 1/ас, где X — теплопроводность магмы) и рассматривать задачу в квазиодномерном приближении; 2) физические параметры среды равны соответствующим средним значениям в рассматриваемых интервалах изменения величин; 3) движение расплава обусловлено лишь усадочными явлениями при кристаллизации, что позволяет им пренебречь.

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

эт __д.хдТ дt дz I дz

CPliq ^=^| ^ I- — (T - Tc ) + AH 3/so1

a c

dt ’

(1)

d( /liq C ) dt

d/sol =

I(D“q/‘q I

-£jkC

d/s

sol

dt

(2)

-F(AT), /sol - 1 /liq,

dt d1(vT)

/liq = 1 При T > Tf, /iq = 0 При T < 7^

(3)

(4)

где T — температура; С—массовая концентрация второго компонента в расплаве, %; f— доля фазы; индексами sol, liq обозначены величины, относящиеся к твердой и жидкой фазам соответственно; с — теплоемкость; Dliq — коэффициент диффузии в жидкой фазе (диффузией в твердой фазе пренебрегаем, поскольку Dsol << Dliq); р — плотность; k— равновесный коэффициент распределения второго компонента в первом; AH = psol к0; к0 — теплота плавления рассматриваемой псевдобинарной системы; е1 - 1 + е; е- (psol - pliq)X х1/psol; d1 — расстояние между первичными осями ден-дритов; vT — скорость охлаждения; F(AT) — феноменологический закон, характеризующий зависимость линейной скорости роста кристаллов от среднего переохлаждения AT расплава в зоне; индексами f, e обозначены значения величин на фронте и в конце гете-рофазной зоны соответственно. При осреднении уравнения массопереноса принималось, что массообмен на границе с вмещающей породой отсутствует. Решение системы уравнений (1)- (4) должно удовлетворять следующим начальным и граничным условиям:

Tt-0 = T~ ’ Ct -0 = C0> /liq It-0 =1;

TL^T~> Cz^ C0> TL=0 = Tc;

(5)

(6)

dz

-A

liq

Ze -0

dC

dz

-xdT

dz

+ psolK^/liqe ze(t)?

e+0

-(1 - k) ze( t);

e+0

-xdT

dz

-Du

= af (Tf Tv°),

1600

о

° 1400

1200

дс liq я dz z=zf = M-0 (Cf C0 ),

. Di + ж"° ~- An +Ж

Di + An

(7)

(8) (9)

Ш

CaMgSi206 20 Диопсид

40 60

Вес. %

80 CaAI2Si208 Анортит

Рис. 1. Схема физической модели двухфазной зоны (а): I — гомогенная жидкость, II — активный район двухфазной зоны, III — область роста твердой фазы в гетерогенной зоне, IV—затвердевшая изверженная порода (Di +Ап), 1 — кристаллическая фаза диопсида Di, 2 — кристаллическая фаза анортита Ап; рассматриваемая диаграмма состояния псевдобинарного расплава Di +Ап (б). Ж — расплав

где 7^ — средняя температура жидкого ядра интрузива; С0 — исходная концентрация второго компонента; af — коэффициент теплоотдачи на фронте; Ц0 — коэффициент массообмена; точка над функцией означает производную по времени.

Система уравнений (1)-(4) с начальными и граничными условиями (5) - (9) представляет собой общую постановку задачи о направленной кристаллизации бинарного сплава [10], решение которой требует применения численных методов. В данной работе представлено приближенное решение рассматриваемой задачи, основанное на некоторых упрощающих допущениях, что, тем не менее, позволяет выявить основные условия возникновения колебательного режима кристаллизации.

Примем, что течением расплава, обусловленным усадочными явлениями в гетерогенной зоне, можно пренебречь. Размер дендритной ячейки dl (рис. 1, а) слабо зависит от изменения скорости охлаждения vT , так что dl = const. Считая далее, что кристаллы имеют форму пластин, а их рост в зоне подчиняется нормальному механизму, запишем:

F(AT) = KAT, (10)

где K — полуэмпирическая константа роста твердой фазы. Диаграмму состояния (рис. 1, б) рассматриваемого псевдобинарного магматического расплава [15] можно аппроксимировать линейной зависимостью. Тогда для величины AT имеем:

AT -ГА-Р0С-T, (11)

где TA — температура плавления чистого компонента; в0 — модуль коэффициента наклона линии ликвидуса. С учетом принятых допущений и соотношений (10), (11) уравнение (3) примет вид:

д/soi = 2K d1

so1 -^ (Ta-Р0С - T).

dt

(12)

Распределение концентрации С определим из решения уравнения массопереноса (2) в квазистацио-нарном приближении [10], рассматривая двухфазную зону как квазиравновесную систему (АТ =0, Т = _ ТА - в0 С) с учетом граничного условия (9). Допуская, что градиент температуры G слабо меняется вдоль зоны, определим его как некоторую среднюю величину <0) _ _ (Gf + 0^/2, где — градиент температуры на фронте зоны; 0е — градиент температуры в конце двухфазной зоны. Тогда Dщ <0)

C-

Р0(1- k) zf

с+

DliqGf

1 --

^0

(1- k) zf

•(1 - /sol)

k -1

(13)

Полученное уравнение при Dliq _ 0 совпадает с известным уравнением неравновесного рычага [16] и является его обобщением на случай учета процессов конвек-

тивного тепло- и массообмена на переднем фронте двухфазной зоны. Принимая во внимание, что для рассматриваемой системы коэффициент распределения к << 1, полагаем далее величину k = 0.

При направленном затвердевании расплава с образованием двухфазной зоны вначале в ней выделяются кристаллы первого компонента, а в конце зоны при достижении расплавом эвтектического состава Се происходит одновременное выделение второй фазы. Безразмерное значение среднего состава затвердевшей породы определим с учетом (13) при k = 0 как

(С> = Ceii(1 -ц0/f

(14)

Се *1^0/ Zf

где Се _ Се /С0; Ь _ ^^/^0^0); *1 _ D1iq <0) х

Х(р0С0^0) .

Выделим внутри двухфазной зоны некоторую узкую область zk < z < zf (рис. 1, а), прилегающую к фронту и характеризующуюся наибольшим для рассматриваемого процесса значением переохлаждения АТ (активный район двухфазной зоны [6]). Ширина этой области может составлять порядка 0.1 размера двухфазной зоны [6, 10]. В данной области протекают основные процессы образования структуры затвердевающего расплава. Тогда как в остальной (квазиравновесной) части зоны, где переохлаждение пренебрежимо мало, осуществляется лишь рост возникших структурных составляющих. Выделенную часть гетерофазной зоны будем рассматривать как динамическую систему, на границах которой происходит тепломассообмен с окружающей средой. При этом считаем, что на фронте z _ zf имеет место конвективный теплообмен с расплавом по закону Нью-тона-Рихмана (8). В конце активного района двухфазной зоны z _ zk _ Zf -10 (10 — его характерный размер) задан удельный поток тепла дк ^) в квазиравновес-ную часть зоны:

^ . (15)

&

dz

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

Zf _ К0(ТА -р0Cf -Tf), (16)

где К0 — кинетический коэффициент роста кристаллов вдоль направления движения фронта.

Переходя далее в уравнениях (1), (8), (12), (13), (15), (16) к безразмерным величинам

0_ (Т -Тс)/АТк , С _ С/С,, zf _ if /V*,

7_ ^1т , т_^и, к_К0/(сАТк), р_р0С^АТк, Ц_Ц^[(1 - к)%], получим

Э0_ 3^0 дз

дт Э72 дт

sol

- Stc6,

(17)

д/

БОІ

= Л[Єііч0 -0-р(С-1)],

Эт

С = ЬР/+ [1 + Ь(1 - р/%)] (1 - /о )

к-1

Э0

дї

_ д0

^ С f = ЗД„-0Д —

2 = Zf

% =01ідО -0f -вЬ

дї

= С к,

(18)

(19)

(20)

(21)

где 0Цч() _ (ТИЧ0 - Тс)/АТк ; А _ 2аК/(К02АТкdx);

0к _ qk 1т/(АТкХ); 1Т _ а/^

^ _ 1т Iа; Йс _ас/(^ичс^х0); _%/(^цчcv*);

х0 _ х0/ 1Т; v* _ К0АТк; АТк _ ТЦц0 - Те; ТНц0 температура ликвидуса; а — температуропроводность. Учитывая, что /,о1 < 1, выражение для С в (19) можно записать в линеаризованном виде:

С _ 1 + Ь + (*1 - Ь)р/Zf +

+ [1 + Ь(1 -р/7£)](1 - к)/оХ. (22)

Из (18) с учетом (22) следует:

дт

= А(0ііЧо -0-вЬ2 -В/яо1),

(23)

В = Р(1 -к )[1 + Ь(1 -р/ tf)],

где Ь2 = Ь + (Ь1 - Ь)р/7Г.

Систему уравнений (17), (23) можно свести к одному уравнению относительно 0, исключив величину /,о1. В результате оно будет иметь следующий вид:

д 2 0 Эт 2

Хг

Э 2 0

д30 Э0

-----з + Х2 — + Хз0 = Ф,

(24)

д!2 дтд!2 ^ дт

где Х1 _ АВ - В/В; Х2 _ Aк + Stc + Х Хз _ StcХl -- АкВ/В; ф_-АВкЬ2 - АкВ (0Кч0-РЬ2)/В.

Из первого граничного условия (20) с учетом выражения для 0f, следующего из (21), имеем:

Оf _ ^и!(1 -у), (25)

где и _0го - 0ич0 + !(■; у _ afД^ДХц,,). Безразмерный параметр у можно представить в виде: у _ Bi0/Sh, где В^ _ ^ 10/X — критерий Био; БИ _ц010/D1iq — критерий Шервуда. Следовательно, параметр у является сложным комплексом, характеризующим соотношения интенсивностей конвективного и молекулярного процессов тепломассопереноса. Выражение для Ь с учетом (25) примет вид:

Ь _у^/в, У1 _у/(1 -у). (26)

Из (21) и (26) найдем соотношение, связывающее температуру среды на переднем фронте двухфазной зоны со скоростью движения фронта:

0f _0го-Ч(1 -у). (27)

Для решения уравнения (24) с граничными условиями (20), (27) воспользуемся продольной разностной

схемой метода прямых [17]. Для этого исходную область 7к < ! < заменим сеточной областью, проведя прямые !■ _ 7к + 1И (г = 1, 2, ..., N - 1), где N _ 10/(МТ) — число отрезков на оси г; h — безразмерный шаг разбиения. На каждой из внутренних прямых ! _ производные по г аппроксимируем значениями функции 0 на нескольких соседних линиях. Для простоты и наглядности исследования ограничимся разбиением области !к < ! < на две ячейки ^ = 2). В результате получим обыкновенное дифференциальное уравнение второго порядка:

01 + (х2 + 2/ ^2 )01 + (хз + 2% / ^2 )01 _

= к 2 (02 + 0О ) + к 2х1 (02 + 0О ) -

- Лк(0ичо -рЙ2)В/В - Ар^,

(28)

(29)

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

0к _01 - О к И, 01 _00 - (1 + В^/(1 -у) и,

Bi _ а И1Т /X.

Уравнения (28), (29) определяют распределение температуры в области !к < ! < '!,1, а также скорость движения фронта кристаллизации. Для единственности решения системы (28), (29) необходимо также задать значения 0г- и Щ при т = 0, то есть начальное распределение температуры и начальное положение фронта.

Подставив в (28) выражения для Ь, В, 02 = 0f, 00 = 0к и 01 из (26), (23), (27), (29), получим дифференциальное уравнение второго порядка относительно переменной и (безразмерной скорости кристаллизации):

и + ^1(и, и)и + gт(и)и _ go + °кgз> (30)

где

Ак

gl _-

<1 + ві + В [(1-у )(р*2 +А0^) - (Ві + 1)и(-

і) I В

(1 + Ві)

(31)

g2 = А(1 -к){ух [1 - р/(и - А0^)](]и --) + Р-0}; gо = А(1 - к)р-;

— = >Віс + Ві/(1 + Ві)]к_2;

-1 = (1 - у) (Ск к + Віс0^) к~21{1 + Ві);

В = (1 - к) {р + ^и [1 -р/ (и-А0^) ]}^

Вх =Ух (1 - к) [1 + рА0со/ (и - А0СО )2 ];

Віс = ^ск2; А0^=0^-0!іЧ0;

gз = (1 -У )/[(1 + Ві) к];

( Р )

+ (Ь - Ьху

Ь21 =

1 --

2(и -А0^)

(и -А0СО)

2

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

Температурный градиент в конце активного района двухфазной зоны можно практически считать постоянным: Gk - const. Поскольку уравнение (30) не содержит явно времени, его можно представить в виде автономной системы

du . du .

■ = u — = g0 - g1(u> u)u - g2(u) u

, , . , „.V , , (32)

dт йт

чьи фазовые траектории удовлетворяют дифференциальному уравнению первого порядка в фазовой плоскости (и, и):

Ом _ g0 - glu - gти йи и

Особые точки будут соответствовать условиям

и _ 0, go - g2(м)и _ 0. (33)

Подставив g0 и g2 из (31) в (33), получим уравнение

(R1 R0u)

P + Y1U

1-

u -A90

- 0,

решения которого выражаются следующим образом:

u1 -

Rl

R0

-P + Y1(A900+p) + ^[P -Y1 (A9^ +Ц)]2 + 4PA900Y1 2Y1

-P + Y1(A9co+P) ^[P-Y1 (A9^ +Ц) ]2 + 4PA9coY1

2y1

В соответствии с теорией устойчивости по Ляпунову [18] для уравнения второго порядка (30) (при условии Gk - const) вида u - /(u,u) и особой точки (u0, u) можно записать характеристическое уравнение

(34)

Точка покоя (особая точка) является:

- устойчивым или неустойчивым узлом, если оба корня ^1 и ^т уравнения (34) действительны и соответственно отрицательны или положительны;

- седлом, если ^1 и ^т действительны и имеют разные знаки;

- устойчивым или неустойчивым фокусом, если ^1 и S2 комплексно сопряжены и имеют соответственно отрицательную или положительную действительную часть;

- центром, если корни чисто мнимы.

В нашем случае при u - 0

du = -^ <1 +Bi + ^[(1 -Y)(P*2 +A9.) -du 1 + Bi [ B

- (Bi + 1)u ] - (1 -Y)pfo21 }-AB - R + Bl( R^u-),

du

ARB + A(1 - k)Y1 Ru - R) X

(

pu

-+1 --

Y

u -A9„

(и -А0го )2

\

Как показал проведенный анализ, автономной системе (32) соответствуют две особые точки типа «седло»: (и0 _ и2; и0 _ 0) и (и0 _ и3; и0 _ 0). Еще одна особая точка (и0 _ и1, и0 _ 0) в зависимости от значений параметров затвердевания и начальных условий может быть центром, устойчивым или неустойчивым фокусом. Соответственно этому в нашей физической системе будут наблюдаться периодические или затухающие колебания, а в случае седла — монотонный характер поведения. Естественно, что особый интерес представляют первые два случая.

Запишем коэффициенты характеристического уравнения (34) в особой точке (и0 _ и1, и0 _ 0), обозначив их 81 и 82:

Ак

1\-dL=-

du 1 + Bi

1 + Bi +

(1 -Y)(Pb +A9.) - (Bi +1)R

R0

-(1 -y)pb211-AB -R0,

82 -^ = - AR0 B. du

(35)

(36)

Тогда условие существования колебательного режима кристаллизации будет определяться соотношением

82 + 482 < 0 при 82 < 0,

(37)

что позволяет установить диапазон варьирования параметров задачи, при которых в нашей физической системе такой периодический процесс будет возможен. Определив значения коэффициентов 81 и 82, можно рассчитать частоту колебаний ю _ 7-§т -8?/4, их период, а также коэффициент затухания 8_-81/2.

3. Обсуждение результатов

Нелинейное дифференциальное уравнение (30) решалось численно методом Адамса. Расчеты были проведены для квазибинарного магматического расплава с составом компонентов Di + 10 % мас. Ап, где диопсид Di является первым компонентом, анортит Ап—вторым

и, 10"'

б

Рис. 2. Фазовые кривые динамической системы при различных параметрах задачи: итч = 1.00045«0, Ок = 265 К/м, /0 = 0.2 м, Sh = 50, Bi = 2.5,

В^ = 10-5, ^0= 5 ■ 10-7м/с, и0 = 8.26 • 10-2, _ 7.46 • 10-7 м/с, af= 52.5 Вт/(м2-К), ас = 2.110-3 Вт/(м2- К) (а); ц^ч = 0.998^, Ок = 130 К/м, /0 =

1и _ и0 .

= 0.15 м, Sh = 70, Bi = 1, Bic = 3.5 • 10-4, ц0 = 9.3 • 10-7м/с, и0 = 8.22 • 10-2, !{\ _ 3.56 -10-7 м/с, af= 28 Вт/(м2-К), ас= 0.13 Вт/(м2-К) (б)

(рис. 1, б). Физические параметры принимались следующими [15]: Т1к10 _ 1665 К, Те _ 1543 К, в0 _ 2.8 К/%, С0 _ 10 % мас., Се _ 42.5 % мас., X = 2.1 Вт/(м-К), с = 1.2-103 Дж/(кг-К), Pliq _ 2.65 • 103 кг/м3, Dlk. _

= 10-9 м2/с, к0 _ 3.2 • 105 Дж/кг, d1 _ 0.005 м, К0 _

_ 10 5 м/(с^К). Расчеты проводились при условиях: АТ^ _ 10 К, Тс _ 423 К, х0 _ 1-2 м; 10 _ 0.1-0.2 м; 0к _ 100-300 К/м, а _ 10-50 Вт/(м2^К), ас _ 0.010.2 Вт/(м2^К), ц0 _ 10-7-10-6 м/с, К= 10-10-10-9 м/с2, БЬ = 50-100, Bi= 1-5, _ 10^-Ю-4. В качестве на-

чальных условий принимались значения и и и, лежащие в окрестности особой точки на фазовой плоскости, характеризующей центр или фокус (и0 _ и1, и0 _ 0).

Результаты численных расчетов представлены на рис. 2 - 5. Фазовая картина поведения динамической системы (рис. 2) показывает, что интересующая нас с точки зрения возникновения колебательных процессов особая точка, соответствующая значениям и0 _ и1, и0 _ 0, в зависимости от начальных условий и значений крис-

таллизационных параметров может быть устойчивым (рис. 2, а) или неустойчивым (рис. 2, б) фокусом. Как видно из уравнения (35), коэффициент затухания 81 определяется целым комплексом параметров. Варьируя эти параметры, например 0к, %, ас, Ц0, К, БЪ, В^ Bic и другие, в допустимых диапазонах (диапазоны варьирования приведены выше), можно получить значение 81 = 0. Значит, особая точка (и0 _ и1, и0 _ 0) может быть и центром, в случае чего в системе возникают гармонические колебания. Кроме того, фазовые портреты системы позволяют предположить существование предельных циклов, что также является условием возникновения гармонических колебаний.

На рис. 3 и 4 представлены временные зависимости соответственно переохлаждения на фронте зоны и скорости продвижения фронта кристаллизации. Как видно из графиков (рис. 3, а, 4, а), при указанных начальных условиях и принятых физических параметрах системы (им соответствует фазовая кривая на рис. 2, а) возмуще-

Рис. 3. Изменение переохлаждения на фронте зоны со временем при параметрах задачи, что и для рис. 2, а (а); при параметрах задачи, что и для рис. 2, б (б)

Рис. 4. Зависимость скорости фронта кристаллизации от времени при тех же параметрах задачи, что на рис. 2

ния на фронте зоны имеют характер слабо затухающих колебаний с периодом порядка 75 ч. В случае иных условий (им соответствует фазовая кривая на рис. 2, б) можно ожидать установления колебаний со слабо нарастающей амплитудой и периодом около 170 ч (рис. 3, б, 4, б). Аналогичная периодичность наблюдается и для среднего состава затвердевшей породы (рис. 5). Надо отметить, что в рассмотренном диапазоне изменения параметров тепломассопереноса характерные значения скорости движения фронта кристаллизации составляют ~10-7м/с. Это соответствует случаю кристаллизации базитовых расплавов в лавовых озерах. Зная период колебательного процесса и скорость движения фронта кристаллизации, можно определять ширину слоев в затвердевших полосчатых структурах.

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

при направленной кристаллизации можно представить следующим образом. При достаточно высоких скоростях кристаллизации (| г | >> DИq/10 ) концентрационное уплотнение перед вершинами растущих кристаллов практически отсутствует («примесь» отводится по бокам кристаллитов). Вершины кристаллов, в соответствии с интенсивностью теплоотвода, продвигаются вглубь перегретого расплава, пока не достигнут точки, где температура станет близкой к локальной температуре ликвидуса. Переохлаждение на фронте двухфазной зоны при этом убывает (рис. 3), чему способствует и диффузионный поток, направленный к фронту зоны и понижающий температуру ликвидуса. Вследствие этого уменьшается температурный градиент на фронте кристаллизации, и скорость перемещения фронта падает (рис. 4). Градиент температуры в конце активного района двухфазной зоны при этом оказывается больше градиента температуры на фронте, т.е. тепловой поток в квазиравновесную часть слитка превышает тепловой поток на фронте. В результате температура расплава на

Рис. 5. Изменение средней безразмерной концентрации Ап в затвердевшей породе в процессе кристаллизации расплава при тех же параметрах задачи, что на рис. 2

фронте зоны начнет снижаться, а концентрационное переохлаждение — увеличиваться. Поэтому скорость продвижения фронта зоны вновь станет возрастать и т.д. (рис. 4). Таким образом, весь цикл периодически будет повторяться. При этом следует отметить, что колебательный режим кристаллизации осуществляется только при определенных соотношениях между физикохимическими параметрами расплава и условиями его охлаждения, на что указывают фазовые кривые, характеризующие режим кристаллизации. Количественно такое условие выражается формулой (37), причем начальные условия задачи должны соответствовать u0 -

- u1 - R-у/R), u0 - 0. В зависимости от тепломассообменных и других кристаллизационных параметров возможны и такие ситуации, когда возникшие случайные возмущения в системе имеют неустановившийся колебательный характер. В этом случае система приходит в прежнее состояние (устойчивый режим) за время релаксации, либо перейдет в новое стационарное состояние (неустойчивый режим).

4. Заключение

Исследован процесс направленного затвердевания псевдобинарных систем с учетом неравновесной кинетики кристаллизации. Для этого использован подход, основанный на предположении о том, что основные процессы формирования структуры кристаллизующегося сплава происходят в активном районе двухфазной зоны, прилегающем к фронту затвердевания. Проведенное исследование позволило выявить основные условия развития ритмического режима кристаллизации. Показано, что характер автоколебаний в системе существенно зависит от соотношения тепло- и массообменных параметров, степени неравновесности процесса, характеризующейся значениями кинетических констант K, К0, скорости кристаллизации и градиентов температур в зоне. Предложенное решение позволяет определять пределы варьирования этих величин. Кроме того, численные эксперименты, проведенные для псевдобинар-ной магматической системы Di-An, демонстрируют

возможность количественного описания полосчатой

неоднородности в затвердевших магматических расплавах.

Литература

1. Борисов В.Т. Кристаллизация бинарного сплава при сохранении устойчивости // ДАН СССР. - 1961. - Т. 136. - № 3. - С. 583-587.

2. Моризейн К., Витт А., Гейтес Х. Распределение примесей в монокристаллах. - М.: Мир, 1968. - С. 25-261.

3. Херл Д. Колебания температуры в расплавленные металлах и их связь со слоистым распределением примесей в кристаллах, выращиваемых из расплавов // Проблемы роста кристаллов. - М.: Мир, 1968. - С. 200-215.

4. Сулимцев И.И., Матвеев Ю.Е., Борисов В.Т. Исследование диффузионного (концентрационного) переохлаждения в двухфазной зоне сплавов олово-цинк // Изв. АН СССР. Металлы. - 1974. - № 6. -С. 206-210.

5. Самойлович Ю.А. О возможности кристаллизации расплава в режиме автоколебаний // Теплофизика высоких температур. -1979. - Т. 17. - № 5. - С. 992-996.

6. Самойлович Ю.А. Системный анализ кристаллизации слитка. -Киев: Наукова думка, 1983. - 246 с.

7. Ландау А.И. К вопросу о волнообразном характере распределения

примеси вдоль длины растущего монокристалла // ФММ. - 1958. -Т. 6. - № 1. - С. 148-155.

8. Anderson A.T Probable relation between plagioclase zoning and magma

dynamics, Fuego Volcano, California // Amer. Mineral. - 1984. -V. 69. - P. 660-676.

9. Tsushyama A. Dissolution kinetics of plagioclase in melt of the system diopside-albite-anorthite and the origin of dusty plagioclase in ande-sites // Contrib. Miner. Petrol. - 1985. - V. 89. - P. 1-15.

10. Шарапов В.Н., Черепанов А.Н. Динамика дифференциации магм. - Новосибирск: Наука, 1986. - 186 с.

11. McBirney A.R. Mechanisms of differentiation in the Scaergaard intrusion // J. Geol. Soc. London. - 1995. - V. 152. - P. 421-554.

12. Френкель М.Я. Тепловая и химическая динамика дифференциации базитовых магм. - М.: Наука, 1995. - 216 с.

13. Vistelius A.B., Pavlov V.M. On rhythmical layering of rocks formed from basaltic magma // Mathem. Geol. - 2003. - V. 35. - No. 4. -P. 399-404.

14. Ярошевский А.А. «Псевдоритмичность» как результат случайные событий (к вопросу о ритмичности расслоенных магматических комплексов) // Геохимия. - 2007. - № 2. - С. 224-228.

15. Минералы. Т. 2. - М.: Наука, 1974. - 415 c.

16. ФлемингсМ.С. Процессы затвердевания. - М.: Мир, 1977. - 423 с.

17. Крыглов В.И., Бобков В.В., Монастыгрныгй П.И. В^гчислительные методы. Т. 2. - М.: Наука, 1976. - 399 с.

18. Эльсгольц Л.Э. Дифференциальные уравнения и вариационное исчисление. - М.: Наука, 1969. - 424 с.

Поступила в редакцию 20.05.2008 г., после переработки 01.06.2009 г.

Сведения об авторе

Черепанова Вера Корнилиевна, к.ф.-м.н., доц. НГТУ, vera_cherepanova@mail.ru