Gravitational differentiation in the regimes from stokes settling to Rayleigh–Raylor flows

Full Text

1. ВВЕДЕНИЕ

1.1. Состояние остывающей Земли
на этапе магматического океана

Земля образовалась за первые 50 млн лет существования Солнечной системы в результате соударения тел. Ядро Земли в основном сформировалось уже на этом первом этапе во время аккреции вещества в процессе сепарации железа от силикатов. На начальной стадии, благодаря энергии соударений тел, энергии радиоактивного распада короткоживущих элементов и гравитационной энергии, большая часть Земли была расплавленной. Несмотря на то, что в ней происходила интенсивная турбулентная конвекция, размешивающая вещество, благодаря большой разнице в плотности и малой вязкости железо все-таки могло эффективно осаждаться. По­этому ядро росло одновременно с образованием Земли.

По мере охлаждения Земли в последующие 100 млн лет на втором этапе Земля и ядро продолжали еще несколько расти. На этом этапе кардинально менялись гидродинамические, физические и химические процессы в мантии Земли, пока не установилась регулярная тепловая конвекция с литосферными плитами и мантийными плюмами [Solomatov, 2015].

При конвекции в мантии устанавливается распределение температуры, близкое к адиабатическому. Адиабатическая температура в мантии растет с глубиной медленнее по сравнению с температурой плавления, поэтому затвердевание мантии происходило снизу. Поскольку силикаты являются химически многокомпонент­ными, то граница затвердевания не резкая, а имеется переходная зона частичного плавления, в которой доля кристаллов растет от 0% до 100%. Однако внутри этой зоны имеется рео­логическая граница при доле кристаллов 60%. При таком содержании расплава кристаллы начинают образовывать каркас, в результате чего вязкость скачком возрастает на несколько порядков.

При остывании Земли в низах мантии образовывался полутвердой слой, вязкость и толщина которого росли по мере остывания Земли. Поэтому при остывании Земли мантия на этом втором этапе была двухслойной. Нижняя высоковязкая часть мантии отделялась от верхней маловязкой части мантии реологической границей. Верхняя часть мантии представляла собой магматический океан, в который сверху продолжало поступать вещество планетозималей, содержащее силикаты и железо. В магматическом океане доля расплавленных силикатов превышала 40%, вязкость была малой (порядка 0.1 Пас), и имела место интенсивная конвекция, которая размешивала и замедляла осаждение как мелких капель железа, так и кристаллов силикатов.

Эта конвекция в магматическом океане обеспечивала быстрый вынос тепла из нижней части мантии. Поэтому охлаждение и повышение вязкости нижней части мантии происходили очень быстро, за время порядка 1000 лет. Вынос же тепла из самой верхней части мантии тормозился толстой первичной корой, поэтому ее затвердевание длилось значительно дольше, порядка 100 млн лет [Solomatov, 2015]. Поскольку вязкость нижней части мантии имела порядок 1014 Пас, то на ее верхней границе скапливалось железо, опускавшееся в виде мелких капель из верхней части мантии. В результате на реологической границе верхней и нижней частей мантии возникал слой жидкого железа. Поскольку плотность железа больше плотности силикатов, то благодаря неустойчивости Релея–Тейлора жидкое железо из этого слоя проникало в нижнюю мантию. Оно эффективно осаждалось в виде каплевидных диапиров, скорость осаждения которых превышала скорость конвективных течений.

Вероятное состояние Земли на втором этапе эволюции Земли в присутствии магматического океана по данным [Walter, Trønnes, 2004; Rubie et al., 2015; Solomatov, 2015; Halliday, Wood, 2015] схематически представлено на рис. 1.

Рис. 1. Схематическая картина состояния ранней остывающей Земли на второй стадии роста ядра и затвердевания магматического океана в первые сотни млн лет: T0 – современная геотерма; T – геотерма ранней Земли; Tm (40%) – температура частичного плавления перидотита с 40%-ным расплавом.

Процессы осаждения железа, имевшие место на второй стадии остывания Земли с магматическим океаном и завершающимся этапе роста ядра, очень сложные и изучены пока лишь частично [Solomatov, 2015]. Интенсивная тепловая конвекция перемешивала кристаллы силикатов и капли железа. Однако и железное ядро быстро росло, и силикатные кристаллы осаждались. Это было возможным благодаря тому, что скорость осаждения превышала скорость конвективных течений. При этом гравитационная дифференциация в разных слоях (различающихся физико-химическими свойствами) происходила в резко отличающихся режимах. В магматическом океане осаждались капли железа размером порядка 1 см. На границе слоя железа и нижней мантии в результате неустойчивости Рэлея-Тейлора возникали струи железа, которые по мере осаждения разбивались на отдельные каплевидные диапиры размером порядка 1 км. Благодаря гидродинамическому взаимодействию оседающие капли и диапиры могли группироваться.

1.2. Седиментационная конвекция в суспензии

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

При инверсии плотности в жидкости может возникать конвективное (циркуляционное) движение. Конвекция в однокомпонентной нагреваемой жидкости называется тепловой. В лабораторном эксперименте она была открыта Бенаром, а теория ее возникновения (теория неустойчивости) была развита Рэлеем. Конвекцию, вызванную инверсией плотности в изотермической многокомпонентной жидкости, в общем случае называют химической. В частности, в отсутствие химических реакций для механической смеси ее называют гравитационной или концентрационной, или композиционной. При диффузионном механизме переноса массы уравнения композиционной (диффузионной) конвекции совпадают с уравнениями тепловой конвекции с точностью до замены диффузионного потока тепла q_T = –k_T · ∇T на диффузионный поток массы q_D = –k_D · r∇C, где: k_T – коэффициент теплопроводности; T – температура; k_D – коэффициент диффузии; C – объемное содержание (концентрация) примесной компоненты; r – ее плотность. Конвекцию, вызванную градиентами температуры и концентрации с диффузионным механизмом, называют двойной диффузионной конвекций.

