CC BY
108
численный расчет продольного демпфирования тела вращения малого удлинения при сверхзвуковом обтекании
А.Ю. ГАЛАКТИОНОВ, доц., МГУЛ, канд. техн. наук(1)
(1)ФГБОУ ВПО «Московский государственный университет леса» 141005, Московская обл., г. Мытищи-5, ул. 1-я институтская, д. 1, МГУЛ
Представлена математическая постановка задачи о свободных колебаниях тела в сверхзвуковом потоке воздуха и численный метод ее решения, позволяющий определить как стационарные, так и нестационарные аэродинамические характеристики. Основное внимание сосредоточено на модели спускаемого марсианского аппарата как на перспективном объекте исследования. Предметом исследования стали нестационарные характеристики модели тела вращения малого удлинения, которые были получены численно на ЭВМ с использованием программы решения полных уравнений Навье-Стокса (уравнения были записаны в подвижной системе координат, что потребовало рассмотрения вектора-столбца источниковых членов), разработанной автором. Зависимости производных демпфирования аэродинамического коэффициента момента тангажа по угловой скорости от угла атаки и числа Маха (интервал изменения 1,1 < М<» < 1,5) определены для ламинарного режима обтекания (ReD = 1,5-106) тела вращения малого удлинения. Достоверность полученных на персональной ЭВМ данных, представляющих интерес для десантируемых экспедиций к Марсу, подтверждена расчетами на различных сетках (в расчетах использовались структурированные расчетные сетки, координаты узлов которых рассчитывались заранее по алгебраическим зависимостям) и сравнением с известными экспериментальными данными. Показано влияние граничных условий на результаты численных расчетов, а также выход решения на устойчивый предельный цикл, соответствующий автоколебательному режиму. Для различных значений числа Маха набегающего потока, как для определяющего параметра задачи, определены количественные характеристики негативного автоколебательного режима, соответствующего случаю потери устойчивости.
Ключевые слова: уравнения Навье-Стокса, динамика вращательного движения, сопряженная задача, численные методы, аэродинамическое демпфирование, тело вращения малого удлинения, сверхзвуковое обтекание
Создание новых ракетно-космических носителей (РКН) может существенно снизить стоимость вывода полезной нагрузки на околоземную космическую орбиту и сделать более доступными различные научно-технические проекты, в частности экспедиции к Луне, Марсу и Венере. Учитывая высокую разряженность атмосферы Марса, стремление к максимизации тормозного аэродинамического импульса СА (С^^-тах) и аэродинамической устойчивости (как статической, так и динамической), а также ограниченные возможности по физическому воспроизведению марсианского обтекания СА в наземных аэродинамических установках [1], важной и актуальной проблемой становится разработка численных методов и математических моделей, позволяющих определять упомянутые газодинамические (аэродинамические) характеристики СА на ЭВМ.
Наряду с отмеченными требованиями (критериями выбора) к аэродинамическим схемам (аэродинамическим коэффициентам) СА, как десантируемым, так и проникающим в поверхность планеты, при
рассмотрении задачи об определении аэродинамического демпфирования необходимо отметить ограничение на амплитуду колебаний СА (|аА| < 5 - 10° [2]) и физическую возможность возникновения антидемпфирования у тел вращения малого удлинения на отмеченных режимах [3].
В качестве объекта настоящего исследования было выбрано тело вращения малого удлинения [3-5] с центровкой хцт~0,5 (рис. 1), совершающее свободные колебания в сверхзвуковом потоке воздуха при малых сверхзвуковых числах Маха (1,1 < Мш < 1,5). Предметом исследования стали аэродинамические характеристики, определяющие устойчивость СА, и в первую очередь его производные коэффициента момента тангажа по угловой скорости: т а+Ш2 = т в.
Целью работы стала демонстрация возможностей адаптированной к ЭВМ математической модели, позволяющей определять как статические, так и динамические аэродинамические характеристики летательных аппаратов.
Задачами исследования, направленными к достижению упомянутой цели, были:
1. Формализация математической постановки сопряженной задачи динамики полета и нестационарной аэродинамики колеблющегося в сверхзвуковом потоке тела вращения малого удлинения.
2. Численный метод и алгоритм решения отмеченной сопряженной задачи на ЭВМ как результат интегрирования уравнений На-вье-Стокса в подвижной системе координат.
3. Обоснование достоверности полученных данных (как внутренняя проверка расчетом на различных сетках - верификация; так и сравнение с известными экспериментальными данными - валидация).
4. Описание особенностей, полученных численным расчетом аэродинамических характеристик выбранной модели СА, связанных с его аэродинамической устойчивостью.
Численный расчет нестационарных аэродинамических характеристик методом свободных колебаний рассматривался как решение аэродинамически сопряженных задач:
- задачи решения уравнений Навье-
Стокса в нестационарной постановке
5ст да дЬ 8с „ /1Ч
— + — + —I-— = Н; (1) дt дх ду дг
- задача математического моделирования колебаний модели СА в продольной плоскости симметрии относительно центра тяжести под действием сил инерции и аэродинамических сил, в общем случае переменных по времени и по углу атаки
,_с12 Мг
¿Й2 ~ Л ~ Ь
ф = ф+ ф ш
ф = ф + фг# Мг= ^(р + хУЦ
ц
1г-момент инерции.
Уравнения Навье-Стокса (1) решались методом конечного объема [7] с расщеплением потоков от вязких и невязких членов. Поток от невязких членов определялся по методу С.К. Годунова как решение задачи Римана о распаде произвольного разрыва. Поток от вязких членов рассчитывался как решение системы линейных алгебраических уравнений, полученных разложением ком-
(2)
понент тензора скоростей деформации в ряд Тейлора.
Уравнения Навье-Стокса решались в подвижной системе координат (ОХУТ) [6], связанной с телом (рис. 1). Неподвижная система координат (ОХ0У010) была связана с подвижной через угол ф, совпадающий с углом атаки (а). Кориолисовы и переносные ускорения были представлены как источниковые члены в уравнениях (1) Нт = (0 pwx pwy pwz 0), w = + вт + ш-шт
Знание распределения давления и тензора вязких напряжений по поверхности тела (2) позволило установить суммарные силы и моменты, действующие на ЛА, что дало возможность определить угловое ускорение по известному моменту инерции. Далее численно рассчитывались угловая скорость и угол атаки.
Выполняя последовательно представленный алгоритм, можно итерационно проводить вычисления угла атаки и аэродинамических характеристик по времени, по аналогии с методом свободных колебаний в аэродинамической трубе.
На рис. 2(а) представлены зависимости угла атаки прямого и обратного конусов от времени, полученные описанным методом свободных колебаний. Прямой конус колебался относительно центра тяжести, расположенного на 0,584 доли длины модели от ее носка. Колебания обратного конуса математически моделировались при центровке 0,3. Начальный угол атаки для прямого конуса был выбран порядка 10°, начальный угол атаки обратного конуса соответствовал значениям 7 и 20°. При различных значениях начального угла атаки обратный конус выходил на угол атаки порядка 15°, что являлось косвенной внутренней проверкой (верификацией) работоспособности численного метода. Прямой конус выходил на угол порядка 0° в процессе колебаний.
Интересной особенностью рассмотренных на рис. 2(а) колебательных процессов было различие примерно в 10 раз демпфирующих характеристик прямого и обратного конусов.
Учитывая, что при малых сверхзвуковых скоростях существенным может оказать-
18
13
Условная граница передней и задней
Рис. 1. Расчетная область. Внешняя (Q^ и внутренняя (Q2) границы. Подвижная (OXYZ) и неподвижная (OX0Y0Z0) системы координат Fig. 1. The computational domain. External (Qj) and internal (Q2) borders. Mobile (OXYZ) and fixed (OX0Y0Z0) coordinate system
5 4 3 2
а) б)
Рис. 2. Зависимости углов атаки от математического времени (числа итераций) осесимметричных моделей: а) для конических моделей, полученные для различных начальных условий; б) для тела вращения малого удлинения на различных сетках (alfa_1 и alfa_2 - 15 точек в окружном направлении; alfa_3 - 30 точек в окружном направлении) и с различными граничными условиями (alfa_1 - «10 точек» на переднюю полусферу; alfa_2 и alfa_3 - «20 точек» на переднюю полусферу) Fig. 2. Dependence of the angle of attack of the mathematical time (number of iterations) of axisymmetric models: а) for bevel models obtained for different initial conditions; б) For a body small elongation rotating at different grids (alfa_1 and alfa_2 - 15 points in the circumferential direction; alfa_3 - 30 points in the circumferential direction) and with different boundary conditions (alfa_1 - «10 points» on the forward hemisphere; alfa_2 and alfa_3 -«20 points» on the front hemisphere)
8
3
0,30 0,25 0,20 0,15 0,10 0,05 0
0 2 4 6 8 10 12 14 16 18 M
■ - результаты экспериментальных исследований в аэродинамической трубе У-6 ЦНИИмаш [8]; — - результаты численных расчетов
Рис. 3. Зависимости аэродинамического коэффициента демпфирования прямого конуса от числа Маха набегающего потока: ■ - результаты экспериментальных исследований в аэродинамической трубе У-6 ЦНИИмаш [8]; ◊ - результаты численных расчетов Fig. 3. Dependence of the aerodynamic damping ratio on a direct cone on the Mach number of the incident flow: ■ - the results of experimental research in the wind tunnel У-6 TsNIIMash [8]; ◊ - the results of the numerical calculations
ся влияние граничных условия, на рис. 2(б) представлены зависимости углов атаки тела вращения малого удлинения (рис. 1):
- для случая воспроизведения невозмущенного потока в части передней полусферы, при постановке т.н. мягких граничных условий (дАдп = 0) во всей остальной внешней границе (а^а_1 - 10 точек вдоль продольного контура внешней границы);
- для случая воспроизведения невозмущенного потока фактически вдоль всего контура внешней границы передней полусферы (а^а_2 и а^а_3 - 10 точек - 10 точек).
Здесь можно отметить, что в процессе решения нестационарной задачи возмущения от внешней границы могут существенно снижать демпфирующие характеристики исследуемой модели и частоту моделируемого процесса. Это явление может быть в некоторой степени аналогично влиянию стенок и возмущений аэродинамической трубы при проведении соответствующих физических экспериментов.
С другой стороны, сравнение переходных процессов на рис. 2(б), полученных на различных расчетных сетках (15 точек (а^а_2) и 30 точек (а^а_3) в окружном, поперечном оси симметрии модели направлении соответственно) частота колебаний остается примерно постоянной, хотя с измельчением сетки демпфирование снижается и наблюдаются
тенденции выхода на устойчивый предельный цикл (рис. 4). В части отмеченного процесса упомянутое явление может быть объяснено более высоким значением демпфирования у граненых и плоских (например, крылья) тел, чем у их осесимметричных аналогов. Данное замечание особенно важно при оценке использования теории Ньютона (где конус и клин имеют одинаковые аэродинамические коэффициенты) и метода искривленных тел в случае определения демпфирующих характеристик [8].
В дальнейшем для подтверждения достоверности результатов расчетов были определены значения коэффициента т р в зависимости от числа Маха набегающего потока от 2 до 20. На рис. 3 приведены значения отмеченной расчетной зависимости и экспериментальные данные [8], полученные в аэродинамической трубе У-6 ЦНИИмаш методом свободных колебаний. Максимальное расхождение расчетных и экспериментальных данных не превосходит 25 %, что является хорошим результатом, учитывая сложность и систематические погрешности определения нестационарных характеристик ЛА расчетными и экспериментальными методами.
Расчетная зависимость (рис. 3) имеет максимум в интервале чисел Маха от 4 до 6 и монотонный нисходящий участок на гиперзвуковом интервале от 8 до 20, что поз-
alfa=f(w) 4 alla
3 ---\
/у 1 \ \
fi 0 1 w
-0J03I -0,002 -0,001 0,001 0,002 р,003/
\\ -2 y y
____
0
-0,002 -0,004 -0,006 -0,008 -0,010 -0,012 -0,014 -0,016
а)
б)
Рис. 4. Графические зависимости, демонстрирующие выход численного решения нестационарной задачи на режим «устойчивый предельный цикл»: а) зависимость угла атаки от угловой скорости движения тела вращения малого удлинения (в окружном направлении расчетная сетка содержит 30 точек, на переднюю полусферу отведено - 20 точек); б) зависимости коэффициента аэродинамического момента тангажа от угла атаки, полученные для стационарной и нестационарной задачи (М„ = 1,5) Fig. 4. Graphics showing the output depending on the numerical solution of the problem of non-stationary mode «stable limit cycle»: а) the dependence of the angle of attack of the angular speed of rotation of the body of small elongation (in the circumferential direction of the computational grid contains 30 points allocated to the forward hemisphere - 20 points); б) dependence of the coefficient of aerodynamic pitch moment on the angle of attack, prepared for the steady and unsteady problem (Мю = 1,5)
0,08 i
0,04
-0,04 J
mz_b —*— mz_b_15 -■- mz_b_12 mz_b_12
3 4
1 alfa
а)
48 a, deg
б)
12 16
Рис. 5. Зависимости коэффициентов демпфирования от угла атаки, полученные методом свободных колебаний: а) численным расчетом для малых сверхзвуковых чисел Маха (1,1, 1,2 и 1,5); б) экспериментально (М = 1,76) [3] Fig. 5. Dependencies damping coefficients on the angle of attack, obtained by the method of free oscillations: а) numerical calculations for low supersonic Mach numbers (1,1, 1,2 and 1,5); б) experimentally (M = 1,76). [3]
воляет распространять результаты трубных испытаний, полученные в гиперзвуковых установках, на участок натурного полета ЛА. При этом необходимо отметить, что с ростом гиперзвуковых чисел Маха наблюдается уменьшение демпфирующих характеристик ЛА, а это может привести к потере
динамической устойчивости на отдельных режимах.
В качестве самостоятельного подтверждения результатов расчета в части возможности выхода решения на устойчивый предельный цикл, соответствующий автоколебаниям, на рис. 4а приведена зависимость угла атаки
0
0
от угловой скорости (первой производной угла атаки по времени), а также зависимость коэффициента аэродинамического момента от угла атаки (рис. 4б). На отмеченном рисунке мо-ментные зависимости представлены для стационарного и нестационарного случаев. Для нестационарного случая можно выделить два участка цели аэродинамического гистерезиса: от 0 до 1,5° - антидемпфирование и от 1,5 до 3° - демпфирование.
Для количественной оценки частных производных момента тангажа по угловой скорости вращения использовалась разница между значением момента тангажа в стационарном расчете и в нестационарном. Таким образом, описанная методика позволила определить зависимости коэффициентов демпфирования от угла атаки для чисел Маха: 1,1, 1,2 и 1,5, соответственно (рис. 5а). Определенным качественным подтверждением правильности полученных результатов является зависимость коэффициентов демпфирования от угла атаки, полученная экспериментально в Центре Арнольда (США) и Институте фон Кармана (Бельгия) для СА «Викинг» [3], имевшего близкую к рассматриваемой модели компоновку и не имеющего заднего цилиндрического участка. В отечественной литературе [8] принято обозначение аэродинамического коэффициента демпфирования -тгр, в странах НАТО [3] - Cq, иногда с отличным знаком (рис. 5).
На рис. 5(а) можно отметить, что по мере приближения к околозвуковому режиму демпфирующие характеристики рассматриваемой модели в окрестности нулевого угла атаки уменьшаются, здесь преобладает антидемпфирование. Последнее объясняет режим автоколебаний с существенно нелинейным коэффициентом «демпфирования», зависящим от параметра (в данном случае от угла атаки).
Дальнейшим направлением настоящего исследования могут быть
- оценки влияния числа Рейнольса как параметра, определяющего перенос результатов трубных испытаний на режим натурного, например, марсианского полета;
- определение величины вклада различных участков рассмотренной модели и
их относительных удлинений в демпфирующие характеристики СА рассмотренной формы.
Заключение
Описаны элементы математической методики решения сопряженной аэродинамической задачи о свободных колебаниях тела вращения в сверхзвуковом потоке воздуха, что дало возможность провести математическое моделирование нелинейных колебательных процессов. Для описанного тела вращения определены демпфирующие аэродинамические характеристики, позволившие установить возникновение предельного цикла при числах Маха М^ = 1,1, 1,2 и 1,5, амплитуда которых растет по мере приближения к околозвуковому режиму. Достоверность полученных данных подтверждена сравнением с экспериментальными результатами, расчетами на различных сетках, с различными граничными условиями и с различными начальными условиями («сброс с различных установочных углов атаки»).
Библиографический список
1. Лукашевич, В.П. Космические крылья / В.П. Лукашевич, И.Б. Афанасьев. - М.: Лента странствий, 2009. - 496 с.
2. Козловский, В.А. Экспериментальное определение в аэродинамических трубах методом свободных колебаний характеристик демпфирования спускаемых в атмосфере планет аппаратов / В.А. Козловский // Космонавтика и ракетостроение, 2005. - Вып. 1(38).
3. Bob L. Useltom and Arthur R. Wallace, Damping-in-pitch and Drag Characteristics of the Viking Configuration at Mach Number from 1.6 through 3, AEDC-TR-72-56, May, 1972.
4. Краснов, Н.Ф.. Аэродинамика тел вращения / Н.Ф. Краснов. - М.: Машиностроение, 1964. - 572 с.
5. Петров, К.П.. Аэродинамика тел простейшей формы / К.П. Петров. - М.: Факториал, 1998. - 432 с.
6. Галактионов, А.Ю. Численное моделирование пространственного взаимодействия боковой струи с набегающим сверхзвуковым потоком / А.Ю. Галактионов // Космонавтика и ракетостроение. - № 2(39). - 2005. -С. 49-58.
7. Липницкий, Ю.М. Численное моделирование нестационарных аэродинамических характеристик затупленного конуса в рамках полных уравнений Навье-Стокса / Ю.М. Липницкий, А.Ю. Галактионов // Космонавтика и ракетостроение, 2006. - № 3. - С. 23-28.
8. Липницкий, Ю.М. Нестационарная аэродинамика баллистического полета / Ю.М. Липницкий, А.В. Красильщиков, А.Н. Покровский и др. - М.: ФИЗМАТЛИТ, 2003. - 176 с.
NUMERICAL CALCULATION OF THE PITCH DAMPING FOR THE ROTATION BODY OF SMALL LENGTH AT THE SUPERSONIC FLOW
Galaktionov A.Yu., Assoc. Prof. MSFU, Ph.D. (Tech.) (1)
(1)Moscow Forest State University (MSFU), 1st Institutskaya st., 1, 141005 Mytishchi, Moscow reg., Russia
The mathematical task of the free fluctuations of the body in the supersonic air flow has been shown and a numerical method of its decision has been given to define stationary and non-stationary aerodynamic characteristics.The main attention was given to the model of the Mars mission as the object of research. The subject of the research were non-stationary aerodynamic characteristics of small length which were obtained by the computer technologies with the author's program of the solution of Navier-Stoks equations (they were written in the moving coordinate system with the Vector-column of the source parts). The equations (with Mach number in the interval from 1,1 to 1,5) for the laminar case (ReD=1,5-106) have been obtained. The reliability of the data received on the personal computers for the Mars expeditions is confirmed by accounts on various grids (algebraic grids of different structure were used for the solution) and comparison with the known experimental data. It was shown that the influence of the boundary conditions on the results of the numerical calculations is significant, and the solution can be found by a limited cycle which corresponds to the vibration case. For the different Mach numbers, as one of the main factors of the tasks, quantitative characteristics were taken.
Keywords: Navier-Stokes equations, dynamic of rotation, aerodynamic damping, supersonic flow, the rotation bodies of small lengthening, connected tasks
References
1. Lukashevich V.P., Afanas'ev I.B. Kosmicheskie kryl'ya, [Space wings]. Moscow: Lenta stranstviy, 2009, p. 496.
2. Kozlovskiy V.A. Eksperimental'noe opredelenie v aerodinamicheskikh trubakh metodom svobodnykh kolebaniy kharakteristik dempfirovaniya spuskaemykh v atmosfere planet apparatov [Experimental definition of unsteady aerodynamic parameters of reentry body for atmosphere planets]. Kosmonavtika i raketostroenie, 2005, V. 1(38).
3. Bob L. Useltom and Arthur R. Wallace, Damping-in-pitch and Drag Characteristics of the Viking Configuration at Mach Number from 1.6 through 3, AEDC-TR-72-56, May, 1972.
4. Krasnov N.F. Aerodinamika tel vrashcheniya [Aerodynamic of rotation objects]. Moscow: Mashinostroenie, 1964. p 572.
5. Petrov K.P. Aerodinamika telprosteyshey formy [Aerodynamic of simple objects]. Moscow: Faktorial, 1998. p. 432.
6. Galaktionov A.Yu. Chislennoe modelirovanie prostranstvennogo vzaimodeystviya bokovoy strui s nabegayushchim sverkhzvukovym potokom [Numerical simulation of 3D-interaction size-jet and supersonic flow]. Kosmonavtika i raketostroenie. № 2(39). 2005. p. 49-58.
7. Lipnitskiy Yu.M., Galaktionov A.Yu. Chislennoe modelirovanie nestatsionarnykh aerodinamicheskikh kharakteristik zatuplennogo konusa v ramkakh polnykh uravneniy Nav'e-Stoksa [Numerical model of unsteady aerodynamic of cone with Navier-Stoks equations]. Kosmonavtika i raketostroenie, 2006. № 3. pp. 23-28.
8. Lipnitskiy Yu.M., Krasil'shchikov A.V., Pokrovskiy A.N., Shmanenkov V.N. Nestatsionarnaya aerodinamika ballisticheskogo poleta [Unsteady aerodynamic of ballistic flight]. Moscow: FIZMATLIT, 2003. p 176.
