Теплофизика высоких температур, 2019, T. 57, № 4, стр. 509-517

ВВЕДЕНИЕ

Безэлектродные ВЧ-плазмотроны [1–7] предоставляют широкие возможности для экспериментальных исследований теплообмена высокоэнтальпийных потоков молекулярных газов и испытаний теплозащитных материалов на термохимическую стойкость в условиях нагрева, соответствующих параметрам входа затупленных тел в атмосферу. Условия теплообмена в высокоэнтальпийных струях воздуха можно менять за счет варьирования энтальпии (5–50 МДж/кг), давления торможения (0.01–1 атм) и геометрии моделей. Представляется важным и интересным вопрос о том, как изменяются параметры теплообмена при варьировании геометрии разрядного канала ВЧ-плазмотрона и неизменной форме модели.

В данной работе экспериментальное исследование теплообмена в дозвуковых струях высокоэнтальпийного воздуха проведено на 100-киловаттном индукционном плазмотроне ВГУ-4 (ИПМех РАН) для трех конфигураций секционированного разрядного канала, состоящего из кварцевой трубки диаметром 80 мм и сменных водоохлаждаемых конических сопел в конце канала с углом полураствора 30°. Диаметры выходных сечений использованных насадков D = 30, 40 и 50 мм. Для всех трех конфигураций канала проведена экспериментальная диагностика в широком диапазоне мощности ВЧ-генератора плазмотрона. Для условий экспериментов проведено численное моделирование течения плазмы в разрядном канале, обтекания модели дозвуковой струей высокоэнтальпийного воздуха и химически неравновесного течения в пограничном слое на оси симметрии потока вблизи критической точки модели. В результате для условий экспериментов из сопоставления экспериментальных и расчетных данных определены значения энтальпии на внешней границе пограничного слоя и скорости набегающего потока на оси симметрии.

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

ЭКСПЕРИМЕНТ НА ВЧ-ПЛАЗМОТРОНЕ

Экспериментальное исследование теплообмена в дозвуковых струях воздуха проведено на 100‑киловаттном индукционном плазмотроне ВГУ-4 для трех конфигураций секционированного разрядного канала с диаметром выходного сечения D = 30, 40 и 50 мм. Исследовалось обтекание модели из меди – водоохлаждаемого цилиндра с плоским торцом диаметром d_m = 20 мм, причем расстояние от среза канала до модели равно диаметру выходного сечения сопла L_z = D. Во всех экспериментах статическое давление в испытательной камере p = 100 гПа, расход воздуха G = 2.4 г/с, мощность генератора по анодному питанию N_ap менялась от 20 до 75 кВт.

Экспериментальная диагностика включала измерения статического давления p_∝ у стенки испытательной камеры, определение с помощью трубки Пито динамического давления p_dyn на оси струи, натекающей на модель, а также измерение теплового потока q_w к медному водоохлаждаемому калориметру в окрестности критической точки на торцевой поверхности модели. Для определения N_pl – мощности, вкладываемой в плазму, – использовались экспериментальные данные для к.п.д. η_eff, полученные ранее для соответствующих режимов работы ВЧ-плазмотрона на воздухе.

Значения полученных в экспериментах параметров q_w и p_dyn приведены в табл. 1–3 для трех серий экспериментов с использованием конических водоохлаждаемых сопел с диаметрами D = = 30, 40 и 50 мм. Эти значения использовались далее как исходные данные при численном моделировании течения в плазмотроне, расчетном восстановлении энтальпии на внешней границе пограничного слоя h_e и характерной скорости дозвукового потока V_s для условий экспериментов. Значения h_e, V_s, а также вычисленный эффективный безразмерный радиус r_eff цилиндрической модели с плоским носком радиусом 10 мм приведены в табл. 1–3. Отметим, что радиусы r_eff для данной модели в случае сопел с D = 40 и 50 мм практически совпадают. В случае сопла с D = 30 мм эффективный радиус r_eff до 9% превышает соответствующие значения из табл. 2 и 3.

