О связывании ближнего и дальнего полей в задаче о звуковом ударе

_________УЧЕНЫЕ ЗАПИСКИ ЦАГИ

_______ __

№3-4

УДК 534.83:629.7.015.016.54

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

Ю. Л. Жилин , В. В. Коваленко

Представлены уточненные методы расчета начальных данных для решения задачи о распространении волны звукового удара от произвольной компоновки летательного аппарата.

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

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

Форма эпюры (распределения возмущенного давления в потоке) в ближнем поле зависит от формы тела и режима полета, а процесс эволюции эпюры — от состояния атмосферы. Все это учитывает теория звукового удара.

Согласно работе [1], зависимость избыточного давления Р от времени / при звуковом ударе определяется уравнениями

Р(0 = кіРіч),

і = к2[ц-кР (т\)],

І2-Ц 1

Л2 - Лі к ’

_ _ ЇЇ2_ (^2+^1 )( Л2 ~ Лі ) = 2 | ^ Л-

ЛІ і

. Функция ^(л) пропорциональна избыточному давлению и задает распределение возмущений около тела. Она должна быть получена из решения задачи о ближнем поле. Функция РЩ — это искомая зависимость давления от времени в волне звукового удара на земле.

В уравнениях (1) переменная п ~ расстояние вдоль оси тела, отнесенное к характерной длине; ? — время (рис. 1).

Уравнения (1) учитывают в первом приближении нелинейные процессы, происходящие при распространении слабых ударных волн в неоднородной атмосфере.

Условия их применимости: возмущения малы по отношению к невозмущенному потоку;

радиус кривизны фронта волны намного больше характерного размера пакета возмущений (длины тела); .

характерный масштаб атмосферных неоднородностей намного больше длины тела.

Первые два уравнения системы (1) описывают затухание возмущений на пути от летательного аппарата до Земли, а следующие два — процесс образования и перемещения ударных волн в пакете возмущений. Величины с индексами (7^, /2, гц, т^) указывают тот отрезок эпюры, который «схлопнется» в ударную волну. Коэффициенты к, к], &2 зависят только от состояния атмосферы и режима полета. Формулы для вычисления этих параметров представлены в работе [1].

Система уравнений (1) может быть решена численно. Алгоритм решения изложен в работе [2].

2. Для решения системы (1) необходимо знать распределение Р(ц). Как указывалось, оно должно быть получено из решения задачи об обтекании тела. Согласно работе [3], для осесимметричных тел указанное распределение вычисляется по так называемой функции Уитхема Ф(х):

Эона сращивания

Дальнее поле Теория звукового удара

Профиль Волны давления^, на поверхности Земли

'"шг-

Рис. 1. Схема решения задачи о звуковом ударе

где 5(х) — распределение по: длине площади проекции на плоскость х= 0 сечений тела плоскостью х = х, - гр (р = а/м2 - 1).

Из асимптотики линейной теории [4]—[6] следует, что на больших расстояниях от произвольного трехмерного тела затухание возмущений в каждой меридиональной плоскости, проходящей через ось тела, происходит так же, как от некоторого эквивалентного тела вращения. Это позволяет свести задачу расчета пространственных возмущений в дальнем поле к двумерной (или же одномерной нестационарной) и воспользоваться формулами (2). Для этого, однако, необходимо определить краевые (начальные) условия для каждой из упомянутых меридиональных плоскостей. В этом случае

/Чл, S) = АфОь Э), Ф

где S’(x, 9) = ~^LS(x, 9), а площадь эквивалентного тела вращения S(x, 9) вычисляется с учетом подъемной и боковой сил:

^

2 q J

где Sb(x, 9) — площадь проекции на плоскость х = 0 сечения поверхности тела плоскостью х = Xj - p(^i sin 9 + у\ cos 9); Y(x, 9) и Z(x, 9) — подъемная и боковая силы, действующие на тело до сечения его указанной плоскостью; q — скоростной напор.

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

