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

УЧЕНЫЕ ЗАПУСКИ Ц АГ И Т о м VI 197 5

№ 3

УДК 533.6.013.2.011.35:629.7.025.73

РАСЧЕТ ОБТЕКАНИЯ КОЛЕБЛЮЩЕГОСЯ ПРОФИЛЯ ТРАНСЗВУКОВЫМ ПОТОКОМ ГАЗА

С. И. Кузьмина

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

Методы установления успешно применяются для расчета обтекания осесимметричных тел и профилей трансзвуковым потоком газа [1, 2]. Так, в работе [1] методом Годунова получено решение стационарной задачи трансзвукового обтекания профиля под углом атаки; результаты расчетов хорошо согласуются с экспериментальными данными. В работе [3] этот метод применяется для расчета переходных процессов.

В настоящей работе конечноразностный метод Годунова с подвижной сеткой применяется к решению нестационарной задачи, в частности к задаче демпфирования колебаний профиля в околозвуковом потоке.

Постановка задачи. Рассматривается задача обтекания симметричного 8%-ного профиля с элероном невязким, нетеплопроводным совершенным газом. Предполагается, что профиль представляет собой жесткое тело, имеющее три степени свободы: вертикальное смещение центра жесткости h, поворот вокруг центра жесткости а и отклонение элерона относительно профиля р. Полагая h = h{t)=0, а = а(t) и заставляя элерон колебаться по гармоническому закону р = р0 sin с заданной частотой со, определим работу аэродинамических сил на движущемся элероне за полный период колебаний. Это позволит судить о том, демпфируются ли колебания при заданном режиме обтекания. Расчеты проводились для нескольких значений числа Моо набегающего потока: М«,=,0,845; 1,18; 1,69; 2,53; во всех рассчитываемых случаях полагалось [30=0,1;

(в = 0,3; Sh = w = -^- , где ш — частота в рад/с, b — хорда элерона,

«со

«оо—скорость звука, Sh — число Струхаля.

Метод расчета. Поставленную задачу будем решать методом Годунова [4], основанным на интегральных законах сохранения для нестационарных уравнений Эйлера. При рассмотрении задачи об обтекании колеблющегося тела воспользуемся расчетной сеткой, связанной с профилем [1] и перемещающейся вместе с ним согласно заданным законам h(t), a(t) и р(£) (фиг. 1,а). Введем две прямоугольные системы координат: неподвижную хОу и жестко связанную с профилем х0О'у0. Профиль с неотклоненным элероном расположим вдоль оси О' хо, начало координат О' поместим в носовой

Фиг. I

части профиля и будем строить расчетную сетку в этой системе координат. Способ построения расчетной сетки для профиля с закругленной передней и острой задней кромками показан на фиг.1, б; для профиля с острыми передней и задней кромками расчетная сетка изображена на фиг. \,г. Проведем N лучей и перенумеруем их, как показано на фиг. 1, б (г = 0, 1, 2, 3, . . . , 1\Г). Область позади хвостовой части элерона разбивается на частей, причем расстояния между лучами увеличиваются в геометрической прогрессии со знаменателем по мере удаления от элерона. Элерон и часть профиля до луча I = А/3 разбиваются соответственно на N^ — N1 и УУ3 — N2 равных частей. В случае закругленной передней кромки на носовой части профиля строится „веер“ лучей с углами наклона 0ог к оси О'х0, согласно формулам:

чт6 _ 3(/~лг4) д-1у4)з

&ш 0 1 2(ЛГ4-ЛУ 2 (Л^4—Л^3)3 ’

cos0O;= — sin2 80. при N3 < / < N — N3.

В остальной области расчета

sin 6,

01 —

- 1

при 0 < і < N3",

sin 60 і = 1 cos60 і = О

при N — yv3 і N.

cos 60 і == 0

При острой передней кромке луч i=Ns смещается против течения в отрицательную область оси О'хй. Это позволяет устанавливать желаемую величину минимального шага по времени х, так как при очень малом t требуется чрезмерно большое машинное время для расчета периода колебаний в трансзвуковом режиме. Каждый луч разбивается на М частей (/= 1, 2, . . . , М). Соединяя точки на лучах при ] ■= 1, получаем ломаную линию, которая аппроксимирует образующую профиля. Расстояния между точками разбиения на каждом луче увеличиваются по закону геометрической прогрессии со знаменателем q2 по мере увеличения номера /'.

Пусть в момент времени t0 обе системы координат совпадают, тогда в другой момент tф t0 координаты узлов сетки можно вычислить, согласно формулам переноса начала и поворота системы координат. Так, если x0lj, у0 i;- — координаты ij узла сетки в системе х0О'у0, то их значения в системе хОу будут

xl} = а + (b — a) cos а 4- (х0 и — Ъ) cos (а — Р) -bj/0 ц sin а;

УU = h + (а — b) sin а — (х0 и — b) sin (а — р) + уй и cos а, где а — расстояние от передней кромки профиля до центра жест-

Аналогично выпишем формулы преобразования углов в системе координат хОу:

Площади ячеек при преобразовании системы координат изменяются только для 0<г<;Д/2 и N—Л^<г<М Так, если 50^— площадь некоторой ячейки в системе х0О'у0 [т. е. при Р(0 = 0], то в системе хОу (в момент времени /) площадь этой ячейки будет

- Зная закон движения сетки в неподвижной системе координат» можно найти скорости движения узлов и границ ячеек сетки.. Рассмотрим ячейку, образованную I—1 и I лучами между у — 1 и у слоями. Границы ячейки пронумеруем так, как показано на фиг. 1, в. Обозначим <р,-, — угол между осью Ох и отрезком ломаной

в слое у — 1 между лучами г—1 и г; таким образом, согласно принятым обозначениям <р/, у_ 1 — угол между осью Ох и третьей

кости.

при

при N3 < і < N — N3;

sin 0г = cos a cos 6,- = — sin a

S(£) = S0cosfi(^).

границей ячейки. Скорость га движения границы по направлению, перпендикулярному к ней, можно вычислить по формулам

Скорости ^2 и вычисляются аналогично. При расчете распада разрыва и при написании законов сохранения для каждой ячейки необходимо учитывать скорость движения границ ячейки и изменение ее площади. Кроме того, вместо и и V — проекций скорости потока на оси Ох и Оу, в расчетах берем составляющие скорости вдоль границы ячейки и по нормали к ней.

Граничные условия- задаем на внешней границе расчетной области и на поверхности тела. Параметры сетки подбираем таким образом, чтобы обеспечить достаточное удаление границ области от профиля. Тогда на внешней границе при —Л^3 задаем

условия в набегающем потоке (рсо=1, рте=1, ^00=О, и.^— Моэаоэсоз а), на всех же остальных границах считаем, что значения газодинамических параметров в соседних ячейках, прилегающих с двух сторон к границе, равны. Расчеты проводились с такими значениями параметров сетки М, М, ци <72, что границы области располагались на расстоянии от профиля, равном трем его хордам. В такой расчетной области, как показывают численные результаты, влияние внешних границ на течение около профиля очень мало. На границе, которая является образующей тела, задаем условие равенства нулю нормального к границе компонента скорости с учетом скорости движения самой образующей. На линии у—1 позади элерона граничные условия специально не задаются; ячейки, прилегающие к этой линии, рассчитываем как внутренние ячейки.

Расчет трансзвукового обтекания профиля. Задачу об обтекании колеблющегося профиля будем решать в два этапа. Сначала определим распределение давления на жестком профиле в момент времени при заданном режиме течения. Используя полученное решение стационарной задачи в качестве начальных условий, найдем распределение давления на движущемся профиле к моменту времени t = где х —шаг повремени, определяемый из усло-

вия устойчивости разностной схемы. Учитывая скорости перемещения точек профиля, вычислим работу сил давления за время х. Полученные на этом шаге расчета значения давления р, плотности р и скоростей и я V снова используем в качестве начальных условий и определяем работу сил давления за время £, = £-(-т. Следующий шаг расчета дает значение работы за время = и т. д. Расчет продолжаем до тех пор, пока не получим интересующую нас величину работы за полный период Т колебания элерона.

В ходе решения задачи было исследовано влияние выбора расчетной схемы на точность результатов. Были построены три расчетные сетки, различающиеся общим числом ячеек, распределением лучей по профилю и удалением внешних границ от профиля. При Моо = 0,845 в этих трех случаях рассчитывалось течение

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

0,16 -дс (1 ~х).

Параметры сеток для этих расчетов приведены ниже, причем в первых двух случаях использовалась сетка, подобная изображенной на фиг. 1, б, а в третьем случае такая, как на фиг. 1, г. Здесь же

приведена величина шага по времени, выбираемого в расчетах с различными сетками.

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

Исследование демпфирования колебаний. Вопрос о демпфировании колебаний элерона при заданных числах БЬ и Моо рассматриваем с помощью следующей модели. Пусть элерон как жесткое тело совершает вынужденные гармонические колебания р (£) = розт*ш£.

м N М Л/з К т

20 104 12 27 47 52 0,22*10-2

25 130 18 38 63 65 0,45-10-2

25 130 15 29 65 - 0,83.10-2

Р/Рс

ЬО

\

У

V 1,

1

> N /

У

N Г

5ч

ч (,

/1

Т=0,22-10 1 г

~~ -- 0,45-ю; — 0,83-Ю

I I I I

0,2 0,4 0,6' 0,8 х

Фиг. 2

Предположим, что зависимость аэродинамических сил давления от времени может быть представлена в виде:

/?(/) = /> з1п (о>* + <Ю,

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

Известно, что работа, совершаемая потоком газа за каждый цикл гармонического движения элерона, положительна при sin']><0 и отрицательна при sin^>0

2Л Ь //Ч

A = dxdt = — b% mizp0 sin 'b.

о 0 ,

Это значит, что при sin <!О демпфирование отрицательно, а при sin положительно. На фиг. 3—7 приведены результаты численных расчетов. Характер изменения по времени давления на профиле при Моо = 0,845 показан на фиг. 3. Расчет проводился по сетке № 2, сплошной кривой показано давление на верхней поверхности элерона в точке х = 0,7 (ячейка 98), пунктирной линией — давление на нижней поверхности при л; = 0,7 (ячейка 33). Видно, что полученные кривые представляют собой периодические функции с периодом Т, но с некоторым запаздыванием относительно гармонического колебания, изображенного штрихпунктирной линией. Сообра-

жения о связи фазы ф со знаком демпфирования указывают на наличие отрицательного демпфирования в рассматриваемом случае. Однако поскольку функция давления носит явно негармонический характер, правильный ответ на вопрос об устойчивости колебаний можно получить, только интегрируя работу аэродинамических сил на профиле по времени. Такое интегрирование приводит к положительной величине работы сил давления (Л) на движущемся элероне за период колебания (кривая с обозначением М = 0,845 на фиг. 4). Зависимость величины работы от времени на движущемся с той же частотой элероне при М = 1,18, 1,69 и 2,53 изображена на фиг. 4. Здесь же для сравнения пунктирной линией нанесена работа при М = 2,5, вычисленная по линейной поршневой теории. Небольшое различие сравниваемых значений работы за период колебаний, вероятно, можно объяснить ошибкой численного метода, имеющего первый порядок точности. Расчеты обтекания профиля

при небольших сверхзвуковых скоростях (Moo=l,l8 и 1,69) дают положительные значения работы за период, но на порядок меньшие, чем при Моо = 0,845. В этом случае нужно применить более точные численные методы расчета.

Приведем некоторые результаты расчета трансзвукового (Моо —0,845) обтекания профиля. На фиг. 5 показано положение скачков уплотнения на нижней (сплошная линия) и верхней (пунктирная) поверхностях профиля в течение двух периодов колебания. Скачки уплотнения и связанные с ними зоны повышенного и пониженного давления перемещаются по поверхности тела таким образом, что восстанавливающий момент относительно шарнира элерона запаздывает по фазе. Величина запаздывания определяется, в основном, режимом обтекания (числами Ми и Sh). На фиг. 6, а изображены линии равных давлений около профиля в поле течения при Моо = 0,845, Р = 5,7° в момент времени Т/4, пунктирной линией показаны границы сверхзвуковых зон; распределение давления на верхней (сплошная линия) и нижней (штрихпунктирная линия) поверхностях профиля представлено на фиг. 6, б. Расположение сверхзвуковых зон на поверхности профиля в последовательные

Т Т 3 т

моменты времени 0, j,y, уи Г представлено на фиг. 7.

Полученные картины течения в каждый момент времени £>0 де совпадают с картинами течения около профиля с неподвижным элероном, отклоненным на соответствующий угол. Из фиг. 7 видно, например, что при 7= 7у2 и (=Т, когда Р(£) = 0, картины течения существенно отличаются от поля течения в стационарном случае при 1 = 0. Интенсивность скачка уплотнения зависит от скорости его перемещения по поверхности профиля [5]. При движении скачка вверх по потоку он будет иметь большую интенсивность, а при движении вниз по потоку — меньшую, чем интенсивность соответствующего неподвижного скачка.

Таким образом, результаты расчета подтверждают возможность возникновения автоколебаний элерона с одной степенью свободы в околозвуковом потоке. Результаты экспериментального исследования автоколебаний профиля, аналогичного рассчитанному9 приведены в работе [6]. Автоколебания элерона были получены в диапазоне чисел М = 0,91-*-0,95. Так как расчеты проводились при других значениях числа М и без учета вязкости, то можно говорить лишь о качественном совпадении полученной картины обтекания с экспериментом.

ЛИТЕРАТУРА

.1. Фонарев А. С. Расчет обтекания осесимметричных тел и несущих крыловых профилей трансзвуковым потоком газа. .Ученые записки ЦАГИ“, т. IV, № 3, 1973.

2. Б е л о ц« р к о в с к и й О. М., Давыдов Ю. М. Нестационарный метод крупных частиц для газодинамических расчетов. „Журн. вычисл. матем. и матем. физ.“, т. И, № 1, 1971.

3. Фонарев А. С. Расчет дифракции ударной волны на профиле с последующим установлением стационарного сверхзвукового и трансзвукового обтекания. ,Изв. АН СССР, МЖГ“, 1971, № 4.

4. Годунов С. К., Забродин А. В., П р о к о п о в Г. П. Разностная схема для двумерных нестационарных задач газовой динамики и расчет обтекания с отошедшей ударной волной. „Журн. вычисл. матем. и матем. физ.", т. 1, № 6, 1961.

5. Nakamura Y. Some contributions on a control surface buzz at high subsonic speeds. Journal of Aircraft, vol. 5, N 2, 1968.

6. Назаренко В. В. Экспериментальное определение переменного давления на элероне при установившихся колебаниях. Труды ЦАГИ, вып. 1541, 1973.

Рукопись поступила 20\11 1974 г.

2—Ученые записки № 3