Труды МАИ. Выпуск № 83
www.mai.ru/science/trudy/
УДК 621.8: 629.78
Математические модели гистерезиса, описывающие деформации механизмов для стыковки космических аппаратов
Яскевич А. В.
Ракетно-космическая корпорация «Энергия», ул. Ленина, 4а, Королев, Московская область, 141070,Россия e-mail: Andrey. Yaskevich@rsce.ru
Аннотация
Механическое соединение космических аппаратов при стыковке обеспечивается специализированными механизмами стыковочного агрегата. Пружины, звенья и отдельные части этих механизмов могут деформироваться с гистерезисом, обусловленным потерями энергии вследствие внутреннего трения. Это явление означает, что силы или моменты реакций зависят не только от текущих деформаций, но и от их величин в предшествующие моменты времени. Математические модели стыковки должны учитывать упругие конструктивные деформации с гистерезисом звеньев и отдельных частей механизмов для корректного расчета нагрузок на стыкуемые объекты. В данной работе рассматриваются простые модели таких деформаций, основанные на использовании экспериментально полученных данных.
Ключевые слова: деформации, математическая модель гистерезиса, стыковка
Введение
Стыковка представляет собой управляемый механический процесс соединения космических аппаратов на орбите с использованием стыковочных агрегатов, конструкция которых впервые рассмотрена в [1]. Три типа механизмов этих агрегатов - механизм сцепки, стыковочный механизм и механизм жесткого соединения обеспечивают выполнение этого процесса. Результаты математического моделирования стыковки должны соответствовать данным динамических испытаний. Практика показывает, что это возможно только при учете деформаций с гистерезисом упомянутых выше механизмов.
Теоретические модели гистерезиса описывают причины возникновения этого явления на основе учета физических свойств различных сплавов и материалов [2-7] или особенностей функционирования биологических объектов [8], органов человека [9] и элементов технических устройств [10].
В механизмах стыковочных агрегатов гистерезис обусловлен деформацией передач вращения с люфтами и сложных сборок, состоящих из большого числа неоднородных по материалу и размерам деталей. Разработка теоретической модели потерь механической энергии в таких механизмах практически невозможна. В данной работе модель гистерезиса описывает не причинно-следственные связи этого явления, а результаты его наблюдения при статических испытаниях. Она предназначена для оценки влияния потерь механической энергии в механизмах на динамический процесс стыковки космических аппаратов. Для этого она должна
обеспечивать значительно более точное описание ветвей гистерезиса, чем то, что обычно используется в теории систем автоматического регулирования [11].
Динамические нагрузки при стыковке определяются в основном стыковочным механизмом. Он обеспечивает устранение рассогласований между направляющими элементами для достижения сцепки, поглощение кинетической энергии относительного движения космических аппаратов и ограничение их максимальных относительных перемещений, выравнивание и совмещение стыковочных плоскостей. Поэтому в данной работе модели деформаций с гистерезисом рассматриваются на примере стыковочного механизма. Модели для других типов могут быть получены аналогично.
Влияние гистерезиса наиболее наглядно проявляется при функционировании стыковочного механизма простейшей системы стыковки типа «штырь-конус».
Деформации стыковочного механизма
Относительная простота стыковочного механизма системы «штырь-конус» состоит в том, что его звенья, обеспечивающие продольные осевые и боковые/угловые перемещения контактирующих буферных элементов могут двигаться независимо. Это делает более легким сопоставление упругих свойств отдельных звеньев и механизма в целом. Упрощенная кинематическая схема механизма представлена на рис. 1.
Рис. 1
Буферным звеном является штырь с полусферической головкой. Он может контактировать с приемным конусом и гнездом пассивного агрегата. На
поверхности головки расположены контактные датчики, а по периметру - четыре
защелки и два датчика сцепки.
Продольная компонента контактной силы, действующей на головку, вызывает
поступательное движение штыря относительно оси подвижного корпуса, связанного с основанием универсальным шарниром. Это перемещение преобразуется шарико-винтовым шарниром во вращательное движение устройств осевого амортизатора -пружинного механизма ПМх, фрикционного тормоза ФрТ, осевого магнитного демпфера Дх.
Боковая компонента контактной силы, действующей на головку, изгибает штырь и разворачивает подвижный корпус. Два поступательных пружинных
механизма ПМ1 и ПМ2 и два боковых магнитных демпфера БД1, БД2
противодействуют этому развороту. ПМ1 и ПМ2 располагаются во взаимно перпендикулярных плоскостях, проходящих через продольную ось стыковочного механизма. Кинематические цепи магнитных демпферов также располагаются во взаимно перпендикулярных плоскостях, развернутых на 45 градусов относительно плоскостей пружинных механизмов.
После поглощения основной части энергии относительного движения стыкуемых аппаратов электродвигатель ЭД привода стыковочного механизма втягивает штырь в подвижный корпус, обеспечивая постепенное угловое выравнивание агрегатов по тангажу и рысканию. При этом боковой амортизатор демпфирует угловые колебания стыкующегося объекта относительно шарнира стыковочного механизма. Выравнивание агрегатов по крену (вращение относительно продольной оси) осуществляется при контактах защелок головки с направляющими пазами приемного гнезда. Демпфирование такого движения обеспечивается конструктивной жесткостью стыковочного механизма. До стягивания стыковочных агрегатов вращение кинематической цепи между фрикционной муфтой и электродвигателем Mot блокируется стопорной муфтой СМ.
Деформации упомянутых выше элементов конструкции стыковочного механизма имеют следующие особенности.
1. Деформации имеют одну или две степени свободы относительно продольной оси деформируемого элемента. Одну степень свободы имеют деформации сжатия-растяжения и кручения. Двумя независимыми степенями свободы может быть приближенно описан пространственный изгиб штыря.
Сферическая форма головки позволяет не учитывать распределение угла прогиба упругого стержня по его длине и использовать только два угла по тангажу и рысканию, определяющие отклонение головки от продольной оси механизма.
2. Все упругие деформации обладают гистерезисом, то есть сила или момент реакции неоднозначно зависят от величины деформации и определяются значением последней не только в текущий, но и в предшествующие моменты времени.
Полученные экспериментально характеристики жесткости всего стыковочного механизма системы типа «штырь-конус» при его угловых деформациях по крену (вращение головки штанги относительно продольной оси) и осевого редуктора с пружинным механизмом при передаче движения на фрикционный тормоз имеют вид, приведенный на рис. 2а и 2б.
а б
Рис. 2 6
Прямая и обратная ветви гистерезиса на этих графиках отмечены буквами Б и В. Направления изменения деформаций показаны стрелками. Переход между прямой и обратной ветвями гистерезиса при смене знака скорости деформации не является мгновенным, и осуществляется по переходной ветви, параметры которой зависят от величины деформации.
Сила реакции ^ на боковое смещение головки определяется
конструктивной изгибной жесткостью штанги и параметрами 2-ступенчатых пружинных механизмов БМ1 и БМ2 (рис. 3).
а б
Рис. 3
Идеализированная, расчетная характеристика жесткости бокового амортизатора стыковочного механизма, полученная без учета деформации штанги и гистерезиса, имеет вид, показанный на рис. 3а, а экспериментальная с учетом этих
факторов - на рис. 3б. При этом к - му линейному участку расчетной характеристики (рис. 3а), соответствуют два линейных участка прямой кР и и обратной кв ветви характеристики с гистерезисом.
Рассмотренные характеристики жесткости с гистерезисом типа «двойная петля» (по классификации [3]) отражают потери механической энергии за счет сил внутреннего конструкционного трения. Учет этих потерь позволяет адекватно отразить форму и периоды следования контактных сил на всех этапах стыковки - от первого контакта стыковочных агрегатов, до завершения их стягивания. В результате корректно определяются нагрузки, действующие на стыкуемые космические аппараты.
Кусочно-линейная модель гистерезиса с постоянными параметрами
Математический анализ стыковки основан на численном интегрировании уравнений движения. Это предъявляет высокие требования к вычислительной эффективности используемых моделей гистерезиса и обуславливает использование кусочно-линейной аппроксимации экспериментальных данных. Такая модель предназначена для описания деформаций, как отдельных звеньев, кинематических цепей так и конструктивных сборок. Ее параметры не изменяются при движении механизма. Эта модель основана на следующих упрощающих допущениях.
1. Все ветви, описывающие характеристики жесткости с гистерезисом, с любой необходимой точностью представляются в виде совокупности отрезков
прямых линий, которые нумеруются по мере возрастания величины деформации от нуля до границы допустимого рабочего диапазона.
2. Силы или моменты реакции примерно равны нулю при малых деформациях, поэтому графики жесткости пересекают начало системы координат. Люфт моделируется отрезком с нулевым номером. Этот отрезок пересекает начало системы координат, имеет нулевой угол наклона и нулевую величину начальной обобщенной силы. В частном случае люфт может быть равен нулю.
3. Прямые и обратные ветви гистерезиса совпадают на первом линейном отрезке, следующем за люфтом, и на последнем, который соответствует деформациям, выходящим за границы рабочего диапазона модели, и имеет максимальный угол наклона (коэффициент жесткости).
4. Переходная ветвь между прямой и обратной ветвями гистерезиса определяется как отрезок прямой линии с коэффициентом наклона таким же, как и на последнем участке с максимальной жесткостью за границей рабочего диапазона модели.
5. Последний по номеру участок обратной ветви определяется линейным отрезком с максимальным углом наклона, как и переходная ветвь.
Общий для обеих ветвей первый линейный отрезок имеет малый размер, и незначительно влияет на точность представления гистерезиса. Он используется для исключения влияния погрешностей численных методов интегрирования на логику изменения состояния модели гистерезиса при малых деформациях Описание переходной ветви линейной функцией (допущение 4) примерно соответствует
экспериментальным данным. Последнее допущение гарантирует переход из состояния за пределами рабочего диапазона модели гистерезиса на последний участок обратной ветви при изменении знака скорости деформации. Необходимая точность представления экспериментальных данных легко обеспечивается выбором числа и параметров линейных участков. Используемые упрощения обеспечивают высокую вычислительную эффективность моделей при практически незначительном снижении их точности. Далее для краткости рассматриваются модели гистерезиса только для положительных деформаций.
Модели гистерезиса с постоянными параметрами для линейных или угловых упругих деформаций различных элементов конструкции могут быть рассчитаны по одному алгоритму, использующему зависимость обобщенной силы реакции Q от обобщенного перемещения (деформации) д и скорости ¿д. В системах с гистерезисом текущий отклик зависит от их текущего состояния. В предлагаемой кусочно-линейной модели ветви гистерезиса нумеруются, и состояние деформируемой системы определяется этим номером. В частности, общему линейному участку прямой и обратной ветвей с малыми деформациями соответствует номер и состояние 0, прямой и обратной ветвям - номера и состояния 1 и 2 соответственно. Линейному участку перехода между прямой и обратной положительными ветвями соответствует номер и состояние 3. Линейному участку, общему для двух ветвей за пределами рабочего диапазона, соответствует номер и состояние 10. Для отрицательных деформаций номера ветвей и состояния модели имеют значения -1, -2, -3 и -10. Каждый / - й линейный участок прямой
положительной ветви характеризуется значениями начальной деформации (г), обобщенной силы QFP (г) и коэффициентом жесткости кРР (г). Для г - го участка обратной положительной ветви аналогичные параметры обозначаются как двр (г),
Qвp (1) и квр (г), г = 1, N (рис. 4).
Рис. 4
Участки прямой и обратной отрицательных ветвей имеют параметры д ш (г),
QFN(г), kFN(г) и двЫ(г), QвN(г), квЫ(г), г = 1, Ыы соответственно. Здесь , Nлr - число точек (концов отрезков), используемых в кусочно-линейной аппроксимации при положительных и отрицательных перемещениях.
Условия изменения состояния модели, то есть перехода между прямой и
Состояние 0 - отсутствие гистерезиса
обратной положительными ветвями гистерезиса, представлены в таблице 1.
11
Таблица 1.
Текущее состояние Условия изменения состояния Новое состояние
Любое ¿вы(2) = ¿ры(2) < д < дрр(2) = ¿вр(2) 0
0 д > 0, д > ¿РР(2) 1
¿1 < ^ д < ¿ш(2) -1
1 д > 0, ¿рр(2) < д < ¿рр (ы) 1
д > 0, д > ¿рр(n) 10
д < 0, ¿рр(2) < д < ¿рр (ы) 3*)
10 д > 0, д > ¿рр(n) 10
д < 0 2
2 1 < 0, д >¿вр(2) 2
1 > 0, д > ¿рр(N) 10
1 > 0, д < ¿рр(Ы) 3*)
3 д зв < д < д зе 3
д > ^ д > д зе 1
1 < ^ д < д зв 2
^т-1-
) Определяются параметры переходного отрезка
Аналогично определяется смена состояний модели при отрицательных обобщенных перемещениях.
Параметры переходного отрезка вычисляются при переходе в состояние 3 из состояний 1 или 2. Этот отрезок прямой имеет единственный постоянный параметр - коэффициент жесткости к3 = кР тах. В начале очередного перехода одна из его
точек является текущей точкой на отрезке прямой или обратной ветви. Другая является точкой пересечения переходного отрезка с каким-либо отрезком соответственно обратной или прямой ветви и является решением системы двух линейных уравнений, в которых первое является уравнением переходного отрезка, а второе - уравнением одного из отрезков той ветви, на которую осуществляется переход.
При переходе от прямой ветви к обратной, то есть из состояния 1 в состояние 2 (рис. 5а), значение деформации в процессе перехода убывает. Поэтому точка (хе з ,УЕ3) на отрезке прямой ветви, соответствующая текущей деформации,
является также конечной точкой переходного отрезка. Ее абсцисса хе 3
соответствует значению деформации на момент смены знака скорости, а ордината У определяется параметрами отрезка прямой ветви на этот момент. Из этих
координат определяется вторая точка ( х 0 3 = хе з - уЕ 3 / к3, у0 3 = 0) линии
переходного отрезка, расположенная на оси абсцисс, и параметры этой прямой
А3 = т 3 = У Е ,3 , В 3 =— 13 =-уЕ, 3/к 3 , С 3 =— А3 Х 0,3 .
а
б
Рис. 5
Далее для каждого очередного г — го отрезка, А 2 х + В2у + С 2 = 0 обратной
ветви, где А 2 = т 2, г = У е,2, г — У В, 2, г, В 2 =— 1 2, г =— ХЕ,2, г + ХВ,2, г, С2 = — А 2 ХВ,2, г — В 2 уВ,2, г
решается система уравнений
А 2 Х + В 2 У = — С2
А 3 Х + В 3 У = — С3
и определяется точка его пересечения ( х , У ) с линией переходного отрезка. При
выполнении условий
хв,2,г < хс < хе,2, г, УВ,2,г < Ус < Уе,2, г - г = 1, ыр
она является начальной точкой переходного отрезка, то есть (хрз = хс , урз = ус)
При переходе от обратной ветви к прямой, то есть из состояния 2 в состояние 1 (рис. 5б), значение деформации возрастает в процессе перехода. Поэтому точка (^з, У в з) на отрезке обратной ветви, соответствующая текущей деформации,
является начальной точкой переходного отрезка. Ее абсцисса хв 3 соответствует
значению деформации на момент смены знака скорости, а ордината ув 3
определяется параметрами отрезка обратной ветви на этот момент. Из этих координат определяется вторая точка (х03 = хвз -увз/к3, у03 = 0) линии
переходного отрезка, расположенная на оси абсцисс, и параметры этой прямой
Аз = т з = ув,з, вз =- Iз =-у в,з/кз, с3 =-А3 х0,3. Далее для каждого очередного г - го отрезка А1 х + в1 у+С1 = 0 прямой ветви,
где а1 = т 1, г = УЕХ г -Ув,1, г> в1 =- 11, г =- ХК,\, < + Х ' С1 =-А1 Хв,1, г -в1 У В,\1 решается
система уравнений
Ах+в у = - С
А 3 х + в 3 У = -Сз
и определяется точка его пересечения ( х , У ) с линией переходного отрезка. При выполнении условий
Хв,1,г < ХС < ХЕ,1, г и Ув,1,г < Ус < УЕ,1, г > 1 = 1, NN она является конечной точкой переходного отрезка, то есть (Х£3 = Хс, уЕЪ = ус).
Аналогично определяются параметры переходного отрезка между ветвями при отрицательных обобщенных перемещениях.
Таблица 1 определяет допустимую последовательность смены состояния модели гистерезиса. Нарушение этой последовательности вызывает остановку процесса моделирования с выдачей диагностического сообщения. Такое событие может быть вызвано влиянием вычислительной погрешности процесса численного интегрирования при очень малых деформациях. Для устранения этого влияния в модели используется первый линейный отрезок малого размера, общий для обеих ветвей гистерезиса. При нарушении логики изменения состояния модели необходимо либо увеличить размер этого участка, либо уменьшить шаг интегрирования.
Модель гистерезиса с постоянными параметрами при моделировании стыковочного механизма системы типа «штырь-конус» описывает деформации штанги и боковых пружинных механизмов при боковых отклонениях головки, деформации сжатия и растяжения крышки, на которой установлен механизм, деформации по крену механизма в целом. При моделировании периферийного стыковочного механизма такая модель описывает деформации сжатия и изгиба кольца, растяжения направляющих выступов, а также редукторов, передающих движение от штанг к элементами системы амортизации - пружинам и демпферам.
Модель гистерезиса с переменными параметрами для передач с фрикционными тормозами
Вращательные фрикционные тормозы (ФрТ) используются в стыковочных механизмах для поглощения основной части энергии сближения космических
аппаратов. Передача осевого поступательного движения штанги к ФрТ в стыковочном механизме системы «штырь-конус» осуществляется через гайку шарико-винтового шарнира, редуктор и ПМх (рис. 1). Статическая характеристика такого амортизатора при проскальзывании ФрТ, то есть при повороте входного вала редуктора на угол Лу после окончания вращения ПМх, показана на рис. 6. Она может быть описана моделью гистерезиса с переменными параметрами.
Основные отличиям такой модели от описанной выше состоят в следующем. 1. Коэффициенты жесткости за пределами последних линейных отрезков аппроксимации (в состояниях 10 или -10) равны нулю, то есть
Дф
Мр
Рис. 6
кРР,N, N = кБР,N = N = NР , ^^,N = = N = ^^ •
2. Вращение входного вала фрикционной муфты являются необратимым, так как это устройство не накапливает, а рассеивает подводимую кинетическую энергию. Поэтому при смене знака скорости необходимо скорректировать аргументы табличного описания осевого амортизатора в соответствии с необратимым угловым перемещением.
Необратимое угловое перемещение амортизатора после завершения вращения муфты в положительном направлении (после поглощения энергии сближения космических аппаратов) равно
А^ = ^Р,шах (NP )
а после завершения вращения муфты в отрицательном направлении (после завершения стягивания) -
,тт -Фт (NN )
где ^Р,тах >ФРР (NP) и ^,тах (NN) - максимальные достигнутые углы
входного вала редуктора в положительном и отрицательном направлении.
Аргументы табличных описаний модели деформаций при смене знака скорости вращения муфты корректируются следующим образом
<рхр(1) = рхр (1) + > <Рхрл(*) = рхря (*) + > ' = 1 ^ ; ^( * ) = Фхя ( * ) + > ФхКК( * ) = ФхЖ ( * ) + А^х > ' = 1 NN .
Коррекция параметров модели гистерезиса может быть многократной. Она определяется числом срабатываний фрикционного тормоза. Смена состояний модели гистерезиса при положительных деформациях описывается таблицей 1.
Изменение параметров модели гистерезиса при отрицательных деформациях и срабатывании фрикционной муфты в обратном направлении описывается аналогично. Такой же метод может быть использован для моделирования деформаций редуктора, передающего вращение к 2-ступенчатой фрикционной муфте.
Модели деформаций упругих передач с фрикционными муфтами используются для описания процессов поглощения энергии в различных стыковочных механизмах центрального («штырь-конус») и периферийного типа.
Результаты моделирования динамики стыковки с учетом гистерезиса.
Описанные выше модели деформаций с гистерезисом используются для расчета реакций механизмов в математических моделях стыковки космических аппаратов. Учет гистерезиса позволяет получить более точное распределение во времени расчетных нагрузок на стыковочные интерфейсы и лучшее соответствие экспериментальным данным, получаемым на 6-степенном динамическом стенде для полунатурного моделирования стыковки. Такой стенд описан упрощенно в [1].
Для системы типа «штырь-конус» уравнения движения, описывающие динамику стыковки, описаны в [12], а алгоритмы определения сил контактного взаимодействия - в [13]. Графики полученных экспериментально и расчетных
контактных сил на начальной фазе стыковки показаны на рис. 7 и 8. Первый из них отражает процесс от момента первого контакта головки штанги и приемного конуса до завершения первого удара головки в дно приемного гнезда, а второй -последовательность ударов головки в дно и упоры приемного гнезда при поглощении кинетической энергии сближения активного корабля.
Рис. 7
Жесткость бокового амортизатора с гистерезисом, имеющим постоянные
параметры (рис. 3б), влияет на форму контактных сил при ударах головки в конус и
стенки его приемного гнезда (рис. 7). Момент сопротивления осевого амортизатора с гистерезисом, имеющим переменные параметры (рис. 2б и 6), а также деформации с гистерезисом переходного люка, на который установлен стыковочный механизм, определяют амплитуду и форму ударных импульсов осевых контактных сил, а также временные интервалы между ними (рис. 8).
Данные испытаний на 6-степенном стенде
-РХ-Ру-р2
1000.0
500.0 0.0 500.0 1000.0
г\
~~V/—
V
0.00 1.00 2.00 3.00 4.00 5.00 6.00 7.00 8.00 9.00 10.00
Время, сек
Интерфейсная контактная сила
Результат математического моделирования —-Р2Ру-?2Ъг
к ГС Активизация 1000.0 т-ФеТ
500.0
0.0 '
-500.0
1000.0
Сопротивление Дх и восстановление ПМх с гистерезисом
Удар головки в дно гнезда - обратный удар -сопротивление ПМх и Дх V деформации с ГИСТерезисом
0.0 1.0 2.0 3.0 4.0 5.0 6.0 7.0 8.0 9.0 10.0
Время, сек
Интерфейсная контактная сила
Результат математического моделирования -Р2Рх
кгс 0.0
200.0
400.0
600.0
800.0
1 0.11 се ^^ 1
-л 0.17 с ек >--<- к 1 -
1 1 1
5.56 5.59 5.63 5.66 5.70 5.73 5.76 5.80 5.83 5.87 5.90
Время, сек
Интерфейсная контактная сила
Рис. 8
Достаточно хорошее совпадение экспериментальных данных и результатов
моделирования обусловлено использованием моделей деформаций с гистерезисом.
21
Расчетные нагрузки на стыкуемые космические аппараты, в частности амплитуда контактных сил, число ударов и частота их следования были бы выше без учета этого фактора.
Заключение.
Разработанные кусочно-линейные модели гистерезиса при упругих деформациях обладают высокой вычислительной эффективностью и приемлемой инженерной точностью. Они могут быть использованы в различных технических приложениях, в которых необходимо не объяснять причину гистерезиса, а оценивать его влияние на динамический процесс.
Библиографический список
1. Сыромятников В.С. Стыковочные устройства космических аппаратов - М.: Машиностроение, 1984. - 216 с.
2. Augusto Visintin, Differential Models of Hysteresis (Applied Mathematical Sciences, vol. 111), Springer, 1994, 407 p.
3. R.V. Lapshin, Analytical model for the approximation of hysteresis loop and its application to scanning tunneling microscope. //Review of Scientific Instruments, volume 66, number 9, 1995, pp. 4718-4730.
4. Mielke, T. Roubicek. A Rate-Independent Model for Inelastic Behavior of Shape-Memory Alloys // Multiscale Modeling & Simulation, Volume 1, Issue 4, 2003, pp. 571-597.
5. I.D. Mayergoyz. Mathematical Models of Hysteresis and their Applications: Second Edition (Electromagnetism). Elsevier, 2003, 499p.
6. R.Smith. Smart material systems: model development, SIAM, 2005, 501 p.
7. Лукичев А.А., Ильина В.В. Простая математическая модель петли гистерезиса для нелинейных материалов. // Известия Самарского научного центра РАН. 2011. Т. 13. №4. С. 39-44.
8. G. A. Tompsett , L. Krogh , D. W. Griffin, W. C. Conner. Hysteresis and Scanning Behavior of Mesoporous Molecular Sieves. //Langmuir, 2005, 21 (18), pp. 82148225.
9. J.D. Escolar, A. Escolar. Lung hysteresis: a morphological view. //Histology and hystopathalogy, 2004, N. 19, pp. 159-166.
10. Пономарев Ю.К. Особенности гистерезиса пар сухого трения при круговых вращениях вибратора // Известия Самарского научного центра РАН. 2011. Т. 13. № 4(3). 1192-1195.
11. Попов Е.П. Теория нелинейных систем автоматического регулирования и управления: Учеб. Пособие. - М.: Наука, 1988. - 256 с.
12. Яскевич А.В. Комбинированные уравнения движения для описания динамики стыковки космических аппаратов с помощью системы «штырь-конус» // Известия РАН. Космические исследования. 2007. Т. 45. №4. С. 325-336.
13. A Yaskevich. Real time simulation of contact interaction during spacecraft docking and berthing. //Journal of Mechanical Engineering and Automation, vo.4, Number 1, January 2014, pp. 1-15.