CC BY
113
МОДЕЛИРОВАНИЕ ПРОЦЕССОВ |
УДК 629.78
В. В. Корянов
МЕТОДИКА РАСЧЕТА ПАРАМЕТРОВ ДИНАМИКИ ДВИЖЕНИЯ СПУСКАЕМОГО АППАРАТА ПРИ ЖЕСТКОЙ ПОСАДКЕ НА ПОВЕРХНОСТЬ ПЛАНЕТЫ
Сформулирована математическая модель посадки спускаемого аппарата на поверхность планеты с учетом относительного движения полезного груза внутри аппарата. Разработан программный комплекс решения уравнений, описывающих этот процесс. Показана необходимость применения демпфирующих материалов для полезного груза. Приведены результаты расчетов перегрузок корпуса спускаемого аппарата и полезного груза для различных параметров подхода к поверхности и видов грунта.
Основная цель работы — сформировать методику расчета динамики движения спускаемого аппарата (СА) при жесткой посадке на поверхность планеты, учитывая относительное движение полезного груза внутри СА.
Работа посвящена заключительному этапу экспедиции "Фобос-Грунт" — посадке СА с образцами грунта на поверхность Земли. Этот процесс сопровождается большими перегрузками, которые могут превышать значения, допустимые для контейнера с полезной нагрузкой. Возможным путем снижения перегрузок, действующих на контейнер, служит применение различного рода демпфирующих устройств.
Математическая модель. Системы координат. Все рассматриваемые четыре системы координат являются прямоугольными и правыми (рис. 1, 2).
Нормальная земная система координат ОиХиУиZи является инер-циальной (см. рис. 1). Начало координат Ои находится на поверхности Земли. Оси ориентированы следующим образом: ОиХи — местная горизонтальная ось, направленная по касательной к поверхности; ОиУи — местная вертикальная ось, направленная вверх по нормали к поверхности; О^и — местная фронтальная ось, перпендикулярная осям ОиХи и ОиУи, дополняет систему координат до правой.
Наклонная система координат ОнХнУнZн (см. рис. 1) используется для задания угла наклона поверхности в районе посадки и характеризует макрорельеф поверхности. Начало координат Он совпадает с началом инерциальной системы координат Ои, плоскость ОнXнZн образует посадочную поверхность. Ось ОнХн наклонена относительно оси ОиХи; ось ОнУн — перпендикулярна к поверхности посадки и
Рис. 1. Расположение наклон- динат, связанных с корпусом СА
ной системы координат и контейнером
направлена вверх от нее, Он— ось, дополняющая систему координат до правой.
Система координат OXYZ (см. рис. 1, 2) жестко связана с твердым недеформируемым корпусом СА. Начало координат О совпадает с центром масс корпуса СА. Оси связанной системы координат ориентированы следующим образом: ОY направлена по продольной оси к хвостовой части; ОX расположена в поперечном направлении относительно корпуса СА; OZ дополняет систему координат до правой. В этой системе координат рассматривается относительное движение контейнера.
Система координат ОРХРYpZp (см. рис. 2) жестко связана с твердым недеформируемом корпусом подвижного контейнера. Начало координат ОР совпадает с центром масс контейнера Р. Оси связанной с контейнером системы координат ориентированны следующим образом: ОРYP направлена по продольной оси контейнера в направлении хвостовой части СА; ОРХР расположена перпендикулярно продольной оси контейнера; ОРZP дополняет систему координат до правой.
Уравнения движения СА. Аппарат конструктивно можно представить как систему двух абсолютно жестких тел: корпуса с центром масс в точке О и подвижного контейнера с центром масс в точке Р (рис. 3).
Поскольку положение центра масс контейнера изменяется во времени, то положение центра масс системы (точка С) также будет изменяться. Принимается шесть степеней свободы контейнера в относительном движении.
Пространственное движение твердого тела в преграде описывается следующими уравнениями.
Уравнение движения центра масс системы имеет вид
ШУОС = Р, (1)
где т — масса системы; г£с — ускорение центра масс системы; Р — главный вектор сил сопротивления со стороны преграды и силы тяжести.
Ускорение центра масс системы определяется по формуле
г£е = г£о + £ х £е + ш х (ш х £е) + г£ег + 2ш х ге,
где (в системе координат ОиХиУ^и) г£о — ускорение точки О; £ — угловое ускорение корпуса; ш — угловая скорость корпуса; в системе OXУZ — £е — радиус-вектор точки С; г£ег — относительное ускорение точки С; £сг — относительная скорость той же точки.
Уравнение вращательного движения аппарата можно представить
как
—К
-ко = т(£е х £о) + М, (2)
где Ко — кинетический момент системы относительно центра О; М — главный момент всех внешних сил относительно центра О; гйе — вектор скорости центра масс системы; £о — вектор скорости полюса.
Уравнение относительного движения контейнера (в связанной системе координат) имеет вид
тР трг = Ё — тР г£р е — тР гиоре, (3)
где тР — масса подвижного контейнера; г£Рг — относительное ускорение центра масс контейнера (точки Р); Ё — равнодействующая всех сил, действующих на контейнер; г£Ре и г£Ре — переносное и кориоли-сово ускорения точки Р.
Уравнение вращательного движения контейнера можно записать
как
АКр -
-ж = МР,
где К£Р — кинетический момент контейнера относительно центра Р; МР — главный момент сил сопротивления со стороны амортизации на подвижный контейнер.
Для определения параметров движения СА при взаимодействии с преградой решается система уравнений [1], представляющая проекции уравнений (1)-(3) на оси связанной с корпусом СА системы координат OXУZ (см. рис. 2), проекции уравнения вращательного движения контейнера на оси связанной с ним системы координат ОРХРУРZP, кинематические уравнения Эйлера и уравнения, определяющие переход от связанной с корпусом системы координат к неподвижной °иХиу^и.
Рассмотрим вывод уравнений движения для СА, представленного системой из двух жестких тел (см. рис. 3). Система имеет следующие
параметры: шк — масса корпуса; шР — масса контейнера; /Х,/у ,/^ — моменты инерции корпуса относительно осей OXYZ; /Р,/Р, /Р — моменты инерции подвижного контейнера относительно осей системы, связанной с центром масс контейнера Р; гР (хР ,уР ,гР) — радиус-вектор точки Р; гс (хс , ус, гс) — радиус-вектор точки С. Для гР, гс — полюс в точке О.
Координаты точки С можно выразить Рис. 3. Конструктивная схема через координаты точки Р следующим аппарата образом:
ШР ШР ШР
xC = ,
mk + mP Учитывая, что
xp; Ус =
mk + mp
yp; zc =
mk + mp
zp.
mk + mp = m,
проекции скорости и ускорения системы в связанной системе координат OXYZ можно представить в виде:
mp
x C =-x p; x p = Vprx;
m mp
У с =-У p; Ус p = vpry;
m mp
Zc =-Zp; Zp = vprz;
m
(4)
mp
xc =-xp;
m mp
Ус = — yp;
m mp
Zc =-Zp;
m
xp = Wprx = vc prx; yp = WprY = V pry; Zp = WprZ = v prZ,
(5)
где г>Рг (^РгХ, г>РгУ, ) и ^Рг (^РгХ, ^РгУ, ) — векторы относительных скорости и ускорения точки Р относительно корпуса.
Здесь и далее все указанные векторы задаются проекциями на оси связанной системы координат OXYZ. Движение контейнера вокруг центра масс Р задается угловой скоростью (^РХ, ^Ру). Этот вектор задан проекциями на оси системы координат ОРXPYPZP.
Для производной вектора скорости по времени имеем [2] =
dt
dt
+ ^ х vO.
В проекциях на оси связанной системы координат это уравнение можно представить в виде
Wox = vOX + UyVoz - UzVoy ;
WOY = vOY + uzVox - ихVoz ; (6)
Woz = vOZ + Ux Voy - uy Vox .
Подставив выражения (6) в уравнение для Wjc и разложив на составляющие, ускорение центра масс системы можно представить как
Wcx = Vox + иуVoz - uzVoy - (up + u|) xc-
- (Uz - ШхШу)yc + (Uy + ШхUz)zc + 2(uyzc - Uzyc) + Xc;
WCY = Voy + uz Vox - их Voz + (UJ z + их иу )xc -
- (uZ + и2х) yc - (cJx - иу uz )zc + 2(uz X c - их Zc) + yc;
wcz = Voz + их Voy - иу Vox - (uJy - их Uz )xc+
+ (cJ X + Шу UZ )yc - (u1 + Ш1) Zc + 2(ux yc - иу x c) + Zc.
Подставив в выражения (7) соотношения (4) и (5), уравнение движения (1) в проекциях на оси связанной с корпусом системы координат можно представить в окончательном виде.
Аналогично выведем уравнения вращательного движения СА.
Кинетический момент системы можно представить в виде
KO = К0орп + Коонт,
где К0орп и К0онт — кинетические моменты корпуса и подвижного контейнера относительно центра О.
Суммарный кинетический момент системы в проекциях на оси связанной системы координат можно записать как
Кох = (Ix + IX + mP (Ур + z1 ))ux;
TD /О О N N
Koy = (Iy + Ip + mp (xP + Zp ))иу;
Koz = (Iz + IZ + mp (xp + yp ))uz .
Уравнение вращательного движения СА (2) и уравнение вращательного движения контейнера в проекциях на оси системы координат, связанной с корпусом и контейнером соответственно можно представить в окончательном виде.
Дополняя полученные ранее уравнения кинематическими, получаем полную систему дифференциальных уравнений.
Уравнения движения центра масс системы:
m(V ох + иу Voz - uz Voy ) - mp [(up + up)xp + (cj z - их иу )yp --(cjY + UxUz)zp - 2(иуVprZ - UzVprY) - VprX] = Gx + Px;
m(v OY + Uz Vox - их Voz ) + mp [(cJ z + их иу )xp - (u2Z + up)yp -
(uJ X + Uy UZ )zp + 2(uz VprX - UX VprZ) - VprY ] = Gy + Py;
m(vVozVoy-шу vox )-mp [(ш y)xp -(ш x +шу )yp+
+ + шУ )zp - 2(шхWpry - шу VprX) - vVPrZ] = Gz + Pz ,
где Gx, Gy, GZ — проекции силы тяжести на оси OXYZ, Px, Py, PZ — проекции главного вектора сил сопротивления со стороны преграды P. Уравнения вращательного движения корпуса:
[Ix + If + mp + zp)]шx = Mj + mp [(шгxp - шхzp + vpry)voz--(шхyp - Шуxp + VprZ)voy - 2шх(vpryyp + VprZZp)]--[Iz + Ip + Iy - Ip + mp (yp - zp)]шу шг;
[Iy + Ip + mp (xp + zp )]ш y = My + mp [(шх yp - шу xp + vprz )vox --(шуZp - шгyp + vprx)voz - 2шу(vprxxp + VprZZp)]--[Ix + Ip - Iz - Ip + mp (zp - xp)^xшг;
[Iz + Ip + mp (yp + xp )]ш z = Mz + mp [(шу zp - шг yp + vprx )voy --(шгxp - шxzp + vpry)vox - 2шz(vprxxp + vpryyp)]--[Iy + Ip - Ix - Ip + mp (xp - y! )^x шу.
Уравнения движения груза:
mp vv prx = Fx - mp [vv ox + шу Voz - шг Voy --(ш2 + шг)xp - (шz - шушx)yp + (шу + шxшг)zp+ +2(шуvprz - шгvpry)];
xp = vprx;
mp vVprY = Fy - mp [vVoy + шг vox - шx voz+ +(ш z + шx шу )xp - (шг + ш2 )yp - (iüx - шу шг )zp+ +2(шгvprx - шxvprz)]; yp = vpry;
mp vv prz = Fz - mp [vv oz + шx voy - шу vox -(ш y - Шx шг )xp + (ÚJx + шу шг )yp - (ш2 + ш2 )zp+ +2^xvpry - шуvprx)];
zp = VprZ.
Переход от связанной системы координат к неподвижной:
xO = vox cos $ cos ф + vOY(- cos y sin $ cos ф + sin y sin ф) + +vOZ (sin y sin $ cos ф + cos y sin ф);
yO = vox sin $ + vOY cos y cos $ + vOZ(- sin y cos $);
zO = vox (- cos $ sin ф) + vOY (cos y sin $ sin ф + sin y cos ф) + +vOZ(- sin y sin $ sin ф + cos y cos ф),
где $, ф и y — углы тангажа, рыскания и крена соответственно. Кинематические уравнения Эйлера:
ф = (шу cos y - шг sin y)/ cos $;
$ = uY sin y + uZ cos y; Y = uX - tg $(uY cos y - uZ sin y).
Уравнения вращательного движения груза:
IXúxp + (Ip - Ip)uypuzp = mxp; Ipшуp + (IX - Ip)uxpuzp = myp ; Ipüzp + (Ip - IX)uxpuyp = mzp .
Кинематические уравнения для груза:
фp = (uYp cos yp - иzp sin yp ) / cos $p ; фp = uYp sin yp + uzp cos yp ;
Yp = uxp - tg $p (uyp cos yp - uzp sin yp ).
В приведенной системе уравнений xo,y0,zo — координаты центра О в неподвижной системе координат; xp,yp,zp — координаты центра P в связанной с корпусом системе координат.
Интегрирование системы дифференциальных уравнений осуществляется методом Рунге-Кутты четвертого-пятого порядка [3].
Математическая модель влияющих силовых факторов. Для определения значений P и M главного вектора и момента сил сопротивления в каждый момент времени необходимо знать распределение "нормального" (an) и "касательного" (аТ) удельных сопротивлений на поверхности проникающего тела:
Un = Av2p + BVn + C, от = (8)
где A, B, С — коэффициенты, характеризующие свойства преграды (грунта) оказывать сопротивление; ц, — коэффициент трения материала преграды (грунта) о корпус тела.
При вычислении силовых факторов на боковой поверхности тела, которая при проникании взаимодействует с инерционно расширяющейся средой, приводимой в движение лобовой частью, скорость движения среды в окрестности рассматриваемых точек боковой поверхности учитывается через vn в указанной выше зависимости.
В этих точках vn находится как vn = РПотн = Vo-гнП = (V - Ve)n, где vn — проекция скорости рассматриваемой точки на нормаль к поверхности в этой точке; vnOTH — проекция скорости рассматриваемой точки боковой поверхности относительно инерционно расширяющейся среды на нормаль к боковой поверхности тела в этой же точке; V — скорость рассматриваемой точки тела; Ve — скорость инерционно расширяющейся среды в окрестности рассматриваемой точки; n — вектор нормали в рассматриваемой точке боковой поверхности тела (рис. 4).
Рис. 4. Графическая интерпретация разложения вектора скорости на составляющие
Скорость рассматриваемой точки тела находим по формуле Эйлера
ус х - г у ^ + г^ Оу
У = (ус X ,ус у ,ус я) - Т X (ох )= ( ус у - г^ ^Х + Гх ^
У СЯ - гхОу + Гуох
Проекции скорости на нормаль к поверхности в каждой точке определяют как
Уп = УпХ + Упу + Уп^ = Ух Пх + Уу Пу + у^ .
Проекции скорости на касательную плоскость к поверхности в каждой точке находят по формулам:
УтХ = Ух - УпХ, Уту = Уу - Упу, Ут^ = У^ - Уп^.
Итак, получаем касательную составляющую скорости
отсюда имеем
Vt = \/ v2X + VrY + v2Z >
VTX VTy Vtz
= -, TY = -, = -.
vT vT vT
Подставляя значение уп в уравнения (8), получаем значения для
Сила сопротивления и ее момент определяются из следующих
уравнений, спроецированных на оси связанной системы координат:
N „
Р = - ^ / (^пПг + Т, г=1 У
и aT.
N
M = - ^ / Pi X (OnUi + GrTi)dSi,
где i — индекс, показывающий на какой характерной поверхности тела (сферической, цилиндрической и т.п.) определяются силовые факторы; Ski — часть i-й поверхности, находящейся в контакте с преградой; иг — единичный вектор нормали; тг — единичный вектор касательной к поверхности; dSi — элементарная площадка на i-й поверхности; рг — радиус-вектор рассматриваемой точки тела в связанной системе координат.
Выражения для щ, dSi, рг целиком определяются геометрией тела, и для каждой поверхности записываются по-своему.
Для каждого элемента можно записать:
dPX = — (onuX + oTтх) dS = —ахdS; dMX = — (pYoz — Pzoy) dS; dPY = — (onuy + oTтy) dS = —oydS; dMy = — (pzoX — pxoZ) dS;
dPz = — (onuz + Ottz) dS = —ozdS; dMz = — (pxoy — pyах) dS.
Окончательно имеем для P и M:
Поверхность Бк, находящаяся в контакте с преградой, определяется в процессе интегрирования по всей поверхности проникающего тела. Для этого вся его поверхность разбивается на элементарные площадки с узловыми точками, в каждой из которых проверяются условия контакта поверхности тела с преградой.
Алгоритм моделирования. Влияние формы тела и его конструкции на силовые факторы сопротивления среды сводится к определению радиуса г сечения тела в рассматриваемой точке, элементарной площадки dSi в окрестностях этой точки и проекции вектора внешней нормали на поверхности тела в этой же точке пхи пУг, пг%.
При взаимодействии корпуса СА с грунтом его начальное положение определяется в неподвижной системе координат ОиХиУи Zи, а скорости — в связанной системе координат ОХYZ. Начальное положение корпуса задается координатами полюса О в неподвижной системе координат.
Составляющие начальной скорости центра О представлены также в неподвижной системе координат, при помощи матрицы перехода они переводятся в связанную с корпусом систему координат.
Px = dPx; Py = ^ dPY; Pz = ^ dPz; Mx = ^ dMx; My = ^ dMY; Mz = ^ dMz
z.
Положение подвижного контейнера при взаимодействии с амортизатором определяется в системе координат, связанной с корпусом СА, а скорости — в системе координат 0РХРУРZP, связанной с контейнером. Начальное положение контейнера задано координатами центра масс Р в системе координат 0XYZ.
Решение сводится к параллельному выполнению двух численных процессов: на каждом шаге интегрирования по времени определяются силовые факторы для корпуса СА при взаимодействии его с грунтом и для подвижного контейнера при взаимодействии с амортизатором и проводится непосредственно интегрирование дифференциальных уравнений.
При определении силовых факторов наиболее целесообразно применение обычных кубатурных формул. Расчет по ним предусматривает разделение поверхности на малые площадки, определение силовых факторов йР и йМ, действующих на каждую из них, по значениям подынтегральных функций в "узловых" точках площадки и последующее их суммирование для получения главных векторов Р и М. При расчете в качестве области интегрирования берется вся поверхность тела, а подынтегральная функция определяется из условий контакта точек тела со средой.
В расчетах сферу разбивают на элементарные площадки следующим образом (рис. 5):
1. Первый параметр разбиения а — деление сферы вдоль диаметра большого круга с шагом йа (рис. 5, а).
2. Каждое полученное сечение-окружность разбивается на равные дуги с шагом йЬ (рис. 5, б).
Тогда получаем: тх = а +--, Ях = у/— тХ, йв = Ях йа йЬ.
2
Разработанная программа универсальна (рис. 6). Она позволяет проводить расчеты по динамике посадки СА [4] различного конструктивного исполнения: корпус СА и подвижный контейнер могут быть представлены набором различных поверхностей (сфера, конус, тор, цилиндр и др.).
Рис. 5. Разбиение сферы на равное число шагов вдоль диаметра большого круга (а) и разбиение окружности сечения на равные дуговые участки (б)
Рис. 6. Принципиальная схема решения задачи
Рассматриваемая форма СА (см. рис. 3) относится к классу сегментально-цилиндрических тел. Поверхность корпуса аппарата можно представить сочетанием трех поверхностей: сферической, цилиндрической и сферической в хвостовой части. Поверхность СА подвижного контейнера представлена притуплением и тороидальной поверхностью. Конструкция лобовой части представляет собой сферу, переходящую в боковую цилиндрическую поверхность. Цилиндрическая часть корпуса в хвостовой части гладко сопрягается со сферой.
Проведение тестовых расчетов. Параметры подхода СА к поверхности следующие: вертикальная составляющая скорости 40 м/с, горизонтальная составляющая скорости 5 м/с, масса спускаемого аппарата 8,7 кг, характерный размер системы 0,5 м.
Виды грунтов: сыпучий и жесткий с наклоном площадки 15°.
Для нахождения коэффициентов А, В, С, А1, В1, Сь характеризующих свойства преграды оказывать сопротивление, необходимо наличие данных об удельном сопротивлении ап для разных скоростей.
Результаты расчетов представлены графическими зависимостями (рис. 7). Проведенные расчеты показали, что в случае использования
Рис. 7. Зависимость перегрузки от времени для корпуса (а) и груза (б). Два вида грунта — жесткий (1, 5) и сыпучий (2, 4) и два вида амортизатора — жесткий (1, 2) и мягкий (3, 4)
мягкого амортизатора перегрузки для груза имеют меньшие значения (см. рис. 7, б), однако перемещение груза внутри корпуса имеет критические значения и груз достигает корпуса СА (см. рис. 8, б). В случае использования жесткого амортизатора перегрузки для груза больше (см. рис. 7, б), перемещение меньше и не критическое (рис. 8, в).
Выводы. 1. Сформированы следующие физико-математические модели процесса:
а) посадка СА с учетом возникающих силовых воздействий от поверхности планеты;
б) относительное перемещение полезного груза внутри СА с учетом возникающих силовых воздействий на корпус аппарата при работе системы амортизации.
2. Разработан программный комплекс решения уравнений, описывающих этот процесс, который позволяет использовать полученные данные в дальнейшем при анализе результатов расчета. Возможно-
Рис. 8. Ключевые моменты процесса внедрения для различных типов амортизатора и грунта при одинаковых условиях подхода к поверхности:
а — момент подхода к поверхности, t = 0 c; б — мягкий амортизатор, сыпучий грунт, t = 0,01 c; в — жесткий амортизатор, сыпучий грунт, t = 0,01 c
сти программного комплекса проиллюстрированы примером расчета (рис. 8). Показана принципиальная возможность снижения значений составляющих перегрузки для полезного груза при работе системы амортизации.
3. Анализ графиков перегрузки при посадке СА на различные грунты показал, что в случае посадки на жесткую поверхность перегрузка для контейнера имеет двойной максимум. В то же время при посадке на сыпучую поверхность кривая перегрузки для контейнера имеет только один максимум.
Работа выполнена совместно со специалистами НПО им. С.А. Лавочкина.
СПИСОК ЛИТЕРАТУРЫ
1. Дмитриевский А. А., Лысенко Л. Н., Богодистов С. С. Внешняя баллистика. - М.: Машиностроение, 1991. - 638 с.
2. Курс теоретической механики: Учеб. для вузов / В.И. Дронг, В.В. Дубинин, М.М.Ильин и др.; Под общ. ред. К.С.Колесникова. - M.: Изд-во МГТУ им. Н.Э. Баумана, 2005. - 736 с.
3.Разоренов Г. Н. Теоретические основы управления полетом баллистических ракет и головных частей. - М.: Военная академия РВСН имени Петра Великого, 2001. - 405 с.
4. Б а ж е н о в В. И., О с и н М. И. Посадка космических аппаратов на планеты. - М.: Машиностроение, 1978. - 158 с.
Статья поступила в редакцию 25.03.2008
Всеволод Владимирович Корянов родился в 1982 г., окончил МГТУ им. Н.Э. Баумана в 2006 г. Аспирант кафедры "Баллистика и аэродинамика" МГТУ им. Н.Э. Баумана. Специализируется в области моделирования динамики движения спускаемых аппаратов.
V.V. Koryanov (b. 1982) graduated from the Bauman Moscow State Technical University in 2006. Post-graduate of "Ballistics and Aerodynamics" department of the Bauman Moscow State Technical University. Specializes in the field of simulation of dynamics of descent vehicle movement.
