Кристаллография, 2019, T. 64, № 6, стр. 973-978

Гидродинамика и массоперенос при выращивании смешанных кристаллов из раствора

Н. А. Верезуб 1, А. Э. Волошин 2, А. И. Простомолотов 1*

1 Институт проблем механики им. А.Ю. Ишлинского РАН
Москва, Россия

2 Институт кристаллографии им. А.В. Шубникова ФНИЦ “Кристаллография и фотоника” РАН
Москва, Россия

* E-mail: prosto@ipmnet.ru

Поступила в редакцию 03.06.2019
После доработки 07.06.2019
Принята к публикации 10.06.2019

Полный текст (PDF)

Аннотация

Рассмотрены результаты математического моделирования оригинального процесса выращивания смешанных кристаллов K2(Co, Ni)(SO4)2 · 6H2O. Обсуждаются особенности гидродинамики и массопереноса при гладкой и шероховатой поверхностях кристаллизации.

ВВЕДЕНИЕ

Из анализа литературных данных следует, что конвекция в растворе приводит к ослаблению процесса образования включений и позволяет увеличить скорость роста без ухудшения качества кристалла [1].

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

Экспериментальное определение распределения поверхностного пересыщения затруднительно, поэтому большую роль играет численное моделирование течения и массопереноса при росте кристаллов из раствора. Например, в [3] исследовались зависящие от времени трехмерные течения при росте кристаллов KDP и было установлено их существенное влияние на рост кристаллов.

В [4] проведено двухмерное моделирование совместного действия вынужденной и естественной конвекции при росте KDP кристаллов, которое показало, что для подавления естественной конвекции требуется значительная интенсивность вынужденного течения в кристаллизационной камере. В [5] были выполнены трехмерные, зависящие от времени расчеты турбулентного течения применительно к условиям высокоскоростного выращивания кристаллов KDP. Они показали, что динамика течения и распределение пересыщения сильно зависят от размеров кристалла, скорости роста и скорости вращения кристалла. В [6] была предложена “самосогласованная” модель роста кристаллов KDP, в которой учитывалась как объемная диффузия, так и реакция на ростовой поверхности кристалла для определения толщины диффузионного слоя вокруг кристалла.

Оригинальная конструкция проточного кристаллизатора применяется для выращивания смешанных кристаллов K2(Co,Ni)(SO4)2 · 6H2O (KCNSH) из смеси двух водно-солевых растворов (кобальтовой K2Co(SO4)2 · 6H2O (KCSH) и никелевой K2Ni(SO4)2 · 6H2O (KNSH) солей) с разным мольным соотношением [7, 8]. Различие транспортных и физико-химических свойств разных солей может приводить к их неравномерному распределению в растворе вдоль грани кристалла, что, в свою очередь, будет приводить к неоднородности растущего смешанного кристалла KCNSH. Такой результат наблюдался в виде чередования более темных и светлых пятен на рентгеновском топографическом снимке поперечного сечения кристалла KCNSH, выращенного из раствора с равным мольным соотношением солей [8].

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

СХЕМА РОСТОВОГО ПРОЦЕССА

Суть ростового процесса состоит в следующем (рис. 1а): раствор (3) поднимается вертикально вверх из трубки (1) на дне цилиндрического сосуда (2) с заданной скоростью Vjet и обтекает поверхность растущего кристалла (4). Размеры кристаллизатора следующие: высота H = 0.04 м, диаметр D = 0.03 м. Диаметр трубки (1) с втекающим раствором d = 0.003 м. Осаждение солей на поверхность кристалла обеспечивает его наращивание сверху вниз по всему диаметру сосуда. Затем “отработавший” раствор вытекает через открытую часть дна сосуда (5).

Рис. 1.

Расчетная модель проточного кристаллизатора: а – трубка (1) для втекания раствора (3) с заданной скоростью Vjet в цилиндрический сосуд (2), поверхность растущего кристалла (4), вытекание раствора (5); б – форма треугольной шероховатости на поверхности кристалла.

Рассматривалась как гладкая поверхность кристалла, так и поверхность с треугольной формой шероховатости. Треугольники имеют одинаковую высоту (δy = 0.4 мм) и два радиальных размера: с большой шириной (δr1 = 2.7 мм) и меньшей (δr2 = 0.6 мм) вблизи оси (рис. 1б). Выбор размеров треугольников меньшего размера лимитировался вычислительной точностью при заданном числе узлов расчетной сетки.

