1. На главную
  2. Электронные версии
  3. Журнал вычислительной математики и математической физики
  4. T. 62, № 4, 2022

Журнал вычислительной математики и математической физики, 2022, T. 62, № 4, стр. 642-658

Об упрощении численных и аналитических “инструментов” описания “звукового удара”

Х. Ф. Валиев 1, А. Н. Крайко 1, *, Н. И. Тилляева 1

1 ЦИАМ
111116 Москва, ул. Авиамоторная, 2, Россия

* E-mail: akraiko@ciam.ru

Поступила в редакцию 27.06.2021
После доработки 28.07.2021
Принята к публикации 17.11.2021

Полный текст (PDF)

Аннотация

В свете современного состояния и тенденций развития методов математического моделирования звукового удара от сверхзвуковых летательных аппаратов обсуждается место численных и аналитических инструментов его описания. Отмечена возрастающая роль численных расчетов на адаптированных к особенностям течения сетках в рамках стационарных уравнений Эйлера (в координатах летательных аппаратов на “крейсерском” режиме полета) до расстояний от нескольких до пары десятков его длин. Другая важная тенденция – замена развивавшегося с середины ХХ в. сложного численно-аналитического аппарата описания “среднего” и “дальнего” полей звукового удара более простыми подходами, в том числе, без обращения к функции Уизема. В развитие этих тенденций в рамках уравнений Эйлера продемонстрирована возможность численного расчета типичных для звукового удара волновых структур без ограничений на расстояния и на интенсивности ударных волн, включая крайне малые. Описание эволюции звукового удара с удаления в 15–20 длин летательных аппаратов и до Земли сведено к мгновенному решению следующих из осесимметричных уравнений Эйлера задач Коши для обыкновенных дифференциальных уравнений. Вязкое размазывание слабых ударных волн описывает известное одномерное стационарное решение уравнений Навье–Стокса. Библ. 60. Фиг. 8.

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

ВВЕДЕНИЕ

Достаточно полное представление о развитии и современном состоянии инструментов математического моделирования звукового удара, возникающего при полете сверхзвуковых летательных аппаратов (ЛА), и аналогичных нестационарных ударно-волновых структур дают статьи и монографии [1]–[45 ] и доклад [46]. Важную роль в этих интенсивно развивавшихся с середины ХХ в. инструментах играют приемы и формулы линейной теории и нелинейные поправки к ним.

Если возмущения течения и интенсивность ударных волн (УВ) превышают допускаемые линейной теорией, то необходимы расчеты ближнего и, возможно, среднего полей обтекания ЛА в рамках более полных уравнений. Часто при моделировании звукового удара обтекание ЛА описывают уравнениями течения идеального (невязкого и нетеплопроводного) газа – уравнениями Эйлера. В этом приближении на “крейсерском” режиме полета в системе координат ЛА течение стационарное и в основном сверхзвуковое. Его расчет возможен почти всюду – маршем по координате х, направленной по скорости V0 невозмущенного (“набегающего”) потока. Программы численного решения уравнений Эйлера или RANS-приближения уравнений Навье–Стокса стали еще одним инструментом, правда, как правило, для расчета ближнего поля с размерами от нескольких до 8–10 длин ЛА (см. [23], [25], [28], [33], [38], [41], [44], [46]). Известное авторам исключение – статья [37] с расчетами до r = 20. Здесь и далее линейный масштаб L° – длина ЛА или заменяющего его в разд. 2 источника возмущений (величины с верхним “градусом” размерные), r – отнесенная к L° радиальная переменная цилиндрических координат xrϕ с ϕ = 0 на полуплоскости xy с направленной вниз (к Земле) осью y декартовых координат xyz с началом в передней точке ЛА. В сравнении с расчетами обтекания ЛА в RANS-приближении расчет по уравнениям Эйлера без сгущения сетки у обтекаемых без пограничных слоев поверхностей проще и быстрее.

Возникает вопрос, что препятствует широкому применению чисто численных инструментов описания звукового удара до r $ \gg $ 1. Существует мнение, что это невозможно из-за погрешностей счета малых возмущений параметров набегающего потока на больших расстояниях от ЛА. В действительности, однако, такое снижение точности присуще разностным сеткам, неадаптированным к особенностям рассчитываемых течений. В [47]–[49] это показали расчеты ударно-волновых структур перед сверхзвуковыми решетками. В связи со звуковым ударом достоинства адаптированных сеток и бесперспективность неадаптированных прямоугольных продемонстрированы в [35]. Еще раньше адаптированные сетки при расчете ближнего поля осесимметричных и простых пространственных тел применили авторы [23]. Кстати, сравнения выполненных в [35] расчетов с экспериментами из [12] по обтеканию тел вращения вопреки утверждениям [35] показали, что при определении ударно-волновой структуры звукового удара уравнения Эйлера не уступают уравнениям Навье–Стокса и их RANS-приближению. Это принципиально, особенно в свете данных [46] о несопоставимых временах RANS-расчетов ближних полей ЛА (при r ≤ 1/6) с расчетами маршем на адаптированных сетках в рамках уравнений Эйлера их средних полей (при 1/6 ≤ r ≤ 10). В [37] неструктурированные адаптированные сетки строились при расчете в рамках уравнений Эйлера сверхзвукового обтекания тел вращения (в том числе, из [12]) и пространственных моделей ЛА.

С учетом опыта расчета плоских (см. [47]–[49]), осесимметричных (см. [23], [35]) и пространственных (см. [23], [35], [37], [46]) ударно-волновых структур для правильного определения параметров течения в конически-кольцевой зоне звукового удара между “головной” и “замыкающей” УВ нужно явно или неявно (сгущением сетки, адаптированной к структуре звукового удара) выделять эти УВ. С удалением от ЛА ширина таких кольцевых зон растет много медленнеe чем r.

В приводимых ниже примерах марш в направлении оси х по стационарному сверхзвуковому аналогу схемы Годунова (см. [50]), модернизированной согласно предложениям (см. [51]–[54]), с высокой точностью рассчитывает звуковой удар до x $ \gg $ 1. Кроме того, при r $ \gg $ 1 в силу известной (см., например, [25]) осесимметричной асимптотики изучаемые пространственные течения в фиксированных меридиональных плоскостях, отвечая обтеканию разных тел вращения, удовлетворяют осесимметричным уравнениям Эйлера. Благодаря этому с r = 15–20 их описание при ϕ = const и при однородной, и при стратифицированной атмосфере сведено к практически мгновенному решению задачи Коши для обыкновенных дифференциальных уравнений – следствий осесимметричных уравнений Эйлера для слабо возмущенных коротких волн. Развитый подход аналогичен, но не тождественен подходу, представленному в [15], [16].

1. УРАВНЕНИЯ ОСЕСИММЕТРИЧНОГО СТАЦИОНАРНОГО СВЕРХЗВУКОВОГО ТЕЧЕНИЯ В СТРАТИФИЦИРОВАННОЙ АТМОСФЕРЕ. ПРИБЛИЖЕНИЕ “КОРОТКОЙ” ВОЛНЫ

При удалении от ЛА интенсивность возникших при его сверхзвуковом обтекании УВ, измеряемая отношением давлений p+/p– за и перед ними, почти монотонно приближается к единице. В слабых УВ приращениями энтропии и одного из инвариантов Римана – величинами порядка (p+/p– – 1)3 (см. [55]–[58]) – будем пренебрегать. Это, тем не менее, не означает изэнтропичности рассматриваемых далее течений. В общем случае удельная энтропия невозмущенного потока s0 переменна из-за неоднородности атмосферы по координате y = $r\cos \phi $. Масштаб изменения s0 много больше длины ЛА.

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

(1.1)
$\frac{{d{{p}_{0}}}}{{dy}} = {{\rho }_{0}}g \Leftrightarrow \frac{{\partial {{p}_{0}}}}{{\partial r}} = {{\rho }_{0}}g\cos \phi .$
Здесь и далее ρ – плотность газа, g – ускорение свободного падения, индекс “нуль”, как и выше, отмечает параметры набегающего потока, а при стратифицированной по y = $r\cos \phi $ атмосфере p0 = p0(y) и ρ0 = ρ0(y). В используемых ниже координатах ЛА скорость набегающего потока $V_{0}^{ \circ }$ – константа. Взяв $V_{0}^{ \circ }$ за масштаб скорости, получим u0 = V0 = 1, ${{{v}}_{0}}$ = θ0 = 0, a0 = 1/M0, где u и ${v}$ – x- и r-компоненты вектора скорости V, θ – угол его наклона к оси х, а – скорость звука и М = V/a – число Маха.

Важный результат исследований звукового удара – упомянутая выше осесимметричная асимптотика. Согласно ей, при r $ \gg $ 1 изучаемое пространственное стационарное течение в фиксированных меридиональных плоскостях описывается осесимметричными уравнениями Эйлера (h = h(р, s) – удельная энтальпия)