ЧИСЛЕННОЕ МОДЕЛИРОВАНИЕ ДОЗВУКОВЫХ ТЕЧЕНИЙ ПЛАЗМЫ И ВЫСОКОЭНТАЛЬПИЙНОГО ВОЗДУХА В РАЗРЯДНОМ КАНАЛЕ И В СТРУЕ, ОБТЕКАЮЩЕЙ МОДЕЛЬ

Численное моделирование течений низкотемпературной плазмы молекулярных газов в установках с индукционным нагревом является неотъемлемой частью диагностики в методах восстановления аэротермодинамических полей для условий экспериментов и определения каталитических характеристик материалов при высоких температурах поверхности [8–17]. Настоящая методика численного моделирования основана на разбиении поля течения в плазмотроне ВГУ-4 на три характерные области [9]: 1) течение равновесной воздушной плазмы в разрядном канале (расчет по программе Alpha на основе уравнений Навье–Стокса); 2) течение дозвуковой струи диссоциированного равновесного воздуха, истекающей из канала и обтекающей испытываемую модель (расчет по программе Beta на основе уравнений Навье–Стокса); 3) течение в неравновесном пограничном слое на оси симметрии у поверхности вблизи критической точки модели (расчет по программе Gamma на основе уравнений пограничного слоя конечной толщины).

Численное моделирование дозвукового течения воздушной плазмы в разрядном канале проводится в предположении, что плазма является термически и химически равновесной, это справедливо при давлении в разрядном канале более 50 гПа. Предполагается также, что значения мощности N_pl, вкладываемой в плазму, не слишком велики, так что можно пренебречь потерями энергии, связанными с излучением плазмы [17, 18]. Течение считается стационарным, ламинарным и осесимметричным с закруткой потока в азимутальном направлении. Для расчета используются полные уравнения Навье–Стокса в цилиндрической системе координат (z, r) с учетом трех компонент скорости: аксиальной, радиальной и азимутальной [8–10]. Уравнения включают источниковые члены, учитывающие влияние электромагнитного поля: аксиальную и радиальную компоненты силы Лоренца и джоулево тепло, которые определяются усредненным по времени вихревым электромагнитным полем и выражаются через тангенциальную компоненту комплексной амплитуды вихревого электрического поля E_θ.

Граничные условия для данной задачи: во входном сечении задаются все параметры газа, включая три компоненты скорости, определяемые условиями эксперимента (заданный массовый расход газа и закрутка потока на входе в канал 45°); нулевые значения для компонент скорости и определенные значения температуры задаются на всех твердых поверхностях (стенки кварцевой трубки и сопла, торец входного интерфейса), температура стенки разрядного канала T_w = 300 K; условия симметрии используются на оси канала; известные “мягкие” граничные условия (равенство нулю производных от газодинамических параметров) применяются в выходном сечении канала.

Считается, что индуктор состоит из отдельных бесконечно тонких круговых витков, перпендикулярных оси симметрии, с одинаковым током в каждом витке. Переменный внешний ток порождает монохроматическое электрическое поле с комплексной амплитудой E_θ(z, r). С использованием предположений, что плазма квазинейтральна и током смещения можно пренебречь, электромагнитное поле внутри индуктора определяется двумерным уравнением Максвелла для комплексной тангенциальной компоненты амплитуды электрического поля E_θ [8, 9]. С учетом особенности разрядного канала плазмотрона ВГУ-4, который является относительно тонким, предполагается также, что электрическое поле меняется в продольном направлении значительно меньше, чем поперек канала. В этом случае двумерные уравнения существенно упрощаются и сводятся к квазиодномерному приближению для уравнений Максвелла [8, 9], которое позволяет определить комплексную амплитуду электрического поля E_θ(z, r). Упрощенная квазиодномерная задача существенно более простая, чем двумерная задача, в то же время она является достаточно точной для относительно тонкого канала ВГУ-4. Сравнение численных решений для упрощенного уравнения и решений с точными расчетами по двумерной модели для электрического поля показало, что отличие параметров плазмы в области индуктора составляет всего несколько процентов, а на срезе канала результаты практически совпадают [16].

