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

МОДЕЛИРОВАНИЕ ПРОЦЕССОВ |

УДК УДК 533.6(075)

Н. И. Сидняев

О ЗАДАЧЕ КОЛЕБАТЕЛЬНОГО ДВИЖЕНИЯ КОНУСА В СВЕРХЗВУКОВОМ ПОТОКЕ С УЧЕТОМ ПРЕДЫСТОРИИ

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

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

Пусть движение модели летательного аппарата конической формы описывается системой

= 0(0) = 00, 0(0) = 01, (1)

где 0(£) — угловое отклонение от некоторого установившегося состояния; I— момент инерции летательного аппарата относительно оси вращения [2]; М(£) — момент тангажа от аэродинамических сил. Аэродинамические силы, рассматриваемые в невязкой среде, можно определять из уравнений газовой динамики.

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

др д ди

дР + дХ (ри) + эу (рм) = 0,

д д 2 д т(ри) + дХ(р + Ри ) + ду(риу) = 0

д д д - м + - (рН + ду(р + ) = 0,

д Е д д

■Ш + дХ((Е + Р)и) + ду((Е + р*" = °,

(2)

где и, V — проекции вектора скорости V на оси декартовой системы координат; р, р, Т — соответственно давление, плотность и температура; Я — газовая постоянная. Здесь все параметры — безразмерные величины, отнесенные к соответствующим параметрам набегающего

Р 7 — 1

потока. Полная энергия Е = р(Н/ч + V2/2), где - =-к, р = рЯТ,

Р 7

7 = 1,4.

Введем систему координат, связанную с телом (г, ф,г), где ось г параллельна оси вращения конуса. Предположим, что уравнение поверхности конуса в момент времени £ = 0 имеет вид ф = ф(г, г). Тогда при £ > 0 положение границы образующей поверхности конуса определяется из уравнения

Е(Х,£) = ф(г,г) + 0(£) — ф = 0 (3)

с соблюдением граничного условия

БЕ

Dt

= 0, B(x, t) = 0; (4)

Б д

здесь Б = д£ + ™

В рассматриваемой задаче уравнения необходимо дополнить условиями перехода через скачок уплотнения (условия Гюгонио). Пусть для головного скачка уплотнения справедливо соотношение V(X, £) = 0, тогда для условий перехода на ударной волне должны выполняться условия сохранения массы, импульса и энергии:

= К2;

Р1^П1 — N ) = р2(К,2 — N);

Р1 + Р1(^1 — N )2 = Р2 + Р2(К,2 — N )2; (5)

^г Р1 + — N )2 = Л- Р2 + (К,2 — N )2, 7 — 1 Р1 7 — 1 Р2

где индекс 1 относится к параметрам до ударной волны, а 2 — к параметрам за ударной волной; Vn — нормальная составляющая скорости к поверхности ударной волны; Vт — тангенциальная составляющая скорости к поверхности ударной волны; N — скорость движения ударной волны по нормали к ее поверхности [4]. Введем следующие обозначения:

V,! = V,! — N, ^2 = К.2 — N, (6)

7 Р1 7 (пл

к1 =--г, (7)

Р1 7 — 1

с учетом которых запишем условия

' VT1 = VT2;

PVrn = PVn2; < _ _ (8) pi + PiVn2 = P2 + Р2^422;

, hi + К2 = h2 + К2,

где Vni и Vn2 — проекции вектора скорости V на нормаль к головной ударной волне.

Для решения задачи об обтекании острого конуса сверхзвуковым потоком был использован стационарный конечно-разностный метод. Расчетным путем определяли область течения между ударной волной и телом при нулевом положении угла атаки [4], а также рассчитывали параметры потока вдоль образующей конуса.

Момент тангажа в системе (1) связан с давлением р на поверхности острого конуса:

M(t)= УУ (x - xc) X p ds, (9)

B(x,t)=0

где Xc — вектор положения оси вращения в плоскости z = 0.

Из полученных уравнений (5)-(9) следует, что момент М(t) зависит не только от текущего состояния 0(t) и 0(t), но должен определяться предысторией движения #(£) (0 < £ < t) [2]. Очевидно, что в данной постановке момент M(t) является функционалом от 0: M(t) = M [#(£)]. Введем параметры установившегося потока вдоль конуса в момент времени t = 0: u0 — скорость потока, в — угол скачка уплотнения, М0 — число Маха, p0 — давление, р0 — плотность, а0 — скорость звука.

Геометрические параметры конуса приведены на рис. 1: а — угол полураствора конуса, АО = А'О = h — расстояние до оси вращения, АС = l — длина конуса.

Для установившегося потока до появления возмущения (t = 0) параметры p0, М0, р0, U0, а0 определяются из выражений:

(l + ^M^ tg3 в - (M2 - 1) ctg а tg2 в+

+ (1 + ^M2) tg в + ctg а = 0, (10)

2

Po =

Y (y + 1) M2

Y Ml sin2 ß -

Y - 1

(11)

Po =

(Y +1) M^ sin2 в (y - 1) MI sin2 в + 2'

Y - 1

здесь ф = в — a. Вводя 5 = ——г +

2 YPo

uo = aoMo, a2 =-, (12)

1

Po

-, принимая M^ ^ то,

Y + 1 ' M2 sin2 a' Y ^ 1 и используя разложения в степенной ряд, получим ф = 5фь

—— = 5^1, y = 1 + 5yi, после чего из уравнений (10)-(12) можно

M2

оо

получить упрощенные выражения:

Фг =

, Y1 . 2

+ —- sin a —:—2-, Mo2 =

sin a cos a

8фгtga

+ 0(1), po = sin2 a + O(S),

tg a

(13)

Po = ti—+ 0(1), uo = cos a + 0(¿). 0Ф1

Все параметры в дальнейшем приводятся к безразмерному виду. Здесь 1/и0 — параметр времени. Следует отметить, что форма ударной волны [3, 4] в данном случае считается присоединенной к передней кромке и задается зависимостью у = х tg ф.

Для дальнейших рассуждений необходимо получить уравнение поверхности конуса в зависимости от времени. Для решения задачи в такой постановке упростим рассматриваемую модель. Рассмотрим только половину конуса (см. рис. 1) и найдем уравнение прямой Л'В'. Для

Рис. 1. Геометрические параметры полуконуса (О — центр вращения)

1

этого будем использовать следующий метод. В общем виде уравнение прямой записывается как y = kx + b . Величина k равна тангенсу угла наклона прямой к оси x. В нашем случае k = tg(9). Следовательно, для нахождения уравнения необходимо получить значение величины b, для этого найдем координаты точки пересечения этой прямой с осью y. Обозначим эту точку как M, тогда А'О = h, а sin(180° — 9 — а) = sin(9 + а), и из этого уравнения, используя теорему синусов, можно определить

h sin а , h sin а KO = -—-—, тогда AK = h — KO = h —

sin(a + в)' sin(a + в)'

Из треугольника АКМ в соответствии с теоремой синусов найдем

h sin(a + в) h sin(a)

AM =

cos 9 cos 9

h sin a cos 9 + h sin 9 cos a — h sin(a)

cos 9

h sin a(cos 9 — 1) cos 9

+ h tg 9 cos a.

/ h sin a(cos в — 1) \ Найдем координаты точки М: М 0;-----h tg в cos a .

cos в

Подставим координаты точки М в уравнение прямой y = tg(0)x + b и выразим b:

h sin a(cos в — 1)

-----h tg в cos а = b.

cos в

Уравнение прямой будет иметь вид

,.w , . h sin а(1 — cos в)

y = tg^)^ — h cos а) +---—---. (14)

cos в

Следовательно, под действием возмущения граница поверхности конуса находится в движении

B (ж, y, t) = (ж — h cos a) tg в(^) — y + h sin а-^) = 0. (15)

