Математическое моделирование многомерных квазистационарных электромагнитных полей в канале электродинамического ускорителя

Размещено на http://www.allbest.ru/

На правах рукописи

Специальность 05.13.18 - математическое моделирование, численные методы и комплексы программ

АВТОРЕФЕРАТ

диссертации на соискание ученой степени кандидата физико-математических наук

Математическое моделирование многомерных квазистационарных электромагнитных полей в канале электродинамического ускорителя

Уразов Сергей Сергеевич

Москва - 2007

Работа выполнена в Ордена Ленина Институте прикладной математики им. М.В. Келдыша Российской Академии Наук

Научный руководитель: доктор физико-математических наук, профессор Галанин Михаил Павлович

Официальные оппоненты: доктор технических наук, ст.н.с. Станкевич Игорь Васильевич доктор физико-математических наук, профессор Фаворский Антон Павлович

Ведущая организация: Московский инженерно-физический институт (государственный университет).

Защита диссертации состоится «__» _______ 200_ г. на заседании диссертационного совета № Д002.024.02 при Институте прикладной математики им. М.В. Келдыша РАН по адресу: 125047, Москва, Миусская пл., д. 4.

С диссертацией можно ознакомиться в библиотеке ИПМ им. М.В. Келдыша РАН.

Автореферат разослан «__» ________ 200_ г.

Ученый секретарь диссертационного совета кандидат физико-математических наук О.В. Щерица

1. ОБЩАЯ ХАРАКТЕРИСТИКА РАБОТЫ

Актуальность темы

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

Цели работы

Целями работы являются:

- развитие математических моделей для описания электромагнитных полей в трехмерных физически и геометрически неоднородных областях канала ускорителей (в том числе с несвязными, негладкими или изменяющимися во времени границами подобластей);

- разработка вычислительных алгоритмов для моделирования явлений в указанных областях;

- проведение комплекса расчетов с использованием разработанного программного обеспечения для исследования разгона макротел в ускорителях различной конфигурации.

Научная новизна

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

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

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

Практическая ценность

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

Апробация работы

Основные результаты работы докладывались на семинарах в Институте прикладной математики им. М.В. Келдыша РАН; на Двенадцатой Международной конференции по вычислительной механике и современным прикладным программным системам, Владимир, 30 июня - 5 июля 2003 г.; X-й научной конференции “Современные проблемы вычислительной математики и математической физики”, посвященной восьмидесятипятилетию академика А.А. Самарского, 24 - 25 февраля 2004 г., ВМ и К МГУ им. М.В. Ломоносова; Международной конференции "Проблемы численного анализа и прикладной математики", Львов, 13 - 16 сентября 2004 г.; Научной сессии МИФИ - 2005 ”Физико - технические проблемы нетрадиционной энергетики и мощная импульсная электрофизика”; Международном научном симпозиуме "Теплофизика и термодинамика ракетно-космических систем" (к 90 - летию профессора Синярева Г.Б.), 3 - 5 октября 2005 г., МГТУ им. Н.Э. Баумана; Международной конференции "Тихонов и современная математика", Москва, 19 - 25 июня 2006 г., ВМ и К МГУ им. М.В. Ломоносова.

Публикации

Результаты выполненной работы представлены в 11 печатных работах. Среди них 3 журнальные статьи, 4 препринта, 4 публикации тезисов докладов на конференциях.

Структура диссертации

Диссертация состоит из введения, четырех глав, заключения и списка использованной литературы. Работа изложена на 120 страницах и содержит 57 рисунков, а также библиографию из 67 ссылок.

2. СОДЕРЖАНИЕ РАБОТЫ

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

Приведены основные положения математической модели, взятой за основу диссертации, и методов ее численной реализации.

Рис. 1. Принципиальная схема рельсотрона: 1 -- направляющий и токоподводящий рельс; 2 -- ускоряемое тело (якорь или иная токовая арматура), 3 -- силовой бандаж канала, 4 -- изолятор.

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

Основное преимущество рельсотрона над распространенными пороховыми или легкогазовыми ускорителями заключается в отсутствии фундаментальных ограничений на величину скорости метания. Диапазон воспроизводимых результатов, полученных при разгоне твердых тел в рельсотроне, простирается до 8.5 км/с. Тем не менее при ускорении металлического якоря силой Лоренца воспроизводимость результатов ограничена скоростями 2.5--3 км/с. В этом случае при достижении скорости порядка 1 км/с наблюдается нарушение металлической проводимости в скользящем контакте между рельсом и якорем. Это приводит к эрозии контакта и ограничению ресурса работы ускорительного канала. Исследование динамики разрушения контакта позволяет разработать соответствующие меры по увеличению скорости макротел в таких устройствах, имеющих максимально простую конструкцию.

