Экспериментально-расчетный метод определения нелинейных аэродинаических характеристик острого конуса Текст научной статьи по специальности «Математика»

2013

ВЕСТНИК ТОМСКОГО ГОСУДАРСТВЕННОГО УНИВЕРСИТЕТА Математика и механика

№ 4(24)

МЕХАНИКА

УДК 531.554

В.И. Биматов, Н.В. Савкина

ЭКСПЕРИМЕНТАЛЬНО-РАСЧЕТНЫЙ МЕТОД ОПРЕДЕЛЕНИЯ НЕЛИНЕЙНЫХ АЭРОДИНАИЧЕСКИХ ХАРАКТЕРИСТИК ОСТРОГО КОНУСА

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

Ключевые слова: траектория движения, аэродинамические характеристики, летательный аппарат.

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

В определенном (узком) смысле целью решения аэродинамической задачи является получение зависимостей коэффициентов аэродинамических сил и моментов от формы ЛА и условий его обтекания потоком газа. Решение этой задачи осуществляется экспериментальными и теоретическими методами. Экспериментальные методы, в свою очередь, разделяются на исследования в аэродинамических трубах и на аэробаллистических установках [1].

За десятилетия накоплен большой опыт исследования АДХ по данным аэробаллистических испытаний, однако лишь к 80-м - 90-м годам для решения этой задачи стали применять методы теории построения нелинейных моделей [2-5]. Совсем неудовлетворительно обстоит дело с привлечением методов решения обратных задач, широко используемых в ядерной физике, тепломассообмене [6, 7]. Для задач оценивания АДХ летательных аппаратов, как будет показано ниже, такой подход является единственно возможным.

Стремительное развитие вычислительной техники в последнее десятилетие позволяет эффективно реализовать методы решения обратных задач и внедрять их в практику экспериментальных исследований. Одной из проблем динамики летательных аппаратов является эффективное оценивание аэродинамических характе-

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

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

В настоящей работе анализируются численные подходы к реализации обратной задачи траекторной баллистики.

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

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

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

dx

— = V cos 9 cos Y, x(0) = x0,

dt 0

= V sin 9, y (0) = yo,

dt

— = -V cos 9 sin Y, z(0) = z0; (1)

dt

^ = -KiV2Cx - gsin9, V(0) = V,, dt

d9

— = K1VCY - gcos 9, 9(0) = 90, dt

d Y

---= K1VCZ / cos 9, Y(0) = Y0; (2)

dt

d ю71

—— = K2mZ1, raZ1(0) = raZ10, dt

d roY1

- = K2mY—, ®Y 1(0) = ®Y 10 ; ()

dt

ddr=aZi’ Q(°)=v dd!r=“"Sñ, v(0)=v°; ()

dt dt cos 0

sin S = sin 0 cos a cos p + cos 0 sin a,

sin y = sin Y cos p + cos Y sin p cos 0. (5)

Здесь K =pS /2m, K2 = pSl / 2 Jz1. Обозначения приведены в конце текста.

Аэродинамические коэффициенты CX, CY, CZ, mZ1, mY1 будем считать функциями углов a, р, 5, угловых скоростей "Z1, "Y1, числа Маха M. Зависимостью от

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

Cx = Cx °(M) + Cx 2 (M )52 + Cx 4 (M )54 +... (6)

CY = CY1(M)a + CY3(M)a3 +... (7)

CZ = CZ1(M)P + CZ3(M)P3 +... (8)

mZ = mZ1 (M)a + mz3 (M)a3 + m""z ; (9)

mY = mY1(M)P + mY3(M)p3 + m""Y . (10)

Коэффициенты разложений (6) - (10) будем считать линейными функциями числа Маха, что вполне допустимо на баллистических трассах небольшой протяженности для начальных значений числа Маха М° > 1,5:

Cx o (M) = Cx ° (M°) + CM, (M° - M). (11)

2. Алгоритм расчета аэродинамических коэффициентов

Прежде чем приступить к обсуждению технологии расчета параметров CJXІ, С^, mJZІ, сформулируем общий подход, учитывающий специфику задачи.

В векторной форме система уравнений (1) - (5) имеет вид

