Введение. При эксплуатации ракеты-носителя (РН) в элементах ее конструкции возникают нагрузки. Речь идет о продольных и перерезывающих силах, изгибающих и крутящих моментах. Данные об указанных силовых факторах задействуют для прочностного анализа при проектировании новых изделий, экспериментальной отработке конструкции [1] и адаптации средств выведения под конкретный пуск [2]. Нагрузки подразделяются на квазистатические и динамические. Квазистатические, возникающие в полете, вычисляются из условия динамического равновесия РН как твердого тела с учетом допустимых параметров траектории выведения.

Для расчета динамических нагрузок обычно применяются обобщенные балочные модели. Такие нагрузки определяются по результатам решения уравнения движения упругой РН, которое в общем случае представляет собой уравнение в частных производных. Методы, основанные на использовании разложений решений по тонам собственных колебаний конструкции, способны обеспечить высокую скорость вычисления в сочетании с достаточной точностью результатов [3]. Как показано в [4], с помощью метода разложения по формам собственных колебаний можно перейти к системе независимых обыкновенных дифференциальных уравнений. Они описывают:

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

С целью расчета нагрузок в полете рассматривается не вся траектория выведения, а только некоторые ее точки, так называемые случаи нагружения, характеризующиеся экстремальным значением отдельных параметров, влияющих на нагружение, либо максимальным значением нагрузок на отдельные элементы конструкции. Один из наиболее важных случаев нагружения — полет РН в плотных слоях атмосферы [1, 5]. Влияние неспокойной атмосферы на нагружение РН можно определить статистическими методами [6, 7] либо в рамках консервативного подхода, когда учитываются максимально возможные (с некоторым уровнем вероятности) характеристики ветра. В данной статье рассматривается второй подход. В качестве внешнего динамического воздействия принимается однократный нормативный порыв ветра. Профиль нормативного порыва, характеризующий изменение скорости ветра с течением времени, может быть задан в трапецеидальном [8], косинусоидальном или синусоидальном виде [9]. В данной статье рассмотрим движение ракеты-носителя под действием трапецеидального порыва. Продолжительность нормативного порыва обычно выбирается сопоставимой с периодом низшего поперечного тона колебаний РН1. При этом часто выдвигаются требования ее варьирования для достижения максимальных усилий в сечениях РН [8, 10]. Трудоемкость расчетов с использованием стандартных программ конечноэлементного (КЭ) анализа обусловлена:

Цель настоящей статьи — разработка методики выбора продолжительности действия нормативного порыва с использованием аналитических решений, полученных для упрощенной динамической модели РН, представленной в виде эквивалентной системы осцилляторов.

Полуаналитический подход с использованием интеграла Дюамеля был успешно применен в [11] для проведения гидроупругого анализа судов. В [12] и ряде других работ интеграл Дюамеля используется в рамках решения задачи нагружения мостов подвижными нагрузками. В настоящей статье интеграл Дюамеля задействовали для аналитических решений реакции ракеты-носителя на кратковременное воздействие порыва ветра в полете.

Материалы и методы. На этапах эскизного проектирования целесообразно использовать плоские расчетные схемы для балочных моделей. При простоте и скорости решения они позволяют определить параметры движения и внутренние усилия (с приемлемой для этого этапа проектирования точностью) [13]. Представим РН в виде упругой балки с переменными по длине массой и жесткостью. Примем обычные для сопротивления материалов допущения, в том числе гипотезу о малости упругих деформаций. Для определения внутренних силовых факторов в сечении РН используем метод ускорений (перегрузок) [1, 4], который можно интерпретировать как метод сечений, адаптированный для динамического расчета. В этом случае внутренние усилия находят из условий статического равновесия мысленно отсеченных частей конструкции под действием внешних распределенных нагрузок, дополненных силами инерции Д’Аламбера, и искомых внутренних усилий. Квазистатические и динамические значения силовых факторов определяются отдельно на основе предварительно вычисленных ускорений, а затем суммируются [4].

В данной работе рассматривается вопрос определения динамического нагружения РН в поперечном направлении под действием порыва ветра, скорость которого направлена перпендикулярно продольной оси РН. Принимается, что нагружение в продольном направлении можно рассчитать независимо. В данной работе оно не рассматривается.

