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

Полный текст

Введение

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

Сценарий, по которому развивается ЛТП пограничного слоя, существенно зависит от уровня возмущений в набегающем потоке [1]. В сверхзвуковых аэродинамических трубах (АДТ) имеются акустические пульсации в рабочей части [2–4], обусловленные шумом турбулентного пограничного слоя на стенках сопла и рабочей части, а также различными неровностями на них. На гладких моделях ЛТП может сопровождаться появлением турбулентных клиньев [5], которые возникают в некоторых местах над обтекаемой поверхностью и не меняют своего расположения от пуска к пуску АДТ.

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

С N-волнами экспериментаторы связывают аномально высокие уровни пульсаций массового расхода (до 20% от среднего удельного массового расхода в набегающем потоке) в пограничном слое плоского треугольного крыла с притупленной передней кромкой [6]. Дело в том, что взаимодействие N-волны и пограничного слоя вблизи притупления приводит к появлению около линии растекания возмущений, которые распространяются в пограничном слое вниз по потоку. Наиболее энергоемкой является низкочастотная (до 3 кГц) часть наблюдаемых пульсаций. Аналогичные эксперименты [7] по возбуждению N-волной высокоинтенсивных возмущений в пограничном слое на модели плоского дельта-крыла показали, что максимальная величина пульсаций составляет 12–15% и слабо зависит от условий обтекания. Различные аспекты взаимодействия N-волн и сверхзвукового пограничного слоя рассмотрены в последовавших экспериментальных работах [8–10].

Численному моделированию взаимодействия N-волн и сверхзвукового пограничного слоя на заостренной пластине посвящены работы [11–13]. В частности, в труде [11] разработана оригинальная методика моделирования тонкой неровности на стенке рабочей части АДТ.

Наша статья является развитием работ [11–13] и направлена на исследование случая пластины с цилиндрически притупленной передней кромкой, обтекаемой воздушным потоком при числе Маха М = 2.5. Численное исследование проводится в рамках полных уравнений Навье–Стокса, а стационарные решения получаются методом установления по времени. Результаты численного моделирования сопоставляются с данными аэродинамического эксперимента [8].

Следует отметить, что в отличие от случая заостренной пластины расчет сверхзвукового обтекания затупленной пластины имеет ряд вычислительных особенностей, связанных с взаимодействием N-волны с головной ударной волной перед притуплением. Такие особенности обсуждаются в нашей работе.

Постановка задачи и численное моделирование

Основные условия задачи и параметры потока соответствуют экспериментальному исследованию [8], которое выполнено в малотурбулентной сверхзвуковой АДТ Т-325 Института теоретической и прикладной механики Сибирского отделения Российской академии наук. Общая конфигурация задачи схематично представлена на рис. 1. Рассматривается течение воздуха в рамках модели совершенного двухатомного газа.

 

Рис. 1. Схема (а) взаимодействия N-волны, возбуждаемой двумерной неровностью на боковой стенке аэродинамической трубы (вид сверху) и возмущение давления p′ от N-волны (б): M_inf, R_e1,inf, T_inf – число Маха, единичное число Рейнольдса, статическая температура набегающего потока соответственно.

 

В работе [8] N-волна возбуждалась на стенке рабочей части с помощью изоленты, которая наклеивалась на расстоянии x^* = -233.5 мм вверх по потоку от передней кромки пластины и формировала малую двухмерную неровность (координата x^* = 0 соответствует передней кромке пластины; звездочка обозначает размерную величину).

При численном моделировании N-волны использована эквивалентная модель [11], в которой тонкая двухмерная прямоугольная неровность заменяется параболической дужкой $z_1=\epsilon F\left(x_1\right)$ (рис. 2):

$\begin{array}{c}F\left(x_1\right)=4x_1(1-x_1),\\ x_1=(x^{*}-x_0^{*})/a^{*},\\ z_1=(z^{*}-z_{\min }^{*})/a^{*},\\ \epsilon =h^{*}\text{}/a^{*},\end{array}$

где h^* – толщина прямоугольной неровности; a^* – ее ширина.

 

Рис. 2. Тонкая лента (а) и ее эквивалентная математическая модель (б).

 

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