Математическая модель сформулирована в соответствии с общепринятой методикой расчета процессов гидродинамики и массопереноса при выращивании кристаллов из водно-солевых растворов [3]. Данные, использованные в расчетах, приведены в табл. 1 и 2. Все известные для KCSH, KNSH и KCNSH данные взяты из [9]. Однако для рассматриваемых кристаллов неизвестна удельная энергия ступени α, в этом случае было взято значение для грани (100) кристалла KDP [10]. Также в литературе отсутствуют данные о коэффициентах диффузии солей KCSH и KNSH в водных растворах, поэтому в расчетах использованы значения для CoCl2 [11] и NiCl2 [12] соответственно. Различие в распределении концентраций никелевой и кобальтовой солей вдоль грани кристалла в значительной степени зависит от соотношения их коэффициентов диффузии. Можно полагать, что разница коэффициентов диффузии CoCl2 и NiCl2 близка к разнице коэффициентов диффузии KCSH и KNSH.

Таблица 1.  

Параметры гидродинамической модели

Кинематическая вязкость раствора ν, м2 ν = 5.06 × 10–7
Коэффициенты диффузии солей Di, м2 D1 = 1.217 × 10–9
D2 = 1.075 × 10–9
Плотность раствора ρ0, кг/м3 ρ0= 1115
Равновесные концентрации солей Cei (при T = 316 K), кг/кг р-ра Ce1= 0.044
Ce2= 0.089
Таблица 2.  

Параметры для задания скорости кристаллизации

T – температура [K] 316
σ0 – пересыщение солью при переохлаждении (ΔT = 0.5 K) 0.09
w – объем молекулы [м3] 3.24 × 10–28
k – константа Больцмана [Дж/K] 1.38 × 10–23
β – кинетический коэффициент ступени на грани (110) [м/с] 5 × 10–4
α – удельная энергия ступени [Дж/м2] 19.5 × 10–3
H – высота ступени на грани (110) [м] 5.3 × 10–10

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

МАТЕМАТИЧЕСКАЯ ФОРМУЛИРОВКА ЗАДАЧИ

Для определения вектора скорости V = (Vr, Vy) и давления P в растворе (3) (рис. 1) решаются уравнения Навье–Стокса и неразрывности, записываемые в векторном виде следующим образом:

(1)
$\frac{{\partial {\mathbf{V}}}}{{\partial t}} + ({\mathbf{V}}\nabla ){\mathbf{V}} = - \frac{1}{{{{\rho }_{0}}}}\nabla {\mathbf{P}} + \nu \Delta {\mathbf{V}},\quad \operatorname{div} {\mathbf{V}} = 0,$
где t – время, ρ0 – плотность и ν – кинематическая вязкость раствора.

Для рассчитанного поля скорости V решаются уравнения конвективного переноса в растворе для двух солей:

(2)
$\frac{{\partial {{M}_{i}}}}{{\partial t}} + {\mathbf{V}}\nabla {{M}_{i}} = {{D}_{i}}\Delta {{M}_{i}},$
где Mi = ρ0Ci – концентрации солей в растворе [кг/м3], ρ0 плотность раствора и Ci – относительные массы солей на 1 кг раствора; i соответствует номеру специи: 1 – KCSH, 2 – KNSH. Параметры гидродинамической модели задаются в соответствии с табл. 1.

Для расчета нормальной скорости роста кристалла из водно-солевого раствора использована формула для простого дислокационного источника [13]:

(3)
$R = \frac{{\beta {{C}_{e}}hkT\sigma _{0}^{2}}}{{19\alpha }}.$

По параметрам табл. 2 оценочная скорость кристаллизации составляет R = 5.66 × 10–9 м/с (или 0.489 мм/сут).

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

на участке втекания раствора (1) задается скорость струи и концентрации солей: Vy = Vjet, Ci = C0i, где С0i = ρ0exp(σ0)Cei. Пересыщение водно-солевого раствора рассчитывается как σi = = ln(Ci/Cei) [10], где Сi – текущие и Сei – равновесные концентрации солей. Взяв из табл. 2 начальное пересыщение σ0 = 0.09, вычисляем: C01 = = 53.68 кг/м3, C02 = 108.58 кг/м3;

на границе вытекания (5) градиент скорости и солевые потоки равны нулю;

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

на поверхности кристаллизации (4) задается соотношение баланса масс для каждой соли с учетом величины скорости кристаллизации, рассчитываемой по формуле (3):

