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

МОСКОВСКИЙ ГОСУДАРСТВЕННЫЙ УНИВЕРСИТЕТ им. №В Ломоносова

ФАКУЛЬТЕТ ВЫЧИСЛИТЕЛЬНОЙ МАТЕМАТИКИ И КИБЕРНЕТИКИ Кафедра общей математики

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

Сазонов Василий Викторович

ПРИМЕНЕНИЕ ВЫЧИСЛИТЕЛЬНОЙ ГЕОМЕТРИИ В ЗАДАЧАХ МОДЕЛИРОВАНИЯ ВРАЩАТЕЛЬНОГО ДВИЖЕНИЯ КОСМИЧЕСКИХ АППАРАТОВ

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

АВТОРЕФЕРАТ диссертации на соискание ученой степени кандидата физико-математических наук

□03065345

Москва 2007

003065345

Работа выполнена на кафедре общей математики факультета вычислительной математики и кибернетики Московского государственного университета им М В Ломоносова.

Научный руководитель.

доктор физико-математических наук, профессор ТТТсттга Евгений Викторович

Официальные оппоненты:

доктор физико-математических наук Галактионов Владимир Александрович

кандидат физико-математических наук Бсрезин Сергей Борисович

Ведущая организация

Институт проблем механики РАН

Защита диссертации состоится 12 октября 2007 года в 14.30 на заседании диссертационного совета К 501 001.07 при Московском государственном университете им M Л Ломоносова по адресу: 119992, Москва, Ленинские горы, МГУ, 2-й учебный корпус, факультет вычислительной математики и кибернетики, яуд. 685

С диссертацией можно ознакомиться в библиотеке факультета ВМК МГУ Автореферат разослан « » сентября 2007 г

Ученый секретарь диссертационного совета, кандидат физико-математических наук, , . Л-~г— доцент

В.М Говоров

Общая характеристика работы

Актуальность темы. Для анализа орбитального и вращательного движения космических аппаратов (КА) необходимо уметь вычислять главный вектор и главный момент действующих на аппарат внешних сил Среди этих сил важное место занимают силы аэродинамического сопротивления и светового давления, действие которых на движение КА в ряде случаев является определяющим Вычисление главного вектора и главного момента сил светового давления (аэродинамического сопротивления) сводится к вычислению интегралов по части поверхности КА, освещаемой Солнцем (омываемой воздушным потоком) Обычно КА — объект сложной формы, и вычисление таких поверхностных интегралов является весьма трудоемким процессом Основная сложность здесь состоит в отыскании части поверхности, по которой следует выполнить интегрирование Существует ряд приемов упрощения расчетов подобного рода, например, приближение поверхности КА какой-либо простой поверхностью, но это увеличивает ошибку результата В такой ситуации применение геометрического моделирования и методов вычислительной геометрии для решения этой задачи представляется целесообразным

Первые работы по применению геометрического моделирования для расчета сил и моментов светового давления и аэродинамического сопротивления, действующих на КА сложной формы, были опубликованы в начале 1990-х годов В работах M M Комарова, Е Ю Зуевой и др поверхность аппарата аппроксимировалась многогранником с треугольными гранями Многогранник строился как оболочка, натянутая на набор контуров - замкнутых пространственных ломаных Поиск освещенных участков (для определенности будем рассматривать вычисление сил и моментов светового давления) использовал частичное упорядочение граней и анализ их взаимного затенения Этот алгоритм эффективен для моделей с небольшим числом граней С ростом числа граней время работы алгоритма быстро возрастает, поэтому его применение к сложным моделям нецелесообразно

В работах M M Комарова и И Ю Балабана для определения освещенных участков моделируемой поверхности применяется метод, основанный на алгоритме удаления невидимых линий Аппеля В принципе, этот метод позволяет работать с моделями большой сложности, но его конкретная реализация не учитывала особых и вырожденных случаев, возникающих при проектировании многогранной поверхности на плоскость вдоль светового потока Отсутствие обработки особенностей проектирования делает алгоритм неустойчивым и ограничивает его применимость в рассматриваемых задачах

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

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

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

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

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

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