С целью определения динамических ускорений РН представляется в виде свободной упругой балки. Ее движение изучается в окрестности момента времени, соответствующего рассматриваемому случаю нагружения, и описывается в отклонениях от состояния динамического равновесия, в котором находилась РН до ветрового порыва, двигаясь по номинальной (невозмущенной) траектории. При этом как бы «замораживаются», т. е. принимаются постоянными и равными характеристикам рассматриваемой точки номинальной траектории такие параметры, как масса и момент инерции РН, тяга двигателя, угол наклона траектории. Возмущенное движение упругой РН исследуется в неподвижной системе координат, связанной с положением, которое РН занимала на момент начала расчета. Возмущенное движение будет представлять собой совокупность плоскопараллельного движения РН как твердого тела в плоскости, в которой приложена динамическая нагрузка, и упругих движений корпуса. Реакцию системы управления не учитываем, т. е. для автомата стабилизации предполагаем большое время запаздывания по сравнению со временем приложения динамической нагрузки. Вообще, под воздействием порыва ветра вместе с упругими колебаниями корпус РН начинает движение как твердое тело. Поперечная составляющая аэродинамической силы, которая при малом угле атаки считается ему пропорциональной, меняет значение за счет смещения корпуса в направлении действия порыва ветра и поворота относительно вектора набегающего потока. Проекция силы тяжести на поперечную ось РН также меняется. Для правильного учета указанных изменений уравнения движения РН необходимо интегрировать с уравнениями, описывающими логику работы автомата стабилизации, что невозможно на ранних этапах проектирования. Учитывая к тому же значительную массу и момент инерции РН, будем считать малыми:

Это позволяет не учитывать влияние указанных выше изменений. Приращение угла атаки (а следовательно, и поперечной составляющей аэродинамической силы) считается зависящим только от величины скорости порыва ветра, заданной в виде функции времени. Таким образом, с учетом принятых допущений аэродинамическая нагрузка в поперечном направлении представляет собой распределенную по длине балки нагрузку с коэффициентом пропорциональности, зависящим от времени. Закон распределения аэродинамической нагрузки по длине РН определяется экспериментально и считается заранее известным. Закон изменения коэффициента пропорциональности (угла атаки) от времени определяется выбором профиля нормативного порыва ветра.

Движение РН, смоделированной в виде упругой балки, можно описать, используя известное уравнение вынужденных поперечных колебаний балки, записанное с учетом гипотезы Фойгта:

(1)

где 𝑚(𝑥) — погонная масса; 𝐵(𝑥) — изгибная жесткость; 𝑞𝑞(𝑥𝑥,𝑡𝑡) — распределенная внешняя нагрузка; ℎ — коэффициент трения.

Данное уравнение должно быть дополнено граничным условием: внутренние усилия в начальном и конечном сечении равны нулю. Значит:

(2)

где 𝐿 — длина РН.

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

𝑞(𝑥,𝑡) = 𝑅(𝑡)𝑌𝑎(𝑥), (3)

где 𝑅(𝑡) — функция, определяющая временную изменчивость аэродинамической силы и изменяющаяся по трапецеидальному закону в соответствии с принятой в настоящей работе моделью порыва ветра; 𝑌𝑎(𝑥) — функция распределения аэродинамической силы по длине РН.

Рассмотрим свободные колебания РН без учета сил трения (при 𝑞(𝑥,𝑡) = 0, ℎ = 0). Подставим метод разделения переменных 𝑦(𝑥,𝑡) = 𝑓(𝑥) ∙ 𝑞(𝑡). В этом случае от уравнения (1) с граничными условиями (2) можно прийти к обыкновенному дифференциальному уравнению

 (4)

с краевыми условиями:

 (5)

Решение (4) с условиями (5) представляет собой классическую задачу Штурма — Лиувилля. Решая ее, можно найти набор собственных форм 𝑓𝑗(𝑥) и собственных частот 𝑝𝑗 рассматриваемой балки (𝑗 = 1, 2, … ). Известно2, что некоторые решения системы (4) соответствуют нулевым собственным частотам. Формы, соответствующие нулевым собственным частотам, определяют поступательное движение РН как твердого тела вместе с центром масс и вращение вокруг центра масс: 𝑓−1 = 1, 𝑓0 = 𝑥 − 𝑥𝐶, где 𝑥𝐶 — координата центра масс РН.

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

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

 (6)