$\Delta p=\frac{p-p_{\infty }}{\rho u_{\infty }^2}=\epsilon \frac{F^{\text{'}}\left(x_1\right)}{\sqrt{M_{\infty }^2-1}}+{\rm O}\left({\epsilon }^2\right),$$\Delta u=\frac{u-u_{\infty }}{u_{\infty }}=-\epsilon \frac{F^{\text{'}}\left(x_1\right)}{\sqrt{M_{\infty }^2-1}}+{\rm O}\left({\epsilon }^2\right),$

$\Delta w=\frac{w}{u_{\infty }}=\epsilon F^{\text{'}}\left(x_1\right)+{\rm O}\left({\epsilon }^2\right),$$\Delta T=\frac{T-T_{\infty }}{T_{\infty }}=\epsilon \frac{\left(\gamma -1\right)M_{\infty }^2{F}^{\text{'}}\left(x_1\right)}{\sqrt{M_{\infty }^2-1}}+{\rm O}\left({\epsilon }^2\right).$

В нашей работе рассматривается N-волна с амплитудой 1%, которая соответствует геометрическим параметрам ленты a^* = 20 мм, h^* = 0.2 мм, x₀^* = -251 мм.

Численный метод

Для решения численной задачи использован оригинальный пакет программ [14], с помощью которого проводится численное моделирование на основе полных уравнений Навье–Стокса. Уравнения Навье–Стокса интегрируются в безразмерных переменных. Декартовы координаты x^* = xL, y^* = yL, z^* = zL отнесены к характерному масштабу длины L, время t ^* = tL/V_∞ – к характерному масштабу времени L/V_∞, компоненты вектора скорости u ^* = u/V_∞, v ^* = v/V_∞, w ^* = w/V_∞ – к модулю вектора скорости набегающего потока V_∞, давление p^* = p(p_∞V_∞²) – к удвоенному скоростному напору набегающего потока, остальные газодинамические переменные – к их значениям в набегающем потоке. Звездочка в верхнем индексе обозначает размерную величину; если звездочка отсутствует, переменная предполагается безразмерной. Символ $\infty$ соответствует значению переменной в набегающем потоке. В нашей работе L = 1 м.

Динамический коэффициент молекулярной вязкости вычисляется по формуле Сазерленда $\mu =(1+T_{\mu })/ (T+T_{\mu })T^{3/2},$ где $T_{\mu }=T_{\mu }^{*}/T_{\infty }^{*}=110.4 K / 103.6 K\approx \mathrm{1.07.}$ Газ является совершенным с постоянными показателем адиабаты g = 1.4 и числом Прандтля $\mathrm{Pr}={\mu }^{*}{c}_p^{*}/{\lambda }^{*}=\mathrm{0.72.}$ Для замыкания задачи использовано уравнение состояния совершенного газа $\gamma M_{\infty }^{2}p=\rho T.$

Набегающий поток считается ламинарным, а фоновые пульсации массового расхода, которые в эксперименте [8] достигали величины 0.4%, отсутствуют.

Расчетная сетка и граничные условия

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

 

Рис. 3. Расчетная подобласть 1 для пластины с острой (а) и с затупленной (б) передней кромкой

 

На рис. 3 также построены сечения с полями продольной компоненты вектора скорости и одна из сеточных плоскостей. Левые границы рис. 3а и 3б соответствуют стенке АДТ; на них накладывается условие прилипания потока и теплоизолированной поверхности, а также условие возбуждения N-волны, описанное выше. N-волна распространяется вниз по потоку. Из приведенных данных можно заключить, что максимальное относительное возмущение продольной компоненты вектора скорости, вызванное неровностью, составляет около 1.4%.

Боковые границы подобласти 1 не визуализированы на рис. 3; на ближней к читателю границе накладывается условие первого рода, соответствующее набегающему вдоль оси x газовому потоку. На другой границе накладывается мягкое граничное условие в виде экстраполяции всех зависимых газодинамических переменных задачи изнутри расчетной области вдоль сеточных линий. На верхней и нижней границах по оси y накладываются условия симметрии.

Установившееся течение с N-волной в подобласти 1 используется для формулирования стационарного граничного условия в подобласти 2. Для этого все зависимые переменные с выходной (правой) границы подобласти 1 накладываются в виде условия первого рода во входном (левом) сечении подобласти 2. При этом выходное сечение подобласти 1 и входное сечение подобласти 2 совпадают узел в узел. Расчетная сетка в области 2 показана на рис. 4.

 

