Введение. Постановка задачи

 

Настоящая работа посвящена исследованию тепловых процессов на контактах электродинамического ускорителя макротел типа рельсотрон [1,2], типичная схема которого представлена на рис. 0.1. Интерес к подобной задаче вызван следующими обстоятельствами. Современные ускорители макротел массой 10^-3 – 1 кг, типичными представителями которых являются различные виды метателей, использующих газодинамическое давление продуктов сгорания химических соединений, обладают скоростью снаряда, ограниченной величиной 1-1,8 км/с. Электродинамический рельсовый ускоритель свободен от этих ограничений. Здесь ускоряющая сила обеспечивается давлением магнитного поля или, иначе, силой Лоренца. Величина силы не имеет явной зависимости от скорости тела, которая в этом случае могла бы иметь верхний предел, ограниченный лишь релятивистскими эффектами. Так, плазменные сгустки удается разгонять до скоростей вплоть до 1000 км/с. К сожалению, реальная ситуация с разгоном твердых тел оказалась много сложнее, чем представлялось на первый взгляд.

Рис. 0.1. 1- рельс, 2- якорь.

 

Уже первые удачные эксперименты показали, что металлический контакт между якорем и рельсами разрушается, а переход к плазменной токовой перемычке сопровождается резким падением КПД преобразования электрической энергии в кинетическую. С ростом скорости контакт якорь/рельс теряет свойства чисто металлического и превращается в смешанный (появляется дуга, т.е. плазма). Это приводит к разрушению как якоря, так и рельса, и не позволяет добиться желаемого результата (электрическая энергия, запасенная в источнике, не преобразуется в кинетическую энергию якоря). Разрушение металлического электрического контакта является результатом действия многих физических факторов, действующих в наиболее напряженном (во всех отношениях) месте ускорителя – месте контакта. Главными из них являются электромагнитные явления и вызванный ими джоулев нагрев. Однако тепловыделение за счет трения также вносит значительный вклад в плавление и испарение материала контактирующих тел. Подробнее физические явления на контактах электромагнитного ускорителя описаны в [3]. Целью данной работы является изучение вклада трения в процесс разрушения электрического металлического контакта.

Рис. 0.2. Геометрия задачи.

 

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

Достаточно полная математическая модель явлений, протекающих при электромагнитном ускорении проводящих макротел, представлена в [4]. Для определения роли одного лишь трения в отдельности ток в модели [4] положен равным нулю. Скорость якоря будем считать заданной. Мы ограничимся пространственно двумерным приближением. Соответствующая геометрия представлена на рис. 0.2.

Поскольку электромагнитное поле отсутствует, то фактически от модели [4] остается одно уравнение энергии - уравнение теплопроводности. С учетом движения проводящего якоря оно принимает следующий вид:

                              (0.1)

Уравнение (0.1) записано с применением смешанных эйлерово - лагранжевых (СЭЛ) переменных. Производная D/Dt означает производную по времени при фиксированных СЭЛ - переменных. В (0.1) r - плотность, e - удельная внутренняя энергия, u - скорость движения вещества, v - скорость движения якоря, k  - удельная теплопроводность, T - температура,  Q_tr - мощность тепловыделения, c_V - удельная теплоемкость при постоянном объеме, k_t - коэффициент трения, s_n - нормальное механическое давление на контакте, d - дельта – функция Дирака.

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

Рис. 1.1. Треугольник и части трех ячеек Дирихле.

 

§ 1. Численный алгоритм.

 

Подробно использованный вычислительный алгоритм описан в работах [4-6]. Укажем здесь лишь некоторые его особенности. Алгоритм для реализации модели на ЭВМ основан на методе конечных элементов [7,8] и удовлетворяет требованиям полной консервативности и квазимонотонности. Для аппроксимации величин применялись два семейства функций, соответствующих разбиению области на треугольники и ячейки Дирихле. На рис. 1.1 представлен один такой треугольник и его разбиение на три части перпендикулярами к серединам сторон, образующими стороны ячеек Дирихле, пересекающихся с данным треугольником. Первое семейство – это плоскость в пределах треугольника, проходящая через 1 в одной из вершин и 0 в двух других. Второе – это характеристическая функция ячеек Дирихле, т.е. 1 внутри ячейки и 0 вне ее. Система линейных алгебраических уравнений для амплитуд конечных элементов получена с помощью алгоритмов Галеркина и Галеркина – Петрова [7].