где 𝑦𝐶(𝑡) — перемещение центра масс балки; ϑ(𝑡) — угол поворота оси недеформированной балки; 𝑓𝑗 — форма собственных колебаний балки, соответствующая тону с номером 𝑗; 𝑞𝑗(𝑡) — обобщенная координата, соответствующая тону с номером 𝑗; 𝑁 — число учитываемых упругих тонов.

После подстановки (6) в (1) и применения процедуры Бубнова — Галеркина можно прийти к системеобыкновенных дифференциальных уравнений с постоянными коэффициентами:

Здесь 𝑚 — масса РН; 𝐼 — момент инерции относительно оси, проходящей через центр масс РН перпендикулярно плоскости вращения; 𝑚𝑚𝑗𝑗 — приведенная (обобщенная) масса для -го тона колебаний и определяемая формулой: 

Обобщенные силы в выражении (7) с учетом (3) определяются следующим образом:

(8-10)

Здесь 𝑄𝑎 — максимальное значение главного вектора поперечной аэродинамической нагрузки; 𝑀𝑎 — максимальное значение главного момента поперечной аэродинамической нагрузки, приведенного к центру масс РН; 𝑄0𝑗 — максимальное значение обобщенной силы, соответствующей обобщенной координате 𝑞𝑗.

Первые два уравнения в (7) определяют закон изменения ускорений точек РН в процессе поступательного и вращательного движения РН как твердого тела. Последнее уравнение в (7) определяет закон движения эквивалентной системы осцилляторов.

Рассмотрим движение одного осциллятора под действием трапецеидальной внешней нагрузки, которая:

Для удобства опустим индексы, характеризующие номер тона. Тогда дифференциальное уравнение движения осциллятора с учетом (10) перепишем в виде:

(11)

где 𝑚 — масса осциллятора; 𝑝 — угловая частота собственных колебаний, выраженная в радианах в секунду; n — коэффициент затухания (определяет демпфирующие осциллятора); 𝑄(𝑡) = 𝑄0𝑅(𝑡) — закон изменения внешней нагрузки.

Представим функцию 𝑅𝑅(𝑡𝑡) как совокупность четырех линейных функций:

Соответственно, внешняя нагрузка 𝑄(𝑡) — это совокупность четырех линейных нагрузок 𝑄𝑖(𝑡) = 𝑄0𝑅𝑖(𝑡), (𝑖 = 1, 2, 3, 4). Нагрузка 𝑄1(𝑡) прикладывается с момента 𝑡𝑡 = 0; 𝑄2(𝑡) — с момента 𝑡 = δ; 𝑄3(𝑡) — с момента 𝑡𝑡 = θ + δ; 𝑄4(𝑡) — с момента 𝑡 = θ + 2δ.

Для определения реакции системы на внешнее воздействие разобьем все время действия нагрузки на четыре интервала (рис. 1).

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

здесь 𝑃(𝑡 − ϑ) — смещенный временной закон изменения внешнего воздействия, а 𝑌𝑌(ϑ) характеризует реакцию системы на единичное импульсное воздействие. Реакцию механической системы 𝑞𝑙(𝑡) на линейно возрастающую нагрузку 𝑃(𝑡) = 𝑘t можно выразить через реакцию системы на единичный импульс 𝑌(ϑ) и на внезапно приложенную единичную нагрузку 𝑌1(ϑ):

Реакция на единичный импульс для механической системы с одной степенью свободы, имеющей коэффициент затухания 𝑛𝑛, будет иметь вид4:

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

Реакция на линейно возрастающую нагрузку:

(12)

Примем обозначения: ∆1= 𝑡 − δ, ∆2= ∆1 − θ, ∆3= ∆2 − δ.

Введем функцию, содержащую гармонические члены решения (12):

Обозначим 𝑞ст перемещение под действием статически приложенной к системе силы 𝑄0 = 𝑞ст𝑚𝑝2 и учтем, что 𝑘 = 𝑄0/δ. Суммарное перемещение осциллятора под действием комбинации нагрузок 𝑄𝑖(𝑡) на каждом временном участке будет суммой соответствующих решений (12). Реакция системы, описываемой уравнением (11), на внешнее воздействие трапецеидального профиля будет иметь вид:

(13)

Закон изменения ускорений осциллятора может быть получен двойным дифференцированием повремени выражения (13). Введем функцию . Дифференцируем и запишем результат с учетом номера тона колебаний:

(14)

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

(15)

Для определения ускорений точек РН дважды продифференцируем (6):

(16)

