CC BY
33
УДК 51-74; 62-768.3
И. В. Плохов, А. В. Ильин, О. И. Козырева
СТРУКТУРА И АЛГОРИТМЫ ИМИТАЦИОННОГО МОДЕЛИРОВАНИЯ ДИНАМИКИ ЭЛЕКТРОФРИКЦИОННОГО
ВЗАИМОДЕЙСТВИЯ
Рассматривается электрофрикционное взаимодействие в электрических машинах. Приводятся основные алгоритмические действия и математические зависимости выполняемые в процессе имитационного моделирования.
Ключевые слова: электрофрикционное взаимодействие, скользящий электрический контакт, компоненты вектора состояния, динамическая модель, фрактальные кластеры.
Моделирование динамики электрофрикционного взаимодействия (ЭФВ) представляет собой многофакторную задачу исследования процессов формирования проводимости в переходном слое скользящего контакта (СК).
В данной статье уделим внимание заданию структуры и алгоритмов моделирования динамики ЭФВ.
Приведём блоки основных алгоритмических действий, выполняемых в процессе имитационного моделирования [1, 2]. Совокупность указанных блоков опишем списком операций:
Блок А. Задание начальных значений компонент векторов состояния контактного элемента (КЭ).
1. Г енерация двумерных массивов ячеек, задающих микрорельефы контактирующих поверхностей с учетом поверхностных пленок.
2. Задание параметров материалов и поверхностных плёнок.
3. Сближение микрорельефов на расстояние, определяемое по величинам усилия нажатия и контактной жёсткости.
4. Заполнение двумерного массива, определяющего распределение кластеров непосредственной и плёночной проводимости;
Блок Б. Формирование и приложение силового вектора управляющих и возмущающих воздействий.
1. Задание скорости скольжения и дискретизация по времени и пути.
2. Задание перемещения микрорельефов путём пошагового переприсвое-ния значений массивов.
3. Задание усилия нажатия и начального сближения микрорельефов.
4. Задание величины электрического напряжения.
5. Задание начальных температур;
Блок В. Модификация компонент векторов состояния КЭ.
1. Вычисление сопротивлений стягивания контактных кластеров.
2. Расчёт величин токов проводящих кластеров.
3. Расчёт мощностей тепловыделения в кластерах.
4. Расчёт нестационарного температурного поля.
5. Изменение контактной жесткости в зависимости от температуры.
6. Изменение контактного сближения в зависимости от контактной жесткости и контактного давления.
7. Изменение величин и типов проводимости в зависимости от температуры.
8. Изменение проводимости в результате фриттинга.
9. Изменение толщины поверхностных плёнок в результате химических, тепловых и механических процессов;
Блок Г. Вычисление и представление интегральных характеристик.
1. Вычисление суммарного тока и переходного падения напряжения.
2. Определение средней температуры кластеров и всего контакта.
3. Вычисление физической (фактической) площади контакта.
4. Построение динамических вольт-амперных аттракторов, графиков тока, температуры и др.
5. Визуализация микрорельефов, картины температурного поля и контактной проводимости.
6. Отображение пользовательских интерфейсов.
Опишем подробно перечисленные операции.
Задание начальных значений компонент векторов состояния контактных элементов
Первым действием вычислительной модели является генерация фрактальных микрорельефов, которая заключается в заполнении двумерных массивов информацией о превышениях поверхностей относительно опорных плоскостей. Заполнение производится по алгоритму Фосса «случайное сложение» [3], который обобщён на большее количество измерений. Модифицируем данный алгоритм применительно к решаемой задаче.
1. Заполним один столбец матрицы превышений микрорельефа:
\] +1 = a ■ \] + Ь ■ ^, при ^ е ]-1,+1[, (1)
где а, Ь — коэффициенты, определяющие параметры шероховатости.
2. Заполним оставшиеся столбцы по алгоритму:
\j = Ыыгр (Н,-и, \j-1, И,J+Д.+и, И,+1,j-1, И,-и+1, И,-и-1, И,+lJ+1) + Ь ■ j, (2)
где Interp — процедура трёхмерной линейной интерполяции по указанным восьми точкам.
3. Имитируем сближение микрорельефов. При этом в соприкосновение вступает всё большее количество КЭ и через определённое число шагов наступает равновесие усилия нажатия и силы упругого противодействия КЭ. В процедуре используем метод половинного деления.
4. Для получения проводящих кластеров различной природы примем, что изначально вся контактная поверхность покрыта равномерным слоем окисной плёнки и имеет хаотические вкрапления продуктов износа трёх типов — материала щётки, контактного кольца и окислов. Поэтому, после завершения проце-
дуры сближения, те элементы, которые сблизились на расстояние, превышающее толщину плёнки и продуктов износа, наделим непосредственной проводимостью (в матрицу записываем «-1»). Для КЭ, не преодолевших толщину плёнки и продуктов износа, запишем типы и толщину слоёв. Для КЭ, не вступивших в соприкосновение, запишем нули.
5. Матрицу теплового состояния заполним начальным значением температуры и, в соответствии с её величиной, установим скорость изменения толщины окисной пленки.
Примем, что компонентами силового вектора F являются независимые входные воздействия: усилие нажатия Fн, скорость относительного перемещения в контактной паре vk, электрическое напряжение и, температура внешней среды 0Ср и начальная температура КЭ.
Усилие нажатия зададим как постоянную величину или как функцию от времени, которую определим приближенно или по результатам моделирования механики контактного взаимодействия.
Скорость перемещения зададим как отношение интервала времени перехода к интервалу расчётной сетки. При этом на каждом шаге перемещения производим переприсвоение значений в матрице превышений подвижного микрорельефа:
где п---размер квадратной матрицы превышений подвижного рельефа.
Напряжение прикладываем к цепи, состоящей из сопротивления тела щётки и параллельного соединения сопротивлений контактных кластеров. Определяем ток и электрическую мощность для каждого КЭ.
Модификация компонент векторов состояния контактных элементов
Данный блок операций является центральным, он определяет основные динамические процессы и содержит главные вычислительные процедуры. Последовательность операций такова:
1. В матрице переходного слоя выделяем кластеры непосредственной и плёночной проводимости (а-кластеры и в-кластеры). Для этого используем процедуру распознавания, транслирующую элементы каждого выделенного кластера в дополнительную матрицу.
2. Подсчитывается количество элементов, составляющих кластер N и вычисляется характеристический радиус кластера г^ как радиус гирации Rg
Формирование и приложение силового вектора управляющих и возмущающих воздействий
hi у+1 = hi у при і = 1... п, } = 1... п - 1, hi,1 = hi,п пРи і = 1 ••• п ,
(3)
(4)
1 N 1 N
X = — У л., 7=—У у.. (6)
с N1=1 ' с N “Т
3. Определяем фрактальную размерность кластера.
4. Находим сопротивление каждого кластера как сумму сопротивлений стягивания верхнего и нижнего полупространств (СК), сопротивлений окисной
плёнки и продуктов износа ^ + Кпи + ^ + Кн ,
где
р
п • = _ П " = (7)
ст. л ^ 5 ст. л т-ч ^ )
4и • г 4и • г
клг клг .
5. Рассчитываем переходное сопротивление контакта как параллельное соединение сопротивлений всех кластеров
П = (У П™-)"'. (8)
6. Вычисляем суммарный ток контакта — сопротивление тела щётки).
I = и /(Пщ + П). (9)
7. Находим тепловую мощность, выделяющуюся в каждом кластере. Мощность содержит электрическую и механическую составляющие.
Электрическая составляющая:
Р . = Ли~, (10)
ЭЛ1 т-\ ^ ^ '
2 Пкл!
ли = и -1 • Ящ, (11)
где Ди — падение напряжения на контактном переходе; и — напряжение, приложенное к контактной паре; I — ток СК; Rщ — сопротивление щётки. Механическая мощность сил трения:
Рмех! = Ск! -лу1 • ктр • V , (12)
где Ск — жёсткость КЭ, Н/м; ДY — сжатие КЭ, м; ктр — коэффициент трения; V — скорость смещения микрорельефов, м/с. Тогда суммарная мощность в КЭ:
Р . = Р .+ Р (13)
кэ1 эл- мех! • V1
8. Определяем приращение температуры КЭ за один квант модельного времени. Т. к. выделяющаяся мощность затрачивается на нагрев объёма КЭ до температуры 0 за время Д^ то в соответствии с методом электротермической аналогии, составим эквивалентную электрическую схему замещения. Источник тепловой мощности Ркэ представим источником тока 1э, а теплоёмкость КЭ Скэ — ёмкостью Сэ, при этом, значению напряжения иэ поставим в соответствие температуру 0 (рис. 1).
із; (Ркэ)
а
Сз; (Скэ)
из; (Є)
_Х_
Рис. 1. Эквивалентная схема замещения
Запишем в операторной форме:
и, (р)-і, ■ хс (р) -
і
Р ■ Сэ
(14)
гДе X (Р):
Р ■ Сэ
— реактивное сопротивление конденсатора Сэ
В соответствии с правилом прямой аналогии:
Р
в( р )= кэ
Скэ ■ Р
Скэ Скэ ■ ткэ ■
(15)
где скэ — удельная теплоемкость; ткэ = Ах2 • Аh • ркэ — масса КЭ; Ah — высота КЭ; р — средняя плотность материала КЭ.
После перехода во временную область, получим:
в{( )=
р„
і.
Скэ •Ах 2 ^ • Рк
Представим последнее выражение в разностях:
(16)
Ав
■ Аі.
с •Ах •Ah • р
кэ кэ
За один шаг моделирования температура КЭ станет:
(17)
ві+1 - ві + Ав .
(18)
9. Температура КЭ изменяется за счёт теплообмена с контактирующими телами и соседними КЭ. Поэтому, после вычисления приращения температуры, вызванного тепловыделением в КЭ, рассчитываем распределение температуры в трёхмерной системе узлов за тот же интервал времени (рис. 2) явным методом [4]. С учётом неравномерного распределения теплопроводности и теплоёмкости по структуре переходного слоя, запишем:
1
9і, і ,5 (ґ + Лґ ) -
9і, і ,5 (ґ ) + Кі+1, і ,5 (ґ ) ' Л9і+1, і ,5 (ґ ) + К, і-1,5 (ґ ) ' Л9і-1, і ,5 (ґ ) +
+ Кі, і +1,5
(ґ) -Л9
, і +1,5 (Ґ) + К,] -1,5 (Ґ) ' Л9і,] -1,5 (Ґ) +
+ К,] ,5+1 (Ґ) ' Л9і,] ,5+1 (Ґ) + Кі,] ,5-1 (Ґ) ' Л9і,] ,5-1 (Ґ)
кц,в+1, Сі],8+1, ^^1,],з+1 (") ki,j-l,s, Ci,j-l,s, 0у-1,8
кі-1 J,s, ci-1,j,s, 01-1,]
кІ+1 ,1,85 СІ+1,1,85 01+1,],;
kІ,j,S, СЦ,8, 0І]
кІ,]+1,8, 0І,І + 1,8 кІ,І,8-1 , СІ,І,8-1 , 0і,],8-1
Рис. 2. К расчёту нестационарного температурного поля СК
где в9 ]5(ґ) — температура узла расчетной сетки с координатами і, і, 8;
Л9і+1,і ,5 (Ґ) - 9і +1, і ,5 (Ґ) - 9і, і ,5 (Ґ) , Л9і-1, і ,5 (Ґ) - 9і-1, / ,5 (Ґ) - 9і, / ,5 (Ґ) , и т. Д. — ра3ности
температур расчётного
Ц
К
сі ,і ,ЛіЛх Ч Рі
узла и примыкающих к нему узлов; • Дt — число Фурье для текущего расчётного узла;
kiJ■sJ, \С1 я/, \Pijsj — средние значения теплопроводности, удельной теплоёмкости и плотности, вычисленные по узлам фрагмента (рис. 2).
Вычислительная процедура имеет устойчивое решение, если число Фурье < 1/6 [4]. Данный показатель имеет наибольшее значение для зон непосредственной проводимости, поэтому выбор шага Дt и пространственную дискретизацию (Дх) производим по КЭ, относящимся к указанным зонам.
Приращение температуры КЭ от электромеханического тепловыделения, и ее изменение за счет распределения тепла по трехмерной сетке, не объединяем в общую процедуру, т. к. при этом нарушается принцип одновременности тепловыделения в КЭ. Расчёт нестационарного температурного поля СК производим для каждого Дt в два этапа (п.7, п.8), а для анализа используем только заключительные значения температур КЭ второго этапа.
Для упрощения алгоритма и увеличения быстродействия введём допущение о том, что температура близлежащих к контакту слоёв КД определяется как полусумма средней температуры всех расчётных КЭ переходного слоя и температур контактирующих тел. Также пренебрегаем конвективным теплообменом на границах контактного слоя. Для этого задаём нулевые значения теплопроводности связей граничных КЭ с внешней средой.
10. На каждом шаге расчёта температурного поля для всех КЭ проверяем условие превышения критических уровней температур материалов: температу-
8
8
ры рекристаллизации 0р; температуры плавления 0тп; температуры испарения 0и.
В зависимости от того, в какой интервал попадает вычисленная температура конкретного КЭ, принимаем решение об изменении его характеристик и параметров.
1) До температуры 0р механическую жёсткость КЭ считаем постоянной. В интервале 0е [0р, 0тп] изменяем жёсткость КЭ. При достижении температуры 0и происходит испарение наименее теплостойкого материала КЭ, и увеличение расстояния между контактирующими поверхностями в данной точке. Для этого вносится коррекция в матрицы микрорельефов. Сближение в КЭ изменяем на величину Дк В результате при неизменном усилии нажатия Fн изменяется эквивалентная жёсткость СК Скэ и сближение контактирующих поверхностей.
ДГ = Сл., где Скэ = 1 Ск,. (20)
Скэ
2) При увеличении температуры повышается проводимость поверхностных плёнок, а при достижении КЭ температуры плавления 0тп наступает тепловой пробой, который определяем как переход КЭ в состояние непосредственной проводимости. При моделировании используем эмпирические зависимости.
11. Фриттинг оксидных плёнок интерпретируем как переход КЭ в состояние непосредственной проводимости при достижении электрической напряжённостью порогового значения Еф, зависящего от напряжения на контактном переходе, и толщины пленки
Е = Ди, (21)
а
где а — толщина оксидной пленки. Если Е > Еф, то КЭ переходит в состояние
непосредственной проводимости.
12. Учитываем зависимости от температуры, удельного сопротивления, удельной теплоёмкости, теплопроводности и мощности силы трения.
Рэл = РэлС •(1 + «-ДТ), (22)
с = сС -(1 + ас -ДТ), (23)
к = к с -(1 + ак -ДТ), (24)
Р (т) = Р
мех\ ! м
мех С
с с ДТ л п ^
1 -
V V Т* - Т0 у )
(25)
где рэл0 — удельное электрическое сопротивление при начальной температуре; а — температурный коэффициент сопротивления; АТ — разница между текущей и начальной температурами; сС и кС — удельная теплоёмкость и теплопроводность при начальной температуре; ас и ак — температурные коэффициенты удельной теплоёмкости и теплопроводности. Коэффициент ак — может быть, как положительным, так и отрицательным, в зависимости от мате-
риалов СК. Рмех 0 — мощность силы трения при начальной температуре; T — текущая температура; n и T* — параметры, подбираемые из условия наилучшей аппроксимации эмпирической зависимости Рмех (T) в интересуемом интервале температур. В большинстве случаев можно ограничиться линейной моделью n = 1, а в качестве T* использовать температуру плавления наиболее легкоплавкого материала СК.
Выводы
Таким образом, разработаны основные алгоритмические действия и математические зависимости, необходимые для процесса имитационного моделирования скользящего электрического контакта: а) задание начальных значений компонент векторов состояния контактного элемента; б) формирование и приложение силового вектора управляющих и возмущающих воздействий; в) модификация компонент векторов состояния контактного элемента; г) вычисление и представление интегральных характеристик. Дальнейшие исследования предполагают программную реализацию описанных алгоритмов и проведение вычислительных экспериментов с последующей их верификацией на практике.
Литература
1. Плохов И. В. Комплексная диагностика и прогнозирование технического состояния узлов скользящего токосъёма турбогенераторов. Диссертация д-ра. техн. наук. С.-Петербург, СПбГПУ, 2002.
2. Плохов И. В. Модель динамики токопередачи через скользящий контакт Электротехника. № 2. 2005. С. 28-33.
3. Voss R. F. Random fractal forgeries / in: Fundamental Algorithms in Computer Graphics (ed R.A.Earnshaw, Springer-Verlag, Berlin, pp. 805-835. 1985. Цветные фотографии на С. 13-16.
4. Крейт Ф., Блэк У. Основы теплопередачи / Пер. с англ. М.: Мир. 1983.
Об авторах
Плохов Игорь Владимирович — заведующий кафедрой «Электропривод и системы автоматизации» ФГБОУ ВПО ПсковГУ, д-р техн. наук, профессор.
E-mail: igor_plohov@list.ru
Ильин Александр Викторович — старший преподаватель кафедры «Электропривод и системы автоматизации» ФГБОУ ВПО ПсковГУ, аспирант.
E-mail: al.ilyin@yandex.ru
Козырева Оксана Игоревна — старший лаборант кафедры «Электропривод и системы автоматизации» ФГБОУ ВПО ПсковГУ, аспирант.
E-mail: ks_33n@rambler.ru
I. V. Plohov, A. V. Iliin, O. I. Kozyreva
STRUCTURE AND ALGORITHMIC STEPS SIMULATION OF THE DYNAMICS ELECTRO FRICTION ENGAGEMENT
Electro friction engagement in electric machines is considered. Basic algorithmic steps and mathematical relationships carried out in the modeling process are offered.
Keywords: electro friction engagement, slide electric contact, component of state vector, dynamic simulator, fractal clusters.