Численное решение уравнений Навье–Стокса основано на методе конечных объемов и алгоритме SIMPLE Патанкара–Сполдинга с применением шахматной сетки [19], сильно неравномерной как в аксиальном, так и в радиальном направлении. Коэффициенты переноса (вязкость, теплопроводность, электропроводность) определялись в процессе расчета с помощью интерполяции по заранее рассчитанным таблицам. Таблицы коэффициентов переноса рассчитывались по программе SoVA [20].

С помощью описанного метода по программе Alpha проведены расчеты течения воздушной плазмы в трех разрядных каналах с коническими соплами (D = 30, 40 и 50 мм) для следующих условий экспериментов на плазмотроне ВГУ-4: давление p = 100 гПа; мощность, вкладываемая в плазму, N_pl = 29 кВт; массовый расход воздуха G = 2.4 г/с. На рис. 1 показаны результаты расчета – изотермы для случая сопла с D = 50 мм.

Расчет течения во второй области – в испытательной камере плазмотрона, где истекающая из разрядного канала дозвуковая струя обтекает модель, – проводится по программе Beta также на основе полных уравнений Навье–Стокса. Влияние электромагнитного поля в этой области не учитывается. Основное внимание при расчете течения во второй области уделяется ядру струи вблизи передней критической точки модели, включая пограничный слой; здесь используется более мелкая подробная сетка. Для численного решения уравнений Навье–Стокса во второй области применяется тот же численный метод [19], что и в канале плазмотрона. В качестве граничных условий используются профили параметров газа на срезе разрядного канала, полученные на предыдущем этапе из расчета по программе Alpha. На твердых поверхностях скорость газа равна нулю и задается температура T_w = 300 K. На выходной границе области используются “мягкие” граничные условия, которые можно применять, если длина области течения существенно превышает ее ширину, что накладывает ограничение на длину расчетной области L_V ≥ 1.5D_V, она должна быть больше, чем 1.5 диаметра камеры D_V.

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

По описанному методу с помощью программы Beta проведены расчеты обтекания модели дозвуковой высокоэнтальпийной струей воздуха для экспериментальных условий, указанных выше: давление p = 100 гПа; мощность, вкладываемая в плазму, N_pl = 29 кВт; массовый расход воздуха G = = 2.4 г/с. Геометрия модели, расположенной на оси симметрии, – цилиндр с плоским торцом диаметром d_m = 20 мм; температура его поверхности T_w = 300 К. На рис. 2 показаны рассчитанные изотермы течения, при этом представлена область ядра струи, включая относительно тонкий пограничный слой вблизи переднего торца модели.

Основными результатами расчета по программе Beta являются не только параметры потока газа в ядре струи перед моделью, но и вычисляемые по этим параметрам безразмерные величины, характеризующие толщину пограничного слоя и завихренность потока на внешней границе пограничного слоя перед моделью. Четыре безразмерные величины определяются в результате расчета по программе Beta по следующим формулам: δ/R_m – безразмерная толщина пограничного слоя; u(Z_e, 0)/U₀ – безразмерная скорость; u_1e = u₁(Z_e) – безразмерная производная радиальной скорости; ${{R}_{m}}{{\left. {\frac{{d{{u}_{1}}(z)}}{{dz}}} \right|}_{{z = {{Z}_{e}};r = 0}}}$ – параметр, характеризующий завихренность потока. Здесь нижний индекс e обозначает величины на внешней границе пограничного слоя; δ – размерная толщина пограничного слоя; R_m – радиус модели; Z_e – продольная координата, соответствующая положению внешней границы пограничного слоя; U₀ – характерная скорость, равная скорости потока на срезе разрядного канала на оси симметрии; u(z, r) – компонента скорости струи в продольном направлении вдоль оси симметрии; z – координата вдоль оси симметрии; r – радиальная координата, т.е. расстояние от оси симметрии струи. Четыре указанных выше безразмерных параметра определяются на внешней границе пограничного слоя на оси струи: Z = Z_e, r = 0. Безразмерная производная радиальной скорости u₁(z) рассчитывается по следующей формуле:

где $\text{v}$(z, r) – компонента скорости струи в радиальном направлении. Внешняя граница пограничного слоя z = Z_e определяется по результатам расчета обтекания модели дозвуковой струей воздуха (программа Beta) из условия максимума производной от u₁(z) на оси струи. Касательная к u₁(z), проведенная через точку (Z_e, u₁(Z_e)), имеет максимальный наклон к оси z (см. рис. 3). Эта касательная является локальной асимптотикой “невязкого” (слабо вязкого) течения, на которую выходит расчет по программе Gamma вблизи внешней границы пограничного слоя при постановке соответствующих граничных условий. На рис. 3 (квадратами) показаны значение u₁(z) на внешней границе пограничного слоя (u_1e) и экстраполированное по касательной на стенку значение u_1w.

Указанные четыре безразмерные величины рассчитаны для всех вариантов, т.е. для трех значений диаметра выходного сечения конического сопла и для всех значений мощности плазмотрона.

ЧИСЛЕННОЕ МОДЕЛИРОВАНИЕ ТЕЧЕНИЯ И ТЕПЛОПЕРЕДАЧИ В НЕРАВНОВЕСНОМ ПОГРАНИЧНОМ СЛОЕ ВБЛИЗИ КРИТИЧЕСКОЙ ТОЧКИ МОДЕЛИ

Пограничный слой для условий эксперимента на ВЧ-плазмотроне ВГУ-4 обычно существенно неравновесный, а при низких давлениях – замороженный. Поверхность модели имеет заданную температуру T_w и характеризуется эффективным коэффициентом каталитической рекомбинации атомов γ_w. Численное моделирование течения и теплопередачи от химически неравновесного пограничного слоя к критической точке модели реализовано в компьютерной программе Gamma [9, 10]. Важная особенность этой программы, позволяющая получать достаточно точные результаты, – тесная интеграция с надежными экспериментальными данными: тепловым потоком q_w, измеренным в критической точке водоохлаждаемой модели медным калориметром, и динамическим давлением p_dyn, измеренным с помощью водоохлаждаемой трубки Пито. Подробное описание датчиков и методов измерения теплового потока и давления приведено в [21, 22].

Программа Gamma основана на численном решении уравнений пограничного слоя на оси симметрии потока с учетом каталитических реакций рекомбинации атомов O и N на поверхности модели, конечной толщины пограничного слоя и завихренности потока на его внешней границе. Модель течения в третьей области построена на основе системы одномерных уравнений пограничного слоя конечной толщины, в которой численное решение уравнений Навье–Стокса используется как внешнее решение [9]. Применяется 5-компонентная модель диссоциированного воздуха (N₂, O₂, NO, N, O).

Отметим, что некоторые обозначения в данном разделе отличаются от использовавшихся в предыдущих разделах в силу специфики третьей расчетной области и используемых здесь уравнений пограничного слоя. Например, компонента скорости V в обозначениях данного раздела соответствует обозначению u из предыдущих разделов, V_s – характерное значение скорости потока вдоль оси симметрии. В обозначениях данного раздела компонента скорости U в направлении x вдоль поверхности модели (т.е. перпендикулярно оси симметрии) равна нулю на оси симметрии, и в уравнения входит градиент этой компоненты скорости u₁. Уравнения пограничного слоя включают безразмерные параметры, полученные ранее из расчета обтекания модели с помощью программы Beta. Для уравнений сохранения импульса и энергии использовались граничные условия прилипания и заданная температура поверхности T_w. При постановке граничных условий для уравнений диффузии компонент на внешней границе пограничного слоя предполагается, что течение во внешнем потоке равновесное. Равновесные концентрации компонент C_ie вычисляются с помощью уравнений химического равновесия для заданных значений давления и температуры (или энтальпии) на внешней границе пограничного слоя с помощью программы SoVA [20].

Чтобы получить граничные условия для уравнений диффузии компонентов на поверхности тела используется простая однопараметрическая модель рекомбинации атомов N и O, полученная в рамках следующих предположений: 1) каталитическая рекомбинация атомов N и O на поверхности (S) описывается реакциями первого порядка

2) концентрация молекул NO мала, и они не участвуют в каталитических реакциях на поверхности; 3) вероятность рекомбинации атомов N и O на поверхности одинакова и равна эффективному коэффициенту γ_wN = γ_wO = γ_w; 4) диффузионные потоки к поверхности пропорциональны только концентрации соответствующей компоненты (модель гулардовского типа), что в рассматриваемом случае дает следующее выражение для массовых диффузионных потоков атомов J_i:

где C_i – массовая концентрация i-й компоненты, ρ – плотность смеси, K_wi – константа скорости гетерогенной рекомбинации i-й компоненты, которая определяется через коэффициенты γ_wN и γ_wO [9, 10, 21, 22]:

Численное решение уравнений химически неравновесного пограничного слоя с учетом его конечной толщины находится с помощью разностной схемы [23] с четвертым порядком аппроксимации. В результате расчета определяются профили параметров поперек пограничного слоя: температуры, градиента скорости, концентраций компонент. Основным результатом расчета является значение теплового потока в критической точке модели q_w, зависящее от температуры стенки T_w и от заданного значения эффективного коэффициента рекомбинации γ_w. Расчет q_w для заданных экспериментальных условий является основой метода восстановления энтальпии набегающего потока.

ОПРЕДЕЛЕНИЕ ЭНТАЛЬПИИ И СКОРОСТИ ДОЗВУКОВОГО ПОТОКА ВЫСОКОЭНТАЛЬПИЙНОГО ВОЗДУХА

Расчетное определение энтальпии с помощью программы Gamma базируется на заданных значениях теплового потока q_w к медному водоохлаждаемому калориметру и динамического давления p_dyn, измеренного с помощью водоохлаждаемой трубки Пито, имеющей ту же геометрию, что и модель с калориметром для измерения теплового потока. Энтальпия определяется из условия равенства измеренного теплового потока q_w и теплового потока q_wc, полученного из расчета по программе Gamma для условий эксперимента: q_w = q_wc(h_e, T_w, γ_w). Рассчитанный тепловой поток q_wc зависит от температуры поверхности модели T_w, коэффициента каталитичности γ_w и энтальпии h_e. При расчетном определении энтальпии струи высокоэнтальпийного воздуха используется предположение, что поверхность химически чистой меди, из которой сделан водоохлаждаемый калориметр, является высококаталитичной и для нее можно принять γ_w = 1. Реальное значение γ_w для меди несколько меньше, чем 1.0, но при испытаниях в дозвуковых струях зависимость теплового потока q_w от γ_w является достаточно слабой в диапазоне 0.1 ≤ γ_w ≤ 1 (см. рис. 2 в работе [21]). Поэтому предположение γ_w = 1 для медного калориметра вносит погрешность в определение энтальпии h_e, сравнимую с погрешностью экспериментальных данных и расчета, если испытания проводятся в дозвуковой струе.

В указанных предположениях определение энтальпии сводится к решению методом итераций обратной задачи

где q_w – тепловой поток, измеренный медным водоохлаждаемым калориметром. При решении задачи (1) на каждой итерации определяется не только h_e, но также и необходимое для расчета значение характерной скорости дозвуковой струи V_s. Оно вычисляется по предварительно измеренному динамическому давлению p_dyn с использованием полуэмпирической формулы где ρ_e – плотность на внешней границе пограничного слоя, k_P = k_P(Re) –коэффициент, зависящий от числа Рейнольдса. В данных расчетах полагалось k_P = 1.1. Значения h_e и V_s, рассчитанные по программе Gamma для трех вариантов разрядного канала D = 30, 40 и 50 мм, приведены в табл. 1–3 и показаны на рис. 4 и 5. Как видно из этих рисунков, различная геометрия разрядного канала (D = = 30, 40, 50 мм) ожидаемо приводит к большому различию в скоростях потоков на оси струи, но при этом зависимости энтальпии от мощности h_e(N_ap) различаются незначительно для разных D.

ПЕРЕСЧЕТ ПАРАМЕТРОВ ТЕПЛООБМЕНА В ДОЗВУКОВЫХ ВЫСОКОЭНТАЛЬПИЙНЫХ СТРУЯХ ВОЗДУХА НА УСЛОВИЯ ВХОДА ЗАТУПЛЕННЫХ ТЕЛ В АТМОСФЕРУ ЗЕМЛИ