Выразим ускорения обобщенных координат из первых двух уравнений системы (7) с учетом (8) и (9). Представим ускорения точек оси РН, которые она приобретает, двигаясь как твердое тело:

(17)

Кроме того, учтем закон изменения функции 𝑅(𝑡) от времени:

Получим:

(18)

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

(19)

С учетом (16) и (17) закон изменения ускорений точек упругой оси РН под действием порыва ветра трапецеидального профиля будет иметь вид:

(20)

где 𝑎(𝑡) определяется выражением (18), 𝑞𝑗(𝑡𝑡) — выражением (19).

При исследовании упругих колебаний РН демпфирование традиционно принимается на основании данных, полученных по результатам натурных динамических испытаний и представленных в виде значений логарифмических декрементов 𝐷𝑗. Тогда коэффициент, определяющий параметр демпфирования и входящий в выражения (14) и (15), можно вычислить по формуле:

(21)

Входящее в формулу (19) изменение обобщенной координаты 𝑞𝑞𝑗𝑗 под действием статически приложенной обобщенной силы 𝑄0𝑗 определяется выражением:

(22)

Принимая во внимание (14), (15), (18), (19), (21), (22), формула (20) — это аналитическое выражение, определяющее функции изменения ускорений по времени для точек оси РН под действием внешней поперечной аэродинамической силы, изменяющейся по трапецеидальному закону с учетом влияния диссипативных сил.

Зная закон изменения ускорений, можно по известным методикам [4] определить динамические, а затем и суммарные корпусные нагрузки, действующие в сечениях РН. Вычислим изгибающий момент, обусловленный инерционными силами, возникающими за счет упругих колебаний корпуса РН:

(23)

Здесь 𝑀𝑗(𝑥) — функция распределения по длине РН единичного (при ускорении 𝑞𝑗(𝑡), равном единице) изгибающего момента для -го тона колебаний. Ее можно найти по формуле:

(24)

Как отмечалось выше, для проведения модального анализа в данной работе использовался метод КЭ. Данный подход обусловлен тем, что на практике у динамической модели РН достаточно сложная структура. Она включает подконструкции, и нельзя пренебрегать их собственной динамикой. Подконструкции могут крепиться к корпусу РН в одном сечении либо располагаться параллельно продольной оси РН и иметь несколько точек крепления. В этом случае расчет динамических характеристик в континуальной постановке является сложной математической задачей. Кроме того, динамические модели отдельных подконструкций представляются предприятиями-разработчиками в конденсированном (матричном) виде в формате Nastran. По этой причине для расчета динамических характеристик конструкции будет оптимальным задействовать пакет программ инженерного анализа Nastran.

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

Для сложной конструкции РН уравнения (7) сохраняют свой вид [3]. Стандартная информация вывода программы Nastran может быть базой для получения параметров эквивалентной системы осцилляторов, массы и момента инерции РН, а также для вычисления обобщенных сил, входящих в третье уравнение системы (7). Для формирования левой части третьего уравнения в (7) необходимы значения собственных частот 𝑝𝑗 (Radians) и обобщенных масс 𝑚𝑗 (Generalized mass). Для определения обобщенных сил в правой части третьего уравнения в (7) необходимы функции собственных форм колебаний 𝑓𝑗(𝑥) (Eigenvector). При вычислении динамических инерционных нагрузок по (23) вместо (24) удобнее воспользоваться единичными инерционными нагрузками (силами и моментами), которые выводятся программой Nastran после стандартного запроса усилий при расчете собственных форм и частот. Единичные инерционные нагрузки выводятся отдельно для каждого тона колебаний, причем помноженные на квадрат собственной частоты, что следует учитывать для их корректного использования.

Результаты исследования. Для проведения тестовых расчетов рассмотрена ракета-носитель среднего класса тандемной схемы. Вид функций погонной изгибной жесткости 𝐵𝐵(𝑥𝑥) и массы 𝑚𝑚(𝑥𝑥), а также распределения сосредоточенных масс 𝑚соср(𝑥) = 𝑚соср𝑟∆(𝑥 − 𝑥𝑟) по длине рассматриваемой РН представлен на рис. 2 (через ∆ обозначена дельта-функция Дирака).