(1.2)
$\begin{gathered} \frac{{dp}}{{dt}} + \rho {{a}^{2}}\left( {\frac{{\partial u}}{{\partial x}} + \frac{{\partial {v}}}{{\partial r}} + \frac{{v}}{r}} \right) = 0,\quad \frac{{du}}{{dt}} + \frac{1}{\rho }\frac{{\partial p}}{{\partial x}} = 0,\quad \frac{{d{v}}}{{dt}} + \frac{1}{\rho }\frac{{\partial p}}{{\partial r}} = g\cos \phi , \\ \frac{{ds}}{{dt}} = 0 \Leftrightarrow \frac{{dh}}{{dt}} = \frac{1}{\rho }\frac{{dp}}{{dt}},\quad \frac{d}{{dt}} = u\frac{\partial }{{\partial x}} + {v}\frac{\partial }{{\partial r}} \\ \end{gathered} $
с нулевыми производными по ϕ и равной нулю ϕ-компонентой скорости.

Следствие второго, третьего и пятого уравнений (1.2) и равенства $dr{\text{/}}dt = {v}$ – уравнение

$\frac{d}{{dt}}\left( {h + \frac{{{{V}^{2}}}}{2} - rg\cos \phi } \right) = 0.$
Введением функции тока ψ оно и четвертое уравнение системы (1.2) интегрируются:
(1.3)
$h + {{V}^{2}}{\text{/}}2 - rg\cos \phi = H(\psi ),\quad s = {{s}_{0}} = S(\psi ),$
с ψ, определяемым дифференциальным равенством (см. [55]–[58])
$d\psi = kr\rho (udr - vdx) = kr\rho V(\cos \theta dr - \sin \theta dx),$
в котором k – произвольный нормирующий множитель. Функции H(ψ) и S(ψ) определяются набегающим потоком: первая – всегда, а вторая – при слабых УВ.

Из четырех независимых дифференциальных уравнений (1.2) проинтегрировались два, для которых линии тока dr/dx = $\operatorname{tg} \theta $ – характеристики. Для получения характеристик, отличных от линий тока, в три первых уравнения (1.2) подставим u = $V\cos \theta $ и ${v}$ = $V\sin \theta $. Производные от V войдут в эти уравнения как производная dV/dt, исключив которую, придем к уравнениям

$\left( {{{M}^{2}} - 1} \right)\left( {\frac{{\partial p}}{{\partial x}} + tg\theta \frac{{\partial p}}{{\partial r}}} \right) - \rho {{V}^{2}}\left( {tg\theta \frac{{\partial \theta }}{{\partial x}} - \frac{{\partial \theta }}{{\partial r}} - \frac{{tg\theta }}{r}} \right) + \rho g\cos \phi \operatorname{tg} \theta = 0,$
$\operatorname{tg} \theta \frac{{\partial p}}{{\partial x}} - \frac{{\partial p}}{{\partial r}} - \rho {{V}^{2}}\left( {\frac{{\partial \theta }}{{\partial x}} + \operatorname{tg} \theta \frac{{\partial \theta }}{{\partial r}}} \right) + \rho g\cos \phi = 0,$
а после добавления к первому из них второго, умноженного на пока неопределенный множитель λ, – к уравнению
$\begin{gathered} (1 - \lambda \operatorname{tg} \theta )\left( {\frac{{\partial \theta }}{{\partial r}} + \frac{{\operatorname{tg} \theta + \lambda }}{{\lambda \operatorname{tg} \theta - 1}}\frac{{\partial \theta }}{{\partial x}}} \right) + \frac{{\left( {{{M}^{2}} - 1} \right)\operatorname{tg} \theta - \lambda }}{{\rho {{V}^{2}}}}\left[ {\frac{{\partial p}}{{\partial r}} + \frac{{{{M}^{2}} - 1 + \lambda \operatorname{tg} \theta }}{{\left( {{{M}^{2}} - 1} \right)\operatorname{tg} \theta - \lambda }}\frac{{\partial p}}{{\partial x}}} \right] + \\ + \;\frac{{\operatorname{tg} \theta }}{r} + \frac{{\operatorname{tg} \theta + \lambda }}{{{{V}^{2}}}}g\cos \phi = 0. \\ \end{gathered} $
Это уравнение будет характеристическим при одинаковых дифференциальных операторах, действующих на θ и р, т.е. при одинаковых множителях при дθ/дх и др/дх. Приравняв их, получим квадратное уравнение для λ, два его решения
$\lambda = \pm \sqrt {{{M}^{2}} - 1} = \pm \operatorname{ctg} \mu $
и уравнения С+- и С–-характеристик с “условиями совместности”
(1.4)
${{C}^{ \pm }}:\frac{{dx}}{{dr}} = \operatorname{ctg} (\theta \pm \mu ),\quad \frac{{d\theta }}{{dr}} \pm \frac{{\operatorname{ctg} \mu }}{{\rho {{V}^{2}}}}\frac{{dp}}{{dr}} + \frac{{\sin \theta \sin \mu }}{{r\sin (\mu \pm \theta )}} \mp \frac{{\operatorname{ctg} (\mu \pm \theta )}}{{{{V}^{2}}}}g\cos \phi = 0.$
Здесь μ – угол Маха, определенный формулой $\operatorname{ctg} \mu = \sqrt {{{M}^{2}} - 1} $ или $\sin \mu $ = 1/M.

Для слабо возмущенных течений p = p0 + δp, ρ = ρ0 + δρ, V = 1 + δV и μ = μ0 + δμ со вторыми слагаемыми, много меньшими первых. Кроме того, |θ| $ \ll $ 1. В силу этого и равенства (1.1) уравнения (1.4) для таких течений примут вид

(1.5)
${{C}^{ \pm }}:\frac{{dx}}{{dr}} = \pm \operatorname{ctg} {{\mu }_{0}} \mp \frac{{\theta \pm \delta \mu }}{{{{{\sin }}^{2}}{{\mu }_{0}}}},\quad \frac{{d\theta }}{{dr}} \pm \frac{{\operatorname{ctg} {{\mu }_{0}}}}{{{{\rho }_{0}}}}\frac{{d\delta p}}{{dr}} + \frac{\theta }{r} + \left( {\theta \mp \frac{{\delta \rho }}{{2{{\rho }_{0}}}}\sin 2{{\mu }_{0}}} \right)\frac{{g\cos \phi }}{{{{{\sin }}^{2}}{{\mu }_{0}}}} = 0.$
Для стратифицированной по y = $r\cos \phi $ атмосферы μ0 = μ0(y).

Из-за малости θ смещение линий тока по r в примыкающем к слабой головной УВ слабо возмущенном течении так мало, что им можно пренебречь. Как следствие этого, можно пренебречь изменениями ψ в правых частях интегралов (1.3). Благодаря чему справедливы формулы

(1.6)
$\delta \rho = \frac{{\delta p}}{{a_{0}^{2}}} = M_{0}^{2}\delta p = \frac{{\delta p}}{{{{{\sin }}^{2}}{{\mu }_{0}}}},\quad \delta \mu = \frac{{\gamma - 1 + 2{{{\sin }}^{2}}{{\mu }_{0}}}}{{{{\rho }_{0}}\sin 2{{\mu }_{0}}}}\delta p.$
Вторая из них получена для совершенного газа с постоянными теплоемкостями и их отношением (показателем адиабаты) γ. Подстановка формулы для δρ в условия совместности из (1.5) приводит их к виду

(1.7)
${{{\text{С}}}^{ \pm }}\,:\quad \frac{{d\theta }}{{dr}} \pm \frac{{\operatorname{ctg} {{\mu }_{0}}}}{{{{\rho }_{0}}}}\frac{{d\delta p}}{{dr}} + \frac{\theta }{r} + \left( {\theta \mp \frac{{\operatorname{ctg} {{\mu }_{0}}}}{{{{\rho }_{0}}}}\delta p} \right)\frac{{g\cos \phi }}{{{{{\sin }}^{2}}{{\mu }_{0}}}} = 0.$

С удалением от ЛА слабо возмущенное течение в кольцевом слое между головной и замыкающей УВ становится “короткой” волной, ибо приращения r при ее пересечениях С–-характеристиками (Δr) малы по сравнению с r, т.е. |Δr|/r $ \ll $ 1. Сказанное поясняет фиг. 1. На ней нарисованы фрагмент короткой волны с тремя УВ (головной SWf, замыкающей SWe и внутренней SWi), приходящими на них С+-характеристиками возмущенного течения, с отрезком С–-характеристики и X = O(1), а также частично залитый цилиндр r = 1, на котором в примерах разд. 2 ставится условие непротекания: ${v}$ = 0. Из-за того, что ось х направлена вверх, а ось r – вправо, система координат xr левая с положительными углами θ при повороте по часовой стрелке.

Фиг. 1.

Схема фрагмента “короткой” волны звукового удара с тремя УВ и приходящими на них С+-характеристиками.

Интегрируя от SWf по отрезку С–-характеристики условие совместности (1.7), пренебрежем изменениями на этом отрезке величин с индексом “нуль” и кубичным по δр/p0 изменением на УВ отвечающего С–-характеристикам инварианта Римана. Если |θ|m – максимум |θ| в примыкающей к SWf ударно-волновой структуре (короткой волне), то в ней с погрешностью, меньшей |θ|mΔr/r $ \ll $ 1, равно нулю возмущение “левого” инварианта Римана

