CC BY
60
Том XXII
УЧЕНЫЕ ЗАПИСКИ ЦАГИ
1991
№ 1
УДК 629.7.015.076.8 + 629.78.015
АЛГОРИТМЫ ОБХОДА ОСОБЕННОСТИ ПРИ ЧИСЛЕННОМ РЕШЕНИИ УРАВНЕНИЙ ДВИЖЕНИЯ ТЕЛ ОТНОСИТЕЛЬНО ЦЕНТРА МАСС В АТМОСФЕРЕ ПРИ ДЕЙСТВИИ ВОЗМУЩЕНИЙ
Предложен эффективный алгоритм обхода особенности в кинематических уравнениях, применимый при их численном решении в комбинированных быстродействующих методах исследования пространственного движения в атмосфере твердых тел вокруг центра масс при наличии возмущений [1 —3]. На основе системы уравнений, использующей в каче1:тве кинематических переменных параметры Родрига —Гамильтона и не .содержащей поэтому особенности, разработан алгоритм перевода параметров относительного движения из кватернионного представления в углы атаки и крена и обратно.
1. Для успешного проведения обширных параметрических исследований задач динамики твердого тела в атмосфере в [
ческая модель, имеющая минимальную трудоемкость численного исследования на ЭВМ (с точки зрения вычислительных затрат и удобства машинной реализации) при заданной точности. Отличительной чертой этой модели является использование для описания ориентации тела относительно инер-циальной системы координат параметров Родрига — Гамильтона/ скомпонованных в кватернион [4].
на точность и быстродействие проводилось ее сравнение с другими математическими моделями движения твердого тела в атмосфере, использующими направляющие косинусы и углы Эйлера. Результаты сравнения приведены в [
кратно проходило вблизи особой точки системы уравнений движения в углах Эйлера. Поэтому для решения задачи в углах Эйлера был разработан специальный алгоритм обхода особой точки, описываемый в настоящей статье.
Система уравнений движения в углах Эйлера, использованная при исследованиях, имеет следующий вид [4] :
Г. М. Лохов, С. И. Подзоров
Т« + » X X« = мст + Ма + Мвозм;
Р = ю* + (ш^соэр — ш^пр)^а; у = (ш^т р — wi/cosP)/siпa,
где векторы г, V, Рв03м — соответственно радиус-вектор центра масс, вектор'его скорости и возмущающей силы — заданы в инерциальной системе координат,
Мст = -—-Б¿тг(а)(0, втр, соэр)1, Мд = (т?, т"*, т?)т -
заданы в связанной системе координат, J — тензор инерции тела в главных осях.
Выбор для исследований системы уравнений (1) обусловлен тем, что при ее . интегрировании углы а, р, у, необходимые для вычисления аэродинамических характеристик тела, определяются непосредственно, без каких-либо дополнительных вычислений.
2. Система уравнений (1) содержит в себе методическую погрешность, вызванную неучетом вращения полускоростной системы координат [5] за счет действия силы гравитации и подъемной силы. Рассмотрим этот вопрос подробнее.
Кинематические уравнения относительного движения в системе (1) определяют углы а, р, у, характеризующие ориентацию связанной системы координат по отношению к полускоростной. Поэтому шх,у,г—это компоненты ш*— вектора относительной угловой скорости вращения связанной системы координат относительно полускоростной. В динамических уравнениях движения вокруг центра масс в системе (1) «—это вектор угловой скорости вращения связанной системы координат относительно инерциальной. Полускоростная система координат вращается, поэтому Разность = ш—обуслов-
лена поворотом вектора скорости за счет действия подъемной с!"лы У (силой гравитации для простоты пренебрегаем). Тогда а — «к = [VXV] /У2, откуда ||ш|| = О (У/тУ), т. е.
(ю*)*,
подставляя это выражение в кинематические уравнения системы (1), получим
<х= ш^пр + шго»р + О (У/тУ), Р = ю* + (юуСОзр-ш^тр) О (У/тУ),
V = (шгsiпр —ШyCOS^ ^та + О (У/тУ),
где Их
координат относительно инерциальноЙ.
Величина -^уг = су(а) для типичных значений параметров тела и внешней среды много меньше значения а.
3. Особыми точками кинематических уравнений системы (1) являются а= О и а = л, Для определенности рассмотрим поведение решения в окрестности а=О, т, е. при аЕ [—а*, а*]. Величина а* выбирается малой, так что численное интегрирование (1) в данной окрестности особой точки затруднительно из-за близости к ней, Сложив последние два уравнения системы (1), получим при соза«1 :Р + у=ш*, а при Шx = coпst соотношение:
р + 7= ро + 7о + ш,*, 2)
где индекс О соответствует моменту входа решения в окрестность а=О, принятому за начальный.
Из (2) находится один из углов р или у по известному другому. Определим приближенную зависимость у (/), использовав аппроксимацию относительного движения тела в окрестности а=О регулярной прецессией при постоянных параметрах движения центра масс. Такое приближение допустимо при малых величинах а*, так как ||МСТ|| пропорционален |а| — малому значению вблизи особой точки а= 0. Дополнительно предположим малость Мв и Мво;>м- Тогда движение тела вокруг центра масс будет представлять собой вращение вектора продольной оси п вокруг вектора кинетического момента К, постоянного в инерциальной системе координат, с постоянной угловой скоростью 'Ф.
Введем инерциальную систему координат Ое\егез с началом в центре планеты и систему координат OxcУcZc с началом в центре масс тела: ось Охс направлена по вектору скорости V, ось Оус — по вектору е2 Х V, ось Ozc дополняет систему до правой [5]. Переход от введенной системы координат к связанной системе осуществляется последовательным поворотом на углы а, р, у, заданные при 1 = О, с помощью матрицы перехода, соответствующей ориентации тела в начальный момент времени. Тогда можно найти в системе координат OxcУcZc векторы Пц_о и К, заданные в связанной системе координат компонентами
п = (1, О, О), К= (
Зависимость п (() формулам конечного поворота:
к? К[К2 К(Кг 0 -к. к2
n (/) =' Ecosi|) + К2К{ кi к2к3 (1 — COS^ + К3 0 -К. sin-ф
К3К{ к3к2 -к2 к 0
где Е—единичная матрица, 'Ф = в системе координат OxcfcZc.
Для вычисления 'Ф разложим вектор (J) на две составляющие, направленные по векторам Кип:
(1) = 'Ф К + 3)
Умножая 3) скалярно на К и n и разрешая полученную систему уравнений относительно 'Ф, найдем ее величину: 'Ф = 11 KII Jz. По зависимости n (()
Y=arccos (zc, [ п]) sign (п, zc).
С учетом того, что в системе координат OxcfcZc V = (V, О, О), получим у = = arccos rt2Csign пзс [2].
Теперь можно определить значения углов ß (1) и а (1) = arccos п^. Знак а (1) остается постоянным на всем расчетном участке [2].
a(t) находим момент времени tF — выхода решения из окрестности особой точки а=О. Далее, численно проинтегрировав с помощью метода Эйлера уравнения движения центра масс твердого тела (первые два уравнения системы (l)) с известными функциями а (t) и ß (t) на отрезке [О, Ifj, находим значения r((
высокого порядка нецелесообразно, так как зависимости а (/) и ß(t) выведены в предположении постоянства векторов r и V в окрестности особой точки.
4. Для завершения построения алгоритма обхода особенности а=О сформулируем оценку величины а*, конкретизировав тем самым требование
малости аэродинамических моментов (примем здесь МВОзМ = О), и оценку шага h интегрирования уравнений движения центра масс методом Эйлера. Обозначим через е, eK допустимые величины ошибок вычисления а и V по разработанному алгоритму. Ошибку вычисления угла атаки оценим как угол поворота вектора К за время tf прохождения решением окрестности особой точки [—а*, а*]. Тогда получим условие AK/K^e, в котором AK«(MCT + + М,,) tf. Величину tf оценим, приблизив a(
==0(1) и oo2«MCT(a)/(a1г) [5]. Отсюда tf~a*/oo, а при MCT</Zoo2a* и К= = О (/zoo) получаем
а*=О (-Ге). 4)
Чтобы оценка (4)
« Ма или Ут;а*« /оот^. В рассматриваемых в данной работе задачах это условие выполняется.
Для оценки шага Л выпишем дифференциальное уравнение для У с учетом приближения аэродинамических коэффициентов зависимостями сх=сХО + ха2, с*, тт" = соп51, справедливыми при малых значениях угла атаки:
(а + а2) + 3ь (5)
где 3 = 5(х + ф, Р1 = сопз1.
Подставляя в (5) приближение для а(() величины погрешности метода Эйлера, получим оценку для шага интегрирования:
При задании допустимых величин ошибок вычисления а и У следует учитывать, что погрешность вносится в вычисления при каждом прохождении решения вблизи особой точки, т. е. накапливается. Поэтому, задав глобальные погрешности вычисления а и У соответственно Е и £у, величины е и еу можно оценить как е« 2лЕ/ (Тоо), еу« 2лЕу/ (Тоо), где Т—продолжительность расчетного участка.
5. Разработанный алгоритм обхода особенности в кинематических уравнениях движения остается применимым и в случае учета дополнительных слагаемых, обусловленных действием подъемной силы (заметим, что в окрестности особой то':ки а= О ||юу|| = О (су (а)) есть величина первого порядка малости по а). Однако учет этих членов необходимо осуществить таким обр-азом, чтобы вид кинематических уравнений в системе (1) не изменился, т. е. ввести дополнительные члены в динамические уравнения.
Подставив 00= оол + Юу в динамические уравнения движения относительно центра масс, получим
J (Ш* + Юу) + (м;. + «у) Х J (»* + ®у) = М};, или, отбрасывая члены порядка а2:
+ [»*Х .Юу + МуХ .!«*] = М}; -.«у,
Величина в квадратных скобках является известной функцией а, р, и параметров движения центра масс твердого тела. Значение выражается через эти же переменные путем дифференцирования по времени выражения
^Х V] /У2
Оценим эффект учета члена .«у. Величина ¡Шу| = О =
= ¡6)*Ц)> т. е. добавочный член пропорционален (1). и является
демпфирующим. Эффект этого демпфирования соответствует добавлению к
(
( II iIIc")
коэффициентам демпфирования величины О (—, которая будет примерно
одного порядка с самим коэффициентом для типичных значений параметров твердого тела. В частности, для случая плоского движения твердого тела
К),фф = m + (Jz/mfK [2].
6. Для исследования движения тела в плотных слоях атмосферы, когда относительное движение является пространственным, а движение центра масс можно считать плоским, при действии кратковременного возмущающего момента М* = М/Jz, перпендикулярного к плоскости траектории, через период колебаний по углу атаки, обычно применяется следующая система уравнений
; = Vsin0, сер = — cose,
ее=.(---)«»*
V --JLsine - ^scx(a), R = — -PJ-SP m; R + М* siny sina, ó = - SI2 TOJ - - R cosa) + —
- S/2m?) R — -^si n 2asin у,
G — Rcosa v pV'2 с / >
У =---. У = -^s-Sc/cO,
sin a ¿
•• . (G — Reasa) (R — Gcosg) . (™
-si^a m
-|I/mV)á + -|£-S/mz(a) = М* cosy
где J O, J z R — соответственно проекции вектора кинетического момента на направление вектора скорости и на продольную ось тела.
При непосредственном численном решении системы уравнений (7) может возникнуть трудность, связанная с наличием сингулярности в уравнении для а, так как при малых углах атаки второе слагаемое в нем есть О (а-1). В момент начала действия возмущений а= О, но G = R, и особенности не возникает, так как
lim <* —Rosal— Iim R -Rcosa =«
а_о sm3a а—О sin2 а 2
В процессе действия возмущения угол атаки а становится ненулевым и появляется разность G — R,*0. В отсутствие возмущений эта разность не обращается в ноль [6],
Отсюда следует постоянство знака а на этих участках [2], не проходит через ноль, и формально особенность в системе (7) на этих участках не возникает. Однако на практике а может принимать весьма малые (хотя и ненулевые) значения, что из-за наличия сингулярности в уравнении для а существенно затрудняет численное интегрирование (7) или даже делает его . невозможным. Что же касается участков действия возмущений, то на них нельзя даже гарантировать отсутствие прохождений а через ноль.
7. Корректный обход особенности в системе уравнений (7) сложной задачей. Поэтому можно попытаться использовать вблизи особых точек системы (
ные с помощью кватернионов [4].
кватернионов состоит в необходимости перевода параметров относительного движения из кватернионного представления в представление системы (7) обратно. Рассмотрим вначале переход от представления (7) т. е. от б, Я, а, а, ' 7 к А, (1). Можно всегда считать, что в начальный момент времени связанная система координат совпадает с полусвязанной [ практически они различаются лишь поворотом на некоторый угол вокруг продольной оси тела, а в силу осесимметричности рассматриваемых тел такой поворот несущественен. Тогда справедливы следующие уравнения динамики [
mV
r = Vх,
Г = vect (л - fl. л) -f. + F х
ita + [ «5x m = м' + ml„,
л = - л - шЕ,
= = Sin а {jD,-dt + J^
где верхний индекс х указывает, что компоненты вектора_ взяты в инер-циальной системе, ; — в связанной, А = А (а, у, u*), u* = Vх/ V, А—сопряженный кватернион;
Сг
спи2 с„и з
ME = -2" Si2 diag (m"",
m,
m,
0
<0
4-sim,
. о
"з
— h
m , =
Кватернион IV можно построить, если известна матрица перехода из связанной системы в инерциальную [8].
ния поворотов вокруг соответствующих осей на углы а и у с учетом того, что вектор скорости лежит в плоскости Ов\в2 инерциальной системы и наклонен к оси Ое\ на известный угол 8 [8].
В итоге получены все необходимые начальные условия для интегрирования уравнений (
8. Несколько сложнее обратный переход. Угол атаки а с точностью до знака определяется по формуле
lal = arccos {п, u).
(9)
Выбор знака угла атаки произволен, так как при Мвозм = О уравнения (7) инвариантны относительно замены а на — а. Примем, что а^О. С целью построения алгоритма определения знака угла крена у введем скоростную систему координат с началом в центре масс твердого тела: ось Ох направлена по вектору' скорости V, ось Ог перпендикулярна к плоскости траектории и образует правую тройку с векторами V и местной вертикали, ось Оу дополняет систему до правой. Твердое тело переводится из положения а= О
(
С, =
sin а
ю
cü
в некоторое положение, характеризуемое углами а ! и у, следующим образом: вначале осуществляется поворот вокруг оси Ог на угол а>О, а затем плоскость (п, и) —плоскость угла атаки — поворачивается вокруг вектора V на угол у; при этом, очевидно, конец вектора п описывает дугу окружности с центром на линии вектора V. Будем считать у>О, если при наблюдении с конца вектора V эта дуга описывается против часовой стрелки. Тогда sigпy = = sigп (п, г0), модуль угла крена |y| = aгccos (г0, [иХп]),
y=aгccos (г0, [иХп])
Дифференцируя формулу (9), найдем
а= - [(и, п) +( (1
где.. = А (-V) = 1' ^Хп = [«Хп].
в (10) не возникает особенности из-за присутствия siпa в знаменателе, так как .в процессе действия возмущения угол атаки становится ненулевым.
9. Эффективный обход особенностей позволил успешно использовать уравнения (1) и (7) для построения регулярных асимптотических методов исследования пространственного относительного движения твердых тел в тех случаях, когда угол атаки может быть близок к О. Матема:гическая модель пространственного движения тела в атмосфере 8) с кинематическими переменными, параметрами Родрига —Гамильтона, позволяет без потери точности существенно уменьшить время прямого счета и применяется в быстродействующих, комбинированных алгоритмах расчета на переходных участках. В заключение отметим, что в [9,
алгоритмы исключения особенности в кинематических уравнениях и оптимальной коррекции матрицы направляющих косинусов, основанные на перестройке координатных осей.
ЛИТЕРАТУРА
1. Ло х о в Г. М., П о д з о р о в С. И. Численное исследование движения твердых тел в атмосфере.—Ж. вычисл. матем. и матем. физ., 1987, т. 27, № 2.
2. Ло х о в Г. М., П о д з о р о в С. И. К исследованию относительного движения
1987, № 2.
3. Л о х о в Г. М., П о д з о р о в С. И. О кинематических переменных в уравнениях движения твердого тела в атмосфере. — Космические исследования, 1988, т. 24, вып. 1.
4. Ло х о в Г. М., П о д з о р о в С. И. Применение кватернионов в задачах динамики твердого тела в атмосфере. — Космические исследования, 1984, т. 22, вып. 4.
5. О х о ц и м с к и й Д. Е., Г о л у б е в Ю. Ф., С и х а р у л и д з е Ю. Г. Алгоритмы управления космическими аппаратами ири входе в атмосферу.— М.: Наука, 1975. .
6. К у з м а к Г. ■ Е. Динамика неуправляемого движения летательных аппаратов при входе в атмосферу. — М.: Наука, 1970. .
7. Ло х о в Г. М., П о д з о р о в С. И., С и д о р е н к о В. В. К исследованию пространственного движения твердого тела в атмосфере при действии возмущений. — М.: 1986. (Препринт № 64. — Ин-т прикл. матем. им. М. В. Келдыша АН СССР). '
8. Лурье А. И. Аналитическая механика. — М.: Физматгиз, 1961.
9. Ш и л о в А. А. Об исключении особенности в общих уравнениях движения летательного ■ аппарата. — Инженерный журнал, 1962, т. 2, вып. 3.
10. Ш и л о в А. А. Оптимальная коррекция матрицы направляющих косинусов при расчетах вращения твердого тела. — Ученые записки ЦАГИ, 1977, т. 8, 5.
Рукопись поступила 18/Х 1989 г.