Однако, кроме диффузионного механизма (на молекулярном уровне), в гетерогенной механической смеси возможен также перенос массы на макроуровне с ненулевой относительной скоростью примеси относительно основной жидкости. Таким примером может служить падение совокупности тяжелых частиц (или всплывание легких) сквозь вязкую жидкость по закону Стокса с относительной скоростью S. При этом распределение совокупности отдельных частиц можно усреднить и рассматривать их как примесную жидкую компоненту, которая может проникать сквозь основную жидкость. Особенностью такого переноса массы является то, что в случае осаждения частиц поток массы будет пропорционален не градиенту концентрации C как в диффузионной конвекции, а самой концентрации и равен qs = S · rC. Здесь коэффициент пропорциональности равен относительной скорости S примесной компоненты. В результате этого дифференциальное уравнение переноса массы в отличие от диффузионного переноса тепла или массы меняется и становится не второго, а первого порядка.

Если инверсия плотности достаточно велика, то, кроме кондуктивного переноса массы (путем стоксового осаждении частиц), может возникнуть более эффективный перенос массы с конвективным движением всей смеси. Эту конвекцию можно назвать седиментационной [Трубицын, Харыбин, 1987], чтобы более подчеркнуть ее отличие от тепловой конвекции или диффузионной, которые также вызываются силой гравитации. В общем случае при массопереносе композиционная конвекция может включать в себя и диффузионную, и седиментационную конвекцию.

Исследование неустойчивости слоев тяжелой и легкой невязкой жидкости при инверсии плотности было проведено Рэлеем в 1900 г. В 1908 г. Харрисон [Chandraskhar, 1961] получил уравнение для роста мод возмущений для слоев вязких жидкостей. В 1950 г. Тейлором теория неустойчивости слоев жидкостей в гравитационном поле была распространена на течения в центрифугах. Чандрасекхаром теория неустойчивости Рэлея–Тейлора была обобщена учетом магнитного поля и вращения [Chandraskhar, 1961]. В последующем [Joseph, 2003] был развит метод потенциального течения, позволяющий приближенно, но более просто рассчитывать растущие моды при неустойчивости Рэлея–Тейлора.

Поскольку в теории неустойчивости Рэлея–Тейлора не учитывается возможность проникновения примесной тяжелой жидкости сквозь основную жидкость, то рассматриваемая в настоящей работе седиментационная конвекция обобщает теорию Рэлея–Тейлора, которая становится предельным частным случаем седиментационной конвекции при стремлении к нулю размеров частиц и, соответственно, скорости их падения сквозь основную жидкость. В этом случае совокупность частиц, вмороженных в основную жидкость, и основную, свободную от частиц, часть жидкости, можно рассматривать как две жидкости с разными плотностями. В работе [Pan et al., 2001] проведено численное моделирование неустойчивости очень плотного слоя частиц, оседающих в вязкой жидкости в предельном случае, в окрестности перехода к неустойчивости Рэлея–Тейлора. В настоящей работе развит метод, позволяющий исследовать весь диапазон режимов осаждения слоя частиц от Стокса до Рэлея–Тейлора, в частности, детально в окрестности перехода к режиму Стокса, с использованием понятия седиментационная конвекция.

Исследование проводилось в приближении не очень быстрых течений, когда инерционными членами в уравнении Стокса можно пренебречь. Постановка проблемы седиментационной конвекции и аналитический анализ конвективной неустойчивости (неустойчивости с возникновением конвекции) (onset of convection problem) в частном случае линейного инверсионного распределения концентрации были проведены еще в работах [Трубицын, Харыбин, 1987]. Авторами были найдены критические числа седиментационного аналога теплового числа Рэлея Rs и критические моды конвективного движения.

В обычной постановке рассматривается неустойчивость покоя или стационарного движения с возникновением конвекции. В настоящей работе рассматривается неустойчивость нестационарного, меняющегося во времени, движения. Такую неустойчивость можно назвать как мгновенная конвективная неустойчивость (instant onset of convection problem). Поскольку при Стоксовом оседании слоя тяжелых частиц общее распределение плотности в суспензии постоянно меняется, то и распределение скоростей смеси постоянно меняется. Однако в любой момент при включении возмущения это движение смеси со Стоксовым осаждением может смениться на конвективное течение.

Изучение конвективного переноса тепла началось в начале прошлого века. При этом для понимания ее природы Рэлеем сначала была исследована конвективная неустойчивость однокомпонентной нагретой жидкости. Интенсивное изучение конвективного массопереноса проводилось позже, когда само явление конвекции уже стало понятным. При этом уже была построена теория неустойчивости Рэлея–Тейлора для двухслойной вязкой жидкости с инверсией плотности, которая, однако, является лишь частным предельным случаем седиментационной конвекции. Со второй половины прошлого века началось интенсивное изучение конечно-амплитудных течений переноса массы на численных моделях. Таким образом, изучение тепловой конвекции сначала началось с построения теории неустойчивости и слабо развитой конвекции, а изучение массопереноса стали сразу изучать с конечно-амплитудного движения, включая процессы осаждения частиц. При этом теория конвективной неустойчивости для массопереноса не была построена. Причиной этого, возможно, были или не востребованность такой теории, или, возможно, необычность постановки задачи и ее решения. Однако процессы осаждения частиц в суспензии, в частности, в режиме перехода от Стоксового осаждения к седиментационной конвекции, наблюдаются не только в экспериментах, но и в технике, и в быту.

К настоящему времени опубликовано очень много работ по численному моделированию развитой гравитационной конвекции в дисперсных системах с учетом осаждения частиц. В первых работах вводились упрощения с пренебрежением объемным содержанием частиц и вариациями эффективной вязкости среды [Hill et al., 1977]. Многочисленные примеры расчетов конвекции в суспензиях с использованием более общей модели, частично учитывающей конечность объема дисперсной фазы, представлены в монографии [McCaffrey, 1997]. Численные решения уравнений развитой композиционной конвекции (при усредненном описании дискретных частиц в виде непрерывного распределения) в применении к осаждению кристаллов в застывающих магматических камерах получены в работе [Rudman, 1992].

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