Вопросы соответствия условий теплообмена в дозвуковой воздушной струе ВЧ-плазмотрона параметрам входа тел в атмосферу Земли в случае осесимметричной конфигурации обтекания цилиндрической модели в эксперименте решаются на основе теории локального термохимического моделирования [24, 25]. Теория разработана для моделирования теплообмена в критической точке затупленного тела, обтекаемого высокоскоростным потоком газа, и основана на требовании одновременного воспроизведения в эксперименте полной энтальпии набегающего потока, давления торможения и градиента скорости. Масштабным фактором моделирования ζ является отношение эффективных радиусов модели и тела, которые выражаются через радиусы затупления поверхностей модели и тела и градиенты скорости в критической точке. При выполнении условий моделирования и одинаковых граничных условиях теплообмена в дозвуковых высокоэнтальпийных потоках на поверхности модели в окрестности критической точки моделируется не только тепловой поток, но локально моделируется и реагирующий неравновесный пограничный слой, включая диффузионные потоки, газофазные и поверхностные реакции.

Выбор режима локального моделирования теплообмена высокоскоростного потока с поверхностью затупленного тела в критической точке определяется критерием, полученным в [24]. Моделирование в дозвуковых высокоэнтальпийных потоках должно применяться в случае, когда масштабный фактор ζ достаточно мал:

где κ_eq, ${{\kappa }_{*}}$ – равновесный и эффективный показатели адиабаты; индекс 2 относится к параметрам высокоэнтальпийного потока в установке. В частности, это случай, когда радиус модели много меньше радиуса затупления тела. Если эффективные радиусы модели и тела – величины одного порядка, то моделирование должно осуществляться в сверхзвуковом режиме обтекания модели высокоэнтальпийным потоком.

В случае М₂ $ \ll $ 1 соотношения между параметрами V_1∞ (скорость) и ρ_1∞ (плотность) высокоскоростного потока, обтекающего тело со сферическим затуплением радиуса R_b, и параметрами h_2e (энтальпия), V_2s (скорость), p_2s (давление торможения) дозвукового потока, обтекающего модель с эффективным радиусом $R_{m}^{*}$, имеют вид [24]

Для заданной геометрии осесимметричной модели эти соотношения позволяют однозначно пересчитать условия эксперимента на ВЧ-плазмотроне на параметры входа тела в атмосферу. С помощью формул (2) на основании результатов численных расчетов течений воздушной плазмы в разрядном канале ВЧ-плазмотрона, расчетов обтекания цилиндрической модели с плоским торцом диаметром 20 мм дозвуковыми высокоэнтальпийными струями воздуха, а также с учетом результатов восстановления энтальпии и скорости дозвуковых потоков воздуха, приведенных в табл. 1–3, параметры экспериментов в ВЧ-плазмотроне пересчитаны на условия движения затупленного тела в атмосфере Земли. Результаты экстраполяции параметров экспериментов на условия входа в атмосферу приведены в табл. 4–6.

Из табл. 4–6 видно, что геометрия сопла влияет главным образом на радиус затупленного тела в критической точке: с увеличением радиуса выходного сечения конического сопла увеличивается соответствующий радиус R_b. При этом моделируемые диапазоны скоростей входа и высот слабо зависят от диаметра D в диапазоне D = 30–50 мм.

ЗАКЛЮЧЕНИЕ

Исследовано обтекание цилиндрической модели с плоским торцом в дозвуковых струях высокоэнтальпийного воздуха, истекающих из разрядного канала плазмотрона, снабженного сменными коническими сопловыми насадками с диаметром выходного сечения D = 30, 40 и 50 мм. Проведена экспериментальная диагностика и численное моделирование течения в канале, обтекания модели, течения и теплообмена в пограничном слое вблизи критической точки модели. Расчетным путем на основе экспериментальных данных определены значения энтальпии и скорости потока на внешней границе пограничного слоя перед критической точкой модели. С помощью теории локального моделирования теплообмена рассчитаны высота полета, скорость полета и радиус кривизны носка тела, соответствующие условиям теплопередачи при обтекании модели струей, истекающей из разрядного канала плазмотрона ВГУ-4 для условий проведенных экспериментов.

Работа выполнена по теме госзадания (№ АААА-А17-117021310383-2) при частичной поддержке РФФИ (грант № 17-01-00054).