3. Применение формул. (3), (4) в их точном виде не является простым делом. Во-первых, следует решить геометрическую задачу об определении сечений аппарата наклонной плоскостью и вычислить распределение площадей проекций этих сечений на плоскость х=0. Далее следует вычислить подъемную силу, действующую на ту часть аппарата, которая лежит впереди (от носка тела до линии пересечения поверхности аппарата с наклонной плоскостью). После этого полученное суммарное распределение надо продифференцировать, подвергнуть интегральному преобразованию (для вычисления Ф(х)) и, наконец, снова продифференцировать для получения искомой функции F{t\).

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

Рассмотрим схему обтекания тела (рис. 2).

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

у = -РосМоо JJ (vnx - ипу )dS, z = -Р«оо JJ (wnx - unz )dS,

s

I = ~ JJ (~P2««x + vny + wnz

где Y и Z — составляющие вдоль осей у и z силы, действующей на находящуюся внутри S часть тела; и, v, w — компоненты возмущенной скорости вдоль осей х, у и z соответственно; пх, пу, nz — составляющие

единичного вектора внешней нормали к S; Б — разность площадей проекции на плоскость х= const концевого и входного сечений тела поверхностью S. Здесь используется прямоугольная система координат, ось х которой направлена вдоль вектора скорости набегающего потока.

^ Выберем контрольную по-

верхность следующим образом. Введем поверхность фронта Jj, отделяющую возмущенный поток от невозмущенного. Проведем плоскость 62 > параллельную оси х и не пересекающую тело (см. рис. 2). На этой плоскости выделим отрезок АВ, где точки А и В лежат на линии пересечения ПЛОСКОСТИ 1S2 и поверхности Si- Проведем характеристические поверхности S3 и 64 (огибающие обратных конусов Маха), содержащие АВ.

Соотношение (5) применим к двум контрольным поверхностям SрЬ и Snb. Эти поверхности содержат часть Si и плоскость S2', при этом Spb содержит внутри тело и замыкается плоскостью S$, a Snb — поверхностью 1S4 и внутри нее нет тела.

В работе [7] было показано, что интегралы (5) по поверхности равны нулю, если на тело набегает невозмущенный поток. Поэтому интегрирование в (5) производится по плоскостям £2, £3 и 1S4, а величина £ равна площади проекции на плоскость х=const сечения тела плоскостью Sj.

Рис. 2. Контрольные поверхности при вычислении интегральных соотношений в сверхзвуковом потоке

На плоскости 62 пх =0, пу = - сое 9, пг = вшв, где 0 — угол между отрицательным направлением оси х и внешней нормалью к 52.

х= const вектора внешней нормали к 53. В дальнейшем будем предполагать, что углы 0 и & фиксированы, а плоскости S2 и £3 перемещаются в пространстве, оставаясь параллельными соответствующим исходным плоскостям. При этом новое положение плоскостей характеризуется величиной лсо — координатой точки пересечения Л3 с осью х, а левые части соотношения (5) являйэтся функциями этой точки. С учетом всего этого (5) можно записать в виде

Упрощение полученных соотношений связано с исключением интегралов по плоскости 63. Для этого первое из соотношений (6) умножается

Соотношение (7) уточняет и обобщает определение эквивалентного тела вращения, принятое в теории звукового удара. Оно позволяет вычислить ^(Хо), если известно поле возмущенных скоростей при обтекании трехмерного тела. Важная особенность соотношения (7) заключается в том, что при фиксированной плоскости 5з можно менять расстояние ПЛОСКОСТИ 52 от оси х. Поэтому функция ^(Хо) может быть вычислена по полю скоростей на любых, в том числе и малых, расстояниях от тела. Применим соотношение (5) к поверхности Так как внутри поверхности нет тела, то левые части (5) равны нулю, а на плоскостях 62

На плоскости 5з пх = пу = -^-совЭ, пг = -^ш 3, где 9 — угол между положительным направлением оси у и проекцией на плоскость

ш: Пу м

Y(х„) = —Рооих cos01| udS2 - (v ~ u$cos$)dS3,

Z(x0) = pxua sin e||udS2 - JJ(w * upsin&)dS3,

на -—j-, а второе — на

Poo^oo

В результате получаем

Pcos9

и все три соотношения складываются.

seq(x 0) = £(Х0) + —&y[y(x0)cos9 + Z (х0) sin 9 ] =

Роо ^00

("Pcos(9 + 9)w - vcosQ + M>sin0)<£y2.

и 5} соответственно имеем пх = 0, Пу = cos0, пг - -sin0 (на ф) и пх = ,

= icos91; = ^sin 9, (на 54), где 9, + 9 = п + 20 (9! — угол между

осью у и проекцией на плоскость х=0 вектора внешней нормали к 64). При этом

0 = собоЦый?^ -•j^'JJ(u - «pcosd,)*/^,

S2 ^4

О = - sin 0 JJ udS2 - JJ

s2 ^4

0 = JJ (,cos0-wsin0)^2 + JJ (-рм + v cos 9, + w sin 9])^54.