2. МОДЕЛЬ, УРАВНЕНИЯ ДВИЖЕНИЯ И МАССОПЕРЕНОСА И ГРАНИЧНЫЕ УСЛОВИЯ ДЛЯ СУСПЕНЗИИ

Рассмотрим слой вязкой двухкомпонентной жидкости, состоящей из основной жидкости с плотностью r0 и небольшого количества тяжелой примеси частиц с плотностью вещества r1. Дискретное распределение частиц будем описывать как усредненное непрерывное распределение тяжелой компоненты с объемной долей (концентрацией) C = Ω₁/(Ω₁ + Ω₀). Здесь Ω₁ – объем элемента смеси, занятый тяжелыми частицами и Ω₀ – объем, занятый основной жидкостью. При этом будем учитывать, что тяжелая компонента может двигаться (проникать) сквозь основную жидкость с относительной скоростью S = V₁ – V₀. Здесь V₁ и V₀ – абсолютные скорости движения тяжелой компоненты и основной жидкости, S – относительная скорость тяжелой компоненты.

При стоксовом осаждении одиночного шара с вязкостью вещества η1 в неподвижной жидкости с вязкостью η0 относительная скорость S приближенно равна [Ландау, Лифшиц, 1986]:

S_s = [(2gΔρr²)/(3η₀)] [(η₀ + η₁)/(2η₀ + 3η₁)]. (1)

При η₁ >> η₀ (для твердого шара) выражение (1) переходит в формулу Стокса Ss = 2gΔrr 2/(9η), а при η1 << η0 (для жидкого шара) выражение (1) переходит в Ss = gΔρr²/(3η₀). Здесь: r – радиус шара, g – ускорение силы тяжести, Δr = r1–r0 – скачок плотности.

В двумерной модели (для бесконечного цилиндра, падающего плашмя), скорость падения цилиндра S_с определяется приближенным соотношением [Слезкин, 1986]:

S_с = ΨgΔρr²/(4η), y(β) = ln β–(β 2–1)/(β 2 + 1), β = β/r, (2)

где: r – радиус цилиндра; g – ускорение силы тяжести; y – множитель, характеризующий влияние боковых стенок с прилипанием; β – расстояние от цилиндра до стенок. Для значений β = 3, 8, 10 функция y принимает, соответственно, значения 0.3, 1.0 и 1.3. При β = 6 скорость падения цилиндра (2) совпадает с (1) для шара.

При осаждении конкретной частицы, с учетом взаимодействия с другими частицами благодаря противотоку жидкости, ее скорость S¢ уменьшается (увеличивается эффективная вязкость), что описывается введением множителя S¢ = xS, зависящего от концентрации x = x(C). Эта зависимость x пока плохо известна и остается предметом многих исследований. По формуле Эйнштейна [Ландау, Лифшиц, 1986] x = 1/(1 + 2.5C). В работе [Richardson, 1954] предлагается более сильная зависимость x от концентрации x = (1–C) 5. Очевидно, что благодаря взаимодействию с другими частицами может меняться и траектория частицы. Изложенная ниже теория седиментационной конвекции показывает, как эти траектории определяются течениями седиментационной конвекции, в пределе переходящими в течения Рэлея–Тейлора. Таким образом, относительная скорость тяжелой компоненты в общем случае может быть функцией концентрации и скоростей течений основной жидкости S = S (C, V0).

Уравнение движения Навье–Стокса для смеси и уравнения массопереноса для смеси и отдельно для тяжелой компоненты можно записать в обычном виде [Ландау, Лифшиц, 1986]

rdV/dt = –∇p+η∇²V + ρg, (3)

δr/δt + div(ρV) = 0, (4)

δC/δt + div(CV₁) = 0. (5)

Здесь: p – давление, η – постоянная эффективная вязкость смеси, ρ – плотность смеси, равная

1/ρ = Cρ₁ + (1–C)ρ₀ = ρ₀(1 + CΔρ/ρ₀). (6)

V – скорость движения центра тяжести элемента смеси, равная

V = C_mV₁ + (1–C_m)V₀, (7)

где: V₁ и ρ₁ – скорость и плотность для тяжелой компоненты V₀ и ρ₀ – для основной жидкости, Δρ = ρ₁–ρ₀, C_m – массовая концентрация. Абсолютная скорость тяжелой компоненты V1 равна векторной сумме относительной скорости S и скорости основной жидкости V₀, V₁ = S + V₀.

Подставим (6) и (7) в (4) и (5). Для малых концентраций C<<1 и C_m<<1, полагая ρ₁ ≈ ρ₀ во всех членах уравнений, кроме члена плавучести, и V₁ = V₀ = V во всех членах уравнений, кроме переноса массы тяжелой компоненты, после простых преобразований получим систему уравнений:

ρdV/dt = –∇p + η∇²V + g(1 + CΔρ/ρ₀)ρ₀, (8)

divV = 0, (9)

δC/δt = –div[C(S + V)]. (10)

Отметим, что в приближении малых концентраций, пренебрегая взаимодействием осаждающихся частиц между собой, скорость осаждения Стокса S можно считать независящей от кон­центрации С. Далее, для медленных течений можно пренебречь и зависимостью скорости осаждения частиц S от распределения скоростей основной жидкости, полагая ее такой же, как и в неподвижной жидкости. Таким образом, в рассматриваемом приближении относительную скорость тяжелой компоненты можно считать постоянной и равной скорости Стокса, известной для рассматриваемых частиц заданной формы.

Система трех уравнений (8)–(10) содержит три неизвестные функции: гидродинамическую скорость смеси V = (U, V, W), распределение объемной доли частиц C и давление p.

Рассмотрим двумерный горизонтальный слой жидкости толщиной D. Введем безразмерные переменные. Примем в качестве единицы измерения для длины – толщину слоя D, для скорости V₀ – заданное значение относительной скорости частицы S₀ (например, скорость падения шаров по Стоксу), для времени – время падения частицы с заданной скоростью S₀ через слой неподвижной жидкости t0 = D/S0, для концентрации частиц – некоторое произвольное значение объемного содержания частиц C = C₀ (например, C₀ = 10% = 0.1), для давления p₀ = ηS₀/D.