Точка с координатами (h cos a; h sin a) — это неподвижная точка на оси конуса.

Условие (4) на границе имеет вид

(9(f)

u tg 9(í) — v +--0 , [ж — h cos a + h sin a sin в(t)] = 0, (16)

cos2 в^)

B (x,y,t) = 0,

(17)

где и, V — компоненты скорости вдоль осей х и у.

Рассмотрим случай, когда уравнение головного скачка уплотнения имеет вид

W(x,y,t) = y — S(x,t) = 0. Условия Гюгонио (5)-(8) в этом случае представим как [3]

2 Sx

(18)

u =

+

(y + 1) M2 Q

u

+

(1 + тгг S-) + T+T(«»- S<)

V = V^ +

2

p =

P =

Y + 1 2Q2

1 + S2

Q 1

1 + S2 M^j

Y - 1

(y + 1)(1 + S2) Y (Y +1) M2' Q2

Y - 1

Q2 +

2

(1 + S2)

(19)

(20) (21) (22)

Y + 1" ' (y + 1) M2 где Q = St + MooSx - Voo.

Положение (xA, yA) вершины конуса определяются из уравнений:

xA = h [cos a — cos(a — #(t))j, = —h [sin a — sin(a — #(t))j .

(23)

Условие присоединения скачка к вершине конуса следует непо-

средственно из выражения

S (x, t) = —h [sin a — sin(a — #(t))j

(24)

при x = h [cos a — cos(a — #(t))] .

Чтобы найти неустановившееся движение [3] передней кромки конуса, необходимо решить систему уравнений (1), (9), (19)-(21), (24).

Момент тангажа M(t) в возмущенном движении с амплитудой е может быть представлен в виде [4]

M (t) =

M

P OÄ12

= —2epo / (x — h) pi| 0dx.

(25)

i

Выражение (25) можно развернуть [3]:

M(t) = - (1 - 2h) 0(t) - - (1 - 3h + h2) 0(t)-

oo

P0M0Ö0 3

л;

tn

tn ^ t-n

+ 8^ЛП í - T- 4 J] (2 - tn) X' ¡ T0(í - T)dT, (26)

л n

- 4 E ~ [h20(t) - (1 - h) (1 - h - tn) 0(t - tn)] +

n=1

t

00

tn / \ / Z—/ t3

n=1 n 0 n=1 n 0

1_н

где Ьп = 1 — тп, т = --— , Н = М0 tg ф. В случае малых возмуще-

1 + Н

ний [4, 5] соотношения можно представить как [2]

H =

2 + (y - 1) K2 27 • K2 - (7 - 1)

_ 2 [1 + (1 + 2H) K2] 0 = (7 +1) K2 '

Л (2H - 1) K2 - 1

Л = (2H +1) K2 — 1' ^ = V = 0'

где K = М^в.

Полученное уточнение поставленной задачи состоит в том, что теперь в системе (1)

1 = / ,2 (27)

p~( — ) l4 V щ )

и момент тангажа М (£)в системе (1) учитывает предысторию движения. В этом случае согласование начальных условий при t = 0 и предыстории достигается через 0(t) = f (t), —1 < t < 0. В особом случае, когда плотность pw = const, из уравнения (27) следует, что

2pW _____/ 1 , 7 2 , < 2

I = — sin a cos а - - h + h2 + tg2 а , (28)

Р^ \3 у

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

Полученная модель, в принципе, позволяет уточнить анализ переходного движения вторым, третьим и другими приближениями, не ограничивая угол полураствора конуса а. Давление на поверхности конуса, находящегося под углом атаки а < в, в ньютоновском потоке во втором приближении представляется выражением

p(£, t) = po (1 + epi + £2p2) + O(e3), (29)

где po = sin2 ф [2].

2

Причем нормированная координата принимает значения £ = х и

£ = mx:

Pi(£,t) = j sin2^

0(t) + (2£ — hi) 0(t) + ±£ (£ — hi) 0(t)

(30)

P2(£, t) = 02(t) cos 2ф + 0(t)0(t) (4£ + 2himi) + 0(t)0(t) (£ — hi) x

+

x (m2£ — hims) + 02(t) cos2 ф [(£ — hi)2 + £hi tg a tgф] + 20(t) + (£ — hi) 0(t)l [£ tg2 ф + him4] cos2 ф0^ — £ sec ф) —

5 sec ф

20(t) sin2 ф cos ф + 0(t) (1 + cos ф) cos2 ф (£ — hi) ] / 0(t — т)dr,

(31)

mi =

cos ф (sin a sin ф — 2 cos a cos ф)

cos a

m2 = 2 + cos ф,

cos ф cos(a + ф)

тз =-,

cos a

m4 = 1 — tg a tg ф.

1 sec <r

Mi(t) = — E / (£ — hi)Pi(£,t)d£. (32)

j=-! о

После подстановки уравнений (28)-(32) в систему (1) получим

2 (J + Ja) 0(t) + L0(t) + So0(t) = 2eM2 [0(t)],

(33)

где

4tg a A 2 ,

L =-cos a--3h2 + 2h2 ,

cos a V3

S0 = 4 tg a cos 2a (1 — 2h2)

Дополнительный момент инерции задается зависимостью

tg(a + a)+tg(a — a) ^ L2

Ja = --- I 4 — 3 h2 + 2 h2 ) ,

cos2 a

Ja > 0,

последний всегда положителен, и оценить влияние предыстории становится возможным через отношение /а//. Если пренебречь членами

второго порядка, то из уравнения (33) видно, что влияние предыстории проявляется только через присутствие момента инерции /а, что

приводит к уменьшению скорости демпфирования с величины — до Ь

величины —-- и к небольшому изменению частоты колебаний.

2 (I + 4) У

Рассмотрим уравнение движения модели при вращательных колебаниях [1]:

_d а da __, ч _ . , ч _

+ ^~dt + K(а + ао) - M(t) = 0;

(34)

здесь Мг (£) — аэродинамический момент, действующий на модель; К — приведенная жесткость измерительной системы. Выражая Мг (£) через коэффициенты аэродинамических производных, из уравнения (34) получим

i - 2 <Z P~si) ^ +