(1.8)
$\theta - \frac{{\operatorname{ctg} {{\mu }_{0}}}}{{{{\rho }_{0}}}}\delta p = 0.$

Насколько погрешность выполнения полученного решения меньше, чем |θ|mΔr/r, зависит от знакопеременного распределения θ поперек короткой волны. Для распределения θ, заданного при r = r0 на фиг. 1, погрешность, обусловленная предпоследним слагаемым в (1.7), меньше примерно втрое. Вклад последнего слагаемого еще меньше. Как видно из (1.4), безразмерное отношение в (1.7) $g{\text{/}}{{\sin }^{2}}{{\mu }_{0}}$ = g°L°/${{(V_{0}^{ \circ }\sin {{\mu }_{0}})}^{2}}$ = g°L°/$a{{_{0}^{ \circ }}^{2}}$. Здесь ускорение свободного падения g отнесено к $V{{_{0}^{ \circ }}^{2}}{\text{/}}L^\circ $, а “градус”, как и ранее, отмечает размерные величины. Отсюда для L° = 30 м и земных условий (g° = 9.8 м/с2, $a_{0}^{ \circ } = 330$ м/с) отношение $g{\text{/}}{{\sin }^{2}}{{\mu }_{0}}$ ≈ 0.003.

Чтобы найти, как эволюционирует при r > r0$ \gg $ 1 осесимметричная (своя на каждой меридиональной плоскости) волна, фрагмент которой изображен на фиг. 1, подставим (1.8) в условие совместности для С+-характеристики из (1.7). Проинтегрировав полученное уравнение, найдем

(1.9)
$\theta (\xi ,r) = \Theta (\xi )\Omega (r,{{r}_{0}}),\quad \delta p(\xi ,r) = \theta (\xi ,r){{\rho }_{0}}\operatorname{tg} {{\mu }_{0}},\quad \Omega (r,{{r}_{0}}) = \sqrt {\frac{{{{r}_{0}}\rho _{0}^{0}\operatorname{ctg} {{\mu }_{0}}}}{{r{{\rho }_{0}}\operatorname{ctg} \mu _{0}^{0}}}.} $
Здесь и далее верхний индекс “нуль” отмечает параметры в общем случае неоднородного (зависящего от y = $r\cos \phi $) невозмущенного набегающего потока при r = r0, ξ – “характеристическая” переменная, постоянная на каждой С+-характеристике возмущенного течения, и Θ(ξ) – известное из расчета среднего поля обтекания ЛА (“начальное”) распределение θ:
(1.10)
$r = {{r}_{0}}:\xi = x({{r}_{0}},\xi ) - {{x}_{0}},\quad \theta ({{r}_{0}},\xi ) = \Theta (\xi ),$
где x0 – начальная координата SWf. На фиг. 1 Θ(ξ) – кусочно-непрерывная кривая.

Подстановка в уравнение С+-характеристик (1.5) δμ из (1.6), а затем θ и δр из (1.9) приводит к уравнению