Направим ось z вверх. Поскольку скорость Стокса S направлена вниз, то div(СS) = –SδС/δz. Для упрощения скорость S принята постоянной. Будем отсчитывать давление от гидростатического распределения, определяемого уравне­нием:

∇p₀ = –gρ₀. (11)

В безразмерных переменных уравнения (8)–(10) примут вид:

(1/Ps)ρdV/dt = –∇p + ∇²V – R_S ∙ C ∙ i_z, (12)

divV = 0, (13)

δC/δt+V ∙ ∇C = δC/δz, (14)

где: R_S – седиментационное число, аналог тепловому числу Рэлея; P_s – седиментационный аналог числа Прандтля P_r или числа Шмидта S_c (диффузионного аналога числа Прандтля).

R_S = ΔρC₀gD²/(S₀η), P_s = ν/(SD). (15)

Седиментационное число характеризует интенсивность течений при массопереносе подобно тепловому числу Рэлея для тепловой конвекции Ra = α∆TgD³/(kη), где: α – коэффициент теплового расширения; k – коэффициент теплопроводности; k = k/(ρc_p) – коэффициент температуропроводности; ∆T – перепад температуры в слое; ν = η/ρ – кинематическая вязкость. Различная зависимость чисел R_S и R_a от толщины слоя обусловлена тем, что поток тепла пропорционален градиенту температуры q_T = k_T ∙ ∇T, а поток тяжелой компоненты при стоксовом падении пропорционален самой концентрации q_s = S ∙ ρC.

Число Прандтля Pr = ν/k характеризует соотношение между интенсивностями диффузии импульса (благодаря вязкости) и диффузии тепла. При двойной диффузии число Шмидта (диффузионный аналог числа Прандтля) S_c = ν/k_D (где k_D – коэффициент диффузии) характеризует соотношение между интенсивностями диффузии импульса и диффузии массы. Появляющийся в теории седиментационной конвекции седиментационный аналог числа Прандтля Ps = ν/(SD) характеризует соотношение между интенсивностями диффузии импульса и седиментационного переноса массы тяжелой компоненты.

Отметим, что уравнение для концентрации (14) можно также получить непосредственно из закона сохранения массы примеси в виде d(ρC)/dt = divq_s,

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

Для исключения давления из системы уравнений (12), (13) применим операцию rotrot к уравнению (12). Тогда с учетом (13) система уравнений (12)–(14) сведется к двум уравнениям для двух переменных W (x, z, t) и C (x, z, t):

∇⁴W = R_Sδ²C/δx², (16)

δC/δt + V ∙ ∇C = δC/δz. (17)

Примем простейшие граничные условия для скорости, взяв верхнюю (z = 1) и нижнюю (z = 0) границы непроницаемыми W = 0 без прилипания, т. е. при нулевом сдвиговом вязком напряжении σ_zx = η(δU/δz + δW/δx) = 0. Так как W равно нулю при всех x, то δW/δx = 0. В результате граничные условия для скорости примут вид:

W = δ² W/δz² = 0 при z = 0, 1. (18)

Принципиальное отличие уравнения седиментационного массопереноса от уравнения теплопереноса выражается в том, что дифференциальное уравнение (17) для объемной доли частиц C (в отличие от уравнения для температуры T) имеет первый порядок. Поэтому для концентрации нужно задавать только одно граничное условие, например, фиксированное значение C₁ на верхней границе:

C = C₁ при z = 1. (19)

Система уравнений седиментационной конвекции (16), (17) для безразмерных переменных содержит всего один параметр – седиментационное число R_S.

Если в качестве начального взять какое-либо распределение тяжелых частиц без инверсии распределения плотности для всей смеси, то частицы будут равномерно осаждаться без возникновения конвективных течений. Противоток вытесняемой вверх основной жидкости не представляет собой конвекцию, так как не образует замкнутых (или квазизамкнутых) траекторий течения жидкости. Однако при инверсии начального распределения плотности может возникать конвекция. Подобно тому, как при нагреве снизу возникает тепловая конвекция, так при возрастании вверх распределения тяжелых частиц может возникнуть конвекция во всей смеси, которую можно назвать седиментационной конвекцией [Трубицын, Харыбин, 1987].

3. ЧИСЛЕННОЕ МОДЕЛИРОВАНИЕ ПРОЦЕССА ОСАЖДЕНИЯ ОДНОРОДНОГО СЛОЯ ТЯЖЕЛЫХ ЧАСТИЦ С ВОЗНИКНОВЕНИЕМ СЕДИМЕНТАЦИОННОЙ КОНВЕКЦИИ

При осаждении слоя частиц в вязкой жидкости мгновенное распределение концентрации частиц имеет вид ступенчатой функции. Аналитическое решение задачи на неустойчивость ниже будет проведено для случая инверсионного линейного распределения оседающих частиц, когда уравнения для возмущений очень упрощаются. При численном решении можно рассмотреть осаждение частиц с любыми распределениями концентрации тяжелых частиц. Рассмотрим наиболее часто встречающийся случай осаждения однородного слоя частиц с фиксированной концентрацией. Значение безразмерной объемной доли частиц было взято равным C = 1 (или C = C₀ в размерном виде) внутри слоя и C = 0 вне слоя. Боковые стенки взяты непроницаемыми и без прилипания. Одна из границ слоя частиц (нижняя) взята неровной, чтобы создать начальное возмущение для развития конвекции, а также продемонстрировать зависимость скорости нарастания возмущения от номера волнового числа k.

Численные решения уравнений седиментационной конвекции (12)–(14) были получены для различных значений седментационных чисел как для малых Rs (в случае, когда слой, не деформируясь, оседает по Стоку), так и для больших Rs, когда нелинейный член V ∙ ∇C, учитывающий обратное влияние изменяющегося распределения тяжелых частиц на скорости течения, становится большим и возникает развитая седиментационная конвекция, деформирующая слой тяжелых частиц.