Для проверки полученных аналитических решений при воздействии трапецеидальной внешней нагрузки рассчитаны динамические ускорения РН по формуле (20) с учетом (14), (15), (18), (19), (21), (22). При этом учтены 5 упругих тонов собственных поперечных колебаний РН. Характеристики 𝑝𝑗, 𝑚𝑗 и 𝑓𝑗(𝑥) были получены по результатам расчета в программном комплексе MSC Nastran с использованием последовательности решения для модального анализа собственных колебаний (SOL 103). Динамическая конечно-элементная модель РН представлена в виде набора балочных элементов с различными инерционными и жесткостными характеристиками. К балочным упруго или жестко присоединены:

Конечноэлементная модель РН включает около 1000 составляющих. Элементы, моделирующие закрепления, не использовались для сохранения способности движения РН как твердого тела.

По той же конечноэлементной модели РН в программном комплексе MSC Nastran рассчитаны ускорения с использованием последовательности решения SOL 119, применяемой для модального анализа переходных процессов. При этом в модальном разложении были учтены все тона собственных колебаний в диапазоне до 100 Гц. Аэродинамическая нагрузка представлена поперечными погонными нагрузками, распределенными по всем балочным элементам, моделирующим корпус РН. Кроме того, задан трапецеидальный закон изменения аэродинамической нагрузки по времени.

На рис. 3 представлены результаты сравнительного анализа полученных с помощью двух различных подходов ускорений 𝑎0(𝑡) = 𝑦(𝑡, 𝑥0) некоторой точки оси РН с координатой 𝑥𝑥 = 𝑥𝑥0. Видно, что хорошо согласуются два решения:

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

На основе значений динамических ускорений получены зависимости (рис. 4). Они показывают, как изменение пиковых значений изгибающего момента 𝑀 для различных сечений РН связано с параметром θ, характеризующим продолжительность действия порыва. На рис. 4 значения изгибающего момента представлены в виде безразмерных величин 𝑀∗. Они рассчитаны делением размерного изгибающего момента на максимальное значение для данного сечения (например для сечения 𝑥 = 0,3𝐿 максимальное значение 𝑀𝑚ax = 4,6 ∙ 105 Н ∙ м), найденное при варьировании величины θ во всем рассматриваемом диапазоне.

Итак, продолжительность действия порыва на графиках представлена безразмерной величиной θ∗, полученной путем деления параметра θ на период 𝑇1 первого тона колебаний РН. Каждая линия графика соответствует одному сечению РН. Все сечения РН сгруппированы по характеру функции 𝑀∗(θ∗) и приведены на различных графиках, а длина РН соответствующим образом разделена на зоны (рис. 5).

Из рис. 4 видно, что для первой зоны максимум изгибающего момента достигается уже при значениях θ не менее 15 % от периода 𝑇𝑇1 первого тона колебаний. Для второй и третьей зон — несколько больше половины периода 𝑇𝑇1. Дальнейшее увеличение значений θ не влияет на величины максимальных значений изгибающего момента. Для зон 4–6 максимум значения изгибающего момента оказывается локальным и располагается в зоне значений θ, близких к величине половины периода 𝑇1 первого тона колебаний РН. Полученные результаты полностью соответствуют итогам конечноэлементного моделирования, проведенного ранее6.

Таким образом, для получения максимальных значений изгибающего момента в сечениях рассматриваемой РН необходим динамический анализ поведения конструкции с внешним воздействием в виде порыва ветра, продолжительность которого определяется параметром θ, близким по значению к половине периода первого тона колебаний РН.

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

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

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

1. Лиходед А. И. Динамика конструкций и определение нагрузок. Королев : Изд-во АО ЦНИИмаш, 2020. 239 с.2. Колесников К. С. Динамика ракет. М. : Машиностроение, 2003. 520 с.3. Бидерман В. Л. Теория механических колебаний. М. : URSS, 2017. 416 с.4. Яблонский А. А. Курс теоретической механики. М. : Интеграл-Пресс. 2007. 603 с.5. Малыхина О. И. Автоматизация обработки результатов конечноэлементного анализа нагружения конструкций ракетно-космической техники : сб. ст. VII науч.-тех. конф. молодых ученых и специалистов Центра управления полетами. Королев : ЦНИИМаш, 2017. С. 427– 434.6. Малыхина О. И., Глуговский М. С. Анализ влияния профиля порыва ветра на величину корпусных нагрузок ракеты-носителя в полетных случаях нагружения // Актуальные проблемы ракетно-космической техники (VI Козловские чтения) : мат-лы VI всерос. науч.-тех. конф. Самара, 2019. Т. 1. С. 133–138.