р - [та + ш^) 2 PooMol S

1

da dt+

+ (K - m?P0U2Sl^a + Kao - шгоPooU0Sl = 0. (35)

Видно, что при вращательных колебаниях т^г является присоединенным моментом инерции, т^ + т^ г — коэффициентом демпфирования и, наконец, т^ соответствует аэродинамической жесткости, Б — миделево сечение. Как известно [1] , общим решением уравнения (35) будет

а = а0 + еоэ^^ — 0). (36)

Здесь а0 — координата среднего положения и средний угол атаки, относительно которого совершаются колебания, а второе слагаемое выражает закон свободных колебаний. Кроме того, р1 — круговая частота свободных колебаний системы, £ — декремент их затухания, а1 и 0 — постоянные интегрирования, зависящие от начальных условий движения. Эти величины следующим образом выражаются через коэффициенты дифференциального уравнения движения (35):

= р - Popupsi'^mq + щдZ) /| 4 I - mZZPoSl3/2

2 _ K - maPo u2Sl/2

Pi =

I - mZZPoSl3/2

— а

(37)

где частота колебаний р1 и величина декремента затуханий £ определяются зависимостями

2п . 1 а Pi = , с = — In —; 1 Т' ц Т а1

Т — период колебаний.

По найденным величинам р и £, учитывая выражение (37), получаем аэродинамические производные момента:

Ка0

тго =

та -

р^ u2Sl/2'

p

о

2п

K - I(p2 + с2)

р^« Sl/2

та + =

р - 2С(1 - p^l3S/2)

р^ «oS/2/2

р - 2С1

р^ MoS/2/2'y

(39)

При этом величиной (1/2)ш^гр^£/3 при испытаниях в воздушном потоке можно пренебречь по сравнению с приведенным моментом инерции системы I.

Для определения приведенной массы т (или приведенного массового момента инерции I) проводится следующий эксперимент. Так, например, в испытаниях при и0 = 0, т.е. без потока, определяется частота р0 свободных колебаний системы [1, 6, 7]. Затем в систему вводится дополнительный груз, масса которого Дт или момент инерции относительно оси колебаний Д!гг известны, и определяется частота колебаний р системы с дополнительным грузом [5-7]. Тогда приведенные масса и момент инерции находятся из соотношений

m = Am

p2

p0 - p2:

I = AI

p2

ZZ 2 о '

p0 - p2

(40)

Коэффициент момента тангажа и его производные для острого конуса с образующим вк определяются из соотношений:

2п l

mZ = 1 J J p(x,t,ae)/0 sin ^d^dx, 00

та =

dmz да

l 2п

1

а=0,

~Q] I I sec2 ßk - Х0) x tg ßk sin ^d^dx,

(41)

Sl J J да

00

ае—а

mt = m? + =

2u

l

2uœ

SÏ2"

I 2n

dp

x sec ßk — xq) x tg ßk sin ^d^dx,

0 0