Численное решение проводилось с помощью программы в конечных элементах СиткомCU (CITCOM-CU), усовершенствованной в статьях [Moresi, Solomatov, 1995; Moresi et al., 1996], переработанной и дополненной, в частности, автоматизированной графикой, [Евсеев, 2008].

На рис. 2 в левой колонке представлены результаты расчета для осаждения крупных частиц со скоростью Стокса S = S₁ при значении седиментационного числа R_S1 = 80. При таком подкритическом значении числа Рэлея в смеси жидкой смеси даже при инверсии плотности конвекция не возникает. Слой частиц опускается как целое, не деформируясь, со скоростью осаждения отдельной частицы в безразмерных единицах S₁.

Во второй колонке приведены расчеты эволюции жидкой смеси с мелкими частицами, безразмерная скорость оседания Стокса для которых равна S₂ = 0.08, что соответствует значению седиментационного числа R_S2 = 1000. Как видно на рис. 2, несмотря на то, что скорость Стокса для частиц в правой колонке в 12.5 меньше, чем в левой колонке, слой достигает дна всего в три раза дольше. При этом возникают течения седиментационной конвекции, характеризуемые числом R_S = 1000, структура которых уже близка к структуре течений Рэлея–Тейлора.

Таким образом, седиментационная конвекция убыстряет перенос массы, особенно для мелких частиц. В пределе одинокие очень мелкие частицы, имеющие очень малую скорость стоксового осаждения, висят как бы вмороженными в основную жидкость. Но, поскольку эти частицы все-таки тяжелее основной жидкости, то в совокупности они эквивалентны слою тяжелой жидкости. Такой предельный случай эквивалентен классической задачи Рэлея–Тейлора. Таким образом, развиваемая в работе теория седиментационной конвекции включает теорию Рэлея–Тейлора как частный, предельный случай.

Рис. 2. Результаты численного моделирования осаждения слоя частиц в вязкой жидкости при различных скоростях Стокса и, соответственно, разных значений седиментационного числа. В левой колонке слой осаждается как целое со скоростью Стокса. В правой колонке при малой скорости осаждения частиц возникают течения, переходящие в течения Рлея–Тейлора.

Как видно на рис. 2, по данным этого моделирования критическое значение седиментационного числа для перехода режима осаждения слоя частиц по Стоксу к режиму конвекции находится в интервале значений от 80 до 1000.

4. ЗАВИСИМОСТЬ СКОРОСТИ РОСТА МОД ВОЗМУЩЕНИЙ ОТ ВОЛНОВОГО ЧИСЛА ДЛЯ ОСЕДАЮЩЕГО СЛОЯ В КРИТИЧЕСКОМ РЕЖИМЕ

Чтобы более точно определить значение критического седиментационного числа для слоя осаждающихся частиц и выявить наиболее быстро растущую моду, далее было проведено численное решение уравнений (12)–(14) для R_s = 100 и R_s = 200, при этом для различных начальных мод возмущений.

На рис. 3 приведены результаты решения для осаждения слоя при значении седиментационного числа R_S = 100. В левой колонке начальное положение возмущенной нижней границы слоя взято в виде z = 0.2 + 0.03sin(πx), а в правой z = 0.2 + 0.03sin(3πx). Как видно на рис. 3, при осаждении слоя частиц возмущения при k = π и при k = 3π почти не растут, так что слой практически не деформируется во время падения. Появляются скорости течений, обусловленные неровностью границы, но эти скорости не растут, и деформации слоя не увеличиваются.

На рис. 4 приведены результаты решения для осаждения слоя при значении седиментационного числа R_S = 200 при различных модах возмущения с волновыми числами k = π, k = 3π, k = 6π. При числе R_S = 200 возмущения явно растут. При этом быстрее растет мода k = π и начинает несколько расти мода k = 3π. Для роста более высоких мод нужно более высокое, свое критическое значение седиментационного числа.

Таким образом, из сравнения рис. 3 и рис. 4, для оседающего слоя (что соответствует ступенчатому распределению концентрации) седиментационная конвекция начинается при значении седиментационного числа, примерно равном R_S* = 100–200. Оценку в качестве критерия возникновения седиментационной конвекции при осаждении тяжелых частиц в суспензии, когда начинают явно расти первые моды возмущений, можно принять приближенно в виде R_S*≈200.

Рис. 3. Осаждение слоя тяжелых частиц при Rs = 100 для двух мод начального возмущения нижней границы при k = π и k = 3π.

Рис. 4. Осаждение слоя тяжелых частиц при Rs = 200 для трех мод начального возмущения нижней границы при k = π, k = 2π и k = 6π.

5. АНАЛИТИЧЕСКОЕ РЕШЕНИЕ ДЛЯ МГНОВЕННОЙ КОНВЕКТИВНОЙ НЕУСТОЙЧИВОСТИ ВЯЗКОЙ ЖИДКОСТИ С ЛИНЕЙНЫМ РАСПРЕДЕЛЕНИЕМ КОНЦЕНТРАЦИИ ОСАЖДАЮЩИХСЯ ТЯЖЕЛЫХ ЧАСТИЦ

В рассмотренной выше модели осаждения однородного слоя тяжелых частиц распределение концентрации описывается сдвигающейся во времени ступенчатой функцией. При этом в результате мгновенной неустойчивости возникает конвекция, при которой в каждый момент меняется ее структура по мере осаждения слоя. Задача аналитического определения критических условий возникновения седиментационной конвекции очень упрощается, если в качестве начального инверсионного распределения тяжелых частиц взять не ступенчатую, а простейшую линейно возрастающую кверху функцию C (z, t = 0) = C₁ – γ ∙ (1–z), где γ – наклон прямой. Эта функция фактически размазывает скачок концентрации ступенчатой функции равномерно по всей глубине области.