Рис. 4. Расчетная подобласть 2 и фрагменты расчетной сетки: (а) – сечение xy, затупленная кромка; (б) – сечение xy, острая кромка; (в) – сетка в сечении xy, затупленная кромка; (г) – сетка в сечении yz, острая кромка.

 

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

В подобласти 2 на поверхности пластины накладываются условия прилипания и теплоизолированной поверхности. На выходной границе подобласти 2 задается мягкое граничное условие. На границах в боковом направлении, а также на границе x < 0, y = 0 накладывается условие симметрии.

Размер ячейки расчетной сетки по направлению x показан на рис. 4а, по направлению z – на рис. 4г. По нормальному к пластине направлению y сетка сгущена к стенке. По опыту расчетных исследований, мелкость сетки в пограничном слое оказывается избыточной для получаемых установившихся ламинарных полей течения, в том числе, возмущенных падающей слабой ударной волной.

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

Особенность моделирования взаимодействия N-волны и ударной волны

При взаимодействии N-волны с головной ударной волной, возникающей перед притуплением, появляются вычислительные трудности, которые могут приводить к численной неустойчивости процесса расчета. Причины этого связаны с особенностями моделирования скачков уплотнения с применением схем сквозного счета [15, 16]. Одним из способов избежать численную неустойчивость такого рода является выравнивание сеточных плоскостей вдоль скачка и измельчение расчетной сетки в его окрестности в целях разрешения профиля ударной волны. Выбор оптимального значения сеточного разрешения на ударной волне учитывает выводы работы [16] и опирается на тестовые расчеты, проведенные в рамках настоящей работы.

На рис. 5 показана структура ударной волны вдоль линии симметрии s₁ (рис. 6а), полученная при решении уравнений Навье–Стокса. Толщина ударной волны, обусловленная конечными значениями вязкости и теплопроводности, составляет примерно 4 мкм. Данные рис. 5 свидетельствуют о сходимости расчетных данных на сетках с различным сгущением.

 

Рис. 5. Распределение давления поперек ударной волны в сечении s1 при y = 0.

 

Чтобы изучить особенности полученных решений, на различных сетках рассматриваются распределения вдоль линий (x = const, y = const), где точка (х^*, у^*) = (–0.0002, 0) находится на линии торможения (s1), а (х^*, у^*) = (0.00245, 0.0026) – в пограничном слое вблизи сопряжения притупления с плоской поверхностью. Разрешение ударной волны анализируется в зависимости от шага расчетной сетки Dx поперек нее.

Из рис. 6б, на которой проиллюстрировано поле давление в сечении y = 0 и расчетная сетка в данном сечении, следует, что ударная волна попадает в область повышенного сеточного разрешения и не выходит за ее пределы. Из рис. 7–9, на которых представлены распределения давления, температуры и скорости вдоль линий s₁ и s₂ в зависимости от Dx поперек ударной волны, следует, что в области ударной волны s₁ осцилляции существенно сильнее, чем в пограничном слое s₂.

По мере уменьшения величины Dx расчет становится более устойчивым, а распределения характеристик вдоль линий s₁ и s₂ остаются вблизи некоторого среднего распределения с погрешностью, не превосходящей 10^–3. Это наблюдение соответствует выводам других исследователей о том, что схемы сквозного счета корректно воспроизводят параметры потока до и после скачка, а проблемы с устойчивостью счета связаны с внутренней структурой скачка и исчезают по мере разрешения вязкой структуры скачка [16].

По результатам данного анализа выбран оптимальный размер ячейки, который зафиксирован для последующих расчетов: Dx = 6 × 10^-7.

Результаты

Визуализация течения

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

 

 

Рис. 6. Поле продольной компоненты вектора скорости в сечении xy (а) и поле давления в сечении xz при y = 0 (Dx* = 8 × 10–7м) (б).

 

Рис. 7. Давление вдоль линии s1 (а) и линии s2 (б).

 

Рис. 8. Температура вдоль линии s1 (а) и линии s2 (б).

 

Рис. 9. Продольная компонента скорости вдоль линии s1 (а) и линии s2 (б).

 