(4)
${{\rho }_{{0i}}}{{D}_{i}}\frac{{\partial {{C}_{i}}}}{{\partial n}} = R({{C}_{{si}}}{{\rho }_{s}} - {{C}_{{ei}}}{{\rho }_{0}}).$
Здесь плотность смешанного кристалла ρs = = 2240 кг/м3. Величины Сsi задаются с учетом известного из эксперимента коэффициента распределения K для этих солей. Он связывает их концентрации в кристалле и растворе:
$K = \frac{{{{C}_{{s2}}}{{C}_{1}}}}{{{{C}_{{s1}}}{{C}_{2}}}} = 3.7,$
где Cs1 + Cs2 = 1. В результате получаем требуемые формулы для расчета Csi:
${{C}_{{s1}}} = \frac{{{{C}_{1}}}}{{{{C}_{1}} + K{{C}_{2}}}},\quad {{C}_{{s2}}} = \frac{{K{{C}_{2}}}}{{{{C}_{1}} + K{{C}_{2}}}},$
где Ci – концентрации солей в растворе на границе фаз, которые определяются в процессе итерационного счета. В данной работе расчеты ограничены их первым приближением: Ci = Cei.

Таким образом, формула (4) сопрягает гидродинамическую и термодинамическую модели.

Программные реализации моделей разработаны на основе пакета программ AnsysFluent [14], дополненного пользовательскими UDF-подпрограммами на языке Cи.

ОСОБЕННОСТИ ГИДРОДИНАМИКИ РАСТВОРА ДЛЯ ГЛАДКОЙ И ШЕРОХОВАТОЙ ПОВЕРХНОСТЕЙ КРИСТАЛЛИЗАЦИИ

Расчеты проведены для гладкой и треугольно-шероховатой поверхностей кристалла (согласно схемам на рис. 1).

Картина циркуляции раствора в кристаллизаторе с гладкой поверхностью кристаллизации показана на рис. 2а при скорости втекания Vjet = = 0.9 м/с, что соответствует скорости потока, создаваемой в реальных экспериментах по выращиванию кристаллов K2(Co,Ni)(SO4)2 · 6H2O. Втекающая струя сталкивается с твердой преградой в виде поверхности кристалла, в результате чего меняет осевое направление течения на радиальное. Ввиду большой скорости центробежной струи вблизи поверхности кристалла возникает интенсивная вихревая циркуляция при ее столкновении с боковой стенкой сосуда. Раствор интенсивно стекает вниз в центральной части кристаллизатора (вне области втекания), причем в остальной радиально-периферийной области возникает слабая замкнутая (запирающая течение) циркуляция.

Рис. 2.

Линии тока в растворе при Vjet = 0.9 м/с при гладкой (а) и шероховатой поверхности кристалла (б).

В случае треугольно-шероховатой поверхности кристалла наблюдается аналогичная структура циркуляции раствора в объеме кристаллизатора. Однако вблизи поверхности кристаллизации проявляется влияние шероховатости. Так, между треугольными выступами возникают слабые вихри, показанные на рис. 2б. Их наличие влияет на распределение скорости обтекания поверхности кристаллизации. Это можно видеть из сравнения радиальных профилей радиальной компоненты скорости Vr(r), показанных на рис. 3.

Рис. 3.

Профиль радиальной скорости вблизи поверхности кристалла y = 0.0395 м при Vjet = 0.9 м/с. Сплошная линия – при шероховатой, пунктир – при гладкой поверхности кристалла.

Заданная высота треугольных выступов δy = = 0.4 мм соответствует координате y = 0.0396 м. Между выступами циркулируют слабые вторичные вихри, а вне − на малом удалении (y = 0.0395 м) существует качественно такая же циркуляция раствора, что и для гладкой поверхности кристалла (рис. 2). В сравнении профили радиальной скорости для гладкой и треугольно-шероховатой поверхностей показаны на рис. 3 при y = 0.0395 м.

При гладкой поверхности кристалла Vr(r) резко возрастает за счет радиального поворота струи от нуля на оси до 0.6 м/с, затем резко падает до 0.2 м/с и остается таковой на большей части обтекаемого кристалла.

Однако наличие шероховатости существенно изменяет профиль скорости Vr(r): из-за дополнительных вихрей между треугольными выступами увеличивается величина и несколько изменяется профиль Vr(r).