где Iq = x sec2 ßk — xo — плечо сил; I — длина острого конуса, колеблющегося около угла ае = Ф — (x sec ßk + xo)/uœ; $ — угол тангажа; xo — продольная координата оси колебаний тела. На рис. 2 показана зависимость коэффициентов момента и давления от колебания в плоскости угла атаки а.

Для определения точности вычислительной модели проводили оценку полученных данных по экспериментальным результатам при

Mœ = 6 по критерию Струхаля Sh = — = 3 • 10-3 и коэффициенту

демпфирующего момента + = 0,12 [1, 4, 5].

Таким образом, следует отметить, что толщина ударного слоя зависит не только от угла атаки и скорости его уменьшения, но и от знака производной угла атаки по времени da/dt. Эти обстоятельства предыстории, означающие запаздывание по времени формирования отмеченных особенностей течения в процессе колебательного движения, существенно влияют на нестационарные аэродинамические характеристики острого конуса. При несимметричном обтекании острого конуса происходит уменьшение разности давлений между наветренной и подветренной его сторонами, что приводит к снижению аэродинамического коэффициента нормальной силы, а это, в свою очередь, приводит к увеличению амплитуды колебаний. Степень влияния угла атаки на величину динамической производной устойчивости зависит от положения оси колебаний. Сдвиг оси колебаний вперед вызывает увеличение динамической устойчивости затупленного тела вращения. При более близком к корме положении оси колебаний с увеличением угла атаки острый конус становится динамически

а

0,44

0,28 - пг.

" о- 0,2- 0,4' - 0,6" Л, р

» Ï а

0 Ï

1 Л

11) L /^v

ml

0,12 -пл.-

■0.04-

0

Рис. 2. Зависимости коэффициентов момента тангажа, давления и угла атаки от времени при затухающих колебаниях конуса

неустойчивым. Экспериментально подтверждаются также предсказанные расчетом тенденции изменения динамической устойчивости при других хь//.

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

СПИСОК ЛИТЕРАТУРЫ

1. Сидняев Н. И. Методика численного расчета сверхзвукового обтекания колеблющегося осесимметричного тела вращения в условиях интенсивного поверхностного массообмена // Вестник МГТУ. Серия "Естественные науки". -2003. - № 1 (10). - С. 71-87.

2. Transient motion of a hypersonic wedre including time history effects: Hiu W.H., Van Roessei H.J., "J. Guid., Contr. And Dyn.", 1986, 9, No. 2, 205-212.

3. Теленин Г. Ф., Липницкий Ю. М. Нестационарное сверхзвуковое обтекание затупленных тел с отошедшей ударной волной // МЖГ. - 1966. - № 4.

4. Сидняев Н. И. Метод расчета нестационарного обтекания тела вращения с поверхностным массообменом в рамках параболизированных уравнений Навье-Стокса // Математическое моделирование. - 2004. - Т. 16, № 5. - С. 55-65.

5. Сидняев Н. И. О методике исследования коэффициентов вращательных производных аэродинамического момента конуса с поверхностным массообменом // Изв. вузов. Авиационная техника. -2004. - № 2. - С. 30-33.

6. Белоцерковский С. М. О коэффициентах вращательных производных. Аэродинамика неустановившихся движений. - М.: Оборонгиз, 1958. - 106 с.

7. К р а м е р В. В. Экспериментальное определение коэффициентов вращательных производных кинематическим способом. Аэродинамика неустановившихся движений // Тр. ЦАГИ. 1958. - Вып. 725. - С. 81-98.

Статья поступила в редакцию 18.02.2005

Николай Иванович Сидняев родился в 1955 г., окончил в 1981г. МВТУ им. Н.Э.Баумана и в 1985г. МГУ им. М.В.Ломоносова. Канд. техн. наук, доцент кафедры "Высшая математика" МГТУ им. Н.Э. Баумана. Автор около 100 научных работ в области прикладной математики и механики.

N.I. Sidnyaev (b. 1955) graduated from the Bauman Moscow Higher Technical School in 1981 and the Moscow State University n.a. M.V. Lomonosov in 1985. Ph. D. (Eng.), assoc. professor of "Higher Mathematics" department of the Bauman Moscow State Technical University. Author of about 100 publications in the field of applied mathematics and mechanics.