По рис. 10а и 11а видно, что влияние N-волны на пограничный слой на заостренной пластине проявляется в виде слабого утолщения пограничного слоя в области переднего фронта N-волны при z ≈ 0.13, который приводит к появлению дефицита скорости в пограничном слое. Аналогично наблюдается слабое утончение пограничного слоя в области заднего фронта N-волны при z ≈ 0, которое приводит к появлению избытка скорости в пограничном слое.

В случае притупленной передней кромки (рис. 10б и 11б) влияние N-волны проявляется иначе на качественном уровне. В частности, наблюдаются две изолированные друг от друга локальные области утолщения и утончения пограничного слоя при z ≈ –0.03 и при z ≈ 0.16, причем изменение толщины пограничного слоя выражено заметно сильнее.

 

Рис. 10. Продольная компонента скорости в сечении x* = 0.06 м: острая (а) и затупленная (б) пластина.

 

Рис. 11. Продольная компонента скорости в сечении x* = 0.2 м: острая (а) и затупленная (б) пластина.

 

Рис. 12. Поле продольной компоненты скорости в сечении y* = 0.0029 м: острая (а) и затупленная (б) пластина.

 

На рис. 12 показаны соответствующие поля продольной компоненты скорости показаны на виде сверху в сечении y ^* = 0.0029 м. За острой кромкой образуется широкий след из пары вихрей (в соответствии с исходной терминологией [6–10]), а за затупленной кромкой образуются два изолированные друг от друга следа меньшего поперечного размера.

Более наглядно качественное отличие между случаями острой и притупленной кромки можно проследить на рис. 13–14, где в различных сечениях показаны поля возмущения продольной составляющей вектора скорости относительно течения без N-волны. Как отмечено ранее, за острой передней кромкой образуется след, соответствующий дефициту и избытку скорости, соответственно, на переднем и заднем фронтах N-волны. В этом случае след имеет большой размер, но амплитуда возмущений в нем относительно мала.

В случае затупленной пластины образуются два таких следа: один соответствует переднему фронту N-волны, второй – заднему. Эти следы имеют меньший размер, но амплитуда соответствующих возмущений примерно втрое выше по сравнению со случаем острой передней кромки.

Сравнение с экспериментальными данными

Для валидации расчетных результатов задачи взаимодействия N-волны с затупленной пластиной проведено сравнение среднего значения массового расхода и среднеквадратичного значения пульсаций массового расхода с экспериментальными данными из работы [8].

Анализ расчетных данных указывает на то, что нестационарность, связанная с N-волной, является слабой. Поэтому при обтекании затупленной пластины взаимодействие N-волны с пограничным слоем можно в главном приближении считать стационарным. В этом случае справедлива процедура численного анализа среднеквадратичного значения пульсаций массового расхода, приведенная в [13]. Так как пульсации N-волны наблюдаются в диапазоне малых частот, предположим, что их можно считать квазистационарными. В этом приближении N-волна и возбуждаемые ею возмущения осциллируют в направлении оси z как целое. Тогда перемещение таких пульсаций есть