Предварительные расчеты температуры контакта показали, что в некоторых случаях численное решение может быть немонотонным по пространству – см. рис. 1.2.а. Это связано с использованием при дискретизации пространственной части уравнения теплопроводности алгоритма Галеркина - Петрова и интерполяции первого вида. Данный подход приводит к аппроксимации первой производной по пространству, которая похожа на схему типа центральной разности. Известно [9], что такие схемы дают немонотонные решения.

Чтобы избежать немонотонности решения, была введена искусственная теплопроводность (ИТ). В правую часть уравнения теплопроводности при этом добавляется член div (|w|rc_vd_T grad T). Тогда для обеспечения квазимонотонности решения (т.е удовлетворения схемой условиям принципа максимума) в каждой из частей ячейки Дирихле имеем следующие условия:

                                      (1.1)

Они выписаны на примере части ячейки Дирихле, прилегающей к 1 - ой точке в рассматриваемом треугольнике с рис. 1.1.

Здесь S₁ – площадь части рассматриваемой ячейки Дирихле в данном треугольнике, n₁ и n₂ – внешние нормали к ее граням в треугольнике. l₁ и l₃ - длины этих граней, j – интерполяционные функции конечных элементов, привязанные к треугольникам, b_i – коэффициенты, с помощью которых задаются интерполянты j (приведены ниже). Условия (1.1) выписаны для каждого коэффициента перед амплитудой конечного элемента по отдельности. В каждом из треугольников функции j имеют вид:

                (1.2)

где w_j – треугольник номера j, - множество треугольников, имеющих точку i₁ своей вершиной. Необходимо добавить еще две тройки неравенств для частей ячеек Дирихле, прилегающих к точкам 2 и 3. Минимальное неотрицательное решение получаемой системы девяти уравнений и дает значение искомой ИТ в данном треугольнике.

a)	б)
в)

 

Рис. 1.2. Немонотонность распределения температуры по пространству (изотермы).

 

Необходимость во введении ИТ возникает, если условие h < k/r cv, где h – характерный линейный размер ячейки, не выполняется. Минимальное значение h ограничено методом генерирования ячеек (т.е. возможностями используемой ЭВМ). В описываемых расчетах оно составляло около 4^.10^-6 м, поэтому необходимо применять ИТ (правая часть последнего неравенства в данном случае около 10^-7 м).

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

Однако возможно выбрать нелинейную ИТ, величина которой зависит от решения в данном треугольнике, аналогично схеме с "лимитерами" из [10]. При этом принцип максимума в каждой ячейке остается справедливым. Данный подход, более подробно описанный в [4-6], позволяет избавиться от немонотонности решения корректным образом. Результат расчета решения с нелинейной ИТ представлен на рис. 1.2.в. Видно, что, с одной стороны, решение является монотонным, с другой – полоса распространения тепла соответствует расчету без ИТ.

 

§ 2. Приближенные аналитические решения.

 

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

2.1. Распределение температуры вдоль поверхности контакта.

В предположении о том, что теплопроводностью вдоль контакта можно пренебречь, можно найти аналитическое решение уравнения (0.1).

При таком допущении (0.1) может быть переписано в виде следующей системы:

                                                                (2.1)

Применение интегрального преобразования Лапласа позволяет найти точное решение (2.1) для температуры вдоль поверхности контакта (y = 0). Оказывается, что такая температура равна:

                        (2.2)

Линейная координата x (расстояние от начала якоря по направлению движения) показана на рис. 0.2.

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

2.2. Решение при наличии скачка теплового потока на контактной границе двух бесконечных полупространств.

Данная постановка задача означает отсутствие в уравнении теплопроводности (0.1) производных температуры по координате вдоль поверхности контакта. В этом случае также можно получить точное решение:

                                      (2.3)

2.3. Снос тепла от источника.

Для выделения роли конвективного переноса полезно рассмотреть исходную задачу в следующей постановке. Опустим в (0.1) слагаемые с производными по x (поперек контакта). При этом фактически имеем задачу о нахождении распределения тепла от источника мощностью q, находящегося на участке (0, l) контакта между двумя полупространствами. Иллюстрация рассматриваемой задачи приведена на рис. 2.1.

Рис. 2.1.

                                                                             (2.4)

Точное решение (2.4) записывается следующим образом:

                                (2.5)

Для его исследования рассмотрим следующие упрощения.

А) Случай нулевой теплопроводности (a = 0).

Тогда имеем решение вида:

                               (2.6)

Здесь H(x) - функция Хевисайда.

Рис. 2.2 иллюстрирует пространственные распределения температуры для различных моментов времени в соответствии с приведенным решением. На данном рис. представлено решение (2.6) в двух пространственных точках в зависимости от времени (два верхних рисунка) и две картины распределения температуры по пространству в два характерных момента времени (два нижних рисунка).