Исследования прикладного характера в области электромагнитного ускорения тел начаты в середине 80-х годов. В Институте прикладной математики им. М.В. Келдыша РАН такие исследования были инициированы академиком А.Н. Тихоновым.

Наиболее важные результаты теоретических и практических исследований в России получены группами Ю.П. Бабакова, М.П. Галанина и Ю.П. Попова, Э.М. Дробышевского, А.Д. Лебедева, А.П. Лотоцкого, В.В. Савичева, Б.Д. Урюкова, В.Е. Фортова, Э.Я. Школьникова, Г.А. Шнеерсона и многих других.

Одним из мощных инструментов исследования здесь является численное моделирование. Для описания электромагнитных полей в диссертации использовано так называемое квазистационарное или МГД- приближение [1] уравнений Максвелла. Замкнутая пространственно трехмерная нелинейная и нестационарная во времени модель [1, 4 - 6] объединяет расчет распределения тока в проводниках и локального тепловыделения с непосредственным расчетом ускоряющей силы. При этом модель однородна по различным подобластям с резко различающейся электропроводностью: типа проводник или диэлектрик.

C математической точки зрения электромагнитные поля в МГД - приближении описываются следующей системой уравнений:

rot H = 4 ?? E,

rot E - rot [u H] = - ,

div Н = 0, j = ?E.

Здесь и далее E и H -- векторы напряженности электрического и магнитного полей соответственно, j -- вектор плотности тока, ? -- электропроводность, u -- вектор скорости движения вещества, r = (x, y, z) -- радиус-вектор, t -- время. Система уравнений (0.1) записана в безразмерном виде. В ней E -- напряженность электрического поля в системе координат, в которой вещество покоится. Будем обозначать через E* напряженность электрического поля в неподвижной (лабораторной) системе координат (). Здесь и всюду ниже в формульных выражениях все величины даются в безразмерном виде (в частности, в (1) H = B, где B -- вектор магнитной индукции).

Приведем постановку задачи [1, 4 - 6] для определения электромагнитных полей внутри области после введения векторного потенциала A:

где векторный потенциал А есть решение следующей задачи:

Здесь учтена неоднородность задачи по пространству: () = 0 в G1 и () = 1 в G2. В (2) G = G1 G2, G -- рассматриваемая область, G1 = {r G: > 0}, G2 = {r G: = 0}, G1 и G2 -- границы G1 и G2 соответственно, G12 = G1 ? G2, Г1 -- часть общей границы G, на которой задано условие для E* (то есть для A), Г2 -- часть G, на которой задано условие для H ( -- известная вектор-функция), G = Г1 Г2, Г12 = Г1 ? G2, 12 = G12 Г12. В записи (2) использованы смешанные эйлерово-лагранжевые (СЭЛ) переменные: D/Dt = /t + (v,), где /t -- производная при фиксированных эйлеровых переменных, D/Dt -- при фиксированных СЭЛ - переменных; v -- скорость движения точек пространственной области (в нашем случае v -- скорость движения якоря как целого, независящая от координат пространственной точки). Индекс n указывает на нормальную по отношению к границе составляющую вектора, ?? -- тангенциальную. В рассматриваемых задачах в декартовой системе координат движение якоря происходит в положительном направлении оси y.

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

В силу геометрической симметрии достаточно найти решение задачи в правой верхней четверти расчетной области в трехмерном случае (или в верхней половине области в двумерном случае).

При разработке модели использовано резкое различие длины ускорителя (по y) и его поперечных размеров. Учтено также, что единственной заданной извне электромагнитной величиной можно считать полный ток, определяемый источником питания. При таком подходе возникает проблема задания граничных условий на передней и задней границах исследуемой области. На боковых границах этой проблемы нет, так как канал рельсотрона обычно заключен в проводящий силовой бандаж. Поэтому естественно рассмотреть модель [4 - 6], в которой на торцах расчетной области заданы тангенциальные компоненты магнитного поля, соответствующие бесконечно длинной (вдоль оси y) системе проводников, для каждого из которых задан полный ток. Согласно этой модели в области G имеется N проводников, по которым протекают заданные токи , где N , т.е. , - область, занимаемая k - ым проводником, - его боковая поверхность, . Поле, необходимое для определения граничных тангенциальных компонент напряженности магнитного поля , является решением соответствующей пространственно двумерной задачи:

в Sk , k=1, 2 … N,

A = 0 в G2 ,

Здесь [A]- скачок вектора при переходе через границу.

В результате для постановки граничных условий необходимо решить две специальные задачи для интегро-дифференциальных уравнений на торцах [4 - 6]. Решение же трехмерной задачи получается путем использования алгоритма [1, 4 - 6] по заданным тангенциальным компонентам магнитного поля. При решении задач использовался метод конечных разностей. Разностная задача формируется при помощи метода опорных операторов [1, 7, 8]. Разностные операторы, аппроксимирующие основные операторы векторного анализа (rоt, div, grad) [1], строятся на основе инвариантных определений, не зависящих от выбора системы координат. В итоге использованная для расчета векторного потенциала разностная схема [1] в общем случае имеет вид:

.

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

Для решения возникающей при моделировании системы линейных алгебраических уравнений использован метод сопряженных градиентов совместно с неполным разложением Xолесского [9 - 11].

Для расчета температурного поля в проводящей области применяется математическая модель[1] в СЭЛ переменных:

.

Здесь w = u - v -- вектор относительной скорости вещества в движущемся со скоростью v объеме, с -- плотность вещества,? -- удельная внутренняя энергия, cv -- удельная теплоемкость, W = - ? grad T -- вектор теплового потока, ? -- коэффициент теплопроводности, T -- температура, (j, E) -- мощность тепловыделения за счет джоулева нагрева. Расчет температурного поля ведется параллельно с расчетом других полей (E, j, H, A).

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

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

Представлен универсальный для различных конфигураций объекта исследования метод описания расчетной области, разработаны новые методы описания границ расчетной области и преобразованные в соответствии с ними вычислительные алгоритмы. С их использованием проведено исследование температурного режима электрического контакта металлического типа с переходом к гибридному и далее к плазменному (см. [16 - 19, 21, 26]).

В экспериментах, выполненных группой А. П. Лотоцкого [12, 17], показано, что при скоростях разгона металлического якоря, превышающих 1.5 км/с, тепловыделение в зоне контакта якоря и рельса приводит к разогреву и испарению задней кромки якоря с образованием плазменной перемычки. Этому способствует эффект скоростного скинирования (явление концентрации тока на границе контакта), приводящий к увеличению плотности тока на задних кромках контактов. Точка перехода к дуговому разряду идентифицировалась исследователями как "кризис" металлического контакта, а разрушение поверхностей контактирующих тел - как эрозия.

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

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

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

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

Отмечено, что зона, в которой наблюдался "кризис", практически недоступна для проведения измерений в силу малых размеров и высокой скорости перемещения якоря. Из-за многочисленных трудностей экспериментального исследования процессов в высокоскоростном электрическом контакте анализ механизмов эрозии может быть выполнен только с помощью математического моделирования, учитывающего вклад различных процессов (распределение токов, нагрев материала и фазовые переходы) в существенно многомерной перемещающейся области. Для исследования динамики этих процессов в данной главе выполнено математическое моделирование процесса разгона якоря и теплового режима скользящих проводящих поверхностей. Результаты математического моделирования сопоставлены с экспериментальными данными.

	Рис. 2. Расчетная область: 1 -- рельс, 2 -- U-образный якорь

В расчетах и экспериментах рассматривался U-образный якорь (см. рис. 2). В соответствии с натурными экспериментами расчеты проведены для двух вариантов процесса ускорения: 1 - ускоряется тело массой 5 г с начальной скоростью 250 м/с; 2 - ускоряется тело массой 2.5 г с начальной скоростью 500 м/с. Параметры расчетной области и материалов совпадают с параметрами в эксперименте.

Данные экспериментов для первого варианта свидетельствуют о существовании постоянного металлического контакта якоря с рельсами при ускорении якоря до скорости около 1 км/с. После достижения якорем скорости 650 - 700 м/с задний фронт токового распределения начинает перемещаться по направлению к передней кромке якоря.

В следующей экспериментальной серии (при ускорении якоря до 2.3 - 2.7 км/с) данные магнитозондовых измерений показали, что надежный металлический контакт существовал в течение 100 - 150 мкс, а к моменту времени 350 мкс от старта разряда наступала дуговая стадия. Наблюдалось разрушение якоря в процессе ускорения.

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