1(0 = f (Р, с, t), p(0) = р,, t 6 [0, tk ]. (12)

Здесь р(0 = {х(0, 7(t), 2(0, К(t), 9^), Y(0, ю7 (t), (t), ), у(0} - вектор решений

системы (1) - (5), с = {Сх 0, С^0,..., т(, х0,..., у 0} - вектор неизвестных параметров (коэффициентов разложений АДХ и начальных условий). Исходной информацией для оценки неизвестных коэффициентов вектора с является система измеренных в опыте дискретных значений кинематических параметров движения тела. При таких посылках математическая постановка задачи определения АДХ представляет собой задачу оценивания вектора с коэффициентов их полиномиальных разложений по известным из эксперимента значениям параметров движения р^) и является примером классической обратной задачи для дифференциальных уравнений [6].

Наиболее употребительные методы решения обратных задач по своей природе являются вариационными. Для задач вида (12) в качестве решения принимают вектор с , доставляющий минимум одному из функционалов:

р(p, f (p,с,t) = minip(x) - f (p,c,t)|; (13)

P( p, g (p, с, t) = min

X

p(x) - (p(Xo) + { f (p, с, t)dx)

X,

(14)

Р(Ри , Рр ) = №>11 Ри ) - Рр )||, (15)

С, Ро

где ри - совокупность экспериментальных (измеренных) данных о векторе рЦ), рр - рассчитанные значения параметров движения.

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

Предлагаемая методика расчета аэродинамических характеристик основана на использовании системы дифференциальных уравнений пространственного или плоского движения твердого тела. Для поиска минимума функционала (15) используется вещественный генетический алгоритм [9, 10], позволяющий находить глобальный экстремум для мультимодальных функций. Этот метод сочетает в себе естественный отбор среди строчных структур с упорядоченным (хотя в чем-то и случайным) обменом информацией. Будучи вероятностным, генетический алгоритм тем не менее не является просто еще одним вариантом случайного поиска, поскольку при отборе новых точек с ожидаемыми более хорошими возможностями он эффективно используют предыдущую информацию. Для решения задач оптимизации в данной работе предлагается вещественный генетический алгоритм (ВГА), совмещающий в себе детерминистический и вероятностный подходы и основанный на механизмах природной селекции и генетики. Важной особенностью генетических алгоритмов является их робастность: они сходятся к глобальному оптимуму (что очень важно для задач, у которых целевая функция имеет несколько локальных экстремумов) и в отличие от классических градиентных методов оптимизации при их реализации не требуется сильных ограничений на гладкость целевой функции, и они позволяют находить оптимум даже для случая, когда целевая функция является разрывной. Данный алгоритм носит итерационный характер и имеет дело с обработкой популяций индивидуумов (или наборов испытаний), каждый из которых представляет потенциальное решение задачи и является вектором в пространстве поиска. Важным свойством генетического алгоритма является также его сравнительно легкая адаптация к параллельным компьютерам,

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

Остановимся более подробно на трех «генетических операторах» - селекции, скрещивании и мутации.

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

Скрещивание. Наиболее простым является одноточечное скрещивание - каждая выбранная таким образом пара строк скрещивается следующим образом: случайным образом выбирается положение точки сечения (целое число к в промежутке от 1 и Z—1 , где l - длина строки). Затем, путем обмена всеми элементами между позициями к+1 и l включительно, рождаются две новые строки. Пусть A = (у1, y2, у3, y4) и A' = (y'j, у'2, y'3, y'4) являются родителями, выбранными в процессе селекции. Тогда (считая, что случайно выбранная точка сечения находится после первого гена), они производят двух детей B = (yi, у'2, у'3, у'4) и B = (у'і, у2, у3, у4). После этого дети замещают родителей в промежуточной популяции P'.

Мутация. Использовалась неоднородная мутация, определенная Михалеви-чем. Если ген у, подвергается мутации, то его новое измененное значение у', выбирается внутри интервала [Min,, Max,] следующим образом:

у'і= у, + * (M - у, ) -§) ,

где s - случайное число из интервала [0,1], g - номер поколения, G - максимальное число итераций, b - параметр уточнения, M, случайным образом выбирается из множества {Min,, Max,}, где Min, и Max, - нижняя и верхняя границы возможного изменения значения переменной у,. Такая адаптивная мутация позволяла соблюдать в процессе реализации ГА (эволюции) необходимый баланс между двумя разномасштабными изменениями (мутациями) генов, так как на первоначальных шагах алгоритма в основном преобладали крупномасштабные изменения (обеспечивающие широкую область поиска), в то время как на заключительном этапе (за счет уменьшения масштаба мутаций) происходило уточнение решения.

3. Целевая функция

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

N _ _ _ N1

Ф(С) = х (Хизм - Хр )W (Хшзм - Хр)+<£ (p )2, (16)

,=1 j=1

Хр = {Хр, ,7р,-, 2рі, 9р,-, } - вектор рассчитанных значений параметров траекто-

ри^ Хизм = {Хизм, Лзм,, Zизм,, &изм,, Уизм, } - вектор измеренных в опыте значений

параметров траектории, W - диагональная матрица весов измерений с элемента-

ми вида 1/стХ , где стХ - среднеквадратические погрешности измерений, £ -

Х1 Х1

параметр регуляризации.

Поставленная задача оптимизации (минимизации функционала со случайными ошибками в исходных данных) является стохастической. Если в детерминированных задачах отыскания локального экстремума стохастические методы имеет смысл применять тогда, когда задача имеет очень большую размерность, то в задачах отыскания глобального экстремума стохастические методы не имеют альтернатив. Подходы поиска решения задачи, реализуемые у авторов [2,4,5], основаны на анализе локальных свойств изучаемой функции и приводят к нахождению лишь локального экстремума. Для отыскания глобального экстремума нужны принципиально другие методы - стохастические.

Минимизацию функционала (16) будем осуществлять методом прямого поиска с использованием вышеизложенного генетического алгоритма.

Статические характеристики и погрешности оценивания параметров АДХ будем исследовать методом имитационного моделирования, то есть проведем анализ динамического процесса с помощью статистических испытаний по следующему алгоритму. Проводим серию расчетов с различными значениями измеряемых параметров, лежащими в интервале [ Хизм -стизм, Хизм -стизм ]. Из полученной выборки значений коэффициентов АДХ определяем средние значения и выборочную дисперсию. По критерию Фишера находим доверительные интервалы оценки коэффициентов АДХ.

4. Оценка значимости параметров разложения аэродинамических характеристик и адекватности модели

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

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

Строго говоря, в случае нелинейных аэродинамических характеристик адекватность их описания не может быть установлена однозначно только на основе анализа траекторных данных без привлечения априорных физических сведений. Объясняется это тем, что, выбирая модель, исследователи, как правило, контролируют ее пригодность лишь по отклонению расчетной траектории от измеренных координат по тем или иным критериям [2, 4, 5]. Но поскольку параметры траектории являются функционалами от правых частей уравнений движения, вполне очевидно, что условиям выбранного критерия может удовлетворять множество параметров правых частей.

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

5. Результаты оценивания аэродинамических характеристик

Изложенные в предыдущих пунктах приемы были применены при обработке траекторных данных, полученных на баллистической установке ФТИ РАН [11].

С учетом точности юстировки оптических систем, погрешности измерения времени, линейных координат и углов оцениваются величинами = 0,6 мкс,

стх = 0,4 мм, ст7 = 0,26 мм, ста = 0,3°.

Располагая этими оценками, авторы обработали траекторные данные острого конуса с углом при вершине 29^ = 30° с диаметром основания 30 мм, массой и 45 г, квадратом поперечного радиуса инерции около 185 мм2, движущегося с начальной скоростью около 680 м/с (М=2). Эксперименты проводились в воздухе при нормальных условиях. Максимальное значение угла тангажа достигало 40°.

Аэродинамические характеристики были представлены в виде

Учитывая особенности баллистической трассы ФТИ РАН, в системе (1) - (4) оставим только уравнения плоского движения острого конуса, а уравнения (5) представим в виде а = 9 - 9.

На рис. 1 приведены зависимости коэффициента аэродинамического момента и оценки относительной погрешности, найденные по методике ФТИ РАН (кривая 1) и по методике данной работы (кривая 2) для групп из трех и четырех опытов. При обработке трех опытов (а > 28°) наблюдается слабая нелинейность статического момента по углу атаки. Относительная погрешность зависимости Ст(а) меняется слабо, оставаясь ниже 3 %. «Завал» моментной характеристики проявляется лишь в четвертом эксперименте при а > 30°. Оценка погрешностей с использованием дисперсионных матриц (1) и методом ГА (2) дали близкие результаты.

тг = т21а + т2 3а3 + т2т2.

0

10

20

30

а, град.

0

-0,02

§ -0,04

-0,06

16

1 - по дисперсной матрице

2 - с помощью генетического алгоритма

20

Рис. 1. Статический аэродинамический момент острого конуса

На рис. 2 даны зависимости коэффициента подъемной силы в трех и четырех опытах. Поскольку при а < 20° подъемная сила линейна, для описания ее «завала» при а > 30° более приемлемой оказалась пятая степень угла атаки по сравнению с третей (она дала меньшие погрешности). Здесь, так же как и для момента, отмечается существенное возрастание погрешности к концу интервала а.

Рис. 2. Зависимость коэффициента подъемной силы

На рис. 3 показаны зависимости коэффициента сопротивления, также соответствующие двум группам опытов. Как и в случае подъемной силы, более удачной оказалась модель с «пропущенной» (четвертой) степенью угла атаки. Коэффициент сопротивления найден с наилучшей точностью и типичное для конца интервала возрастание погрешности здесь не выглядит катастрофическим (около 3 %).

Рис. 3. Зависимость коэффициента сопротивления

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

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

Список основных обозначений в уравнениях движения

х, y, z - оси лабораторной системы координат,

Сх, CY, cz - коэффициенты лобового сопротивления , подъемной и боковой силы, mz, mY - коэффициенты момента тангажа и рыскания, g - ускорение силы тяжести, р - плотность воздуха, m, l - масса и длина тела,

S - площадь миделя тела,

М - число Маха, t - время,

V - модуль вектора скорости тела,

9, Y - угол наклона траектории и курсовой угол,

S, у - угол тангажа и рыскания,

а, р - угол атаки и скольжения.

5 - пространственный угол атаки: cos 5 = cos а cos p.

ЛИТЕРАТУРА

1. Баллистические установки и их применение в экспериментальных исследованиях / Н.А. Златин, А.П. Красильщиков, Г.И. Мишин, Н.Н. Попов; под ред. Н.А. Златина и Г.И. Мишина. М.: Наука, 1974. 344 с.

2. Chapman G.T., KirkD.B. A Method for extracting aerodynamic coefficients from free-flight data // AIAA J. 1970. 8(4). P. 753-758.

3. Биматов В.И. Определение аэродинамических характеристик тел сложных форм по данным аэробаллистических испытаний // Материалы II Всесоюзной конференции по методам аэрофизических исследований. Ч. 2. Новосибирск, 1979. С. 179-181.

4. Менде Н.П. Об одном методе определения нелинейных аэродинамических сил и моментов // Физико-газодинамические баллистические исследования / под ред. Г.И. Мишина. Л.: Наука, 1980. С. 200-224.

5. Костров А.В. Движение асимметричных баллистических аппаратов. М.: Машиностроение, 1984. 272 с.

6. Тихонов А.Н., Арсенин В.Я. Методы решения некорректных задач. М.: Наука, 1979. 288 с.

7. Алифанов О.М. Идентификация процессов теплообмена летательных аппаратов. М.: Машиностроение, 1979.

8. Рао С.Р. Линейные статистические методы и их применение. М.: Наука, 1968.

9. Казаков В.Ю., Пейгин С.В., Тимченко С.В. Оптимизация траектории входа в атмосферу земли по интегральному тепловому потоку // ПМТФ. 2000. Т. 41. № 4. С. 112-121.

10. Биматов В.И., Тимченко С.В. Применение генетических алгоритмов к решению задач оптимизации в гиперзвуковой аэродинамике // Фундаментальные и прикладные проблемы современной механики: материалы конф. Томск: Изд-во Том. ун-та, 2002. С. 209-210.

11. Басаргин И.В., Дементьев И.М., Мишин Г.И. Полигон для аэродинамических исследований // Аэрофизические исследования сверхзвуковых течений / под ред. Ю.А. Дунаева. М.-Л.: Наука, 1967. С. 168-178.

Статья поступила 13.03.2013 г.

Bimatov V. I, Savkina N. V. EXPERIMENTAL AND COMPUTATIONAL METHOD OF DETERMINING NONLINEAR AERODYNAMIC CHARACTERISTICS OF A SHARP CONE. Nonlinear aerodynamic characteristics of a sharp cone are estimated using trajectory data for bodies freely flying in the atmospheric environment. To construct a calculation algorithm, the problem is analyzed as an inverse problem for differential equations. Results of the estimation of aerodynamic characteristics of a sharp cone at Mach numbers M = 2 are presented.

Keywords: motion trajectory, aerodynamic characteristics, airborne vehicle

BIMATOV Vladimir Ismagilovich (Tomsk State University)

E-mail: vbimatov@mail.ru

SAVKINA Nadezhda Valer’evna (Tomsk State University)

E-mail: pantera@ftf.tsu.ru