Numerical modeling of a wing leading-edge thermal regimes for a reusable space vehicle

Cover Page

Abstract


Throughout the history of human exploration of outer space, work is underway to reduce the cost of bringing cargo into space. One of the technically feasible solutions to achieve this is the use of smallsized reusable aerospace vehicles. As the new thermal protection materials are developed, they are employed for the construction of the reusable aerospace vehicles (RSV). In this paper, the assessment is given of the possibility of making RSV wing leading edge from an Al2O3 fiber based heat-resistant porous ceramic. The main advantages of using such material are its relatively low values of thermal conductivity and density, which makes it possible to improve weight characteristics of the RSV. The material of the support structure is heat-resistant carbon fiber reinforced polymer (CFRP). Due to the porous nature of such thermal protection system (TPS), it is necessary to consider the effect of air pressure on the thermal conductivity of the material. Therefore, a computational mathematical model is proposed that allows one to take into account this dependence of thermal conductivity on temperature and pressure for the wing edge porous TPS of an aerospace vehicle, during its re-entry in the atmosphere. Based on the temperature field inside the leading edge, the minimum thickness of the thermal protection coating was determined so that the support structure temperature stays within its maximum permissible operating range. It is shown that the Al2O3 heat-resistant porous ceramic can provide the required thermal protection, so that the maximum temperature of the composite support structure does not exceed 250 °C on the entire re-entry flight path.


