Научная статья на тему 'Алгоритмы обхода особенности при численном решении уравнений движения тел относительно центра масс в атмосфере при действии возмущений'

Алгоритмы обхода особенности при численном решении уравнений движения тел относительно центра масс в атмосфере при действии возмущений Текст научной статьи по специальности «Физика»

CC BY
219
59
i Надоели баннеры? Вы всегда можете отключить рекламу.
Журнал
Ученые записки ЦАГИ
ВАК
Область наук

Аннотация научной статьи по физике, автор научной работы — Лохов Г. М., Подзоров С. И.

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

i Надоели баннеры? Вы всегда можете отключить рекламу.

Похожие темы научных работ по физике , автор научной работы — Лохов Г. М., Подзоров С. И.

iНе можете найти то, что вам нужно? Попробуйте сервис подбора литературы.
i Надоели баннеры? Вы всегда можете отключить рекламу.

Текст научной работы на тему «Алгоритмы обхода особенности при численном решении уравнений движения тел относительно центра масс в атмосфере при действии возмущений»

Том 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, А—сопряженный кватернион;

Сг

iНе можете найти то, что вам нужно? Попробуйте сервис подбора литературы.

спи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 а

ю

в некоторое положение, характеризуемое углами а ! и у, следующим образом: вначале осуществляется поворот вокруг оси Ог на угол а>О, а затем плоскость (п, и) —плоскость угла атаки — поворачивается вокруг вектора 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 г.

i Надоели баннеры? Вы всегда можете отключить рекламу.