s2 s4

Умножая первое из этих соотношений на pcos9[, второе на р sin 9] и комбинируя с (7), получаем

Seq(x0) = cos20cos(0 - 9|) \[udS2 =

2

= —JJ(-psin20sin(0 - - ucos0 + wsin d)dS2.

°° ^2

Соотношения упрощаются, если плоскость 1S2 и отрезок А В перпендикулярны ОСИ у. При ЭТОМ 0 = 9, $1 =71 и -

Seg(x0) = Е(х0) + SIM = - j_ ff (Р„ + v)dS2 =

РсоК u<» V

„

= ffiftJS, = —— \\vdS2.

«о, {J ‘ J J z

*>2 s2

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

. Выполнив дифференцирование выражения для ^„(xq) и про-

dxQ

ведя необходимые преобразования, можно получить

- _L f

dx0 Ux J4 и v ' x^x'(z') ’ 4 7

z'a

а при 0 = 9, 9j =7t

z'b z’b - z'b

^i-± h+u *. .a j. 1 л*=_jl j„ u *,

. za. za

где интегралы вычисляются вдоль отрезка АВ\ х' = х'(г') — уравнение отрезка АВ.

Из приведенных соотношений следует, что функция ^(х0) и ее производная могут быть вычислены несколькими способами, если известно поле возмущенных скоростей на любом расстоянии от тела. Это обстоятельство может оказаться полезным для проверки теоретических результатов и выбора подходящего алгоритма расчета. Для таких вычислений можно использовать любые программы расчета обтекания трехмерных тел сверхзвуковым потоком газа.

На рис. 3 представлен пример расчета звукового удара от компоновки корпус — крыло (диаметр корпуса 4 м, длина 80 м, масса 300 т, высота полета 16 км, число Маха полета М = 1,5). Поле течения около тела вычислялось при помощи интегрирования полной системы уравнений Эйлера [9]. Интегралы (12) вычислялись для трех расстояний г плоскости ^2 от оси тела (г=0,1; 0,2 и 0,3 полной длины тела /_,). Видно, что результат не зависит от г (небольшие отличия в зависимостях йБ/бх и Дл) являются следствием слабой нелинейности, которая улавливается полной системой Эйлера).

На рис. 4 дано сопоставление расчетных результатов для компоновок с передним горизонтальным оперением (ПГО) и без него. Как видно, использование (12) позволяет легко и естественно отслеживать влияние изменений в форме аппарата.

4. Вычисление функции Ял) по 5ея(х) согласно формулам (4) или (7) есть следствие линейной теории обтекания тел сверхзвуковым потоком. Соответственно, в тех случаях, когда применять линейную теорию для расчета поля течения около летательного аппарата нельзя

&$ед{х)

АХ гГ°’Г

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

*(*)[ 700 Ъ i

Рис. 4. Сопоставление расчетных данных для компоновок с ПГО и без

него

(например, при больших сверхзвуковых скоростях или же когда нельзя пренебрегать интерференцией толщины и подъемной силы), указанный выше способ связывания ближнего и дальнего полей неприменим.

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

Так как из выражения для коэффициента кх следует [1], что кл

обратно пропорционален 4г, то непосредственно из

£><*>-«>.

У V р V L рх v '

В выражении (13) L — длина тела, У — подъемная сила, х — продольная координата. Независимость F от г для трехмерных течений является критерием выхода на квазиосесимметричное асимптотическое решение. Таким образом, расчет ближнего поля необходимо вести до тех пор, по-

ка не установится связь (13), далее решение может быть продолжено как решение системы (1). Надо отметить, что соответствие между переменными X И Т1 меняется при изменении г со вторым уравнением системы (1), т. е. происходит деформация эпюры Р(х) в продольном направлении. Сопоставление кривых при разных г следует проводить с учетом пересчета продольной координаты.

Расчет обтекания тела с использованием полной системы уравнений Эйлера — это процедура стандартная и не вызывающая затруднений; в данной работе использованы методика и программа, описанные в [9].

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

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

В программе [9] уравнения Эйлера интегрируются маршевым конечно-разностным методом в криволинейных координатах в дифференциальной форме. Известно, что решение при этом подвержено влиянию «криволинейности», и надо принимать специальные меры для его устранения [10], [11]. В настоящей работе это выполняется следующим образом.

При переходе от исходных декартовых координат (х, у, I) к криволинейным е=х , /=/(*, у, г), .5 (х, у, г) система уравнений Эйлера мо-

жет быть представлена в следующем виде:

Ее + Ft + Gs — 0, те Е = gE; F = g(txE + tyF + tzG); G = g(sxE + syF + szG);

' р и ' pv £ CL

Е = ри2 + кр , F = puv , G = piw

р UV pv2 + кр pwv

<Puw , vp VW J vpw2 + kp,

k = 2l-, s = det

1 - 1 8

В уравнениях (14) и, v, w — компоненты полной скорости (а не возмущений, как в (5)—(12)); они обезразмерены отнесением к скорости адиабатического истечения в вакуум, а р и р — к параметрам торможения невозмущенного потока. Если поток не возмущен, т. е. и = им, V = W = Woo, Р = Рх> Р = Роо» то Уравнения (14) выполняются тождественно, так как сводятся к соотношениям типа zst - Zts и тому подобным. По-иному дело обстоит с конечно-разностными уравнениями , полученными при разностной аппроксимации уравнений (14) с использова-

нием значений переменных и метрических коэффициентов в узлах сетки. В этом случае ошибки аппроксимации приводят к тому, что Е-const не является точным решением разностных уравнений, т. е. невозмущенный поток подвергается искажению. Этой неприятности можно избежать, если рассматривать дифференциальные уравнения (14) в терминах вариаций относительно Ет, т. е. записать Е = Е00 + Е, F^F^ + F,

G = Gx + G, и с учетом (Ех)е + (F^t + (G^g = 0 вместо (15) конечно-

разностному решению подвергнуть уравнение

Ёе + F, + Gs =0.

Так как уравнения (14) и (15) однородны относительно Е, F и G, то конечно-разностные уравнения обращаются в тождество при Е = F = G = 0

и, следовательно, невозмущенный поток является точным решением разностных уравнений.

Система (14) или (15) замыкается уравнением Бернулли, которое с

2 2 2

учетом принятого обезразмеривания имеет вид р = р(\ - и -v - w ). На поверхности тела ставится условие непротекания '

(Vn) = ufx + vfy + wfz = 0,

где f(x,y,z) = 0 — уравнение поверхности тела, а на головном скачке уплотнения — условия Гюгонио: •

PlKl = P2Vn2’ Pl^nl^nl + kpi = p 2Vn2Vn2 + кРЪ Kl = K2-

В уравнениях (16) величины с индексом 1 относятся к параметрам течения до скачка, а с индексом 2 — за ним; Vn — нормальный к поверхности скачка компонент скорости; Vx — вектор скорости, касательной к скачку.

Одной из неприятных особенностей задачи является то, что при перемещении маршевой плоскости вниз по потоку растет поперечный размер расчетной области (размер отрезка Со — F на рис. 5). Это ведет к укрупнению ячеек разностной сетки, что, в свою очередь, требует дополнительных усилий для поддержания высокой точности в пространственном разрешении особенностей течения. Следует заметить, что для расчета звукового удара вычисления решения во всей зоне от Со до F не требуется. Как видно из характерного вида эпюры в ближнем поле (см. рис. 1), в поле течения содержатся головной и хвостовой скачки уплотнения с некоей зависимостью ^(л) между ними (зона а — Ь), а далее (за хвостовым скачком) следует зона b — с, где происходит окончательное затухание возмущений при с -> оо. Уровень возмущений в зоне b — с намного меньше уровня в основной части эпюры, и поэтому эту зону (или большую ее часть) Можно исключить из рассмотрения и расчета. При конечно-разностном решении уравнений Эйлера это значит, что можно исключить из расчета часть зоны Cq — D (рис. 5). Этого можно достичь введением, внутренней границы расчетной области А — Q — В, которая в настоящей работе выстраивалась как характеристическая поверхность,

Рис. 5. Расчетная область при конечно-разностном решении уравнений Эйлера

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

В качестве примера на рис. 6 представлены результаты расчета звукового удара, порождаемого аппаратом, перемещающимся со скоростью М = 2,5 на высоте 18000 м. Как видно из графиков зависимостей Дл), решение при г> 4Ь следует считать установившимся.

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

поле

5. Таким образом, начальные данные для расчета распространения волны звукового удара в атмосфере (данные ближнего поля) могут быть вычислены различными способами. Для слабовозмущенных течений, когда применима линейная теория, эффективным является использование интегральных соотношений (7) и (12). Для нелинейных течений

п

<

і

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

ЛИТЕРАТУРА

1. Жилин Ю. Л. О звуковом ударе//Ученые записки ЦАГИ.— 1971.

Т. II, № 3.

2. Ж и л и н Ю. Л., Ч е р н ы ш е в С. JT. Алгоритм построения эпюры избыточного давления при звуковом ударе/Друды ЦАГИ.— 1981. Вып. 2110.

3. W h i t h a m G. B. The behavior of supersonic flow past a body of revolution, far from the axis//Proc. Roy. Soc., ser. A.— 1950. Vol. 201, N 1064.

4. W a 1 k d e n F. The shock pattern of a wing-body combination far from the flight path//Aeron. Quart.— 1958. Vol. IX, N 2.

5. Hayes W. D. Linearized supersonic theory//North Amer. Am. Co.,

Re AL-222, Los Angeles.— 1947.

6. G e о rg e A. K. Reduction of sonic boom by azimuthal redistribution of overpressure//AIAA.—1969. Vol. VIII, N 3.

7. Жилин Ю. JI. Крылья минимального сопротивления//ПММ.—

1957. Т. 21, вып. 2.

8. Жилин Ю. Л. Полная подъемная сила, действующая на тела, слабовозмущающие сверхзвуковой поток//Ученые записки ЦАГИ. — 1971.

Т. II, № 4.

9. Коваленко В. В., Минайлос А. Н. Расчет невязкого сверхзвукового течения около комбинации крыло — фюзеляж//Труды ЦАГИ.—

1984. Вып. 2251.

10. Thomas P. D., Lombard С. К. The geometric conservation law a link between finite-difference and finite-volume methods of flow computation on moving grids//ALAA Paper 78-120,— 1978.

11. Hindman R. G. Geometrically induced errors and their relationship to the governing equations and the treatment of generalized mapping//AIAA Paper 81-1008.- 1981.

Рукопись поступила 8/II1997 г.