При таком начальном условии, в отсутствии конвекции, решение уравнений (16)–(19) будет иметь вид (см. рис. 5):

C(z, t) = C₁ при z < 1–S ∙ t

C(z, t) = C₁–γ ∙ (1–S ∙ t–z) при z < 1–S ∙ t. (20)

Это решение соответствует равномерному осаждению облака частиц через неподвижную жидкость. Заметим, что в рассматриваемом приближении малой концентрации частиц мы пренебрегаем небольшим противотоком жидкости.

Однако такой режим осаждения слоя кристаллов при определенных условиях может нарушаться с возникновением конвективного движения всей жидкости [Трубицын, Харыбин, 1987]. Будем искать решение уравнений (16)–(19) в виде суперпозиции решения C (z, t) для осаждения частиц в неподвижной жидкости W = 0 и возмущения, растущего во времени, т. е. применим подход, развитый Рэлеем для анализа возникновения тепловой конвекции. В малом нелинейном члене V ∙ ∇C заменим распределение концентрации C на ее начальное линейное распределение C. Тогда ∇C ≈ ∇C = γi_z. В результате нелинейное дифференциальное уравнение (17) упрощается до линейного уравнения:

δC/δt + gW = С/z. (21)

Рис. 5. Начальное С(0) и мгновенное С(t) распределение осаждающейся тяжелой компоненты по глубине.

Будем искать решения уравнений (16), (21) для скорости движения жидкости W и для распределения тяжелой компоненты C в виде, соответствующем конвективным течениям (ана­логично теории тепловой неустойчивости Рэлея), в виде:

γW(x, z, t) = exp(λt) ∙ φ(z) ∙ cos(kx),

C(x, z, t) = C (z, t) + exp(λt) ∙ f(z) ∙ cos(kx), (22)

где: показатель λ характеризует скорость нарастания возмущений во времени; k – волновое число, характеризующее его длину волны 2π/k с периодичностью по оси x.

После подстановки (22) в уравнения (16), (21) с учетом (20) получим систему уравнений для функций φ(z) и f(z):

(d²/dz²–k²)²φ(z) = –k²R_sf(z),

λf(z) + φ(z) = df/dz. (23)

Исключая j(z) из (23), получим для f(z) дифференциальное уравнение пятого порядка, содержащее три параметра Rs, k и l:

(d²/dz²–k²)² (df/dz–λf)) = –k²Rsf(z). (24)

С граничными условиями:

f|_z=1= df/dz|_z=1= d²f/dz²|_z=1= df/dz|_z=0=d²f/dz²|_z=0=0. (25)

Будем искать частные решения уравнения (24) в виде экспоненциальной функции от z:

f(z) = A · exp(qz). (26)

Подставив (26) в уравнение (24), получим алгебраическое уравнение пятой степени для определения коэффициента q:

(q²–k²)²(q–λ)) = –k²Rs. (27)

Это уравнение имеет пять корней q1, … q5. Поэтому общее решение уравнения (24) можно записать в виде суммы пяти функций:

f(z) = ∑A_n · exp(q_nz). (28)

Пять коэффициентов A_n находятся из пяти граничных (25). Поскольку получающаяся система пяти алгебраических уравнений однородна, то она имеет ненулевые решения при условии равенства нулю ее определителя. В результате получается одно алгебраическое уравнение, связывающее параметры R_s, λ и k. В явном виде его можно записать как Rs=R_S(k, λ) или как λ=λ(R_S, k). По (22) коэффициент λ характеризует скорость нарастания возмущений. Они не затухают, если λ > 0. Для фиксированной моды (фиксированного волнового числа k*) уравнение λ = λ(R_S, k*) определяет критическое значение седиментационного числа R_S, при котором λ становится положительным и возникает седиментационная конвекция с данной модой k*. При этом из всех возможных мод k* должна реализоваться та мода течений k*_min, для которой критическое седиментационное число R_S* минимально.

Вычисления показали [Трубицын, Харыбин, 1987], что минимальное значение критического седиментационного числа равно R_S* = 105 и оно достигается для моды с волновым вектором k*=3.0. На рис. 6 приведена функция λ=λ(R*S, k*) при фиксированном значении k*=3.0. Как видно на рис. 6, этот коэффициент положителен при R_S*≥105. Таким образом, это значение седиментационного числа R_S* является критерием нарастания возмущений и начала седиментационной конвекции.

Минимальное критическое число Рэлея R_S*=105 согласно (25) получено для обеих скользких границ. При этом длина волны наиболее быстро растущей моды равна 2π/k*, а конвективная ячейка имеет аспектное отношение (отношение ширины к высоте), равное a* = π/k* = 1.05. Для обеих шероховатых границ (с полным прилипанием) R_S* = 220,8 и a* = 0.68. Для одной скользкой и другой шероховатой границ R_S* = 157.1 и a* = 0.83. Как известно, при тепловой конвекции с начальным линейным распределением температуры для разных типов границ критические числа Рэлея и аспектные отношения, соответственно, равны: Ra* = 657.5 и a* = 1.41; Ra* = 1707.8 и a* = 1.01; Ra* = 1100.7 и a* = 1.17.

Рис. 6. Зависимость скорости нарастания возмущения λ от седиментационного числа RS при k = k*.

Таким образом, при значениях седиментационного числа R_S, меньших критического R_S*, возмущения вида (22) затухают (λ<0), и совокупность тяжелых частиц с начальным линейным распределением должна оседать равномерно без возникновения конвективных течений. При R_S>R_S* в смеси жидкости и тяжелых частиц конвективные возмущения нарастают (λ>0), и должна возникать седиментационная конвекция, при этом скорость суммарного массопереноса частиц увеличивается, т. к. к стоксовому осаждению частиц сквозь основную вязкую жидкость добавляется их перенос конвективными течениями.