С учетом заданной зависимости полного тока от времени рассчитаны временные зависимости сосредоточенных характеристик процесса ускорения: максимальной температуры в области, скорости и координаты якоря. Ускоряющая сила (используемая при расчете скорости) вычислялась интегрированием вектора плотности силы Лоренца f = [j, H] по объему. Из результатов расчета для первого варианта ускорения следует, что при достижении момента времени t = 0.7 мсек, начинается плавление якоря. Якорь покидает ускоритель со скоростью около 900 м/c. В процессе ускорения температура кипения не достигнута. Для второго варианта ускорения (разгон металлического якоря до скорости более 2 км/с) плавление якоря начинается при достижении момента времени t = 0.21 мс, v = 1 км/c. Далее при достижении момента времени t = 0.35 мс и скорости v = 1.8 км/c начинается кипение материала якоря. Об этом свидетельствует и своеобразный "излом" на зависимости максимальной температуры от времени. Кипение и плавление начинаются в задней (по направлению движения) угловой точке якоря, контактирующей с рельсом. Якорь вылетает со скоростью около 3 км/c.

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

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

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

Распределение Нx (безразмерное).	Распределение T (К) .
Рис. 3. Сечение x = const , плоскость контакта якоря и рельса. Распределение компонент физических полей при t = 0.233 (плавление). По осям y,z см.

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

Джоулев нагрев в рассматриваемом случае может привести к частичному или полному разрушению ускоряемого тела и вызвать срыв процесса ускорения. Рассмотрена возможность сквозного проплавления материала ускоряемого U-образного тела. Расчеты показали, что в конце ускорительного канала якорь может полностью разрушиться (распасться на 2 или 3 части). Проходящий через якорь ток концентрируется прежде всего на поверхности проводящего тела, в результате чего внутренняя граница (по оси y) U-образного якоря нагревается до температур, превышающих температуру плавления.

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

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

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

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

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

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

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

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

Рассмотрены явления, способствующие возникновению сил, способных разрушить соединение. Представлены линии уровня распределений плотности тока в различных сечениях. Определены участки, где плотность тока достигает наибольшего по модулю значения. Проанализирован вклад компонент вектора плотности тока в значение его модуля в различных участках стальной и алюминиевых частей. Наличие выступов и угловых точек приводит к появлению вихрей в распределении тока.

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

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

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

Целью третьей главы является разработка методов моделирования процессов в канале ускорителя в случае изменения границ областей при испарении материала и в случае несвязности границ (см. [22, 23, 25]).

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

Все исследуемые виды областей встречаются при исследовании импульсных электродинамических ускорителей типа рельсотрон.

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

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

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

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

Ранее [1] при моделировании испарения проводника электропроводность в испарившейся части полагалась малой ("фоновой") величиной (на несколько порядков меньшей значения до испарения), что обеспечивало соответствующее перераспределение тока. Но при таком подходе образуются проводящие подобласти с малой электропроводностью, которые должны соответствовать диэлектрическим. Использование для моделирования диэлектрика проводника с малой электропроводностью может серьезно ухудшить устойчивость решения по отношению к возмущениям правой части системы уравнений, что неизбежно приводит к плохой сходимости итераций и резкому уменьшению шага по времени. Если величина "фоновой" электропроводности берется большой, то неверно считается тепловыделение и распределение токов. Однако в идеальном случае она должна стать равной нулю. Для преодоления описанных трудностей в данной главе предложена следующая модель: при превышении температуры кипения в проводящей ячейке такая ячейка заменяется диэлектрической с соответствующей перестройкой логических массивов, используемых для описания расчетной области и разностных схем, соответствующих дифференциальным операторам. Использование такой модели позволяет получить нормальное (обладающее минимальной нормой) решение задачи с нулевой "фоновой" электропроводностью.

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

Для сопоставления различных методов моделирования испарения исследованы границы спектра разностного оператора для таких методов.

Рассмотрены задачи в двумерном и трехмерном приближении (двумерное приближение рассматривается для получения качественных картин, основным является трехмерное).

Для иллюстрации на рис. 4 представлена часть двумерной расчетной области и разностной сетки в месте испарения материала проводника. Ось x направлена вертикально, y -- горизонтально. С испарением ячеек меняются границы G1, G2, G12 и ?12.