$\Delta {m^{\text{'}}}_{N-wave}=f^{\text{'}}\left(z\right)\Delta z\left(t\right),$

 

Рис. 13. Возмущение продольной компоненты скорости в сечении y* = 0.0033 м: острая (а) и затупленная (б) пластина.

 

Рис. 14. Возмущение продольной компоненты скорости: острая (а) и затупленная (б) пластина при разных x *.

 

где f (z) – распределение стационарного массового расхода.

Пусть фоновые пульсации массового расхода в потоке без N-волны равны $\Delta m{\text{'}}_F\left(t\right)$ и не зависят от z. Тогда суммарные пульсации (дрожащая N-волна плюс фон) имеют вид

$\Delta m^{\text{'}}=\Delta {m^{\text{'}}}_{N-wave}+\Delta {m^{\text{'}}}_F.$

Рис. 15. Средние значения массового расхода: (а) – х* = 60 мм, у* ≈ 2.983 мм, (б) – х* = 90 мм, у* ≈ 3.079 мм.

 

Рис. 16. Среднеквадратичные пульсации массового расхода: (а) – х* = 60 мм, у* ≈ 2.983 мм, (б) х* = 90 мм, у* ≈ 3.079 мм.

 

Их среднеквадратичные значения равны

$⟨\Delta m^{\text{'}}⟩={\left[f^{\text{'}}\left(z\right)\right]}^2⟨{\left(\Delta z\right)}^2⟩+⟨{\left(\Delta {m^{\text{'}}}_F\right)}^2⟩.$

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

$f^{\text{'}}\left(z\right)=\frac{3^{*}f\left(z_{j+1}\right)-4^{*}f\left(z_j\right)+f\left(z_{j-1}\right)}{z_{j+1}-z_{j-1}}.$

В результате среднеквадратичные пульсации массового расхода имеют вид

$⟨m^{\text{'}}⟩=\sqrt{\alpha {\left[f\text{'}\left(z\right)\right]}^2+\beta },$

где $\beta =⟨{\left(\Delta {m^{\text{'}}}_F\right)}^2⟩$ определяется по экспериментальным данным и характеризует фоновые пульсации массового расхода, а $\alpha =⟨{\left(\Delta z\right)}^2⟩$ выбирается так, чтобы расчетное распределение совпадало с экспериментальным в точке локального максимума, соответствующего переднему фронту N-волны.

На рис. 15 средние значения массового расхода вдоль линии х^* = 60 мм и х^* = 90 мм сопоставлены с данными эксперимента. Соответствующее сопоставление для среднеквадратичного значения пульсаций массового расхода приведено на рис. 16. Из сопоставления следует, что результаты численного моделирования удовлетворительно согласуются с данными эксперимента. Это свидетельствует о корректности численного моделирования взаимодействия N-волны и затупленной пластины, обтекаемой сверхзвуковым потоком газа.

Обсуждение

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

Каждый из двух следов, по-видимому, будет располагаться там, где за острой кромкой в поперечном направлении располагались линии переднего и заднего фронтов N-волны. Анализ влияния N-волны на пограничный слой в зависимости от величины притупления требует трудоемкого параметрического исследования, которое не было проведено в данной работе, но несомненно представляет интерес.

Местное изменение характеристик пограничного слоя из-за его взаимодействия с падающей N-волной оказывает влияние на устойчивость пограничного слоя и процесс его турбулизации. В работе [13] это было продемонстрировано для случая заостренной пластины. На стационарное течение в пограничном слое с падающей N-волной и без нее накладывались периодические возмущения, которые вводились в пограничный слой с помощью генератора типа вдув–отсос, размещаемого на поверхности пластины. Спектральные и геометрические характеристики генератора подбирались на основе расчета линейной устойчивости пограничного слоя.

Было показано, что в области дефицита скорости, которая формируется в следе от переднего фронта N-волны, линия ЛТП искажается, смещаясь вверх по потоку и формируя подобие турбулентного клина. Аналогично в области избытка скорости, которая формируется в следе от заднего фронта N-волны, линия ЛТП смещается вниз по потоку.

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

В работе [17] были получены первые качественные расчетные данные о влиянии N-волны на ЛТП на затупленной пластине. Пограничный слой возбуждался идентично случаю острой пластины [13].

В частности, использовались характеристики генератора возмущений, определенные из расчета устойчивости пограничного слоя на острой, а не на затупленной пластине. В связи с этим линия ЛТП в отсутствие N-волны оказалась ниже по потоку по сравнению с ее положением на острой пластине.

Однако появление N-волны привело к существенно более значительному относительному смещению положения ЛТП вверх по потоку на обоих следах от N-волны, чем смещение на линии переднего фронта от N-волны в случае заостренной пластины.

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

Выводы

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

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

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

Финансирование

Работа выполнена в МФТИ при финансовой поддержке РНФ по проекту № 21-19-00307 (И. В. Егоров, Н. К. Нгуен – постановка задачи, модельные граничные условия, численное моделирование) и по проекту № 19-19-00470 (П. В. Чувахов – анализ численных результатов; особенности моделирования ударной волны п. 2.3).

Об авторах

И. В. Егоров

Московский физико-технический институт (национальный исследовательский университет); Центральный аэрогидродинамический институт им. проф. Н. Е. Жуковского

Email: pavel_chuvahov@mail.ru
Россия, Долгопрудный; Жуковский

Н. К. Нгуен

Московский физико-технический институт (национальный исследовательский университет)

Email: pavel_chuvahov@mail.ru
Россия, Долгопрудный

П. В. Чувахов

Московский физико-технический институт (национальный исследовательский университет)

Автор, ответственный за переписку.
Email: pavel_chuvahov@mail.ru
Россия, Долгопрудный

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

Дополнительные файлы