Выражение для седиментационного числа (15) можно получить качественно из баланса сил. При возмущении границы раздела между основной жидкостью и вышележащей смеси элемент жидкости, в котором объемная концентрация тяжелых частиц равна C, обладает отрицательной плавучестью f₁=a₁CΔρg и при движении вниз тормозится силой трения f₂ = a₂ηV/D². При установившемся течении обе силы должны быть равны друг другу. Поэтому скорость опускания рассматриваемого элемента будет равна V=(a₁/a₂)CΔρgD²/η. Движение вниз рассматриваемого элемента возмущенной жидкости не затухнет, если его скорость больше, чем скорость S, с которой тяжелые частицы по закону Стокса выпадают из слоя, V>S. В противном случае тяжелые частицы уйдут из рассматриваемого элемента жидкости, и он перестанет быть тяжелым. Таким образом, критерий возник­новения конвективных течений в жидкой смеси принимает вид V=(a₁/a₂)CΔρgD²/η>S или ΔρCgD²/Sη>(a/a1). Таким образом, в соответствии с определением седиментационного числа (15) конвекция не затухает, если ΔρCgD²/Sη=Rs больше некоторого критического значения.

Заметим, что поскольку теплоперенос опи­сывается дифференциальным уравнением второго порядка и существует стационарное состояние теплопереноса, то в теории неустойчивости Рэлея–Бенара в качестве начального состояния можно взять состояние покоя при линейном распределении температуры с кондуктивным потоком тепла. Тепловая конвекция возникает, когда (при малом коэффициенте теплопроводности и вязкости) конвективный перенос тепла более эффективный, чем кондуктивный. При седиментационном потоке массы, поток массы пропорционален не градиенту концентрации (как в тепловой конвекции), а самой концентрации. Поэтому массоперенос описывается дифференциальным уравнением первого порядка, для которого не существует стационарного решения. Поэтому исследуется мгновенная неустойчивость состояния суспензии, при которой рой частиц продолжает опускаться и одновременно начинает менять свою форму, и возникают конвективные течения во всей смеси.

6. ЧИСЛЕННЫЕ РЕШЕНИЯ СЛАБО НАДКРИТИЧНОЙ И РАЗВИТОЙ СЕДИМЕНТАЦИОННОЙ КОНВЕКЦИИ, ВЫЗВАННОЙ ИНВЕРСИОННЫМ ЛИНЕЙНЫМ РАСПРЕДЕЛЕНИЕМ КОНЦЕНТРАЦИИ ТЯЖЕЛЫХ ЧАСТИЦ

Чтобы проверить полученные выше результаты аналитического решения для конвективной неустойчивости, было проведено численное решение уравнений (12), (14) для модели с таким же указанным выше начальным распределением частиц, линейно возрастающим кверху С₀(z, t = 0) = С₁–γ ∙ (1–z) для двух значений седментационного числа R_S1 и R_S2, меньшего критического значения и больше критического значения R*_S.

На рис. 7 в левой колонке представлены результаты численного решения уравнений (12)–(17) для эволюции распределения тяжелых частиц с безразмерной скоростью оседания, равной S₁ = 1 (или S₁ = S₀ в размерном виде) при значении седиментационного числа, равного R_S1 = 95, что меньше критического значения, равного R*_S = 105. В этом случае в соответствии с результатами аналитического решения даже при слабом возмущении распределения конвекции не возникает, и рой тяжелых частиц оседает как целое по Стоксу со скоростью падения одиночной частицы. На рис. 7 в левой колонке приведено вычисленное распределение частиц для моментов безразмерного времени (в единицах D/S) t = 0, t = 0.2 и t = 0.3.

В правой колонке рис. 7 приведены результаты расчета для частиц несколько меньшего размера, для которых скорость Стокса равна S₂ = 0.85 (или S₂ = 0.85S0 в размерном виде). При такой скорости Стокса седиментационное число по (15) оказывается в 1/0.85 = 1.17 раз большим и равным R_S1 = 110, что больше критического значения R*_S = 105. В этом случае кондуктивный массоперенос путем простого стоксового осаждения становится малоэффективным, мгновенное распределение частиц в виде медленно оседающего неоднородного роя становится неустойчивым и подобно Релей–Тейлоровской неустойчивости возникают течения всей смеси. Скорости течений показаны стрелками. Соответственно, деформируется распределение тяжелых частиц. Чтобы легче сравнить положение верхней границы распределения частиц для обоих рассмотренных случаев, в правой колонке рис. 7 для более медленно оседающих частиц показаны моменты времени с несколько более длительными интервалами безразмерного времени t = 0, t = 2/0.85 = 2.35 и t = 3/0.85 = 3.53.

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

Рис. 7. Эволюция линейного распределения частиц, оседающих в вязкой жидкости. Концентрация тяжелых частиц показана темными тонами с интервалами 0.1C1. Левая колонка для частиц, оседающих в неподвижной жидкости в виде недеформирующегося распределения с безразмерной скоростью S₁ = 1 при значении седиментационного числа RS = 95 < RS* = 105 для моментов безразмерного времени t = 0, t = 0.2 и t = 0.3. Правая колонка – для частиц меньшего размера при меньшей скорости осаждения S₂ = 0.85 и, соответственно, большем седиментационном числе RS = 110 > RS* для моментов времени t = 0, t = 2.35 и t = 3.53. Скорости течений седиментационной конвекции показаны белыми стрелками.

7. ЗАКЛЮЧЕНИЕ

Классическими теориями неустойчивости в вязкой жидкости являются теория Рэлея–Бенара и теория Рэлея–Тейлора. При неустойчивости нагретой снизу однокомпонентной жидкости кондуктивный теплоперенос сменяется на регулярную тепловую конвекцию Рэлея–Бенара с ускорением теплопереноса. Неустойчивость слоя тяжелой изотермической жидкости, лежащего на легкой жидкости, приводит к возникновению струйных течений, описываемых теорией Рэлея–Тейлора. С одной стороны, поскольку струйные течения Релея–Тейлора отличаются от течений стационарной тепловой конвекции с регулярными ячейками, может показаться, что эти теории качественно различны. С другой стороны, течения Рэлея–Тейлора похожи на плюмы тепловой конвекции, и может показаться, что эти теории отличаются лишь обозначениями. Действительно, диффузия компонент вещества и диффузия тепла описываются с точностью до обозначений одинаковыми дифференциальными уравнениями второго порядка. Однако в отличие от уравнения теплопереноса, уравнение седиментации содержит только первую производную от концентрации и поэтому не имеет стационарного решения.

