Содержание к диссертации
Введение
Глава 1. Алгоритмы решения двух электрофизических задач 23
1.1 Методы решения задачи индукционного каротажа 24
1.1.1 Постановка задачи 24
1.1.2 Конечнообъемная аппроксимация 26
1.1.3 Сходимость сеточных решений 29
1.1.4 Итерационные алгоритмы 31
1.1.5 Численные исследования 37
1.2 Численное решение трехмерной задачи магнитостатики в области с ферромагнетиками 40
1.2.1 Постановка задачи 40
1.2.2 Конечнообъемная аппроксимация 43
1.2.3 Численная реализация формулы Био-Савара-Лапласа 52
1.2.4 Аппроксимация табличной функции магнитной индукции эрмитовым сплайном 53
1.2.5 Итерационные алгоритмы решения нелинейной системы 54
1.2.6 Численные исследования 56
Глава 2. Численное решение системы трехмерных уравнений электростатических и тепловых задач 63
2.1 Постановка задачи 64
2.2 Конечнообъемная аппроксимация 65
2.3 Итерационные алгоритмы 68
2.4 Сходимость итерационного процесса 71
2.5 Примеры вычислительных экспериментов 73
Глава 3. Комплексное моделирование электромагнитных и тепловых полей алюминиевого электролизера 84
3.1 Физическая модель электролизера 85
3.2 Решение задачи о совместном распределении электрических полей в электролизере и ошиновке 90
3.3 Вычислительная схема комплексного моделирования 93
3.4 Результаты численных расчетов 95
Приложение. Технологические аспекты решения трехмерных краевых задач 104
- Численное решение трехмерной задачи магнитостатики в области с ферромагнетиками
- Аппроксимация табличной функции магнитной индукции эрмитовым сплайном
- Конечнообъемная аппроксимация
- Решение задачи о совместном распределении электрических полей в электролизере и ошиновке
Численное решение трехмерной задачи магнитостатики в области с ферромагнетиками
В последние 5-10 лет резко возросла производительность ЭВМ. Известно, что при необходимости рядовой пользователь может обзавестись компьютером с частотой до 3400Мгц, а оперативной памятью до ЗГб, размеры жестких дисков для хранения данных и результатов расчета достигают порядка 200Гб. Благодаря этому стало возможным проводить более серьезные численные расчеты на персональных компьютерах: не ограничиваясь памятью и временными рамками, подтверждать теоретические предпосылки и теоремы о сходимости численного решения при сгущении расчетной сетки, при наличии аналитического решения проверять сходимость к таковому, а при отсутствии аналитики говорить о гарантированной точности численных экспериментов. Появились широкие возможности по подтверждению и (или) опровержению теоретических результатов о сходимости итерационных методов, сравнению этих методов между собой. Наконец, повышение производительности ЭВМ играет большую роль в уходе от многих одномерных и двумерных задач математической физики, переходе на качественно новый уровень решения трехмерных задач, в комплексном моделировании взаимосвязанных физических полей,
Несмотря на большой прогресс в развитии вычислительной техники, остается открытым вопрос о создании комплексов программ, позволяющих эффективно моделировать различные физические процессы, возникают все новые и новые требования к таким программам. Цель диссертации: используя эффективные методы решения краевой задачи Пуассона создать и программно реализовать алгоритмы, структуры данных для решения задачи индукционного каротажа, комплексного решения нелинейных трехмерных задач магнитостатики, электростатики и теплопроводности; применить разработанные методы и программы для комплексных расчетов физических полей в алюминиевом электролизере. В процессе работы над диссертацией создан также алгоритлі совместно му го решения дифференциальной задачи электропроводности в трехмерной расчетной области и уравнений Кирхгофа для электрических цепей, соеди ненных непосредственно с этой областью и влияющих на общее распреде ление электрического тока. Разработаны алгоритмы обработки численных результатов с их представлением в виде изолиний, силовых линий, вектор N ных полей и т.д. v - Индукционный каротаж (ИК) является одним из способов исследования верхнего слоя земной коры, толщиной до нескольких километров, основан ный на использовании переменного тока. Существует множество геофи зических методов изучения свойств земной коры [3], [20], [64], [6 5], [66], [70]. Часть из них основана на электрометрии скважин, другие являются ядер ными, акустическими, сейсмическими методами исследования. ( Сущность индукционного каротажа состоит в измерении и анализе вто ричного индукционного магнитного поля, созданного в горной породе под действием первичного переменного тока, а также в решении обратной задачи ИК. Обратная задача заключается в идентификации геометрических и материальных параметров модели по измеренным магнитным полям. Задача ИК широко используется в нефтеразведке [22]. Для решения обратной задачи [25] требуются многочисленные расчеты прямых задач ИК, т.е. моделирование электромагнитного поля по заданным свойствам среды. Впервые приближенная теория индукционного каротажа была предложена Г.Г.Доллєм [74]. Теория ИК, методы решения прямой задачи, задачи ИК с различным числом зондов подробно описаны в работах А.А.Кауфмана [33],[34], М.И.Эпова, Л.А.Табаровского [3],[37],[58],[69] и других зарубежных авторов [92], [93]. Если в горизонтальном и вертикальном зондировании можно ограни читься одномерным [36] и двумерным математическим моделированием [37], то в наклонном зондировании перспективным является трехмерное мо делирование [72],[75],[88]. Предложено множество различных методов для решения задач электромагнитной индукции, например [58],[92],[93]. Одним из методов для решения трехмерных задач геоэлектрики является метод Щ вторичных источников, предложенным О.В.Тозони [60],[61],[62]. Этот ме тод сводится к методу интегральных уравнений [4],[16],[71], который в свою fr очередь основывается на понижении размерности задачи. Однако, матрица, полученная в результате аппроксимации этим методом, является плотной, тем самым необходимо больше места для хранения ее элементов. К примеру, в трехмерном случае, при использования методов конечных элементов или конечных объемов, число неизвестных 0(/i 3), число ненулевых элементов матрицы 0(/і_3). При использовании интегральных методов число неизвестных сокращается до 0(/г 2), однако число ненулевых элементов в матрице возрастает до 0(h {). В результате, из-за большой ресурсоемкое (по памяти, скорости расчетов) этот метод редко используется. В последнее время для решения трехмерных задач ИК широко используются методы конечных элементов (МКЭ), а также векторные МКЭ, предложенные в работах Э.П.Шуршюй [68(,(89] и Ю.Г.Соловейчика [49],[50,[51,[52]. Большое внимание при моделировании магнитных полей уделяется дву ( мерным задачам [83]. Эти задачи возникают в области с осесимметричными свойствами среды, а также кольцеобразными токами (например, катушки), возбуждающими магнитное поле. Методы расчета магнитного поля в ускорителях заряженных частиц с использованием треугольных сеток описаны в работах С.Б.Ворожцова [15],[13],[12],[14]. В осесимметричных моделях задача о нахождении напряженности Н ноля сводится к поиску скалярного потенциала р [21]. В общем случае необходимо рассматривать трехмерную математическую модель. В трехмерном случае, помимо постановки задач со скалярным потенциалом, широко применяется смешанный метод, в котором искомыми функциями являются потенциал tp и поток j [86] ,[85], [84]. При невозможности использования скалярного потенциала, для вектора магнитной индукции В задают векторный потенциал X If = rotA [91]. В однородной среде, в которой всюду магнитная проницаемость \i — const, распределение магнитного поля описывается законом Био-Савара Лапласа. Области, в которых распространяются магнитные поля, по ти пу \i разделяются на парамагнетики pi 1, диамагнетики р, 1 и фер ромагнетики с нелинейной зависимостью функции магнитной проницае мости д = fi{H). Функция //(Я) может задаваться как аналитически, Щь так и как сплайн табличной функции, построенной на основании опытных данных [28].
Аппроксимация табличной функции магнитной индукции эрмитовым сплайном
Примеры 2 и 3 продемонстрировали второй порядок сходимости к аналитическому решению, что соответствует теоретическим предпосылкам.
Так как задача (2.1.1)-(2.1.4) актуальна при моделировании алюминиевых электролизеров [39], то геометрию и физические параметры дальнейших примеров возьмем соответствующие. Функции электропроводности материалов изображены на рис.2.2. Для следующих двух примеров расчетные области показаны на рис.2.3
Пример 4. Рассматривается область П, состоящая из 3-х параллелепипедов (рис.2.3.а) - одного внешнего 0,\ и двух внутренних 1 и Пз-Внешний Qj = [0,0.8] х [0,0.4] х [0,0.4] состоит из графита с коэффициентом теплопроводности А = 20. Два внутренних параллелепипеда П2 = [0.2,0.3] х [0.15,0.25] х [0,0.4] и П3 = [0.5,0.6] х [0.15,0.25] х [0,0.4] состоят из стали с коэффициентом теплопроводности А = 40.
На верхней грани задаем условие втекания тока силой 16 кА, а на торцевых гранях внутренних параллелепипедов — условие вытекания тока силой 4 кА с каждой грани. На остальных границах с внешней средой условие непротекания тока. Для тепла на нижней грани задаем условие постоянной температуры Т — ЗООК, а на остальных внешних гранях — условие не протекания теплового потока — 0. Расчеты проводятся на трех равномерных сетках (табл.2.5). В качестве следующего примера берем модель электролизера, рассчитанного пакетом программ Пример 5. Рассматривается область, состоящая из 14-ти параллелепипедов (рис.2.3.Ь). Из-за сложности геометрии ее точные размеры приводить не будем. Четыре вертикальных и три горизонтальных штыря состоят из стали с коэффициентом теплопроводности Л = 40. Нижний горизонтальной слой -графитосодержащая подина (Л = 20), 2-ой горизонтальный слой - алюминий (Л = 110), 3-ий слой - электролит (Л = 2). Четыре блока, находящиеся сверху электролита, состоят из графита (Л = 20). Функции электропроводности задаются в соответствии с рисунком (2.2). Сверху четырех вертикальных стержней задаем условие втекания тока силой 6 кА в каждую грань, с торцов трех горизонтальных стержней условие вытекания тока силой 4 кА с каждой грани, а на остальных внешних гранях условие непротекания тока. Для тепла задаем условие постоянной температуры на всей границе области, равной 300К. В таблицах 2.5 и 2.6 продемонстрированы результаты выполнения примеров 4 и 5. Из табличных значений ( и 8 видна сходимость приближенного решения с порядком h. Во всех примерах число внешних итераций nout практически не увеличивается с дроблением сетки. Можно предположить, что это связано с тем, что внешний итерационный процесс связывает два дифференциальных уравнения из (2.1.1). Таким образом, при сгущении расчетной сетки, число внешних итераций приближается к некоторому числу, зависящему от нелинейности коэффициента а. Легко заметить, что если a = const, то W = 1. По значениям из таблиц заметно, что число nout больше в 2,5-ой и 2.6-ой таблицах, чем в 2.3-ой и 2.4-ой. Это подтверждает наше предположение о числе внешних итераций, т.к. в 4-ом и 5-ом примерах используются среды с различными по свойствам нелинейности коэффициентами а. Для примеров 4,5 на рисунках 2.4 - 2.5 приведены изолинии скалярного поля электрического потенциала, а также изотермы теплового поля. Помимо численных результатов, важным критерием оценки правильности расчетов является анализ графического представления данных. По изолиниям электростатического поля рисунка 2.4 можно судить о том, что электрический ток втекает сверху, что соответствует адекватному учету граничных условий из примера 4. По изотермам можно судить о распределении температуры в расчетной области. Максимальная температура достигается на верхней грани. Это предсказуемо, так как на всех гранях, кроме нижней, задано условие ненротекания тепла. Тепло, выделяясь внутри расчетной области,выходит снизу, где задано условие охлаждения (условие Дирихле). Рис. 2.4. Изолинии электрического потенциала и изотермы теплового поля для примера 4
По характеру изолиний и изотерм на рисунке 2.5 также можно говорить о правильности задания граничных условий в примере 5. Изолинии перпендикулярны границам - ток не вытекает, а для поля температур на границе области выполняется заданное условие постоянной температуры Т = 300К. Максимальная температура достигается в верхней части вертикальных штырей. Это объясняется тем, что сверху штырей втекают токи, а затем эти токи распределяются по всей области. Плотность токов на верхней грани штырей максимальна, что влияет на интенсивность выделение тепла. Также отметим, что важным критерием правильности расчетов является симметрия на приведенных рисунках.
Конечнообъемная аппроксимация
Как отсюда видно, характеристики неравномерности токораспределения но штырям для данных ВАМИ и НКАЗ совпадают с очень высокой точностью, а токорасіїределения для блюмсов отличаются значительно.
Отметим, что данные расчеты проводились на кусочно-равномерной сетке с числом узлов около 500000. Контрольные вычисления для проверки достаточной точности моделирования осуществлялись на более густых сетках с количеством узлов около 4000000 и 9100000. Приведенные результаты получены путем решения системы (2.3.9) с общим количеством итераций 5, после которых последовательные приближения практически не менялись (с точностью до 5 десятичных знаков).
Более детальные характеристики вычислительного процесса при решении комплексной задачи приведены в табл.3.9. Здесь для трех различных сеток даны количества "внутренних" итераций (итерационная точность бралась є = 10 6) при последовательных расчетах электрического потенциала р и температуры Т па разных внешних нелинейных итерациях. Ошибки численных решений оценивались по асимптотическому решению, оцененному с помощью экстраполяции Ричардсона по результатам расчетов на трех указанных сетках.
Расчеты проводились на компьютере Pentium-IV с частотой процессора 2.4GHz и 512Mb оперативной памяти. Помимо этих расчетов, был проведен ряд других экспериментов по моделированию электротеплового баланса в алюминиевом электролизере. Среди них можно выделить следующие: моделирование нагрева электролизера при изменения межполюсного расстояния от 45мм до 65мм; моделирование динамики распределения теплового поля во время сгорания анода с течением времени; моделирование изменения тока по анодным штырям и теплового поля при перестановки штырей с нижнего уровня на верхний; моделирование изменения распределения тепла при уменьшения теплопроводности подины; моделирование влияния корки сверху электролита, с учетом изменения коэффициента теплоотдачи свободной поверхности электролита; моделирование перегрева стенок электролизера при увеличении неравномерности токов, втекающих в анодные стержни; моделирование отключения стояка путем распределения втекающего тока по трем из четырех стояков. На рисунке (3.8) представлены гистограммы распределения рассчитанных токов по анодным штырям. Продемонстрирован цикл сгорания анода: начальное расположение штырей, а также поочередная перестановка штырей с нижнего уровня на верхний. По гистограммам нетрудно заметить, что максимальный ток протекает по штырям первого уровня, а минимальный — 4-го. Это объясняется тем, что расстояние между подошвой анода и штырями растет с увеличением номера уровня, при этом электропроводность анодных штырей значительно больше электропроводности самого анода. На Рис. 3.9 показана динамика распределения теплового поля в электролизере при уменьшении коэффициента теплоотдачи поверхности электролита (модель образовавшейся теплоизоляционной корки). Значительная часть выделяемого в электролизере тепла выходит через эту поверхность. Поэтому при уменьшении коэффициента теплоотдачи значительно возрастает температура внутри ванны и рабочего пространства электролизера. Моделирование некачественного соединения анодных штырей с ошинов кой производится путем задании неравномерности токов в анодных шты рях. Следствием такой неравномерности токов является перегрев нижней части анода и части подины (рис. 3.10), расположенных под штырями с У 1 большей плотностью тока. Также немаловажным является перегрев лево го борта ванны, что может привести к аварийной ситуации.
На основе рассчитанного электротеплового баланса проводилось моделирование распределения магнитного поля. В качестве токов, возбуждающих магнитное поле, задавалось рассчитанное распределение токов по стоякам, анодной и катодной ошиновке. Расчеты проводились как с учетом ферромагнитного корпуса (см.1.2), так и без него. Для расчета магнитного поля в областях с ферромагнетиками используется алгоритм, описанный в 1.2. Характеристики вычислительного процесса представлены в примере 4, п.1.2.6. На Рис. 3.11—3.12 показано расчитанное распределение напряженности Н магнитного поля на поверхности алюминия.
В ряде случаев необходимо визуально определить границу фазового перехода вещества их одного состояния в другое (затвердевание анодной массы, электролита). На Рис. (3.13) показана изопо верхи ость для некоторого абстрактного значения температуры Т = 600. По рисунку можно сделать некоторые выводы о нагреве различных конструктивных элементов электролизера выше 600 градусов. К примеру, в нижней части анода температура в основном выше 600 градусов, а в бортовой футеровке, из-за своей малой теплопроводности, температура значительно меньше.
Решение задачи о совместном распределении электрических полей в электролизере и ошиновке
В третьих, в разрабатывающемся пакете предусмотрена возможность задания пересечения двух или, в общем случае, нескольких параллелепипедов. Для этого разработан алгоритм заполнения массива med, учитывающий новое свойство параллелепипедов — lay. Каждому параллелепипеду cub пользователь присваивает номер слоя lay(cub), которому он принадлежит. Слои определяют последовательность, в которой рассматриваются параллелепипеды, при заполнении массива med. Сначала "перебираем" первый (нижний) слой, затем второй и т.д. Заканчиваем последним (верхним). Для параллелепипедов из одного слоя условия не пересечения, а также порядок перебора остаются такими же, как и прежде. Для двух параллелепипедов из разных слоев ограничений на взаимное расположение не существует. Однако, если они пересекаются в пространстве, то в процессе формирования массива med параллелепипед из верхнего слоя "затрет" нижний в местах их пересечения. Например, для задания геометрии, изображенной на рис.2, достаточно определить не семь, а четыре параллелепипеда. При этом 2-ой, 3-ий и 4-ый присвоить слою JVH, а 1-ый слою №2. Остается заметить, что прежнее задание геометрии аналогично принадлежности всех параллелепипедов одному слою.
Отметим еще то, что в некоторых частных случаях при создании геометрической модели возможно клонирование (тиражирование) объектов. Это означает, что, если объекты с одинаковыми геометрическими и физическимн параметрами расположены в пространстве по некоторому правилу, то достаточно определить это правило и описать геометрические размеры одного из этих объектов.
Основные алгоритмы обработки численных результатов При моделировании физических процессов и явлений возникает вопрос о способах визуализации данных в процессе выполнения задачи и после окон чания расчетов. Это могут быть как некоторые численные представления, например таблицы, так и визуальное графическое представление резуль татов: одномерные и двумерные графики, гистограммы, изолинии, изопо верхности, векторные поля и т.д. Так как численные данные — это набор некоторых чисел, то вопрос о их представлении в виде чисел напрямую свя зан с понятием таблиц, баз данных. В большинстве случаев пользователя интересует графическое распределение данных, направления потока, области распределения температур, векторные поля и многое другое. При этом, если в 2-мерном случае геометрия расчетной области и координатная сетка могут быть непосредственно отражены па экране, то в 3-х мерном случае расчетная область отображается в виде некоторых проекций или сечений. Существует множество программных пакетов, позволяющих визуализи ровать данные, начиная от простых типа Excel и заканчивая более уни версальными, например Matlab [38]. Но при разработке собственного программного пакета по математическому моделированию желательно иметь встроенный в него постпроцессор данных, если не такой профессиональный, как Matlab, то с функциональными возможностями, отвечающими требованиям решаемых задач. Для этого разработаны и реализованы алгоритмы по построениям изолиний, изоиоверхностей, полей градиентов, линий тока. Введем ряд определений. Изолиния (для функции f(x,y)) — это геометрическое место точек ( : у) — 0 на плоскости Оху% на которой значение функции равно заданной константе, т.е. f{x,y) — const. Для поля температуры обычно употребляется термин "изотерма , для электрического потенциала — эквипотен циаль, а более общее название — линия одного уровня, употребляющееся в геодезии (одинаковая высота на географических картах). В дальнейшем,для краткости, везде будем применять термин изолиния. Изолинии в сечении — это набор изолиний в плоскости, проходя и ірії. перпендикулярно заданной координатной оси через заданную точку. Изоповерхносгпъ (для функции f(x,y,z)) — это геометрическое место точек р(:г, у, ) — 0 в пространстве, на которой значение функции равно заданной константе: f(x, у, z) — const. Пространство поля градиентов (для функции /(a;, y,z)) — множество точек в пространстве, заданных некоторым образом, и векторов, направленных из них в противоположную сторону направления градиента функции f(oc,y, z). Длина векторов пропорциональна норме градиента. Поле градиентов в сечении — это пространство поля градиентов в плос кости проходящей перпендикулярно заданной координатной оси через заданную точку. Линия тока (силовая линия) (для функции f(x,y,z)) — это линия l(lx(t)fly(t),lz(t)) в пространстве, проходящая через заданную точку ро(гг0,2/о,2о)- В каждой точке линии тока градиент функции / является касательным к этой линии I: ( ; - ; Щ-) — V/. Вышеприведенные определения используют некоторые аналитические функции координат. Однако при численных расчетах мы имеем дело с кусочно-полиномиальном представлением сеточных функций. Для сеточ ных функций разработан постпроцессинг — алгоритмы построения изоли ний, изоповерхностей, полей градиентов, а также линий тока для сеточных функций. Полагаем, что известна сеточная функция / = {/г ,к} заданая на некоторой параллеленипедалыюй сетке Qh — ({#i}Ui, {VJ}JLI, izkVk=i)- e зультатом работы алгоритмов являются наборы координат отрезков, треугольников, а также соответствующие им цвета. Для их отображения на экран можно использовать различные программные инструменты — встраиваемые графические библиотеки для языков программирования, графические пакеты типа "Visual С", "Java 3D" и др. Например, в пакете Базис-А [18] для отображения результатов работы иостпроцесспнга используется приложение "Java 3D".
Предположим, что задана двумерная сеточная функция {/ } на прямоугольной сетке Q,h — ({xt}f=i, {yj}fLi) и значение с, которое должна принимать функция / на изолинии. Первый алгоритм построения изолиний состоит в следующем. Перебором всех сеточных элементов Ег } = {(xt % I+I)J ( У У;+і)} находим первый элемент, через который проходит изолиния. Линейно интерполируем функцию / внутри этого элемента и находим координаты начала и конца изолинии в этом элементе, запоминаем эти координаты (рис.б.а). Затем таким же образом, рассматриваем элемент, прилегающий к найденному, в который продолжается изолиния. Тем самым мы найдем набор отрезков, задающих эту изолинию. Однако, может быть еще несколько изолиний для заданной константы с и их тоже необходимо учитывать. Более того, допускается пересечение и неоднозначность при построении изолинии в некотором элементе (рис.б.Ъ).