Рис. 4. а -- расчетная область после испарения материала двух ячеек (испарился материал двух ячеек с общим ребром 2, в задаче нет движения), б -- изменение шаблонов разностной схемы для ребер 1 - 3 (преобразование шаблона разностной схемы для оператора (rot rot) c учетом условий div A = 0 на ?12), в -- изменение шаблона для ребра 4 (к виду шаблона для оператора (rot rot - grad div)). Ребра, для описания которых изменяются шаблоны разностной схемы, обозначены символом Х.

Исследована обусловленность матрицы М системы линейных алгебраических уравнений, получаемых при разностной аппроксимации уравнений Максвелла для различных способов моделирования испарения (см. таб. 1,2).

Расчеты проведены для следующих вариантов: 1 -- в области нет кипения; 2 -- испаряется материал двух крайних ячеек разностной сетки на границе якоря (в двумерном случае) или на ребре якоря (в трехмерном случае), вводится "фоновая" электропроводность ?0 = 10-5 (обычно использовалось в расчетах) и ?0 = 10-10 (вариант 2a), соответственно, при исходной ??>1; 3 -- испаряется материал двух ячеек, и разностные схемы перестраиваются по предложенному алгоритму. Собственные значения приведены для задачи без движения (процессы кипения и испарения проходят в рассматриваемых случаях качественно сходно для различных скоростей движения, в трехмерном случае для моделирования рассматривался U-образный якорь [4 - 6]).

Таблица 1. Собственные значения матрицы M для различных способов моделирования испарения материала в двумерном случае

№ варианта	Наибольшее собственное значение	Наименьшее собственное значение	Число обусловленности
1	52.1784	1.6894·10-2	3.0885·103
2	52.1784	2.5875·10-6	2.0165·107
2a	52.1784	2.5875·10-11	2.0165·1012
3	52.1784	1.6398·10-2	3.1819·103

Таблица 2. Собственные значения матрицы M для различных способов моделирования испарения материала в трехмерном случае

№ варианта	Наибольшее собственное значение	Наименьшее собственное значение	Число обусловленности
1	20.457	5.331·10-5	3.837·105
2	20.457	6.906·10-9	2.963·109
3	20.457	5.326·10-5	3.840·105

Также приведены полученные при численном моделировании изотермы, демонстрирующие картину испарения материала (движения волны кипения) и картины перераспределения магнитного поля при испарении материала в задаче с перестроением разностных схем.

Из приведенных результатов видно, что перестроение разностных схем при испарении материала в двумерном и трехмерном случае позволяет улучшить обусловленность решаемой системы на много порядков. Отметим, что числа обусловленности определяются значением "фоновой" электропроводности. Чем она меньше, тем больше число обусловленности, стремящееся к бесконечности при стремлении "фоновой" величины к нулю. Улучшение обусловленности системы позволяет проводить расчет с большим шагом по времени. Таким образом, при исследовании процессов в ускорителях предпочтительной является модель с перестроением разностных схем и описателей границ при испарении материала в ячейках разностной сетки, относящихся к проводнику.

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

Алгоритм обобщается на другие случаи изменения границ.

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

Интерес представляет исследование единственности решения системы уравнений Максвелла в трехмерной области, в которой граница ?12 несвязна.

Согласно [1] решение системы уравнений Максвелла в квазистационарном приближении в проводнике единственным образом определяется граничными и начальными условиями.

Квазистационарность рассматриваемых полей может привести к неединственности поля E в диэлектрике, причем поле H определяется единственным образом во всей области [1].

Разработаны методы получения единственного квазистационарного решения системы уравнений Максвелла в трехмерной области с многосвязными диэлектрическими и проводящими подобластями. Приведен метод построения системы независимых векторов, представимых в виде градиента скалярной функции для случая m-связной границы ?12. Для двусвязной и трехсвязной областей построены векторы, определяющие подпространство, соответствующее ядру оператора квазистационарного приближения уравнений Максвелла. При этом использовалось представление разности решений уравнений Максвелла в виде градиента потенциала.

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

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