ОСОБЕННОСТИ РАСПРЕДЕЛЕНИЯ ПЕРЕСЫЩЕНИЯ ВДОЛЬ ГЛАДКОЙ И ШЕРОХОВАТОЙ ПОВЕРХНОСТЕЙ КРИСТАЛЛИЗАЦИИ

Рассмотренные выше особенности гидродинамики влияют на распределения концентрации солей в растворе. В целом макрокартины изолиний концентрации солей KCSH и KNSH идентичны для гладкой и шероховатой поверхностей кристаллизации. Поэтому достаточно их проанализировать для одной из солей – KCSH.

Распределение изолиний концентрации KCSH на рис. 4 показывает, что за счет втекающей из трубки струи (рис. 1) в приосевой области существует максимальная концентрация соли, которая достигает поверхности кристалла и растекается по ней за счет радиального поворота струи. Осаждение (кристаллизация) соли на поверхности кристалла снижает ее концентрацию в остальном объеме кристаллизатора, из которого раствор вытекает через его донную поверхность (5) (рис. 1).

Рис. 4.

Картина изолиний концентрации соли KCSH кг/м3 при скорости струи Vjet = 0.9 м/с.

По данным концентраций солей рассчитано их пересыщение σi на фронте кристаллизации по формуле [3]: σi = (СiCei)/Cei.

В соответствии с данными, приведенными на рис. 4, в случае гладкой поверхности кристаллизации пересыщение в центре, на участке натекающей струи ∼9.8% и оно монотонно падает при удалении от оси до ∼8.2%. При этом различие обеих солей не значительное (рис. 5а).

Рис. 5.

Профили пересыщенности солями при Vjet = = 0.9 м/с (сплошная линия – KCSH, пунктир – KNSH): а – на гладкой поверхности кристалла, б – вблизи шероховатой поверхности кристалла при y = 0.0395 м (hole – углубление, hill – выступ на треугольно-шероховатой поверхности).

Однако в случае шероховатой поверхности кристаллизации проявляется осциллирующий характер σi в радиальном направлении (рис. 5б), вызванный треугольными возвышениями (hill) и впадинами (hole), поскольку их наличие обусловливает слабые вторичные вихри и осциллирующее радиальное изменение концентраций солей KCSH и KNSH.

ЗАКЛЮЧЕНИЕ

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

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

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

Работа выполнена при поддержке Российского научного фонда (грант № 15-12-00030).

Список литературы

  1. Dinakaran S., Verma S., Das S.J. et al. // Physica B. 2010. V. 405. № 18. P. 3919.

  2. Vartak B., Yeckel A., Derby J.J. // J. Cryst. Growth. 2005. V. 281. № 2–4. P. 391.

  3. Zhou C., Lin M., Hu Z. et al. // J. Cryst. Growth. 2016. V. 450. P. 103.

  4. Brailovskaya V.A., Zilberberg V.V., Feoktistova L.V. // J. Cryst. Growth. 2000. V. 210. P. 767.

  5. Robey H.F., Maynes D. // J. Cryst. Growth. 2003. V. 259. P. 388.

  6. Liiri M., Enqvist Y. // J. Cryst. Growth. 2006. V. 286. P. 413.

  7. Гребнев В.В., Ковалев С.И., Волошин А.Э. и др. // Кристаллография. 2017. Т. 62. № 6. С. 994.

  8. Voloshin A.E., Manomenova V.L., Rudneva E.B. et al. // J. Cryst. Growth. 2018. V. 500. P. 98.

  9. Воронцов Д.А., Волошин А.Э., Гребенев В.В. и др. // Кристаллография. 2017. Т. 62. № 6. С. 986.

  10. Rashkovich L.N. KDP-family single crystals. N.Y.; Bristol: Adam Hilger, 1991. 200 p.

  11. Ribeiro A.C.F., Lobo V.M.M., Natividade J.J.S. // J. Chem. Eng. Data. 2002. V. 47. P. 539.

  12. Ribeiro A.C.F., Verissimo L.V.M.M., Gomes J.C.S. et al. // Comptes Rendus Mecanique. 2013. V. 341. P. 417.

  13. Современная кристаллография. Т. 3. / Под ред. Вайнштейна Б.К. М.: Наука, 1980. 407 с.

  14. Ansys C.F.D. // Lisence of IPMech RAS. 659778-23-Aug-2011.

Дополнительные материалы отсутствуют.