Rapid Estimation of Supersonic Civil Aircraft Sonic Boom Characteristics in Standard Atmosphere Based on Analytical Solutions. Cruise Flight

Texto integral

ВВЕДЕНИЕ

Начиная с 50-х гг. XX века и по настоящее время в силу разработки сверхзвуковых пассажирских самолетов (СПС) первого и второго поколений актуальной является задача моделирования распространения волн звукового удара (ЗУ) в земной атмосфере. Для моделирования распространения волн ЗУ в стратифицированной атмосфере развиты модели на основе подходов нелинейной геометрической акустики (НГА) (cм., например, [1–3]), реализованные численно во многих программных кодах (см. обзоры [4, 5]). В более ранних и долго сохранявших свою актуальность программных кодах моделирование производилось в рамках модели идеальной среды. Наиболее современные модели учитывают и различные диссипативные факторы, связанные с неидеальностью среды.

Подходы высокочастотной акустики применимы в предположении, что распространяющиеся волны можно считать слабо модулированными, т. е. локально плоскими. Применительно к задаче о ЗУ это предположение может нарушаться главным образом при возникновении фокусировки акустических возмущений, происходящей по двум основным причинам. Во-первых, причиной фокусировки волн ЗУ служит мелкомасштабная атмосферная турбулентность. Для численного моделирования распространения волн ЗУ в турбулентном приземном слое атмосферы в настоящее время применяются эволюционные модели типа Хохлова-Заболотской [6–8], основанные на параболическом приближении теории дифракции, и более точные модели типа HOWARD (Heterogeneous One-Way Approximation for Resolution of Diffraction) [9–12]. Во-вторых, к фокусировке возмущений может приводить маневрирование летательного аппарата (ЛА). Наличие ускорения ЛА приводит к генерации сходящихся волн, фокусирующихся в атмосфере и на поверхности земли. Возникающие при этом интенсивности зависят от конкретных параметров полета ЛА, а амплитуда волн в результате фокусировки в атмосфере может в 2–5, а иногда и в 9–13 раз возрастать по сравнению с номинальными значениями [13]. В наиболее простых случаях маневрирование ЛА приводит, в терминах геометрической акустики, к возникновению простых (гладких) каустик [14]. Для описания акустического поля в окрестности простой каустики в рамках подхода сращиваемых асимптотических разложений в 1960-х гг. были получены модели LTE (Linear Tricomi Equation) [15] и NTE (Nonlinear Tricomi Equation) [16], реализованные с некоторыми модификациями в современных программных кодах [14, 17–19].

Для задач предварительного проектирования элементов СПС и оценок характеристик ЗУ в режиме реального времени целесообразно применение моделей распространения, основанных на приближении НГА идеальной среды, позволяющих получать точные решения, сводя тем самым к минимуму численные расчеты. В рамках подхода НГА идеальной движущейся плоско-слоистой среды задача о распространении акустических волн решена в квадратурах [20, 21]. Для моделирования прохождения волн через турбулентные флуктуации в приземном слое атмосферы может быть применена модель фазового экрана [22–24], для которой получено точное решение трехмерной задачи распространения в рамках НГА идеальной среды [25]. Приближенные оценки в рамках геометрической акустики [14] и ее модификаций [26] возможны также и для волн ЗУ вблизи каустик.

Текущая работа посвящена развитию упрощенных моделей распространения волн ЗУ в стандартной атмосфере, позволяющих проводить оценки характеристик распространяющихся волн ЗУ в режиме реального времени в отсутствие фокусировок. В разделе 1 приведено точное общее решение задачи на акустическое давление в рамках подхода НГА идеальной среды. В разделе 2 показано, что в общем случае задача о лучевой картине в атмосфере с кусочно-линейным профилем температуры может быть решена в элементарных функциях, а площадь лучевой трубки может быть вычислена алгебраически. Модельные профили параметров стандартной атмосферы и постановка начальных условий для задачи о ЗУ в рамках НГА приведены в разделах 3 и 4 соответственно. С помощью предложенного в данной работе подхода в разделе 5 проведен анализ зависимости геометрии распространения волн ЗУ в стандартной атмосфере от параметров крейсерского полета СПС. В условиях третьего семинара SBPW (Sonic Boom Prediction Workshop) 2020 произведен расчет эпюр избыточного давления на земле от демонстратора X-59, а также проведено сравнение с результатами численного моделирования [19].

1. ТОЧНОЕ РЕШЕНИЕ ДЛЯ АКУСТИЧЕСКОГО ДАВЛЕНИЯ

В приближении идеальной неподвижной среды акустическое давление моделируется следующим уравнением в лучевых координатах (см., например, [21]):

$\frac{\partial p}{\partial s}=\frac{1}{\mathcal{B}}\frac{\partial \mathcal{B}}{\partial s}p+\frac{\beta }{{\rho }_0{c}_0^3}p\frac{\partial p}{\partial \tau }$, (1)

где p — акустическое давление, s — переменная длина дуги геометроакустического луча, τ = t – φ — время в сопровождающей системе координат, φ — эйконал, ρ₀ и c₀ — плотность и скорость звука в среде соответственно, являющиеся в общем случае функциями пространственных координат, $\mathcal{B}=\sqrt{\frac{{\rho }_0{c}_0}{A}}$ — коэффициент инварианта Блохинцева [27] (если нелинейностью в задаче можно пренебречь, то вдоль луча сохраняется величина $\frac{p}{\mathcal{B}}$ — инвариант Блохинцева), A — площадь лучевой трубки, $\beta =\frac{\gamma +1}{2}$ — коэффициент нелинейности среды (γ — показатель адиабаты). Данная модель включает два основных механизма эволюции слабых ударных волн в атмосфере — геометрическое рассеяние и нелинейность. Аналогичные приближения долгое время применялись при оценке характеристик ЗУ на земле [4, 5, 28].

Известно, что уравнение (1) можно привести к уравнению простых волн с постоянными коэффициентами. Введем новую функцию  и новую эволюционную координату  соотношениями

