Содержание к диссертации
Введение
Глава 1 Математическая постановка задачи 16
1.1 Физическая модель 16
1.1.1 Устройство алюминиевого электролизёра 16
1.1.2 Основные физические процессы в электролизной ванне 19
1.1.3 Рабочая гипотеза математической модели 20
1.2. Математическая модель 21
1.2.1 Уравнение объёмного баланса смеси 23
1.2.2 Уравнение неразрывности 24
1.2.3 Уравнение движения 25
1.2.4 Моделирование давления 27
1.2.5 Математическое моделирование электродинамических процессов 31
1.2.6 Математическое моделирование электрохимических процессов 35
1.2.7 Математическое моделирование тепловых процессов 38
1.2.8 Полная математическая модель 39
Глава 2 Численный метод решения 45
2.1 Разделение на этапы по физическим процессам 45
2.2 Численный метод решения задачи 50
2.3 Численное исследование разностного метода . 66
2.3.1 Устойчивость разностного метода. 66
2.3.2 Оценка точности метода. 67
Глава 3 Анализ численных экспериментов 71
3.1 Тестовые и верификационные расчёты 71
3.2 Сравнение многоанодного электролизёра и электролизёра Содерберга 89
3.4 Математическое моделирование анодного эффекта 116
Заключение 130
Список литературы 131
- Основные физические процессы в электролизной ванне
- Математическое моделирование электродинамических процессов
- Численное исследование разностного метода
- Сравнение многоанодного электролизёра и электролизёра Содерберга
Основные физические процессы в электролизной ванне
Алюминиевый электролизёр представляет собой сложный нелинейный объект, который включает в себя множество процессов, влияющих друг на друга. Основные процессы можно разбить на несколько групп:
Перечисленные процессы оказывают сильное влияние друг на друга, однако имеют разное характерное время протекания. В связи с этим многие авторы используют усреднение по времени, проводя моделирования без учёта взаимовлияния перечисленных выше процессов. В представленной диссертации усреднений не проводится. Это позволяет провести математическое моделирования всех процессов во взаимосвязи. Остановимся подробнее на роли каждого из процессов при получении первичного алюминия.
Гидродинамические процессы включают в себя изменение конфигурации компонентов смеси внутри электролизной ванны и изменение скоростей отдельных фаз смеси, при этом главной вынуждающей силой является сила Лоренца.
Электромагнитные процессы описывают изменение плотности электрического тока, напряжённости электрического и магнитного полей, определяют характер протекания электролиза алюминия в ванне и обуславливают развитие МГД-нестабильности. Поскольку МГД-нестабильность, развивается довольно быстро, в течение долей секунд, то их можно отнести к быстротекущим процессам.
Химические процессы описывают восстановление первичного алюминия из криолит-глинозёмного расплава, а также окисление анода с выделением газа и протекают очень медленно. Еще одной их отличительной чертой, является протекание в тонком слое под подошвами анодов и на поверхности расплавленного ме талла. Накапливаемые в приэлектродных областях заряды вносят вклад в формирование электромагнитных полей. Однако область, где заряд играет существенную роль, настолько мала по сравнению с объёмом ванны, что в целом рассматриваемая смесь берётся электронейтральной. Химические процессы описывают также дрейф ионов в электролите и являются главной причиной развития анодного эффекта.
Тепловые процессы также оказывают существенное влияние на процесс электролиза алюминия. Повышение температуры выше критической может привести к тому, что плотность алюминия станет меньше плотности электролита, вызывая при этом перемешивание всей смеси. В случае, когда температура становится меньше определенного значения, происходит кристаллизации жидких компонент расплава с последующим оседанием их на боковых стенках ванны и в падине. При повышении до величины критической температуры начинается плавление боковых стенок ванны. В обоих случаях происходит изменение геометрии рабочего пространства ванны, которое существенно влияет на МГД-стабильность процесса электролиза алюминия. Таким образом, математическое моделирование изменения формы рабочего пространства ванны должно учитывать изменение температуры вблизи гарнисажа и настыли.
Рассматриваемая в электролизной ванне смесь состоит из трёх фаз: жидкого алюминия, электролита и газа, таким образом, в основе предлагаемой математической модели находится многофазный и многоскоростной подход [23, c. 14-39].
На интервалах времени, для которых была разработана математическая модель, делаются следующие предположения:
В сделанных предположениях математическое моделирование проводилось как для многоанодных электролизёров с обожжёнными анодами, так и для од-ноанодных электролизёров (электролизёров Содерберга). Вычислительные эксперименты проводились как для упрощённой конфигурации алюминиевого электролизёра, которая представляет собой прямоугольный параллелепипед без внутреннего гарнисажа и настыли, так и для конфигураций с реальной геометрией рабочего пространства промышленного электролизёра. В обоих случаях исследовалось влияние формы рабочего пространства на МГД-стабильность ванны.
Введем в рабочем объёме алюминиевого электролизёра декартову систему координат. Пусть ось ОХ будет направлена вдоль короткой стороны ванны, ось OY - вдоль длинной стороны ванны, а ось OZ направлена вертикально вверх (рисунок 1.3). Начало координат совмещено с углом ванны, к которому ближе первый анод, на уровне подошвы слоя металла. Введенная система координат изображена на рисунке 1.3.
Для описания гидродинамических процессов в электролизной ванне используется система уравнений Навье-Стокса в переменных Эйлера, записанная в средах жидкого металла, электролита и газовой фазы. При моделировании электромагнитных полей и плотности электрического тока за основу бралась система уравнения Максвелла, записанная в системе единиц измерения СИ. Химическая кинетика электролиза алюминия описывается уравнениями Нернста-Планка для основных ионов участвующих в процессе.
В гетерогенной смеси каждая фаза занимает лишь часть элементарного объёма, в связи с этим вводятся объёмные доли а (т = 1,2,3) для каждой фазы, кото рые характеризуют доли, занимаемые каждой фазой смеси, в элементарном объёме [23, с. 14-39]. При этом естественно считать, что в элементарном объёме выполняется уравнение объёмного баланса смеси
Математическое моделирование электродинамических процессов
Система (1.4), (1.11), (1.14), (1.30), (1.31) и (1.34) решается методом расщепления на этапы по физическим процессам [31]. Ниже перечислены этапы решения задачи на одном временном слое. Начальные условия, которые требуются для решения соответствующих уравнений, считаются известными. корка Этап 11 включает в себя также изменение формы рабочего пространства ванны за счёт кристаллизации вещества и плавления гарнисажа или настыли.
На блок схеме представлен алгоритм расщепления по физическим процессам [31][36, c. 246-256].
Сетка Qc состоит из внутренних узлов расчётной области, граничных узлов и внешних узлов. Расчёт всех скалярных величин кроме температуры происходит во внутренних узлах сетки Qc, температура также рассчитывается в граничных узлах, в которых значения потоковых величин скорости и плотности тока принимаются равными нулю. имеет второй порядок аппроксимации по пространству и первый порядок аппроксимации по времени, т.е. iy = 0(Ax2,Ay2,Az2,t).
Потоки на границах аппроксимируются аналогично основному уравнению с первым порядком точности по пространству:
Порядок аппроксимации = 0(Ax2, Ay2, Az2, ґ), так как разностные производные по пространству, аппроксимирующие дивергенцию, определенны на разнесённой сетке, имеют второй порядок.
Разностный аналог уравнения (2.4) имеет второй порядок аппроксимации по пространству. Для вычисления потенциала электрического поля р используется ме тод верхних релаксаций с параметром со є [1.8,1.95]. Для ускорения сходимости в начале итерационного процесса берётся со = 1.95. В дальнейшем итерационный параметр уменьшается до со = 1.8 для получения необходимой точности. Итерационный процесс останавливается когда разность итераций в норме С по модулю меньше 10-6.
Плотность электрического тока находится явно из уравнений (индекс временного слоя опущен, все величины берутся с «-го слоя по времени), которые аппроксимируют (2.5) с i// = 0(Ax2,Ay2,Az2) [45]: Для уравнений (2.6) строятся разностные аналоги по схеме «крест» (индекс временного слоя опущен, все величины берутся с n-го слоя по времени) со вторым порядком аппроксимации по пространству [44,45]:
С учётом введённой сетки граничные условия для всех компонент вектора потенциал магнитного поля имеют следующие разностные аналоги (индекс компоненты опущен):
Разностный аналог уравнения (2.6) имеет второй порядок аппроксимации по пространству. Для вычисления потенциала напряжённости магнитного поля А используется метод верхних релаксаций с параметром со є [1.8,1.95]. Для ускорения сходимости в начале итерационного процесса берётся со = 1.95. В дальнейшем итерационный параметр уменьшается до с? = 1.8 для получения необходимой точности. Итерационный процесс останавливается когда разность итераций в норме С по модулю меньше 10-6.
Компоненты вектора магнитного поля находится явно из уравнений (индекс временного слоя опущен, все величины берутся с «-го слоя по времени), которые аппроксимируют (2.7) со вторым порядком [45]: Я =( --1 ),-(Л -1Ь = (4 -1)z-(4 b Hz,k={A _1jk)x-{A _1k) Учитывая, что величины магнитного поля и векторного потенциала определены в различных узлах сеток, разностная аппроксимация имеет второй порядок.
При численном решении уравнения движения основная сложность связана с аппроксимацией нелинейного конвективного члена уравнения. Воспользуемся методом Мак-Кормака второго порядка [34, c. 199-201]. Рассмотрим применение этого метода на примере решения модельной одномерной задачи: Давление/? и внешние силы/пересчитывается отдельно как для шага предиктора, так и для шага корректора по известным формулам.
В дальнейшем будем пользоваться следующими их покомпонентными обозначениями покомпонентной аппроксимации на п-ом слое по времени шага предиктора или корректора для расчёта вектора вязкого трения Р, коэффициента трения К и сил Лоренца F [45]: Pm
Разностный аналог уравнения (2.8) имеет второй порядок аппроксимации по пространству. Для вычисления давления р на шаге предиктора используется метод верхних релаксаций с параметром со є [1.8,1.95]. Для ускорения сходимости в начале итерационного процесса берётся ю = 1.95. В дальнейшем итерационный параметр уменьшается до й? = 1.8 для получения необходимой точности. Итерационный процесс останавливается когда разность итераций в норме С по модулю меньше 10-6.
Численное исследование разностного метода
Данный тест подтверждает корректность многофазного подхода математической модели. Проверяется корректность описания вязкого трения между жидким металлом и электролитом. В данной модельной задаче токи, напряжённость электромагнитного поля и электромагнитные силы заданы равными нулю. На рисунке 3.16 изображены скорости в электролите в начальный момент времени, скорости в металле в начальный момент времени равны нулю. Поверхность раздела фаз алюминий-электролит берётся горизонтальной, без каких-либо возмущений.
Поскольку уравнение движения учитывает передачу импульса между фазами в смеси, скорости заданные в начальный момент времени в электролите должны передать импульс жидкому алюминию и повлиять на его движение. На рисунках 3.17 и 3.18 показано поле скоростей в среднем слое алюминия, интенсивность которого увеличивается со временем.
Данный вычислительный эксперимент показывает, что слагаемое, отвечающее за передачу импульса между фазами, описано корректно.
В работах [15,39] были представлены нестационарные математические модели для расчёта напряжённости магнитного поля, которые достаточно точно описывают его эволюцию. Однако существенным недостатком таких моделей яв ляется их зависимость от начальных данных, которые в случае алюминиевого электролизёра несут в себе большую погрешность, поэтому сила Лоренца в таких моделях качественно некорректна. Замеры промышленного электролизёра имеют 11 узлов по оси OX и 31 узел по оси OY, что позволяет увидеть только общую структуру вектора напряжённости магнитного поля. Согласно замерам напряжённость магнитного поля в среде электролита имеет один большой вихрь (рисунок 3.19).
Из общих физических соображений, для многоанодного электролизёра в среде электролита следует ожидать небольшие вихри магнитного поля непосредственно под каждым анодом, которые будут формировать один большой обводящий вихрь (рисунок 3.20). Рисунок 3.19. Напряжённость магнитного поля многоанодного электролизёра в среднем слое электролита. Экспериментальные данные
Поскольку в моделях [15,39] использовались эволюционные модели, погрешность начальной напряжённости магнитного поля сохранялась в процессе расчёта. Силы Лоренца, рассчитанные на основе такого распределения напряжённости магнит-85 ного поля, качественно отличались от данных, полученных по представленной модели. При разряжении сетки, напряжённость электромагнитного поля, полученная из расчёта, будет соответствовать данным полученным при замерах. Тест 5 Для проверки работы математической модели в целом проводится численное моделирование развития МГД-нестабильности при моделировании замены крайней пары (11-го и 22-го) анодов (рисунок 3.21).
Согласно наблюдениям технологов, замена этой пары анодов может привести к развитию нежелательных колебаний на границе раздела сред алюминий-электролит. В работах [15,26] проводилось математическое моделирование замены выгоревшей пары крайних анодов, которое показало, что поверхность алюминия поднимается выше критических значений, что говорит о развитии МГД-нестабильности процесса электролиза. На рисунке 3.22 приведён график максимальной высоты поверхности уровня алюминия после извлечения анодов полученный авторами в работе [26]. Из графика видно, что в момент времени t =0.25 с после отключения подачи тока на 11-й и 22-й аноды, высота достигает 30 см.
Как уже было отмечено, напряжённость магнитного поля, найденная с использованием стационарной модели, описанной в пункте 1.2.5 первой главы, качественно отличается от поля, рассчитанного в работе [26]. В связи с этим, поверхность раздела фаз также отличается от раздела фаз алюминий-электролит, однако общий профиль поднятия поверхности алюминия сохраняется (рисунок 3.23).
Динамика высоты алюминия после извлечения 11 и 22 анодов. Шаг по времени – 0.015 с, макс. высота гр. Раздела – 0.311 см График межполюсного расстояния (МПР) (рисунок 3.25), полученный по описанной в главе 1 математической модели, также имеет максимум в момент времени t = 0.26c с максимальной амплитудой 0,341 м. После достижения максимума в обоих случаях идёт затухание колебаний поверхности жидкого металла.
Сравнение численных экспериментов по моделированию технологического процесса замены крайней пары анодов показало, что представленную в работе математическую модель и численный алгоритм ее решения можно использовать для моделирования процесса электролиза алюминия.
Современная промышленность использует при добыче первичного электролиза алюминия два разных типа электролизёров: многоанодный электролизёр с 22-мя и более обожженными анодами, погруженными в электролизную ванну в два ряда; электролизёр Содерберга с одним самообжигающимся анодом.
В работе [6] обсуждаются преимущества и недостатки обоих типов электролизёров с экологической точки зрения, делается вывод, что многоанодные электролизёры производят меньше парниковых газов, чем одноанодные. В данном исследовании проводится сравнение электролизёров с точки зрения их МГД-устойчивости [43].
Рассмотрим подробнее электролизёр Содерберга. Поскольку для данного типа электролизёра используется самообжигающийся анод, необходимо периодически насыпать в него новую анодную массу для поддержания процесса электролиза. Анодная масса в процессе обжига становится неоднородной с точки зрения проводимости электрического тока, что может повлиять на МГД-стабильность процесса.
Сравнение многоанодного электролизёра и электролизёра Содерберга
Подина – область ванны, где потенциал считается равным нулю, занимает всю горизонтальную плоскость дна ванны. Исходя из описанной конфигурации анодов и катода, распределение плотности тока в срезах YZ и XZ принимает вид рисунка 3.44. В представленных графиках, области наибольшей по модулю плотности тока находятся в углах анодов и в углах подины, данные распределения совпадают с исследованиями, представленными в работах по моделированию плотности тока других авторов [26,40]. а) срез плоскостью YZ, ближайший б) срез плоскостью XZ к началу координат ряд анодов пятая пара анодов Области повышенной температуры (рисунок 3.45) соответствуют областям с наибольшей по модулю плотностью тока, что говорит о прямой связи между плотностью тока и температурой в электролизной ванне. а) срез плоскостью YZ, ближайший б) срез плоскостью XZ к началу координат ряд анодов пятая пара анодов При одинаковом распределении потенциала по анодам (примерно 4 В), плотность тока и распределение температуры качественно остаются неизменными.
Численный эксперимент при уменьшенном потенциале на пятой паре анодов. В данном вычислительном эксперименте потенциал на пятой паре анодов был уменьшен на 15%. Это вызвало перераспределение плотности тока в электролизной ванне и изменение температуры. Результирующее рабочее пространство электролизной ванны (рисунок 3.46) имеет направленную внутрь деформацию, которая вызвана нарастанием гарнисажа и натыли, изменение составляет порядка 12 – 16см.
На рисунках 3.47 и 3.48 представлено результирующее рабочее пространство ванны в срезах плоскостями YZ и XZ соответственно. Рисунок 3.48 более подробно иллюстрирует изменение гарнисажа и настыли.
При уменьшении потенциала электрического поля на пятой паре анодов (рисунок 3.49) происходит перераспределение плотности тока, протекающего в рабочем пространстве ванны. В связи с этим уменьшается выделение тепла (рисунок 3.50), связанного с законом Джоуля-Ленца, и происходит остывание области ванны по пятой паре анодов с последующим нарастанием гарнисажа.
На рисунке 3.51 представлен график динамики МПР и колебания границы зоны обратного окисления. Красной линией обозначено максимальное поднятие поверхности металла в процессе электролиза алюминия, синей – максимальная глубина зоны обратного окисления.
Рисунок 3.51. Динамика МПР 3. Численный эксперимент при уменьшенном потенциале на крайней паре анодов.
В данном вычислительном эксперименте потенциал на крайней (первой) паре анодов был уменьшен так же на 15%. В данном случае изменение формы ванны происходит с одного края ванны симметрично относительно оси OY. Результирующая форма ванны представлена на рисунке 3.52.
На рисунках 3.53 и 3.54 представлено результирующее рабочее пространство ванны в срезах плоскостями YZ и XZ соответственно. Рисунок 3.54 более подробно иллюстрирует изменение ФРП у левого края ванны, где настыль заходит за проекцию анодов на дно ванны.
При уменьшении потенциала электрического поля на крайней паре анодов (рисунок 3.55), так же как и в первом численном эксперименте, происходит перераспределение плотности тока, протекающего в рабочем пространстве ванны. Так как пара анодов является крайней, перераспределение тепла отличается от перераспределения тепла для случая с пятой парой анодов, из-за этого рабочая область ванны остывает сильнее (рисунок 3.56), что делает нарастание гарнисажа более интенсивным (в среднем на 6 см). а) срез плоскостью YZ, ближайший б) срез плоскостью XZ к началу координат ряд анодов первая пара анодов Рисунок 3.55. Результирующее распределение плотности тока, t = 20 c Ill а) срез плоскостью YZ, ближайший б) срез плоскостью XZ к началу координат ряд анодов первая пара анодов Рисунок 3.56. Результирующее распределение температуры, t = 20 c На рисунке 3.57 представлен график динамики МПР и колебания границы зоны обратного окисления. Красной линией обозначено максимальное поднятие поверхности металла в процессе электролиза алюминия, синей – максимальная глубина зоны обратного окисления.
Численный эксперимент при увеличенном потенциале на пятой паре анодов. В данном численном эксперименте потенциал на пятой паре анодов был увеличен на 15% по сравнению с исходным значением. Данное увеличение потенциала приводит к увеличению плотности протекающего по ванне тока и как следствие к увеличению выделяемого тепла. Получившаяся форма рабочего пространства представлена на рисунке 3.58.
На рисунок 3.59 и 3.60 представлено результирующее рабочее пространство ванны в срезах плоскостями YZ и XZ соответственно. Рисунок 3.60 более подробно иллюстрирует увеличение гарнисажа, которое является симметричным относительно оси OY. Увеличение ФРП ванны составляет порядка 10 – 13 см. относительно начальной формы.