В теории Рэлея для тепловой конвекции рассматривается неустойчивость стационарного линейного распределения температуры, которое устанавливается в отсутствие конвекции. В настоящей работе развивается теория седиментационной конвекции и рассматривается неустойчивость не покоя, а мгновенная неустойчивость нестационарного, меняющегося во времени распределения концентрации (распределения тяжелой компоненты смеси). При этом акцентируется внимание на то, что неустойчивость Рэлея–Тейлора является частным предельным случаем седиментационной конвекции в суспензии. Течения Рэлея–Тейлора подобны течениям тепловой конвекции при очень высоком значении числа Рэлея. Если резко нагреть нижнюю часть слоя вязкой жидкости с очень малым коэффициентом теплопроводности, то появится инверсия плотности, и вначале возникнут струйные течения, подобные течениям Рэлея–Тейлора. При этом возбуждено очень много высоких мод. Со временем, по мере выравнивания скачка температуры, эти течения заменятся на течения тепловой конвекции Рэлея–Бенара. Аналогично при очень малых размерах частиц в облаке тяжелых частиц и, соответственно, малой скорости осаждения, частицы в суспензии оказываются вмороженными в жидкость и область, занятая частицами, оказывается аналогичной тяжелой жидкости. Поэтому появляется неустойчивость Рэлея–Тейлора, и возникают струйные течения жидкости с вмороженными частицами. При конечных размерах частиц с конечной скоростью стоксового осаждения в общем случае возникает седиментационная конвекция, аналогичная тепловой конвекции.

Следует отметить, что, как видно из выражения для седиментационного числа (15), оно растет с ростом концентрации и уменьшением скорости осаждения Стокса (зависящей от размера частиц). Поэтому предельный случай струйных течений Рэлея–Тейлора, когда Rs ® ∞, может достигаться при стремлении к нулю скорости осаждения частиц Стокса, т. е. при малом размере частиц. При этом, значение концентрации частиц в облаке может оставаться конечной. В частности, оно может быть достаточно малым, чтобы выполнялось приближение малых концентраций, используемое выше при выводе уравнений для седиментационной конвекции.

При указанном сходстве теплой и седиментационной конвекции имеется принципиальное различие между ними. Оно состоит в следующем. Если в нагреваемой жидкости создать скачок температуры, т. е. бесконечный градиент температуры и соответственно бесконечно большое значение локального числа Рэлея, то мгновенная конвективная неустойчивость на этой границе вызовет струйные течения типа Рэлея–Тейлора (до тех пор, пока не сгладится скачок температуры). Однако, как было показано выше на модели оседающего слоя частиц, несмотря на скачок концентрации, возникающие течения седиментационной конвекции являются регулярными и становятся струйными только в предельном случае плотного роя очень мелких частиц. Это обусловлено тем, что поток тепла, пропорциональный градиенту концентрации, стремится к бесконечности при скачке температуры. В то же время седиментационный поток массы, пропорциональный самой концентрации, всегда конечен.

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

Аналитическая теория и численные результаты получены в настоящей работе для случая высоковязких изотермических жидкостей при пренебрежении инерционными членами и в приближении описания роя дискретных частиц усредненным распределением и в приближениии малых концентраций. Поэтому для непосредственного применения к процессам осаждения железа в ранней Земле или в застывающих магматических камерах полученные результаты нужно дополнить учетом неинерциальных членов при немалом седиментационном числе Прандтля, а также учетом процессами перемешивания течениями тепловой конвекции.

About the authors

V. P. Trubitsyn

Author for correspondence.
Email: trub@ifz.ru
Russian Federation, Bolshaya Gruzinskaya str., 10-1, Moscow 123242, Russia

References

Supplementary files

Supplementary Files

Action

1. JATS XML

Download

2. Fig. 1. A schematic picture of the state of the early cooled Earth at the second stage of the growth of the nucleus and solidification of the magmatic ocean in the first hundreds of million years: T0 is a modern geotherm; T is the geotherm of the early Earth; Tm (40%) is the temperature of partial melting of peridotite with 40% melt.

Download (134KB)

Indexing metadata

3. Fig. 2. Results of numerical simulation of the deposition of a layer of particles in a viscous fluid at different Stokes velocities and, accordingly, different values of the sedimentation number. In the left column, the layer precipitates as a whole with the Stokes velocity. In the right column, at a low sedimentation rate of particles, flows occur, passing into the Rleigh – Taylor flows.

Download (243KB)

Indexing metadata

4. Fig. 3. The deposition of a layer of heavy particles at Rs = 100 for two modes of the initial perturbation of the lower boundary at k = π and k = 3π.

Download (570KB)

Indexing metadata

5. Fig. 4. Deposition of a layer of heavy particles at Rs = 200 for three modes of the initial perturbation of the lower boundary at k = π, k = 2π and k = 6π.

Download (574KB)

Indexing metadata

6. Fig. 5. The initial C (0) and the instantaneous C (t) distribution of the deposited heavy component in depth.

Download (8KB)

Indexing metadata

7. Fig. 6. Dependence of the growth rate of the perturbation λ on the sedimentation number RS at k = k *.

Download (7KB)

Indexing metadata

8. Fig. 7. The evolution of the linear distribution of particles deposited in a viscous fluid. The concentration of heavy particles is shown in dark tones with intervals of 0.1C1. The left column for particles deposited in a stationary liquid in the form of a non-deformable distribution with a dimensionless velocity S1 = 1 with a value of sedimentation number RS ​​= 95 RS * for moments of time t = 0, t = 2.35 and t = 3.53. The flow rates of sedimentation convection are shown by white arrows.

Download (190KB)

Indexing metadata