$П=\frac{\overline{\mathcal{B}}}{\mathcal{B}}p, x\left(s\right)={\int }_{\overline{s}}^s\frac{{\overline{\rho }}_0{\overline{c}}_0^3}{{\rho }_0{c}_0^3}\frac{\mathcal{B}}{\overline{\mathcal{B}}}d{s}^{\text{'}}$. (2)

Здесь верхней чертой обозначены значения параметров на некотором расстоянии $\overline{s}>0$ от источника (ЛА), на котором для уравнения (1) ставится начальное условие. Введенная функция П сохраняется вдоль луча, если отбросить нелинейные эффекты. Подставляя замену (2) в уравнение (1), приходим к уравнению простых волн с постоянными коэффициентами

$\frac{\partial П}{\partial x}=\frac{\beta }{{\overline{\rho }}_0{\overline{c}}_0^3}П\frac{\partial П}{\partial \tau }$, (3)

для которого известно неявное точное решение. Пусть начальное условие для уравнения (1) имеет вид

${\left.p\right|}_{s=\overline{s}}=p_0\left(\xi ,\tau \right)$, (4)

где ξ = (ξ₁, ξ₂) — лучевые координаты, дополняющие координату s. Тогда для функции П имеем соответствующее условие ${\left.П\right|}_{x=0}=p_0\left(\xi ,\tau \right)$, и неявное решение (3) запишется в виде (координаты ξ далее опускаются)

$П\left(x, \tau \right)=p_0\left(\eta \right), \eta =\tau +\frac{\beta p_0\left(\eta \right)}{{\overline{\rho }}_0{\overline{c}}_0^3}x$. (5)

Тогда решение задачи (1), (4) в исходных переменных принимает вид:

$p\left(s, \tau \right)=\frac{\overline{\mathcal{B}}}{\mathcal{B}}{p}_0\left(\eta \right), \eta =\tau +\frac{\beta }{\overline{\mathcal{B}}}{p}_0\left(\eta \right){\int }_{\overline{s}}^s\frac{\mathcal{B}}{{\rho }_0{c}_0^3}d{s}^{\text{'}}$. (6)

Непосредственное применение данной формулы для вычислений является затруднительным по следующим причинам:

а) в неоднородной среде входящий в решение интеграл и коэффициент A(s, ξ) требуют численного расчета с учетом зависимости площади лучевой трубки  и характеристик среды от пространственных координат;

б) решение в общем случае является неоднозначным и требует обработки возникающих «перехлестов» интегральной кривой путем определения положения разрывов и их расстановки [2].

В следующем разделе показано, что в рамках модели стандартной атмосферы площадь лучевой трубки A может быть вычислена точно в силу кусочно-линейной структуры профиля температуры, что существенно упрощает анализ. В данном разделе остановимся подробнее на второй проблеме — неоднозначности решений (5) и (6).

Рассмотрим функцию

$Ф={\int }_{-\infty }^{\tau }П\left(x, {\tau }^{\text{'}}\right)d{\tau }^{\text{'}}$,

часто называемую в литературе потенциалом или псевдо-потенциалом в силу своего вида [9, 10, 28]. Она обладает следующим замечательным свойством [29]: решение задачи в терминах Ф всегда непрерывно и может быть непосредственно найдено интегрированием выражения (5) в предположении гладкости начальной функции p₀(τ):

$\begin{array}{l}Ф\left(x, \tau \right)={\int }_{-\infty }^{\tau }П\left(x, {\tau }^{\text{'}}\right)d{\tau }^{\text{'}}=\\ ={Ф}_0\left(\eta \right)-\frac{\beta x}{2{\overline{\rho }}_0{\overline{c}}_0^3}{p}_0^2\left(\eta \right). \end{array}$ (7)

Здесь введено обозначение для начального потенциала ${Ф}_0\left(\eta \right)={\int }_{-\infty }^{\eta }{p}_0\left({\eta }^{\text{'}}\right)d{\eta }^{\text{'}}$. Решение (7) остается справедливым и после возникновения разрывов в решении для П. Оно, так же, как и решение (5), является неоднозначным, но позволяет легко отобрать физически допустимое решение в соответствии с принципом максимума: при любом фиксированном значении x₀ > 0 физически допустимая функция Ф(x₀, τ) есть огибающая максимумов многозначного решения (7). Отслеживание переходов между ветвями неоднозначного решения в процессе построения допустимого решения для Ф позволяет автоматически определить положение разрывов в соответствующем решении для П. Алгоритм построения физически допустимого решения для искомой функции с помощью решения для потенциала, применяемый в данной работе, описан в статьях [29, 30].

Аналогичным образом в том же приближении может быть получено решение для акустического давления при наличии ветра в атмосфере [21].

2. ЛУЧЕВАЯ КАРТИНА В СТАНДАРТНОЙ АТМОСФЕРЕ

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

Уравнения лучей с длиной дуги луча s в качестве эволюционной координаты в неоднородной неподвижной среде имеют вид (см., например, [31]):

$\frac{dr}{ds}=c_{0}q, \frac{dq}{ds}=-\frac{\nabla c_0}{c_0^2}$. (8)

Здесь q = ∇φ — градиент эйконала. В плоско-слоистой атмосфере, где все параметры являются функциями геометрической высоты z, задача на шесть уравнений (8) может быть сведена к системе из двух уравнений и решена в квадратурах путем выбора геометрической высоты z в качестве эволюционной координаты с учетом того, что величина q удовлетворяет уравнению эйконала

${q_x}^2+{q_y}^2+{q_z}^2=\frac{1}{c_0^2}$.

В силу плоской слоистости среды c₀ = c₀(z), и $\nabla c_0=k\frac{d{c}_0}{dz}$, где k — орт вдоль оси Oz. Следовательно, в силу уравнений системы (8) проекции q_x и q_y вектора q на горизонтальные оси земной системы координат сохраняются вдоль луча, т.е. являются функциями только координат ξ, а вертикальная проекция q_z вычисляется из уравнения эйконала:

$q_z\left(z, \xi \right)=\pm \sqrt{\frac{1}{c_0^2\left(z\right)}-Q^2\left(\xi \right)}, Q= \sqrt{{q_x}^2+{q_y}^2}$. (9)

Таким образом, второе векторное уравнение системы (8) исключается, а так как в качестве независимой переменной выбрана координата z, первое векторное уравнение сводится к системе из двух уравнений:

$\frac{dx}{dz}=\frac{q_x}{q_z}, \frac{dy}{dz}=\frac{q_y}{q_z}$. (10)

Решение уравнений (10) имеет форму

$x=x\left(z, \xi \right), y=y\left(z, \xi \right)$. (11)

Известно, что в стандартной атмосфере на высотах ниже 25 км профиль температуры является монотонной кусочно-линейной функцией высоты. Этот факт дает ряд преимуществ при решении задачи о лучевой картине применительно к проблеме первичного ЗУ (т.е. из рассмотрения исключены лучи, отражающиеся от верхних слоев атмосферы). Во-первых, монотонность профиля температуры исключает возможность возникновения точек разворота первичного ЗУ, т.е. выбранная параметризация лучей геометрической высотой  может считаться допустимой (в смысле строгой монотонности z = z(s)). Во-вторых, кусочно-линейная структура профиля температуры позволяет вычислить в элементарных функциях все интегралы, входящие в решение. Фиксируя знак величины q_z < 0 (т.е. выбирая знак минус в решении (9)) и подставляя выражение (9) для q_z в систему (10), получаем общее решение:

$\left\{\begin{array}{c}x\left(z, \xi \right)=x_0\left(\xi \right)-q_x\left(\xi \right)I_0\left(z, \xi \right),\\ y\left(z, \xi \right)=y_0\left(\xi \right)-q_y\left(\xi \right)I_0\left(z, \xi \right),\end{array}\right.$ (12)

где

$I_0\left(z, \xi , \xi \right)={\int }_{z_0\left(\xi \right)}^z\frac{d{z}^{\text{'}}}{\sqrt{\frac{1}{c_0^2\left(z^{\text{'}}\right)}-Q^2\left(\xi \right)}}$, (13)

а отображение $r_0\left(\xi \right)=\left(x_0\left(\xi \right), y_0\left(\xi \right),z_0\left(\xi \right)\right)$ задает множество начальных точек лучей (начальную поверхность). При кусочно-линейном профиле температуры профиль квадрата скорости звука также кусочно-линеен. Пусть на некотором линейном участке профиля квадрат скорости звука задается в виде $c_0^2\left(z\right)=a+bz$. На данном участке профиля неопределенный интеграл F₀, соответствующий (13), может быть вычислен с помощью замены переменных $\zeta \left(z\right)=\sqrt{1-Q^2\left(a+bz\right)}$. В итоге получаем

$\begin{array}{c}{F}_0\left(\zeta \right)=-\frac{2}{Q^{3}b}\int \sqrt{1-{\zeta }^2}d\zeta =\\ =-\frac{1}{Q^{3}b}\left(\zeta \sqrt{1-{\zeta }^2}+\arcsin \left(\zeta \right)\right)+C. \end{array}$ (14)

Здесь и далее в похожих выражениях C — константа интегрирования, своя в каждом конкретном случае.

В лучевых координатах (z, ξ), в отличие от координат (s, ξ), площадь лучевой трубки A'(s, ξ) вырождается, и соответствующий определитель понижает свою размерность (частные производные обозначены нижними индексами):

$\begin{array}{c}{A}^{\text{'}}\left(z, \xi \right)=\left|\begin{array}{ccc}{x}_z& x_{{\xi }_1}& x_{{\xi }_2}\\ y_z& y_{{\xi }_1}& y_{{\xi }_2}\\ z_z& z_{{\xi }_1}& z_{{\xi }_2}\end{array}\right|=\\ =\left|\begin{array}{ccc}{x}_z& x_{{\xi }_1}& x_{{\xi }_2}\\ y_z& y_{{\xi }_1}& y_{{\xi }_2}\\ 1& 0& 0\end{array}\right|=\\ =x_{{\xi }_1}{y}_{{\xi }_2}-x_{{\xi }_2}{y}_{{\xi }_1}. \end{array}$ (15)

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

$\begin{array}{l}A\left(s, \xi \right)=-\left|\begin{array}{ccc}{x}_s& x_{{\xi }_1}& x_{{\xi }_2}\\ y_s& y_{{\xi }_1}& y_{{\xi }_2}\\ z_s& z_{{\xi }_1}& z_{{\xi }_2}\end{array}\right|=\\ =-\frac{\partial z}{\partial s}\left(s, \xi \right)A^{\text{'}}\left(z\left(s, \xi \right), \xi \right), \end{array}$ (16)

где z(s, ξ) — решение уравнения для координаты z системы (8), представляющее собой формулу для замены параметризации лучей. В силу соотношения (9) для ее частной производной справедливо соотношение

$\frac{\partial z}{\partial s}=-\sqrt{1-c_0^2\left(z\right)Q^2\left(\xi \right)}$. (17)

Начальное условие для (17), согласованное с решением (12), (13), можно записать в виде

$z\left(\mathrm{0,}\xi \right)=z_0\left(\xi \right)$. (18)

Решение задачи (17), (18) имеет вид

$s=-{\int }_{z_0\left(\xi \right)}^{z\left(s, \xi \right)}\frac{d{z}^{\text{'}}}{\sqrt{1-c_0^2\left(z^{\text{'}}\right)Q^2\left(\xi \right)}}$. (19)

При $c_0^2\left(z\right)=a+bz$ для соответствующего неопределенного интеграла в терминах введенной ранее переменной ζ получаем

$F_s\left(\zeta \right)=-\frac{2}{b{Q}^2} \zeta +C$. (20)

Частные производные решения (12), (13), входящие в (15), вычисляются с использованием формулы Лейбница для производной интеграла, зависящего от параметра. Применительно к задаче о ЗУ множество начальных точек лучей вырождается в кривую и представляет собой траекторию полета ЛА. Далее будем полагать, что r₀(ξ) = r₀(ξ₁), и обозначим для удобства производную $\frac{d{r}_0}{d{\xi }_1}\equiv {\dot{r}}_0$. Тогда для площади лучевой трубки  можно получить соотношение:

$$\begin{gathered} A^{\text{'}}\left(z, \xi \right)=\frac{\partial q_x}{\partial {\xi }_2}\left\{k_{xx}\left({\dot{y}}_0-{\dot{z}}_0{\left.\frac{q_y}{q_z}\right|}_{z_0}\right)-k_{xy}\left({\dot{x}}_0-{\dot{z}}_0{\left.\frac{q_x}{q_z}\right|}_{z_0}\right)\right\}+ \\ +\frac{\partial q_y}{\partial {\xi }_2}\left\{k_{xy}\left({\dot{y}}_0-{\dot{z}}_0{\left.\frac{q_y}{q_z}\right|}_{z_0}\right)-k_{yy}\left({\dot{x}}_0-{\dot{z}}_0{\left.\frac{q_x}{q_z}\right|}_{z_0}\right)\right\}+ \\ +\left(\frac{\partial q_x}{\partial {\xi }_1}\frac{\partial q_y}{\partial {\xi }_2}-\frac{\partial q_x}{\partial {\xi }_2}\frac{\partial q_y}{\partial {\xi }_1}\right)\left(k_{xx}{k}_{yy}-k_{xy}^2\right), \end{gathered}$$ (21)

где введены обозначения

$k_{xx}=I_2-q_y^2{I}_1, k_{xy}=q_x{q}_y{I}_1, k_{yy}=I_2-q_x^2{I}_1$,

$\begin{array}{l}{I}_1={\int }_{z_0\left(\xi \right)}^z\frac{d{z}^{\text{'}}}{{\left(\frac{1}{c_0^2\left(z^{\text{'}}\right)}-Q^2\left(\xi \right)\right)}^{3/2}}, \\ I_2={\int }_{z_0\left(\xi \right)}^z\frac{d{z}^{\text{'}}}{c_0^2\left(z^{\text{'}}\right){\left(\frac{1}{c_0^2\left(z^{\text{'}}\right)}-Q^2\left(\xi \right)\right)}^{3/2}}. \end{array}$ (22)

Первообразные для неопределенных интегралов, соответствующих I₁ и I₂, при $c_0^2\left(z\right)=a+bz$ имеют вид:

$$\begin{gathered} F_1\left(\zeta \right)=\frac{2}{Q^{5}b}\left\{\frac{\sqrt{1-{\zeta }^2}}{2\zeta }\left({\zeta }^2+2\right)+\frac{3}{2}\arcsin \left(\zeta \right)\right\}+C, \\ {F}_2\left(\zeta \right)=\frac{2}{Q^{3}b}\left\{\frac{\sqrt{1-{\zeta }^2}}{\zeta }+\arcsin \left(\zeta \right)\right\}+C. \end{gathered}$$ (23)

Таким образом, в атмосфере с кусочно-линейным профилем температуры, в частности, в стандартной атмосфере, лучи (формулы (12), (13), (14)) и площадь лучевой трубки (формулы (21), (22), (23)) первичного ЗУ могут быть рассчитаны точно в элементарных функциях. Для применения полученных первообразных к вычислению интегралов необходимо на каждом из линейных кусков профиля температуры задавать соответствующие значения параметров a и b профиля квадрата скорости звука, после чего интегралы могут быть вычислены в виде сумм интегралов по линейным участкам. Отметим, что полученные первообразные обладают особенностью при b = 0. В этом случае все первообразные представляют собой линейные функции высоты z. Первообразные (23) для вычисления площади лучевой трубки также обладают особенностью в точке ζ = 0. Физический смысл данной особенности заключается в возникновении точек разворота лучей (см., например, [31]). Разворот луча проявляется в смене знака величины q_z, связанной с переменной ζ соотношением ζ = –c₀q_z. Условие q_z = 0 нарушает допустимость выбранной параметризации. Как сказано ранее, в силу монотонности профиля температуры в стандартной атмосфере лучи первичного ЗУ не будут обладать точками разворота, а значит, при рассмотрении возмущений первичного ЗУ для всех лучей, достигающих земной поверхности, замена параметра z = z(s) будет допустимой. Те же лучи, которые имеют точку разворота, не могут достичь земли, не отразившись от верхних слоев атмосферы. Следовательно, они относятся к возмущениям вторичного ЗУ, имеющим существенно меньшую интенсивность на земле и потому представляющим меньший интерес. Разворот луча, сопровождающийся его уходом в верхние слои атмосферы, называют эффектом отсечки.

Полученные решения (16), (17), (21)–(23) могут быть применены для вычисления эволюционной координаты x(s) (2) в уравнении для акустического давления. Для этого осуществим переход к переменной z под интегралом, являющейся более естественной в силу плоской слоистости атмосферы. Тогда интеграл принимает вид:

$x\left(z\right)=-\frac{{\overline{\rho }}_0{\overline{c}}_0^3}{\overline{\mathcal{B}}}{\int }_{\overline{z}}^z\frac{1}{\sqrt{{\rho }_0{c}_0^5{\left(1-c_0^2{Q}^2\right)}^{3/2}}} \frac{1}{\sqrt{A^{\text{'}}}}d{z}^{\text{'}}$. (24)

Здесь $\overline{z}=z\left(\overline{s},\xi \right)$. Для вычисления интеграла (24) необходимо задать высотные профили скорости звука и плотности. В следующем разделе описано задание данных профилей для стандартной атмосферы.

3. ПАРАМЕТРЫ СТАНДАРТНОЙ АТМОСФЕРЫ

Высотный профиль температуры в стандартной атмосфере до высоты z^** = 25000 м может быть задан следующим образом [32] (отличием геометрической высоты от геопотенциальной до характерных высот полета ЛА будем пренебрегать):

$T_0\left(z\right)=\left\{\begin{array}{c}{T}_{з}-{\beta }_{T}z, 0\le z\le z^{*},\\ T^{*}, z^{*}<z\le z^{**}.\end{array}\right.$

Здесь Т_з = 288.15 К, β_Т = 0.0065 К/м, z^* = 11000 м, Т^* = 216.65 К, z^** = 25000 м. Соответствующий профиль квадрата скорости звука имеет вид:

$c_0^2\left(z\right)=\left\{\begin{array}{c}\gamma R_{в}\left(T_{з}-{\beta }_{T}z\right), 0\le z\le z^{*},\\ \gamma R_{в}{T}^{*}, z^{*}\le z\le z^{**}.\end{array}\right.$ (25)

Здесь R_в = 287.05287 Дж/(кг К) — газовая постоянная воздуха, γ = 1.4 — показатель адиабаты.

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

$\frac{d{p}_0}{dz}=-{\rho }_{0}g\left(z\right)$,

где g(z) — ускорение свободного падения с учетом зависимости от высоты z. С учетом уравнения состояния идеального газа и высотного профиля температуры высотный профиль давления может быть вычислен в терминах геопотенциальной высоты

$H\left(z\right)=\frac{1}{g_0}{\int }_0^{z}g\left(z^{\text{'}}\right)d{z}^{\text{'}}$,

где g₀ = g(0) = 9.80665 м/с². Пренебрегая отличием геопотенциальной высоты от геометрической, профиль давления зададим следующими соотношениями:

$p_0\left(z\right)=\left\{\begin{array}{c}{p}_{з}{\left(1-\frac{\beta z}{T_{з}}\right)}^{\frac{g_0}{{\beta }_T{R}_{в}}}, 0\le z\le z^{*},\\ p^{*}{e}^{-\frac{g_0}{R_{в}{T}^{*}}\left(z-z^{*}\right)}, z^{*}<z\le z^{**}.\end{array}\right.$

Здесь p_з = 101325 Па, $p^{*}=p_{з}{\left(1-\frac{\beta z^{*}}{T_{з}}\right)}^{\frac{g_0}{{\beta }_T{R}_{в}}}$ = 70152.24 Па. Тогда из уравнения состояния p₀ = ρ₀R_вT₀ получаем профиль плотности:

${\rho }_0\left(z\right)=\left\{\begin{array}{c}{\rho }_{з}{\left(1-\frac{\beta z}{T_{з}}\right)}^{\frac{g_0}{{\beta }_{T}Rв}-1}, 0\le z\le z^{*},\\ {\rho }^{*}{e}^{-\frac{g_0}{R_{в}{T}^{*}}\left(z-z^{*}\right)}, z^{*}\le z\le z^{**}.\end{array}\right.$. (26)

Здесь введены обозначения

${\rho }_{з}=\frac{p_{з}}{R_{в}{T}_{з}}, {\rho }^{*}=\frac{p^{*}}{R_{в}{T}^{*}}$.

Таким образом, подынтегральная функция в выражении (24) в условиях плоско-слоистой среды с кусочно-линейным профилем температуры является композицией элементарных функций. Тем не менее, чтобы избежать громоздких выкладок, его вычисление в данной работе проводилось численным путем.

4. НАЧАЛЬНЫЕ УСЛОВИЯ НА ТРАЕКТОРИИ ПОЛЕТА

При решении задачи о ЗУ в качестве лучевых координат ξ₁ и ξ₂ часто выбирают t_a и φ, где t_a — физическое время, параметризующее траекторию ЛА, φ — азимутальный угол выхода луча из начального конуса Маха (см. рис. 1). Траектория ЛА, представляющая собой множество начальных точек лучей, задается как параметрическая кривая

$r=r_a\left(t_a\right)=\left(x_a\left(t_a\right),y_a\left(t_a\right), z_a\left(t_a\right)\right)$.

 

Рис. 1. (а) — Начальные условия на конусе Маха и (б) — угловые положения вектора воздушной скорости

 

Начальные условия для уравнений (10) запишутся в виде

${\left.x\right|}_{z=z_a\left(t_a\right)}=x_a\left(t_a\right),{\left.y\right|}_{z=z_a\left(t_a\right)}=y_a\left(t_a\right)$. (27)

Значения параметров q_x и q_y на лучах определяются ориентацией осей и углами раствора конусов Маха в точках траектории ЛА, зависящими от числа Маха полета, углов курса и наклона траектории ЛА, а также атмосферного ветра. Ось конуса Маха направлена вдоль вектора V воздушной скорости ЛА, направление которого в земных координатах (x, y, z) задается углами θ₁ и θ₂ (рис. 1б). Вектор воздушной скорости ЛА в данной точке траектории отличается от касательного к траектории вектора земной скорости V_з на вектор скорости ветра в данной точке v₀. В рамках текущей работы атмосферный ветер считается нулевым, и вектора V и V_з совпадают. В рассматриваемой точке траектории введены локальные координаты (x', y', z'), такие что ось x' совпадает по направлению с вектором воздушной скорости, а оставшиеся координатные оси, дополняющие систему до правой тройки, выбраны таким образом, чтобы плоскость (x', z') была параллельна оси z земной системы координат. Азимутальный угол φ вводится как полярный угол в плоскости (y', z') (см. рис. 1а). Начальные направления лучей определяются нормалями к конусам Маха и задаются в виде семейства

$q_0=q_0\left(t_a,\phi \right)$.

Во введенных вспомогательных координатах (x', y', z') в отсутствие ветра компоненты q₀ определяются углом Маха μ и азимутальным углом φ:

$\begin{array}{c}\left({q_0}_{x\text{'}},{q_0}_{y\text{'}},{q_0}_{z\text{'}}\right)=\\ =\frac{1}{c_0\left(z_a\left(t_a\right)\right)}\left(\sin \mu \left(t_a\right),\cos \mu \left(t_a\right)\sin \phi , -\cos \mu \left(t_a\right)\cos \phi \right). \end{array}$ (28)

Для перехода в земные координаты (x, y, z) необходимо совершить два поворота на углы θ₁ и θ₂.

Начальные условия для акустического давления (1) формируются при помощи инвариантов преобразования Галилея, осуществляемого при переходе от системы координат, связанной с ЛА, к земной системе координат [20]. Модели типа (1), основанные на приближении нелинейной геометрической акустики (НГА), справедливы там, где выполнено предположение о независимости распространения возмущений вдоль лучей. Поэтому начальные условия задаются на некотором удалении от ЛА. Предположим, что поле избыточного давления известно в некоторой области с цилиндрической боковой поверхностью с образующими, параллельными оси x', т.е. направлению воздушной скорости ЛА (рис. 2). Избыточное давление на боковых поверхностях, соответствующих различным точкам траектории, задается функцией

$p_0=p_0\left(t_a,\phi ,l\right)$. (29)

 

Рис. 2. К постановке начальных условий для акустического давления

 

Здесь величина l играет роль фазовой переменной. В геометрической акустике в качестве фазовой переменной выступает время в сопровождающей системе τ = t – φ. При фиксированных пространственных координатах τ изменяется с физическим временем t. Чтобы связать две фазовые переменные l и τ, рассмотрим третью фазовую переменную: расстояние N по нормали к волновым фронтам, измеряемое относительно фронта, выходящего из точки O. Данная величина является галилеевым инвариантом при переходе от системы координат, связанной с ЛА, к земной системе.

С одной стороны, в связанной с ЛА системе координат имеем связь

$N=l\sin \mu$.

В земной системе через фиксированную точку в пространстве возмущение проходит расстояние N по нормали к фронту со скоростью c₀(z_a(t_a)) (на масштабах рассматриваемой области атмосфера полагается однородной) за время t = τ:

$N=c_a\tau$.

Таким образом, фазовые переменные l и τ связаны соотношением

$\tau =\frac{l\sin \mu }{c_a}$. (30)

Пусть радиус цилиндрической боковой поверхности равен R. В локальной системе координат, в которой атмосфера покоится, распространение волны до цилиндрической границы происходит за время T, определяемое соотношением

$c_{a}T=\frac{R}{\cos \mu }$.

За то же время в земной системе при отсутствии ветра луч проходит расстояние

$\overline{s}=c_{a}T$.

Тогда для начального значения длины дуги луча $\overline{s}$ имеем:

$\overline{s}\left(t_a\right)=\frac{R}{\cos \mu \left(t_a\right)}$, (31)

а для соответствующего значения $\overline{z}$ справедлива формула

$\overline{z}\left(t_a,\phi \right)=z_a\left(t_a\right)+c_0\left(z_a\left(t_a\right)\right){q_0}_z\left(t_a,\phi \right)\overline{s}\left(t_a\right)$. (32)

Здесь q_0z(t_a, φ) — -компонента начального вектора q₀ в земной системе координат.

Таким образом, начальное условие для уравнения (1) можно записать в виде

${\left.p\right|}_{s=\frac{R}{\cos \mu \left(t_a\right)}}=p_0\left(t_a,\phi ,\frac{\tau c_0\left(z_a\left(t_a\right)\right)}{\sin \mu \left(t_a\right)}\right)=F\left(t_a,\phi ,\tau \right)$. (33)

5. ЗВУКОВОЙ УДАР В КРЕЙСЕРСКОМ РЕЖИМЕ ПОЛЕТА ЛЕТАТЕЛЬНОГО АППАРАТА

5.1. Геометрия распространения волн звукового удара

С использованием точного решения (12), (13), (14) для лучей в стандартной атмосфере можно построить геометрию распространения возмущений ЗУ от сверхзвукового ЛА на различных режимах его полета. Лучи, выходящие из фиксированной точки траектории, образуют обратный коноид Маха. Семейство обратных коноидов Маха на траектории полета ортогонально семейству волновых фронтов — прямых коноидов Маха. Поскольку оба семейства заполняют одну и ту же область пространства, геометрию распространения волн можно характеризовать как прямыми, так и обратными коноидами. При наличии полученного точного решения для лучей распространение волн ЗУ удобнее характеризовать геометрией обратных коноидов.

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

Будем предполагать, что крейсерский полет ЛА происходит в направлении оси x в плоскости y = 0 на высотах $z_a=const\in \left[z^{*},z^{**}\right]$ (величины z^* и z^** определены в разделе 3). Тогда решение (12), (13) с начальными условиями (27), (28) примет вид:

$\begin{array}{l}x\left(z,t_a,\phi \right)=c_{a}M{t}_a-\frac{\sin \mu }{c_a}{I}_0\left(z,\phi \right), \\ y\left(z, t_a,\phi \right)=\frac{\cos \mu \sin \phi }{c_a}{I}_0\left(z, \phi \right), \end{array}$ (34)

$I_0\left(z,\phi \right)=\left\{\begin{array}{l}\frac{c_a}{\cos \mu \cos \phi }\left(z^{*}-z_a\right)+F_0\left(\zeta \left(z\right)\right)-\\ -F_0\left(\zeta \left(z^{*}\right)\right), 0\le z\le z^{*};\\ \frac{c_a}{\cos \mu \cos \phi }\left(z-z_a\right), z^{*}<z\le z_a.\end{array}\right.$ (35)

Здесь I₀(z, φ) — интеграл, заданный формулой (13), функция F₀(ζ) определяется формулой (14), c_a = c₀(z_a), а z^* = 11000 м. В текущей постановке число Маха полета M и скорость звука на высоте полета c_a постоянны, поэтому I₀ не зависит от параметра t_a. Величина t_a определяет лишь координату вершины коноида Маха и потому не является существенной при изучении его геометрии, сохраняющейся во время полета ЛА. Для определенности положим t_a = 0, фиксируя тем самым точку траектории полета (0, 0, z_a) и соответствующий обратный коноид Маха, имеющий вершину в этой точке.

Строго говоря, величина z в решении (34), (35) ограничена не только условием $z\in \left[\mathrm{0,}{z}_a\right]$. Полученное решение также необходимо дополнить условием допустимости параметризации лучей

$\begin{array}{l}\frac{1}{c_0^2\left(z\right)}-Q^2\left(\phi \right)=\frac{1}{c_0^2\left(z\right)}-\\ -\frac{1}{c_a^2}\left({\sin }^2\mu +{\cos }^2\mu {\sin }^2\phi \right)>\mathrm{0,} \end{array}$ (36)

вытекающим из предположения q_z < 0, где величина q_z определяется отрицательной ветвью формулы (9). При z^* < z < z_a, как можно убедиться, обратившись к решению (34), (35), лучи представляют собой прямые линии, поэтому в данном диапазоне высот отсечки не происходит, и величина q_z сохраняет свое значение вдоль луча. Условие (36) для этого случая сводится к неравенству cos²φ > 0, запрещающему лучи, горизонтально выходящие из начального конуса Маха. При 0 ≤ z < z^* с учетом зависимости (25) для квадрата скорости звука в стандартной атмосфере выражение (36) можно переписать в следующем виде:

$z>\frac{1}{\gamma R_{в}{\beta }_T}\left(\frac{c_a}{\sqrt{{\sin }^2\mu +{\cos }^2\mu {\sin }^2\phi }}-c_g\right)$. (37)

Здесь c_g = c₀(0) — скорость звука на поверхности земли. Полученное неравенство следует воспринимать следующим образом: при положительной правой части неравенства величина его правой части дает высоту точки разворота луча z_CO, а при отрицательной правой части луч достигает поверхности земли, и минимальное значение z на нем равно нулю.

Построение обратного коноида Маха производится путем трассировки лучей, выходящих из рассматриваемой точки траектории под различными азимутальными углами φ в соответствии с решением (34), (35) при t_a = 0 с учетом ограничения (37). На рис. 3а изображены обратные коноиды Маха, построенные указанным способом для случая крейсерского полета ЛА на высоте z_a = 15 км с числами Маха M = 1.05, 1.15, 1.5, 3 в положительном направлении оси x земной системы координат. При построении каждого коноида производилась трассировка 281 лучей для значений φ, равномерно распределенных на отрезке [–70°, 70°]. Видим, что для коноидов, соответствующих числам Маха 1.05 и 1.15, все лучи дошли до точки поворота. Для чисел Маха, меньших некоторого порогового значения, волны первичного ЗУ не достигают поверхности земли.

 

Рис. 3. (а) — Обратные коноиды Маха и (б) — линии пересечения обратных коноидов с поверхностью земли в крейсерском полете ЛА на высоте 15 км при различных числах Маха

 

Линия пересечения обратного коноида Маха с поверхностью земли может быть построена как параметрическая кривая x = x(φ), y = y(φ) на основе формул (34), (35) при z = 0 и t_a = 0, т.е. как координатная линия на обратном коноиде. Значения параметра φ на линии пересечения ограничены неравенством (37) при z = 0. Для этого случая неравенство запишется в виде

${\sin }^2\phi <\frac{{\left(\frac{c_a}{c_g}\right)}^2-{\sin }^2\mu }{{\cos }^2\mu }$. (38)

При положительной правой части предельный случай равенства в (38) определяет азимутальный угол отсечки φ_CO > 0, и на линии пересечения |φ| < φ_CO. Заметим, что данное ограничение при выбранных высотах полета не зависит от высоты полета и определяется только величиной числа Маха. В случае отрицательной правой части обратный коноид Маха не пересекает поверхность земли ни при каких значениях φ, т.е. происходит глобальная отсечка волн первичного ЗУ градиентом температуры. Знак правой части определяется числом Маха полета $M=\frac{1}{\sin \mu }$. При выбранных параметрах полета в соответствии с профилем квадрата скорости звука в стандартной атмосфере (25) граничное число Маха отсечки M_CO есть величина

$M_{CO}=\frac{c_g}{c_a}\approx 1.153$. (39)

Семейство линий пересечения для случая крейсерского полета ЛА на высоте z_a = 15 км для ста значений чисел Маха, равномерно распределенных на отрезке [1.05, 3], построено на рис. 3б. Видно, что при изменении числа Маха геометрия линии пересечения существенно изменяется, а ширина возмущенной области на поверхности земли монотонно увеличивается при увеличении числа Маха. Картина линий пересечения с увеличением числа Маха уплотняется, а начиная с некоторого значения числа Маха вовсе перестает существенно меняться, т.е. геометрия линии пересечения «стабилизируется». При M → ∞ азимутальный угол отсечки φ_CO стремится к предельному значению

${{\phi }_{CO}}_{\max }=\arcsin \left(\frac{c_a}{c_g}\right)\approx 60.1°$,

не зависящему от высоты полета ЛА.

Линии пересечения обратных коноидов Маха с поверхностью земли также изменяют свою геометрию при изменении высоты полета ЛА. На рис. 4а построено семейство линий пересечения для ста значений z_a, равномерно распределенных на отрезке [11, 20] км при числе Маха полета M = 1.5. Видим, что при увеличении высоты полета ЛА кривизна линии пересечения уменьшается, а ширина увеличивается. Высокая плотность расположения кривых свидетельствует о более плавной зависимости положения и формы линии пересечения от высоты полета, чем от числа Маха.

 

Рис. 4. (а) — Линии пересечения обратных коноидов Маха с поверхностью земли в крейсерском полете с числом Маха M = 1.5 и (б) — семейство зависимостей ширины возмущенной области от числа Маха при различных высотах крейсерского полета

 

Определим ширину L возмущенной области на земле как удвоенную координату y луча на земле, удовлетворяющего условию отсечки ζ(z = 0) = 0. В силу отсутствия особенностей при ζ = 0 у первообразной (14) для величины L при M ≥ M_CO можно получить:

$L\left(M,z_a\right)=2\frac{\sqrt{1-{\left(\frac{M_{CO}}{M}\right)}^2}}{\sqrt{M_{CO}^2-1}}\left(z_a+z^{*}\frac{\arcsin \left(\sqrt{1-\frac{1}{M_{CO}^2}}\right)}{\sqrt{1-\frac{1}{M_{CO}^2}}}\right)$. (40)

При $M\in \left[\mathrm{1,}{M}_{CO}\right]$ необходимо положить L = 0. Видим, что ширина возмущенной области на земле является линейной функцией высоты полета ЛА. Формула (40) также свидетельствует об эффекте «стабилизации» величины L при возрастании числа Маха полета. На рис. 4б изображено семейство зависимостей ширины L возмущенной области на земле от числа Маха при различных высотах крейсерского полета ЛА. С ростом числа Маха при переходе через данное пороговое значение величина L начинает быстро нарастать, и превышает 20 км уже при числе Маха M = 1.2 для всех значений высоты полета. При этом влияние изменения высоты полета несущественно при числах Маха, близких к M_CO, и существенно возрастает при увеличении числа Маха до тех пор, пока увеличение числа Маха не перестанет ощутимо изменять геометрию благодаря эффекту «стабилизации» и высота полета не станет единственным существенно влияющим на величину L фактором. Таким образом, ширина области, подверженной возмущениям ЗУ на земле, монотонно возрастает и с числом Маха, и с высотой полета ЛА. При выбранных диапазонах изменения числа Маха и высоты полета величина L может превосходить 100 км.

5.2. Пример расчета эпюр избыточного давления на земле

В качестве примера применения полученных решений произведен расчет распространения волн ЗУ от демонстратора X-59 в условиях третьего семинара SBPW (Sonic Boom Prediction Workshop), организованного NASA в 2020 г. [33]. X-59 — экспериментальный сверхзвуковой самолет, разрабатываемый в рамках миссии NASA QueSST, призванный продемонстрировать низкий уровень ЗУ на земле (первый полет запланирован на 2024 г.). Участникам семинара представлялись начальные данные для расчета распространения волн ЗУ, полученные с помощью методов вычислительной газовой динамики. Для случая крейсерского полета демонстратора в стандартной атмосфере задано число Маха полета M = 1.4, высота полета z_a = 16459.2 м, высота уровня земли z_g = 110.011 м, а также начальные условия для расчета акустического давления на земле в форме (29) при радиусе боковой границы расчетной области течения вблизи демонстратора R = 82.296 м. Отражение от земной поверхности предлагается учитывать простым домножением результата распространения на коэффициент r = 1.9.

На рис. 5 демонстрируется работа алгоритма извлечения однозначного разрывного решения для функции П из неоднозначного решения (5) с помощью принципа максимума для потенциала Ф в случае φ = 0°. Видим, что выбор максимального значения потенциала для каждого значения τ дает непрерывную функцию, огибающую многозначное решение. Отслеживание переходов между однозначными «ветвями» при построении данной огибающей и совершение соответствующих переходов в решении для функции П позволяет автоматически расставить разрывы, удовлетворяющие правилу равенства площадей [34].

 

Рис. 5. (а) — Решения для функции П и (б) — потенциала Ф в случае φ = 0°

 

Для заданных выше параметров крейсерского полета демонстратора азимутальный угол отсечки на земле φ_СО равен 44.84°. Расчеты эпюр избыточного давления на земле проведены для 81 значения азимутального угла φ от φ_min = –40° до φ_max = 40° с шагом 2°. Начальные условия в форме (33) и соответствующие результаты расчетов на земле приведены на рис. 6 (в силу симметрии задачи на рисунке изображены только эпюры для неотрицательных значений φ). Можно увидеть, что при увеличении азимутального угла происходит упрощение эпюры давления на земле в результате слияния скачков, чему способствует увеличение расстояния, проходимого возмущением до достижения поверхности земли. Тем не менее, формы волн остаются многоскачковыми и качественно довольно далекими от N-образной. Это обеспечивает большее время нарастания избыточного давления на головном и хвостовом ударных фронтах и меньшую громкость волн ЗУ по сравнению с идеальной N-волной той же амплитуды [35].

 

Рис. 6. (а) — Начальные эпюры избыточного давления и (б) — соответствующие эпюры на земле для различных азимутальных углов выхода лучей

 

В силу важности учета конечности времени нарастания в скачках избыточного давления полученные результаты в представленном виде не позволяют рассчитать громкость с помощью современных метрик из-за наличия в эпюре математических разрывов. Несмотря на то, что амплитуда (максимальное в эпюре значение) в полученных решениях убывает строго монотонно с ростом азимутального угла (рис. 6б), расчет громкости для аналогичных результатов с учетом диссипативных факторов, проведенный в работе [19] с помощью известного программного кода sBOOM (NASA), показывает немонотонное поведение ее значения. Для корректной оценки громкости при расчетах распространения волн ЗУ в приближении идеальной среды ранее применялись эмпирические поправки (см., например, [28]), основанные на результатах статистического анализа экспериментальных данных. В настоящее время применяются модели, учитывающие физические процессы диссипации волн в атмосфере (см., например, обзор [5]). Несмотря на отмеченный недостаток, присущий всем подходам, игнорирующим диссипативные процессы в волнах ЗУ, упрощенные модели воспроизводят основную структуру эпюр избыточного давления на земле, что позволяет судить об успешности проектирования геометрии СПС, основываясь на расположении и интенсивности скачков в волне ЗУ, проводя расчеты в режиме реального времени.

Влияние диссипативных факторов на форму волны наглядно продемонстрировано на рис. 7, где результаты расчетов на основе предложенного подхода сравниваются с результатами, полученными численно в работе [19] на основе современной модели «дополненного уравнения Бюргерса» (augmented Burgers equation) [36, 37] для трех значений азимутального угла φ. Отличие данной модели от рассматриваемой в текущей работе состоит в учете диссипативных факторов — термовязкого затухания и молекулярной релаксации. Видим, что положение «скачков» в эпюре с учетом диссипации изменяется слабо в условиях текущей задачи, однако форма волны существенно сглаживается с понижением абсолютного значения избыточного давления и появлением конечных времен нарастания в «скачках», что для малых скачков приводит к их неразличимости и формированию единого плавно нарастающего ударного фронта.

 

Рис. 7. Эпюры избыточного давления на земле в приближении идеальной среды (сплошные линии) и с учетом диссипативных факторов (пунктирные линии) [19] для азимутальных углов (а) — φ = 0°, (б) — φ = 20° и (в) — φ = 40°

 

ЗАКЛЮЧЕНИЕ

С помощью полученного в данной работе решения задачи о лучевой картине в стандартной атмосфере исследована геометрия распространения волн ЗУ для случая сверхзвукового крейсерского полета СПС. Определена зависимость геометрии возмущенной области на земле от высоты и числа Маха полета. Для ширины возмущенной области на земле получена простая расчетная формула.

Для расчета эпюры избыточного давления на земле предложен подход, основанный на модели НГА идеальной среды, позволяющий проводить расчеты акустического давления на земле с произвольным начальным профилем волны в режиме реального времени путем введения псевдо-потенциала, обеспечивающего автоматическое выполнение правила равенства площадей. В условиях семинара SBPW3 2020 для демонстратора X-59 проведен расчет акустического давления в волнах ЗУ на земле; путем сравнения результатов расчета с результатами, полученными на основе численного решения дополненного уравнения Бюргерса с помощью программного кода sBOOM (NASA), проиллюстрировано, что учет диссипативных факторов не приводит к изменению основной структуры эпюр избыточного давления на земле.

Полученное решение задачи о ЗУ в стандартной атмосфере, помимо быстрых оценок характеристик ЗУ в крейсерском режиме полета СПС, может быть применено для исследования геометрии распространения волн ЗУ в различных режимах полета СПС, в том числе, для определения условий возникновения фокусировки волн на земле при сверхзвуковом маневрировании. Оценка характеристик ЗУ на земле при сверхзвуковом маневрировании ЛА является задачей будущих работ.

Статья подготовлена в рамках реализации Программы создания и развития научного центра мирового уровня «Сверхзвук» на 2020-2025 гг. при финансовой поддержке Минобрнауки России (соглашение от 17 мая 2022 г. № 075-15-2022-1023).

Sobre autores

A. Korunov

Central Aerohydrodynamic Institute

Autor responsável pela correspondência
Email: korunov.ao@phystech.edu
Rússia, Zhukovsky

V. Gusev

Central Aerohydrodynamic Institute; Lomonosov Moscow State University

Email: vgusev@bk.ru
Rússia, Zhukovsky; Moscow

V. Gorbovskoy

Central Aerohydrodynamic Institute

Email: korunov.ao@phystech.edu
Rússia, Zhukovsky

Bibliografia

Arquivos suplementares