(1.11)
$\frac{{dx}}{{dr}} = \operatorname{ctg} (\theta + \mu ) = \operatorname{ctg} {{\mu }_{0}} - 2\Theta (\xi )\Phi (r,{{r}_{0}}),\quad \Phi (r,{{r}_{0}}) = \frac{{\gamma + 1}}{{{{{\sin }}^{2}}2{{\mu }_{0}}}}\Omega (r,{{r}_{0}}),$
справедливому на любой С+-характеристике от r = r0 до ее прихода на одну из УВ. При рассмотрении внутренних УВ, с обеих сторон от которых, как в случае SWi, поток возмущен, переменную ξ будем отмечать нижним индексом “минус” (“плюс”) для С+-характеристик, приходящих на УВ сверху (снизу) по течению. С учетом такого выбора, проинтегрировав уравнение (1.11) при постоянном ξ– от сечения r = r0 до УВ, придем к равенству
${{x}_{{SW}}} = {{x}_{{SW0}}} + {{\xi }_{ - }} + \int\limits_{{{r}_{0}}}^r {\operatorname{ctg} {{\mu }_{0}}(\bar {r}\cos \phi )d\bar {r}} - 2\Theta ({{\xi }_{ - }})\int\limits_{{{r}_{0}}}^r {\Phi (\bar {r},{{r}_{0}})d\bar {r}} ,$
в котором r и xSW – координаты УВ. Продифференцировав его по r, получим
(1.12)
$\frac{{d{{x}_{{SW}}}}}{{dr}} = \operatorname{ctg} {{\mu }_{0}} + \left[ {1 - 2\Theta {\kern 1pt} '({{\xi }_{ - }})\int\limits_{{{r}_{0}}}^r {\Phi (\bar {r},{{r}_{0}})d\bar {r}} } \right]\frac{{d{{\xi }_{ - }}}}{{dr}} - 2\Theta ({{\xi }_{ - }})\Phi (r,{{r}_{0}}).$
Здесь и далее “штрих” означает дифференцирование по ξ– или по ξ+.

Для слабых УВ, догоняющих С+-характеристики и догоняемых характеристиками того же семейства (см. [56], [60]),

(1.13)
$\frac{{d{{x}_{{SW}}}}}{{dr}} = \frac{{\operatorname{ctg} ({{\theta }_{ - }} + {{\mu }_{ - }}) + \operatorname{ctg} ({{\theta }_{ + }} + {{\mu }_{ + }})}}{2},$
а в силу равенства (1.11)
(1.14)
$\frac{{d{{x}_{{SW}}}}}{{dr}} = \operatorname{ctg} {{\mu }_{0}} - \left[ {\Theta ({{\xi }_{ - }}) + \Theta ({{\xi }_{ + }})} \right]\Phi (r,{{r}_{0}}).$
Разность (1.12) и (1.14) дает уравнение
$\left[ {1 - 2\Theta {\kern 1pt} {\text{'}}{\kern 1pt} ({{\xi }_{ - }})\int\limits_{{{r}_{0}}}^r {\Phi (\bar {r},{{r}_{0}})d\bar {r}} } \right]\frac{{d{{\xi }_{ - }}}}{{dr}} = [\Theta ({{\xi }_{ - }}) - \Theta ({{\xi }_{ + }})]\Phi (r,{{r}_{0}}).$
Введя переменную χ с начальным условием χ(r0) = 0, заменим его на два:
(1.15)
$\frac{{d\chi }}{{dr}} = 2\Phi (r,{{r}_{0}}),\quad \frac{{d{{\xi }_{ - }}}}{{dr}} = \frac{{\Theta ({{\xi }_{ - }}) - \Theta ({{\xi }_{ + }})}}{{1 - \chi \Theta '({{\xi }_{ - }})}}\Phi (r,{{r}_{0}}).$
Аналогичным образом с той же функцией χ получим уравнение для ξ+

(1.16)
$\frac{{d{{\xi }_{ + }}}}{{dr}} = \frac{{\Theta ({{\xi }_{ + }}) - \Theta ({{\xi }_{ - }})}}{{1 - \chi \Theta {\text{'}}({{\xi }_{ + }})}}\Phi (r,{{r}_{0}}).$

Итак, каждая внутренняя УВ определяется решением задачи Коши для обыкновенных дифференциальных уравнений первого порядка: (1.15) и (1.16) – для ξ– и ξ+, (1.14) – для xSW и первого уравнения (1.15) – для χ (одного на все УВ) c начальными условиями ξ– = ξ+ = xSW = xSW0 и χ = 0 при r = r0. По найденным ξ–(r) и ξ+(r) величины θ± и δp± с обеих сторон от УВ определят формулы (1.9), а δρ± и δμ± – формулы (1.6). Для головной УВ, поток перед которой не возмущен, остаются три уравнения (1.14) и (1.16) с Θ(ξ–) ≡ 0. Если не возмущен поток за замыкающей УВ, то также остаются три уравнения, но с Θ(ξ+) ≡ 0.

При однородном набегающем потоке уравнение для χ интегрируется и дает

$\chi = \frac{{4(\gamma + 1){{r}_{0}}}}{{{{{\sin }}^{2}}2{{\mu }_{0}}}}\left( {\sqrt {r{\text{/}}{{r}_{0}}} - 1} \right).$
Если при r = r0 за головной УВ Θ(ξ+) = α – βξ+ с положительными константами α и β, то, проинтегрировав уравнения (1.12) и (1.14) с Θ(ξ–) = 0, для SWf придем к решениям и к их следствиям для r/r0$ \gg $ 1:
(1.17)
$\begin{gathered} \frac{{\delta p(r)}}{{\delta p({{r}_{0}})}} = \frac{{\theta (r)}}{{\theta ({{r}_{0}})}} = \frac{{\sqrt {{{r}_{0}}{\text{/}}r} }}{{\sqrt {1 + \beta \chi } }}\;\mathop \approx \limits_{r/{{r}_{0}} \gg 1} \;\frac{{\sin 2{{\mu }_{0}}}}{{2\sqrt {(\gamma + 1)\beta {{r}_{0}}} }}{{\left( {\frac{{{{r}_{0}}}}{r}} \right)}^{{3/4}}}, \\ {{x}_{{SW}}} - {{x}_{{SW0}}} - (r - {{r}_{0}})\operatorname{ctg} {{\mu }_{0}} = - \frac{\alpha }{\beta }(\sqrt {1 + \beta \chi } - 1)\;\mathop \approx \limits_{r/{{r}_{0}} \gg 1} \; - \,{\kern 1pt} \frac{{2\alpha \sqrt {(\gamma + 1){{r}_{0}}} }}{{\sqrt \beta \sin 2{{\mu }_{0}}}}{{\left( {\frac{r}{{{{r}_{0}}}}} \right)}^{{1/4}}}. \\ \end{gathered} $
Первым закон (1.17) затухания головной УВ установил О.С. Рыжов (см. [10]), а справедливую для r/r0$ \gg $ 1 степенную зависимость δp ∼ 1/r3/4 – еще Л.Д. Ландау (см. [1]).

Обычны ситуации, когда одна УВ (“вторая” – SWII) догоняет другую (“первую” – SWI), например, головную. Так как обе УВ слабые, то при их слиянии интенсивность единственной результирующей УВ определят величины θ– перед SWI и θ+ за SWII в момент слияния. При численном решении описанных выше задач Коши момент слияния (значение r = ${{r}_{{{\text{I}} + {\text{II}}}}}$, при котором оно происходит) определит равенство xSWII = xSWI.

Интегрирование уравнения (1.11) для одной-двух “промежуточных” характеристик позволяет построить распределения θ по х между УВ при r = ${{r}_{{{\text{I}} + {\text{II}}}}}$, т.е. новые Θ(ξ) для продолжения счета при r > ${{r}_{{{\text{I}} + {\text{II}}}}}$. Это может понадобиться для представления волн сжатия с их “стремлением опрокинуться” и лишь иногда для еще искривленных волн разрежения. Те же волны разрежения, которые были прямолинейными при r = r0, останутся такими и далее. Для них промежуточных характеристик считать не надо. Описанный подход назовем “приближением короткой волны” (ПКВ). В связи с газодинамическими приложениями термин “короткая волна” введен С.А. Христиановичем и соавт. (см. [7]). Первые результаты по затуханию таких волн, не называя их “короткими”, получил Л.Д. Ландау (см. [1]).

Обыкновенные дифференциальные уравнения (1.14)–(1.16) ПКВ – простые следствия уравнений Эйлера (1.2). Этим ПКВ отличается от развитого в [15], [16] подхода, также сводящего расчет дальнего поля звукового удара к решению обыкновенных дифференциальных уравнений.

2. ПРИМЕРЫ РАСЧЕТА УДАРНО-ВОЛНОВЫХ СТРУКТУР В СТАЦИОНАРНЫХ СВЕРХЗВУКОВЫХ ТЕЧЕНИЯХ

Ниже для однородной атмосферы приводятся результаты расчета стационарных ударно-волновых структур в идеальном совершенном газе с γ = 1.4 и М0 = 2. Для слабо возмущенных осесимметричных сверхзвуковых течений полученные результаты демонстрируют возможность аккуратного расчета таких структур на любых расстояниях от источника возмущений со сколь-угодно малыми возмущениями параметров. Расчет таких структур велся по явной маршевой в направлении оси х схеме с построением по уравнению (1.13) головной и замыкающей УВ-адаптированных к рассчитываемому течению сеточных линий и в рамках ПКВ, развитого в разд. 1 для r $ \gg $ 1.

При маршевом счете на слоях x = const число K = ${{2}^{{k--1}}}$ разностных ячеек между двумя УВ было фиксировано, а между замыкающей УВ и цилиндром r = r0 в начале счета увеличивалось от 1 до K, а затем сохранялось неизменным. В итоге общее число ячеек на слоях x = const становилось равным 2k при равномерном разбиении по r в каждой из двух областей. Решение на границах всех ячеек задачи взаимодействия равномерных сверхзвуковых потоков (фиг. 2) – аналога нестационарной задачи о распаде разрыва здесь упрощается из-за малого отличия параметров взаимодействующих потоков. Подобно нестационарной задаче [59], благодаря этому, условия, справедливые для слабых УВ и центрированных волн разрежения, записываются единообразно через разности р и θ. В приближении слабых УВ на разделительных линиях тока непрерывны все параметры. В точной постановке эти линии тока – тангенциальные разрывы, на которых непрерывны только р и θ.

Фиг. 2.

Схема взаимодействия сверхзвуковых потоков для определения параметров на левой и правой границах (LB и RB) ячейки (LW – левая волна, RW – правая волна, LS – линия тока).

Согласно условиям совместности из (1.4), для совершенного газа на LW и RW – “левой” и “правой” слабых ударных или центрированных волнах – с равной погрешностью (кубичной по приращению давления) выполняются равенства

(2.1)
$\begin{gathered} LW\,:\quad \Theta = {{\theta }_{{II}}} + \left[ {\frac{{G({{M}_{{II}}})}}{{{{p}_{{II}}}}} + \frac{{G(M)}}{P}} \right]\frac{{{{p}_{{II}}} - P}}{{2\gamma }}, \\ RW\,:\quad \Theta = {{\theta }_{I}} - \left[ {\frac{{G({{M}_{I}})}}{{{{p}_{I}}}} + \frac{{G(M)}}{P}} \right]\frac{{{{p}_{I}} - P}}{{2\gamma }},\quad G(M) = \frac{{\sqrt {{{M}^{2}} - 1} }}{{{{M}^{2}}}}. \\ \end{gathered} $
Здесь и ниже индексы “I” и “II”отмечают параметры взаимодействующих потоков, полученные интерполяцией в согласии с [51] на границах, отличных от SWf и SWe, экстраполяцией для SWe и за SWf и с заменой индекса “I” на “0” перед SWf, а “большие величины” без индексов P, Θ, R, A и М – значения p, θ, ρ, a и числа Маха между LW и RW. Исключив Θ из уравнений (2.1), придем к квадратному уравнению для Р и к его корню
(2.2)
$\begin{gathered} P = {{\left( {C + \sqrt {{{C}^{2}} + BF} } \right)} \mathord{\left/ {\vphantom {{\left( {C + \sqrt {{{C}^{2}} + BF} } \right)} B}} \right. \kern-0em} B}, \\ B = \frac{{G({{M}_{I}})}}{{{{p}_{I}}}} + \frac{{G({{M}_{{II}}})}}{{{{p}_{{II}}}}},\quad C = \frac{{G({{M}_{I}}) + G({{M}_{{II}}})}}{2} - G(M) + \gamma ({{\theta }_{{II}}} - {{\theta }_{I}}),\quad F = ({{p}_{I}} + {{p}_{{II}}})G(M). \\ \end{gathered} $
Зависимость от Р числа Маха М, входящего в C и F, определяют формулы

$R = \frac{{{{P}^{{1/\gamma }}}}}{{p_{0}^{{1/\gamma }}}},\quad {{A}^{2}} = \frac{{\gamma P}}{R},\quad {{M}^{2}} = \frac{{{{I}_{0}} - 2{{A}^{2}}}}{{(\gamma - 1){{A}^{2}}}},\quad {{I}_{0}} = \frac{{2 + (\gamma - 1)M_{0}^{2}}}{{M_{0}^{2}}}.$

Уравнение (2.2) решалось итерациями с Р(0) = pII, после сходимости которых R, A и M известны, V = AM, а Θ определяет любое из уравнений (2.1). При решении задач взаимодействия для всех ячеек сначала на полуслое x = xn + Δх/2 находятся координаты r выделяемых УВ. По ним при равномерном разбиении по r отрезков между SWf и SWe и между SWe и цилиндром r = 1 находятся наклоны всех границ и их координаты на полуслое. По наклонам границ устанавливается, куда в задачах взаимодействия они попадают, и определяются параметры течения на этих границах (“большие величины” в законах сохранения, записанных для каждой двумерной ячейки между слоями x = xn и x = xn + Δx/2) и делается полушаг, в результате которого находятся параметры во всех ячейках полуслоя. По найденным параметрам на полуслое определяются параметры с индексами “I” и “II”, решаются новые задачи взаимодействия и делается полный шаг до хn = xn + Δx. Не вдаваясь в дальнейшие детали, отметим, что, согласно расчетам, максимальное число Куранта, с которым созданный алгоритм реализует устойчивый счет, равно $0.8\operatorname{ctg} {{\mu }_{0}}$.

Начальное возмущение (фиг. 3а) однородного потока с М0 = 2, V0 = ρ0 = 1, ${{a}_{0}} = 1{\text{/}}{{M}_{0}} = 0.5$, ${{p}_{0}} = 1{\text{/}}(\gamma M_{0}^{2}) \approx 0.179$ задавалось в кольце x = х0 = 1, 1 ≤ r ≤ r0, что близко к фиг. 1, но с r0 = 1 + $\operatorname{tg} {{\mu }_{0}}$ ≈ ≈ 1.58. Максимальные θm ≈ 0.01 у передней границы кольца (SWf) и на его среднем радиусе (SWi), минимальное θ = –θm на замыкающей границе (SWe) и линейные по r между ними. Начальное возмущение р принималось равным δp = θ$\operatorname{tg} {{\mu }_{0}}$ ≈ 0.58θ, начальная плотность определялась по р условием p/ργ = p0, а модуль скорости газа – условием постоянства полной энтальпии. При θm ≈ 0.01 максимальные |δp|m/p0 = 0.032. Это настолько мало, что энтропия и “левый” инвариант Римана в УВ практически не изменяются. При x ≥ 1 на цилиндре r = 1 выполнялось условие непротекания.

Фиг. 3.

Начальные распределения δp/p0 и θ (см. текст).

На фиг. 4 приведены распределения δp/p0, в которые начальное возмущение превратилось при х = 10. Сравниваются результаты расчетов при разном числе ячеек 2k в радиальном направлении с выделением SWf и SWe – двух границ слегка расширяющейся узкой конически-кольцевой волны между ними. На каждом фрагменте фиг. 4 черные (красные) кривые – результаты расчетов с 2k (${{2}^{{k--1}}}$) ячейками. На фиг. 5 таким же способом сравниваются N-волны, в которые при разном числе ячеек это возмущение превратилось при х = 400. Расчет уже на 32 ячейках в столь дальнем поле дает близкую к точной N-волну за 50 с работы одноядерного ПК. На N-волну здесь приходится всего 16 ячеек, а их размер по r при х = 400 в широкой приосевой зоне больше, чем между SWf и SWe в 100 раз.

Фиг. 4.

Кривые δp/p0 = f(r) при x = 10, рассчитанные маршем на сетках с 2k (черные кривые) и с ${{2}^{{k--1}}}$ (красные кривые) ячейками в слоях x = const.

Фиг. 5.

N-волны при счете маршем сверхзвукового осесимметричного течения до x = 400 с 2k (черные кривые) и ${{2}^{{k--1}}}$ (красные кривые) ячейками в слоях x = const.

Обычно замыкающая УВ возникает не сразу, а на некотором удалении от ЛА. При счете маршем ее заменит УВ нулевой интенсивности, т.е. С+-характеристика, которая с ростом х станет замыкающей УВ. Как это происходит, показывает пример с начальными параметрами без замыкающей УВ (фиг. 3б). Главные их отличия от рассмотренных – при $1 \leqslant r < 1 + 1{\text{/}}(2\sqrt 3 )$ минимум $\theta = - {{\theta }^{m}}{\text{/}}\sqrt 2 $ в вершине $r = 1 + (\sqrt 2 - 1){\text{/}}(2\sqrt 6 )$ равнобедренного треугольника, равного по площади треугольнику, примыкающему к внутренней УВ. Кроме того, в этом примере θm ≈ 0.005, т.е. вдвое меньше, чем в предыдущем примере.

Описание формирования замыкающей УВ требует большего числа ячеек, чем расчет, когда такая УВ есть изначально. Получившиеся в данном случае профили δp/p0 = f(r) для x = 5 и 10 при счете на сетках с k = 6, 7 и 8 представлены на фиг. 6. При х = 10 замыкающая УВ полностью сформировалась, и продолжать счет на столь тонких сетках нет необходимости. Для дальнейшего счета достаточно сетки с k = 5.

Фиг. 6.

Кривые δp/p0 = f(r) для x = 5 и 10, рассчитанные на трех сетках при начальных параметрах без замыкающей УВ.

Последние примеры относятся к описанию звукового удара с удалений в несколько десятков длин ЛА, когда применимо ПКВ. Для смещенных по оси r на Δr ≈ 19 начальных распределений δp/р0 (на фиг. 3а это отвечает r0 ≈ 20.6) сравнивались результаты ПКВ, счета установлением и маршем. В данном примере при M0 = 2 в силу равенства (1.8) δp/р0 ≈ 3.2θ. Сравнение параметров в сечении слияния SWf и SWi (r ≈ 103, x ≈ 141) и N-волн при x = 200, 300 и 400 показало, что для не столь большой величины r0 ≈ 21 погрешности ПКВ в определении θ и δр не превышают 1%, а в определении координат N-волн – 0.005%.

После такого подтверждения точности ПКВ покажем его возможности на начальных распределениях, более сложных, чем на фиг. 3а и 3б. Два таких распределения θ = θ(х, r0) = Θ(ξ) с константами 0 < θm$ \ll $ 1 и Х > 1 приведены на фиг. 3в. Первое образуют отрезки прямых, из которых вертикальные – УВ, а наклонные – волны разрежения. Горизонтальный участок 1+-2– введен из-за интереса к таким участкам, например, в [26], [29], [32], [34], [38], [43]. Во второе распределение введена красная волна сжатия при ξ > Х, а начальная ломаная 1–-1+-2±-3– заменена гладкой красной волной сжатия-разрежения 1-2-3– без головной УВ (с “головной С+-характеристикой” $C_{1}^{ + }$). Нечто похожее на такую волну есть в [41], [42].

Первое начальное распределение не содержит непрерывных волн сжатия, а волны разрежения и горизонтальный участок 1+-2– – отрезки прямых. С ростом r они остаются прямыми. К значениям ξ±, θ± и координат х всех УВ как функций r, которые определяются уравнениями (1.14)–(1.16), в данном случае добавляется х2± – координата С+-характеристики $C_{2}^{ + }$, вышедшей из точки 2±. Координаты х любых чем-то интересных С+-характеристик можно определять либо интегрированием уравнения (1.11) с постоянным ξ = ξ2± (в случае точки 2±), либо для упрощения логики программы – как УВ нулевой интенсивности.

Результаты расчетов для первого распределения θ(х, r0) = Θ(ξ) (фиг. 3в) с r0 = 20, X = 2 и θm = 0.04 приведены на фиг. 7. Наряду с начальным профили θ = θ(х) даны для r, при котором $C_{2}^{ + }$ догнала SW1, для тех r, при которых одна из УВ догнала другую (по координате х – “предыдущую”) и слилась с ней, и для r = 400. За УВ, получившейся при слиянии, сохранялся номер предыдущей. В рассчитанном примере после четырех слияний остались три УВ: SW1, SW6 и SW7. Интенсивность SW6 настолько мала, что она SW1 уже не догонит, т.е. сохраняется слегка отличная от N-волны структура с третьей очень слабой УВ (фиг. 7е).

Фиг. 7.

Профили θ = θ(х): начальный r = r0 = 20 (а), при слияниях SW1 с $C_{2}^{ + }$ r = 28.71 (б), с SW3 r = 38.17 (в), с SW4 r = 38.96 (г) и с SW5 r = 39.85 (д) и при r = 400 (е).

В примере с волной сжатия-разрежения 1-2-3–, головной SW1 нулевой интенсивности (с “головной С+-характеристикой” $C_{1}^{ + }$) и волной сжатия за SW7 уточним некоторые моменты. При описании эволюции волны 1-2-3– находятся координата х1 ее начальной точки, координата х2 и значение θ2 максимума и те же величины в ее концевой точке (х3 и θ3–). Начальная точка и точка максимума принадлежат фиксированным С+-характеристикам, переменные ξ которых изменяются только при пересчетах, вызванных слияниями УВ или образованием “висячей” УВ на $C_{1}^{ + }$. Значение r, при котором это случится, зависит от величины $\Theta '(0)$. Для определения такого r и описания движения SW3 разобьем волну 1-2-3– на два участка и, согласно фиг. 3в, начальные Θ(ξ) на них зададим многочленами

(2.3)
$\begin{gathered} 0 \leqslant \xi \leqslant {{\xi }_{2}} = 0.22X,\quad \Theta (\xi ) = {{\beta }_{1}}\xi + {{\beta }_{2}}{{\xi }^{2}} + {{\beta }_{3}}{{\xi }^{3}},\quad \Theta {\text{'}}(\xi ) = {{\beta }_{1}} + 2{{\beta }_{2}}\xi + 3{{\beta }_{3}}{{\xi }^{2}}; \\ {{\xi }_{2}} \leqslant \xi \leqslant {{\xi }_{3}} = 0.34X,\quad \Theta (\xi ) = {{\beta }_{4}} + {{\beta }_{5}}{{(\xi - {{\xi }_{2}})}^{2}},\quad \Theta {\text{'}}(\xi ) = 2{{\beta }_{5}}(\xi - {{\xi }_{2}}), \\ {{\beta }_{1}} = 3.3\frac{{{{\theta }^{m}}}}{X},\quad {{\beta }_{2}} = - 3.6\frac{{{{\theta }^{m}}}}{{{{X}^{2}}}},\quad {{\beta }_{3}} = - 12\frac{{{{\theta }^{m}}}}{{{{X}^{3}}}},\quad {{\beta }_{4}} = {{\beta }_{1}}{{\xi }_{2}} + {{\beta }_{2}}\xi _{2}^{2} + {{\beta }_{3}}\xi _{2}^{3}, \\ {{\beta }_{5}} = \frac{{{{\beta }_{1}} + {{\beta }_{2}}({{\xi }_{3}} + {{\xi }_{2}}) + {{\beta }_{3}}\left( {\xi _{3}^{2} + {{\xi }_{3}}{{\xi }_{2}} + \xi _{2}^{2}} \right)}}{{{{\xi }_{3}} - {{\xi }_{2}}}}. \\ \end{gathered} $
Поскольку Θ(ξ+ = 0) = 0, то ξ+ = 0 в силу уравнения (1.16) с Θ(ξ–) ≡ 0, пока не обратится в нуль знаменатель правой части этого уравнения

(2.4)
$1 - \Theta {\text{'}}(0)\chi ({{r}_{ * }},{{r}_{0}}) = 1 - {{\beta }_{1}}\chi ({{r}_{ * }},{{r}_{0}}) = 0.$

При x = x(ξ+, r) на $C_{1}^{ + }$ найдем, как на ней с ростом r изменяется производная ∂x/∂ξ+. Продифференцировав уравнение (1.11) по ξ при постоянном r и преобразовав результат с учетом первого уравнения (1.15), придем к задаче

$\frac{{d(\partial x{\text{/}}\partial {{\xi }_{ + }})}}{{d\chi }} = - \Theta {\text{'}}(0),\quad r = {{r}_{0}},\quad \chi = 0,\quad {{(\partial x{\text{/}}\partial {{\xi }_{ + }})}_{{r = {{r}_{0}}}}} = 1$
с начальным условием для ∂x/∂ξ+, вытекающим из определения ξ в (1.10), и с производной по χ при ξ+ = const = 0. Решение этой задачи – равенство
(2.5)
$\partial x{\text{/}}\partial {{\xi }_{ + }} = 1 - \Theta {\text{'}}(0)\chi (r,{{r}_{0}}),$
в котором $\Theta {\text{'}}(0) > 0$, а χ – монотонно растущая функция r. Поэтому при $r = {{r}_{ * }}$ из (2.4) ∂x/∂ξ+ становится нулем на $C_{1}^{ + }$ (если раньше ее не догнала SW3), свидетельствуя о зарождении SW1. Наконец, из (1.9) и (2.5) следует, что на $C_{1}^{ + }$
$\frac{{\partial \theta ({{\xi }_{ + }},r)}}{{\partial x}} = \frac{{\Theta {\text{'}}(0)\Omega (r,{{r}_{0}})}}{{1 - \Theta {\text{'}}(0)\chi (r,{{r}_{0}})}}.$
Таким образом, при том же $r = {{r}_{ * }}$ эта производная обращается в бесконечность.

При $r = {{r}_{ * }}$ пересчитаем ξ и Θ(ξ): на отрезках 1-2 и 2-3– по известным координатам х1, х2 и х3 и значениям θ2 и θ3– найдем ξ1N = х1, ξ2N = х2 и ξ3N = х3, а “старые” формулы (2.3) заменим на “новые”:

${{\xi }_{{1N}}} \leqslant \xi \leqslant {{\xi }_{{2N}}}\,:\quad \Theta (\xi ) = {{\beta }_{{1N}}}\sqrt {\xi - {{\xi }_{{1N}}}} + {{\beta }_{{2N}}}(\xi - {{\xi }_{{1N}}}),\quad \Theta {\text{'}}(\xi ) = \frac{{{{\beta }_{{1N}}}}}{{2\sqrt {\xi - {{\xi }_{{1N}}}} }} + {{\beta }_{{2N}}};$
${{\xi }_{{2N}}} \leqslant \xi \leqslant {{\xi }_{{3N}}}\,:\quad \Theta (\xi ) = {{\theta }_{2}} + {{\beta }_{{3N}}}{{(\xi - {{\xi }_{{2N}}})}^{2}},\quad \Theta {\text{'}}(\xi ) = 2{{\beta }_{{3N}}}(\xi - {{\xi }_{{2N}}}),$
${{\beta }_{{1N}}} = \frac{{2{{\theta }_{2}}}}{{\sqrt {{{\xi }_{{2N}}} - {{\xi }_{{1N}}}} }},\quad {{\beta }_{{2N}}} = \frac{{ - {{\theta }_{2}}}}{{{{\xi }_{{2N}}} - {{\xi }_{{1N}}}}},\quad {{\beta }_{{3N}}} = \frac{{{{\theta }_{{3 - }}} - {{\theta }_{2}}}}{{{{{(\xi - {{\xi }_{{2N}}})}}^{2}}}}.$
Если же точки 2 на 1-3– уже нет (SW3 догнала $C_{2}^{ + }$), то формулы (2.3) заменим на
${{\xi }_{{1N}}} \leqslant \xi \leqslant {{\xi }_{{3N}}}\,:\quad \Theta (\xi ) = {{\beta }_{{1N}}}\sqrt {\xi - {{\xi }_{{1N}}}} ,\quad \Theta {\text{'}}(\xi ) = \frac{{{{\beta }_{{1N}}}}}{{2\sqrt {\xi - {{\xi }_{{1N}}}} }},\quad {{\beta }_{{1N}}} = \frac{{{{\theta }_{{3 - }}}}}{{\sqrt {{{\xi }_{{3N}}} - {{\xi }_{{1N}}}} }}.$
Кроме того, здесь “старое” r0 заменяется на ${{r}_{ * }}$, а при слияниях УВ – на ${{r}_{{{\text{I}} + {\text{II}}}}}$, и интегрирование уравнений (1.15) возобновляется с начальным условием χ(r0,r0) = 0.

При $r = {{r}_{ * }}$ в правой части уравнения (1.16):

$\frac{{d{{\xi }_{ + }}}}{{dr}} = \frac{{\Theta ({{\xi }_{ + }})\Phi (r,{{r}_{0}})}}{{1 - \chi \Theta {\text{'}}({{\xi }_{ + }})}},\quad \Theta ({{\xi }_{ + }} \to {{\xi }_{{1N}}}) \to 0,\quad \chi (r \to {{r}_{0}},{{r}_{0}}) \to 0,\quad \Theta {\text{'}}({{\xi }_{ + }} \to {{\xi }_{{1N}}}) \to \infty ,$
возникли неопределенности. Для их разрешения воспользуемся его следствиями:
$\frac{{d{{\xi }_{ + }}}}{{d\chi }} = \frac{{\Theta ({{\xi }_{ + }})}}{{2 - 2\Theta {\text{'}}({{\xi }_{ + }})\chi }} \Leftrightarrow \frac{{d\chi }}{{d{{\xi }_{ + }}}} + 2\frac{{\Theta {\text{'}}({{\xi }_{ + }})}}{{\Theta ({{\xi }_{ + }})}}\chi - \frac{2}{{\Theta ({{\xi }_{ + }})}} = 0.$
Решив это линейное уравнение, найдем
$\chi = \frac{2}{{{{\Theta }^{2}}({{\xi }_{ + }})}}\int\limits_{{{\xi }_{{1N}}}}^{{{\xi }_{ + }}} {\Theta \left( {{{{\bar {\xi }}}_{ + }}} \right)d{{{\bar {\xi }}}_{ + }}} \approx \frac{2}{{{{\beta }_{{1N}}}({{\xi }_{ + }} - {{\xi }_{{1N}}})}}\int\limits_{{{\xi }_{{1N}}}}^{{{\xi }_{ + }}} {\sqrt {{{{\bar {\xi }}}_{ + }} - {{\xi }_{{1N}}}} d{{{\bar {\xi }}}_{ + }}} = \frac{{4\sqrt {{{\xi }_{ + }} - {{\xi }_{{1N}}}} }}{{3{{\beta }_{{1N}}}}}.$
Таким образом, для r, близких к новому ${{r}_{0}} = {{r}_{ * }}$, имеем
${{\xi }_{ + }} - {{\xi }_{{1N}}} \approx \frac{9}{{16}}\beta _{{1N}}^{2}{{\chi }^{2}}.$
С той же точностью χ ≈ 2Φ(r0, r0)(r – r0), и с учетом формулы для Φ(r0, r0) получим

$r - {{r}_{0}} \ll 1\,:\quad {{\xi }_{ + }} - {{\xi }_{{1N}}} \approx \frac{9}{4}\beta _{{1N}}^{2}{{\Phi }^{2}}({{r}_{0}},{{r}_{0}}){{(r - {{r}_{0}})}^{2}},\quad \frac{{d{{\xi }_{ + }}}}{{dr}} \approx 9\frac{{\beta _{{1N}}^{2}{{{(\gamma + 1)}}^{2}}}}{{2{{{\sin }}^{4}}2{{\mu }_{0}}}}(r - {{r}_{0}}).$

Результаты расчетов для второго распределения θ(х, r0) = Θ(ξ) (фиг. 3в) с r0 = 20, X = 2 и θm = 0.01 приведены на фиг. 8. Наряду с начальным, профили θ = θ(х) даны для значения $r = {{r}_{ * }}$ образования висячей SW1 на $C_{1}^{ + }$ для тех r, при которых SW3, SW4 и SW5 одна за другой догоняли SW1, сливаясь с ней, и для r = 400. Из-за вчетверо меньшего θm слияния УВ происходили при заметно бóльших r, чем в примере фиг. 7. При интегрировании обыкновенных дифференциальных уравнений методом Рунге–Кутты до r = 103 времена счета в этих примерах – доли секунды на ПК.

Фиг. 8.

Профили θ = θ(х): начальный r = r0 = 20 (а), при образовании SW1 r = 30.59 (б), при слияниях SW1 с SW3 r = 86.92 (в), с SW4 r = 106.69 (г), с SW5 r = 118.36 (д) и при r = 400 (е).

УВ, рассматриваемые как поверхности разрыва, на самом деле имеют конечную ширину $\Delta ^\circ $ (“градус” по-прежнему отмечает размерные величины). Для УВ умеренной и большой интенсивности она пренебрежимо мала (равна нескольким длинам свободного пробега). Однако размазывание слабых УВ с $\Delta ^\circ $, обратно пропорциональным разности $p_{ + }^{ \circ } - p_{ - }^{ \circ } = \delta p_{ + }^{ \circ } - \delta p_{ - }^{ \circ }$, может изменить профиль давления в звуковом ударе и, как следствие, – высокочастотную часть его спектра (см. [60]). Решение из [55], определяющее в рамках одномерных стационарных уравнений Навье–Стокса профиль δр размазанной УВ и ее ширину Δ, можно записать так:

(2.6)
$\begin{gathered} \delta p = \frac{{\delta {{p}_{ + }} + \delta {{p}_{ - }}}}{2} + \frac{{\delta {{p}_{ + }} - \delta {{p}_{ - }}}}{2}\operatorname{th} \frac{n}{\Delta },\quad \Delta = \frac{{\Delta ^\circ }}{{{{L}^{ \circ }}}} = \frac{{4\gamma [4Pr + 3(\gamma - 1)]{{M}_{0}}}}{{3(\gamma + 1)(\delta {{p}_{ + }} - \delta {{p}_{ - }})PrRe}}, \\ Re = \frac{{L^\circ V_{0}^{ \circ }}}{{\nu ^\circ }},\quad Pr = \frac{{\nu ^\circ }}{{\kappa ^\circ }}\rho _{0}^{ \circ }c_{p}^{ \circ },\quad \nu ^\circ = \frac{{4\eta ^\circ + 3\zeta ^\circ }}{{4\rho _{0}^{ \circ }}}. \\ \end{gathered} $
Здесь координата n растет в направлении потока от “середины” УВ по нормали к ней, δр отнесены к $p_{0}^{ \circ }$, $\eta ^\circ $ и $\zeta ^\circ $ – коэффициенты “первой” и “второй” вязкости, $\kappa ^\circ $ – теплопроводность воздуха, а Re и Pr – числа Рейнольдса и Прандтля.

При полете с $M_{0}^{0} = 2$ на высоте 20 км в рамках “Стандартной атмосферы” у Земли М0 ≈ 1.73, и формула для Δ при Pr = 0.8 и длине l ЛА в метрах принимает вид Δ ≈ 1.7/[l(δp+ – δp–)106]. Отсюда для l = 30 ширина Δ крайне слабой УВ с приращением давления всего в 1 Па не превысит 6 × 10–3, что много меньше равной нескольким единицам ширины дошедшей до Земли волны звукового удара. Еще меньшее размазывание УВ при распространении от ЛА до Земли практически не может влиять на их эволюцию. При расчете громкости звукового удара эффекты вязкого размазывания слабых УВ правильно и легко учесть с помощью решения (2.6). Воспользовавшись уравнением Ландау–Теллера, еще проще учесть размазывание слабых УВ из-за неравновесности колебательных степеней свободы молекул воздуха. Это на порядки быстрее и много проще, чем численное решение уравнения Бюргерса в том же, что и ПКВ, осесимметричном приближении (см. [39]).

ЗАКЛЮЧЕНИЕ

До настоящего времени инструменты предсказания звукового удара вдали от предполагаемых высоколетящих сверхзвуковых гражданских ЛА во многом опираются на весьма сложные методы и решения линейной и нелинейной акустики движущихся неоднородных сред. Отмеченная сложность проявляется уже в процессе получения применяемых затем уравнений и решений.

Упрощение аналитической части данного исследования обязано интеграции численного и аналитического подходов. “Аналитика” подсказывает правильную адаптацию разностных сеток, а “вычислительный инструмент” ограничивает применение аналитики расстояниями, начиная с которых исследуемое течение – слабая короткая волна при том, что с удалением от ЛА в каждой меридиональной плоскости она близка к осесимметричной. Учет этих особенностей позволил начинать анализ не с пространственных, а с более простых осесимметричных уравнений и затем работать с уравнениями их характеристик и решениями последних для слабых коротких волн. Существенно, что при таком подходе эффекты осевой симметрии, нелинейности и стратификации невозмущенной атмосферы учитываются не последовательно, а одновременно. Еще одно упрощение связано со скоростью слабых УВ как полусуммы скоростей приходящих к ним характеристик. Вместо этого большинство исследователей, начиная с Л.Д. Ландау (см. [1]), пользуются более сложным “правилом площадей” (исключение – К.Е. Губкин, см. [5]).

В свете выполненного исследования и в развитие полученных ранее результатов можно надеяться на реальность многократного сокращения времен расчета в приближении пространственных уравнений Эйлера ближнего и среднего полей стационарного (на крейсерском режиме полета) обтекания ЛА на адаптированных к их особенностям разностных сетках. Если указанные времена станут соизмеримы с временами счета дальнего поля существующими программными комплексами (типа “sBOOM”), замена их развитым выше также в рамках уравнений Эйлера ПКВ будет обусловлена не только его простотой, но и на порядки меньшими временами счета. Привлечение уравнений Навье–Стокса в их одномерном стационарном варианте при этом свелось к известному локальному решению, описывающему размазывание слабых УВ без влияния на их распространение от ЛА до Земли.

Авторы признательны К.С. Пьянкову и В.С. Горбовскому за полезные обсуждения.

Список литературы

  1. Ландау Л.Д. Об ударных волнах на далеких расстояниях от места их возникновения // Прикл. матем. и механ. 1945. Т. 9. Вып. 4. С. 286–292.

  2. Whitham G.B. The Propagation of Spherical Blast // Proceed. of the Royal Soc. Ser. A. 1950. V. 203. № 1075. P. 571–581.

  3. Whitham G.B. The Flow Pattern of a Supersonic Projectile // Comm. Pure Appl. Mathem. 1952. V. 5. № 3. P. 301–348.

  4. Христианович С.А. Ударная волна на значительном расстоянии от места взрыва // Прикл. матем. и механ. 1956. Т. 20. Вып. 5. С. 599–605.

  5. Губкин К.Е. Распространение разрывов в звуковых волнах // Прикл. матем. и механ. 1958. Т. 22. Вып. 4. С. 561–564.

  6. Walkden F. The Shock Pattern of a Wing-Body Combination Far from the Flight Path // Aeronautical Quarterly. 1958. V. 9. P. 164–194.

  7. Гриб А.А., Рыжов О.С., Христианович С.А. Теория коротких волн // Прикл. механ. и техн. физ. 1960. № 1. С. 63–74.

  8. Рыжов О.С. Затухание ударных волн в неоднородных средах // Прикл. механ. и техн. физ. 1961. № 2. С. 15–25.

  9. Полянский О.Ю. О затухании ударных волн в движущейся среде с переменными плотностью и температурой // Прикл. матем. и механ. 1960. Т. 24. Вып. 5. С. 912–915.

  10. Рыжов О.С. Затухание ударных волн в стационарных течениях // Прикл. механ. и техн. физ. 1961. № 6. С. 36–43.

  11. Carlson H.W. Influence of Airplane Configuration on Sonic Boom Characteristics // J. Aircraft. 1964. V. 1. № 2. P. 82–86.

  12. Carlson H.W., Mack R.J., Morris O.A. A Wind-Tunnel Investigation of the Effect of Body Shape on Sonic-Boom Pressure Distributions. NASA TN D-3106. 1965. 36 p.

  13. Seebass R. Sonic Boom Theory // J. Aircraft. 1969. V. 6. № 3. P. 177–184.

  14. Жилин Ю.Л. О звуковом ударе // Уч. зап. ЦАГИ. 1971. Т. 2. № 3. С. 1–11.

  15. Thomas Ch.L. Extrapolation of Wind-Tunnel Sonic Boom Signatures without Use of a Whitham F-Function. NASA SP-255. 1970. P. 205–217.

  16. Thomas Ch.L. Extrapolation of Sonic Boom Pressure Signatures by the Waveform Parameter Method. NASA TN D-6832. 1972. 35 p.

  17. Carlson H.W., Maglieri D.J. Review of Sonic Boom Generation Theory and Prediction Methods // J. Acoustical Soc. Am. 1972. V. 51. P. 675–682.

  18. Жилин Ю.Л. Теория звукового удара // Сборник работ по звуковому удару. Тр. ЦАГИ. 1973. Вып. 1489. С. 3–24.

  19. Уизем Дж.Линейные и нелинейные волны. М.: Мир, 1977. 622 с.

  20. Darden C.M. Minimization of Sonic-Boom Parameters in Real and Isothermal Atmospheres. NASA TN D-7842. 1975. 28 p.

  21. Darden C.M. Sonic Boom Theory: Its Status in Prediction and Minimization // J. Aircraft. 1977. V. 14. № 6. P. 569–576.

  22. Mack R.J., Darden C.M. Wind-Tunnel Investigation of the Validity of a Sonic-Boom-Minimization Concept. NASA Technical Paper 1421. 1979. 44 p.

  23. Siclari M.J., Darden C. M. An Euler Code Prediction of Near Field to Midfield Sonic Boom Pressure Signatures. AIAA-90-4000-CP. 1990. 14 p.

  24. Darden C.M. Limitation of Linear Theory for Sonic Boom Calculations // J. Aircraft. 1993. V. 30. № 3. P. 309–314.

  25. Жилин Ю.Л., Коваленко В.В. О связывании ближнего и дальнего полей в задаче о звуковом ударе // Уч. зап. ЦАГИ. 1998. Т. 29. №. 3–4. С. 111–122.

  26. Makino Y., Aoyama T., Iwamiya T. et al. Numerical Optimization of Fuselage Geometry to Modify Sonic-Boom Signature // J. Aircraft. 1999. V. 36. № 4. P. 668–674.

  27. Plotkin K.J. State of the Art of Sonic Boom Modeling // J. Acoustical Soc. Am. 2002. V. 111. № 1. Pt. 2. P. 530–536.

  28. Makino Y., Suzuki K., Noguchi M. et al. Nonaxisymmetrical Fuselage Shape Modification for Drag Reduction of Low-Sonic-Boom Airplane // AIAA J. 2003. V. 41. № 8. P. 1413–1420.

  29. Жилин Ю.Л., Ивантеева Л.Г., Чернышев С.Л. Эквивалентные тела вращения, обладающие минимальным уровнем звукового удара // Тр. ЦАГИ. 2005. Вып. 2670. С. 12–19.

  30. Жилин Ю.Л., Ивантеева Л.Г., Чернышев С.Л. О телах вращения минимального звукового удара // Уч. зап. ЦАГИ. 2006. Т. 37. № 3. С. 36–45.

  31. Коваленко В.В., Чернышев С.Л. К вопросу о снижении звукового удара // Уч. зап. ЦАГИ. 2006. Т. 37. № 3. С. 53–62.

  32. Farhat C., Maute K., Argrow B. et al. Shape Optimization Methodology for Reducing the Sonic Boom Initial Pressure Rise // AIAA J. 2007. V. 45. № 5. P. 1007–1018.

  33. Choi S., Alonso J., Kroo I. et al. Multifidelity Design Optimization of Low-Boom Supersonic Jets // J. Aircraft. 2008. V. 45. № 1. P. 106–118.

  34. Li W., Shields E., Le D. Interactive Inverse Design Optimization of Fuselage Shape for Low-Boom Supersonic Concepts // J. Aircraft. 2008. V. 45. № 4. P. 1381–1397.

  35. Bui Trong T. CFD Analysis of Nozzle Jet Plume Effects on Sonic Boom Signature, NASA dryden flight research center, California. 2009. 28 p.

  36. Кючул Чо. Анализ звукового удара с использованием полей возмущений, рассчитанных по нелинейной теории // Уч. зап. ЦАГИ. 2009. Т. 40. № 4. С. 45–55.

  37. Alauzet F., Loseille A. High-Order Sonic Boom Modeling Based on Adaptive Methods // J. Comput. Phys. 2010. V. 229. P. 561–593.

  38. Чернышев С.Л. Звуковой удар. М.: Наука, 2011. 351 с.

  39. Rallabhandi S. Advanced Sonic Boom Prediction Using Augmented Burgers Equation // J. Aircraft. 2011. V. 48. № 4. P. 1245–1253.

  40. Chiba K., Makino Y., Takatoya T. Design-Informatics Approach for Intimate Configuration of Silent Supersonic Technology Demonstrator // J. Aircraft. 2012. V. 49. № 5. P. 1200–1211.

  41. Ordaz I., Li W. Integration of Off-Track Sonic Boom Analysis in Conceptual Design of Supersonic Aircraft // J. Aircraft. 2014. V. 51. № 1. P. 23–28.

  42. Li W., Rallabhandi S. Inverse Design of Low-Boom Supersonic Concepts Using Reversed Equivalent-Area Targets // J. Aircraft. 2014. V. 51. № 1. P. 29–36.

  43. Minelli A., Salah el Din I., Carrier G. Inverse Design Approach for Low-Boom Supersonic Configurations // AIAA J. 2014. V. 52. № 10. P. 2198–2212.

  44. Ordaz I., Geiselhart K., Fenbert J. Conceptual Design of Low-Boom Aircraft with Flight Trim Requirement // J. Aircraft. 2015. V. 52. № 3. P. 932–939.

  45. Киселев А.Ф., Коваленко В.В., Притуло Т.М. Исследование звукового удара: расчет и эксперимент // Инженерн. журн.: наука и инновации. 2017. № 8. С. 1–11.

  46. Housman J., Kenway G., Jensen J. et al. Efficient Near-Field to Mid-Field Sonic Boom Propagation Using a High-Order Space Marching Method // AIAA Aviation 2019, Session: APA-25, Commercial Supersonic Technologies III, Dallas, Texas, June 20, 2019. 24 p.

  47. Браилко И.А., Крайко А.Н., Пьянков К.С. и др. Аэродинамические и акустические характеристики сверхзвуковой вентиляторной решетки с дозвуковой осевой компонентой вектора скорости // Аэромеханика и газовая динамика. 2003. № 4. С. 9–22.

  48. Efremov N.L., Kraiko A.N., Pyankov K.S. et al. Mathematical Simulation of Shock-Wave Structures Arising ahead of Fan Plane Cascades and Wheels // In “Turbomachines: Aeroelasticity, Aeroacoustics, and Unsteady Aerodynamics”. Moscow: Torus Press, 2006. P. 281–294.

  49. Ефремов Н.Л., Крайко А.Н., Пьянков К.С. и др. Ударно-волновые структуры перед неоднородной вентиляторной решеткой // Изв. РАН. Механика жидкости и газа. 2010. № 2. С. 135–152.

  50. Годунов С.К., Забродин А.В., Иванов М.Я. и др. Численное решение многомерных задач газовой динамики. М.: Наука, 1976. 400 с.

  51. Колган В.П. Применение принципа минимальных значений производной к построению конечно-разностных схем для расчета разрывных решений газовой динамики // Уч. зап. ЦАГИ. 1972. Т. 3. № 6. С. 68–77.

  52. Тилляева Н.И. Обобщение модифицированной схемы С.К. Годунова на произвольные нерегулярные сетки // Уч. зап. ЦАГИ. 1986. Т. 17. № 2. С. 18–26.

  53. Ганжело А.Н., Крайко А.Н., Макаров В.Е. и др. О повышении точности решения газодинамических задач // Кн. “Современные проблемы аэромеханики”. М.: Машиностр., 1987. С. 87–102.

  54. Родионов А.В. Повышение порядка аппроксимации схемы С.К. Годунова // Ж. вычисл. матем. и матем. физ. 1987. Т. 27. № 12. С. 1853–1860.

  55. Ландау Л.Д., Лифшиц Е.М. Теоретическая физика. Т. VI. Гидродинамика. Изд. третье, переработанное. М.: Наука, 1986. 733 с.

  56. Черный Г.Г. Газовая динамика. М.: Наука, 1988. 424 с.

  57. Овсянников Л.В. Лекции по основам газовой динамики. Москва–Ижевск: Ин-т комп. исслед., 2003. 336 с.

  58. Крайко А.Н. Теоретическая газовая динамика: классика и современность. М.: Торус пресс, 2010. 440 с.

  59. Валиев Х.Ф., Крайко А.Н., Тилляева Н.И. Устойчивость одномерных стационарных течений с детонационной волной в канале переменной площади // Ж. вычисл. матем. и матем. физ. 2020. Т. 60. № 4. С. 711–724.

  60. Горбовской В.С., Кажан А.В., Кажан В.Г. и др. О влиянии формы эпюры избыточного давления на громкость звукового удара // Уч. зап. ЦАГИ. 2020. Т. 51. № 2. С. 3–17.

Дополнительные материалы отсутствуют.

Инструменты

следующая статья выпуска предыдущая статья выпуска содержание выпуска

Журнал вычислительной математики и математической физики

Архивы выпусков Информация о журнале Отправить рукопись в журнал