Представлены алгоритмы получения нормального решения задачи, принадлежащего пространству, ортогональному ядру оператора задачи (в m-связной области), а именно: построены алгоритмы преобразования исходной системы уравнений (разностной схемы) к виду, дающему единственное нормальное решение задачи. Показана положительная определенность матриц преобразованных систем. Для получения одного из решений и последующего нахождения нормального построены простейший вид преобразования системы (с сохранением симметричности матрицы) и оператор, строящий проекцию полученного решения на подпространство, ортогональное ядру разностного оператора M. Предложенные алгоритмы можно применять для нахождения нормального решения и в других случаях вырожденной матрицы M. Алгоритмы решения системы уравнений Максвелла в многосвязной области, позволяющие получить единственное решение, доведены до программной реализации. Получены результаты расчета электромагнитных полей по построенным алгоритмам в областях с несвязными границами в нестационарном случае.

В четвертой главе рассмотрена задача численного моделирования квазистационарных электромагнитных полей в областях с негладкими границами проводящих и диэлектрических подобластей (см. [24]).

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

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

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

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

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

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

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

(E= - DA/Dt + (v,)A + [u, rot A] + grad ? = - DA/Dt + grad (v, A) +[w, rot A] + grad ?).

При выборе кулоновской калибровки векторный потенциал А является решением задачи (2).

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

Напряженность магнитного поля не изменится, если в (1) вместо нулевого взять любое согласующееся с граничными условиями для Е значение скалярного потенциала.

Показан вид потенциала (? = - (v,A)), при использовании которого (в данной задаче) конвективные слагаемые будут входить только в уравнения для A в рельсе.

Тогда поля описываются уравнениями:

4??? (- DA/Dt +[w, rot A]) = rot rot A - ??(?) grad div A

или 4??? (- DA/Dt ) = rot rot A в якоре;

4??? (- DA/Dt - [v, rot A]) = rot rot A в рельсе

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

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

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

Проведено сравнение числа итераций для различных способов моделирования в двумерном и трехмерном случае. На рис. 5 приведен пример полученной в расчете зависимости суммарного количества итераций (внешних и внутренних) от времени для двух вариантов калибровки с ? = 0 и с ? = -(v,A) для двумерного случая.

Рис. 5.Суммарное количество итераций.
Калибровка ? = 0.	Калибровка ? = -(v,A).

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

При явном выделении особенности использована возможность выделения сингулярной части решения в явном виде [15].

Решение представлено в виде суммы негладкой (выделяющей особенность и обозначенной через A0) и гладкой (A*) частей A= A0 + A0*. При таком задании A0 выполняются требования: div A0 = 0 и rot A0 = 0 в диэлектрической подобласти. Получено решение с устраненными экстремумами вблизи угловой точки.

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

Работы 2003 - 2007 годов поддержаны грантами Российского фонда фундаментальных исследований (проекты РФФИ № 03 - 01 - 00461 и РФФИ № 06 - 01 - 00421) и Фонда содействия отечественной науке, а также грантами в рамках Программы № 3 ОМН РАН (проект № 3.2).

В заключении приведены основные результаты диссертации, а именно:

- разработаны методы математического моделирования квазистационарных электромагнитных полей в неоднородных областях канала ускорителя (в том числе с изменяющимися во времени, несвязными и негладкими границами подобластей);

- построены и программно реализованы вычислительные алгоритмы для моделирования процесса электромагнитного ускорения в указанных областях;

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

 ЛИТЕРАТУРА

1. М. П. Галанин, Ю. П. Попов. Квазистационарные электромагнитные поля в неоднородных средах: Математическое моделирование. М.: Наука, 1995. 320 с.

2. Материалы I Всесоюзного семинара по динамике сильноточного дугового разряда в магнитном поле (Новосибирск, 10 - 13 апреля 1990 г.). / Под ред. М. Ф. Жукова Новосибирск: Изд. ИТ СО АН СССР, 1990. 350 с.

3. Материалы II Всесоюзного семинара по динамике сильноточного дугового разряда в магнитном поле (Новосибирск, 4 - 6 декабря 1991 г.). / Под ред. В. Е. Накорякова Новосибирск: Изд. ИТ СО РАН. 1992. 367 с.

4. М. П. Галанин, А. П. Лотоцкий, Ю. П. Попов, С. С. Храмцовский. Численное моделирование пространственно трехмерных явлений при электромагнитном ускорении проводящих макротел // Математическое моделирование. 1999. Т. 11, № 8. С.c. 3 - 22.

5. М. П. Галанин. Компьютерное моделирование в задачах конвертирования электромагнитной и кинетической энергии. Задачи и модели. // Информационные технологии и вычислительные системы. 2002. № 4. С.c. 109 - 123.

...