Введение На протяжении всей истории освоения людьми космического пространства, ведутся работы по снижению стоимости вывода грузов в космос. Одним из технически реализуемых решений для достижения этой цели новых направлений в освоении околоземного космического пространства является использование малоразмерных многоразовых аэрокосмических аппаратов (МКА) «крылатого» типа [1; 2]. Для обеспечения тепловой защиты конструкции данных аппаратов применяются современные термостойкие композитные покрытия. Так, беспилотный аэрокосмический аппарат X-37 (рис. 1), имеет теплозащитную конструкцию передней кромки крыла, частью которой служит материал AETB-8 (Alumina Enhanced Thermal Barrier), представляющий собой высокопористую керамику на основе волокон оксида алюминия (Al2O3), кремния и алюминий-боросиликатных волокон [3; 4]. Теплозащитное покрытие (ТЗП) аппарата способно выдерживать температуры до 1970 К [5]. Рис. 1. Многразовый космический аппарат X-37 [6] [Fig. 1. X-37 reusable space vehicle [6]] Во многих других космических кораблях также используют композитные материалы в качестве ТЗП (табл. 1). Таблица 1 Перечень материалов передних кромок малоразмерных аппаратов «крылатого» типа [List of wing leading-edge materials for “winged” small-sized space vehicles] Название космического аппарата Год запуска Страна Особенности конструкции передней кромки крыла БОР 4 [7; 8] 1982-1984 СССР ПКТ-ФЛ (уносимая) из фенол-формальдегидной ткани, пропитанной смесью фенол-формальдегидных смол; между верхней и нижней поверхностями крыльев также находился материал типа фетр, пропитанный специальным составом на основе воды. Испарение воды обеспечивало эффективное охлаждение во время интенсивного нагрева. Характерный размер кромки: 70-80 мм. Максимальная рабочая температура: до 1700 °С Бор 5 [7; 8] 1984-1988 СССР молибденовый сплав и дополнительное антиокислительное покрытие. К кромке из молибденового сплава крепился ТЗП на основе кварцевого волокна и хром-алюминий-фосфатного связующего. Характерный размер кромки: » 100 мм. Максимальная рабочая температура: до 2000 °С X-33 [8-10] Отсутствует США, НАСА и др. Углерод-углеродный материал X-34 [11] Отсутствует США Многоразовая теплозащитная плитка на основе Кремниевой керамики X-38 [12; 13] Отсутствует NASA, ESA Карбидокремниевая керамика Hermes [8; 14] Отсутствует Европейские страны, Россия, Канада и др. Углеродный носок (С/SiC) - высокопрочный углеродный материал типа «Карбоксил» в виде оболочки толщиной 3 мм, подкрепленной поперечными ребрами толщиной 6 мм Skylon [15; 16] В разработке ESA Углерод - углеродный материал. Композитная внутренняя силовая конструкция крыльев В ряде перспективных аппаратов «крылатого» типа, созданных для совершения орбитальных полетов на низкой околоземной орбите (НОО), также следует отметить аппараты FTB-1 (Italian Aerospace Research Centre) и DreamChaser (Sierra Nevada Corporation) (рис. 2). а б Рис. 2. Космические аппараты: а - FTB-1; б - DreamChaser [17; 18] [Fig. 2. Space vehicles: а - FTB-1; б - DreamChaser [17; 18]] На первом аппарате, для кромки малого радиуса (R = 0,04 м), планируется использование высокотемпературной керамики (UHTC) на основе ZrB2 или SiC [19-21]. Однако ТЗП аппарата DreamChaser, ввиду большего радиуса кромки крыла, также, как и на X-37, выполнено из материала AETB-8 [22]. Как было отмечено, этот материал имеет пористую структуру, а, следовательно, меньшую плотность, что позволяет значительно улучшить весовые характеристики ТЗП по сравнению с UHTC. Материалы и допущения Приведенные сведения доказывают высокий интерес космической индустрии к композитам, а именно к пористой керамике как к современному теплозащитному материалу передней кромки крыла МКА. Цель данного исследования состоит в оценке возможности выполнения передней кромки крыла МКА из термостойкой пористой керамики на основе волокон оксида алюминия (Al2O3) [23] (табл. 2). Основные характеристики пористой керамики на основе Al2O3 [List of properties for porous Al2O3 based ceramics] Таблица 2 Материал Структура Плотность, кг/м3 Теплопроводность, Вт/м·К Теплоемкость, Дж/кг·К ТЗМК1700 [23] Пористая (термические характеристики зависят от давления) Матрица: кремнийорганическое связующее Волокна: Al2O3 250 Давление 10-5 Па: 0,07 (при 20 °С) - 0,43 (при 1700 0С) Давление 105 Па (1 атм.): 0,1 (при 20 °С) - 0,5 (при 1700 °С) 680 (при 20 °С) - 1260 (при 1700 °С) Представлена математическая модель, с помощью которой проводился анализ для наиболее теплонагруженного участка траектории - входа космического аппарата в атмосферу. Основным критерием для оценки ТЗП была максимальная температура силовой конструкции, выполненная из термостойкого углепластика. На протяжении всей траектории полета эта температура не должна превышать значений в пределах 250-300 °С [24]. В связи с пористой структурой материала, его теплопроводность значительно меняется (»30% в интервале от 10-5 Па до 105 (1 атм.) в зависимости от давления воздуха внутри материала (рис. 3). Теплопроводность, Вт/м·К [Thermal conductivity, W/m·K] 0,6 0,5 0,4 0 Па [Pa] 101325 Па [Pa] 0,3 0,2 0,1 0 0 500 1000 1500 2000 Температура, °С [Temperature, °C] Рис. 3. Зависимость теплопроводности материала от температуры и давления [Fig. 3. Temperature and pressure dependence of the porous material] Исходя из особенностей свойств материала и теплового нагружения, в представленной математической модели использовались следующие допущения: 1. давление воздуха внутри теплозащитного материала кромки считалось близким к статическому атмосферному давлению на высоте полета; 2. тепловой поток на всей поверхности кромки считается ламинарным; 3. температура конструкции при t = 0 с (т.е. перед началом спуска с орбиты 150 км) принималось равной 30 °С; 4. излучательная способность внешней поверхности кромки принималась равной ε = 0,8 [25]. Расчетная математическая модель Для определения рациональной конструкции кромки была разработана математическая модель в пакете конечно-элементного (КЭ) анализа ANSYS Transient Thermal 16.2. С помощью программного кода ADPL (Ansys Parametric Design Language), учитывалась зависимость тепловодности материала от температуры и давления, λ = f(T, P), на всей траектории спуска. Функция f(T, P) задавалась эмпирическим уравнением, с помощью кривых теплопроводности для двух крайних значений давления (» 0 Па и 105 Па) по формуле [26] l(T , P ) = l1(T ) l1(T ) - l0 (T ) 1+ 0,656×10-3 × P ⎛1+ , 124 ⎞ ⎝⎜ T + 273⎟⎠ где T - температура в данной точке материала, °C; P - атмосферное давление в данный момент времени, Па; λ1 - теплопроводность при P = 105 Па, Вт/м·К; λ0 - теплопроводность при P » 0 Па, Вт/м·К. Нестационарный тепловой анализ проводился для времени полета tобщ = = 2870 мин = 47,8 мин с шагом по времени в 1 с. Изменение свойств материала, а именно λ = f(T, P) по формуле , задавалось с помощью метода Singleframe Restart в программе КЭ анализа. Программа, сохранив результаты расчета после временного шага, обращалась в препроцессор (PREP7), где перезаписывались кривые теплопроводности материала, λ = f(T), характерные для новой точки на траектории (т.е. нового атмосферного давления) и расчет продолжался для нового значения времени. Граничными условиями модели были конвективный тепловой поток, qконв, излучаемый радиационный поток, qрад, а также адиабатические стенки, которыми считались внутренние части силовой конструкции (рис. 4). Подводимый конвективный тепловой поток, q, считался ламинарным и задавался функцией координат x, y и времени спуска аппарата, t. Для определения q в окрестности критической точки для заданного угла, θ, при ламинарном режиме обтекания использовалась формула (рис. 5) q(q) = 0,55 + 0,45cos(2q) при 0 £ q £ p , q0 2 где q0 - плотность теплового потока в критической точке, Вт/м2, определялась из эмпирической формулы Фэя - Риддела [27; 28]: 2,56×10-5 0 q = ¥ e w r¥ V 3,25(h - h ), he R где ρ¥ - плотность потока, кг/м3; V¥ - скорость аппарата, м/с; he - энтальпия восстановления, Дж/кг; hw - энтальпия газа при температуре стенки, Дж/кг; R - радиус закругления носка, м. q рад = f(e, T, t) q конв = f(x, y, t) z y q = 0 Адиабатические стенки (q = 0) [Adiabatic walls (q = 0)] x z Ðèñ. 4. Схема теплового нагружения передней кромки крыла МКА [Fig. 4. Boundary conditions of the RSV wing leading-edge] Плотность радиационного теплового потока определялась при постоянном значении ε = 0,8 и температуре, T, в определенной точке поверхности в момент времени t (см. рис. 4). q(q) q q 0 Ðèñ. 5. Схема теплового потока на поверхности кромки [Fig. 5. Surface heat flux distribution schematics] Плотности теплового потока были рассчитаны для кромки крыла аппроксимированной уравнением эллипсоида (рис. 6). Геометрия ТЗП представляла собой объемную форму постоянной толщиной в 0,05 м (см. рис. 6). Теплозащитное покрытие крепилось к силовой конструкции из углепластика с теплопроводностью 3,3 Вт/м·К, плотностью 1600 кг/м3 и теплоемкостью 1200 Дж/кг·К [29]. Толщина ТЗП в верхней и нижней зонах (50 мм и 80 мм) определялось исходя из максимально допустимой рабочей температуры термостойкого углепластика в пределах 250-300 °С [24]. Как будет описано далее, меньшая толщина ТЗП на подветренной части кромки связана с меньшим значением теплового потока в данной зоне. Эллипсоид [Ellipse]: 2 æ (x - a)2 ö y 1 b = çè a2 ÷ø 175 мм 200 мм 50 мм 225 мм 4 мм Термостойкий углепластик y [Heat-resistant carbon] 80 мм x z Пористая керамика [Porous ceramics] Рис. 6. Геометрия передней кромки крыла аппарата [Fig. 6. Leading-edge geometry schematic] Результаты моделирования Тепловое нагружение кромки определялось для траектории аэрокосмического аппарата схожей с полетом таких кораблей как Буран, Space Shuttle, БОР 4, БОР 5, X-37 [7; 30]. Законы изменения высоты и скорости по времени спуска показаны на рисунке 7. Были рассмотрены максимальные тепловые нагрузки, которые имеют место быть при углах атаки в 0 и 40° для верхней и нижней поверхностей кромки, соответственно. 150 35 30 Высота, км [Altitude, km] Число Маха [Mach number] 25 100 20 15 50 10 0 0 10 20 30 40 50 Время, мин [Time, min] 5 0 0 10 20 30 40 50 Время, мин [Time, min] Рис. 7. Изменение высоты и скорости аэрокосмического аппарата со временем [Fig. 7. Height and Mach number vs. time plots] Расчеты показали, что наиболее интенсивный нагрев кромки поверхности происходит при t = 26 мин (рис. 8). Максимальная температура в критической точке в этот момент времени достигает значений в 1735 °С. 1 q/q в критической точке [q/q at stagnation point] 0,8 0,6 0,4 0 0 0,2 0 0 10 20 30 40 50 Время, мин [Time, min] Рис. 8. Изменение плотности теплового потока во времени [Fig. 8. Heat flux vs. time variation at stagnation point] В зависимости от угла атаки, координаты критической точки смещаются в сторону наветренной стороны носка. Изменение плотности теплового потока вдоль поверхности при t = 26 мин для разных углов атаки показано на рисунке 9. Однако не зависимо от изменения положения критической точки вдоль кромки крыла, значение ее максимальной температуры не претерпевает существенных изменений (рис. 10). Результаты расчетов свидетельствуют, что изменение угла атаки в большей степени оказывает влияние на распределение температуры внутри конструкции. С его увеличением наблюдается рост температуры на нижней поверхности кромки, что, как следствие, требует большей толщины ТЗП в этой зоне. Поля температур в последний момент времени, t = 47,8 мин, для углов атаки 0°, 25° и 40° приведены на рисунках 11-13, соответственно. 14 [q/q at t = 1560 sec] 1 q/q при t = 1560 с 0,8 АВИАЦИОННАЯ И РАКЕТНО-КОСМИЧЕСКАЯ ТЕХНИКА 0,6 0 0 0,4 0,2 4 5 0 1 + 2 3 1 - 4 5 2 3 [q/q at t = 1560 sec] 1 q/q при t = 1560 с 0,8 0,6 0,4 0 0 4 0,2 5 0 + 2 3 1 1 - 4 5 2 3 [q/q at t = 1560 sec] 1 q/q при t = 1560 с 0,8 0,6 4 0 0,4 0 0,2 5 0 + 2 3 1 - 1 4 5 2 3 -1 -0,75 -0,5 -0,25 0 0,25 0,5 0,75 1 -1 -0,75 -0,5 -0,25 0 0,25 0,5 0,75 1 -1 -0,75 -0,5 -0,25 0 0,25 0,5 0,75 1 Относительное расстояние вдоль поверхности носка [Relative distance along the edge surface] Относительное расстояние вдоль поверхности носка [Relative distance along the edge surface] Относительное расстояние вдоль поверхности носка [Relative distance along the edge surface] à á â Ðèñ. 9. Плотность теплового потока вдоль поверхности при углах атаки 0° (à), 25° (á) è 40° (â) [Fig. 9. Surface heat flux distribution for angles of attach of 0° (à), 25° (á) and 40° (â)] Bodnya I.S., Timoshenko V.P. RUDN Journal of Engineering researches, 2018, 19 (1), 7-21 Температура в критической точке, °С [Stagnation point temperature, °C] 1800 1600 1400 1200 1000 800 600 400 200 0 0 10 20 30 40 50 Время, мин [Time, min] Условные обозначения: - угол атаки 0° [Angle of attack 0°]; - угол атаки 25° [Angle of attack 25°]; - угол атаки 40° [Angle of attack 40°] Рис. 10. Зависимость изменения температуры в критической точке со временем для различных углов атаки [Fig. 10. Temperature vs. time variation for stagnation point for various angles of attack] Рис. 11. Поле температур при t = 47,8 мин для угла атаки 0° [Fig. 11. Temperature field at t = 47,8 min for 0 degrees angle of attack] Рис. 12. Поле температур при t = 47,8 мин для угла атаки 25° [Fig. 12. Temperature field at t = 47,8 min for 25 degress angle of attack] Рис. 13. Поле температур при t = 47,8 мин для угла атаки 40° [Fig. 13. Temperature field at t = 47,8 min for 40 degrees angle of attack] Для определения влияния угла атаки на температуру силовой конструкции были построены зависимости максимальной температуры углепластика от времени полета (рис. 14). Максимальная температура углепластика, °С [Maximum carbon temperature, °C] 350 300 250 200 150 100 50 0 0 10 20 30 40 50 Время, мин [Time, min] Условные обозначения: - угол атаки 0° [Angle of attack 0°]; - угол атаки 25° [Angle of attack 25°]; - угол атаки 40° [Angle of attack 40°] Рис. 14. Изменения максимальной температуры силовой конструкции во времени для углов атаки 0°, 25° и 40° [Fig. 14. Support structure maximum temperature vs. time for 0°, 25° and 40° angles of attack] Выводы Результаты математического моделирования свидетельствуют, что максимальная температура силовой конструкции из углепластика не превышает 225 °С и 250 °С при углах атаки в 42° и 25°, соответственно. Следует отметить, что, как правило, полет МКА совершается при углах атаки в пределах от 42° до 7°. Максимальная тепловая нагрузка же имеет место еще в более узком интервале от 42° до 25° [30]. Из этого следует, что действительная температура силовой конструкции будет иметь максимальное значение, не превышающее »250 °С. Однако в реальных условиях, при наличии излучения во внутреннюю полость крыла, ожидается дополнительное снижение температуры силовой конструкции на 5-10%. Следовательно, максимальная температура силовой конструкции на всей траектории полета будет находиться в пределах рабочих температур термостойкого углепластика до 300 °С [24]. Из полученных результатов сделан вывод, что пористая керамика на основе волокон Al2O3 способна обеспечить необходимую тепловую защиту корабля для заданной траектории и скорости полета.

Ivan S Bodnya

Bauman Moscow State Technical University (National research university of technology)

Author for correspondence.
Email: ivanbodnya@gmail.com
5/1, 2-nd Baumanskaya str., Moscow, 105005, Russian Federation

master student at Bauman Moscow State Technical University in the Mechanical engineering department SM-13 “Space-Rocket Composite Designs”. Research interests: heat transfer, thermal regimes of space vehicles

Valery P Timoshenko

Bauman Moscow State Technical University (National research university of technology)

Email: moltim@yandex.ru
5/1, 2-nd Baumanskaya str., Moscow, 105005, Russian Federation

professor at Bauman Moscow State Technical University in the Mechanical engineering department SM-13 “Space-Rocket Composite Designs”. Research interests: heat transfer, thermal protection of space vehicles, space vehicles tests.

  • Dumbacher D. NASA’s Second Generation Reusable Launch Vehicle Program Introduction, Status and Future Plans. In: 38th AIAA/ASME/SAE/ASEE Joint Propulsion Conference & Exhibit [Internet]. Huntsville, Alabama: AIAA; 2002. Available from: http://dx.doi. org/10.2514/6.2002-3613
  • Wang Z., Huang S., Shen L., Zhou H., Zhi J. Conceptual evaluation of multi-purpose aerospace plane. Russian-American scientific journal: Actual problems of aviation and aerospace systems: processes, models, experiment. Nanjing, China, 2014, 1(38), 45—57.
  • Johnson S.M. Thermal Protection Materials: Development, Characterization and Evaluation. Munich, Germany. 2012.
  • Daryabeigi K., Knutson J.R., Cunnington G.R. Heat Transfer Measurement and Modeling in Rigid High-Temperature Reusable Surface Insulation Tiles. AIAA. 2011; 345: 2011.
  • Nanowick L., Flow C. Lightweight Thermal Protection System for Atmospheric Entry. NASA Tech Briefs. 2007; (October 2007): 20–21.
  • Desarrollo y Defensa [Internet]. [cited 2017 Oct 13]. Available from: http://desarrolloydefensa. blogspot.ru/2017_07_09_archive.html (accessed: 05.10.2017).
  • Lozino-Lozinsky G., Timoshenko V. “Lessons learned from the BOR flight campaign”. In: Proceedings of the 3rd European Symposium on Aerothermodynamics for space vehicles. ESTEC; 1999. Pp. 9.
  • Gofin M. Thermal protection and hot structures of reusable space vehicles. Moscow: MIR publ. 2003. 637. (In Russ.).
  • Glass D.E. Ceramic Matrix Composite (CMC) Thermal Protection Systems (TPS) and Hot Structures for Hypersonic Vehicles. Seminar. 2008; 2682(May): 1—36.
  • Jenkins D.R., Landis T., Miller J. American X-Vehicles: An Inventory-X-1 to X-50. Monographs in Aerospace History. 2003. 65.
  • Palmer G., Polsky S. Aerothermal Analysis of the X-34 Vehicle. Access in Space. 1998;(January 2014): 84—86.
  • Hilfer G. Flight Qualification Testing of X-38 TPS Components Lessons Learned. In: A. Wilson, editor. Hot Structures and Thermal Protection Systems for Space Vehicles, Proceedings of the 4th European Workshop. Palermo, Italy: European Space Agency; p. p. 169.
  • Stewart D.A., Leiser D.B., DiFiore R.R., Katvala V.W. High efficiency tantalum-based ceramic composite structures [Internet]. Vol. 1. 2010. Available from: http://www.google.com/patents/ US7767305%5Cnhttp:// patentimages.storage.googleapis.com/pdfs/US7767305.pdf (accessed: 25.09.2017).
  • Ralf R., I-Wei C. Ceramics Science and Technology. Volume 1: Structures. In Wiley-VCH; 2008. p. 565—566.
  • European Space Agency. Skylon Assessment Report. Noordwijk, Netherlands; 2011.
  • Kuczera H., Sacher P.W. Reusable Space Transportation Systems [Internet]. Berlin: Springer Berlin Heidelberg; 2011. Available from: http://www.springer.com/us/book/9783540891802 (accessed: 01 Oct 2017).
  • Rufolo G., Roncioni P., Marini M. USV FTB-1 Reusable vehicle aerodatabase development. 2007.
  • NASA, SNC. Photo of Dream Chaser [Internet]. Available from: https://www.nasa.gov/sites/ default/files/2013-3230_0.jpg (accessed: 05 Oct 2017).
  • Pezzella G., Battista F., Schettino A., Marini M., Matteis P.De. Hypersonic Aerothermal Environment Preliminary Definition of the Cira Ftb-X Reentry Vehicle. Environment. 2007; (November): 1—25.
  • Viviani A., Pezzella G. Heat Transfer Analysis for a Winged Reentry Flight Test Bed. International Journal of Engineering. 2009; 3(3): 329—345.
  • Ii F., Ingegneria F.DI, Di D., In R., Aerospaziale I. a Study of a High Lift Wing-Body Configuration for Low Earth Orbit Re-Entry. 140.
  • Tatsuki O., Mrityunjay S. Engineered Ceramics: Current Status and Future Prospects. New Jersey: John Wiley & Sons; 2015. 232 p.
  • Gribkov B.N., Mizurina G.T., Shetanov B.V., Lyapin V.V. Possibilities of fibrous thermal protection. In: Proceedings of the First International Aviation Conference ‘Man-Earth-Space’. Moscow: Russian Engineering Academy. Sec. ‘Aerospace’; 1995. p. 223—231.
  • Shalin R.E., Zinoviev S.N., Pomerantsev K.P., Moiseev E.B., Shepeleva L.I. Thermostable carbon plastic KMU-8. Aviation industry, 1987, (5), 53—55. (In Russ.).
  • Stewart D.A., Leiser D.B. Toughened Uni-piece, Fibrous, Reinforced, Oxidization-Resistant Composite. Vol. 1. USA; 7314648, 2008.
  • Kostylev V.M. Thermal conductivity of dispersed bodies at different atmospheric pressure. Thermal physics of high temperatures. 1964. 2(1), 1—18. (In Russ.).
  • Surzhikov S.T., Shuvalov M.P. Analysis of radiation-convective heating of four types of descent space vehicles. Physico-chemical kinetics in gas dynamics. 2014. 15(4). 1—18 (in Russ.).
  • Bobylev A.V., Vaganov A.V., Dmitriev V.G., Zadonsky S.M., Kireev A.Y., Skuratov A.S., et al. Development of aerodynamic configuration and research of aerothermodynamic characteristics of a small-sized winged space vehicle. Scientific notes TsAGI. 2009, XL(3), 3—15. (In Russ.).
  • Denisov O., Minakov D., Kirbay A. Methodical Specifics of Thermal Experiments with Thin Carbon Reinforced Plates. Science and Education of the Bauman MSTU [Internet]. 2015. (7). 171—184. Available from: http://technomag.bmstu.ru/doc/781946.html (accessed: 13.08.2017).
  • Nguyen V., Poladian D., Falangas E., Chaudhary A., Tran H. Dynamics and stability and control characteristics of the X-37. American Institute of Aeronautics and Astronautics. 2001; 1—10.

Views

Abstract - 85

PDF (Russian) - 30

PlumX


Copyright (c) 2018 Bodnya I.S., Timoshenko V.P.

Creative Commons License
This work is licensed under a Creative Commons Attribution 4.0 International License.