Разработанный программный комплекс использовался для расчета аэродинамического момента при моделировании вращательного движения спутника Фотон-Ш и расчета момента сил светового давления, действующего на КА с солнечным парусом

Апробация работы. Основные результаты диссертации докладывались на XIV Международной конференции «Математика Компьютер Образование» (г Пущино, 22-27 января, 2007 г) и на семинаре по компьютерной графике под руководством д ф -м н , проф Е В Шикина, к ф -м н А В Бореско-ва и к ф -м н СБ Березина (факультет ВМК МГУ)

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

Структура и объем диссертации. Диссертация состоит из введения, четырех глав и списка литературы Текст работы изложен на 109 страницах, диссертация содержит 60 рисунков и 4 таблицы Список литературы состоит из 29 наименований

Содержание работы

Во введении обоснована актуальность выбранной темы, приведен краткий обзор полученных ранее результатов, сформулированы основные решаемые задачи

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

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

Пусть КА движется поступательно относительно атмосферы со скоростью у и молекулы воздуха, сталкиваясь с ним, испытывают абсолютно неупругий удар Построим плоскость П, перпендикулярную вектору V Главный вектор и главный момент действующих на аппарат сил аэродинамического сопротивления можно представить следующим образом

¥ = -Р80\у\у, М = р | у | (у х Р), Р = 2\т1Г(Т1,

1

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

центра проекции 1 -го участка относительно точки, являющейся проекцией полюса приведения аэродинамических сил

Введем определение Будем говорить, что точка р множества Q видима в направлении (1, если луч г = р + «(1, а > 0, не содержит никаких других точек множества Q кроме точки р Таким образом, сг,— это площадь проекции 1-го участка поверхности аппарата, видимого в направлении его скорости

V

Главный вектор и главный момент сил светового давления вычисляются по формулам

2

Г = -2р£ е,а, (в,:п, )2п, - -р^УР, (8, п, )п, - (1 - е, (я, п, ),

1 -Э I «

2

М = Б,а,(я,п,У (п, * г()■+ -г,<т,(*>Xй, * О +

I з I

+/**£( 1-г,)ег1(в,п,)г1

Здесь р — световое давление на абсолютно черную пластинку, расположенную перпендикулярно световому потоку, 8 — орт направления «КА - Солнце», <т, — площадь »-го участка поверхности аппарата (такой участок - грань или часть грани), который освещается Солнцем, п, — внешняя нормаль к г-ому участку, г, — радиус-вектор центра масс ;-го участка относительно полюса приведения сил светового давления, е1 и у, — коэффициенты зеркального и диффузного отражения 1-го участка. Считаем, что для каждой грани объекта коэффициенты зеркального и диффузного отражения постоянны

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

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

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

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

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

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

Ребро будем называть лицевым, если оно принадлежит двум лицевым граням, нелицевым, если оно принадлежит двум нелицевым граням, и граничным, если оно является общим ребром лицевой и нелицевой граней

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

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

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

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

Вершина называется особой, если ей инцидентно более двух граничных ребер или она является точкой сборки поверхности

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

Рис. 1. Сборка триангулированной поверхности

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

Анализ задачи в такой постановке показывает, что индекс затенения является кусочно-постоянной функцией, заданной на лицевых поверхностях объектов. Рассмотрим точку Р лицевой поверхности (р и предположим, что проекция граничного ребра другой лицевой поверхности у эту точку содержит (рис. 2). Ясно, что в ^-окрестности точки Р индекс затенения поверхности принимает различные значения и, следовательно, эта точка является точкой разрыва функции индекса затенения. Будем считать, что в самой точке Р эта функция принимает значение, равное максимуму значений индекса затенения из рассматриваемой е-окрестности.

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

Рис. 2. Индекс затенения в окрестности точки проекции 1раничного ребра

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

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

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

Количество контурных циклов, ограничивающих лицевую поверхность, зависит от формы поверхности рассматриваемого объекта и его положения в пространстве. У эллипсоида лицевая поверхность всегда будет ограничена одним контурным циклом (рис. За), а у тора существуют положения, при котором лицевая поверхность ограничена двумя контурными циклами (рис. 36).

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

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

• одна поверхность частично закрывает другую (рис. 4а);

• одна поверхность полностью закрывает другую (рис, 4в, 4е, 4к, 4л, 4м);

• одна поверхность полностью проектируется на другую (рис. 46, 4д);

• поверхности не перекрываются (рис. 4г, 4ж, 4и).

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

• проекции циклов не пересекаю тся, и проекции внутренних точек поверхности не пересекаются;

• проекции циклов не пересекаются, ко проекция одного из циклов лежит в области, ограниченной проекцией другого цикла;

• проекции циклов имеют только точки касания;

в проекции циклов имеют только точки пересечения;

• проекции циклов имеют и точки касания и точки пересечения.

б

Рис. 3. Возможные формы лицевых поверхностей

4).

О © О

к Л И

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

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

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

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

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

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

• точки ребер, соответствующих пересечению проекций граничных ребер,

• точки сборок поверхности

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

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

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

|свет

Рис 6 Стек граней

На рис. 6 изображен стек граней для точки Р. В стек входят грани сие в указанном порядке. Вместо грани е можно включить грань 4. В различных ситуациях удобно выбирать ту или иную грань из смежных.

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

Предлагаемый а,тгоритм поиска освещенных участков поверхностей сцены состоит из четырех этапов:

поиск лицевых поверхностей и ограничивающих их контурных

циклов; основная цель этого этапа — сформировать лицевые поверхности, найти граничные ребра и объединить их в контурные циклы; на э том этапе проводится анализ контурных циклов и выявление особых вершин; поверхности всех объектов сцены могут обрабатываться независимо;

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

3. отыскание освещенных частей граничных ребер, составляющих контурные циклы, и их проекций на грани лицевых поверхностей объектов;

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

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

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

Рис. 7. Результат работы алгоритма

Нетреуюльные грани для хранения и обработки неудобны, поэтому они разбиваются на треугольники

Учитывая оценки времени работы каждого этапа, получаем, что весь алгоритм требует 0{И + К\о%К + + + >1)М) операций, где N —

общее количество граней сцены, К — число граничных ребер сцены, Ка — количество пар пересекающихся граничных ребер, Ктг — среднее число граничных ребер, проекции которых пересекают проекцию одного граничного ребра, « — количество поверхностей сцены, й — глубина стека граней Параметр А зависит от формы поверхностей сцены и положения источника света, но не зависит от количества граней сцены Из оценки сложности следует, что сложность предлагаемого алгоритма линейно зависит от количества граней сцены и почти линейно от количества контурных ребер сцены Количество граничных ребер зависит от количества граней и положения источника света Оценить характер зависимости достаточно сложно В работе Кетт-нера и Велцла (1997 г) проводится анализ этой зависимости и приведены результаты экспериментов Если рассматривать многогранную аппроксимацию сферы, то с уменьшением хаусдорфова расстояния е число ребер в аппроксимации составляет а число граничных ребер не превышает О

при произвольном ортогональном проектировании

В ходе реализации предложенного алгоритма были произведены замеры времени обработки модели космического корабля Фотон-М2 для одного и того же положения источника света с различной степенью детализации Результаты замеров приведены в таблице 1

Таблица 1 Результаты измерений времени работы алгоритма

кол-во граней кол-во ребер кол-во конт ребер макс глуб стека граней время, сек

4528 6792 199 3 0 012019

12928 19392 360 3 0 023803

51712 77568 720 3 0 088785

186112 279168 1364 3 033119

На рис 8 и 9 приведены график зависимости времени работы алгоритма от числа граней сцены и график зависимости числа контурных ребер от общего числа ребер сцены

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

О 35 • 0 301

5 025-^ 0,20,8 015-

60000 100000 150000 Количество граней сцены

Рис 8 График зависимости времени работы алгоритма от количества граней сцены

\ 1200-О

I 900£

О 50000 100000 150000 200000 250000 ЗООООО КоЛ№ес«о ребер сиены

Рис 9 График зависимости количества контурных ребер от общего количества ребер сцены

Третья глава посвящена применению разработанного программного комплекса в задачах моделирования вращательного движения спутников «Фотон» Применяемая модель основана на динамических уравнения Эйлера движения твердого тела и учитывает действие на спутник четырех внешних механических моментов гравитационного момента, восстанавливающего аэродинамического момента, момента с постоянными компонентами в связанной со спутником системе координат и момента, возникающего при взаимодействии магнитного поля Земли с собственным магнитным моментом спутника Описанный в главе 2 программный комплекс использовался для расчета аэродинамического момента С этой целью была разработана геометрическая модель внешней оболочки спутника (рис 10) Именно в разработке этой модели и создании экономичных по вычислительным затратам процедур ее

использования заключалось участие диссертанта в создании модели движения спутника.

Рис. 10. Геометрическая модель спутника Фотон М-2

Детальный вид формул, задающих аэродинамический момент в уравнениях движения спутника, учитывает те факты, что внешняя оболочка имеет ось симметрии и что центр масс спутника лежит на этой оси. Поскольку координата центра масс на оси симметрии может меняться, она явдяется параметром модели. В такой ситуации аэродинамический момент можно задать двумя скалярными функциями угла между осью симметрии и вектором скорости спутника относительно атмосферы. Значения этих функций заранее насчитываются в узлах равномерной сетки на отрезке [0, 180°] и затем в расчетах при произвольных значениях угла интерполируются при помощи конечных рядов Фурье.

Тестирование модели выполнено посредством определения с ее привлечением фактического вращательного движения спутника Фотон М-2 (на-

ходился на орбите 31 05 - 16 06 2005) по данным бортовых измерений напряженности магнитного поля Земли Данные, полученные на интервале времени длиной несколько часов, обрабатывались совместно методом наименьших квадратов путем интегрирования уравнений движения спутника относительно центра масс В результате обработки оценивались начальные условия движения и параметры используемой математической модели, в том числе, координата центра масс спутника на оси симметрии его внешней оболочки Использование новой модели вращательного движения уменьшило среднюю квадратичную погрешность аппроксимации данных измерений на 12-15%

Четвертая глава посвящена сравнению двух способов расчета главного вектора и главного момента сил светового давления, действующих на космический аппарат с солнечным парусом Парус составлен из 8 одинаковых плоских лепестков в форме равнобочных трапеций, расположенных симметрично относительно некоторой оси, которую будем называть осью паруса (рис 11) Форма паруса инвариантна относительно поворотов вокруг этой оси на углы, кратные 45° Первый способ расчета основан на аналитических формулах, полученных без учета затенения одних лепестков паруса другими Эти формулы были получены ранее и существенно использовались при исследовании режима одноосной солнечной ориентации искусственного спутника Земли, несущего такой парус, и закона управления парусом, увеличивающим большую полуось орбиты спутника Второй способ использует детальную геометрическую модель паруса, которая построена при помощи разработанного программного комплекса и позволяет учесть взаимное затенение лепестков Оба способа исходят из того, что парус освещен плоскопараллельным пучком света и что фиксированная часть попавших на него фотонов отражается от его поверхности зеркально, а остальные полностью поглощаются

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

Рис 11 Геометрическая модель КА с солнечным парусом

Основные результаты

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

Программная реализация предложенного алгоритма удаления невидимых поверхностей использована для расчета аэродинамического момента в модели вращательного движения космических аппаратов «Фотон» Построенная модель использована при реконструкции вращательного движения спутника Фотон М-2 по данным бортовых измерений магнитного поля Земли Использование новой модели вращатель-

ного движения уменьшило среднюю квадратичную погрешность аппроксимации данных измерений на 12 - 15% 4 Разработана новая геометрическая модель космического аппарата с солнечным парусом, которая использовалась для вычисления действующих на аппарат главного вектора и главного момента сил светового давления Модель позволила сравнить два способа такого вычисления способа, основанного на упрощенных аналитических формулах, и способа использующего модель паруса Установлена область изменения параметров ориентации паруса относительно падающих на него солнечных лучей, в которой аналитические формулы пригодны для ис-потьзования

Публикации автора по теме диссертации

1 Сазонов Вас В Алгоритм отыскания освещенных участков многогранных поверхностей в плоскопараллельном световом потоке // Математическое моделирование, 2007, июнь, т 19, №6, с 16-30

2 Сазонов Вас В, Комаров ММ Алгоритм построения освещенных участков полигональной поверхности в плоскопараллельном световом потоке // Математика, Компьютер, Образование, сборник научных тезисов Москва, Ижевск РХД, 2007, январь, с 96

Ъ Сазонов Вас В , Сазонов В В Использование уточненной модели аэродинамического момента в задачах исследования вращательно движения спутников «Фотон» - М Институт прикладной математики им М В Келдыша РАН, 2007, июль, препр № 32, 22 с

4 Сазонов Вас В, Сазонов В В Расчет главного вектора и главного момента сил светового давления, действующих на космический аппарат с солнечным парусом - М Институт прикладной математики им М В Келдыша РАН, 2007, август, препр № 46, 25 с

Напечатано с готового оригинал-макета

Издательство ООО "МАКС Пресс" Лицеями ИД N00510 от 01 12 99 г Подписано к печати 04 09.2007 г Формат 60x90 1/16 Усллечл. 1,25 Тираж 85 экз. Заказ 410 Тел 939-3890 Телефакс 939-3891 119992, ГСП-2, Москва, Ленинские горы, МГУ им М В Ломоносов 2-й учебвыв корпус, 627 к.

СОДЕРЖАНИЕ

ВВЕДЕНИЕ

ГЛАВА 1. ВЫЧИСЛЕНИЕ ГЛАВНЫХ ВЕКТОРОВ И ГЛАВНЫХ МОМЕНТОВ СИЛ СВЕТОВОГО ДАВЛЕНИЯ И СИЛ АЭРОДИНАМИЧЕСКОГО СОПРОТИВЛЕНИЯ,

ДЕЙСТВУЮЩИХ НА ПОВЕРХНОСТЬ КОСМИЧЕСКОГО АППАРАТА

1.1. Определение моментов

1.2. Задача удаления невидимых поверхностей

1.2.1. Алгоритм удаления невидимых линий Аппеля

1.2.2. Алгоритм удаления невидимых линий Блинна

1.2.3. Алгоритм Хегера

1.2.4. Алгоритм Балабана

1.2.5. Алгоритм Гудрича 19 1.2.4. Выводы

ГЛАВА 2. АЛГОРИТМ ПОИСКА ОСВЕЩЕННЫХ УЧАСТКОВ МНОГОГРАННЫХ

ПОВЕРХНОСТЕЙ В ПЛОСКОПАРАЛЛЕЛЬНОМ СВЕТОВОМ ПОТОКЕ

2.1. Используемые определения и допущения

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

2.3. Способ решения задачи

2.4. Контурные циклы

2.4.1. Определение контурного цикла

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

2.4.3. Затененность контурных циклов

2.4.4. Проектирование освещенных участков контурных циклов на лицевые поверхности объекта

2.5. Описание алгоритма 31 2.5.1. Организация входных данных

2.5.2.1-й этап: отыскание лицевых поверхностей и контурных циклов

2.5.3.2-й этап: отыскание пересечений проекций контурных ребер

2.5.4.3-й этап: нахождение освещенных участков контурных циклов и их проекций на грани объекта

2.5.5.4-й этап: нахождение освещенных участков поверхностей

2.6. Оценка сложности алгоритма

2.9. Программная реализация алгоритма

2.10. Выводы

ГЛАВА 3. ИСПОЛЬЗОВАНИЕ УТОЧНЕННОЙ МОДЕЛИ АЭРОДИНАМИЧЕСКОГО МОМЕНТА В ЗАДАЧАХ ИССЛЕДОВАНИЯ ВРАЩАТЕЛЬНОГО ДВИЖЕНИЕ

СПУТНИКОВ ФОТОН

3.1. Математическая модель вращательного движения спутника

3.2. Методика определения движения спутника по измерениям бортовых магнитометров

3.3. Результаты определения вращательного движения

3.4. Выводы

ГЛАВА 4. ПРИМЕНЕНИЕ АЛГОРИТМА ПОИСКА ОСВЕЩЕННЫХ ОБЛАСТЕЙ МНОГОГРАННОЙ ПОВЕРХНОСТИ К МОДЕЛИРОВАНИЮ ДВИЖЕНИЯ

КОСМИЧЕСКОГО АППАРАТА С СОЛНЕЧНЫМ ПАРУСОМ

4.1. Космический аппарат с солнечным парусом

4.2. Форма паруса и аналитические выражения для главного вектора и главного момента действующих на него сил светового давления

43. Численный расчет главного вектора и главного момента действующих на КА сил светового давления

4.4. Выводы

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

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

В работе [2] рассматривается задача вычисления сил и моментов, действующих на космический аппарат, который использует в качестве движителя солнечный парус. Метод моделирования поверхности космического аппарата такой же, как и в работе [1], однако для определения освещенных участков моделируемой поверхности применяется алгоритм, основанный на алгоритме удаления невидимых линий Аппеля [3]. Этот метод позволяет работать с моделями большой сложности, однако разработанный алгоритм не учитывает особых и вырожденных случаев, которые связаны с проектированием. Отсутствие обработки особенностей проектирования делает алгоритм неустойчивым и ограничивает его применимость для задачи вычисления моментов.

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

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

Целью работы является разработка эффективного метода поиска освещенных участков многогранных поверхностей и применение этого метода для вычисления главного вектора и главного момента сил светового давления и аэродинамического сопротивления, действующих на космический аппарат произвольной формы с целью использования для моделирования вращательного движения аппарата «Фотон-М2» и аппарата, оснащенного солнечным парусом.

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

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

Разработанный программный комплекс использовался для моделирования вращательного движения космического аппарата «Фотон-М2». Модель основана на динамических уравнениях Эйлера движения твердого тела и учитывает действие на спутник четырех внешних механических моментов: гравитационного момента, восстанавливающего аэродинамического момента, момента с постоянными компонентами в связанной со спутником системе координат и момента, возникающего при взаимодействии магнитного поля Земли с собственным магнитным полем спутника. Использование новой геометрической модели для расчета аэродинамического момента повысило точность определения вращательного движения.

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

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

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

Третья глава посвящена применению разработанного программного комплекса для моделирования вращательного движения корабля «Фотон-М2», который совершал полет с 31.05 по 16.06.2005. Описана технология моделирования и приведены данные расчетов режимов вращательного движения аппарата на временном отрезке в 270 минут, на котором вращательное движение аппарата было неуправляемым. Результаты расчетов показывают, что использование новой геометрической модели поверхности космического аппарата позволило в некоторых случаях уменьшить ошибку определения параметров движения на 15%.

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

Основные результаты опубликованы в четырех печатных работах [5], [6], [7] и [8].

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

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

3. Указанный программный комплекс использован для моделирования вращательного движения космических аппаратов «Фотон». При помощи построенной модели произведена оценка параметров вращательного движения аппарата во время полета с 31.05 по 16.06.2005. Произведено сравнение результатов оценки с результатами, полученными ранее, которое показывает, что использование новой геометрической модели поверхности космического аппарата позволило в некоторых случаях уменьшить ошибку определения параметров движения на 15%.

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

Заключение

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

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

1. Е.Ю. Зуева, М.М. Комаров, В.В.Сазонов. Представление аэродинамического момента в задачах моделирования вращательного движения искусственных спутников Земли // Космические исследования, 1992, т.ЗО, № 6, с. 771-779.

2. И.Ю. Балабан, М.М. Комаров, В.В. Сазонов. Определение освещенных участков поверхности космических объектов при расчете сил и моментов светового давления. -М.: Институт прикладной математики им. М.В.Келдыша РАН, 1995, препр.№57,28 с.

3. A. Appel, The Notion of Quantitative Invisibility and Machine Rendering of Solids // Proc. of ACM National Conference. Thompson Book, 1967, pp. 387-393.

4. Вас.В. Сазонов. Алгоритм отыскания освещенных участков многогранных поверхностей в плоскопараллельном световом потоке // Математическое моделирование, 2007, июнь, т. 19, №6, с. 16-30.

5. Вас.В. Сазонов, М.М. Комаров. Алгоритм построения освещенных участков полигональной поверхности в плоскопараллельном световом потоке // Математика, Компьютер, Образование, сборник научных тезисов. Москва, Ижевск: РХД, 2007, январь, с. 96.

6. Вас. В. Сазонов, В.В. Сазонов. Использование уточненной модели аэродинамического момента в задачах исследования вращательно движения спутников Фотон. М.: Институт прикладной математики им. М.В.Келдыша РАН, 2007, июль, препр.№32,22 с.

7. Вас. В. Сазонов, В.В. Сазонов. Расчет главного вектора и главного момента сил светового давления, действующих на космический аппарат с солнечным парусом. М.: Институт прикладной математики им. М.В.Келдыша РАН, 2007, август, препр.№46,25 с.

8. В.В. Белецкий. Движение искусственного спутника относительно центра масс. М.: Наука, 1965,416 с.

9. Е.В. Шикин, А.В. Боресков. Компьютерная графика. Полигональные модели. М.: ДИАЛОГ-МИФИ, 2000,464 с.

10. Е. Sutherland, R. F. Sproull, R. A. Schumacher. A characterization of ten hidden-surface algorithms // ACM Computing Surveys. New York: ACM Press, 1974, pp. 1-55.

11. M.T. Goodrich. A Polygonal Approach to Hidden-Line and Hidden-Surface Elimination // CVGIP: Graphical Models and Image Processing, 1992, vol. 54, №1, pp. 1-12.

12. J. Blinn. Fractional invisibility// IEEE CG&A, 1988, vol. 8, №6, pp. 77-84.14. http://wheger.tripod.com/vhl/vhl.htm.

13. И.Ю. Балабан. Применение методов символических возмущений для обработки вырожденных ситуаций в геометрических алгоритмах. М.: Институт прикладной математики им. М.В.Келдыша РАН, 1995, препр.№61,30 с.

14. L. Kettner, Е. Welzl. Contour Edge Analysis for Polyhedron Projections // Geometric Modeling: Theory and Practice. New York: Springer-Verlag New York, 1997, pp. 379-394.

15. Ф. Препарата, М. Шеймос. Вычислительная геометрия. Введение. М.: Мир, 1989, 478 с.

16. И.Ю. Балабан. Алгоритм поиска пересечений множества отрезков. М.: Институт прикладной математики им. М.В.Келдыша РАН, 1994, препр.№45,30 с.

17. В.А. Сарычев, В.В. Сазонов, М.Ю. Беляев, Н.И. Ефимов. Повышение точности определения вращательного движения орбитальных станций Салют-6 и Салют-7 по данным измерений //Космические исследования, 1991, т. 29, № 3, с. 375-389.

18. Е.В. Бабкин, М.Ю. Беляев, Н.И. Ефимов, В.В. Сазонов, В.М. Стажков. Неуправляемое вращательное движение орбитальной станции Мир // Космические исследования, 2001, т. 39, №1, с. 27-42.

19. ГОСТ Р 25645.166-2004. Атмосфера Земли верхняя. Модель плотности для баллистического обеспечения полетов искусственных спутников Земли.

20. В.В. Сазонов, С.Ю. Чебуков, В.И. Абрашкин, А.Е. Казакова, А.С. Зайцев. Низкочастотные микроускорения на борту ИСЗ Фотон-11 // Космические исследования, 2004, т. 42, №2, с. 185-200.

21. В.Н. Лихачев, В.В. Сазонов, А.И. Улъяшин. Одноосная солнечная ориентация искусственного спутника Земли // Космические исследования, 2003, т. 41, №2, с. 174-185.

22. В.Н. Лихачев, В.В. Сазонов, А.И. Улъяшин. Эволюция орбиты искусственного спутника Земли с солнечным парусом // Космические исследования, 2004, т. 42, №1, с. 83-87.

23. В.Н. Лихачев, В.В. Сазонов, А.И. Улъяшин. Режим одноосной солнечной ориентации искусственного спутника Земли. М.: Институт прикладной математики им. М.В.Келдыша РАН, 2001, препр.№15,32 с.

24. В.Н. Лихачев, В.В. Сазонов, А.И. Улъяшин. Исследование режима одноосной солнечной ориентации искусственного спутника. М.: Институт прикладной математики им. М.В.Келдыша РАН, 2002, препр.№65,32 с.

25. ММ Комаров, Д.Н. Климович, В.В. Сазонов. Расчет сил и моментов светового давления, действующих на роторный солнечный парус. М.: Институт прикладной математики им. М.В.Келдыша РАН, 1995, препр.№59,24 с.