Б) Случай малой теплопроводности.

Перепишем решение (2.5) в следующем виде:

Его можно представить как сумму  - главной части решения - и двух интегралов I₁ и I₂:

           (2.7)

Для части решения I₁ при малых значениях a может быть получена следующая асимптотическое соотношение:

                                      (2.8)

Рис. 2.2.

 

Для части решения I₂ верно все то же самое с заменой y на y-l.

Видно, что при стремлении a к нулю главная часть решения (2.7) переходит в соотношение (2.6).

Асимптотика (2.8) получена в соответствии с методами [11].

 

§ 3. Численные расчеты

 

1. На рис. 3.1 сопоставлены продольные распределения температуры, полученные численно (с нелинейной ИТ) и аналитически (решение (2.2)), на момент времени 0.3 ms. При этом выбраны реальные физические параметры рельсотрона. Главные из них таковы: материал рельса - медь, материал якоря – алюминий, скорость якоря - 1000 м/с, тепловой поток, вызванный трением – 5^.10¹⁰ Вт/м².

Рис. 3.1.

 

Рис. 3.1 демонстрирует превышение аналитическим решением численного, что свидетельствует о значительности теплопроводности в продольном направлении, которая не учитывается решением (2.2).

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

 

Рис. 3.2

 

2. На рис. 3.3. сопоставлены результаты численных расчетов по двумерной модели с решением (2.6). На этом рисунке показаны распределения температуры вдоль поверхности контакта на два различных момента времени. Скорость движения якоря была постоянной и равнялась 1000 м/с. Теплофизические свойства якоря и рельса в численной модели, а также длина якоря, были взяты как в реальном рельсовом ускорителе. Величина q из уравнения (2.6) определялась через тангенс угла наклона правого фронта численной кривой распределения температуры на рис. 3.3.а. Якорь (зона источника тепла) занимает отрезок от 0.0 до 0.155 м. Зона, для которой численно решалась задача определения температуры, занимает отрезок от –0.004 м до 0.0165 м, поэтому численные кривые обрываются слева.

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

а) t=12,3 мкс, (случай t<l/v)
б) t=22,8 мкс, (t>l/v)

Рис. 3.3. Сопоставление численно полученного распределения температуры на контакте со случаем чисто конвективного уноса тепла из зоны тепловыделения. Сплошные кривые – численный расчет, пунктир – решение (2.6).

 

С другой стороны, отличия расчетного и аналитического решений обусловлены различием численно решаемой задачи (0.1) и задачи (2.4). Главные из этих отличий: двумерность основной задачи (0.1) и присутствие теплопроводности во всех направлениях, различие источников тепла. В задаче (2.4) величина q есть мощность тепловыделения, в то время как в (0.1) Q – скачок теплового потока на границе раздела проводников.

При этом в зоне якоря в численных расчетах наблюдается тенденция к выходу распределения температуры контакта на стационарное распределение (см. численную кривую на рис. 3.3.б).

 

Заключение

 

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

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

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

 

Литература

 

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

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

3. М.П. Галанин, А.Д. Лебедев, А.П. Лотоцкий, К.К. Миляев. Тепловые и электромагнитные процессы на контактах электродинамического ускорителя // Препр. Ин. прикл. матем. им. М.В.Келдыша РАН. 2000. N 42. 32 с.

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

5. Галанин М.П. Нелинейная квазимонотонная конечноэлементная схема для решения двумерных задач с переносом. Задача о скоростном скин–слое // Дифференциальные уравнения. 1996. Т. 32. N 7. С. 935-942.

6. Галанин М.П. Задача о скоростном скин - слое и квазимонотонная конечноэлементная схема для ее численного решения. ускорителя // Препр. Ин. прикл. матем. им. М.В.Келдыша РАН. 1995. N 21. 29 с.

7. Г.И. Марчук, В.И. Агошков. Введение в проекционно - сеточные методы. М., Наука, 1981, 416 с.

8. П. Сильвестер, Р. Феррари. Метод конечных элементов для радиоинженеров и инженеров – электриков. М., Мир, 1986, 229 с.

9. А.А. Самарский. Теория разностных схем. М., Наука, 1978, 656 с.

10. К.В. Вязников, В.Ф. Тишкин, А.П. Фаворский, М.Ю. Шашков. Квазимонотонные разностные схемы повышенного порядка для ее точности // Препр. Ин. прикл. матем. им. М.В.Келдыша АН СССР. 1987. N 36. 27 с.

11. М.В. Федорюк. Асимптотика: Интегралы и ряды. М., Наука, 1987, 544 с.