МОДЕЛИРОВАНИЕ ДИНАМИЧЕСКОЙ УСТОЙЧИВОСТИ ТОНКОСТЕННЫХ КОНСТРУКЦИЙ, ЧАСТИЧНО ЗАПОЛНЕННЫХ ЖИДКОСТЬЮ, ПРИ ГИДРОСТАТИЧЕСКОМ ВОЗДЕЙСТВИИ Текст научной статьи по специальности «Физика»

УДК 534.131.2

DOI: 10.18698/2309-3684-2022-3-317

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

© Сонги Пак, В.Г. Григорьев Московский авиационный институт, Москва, 125993, Россия

В данной работе рассматривается задача устойчивости тонкостенной оболочечной конструкции с двумя полусферическими днищами одинаковой толщины, частично заполненной жидкостью, которая погружена во внешнюю жидкую среду и находится под действием гидростатического давления. Получены динамические характеристики такой конструкции, содержащей ограниченный объем жидкости, находящейся под внутренним давлением и гидростатическим давлением. Разработанная программа для расчета динамических характеристик осесимметричных оболочечных конструкций, содержащих жидкость, основана на методе конечных элементов. Конечные элементы имеют кольцевую форму при вращении вокруг оси симметрии. Программа реализована в среде табличного процессора Excel на языке Visual Basic for Applications (VBA). Она позволяет рассчитывать собственные частоты тонкостенных упругих конструкций, взаимодействующих с произвольным количеством ограниченных объемов жидкостей с учетом влияния статического напряженно-деформированного состояния, вызванного гидростатическими внутренним давлением и прочими внешними силами, не нарушающими осевую симметрию задачи.

При фиксированном значении внутреннего давления выполняется расчет низших собственных частот колебаний с различными числами волн по окружности. Последовательным уточнением определяется критическая толщина оболочки, при которой хотя бы одна из собственных частот достигает нулевого значения. Внутреннее давление р изменяется от 0 до 1 атм. с шагом 0,1 атм. и расчеты повторяются для получения каждого критического значения. При каждом значении давления на графике «количество волн — собственная частота» строятся кривые. Построены графики, иллюстрирующие процесс получения этого критического значения. На координатной плоскости «давление — толщина оболочки» построена граница области неустойчивости

Ключевые слова: динамические характеристики, собственные частоты, осесим-метричная оболочка, метод конечных элементов, подводный аппарат, упругая оболочка с жидкостью

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

Excel на языке Visual Basic for Applications (VBA). Программа основана на использовании метода конечных элементов, теоретические основы ее содержатся в работах [1-4].

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

Модель динамической устойчивости конструкции в трехмерной постановке. Рассматривается общий принцип построения конечных элементов для расчета собственных колебаний упругих систем с учетом влияния статического напряженно-деформированного состояния посредством нелинейности деформаций.

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

S. + S

G

Рис. 1. Конструкция и обозначения

На рис. 1. обозначим Q и Q0 — объемы, занятые упругим телом и жидкостью, £ — поверхность упругого тела, £0 — поверхность полости, смоченная жидкостью объема Q0, £ — поверхность полости, подверженная действию избыточного внутреннего давления

X

X

X

2

газов , £ и £ст — участки поверхности, на которых заданы кинематические и динамические граничные условия. Область над поверхностью жидкости не именуется, поскольку динамика газов, создающих избыточное статическое давление, в рассматриваемом низкочастотном диапазоне не учитывается. Границу жидкой массы составляют свободная поверхность £ и контактная поверхность .

Введены вектор й и II — смещения точек упругого тела и жидкости, С — вектор ускорения гравитационных сил, п0 — вектор нормали к поверхности £.

Конструкция находится в однородном гравитационном поле с

вектором ускорения свободного падения С, который ортогонален свободной поверхности жидкости в невозмущенном состоянии. В соответствии с формулировкой работы [3] поле малых перемещений частиц жидкости представлено градиентом скалярного поля:

и^Ф в (30, (1)

где Ф — потенциал смещений. Граничное условие на смоченной

поверхности имеет вид условия непротекания:

¿ф

— = (и,п0) на (2)

ан

Из условия несжимаемости следует уравнение Лапласа для потенциала смещений:

ДФ = 0 в Q0. (3)

На свободной поверхности жидкости должны быть выполнены условия:

аФ

— = 1 на ^ (4)

ап

Ф + Ог] = 0 на Е, (5)

где Т] — нормальные смещения точек свободной поверхности,

О =

С

Для описания движения деформируемого упругого тела используем нелинейные геометрические соотношения [5-8], записанные в виде:

£(П)=^(й) + £Ж(й,й), (6)

где £1 — линейный дифференциальный оператор, £ж — билинейный

симметричный дифференциальный оператор, а компоненты деформаций представлены векторами-столбцами.

Физическое соотношение полагаем линейным:

а - Бе, (7)

где Б — матрица модулей упругости, а а — вектор компонент тензора напряжений. В настоящей работе упругое тело представлено тонкой упругой оболочкой [6], что не меняет формальную запись (6)-(7).

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

иполн\1) = и + ). (8)

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

£(ио +и) = £^(ио) + г^(и) + £ж(ио,ио) + +2еыь (ио ,и) + еж (и,и).

Выражение для полной потенциальной энергии деформации имеет вид:

й0 +й)тВг(й0 + й)с!Г, (10)

2 е

где е — область пространства, занятая упругим телом, которое в данной работе является упругой оболочкой.

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

Гв = ± | ([(е£ (и) + 2гж (ио, и))]' Б [(е£ (и) + 2гмь (ио, и))"

е

(11)

+2г(ио )т Вгж (и, и))ёУ.

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

перемещении точек упругого тела и, скалярное поле потенциала смещений частиц жидкости Ф и поле вертикальных перемещений частиц на поверхности жидкости Т].

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

V = VQ + ¥,+ Ув + .

(12)

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

VЕ =

1 ¡р0О]2 я.

(13)

Потенциал сил контактного взаимодействия тела с жидкостью определяется формулой:

Зй

Зй

г„ х

--гя х-

др р да

Л

(14)

- р0 (б, й) (й, га X гр )} ёаёр,

где а и /3 — координаты на контактной поверхности (не обязательно ортогональные), г = г (а,/?) — радиус-вектор точки на поверхности и касательные к координатным линиям векторы га = дг/да, Тр = Зг/д/З. Здесь гидростатическое давление на смоченную поверхность

=/'о(С,Г-Х0),

где х0 — произвольная точка свободной поверхности жидкости. Потенциал давления газа равен

V' -11

(

и-

Зй

Зи

ГаХ--ГЯ Х-

а др р да

Л

ёаёр,

(15)

где область интегрирования включает смоченную жидкостью поверхность в соответствии с законом Паскаля.

Кинетическая энергия гидроупругой системы может быть вычислена по формуле:

к

0

т = I\ptidV+11 р0 (уф)2 ¿Г. (16)

2 <2 2<Зо

Однако при попытке использовать функционал Лагранжа в виде Ь — Т - V для конечно-элементного моделирования собственных колебаний конструкции и вычисления ее собственных частот возникнут трудности с обеспечением выполнения кинематических ограничений, налагаемых уравнениями (2), (3) и (4).

В работах [1-3] предложена модификация лагранжиана, в которой изменено выражение кинетической энергии в виде:

о Оо (17)

/?0Ф(й • п0)й?5'+| р0 Ф

&0 2

Выражение (17) в точности равно кинетической энергии конструкции, если присутствующие в нем поля удовлетворяют указанным ограничениям. В то же время, численные решения, полученные с использованием модифицированной формулы лагранжиана

Ь - Т*- V, (18)

автоматически этим ограничениям удовлетворяют.

Осевая симметрия конструкции ведет к важному следствию. А именно, спектр конструкции распадается на множество подспектров по признаку числа волн смещений на окружностях горизонтальных ее сечений [9]. Сказанное означает, что в цилиндрической системе координат От^г, ось которой совпадает с осью симметрии конструкции, любая собственная форма колебаний может быть представлена в виде (условно нумеруем ее п):

(19)

где перемещения u, v, w представляют радиальное, окружное и осевое направления, параметр mn — представляет число волн в окружном направлении, а перед тригонометрическими функциями стоят амплитудные значения соответствующих величин. Особо надо отметить роль параметра гя, который может принимать значения 0 или

1. В случае неосесимметричной формы колебаний (mn > 0), ), когда

каждая собственная частота имеет кратность 2, он позволяет различать две линейно независимые собственные формы, соответствующие одной собственной частоте. В случае осесимметричных колебаний (mn = 0) очевидным образом значение zn = 0 соответствует продольно-радиальным, а /и = 1 крутильным колебаниям конструкции.

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

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

u(x) = N(x)ö, (20)

где 5 — вектор-столбец дискретных узловых параметров элемента, а N(x) — матрица функций формы, которую можно представить в виде совокупности столбцов (к — размерность вектора 5 ):

N(x) = [iij (х) ... п,(х)]. (21)

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

1 т

V = - 5TK, 5, (22)

где K — матрица жесткостей конечного элемента, вычисляемая по формуле:

K = j[B(x)T DB(x) + S(x) ] dV, (23)

Qe

в которой B(x) получается из первого слагаемого подинтегрального выражения в формуле (11):

£L (и) + 2вЖ (ио, и) -

(24)

= [е'(п1) + 2смь(и0,п1) ... в^К) + 2вж(и0,п,)]б^В(х)б, а компоненты матрицы 8 (г) получены из второго слагаемого там же:

(х) - 2£(ио)тD£NL(п,п7). (25)

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

формуле, аналогичной (20).

В формуле (23) учтено наличие статического напряженно-деформированного состояния, и ее можно использовать для расчета малых колебаний в окрестности этого состояния. На выражение для матрицы масс наличие статического напряженно-деформированного состояния влияния не оказывает. Матрица масс получается на основе выражения (17), в котором кроме аппроксимации перемещений аналогично (20) добавляются члены, аппроксимирующие значения потенциала смещений и вертикальных перемещений точек свободной поверхности жидкости. Очевидно, что кроме элементов упругого тела в этом выражении должны участвовать элементы жидкости и элементы свободной поверхности.

Таким образом, после объединения матриц отдельных элементов получаем задачу определения собственных частот и форм колебаний конструкции с матрицей масс М и матрицей жесткостей К в виде:

(К -ЯМ)б - 0, (26)

Собственные значения задачи (26) образуют дискретный ряд:

Я<Я<Я< •••

и, соответственно им, вычисляются собственные вектора:

бр б2, йв —

представляющие собственные формы колебаний.

Угловая собственная частота с определяется как корень квадратный из Я, и тогда справедлива формула

(К -с2М) б - 0, (27)

а собственная частота колебаний в Гц /к определяется из формулы

С - 2ж/к. (28)

Изменение матрицы жесткостей при учете статического напряженно-деформированного состояния конструкции может привести к тому, что один или несколько членов последовательности Хк станут отрицательными. В таком случае статическое состояние конструкции оказывается неустойчивым, и каждый из векторов àk

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

Поскольку уравнения собственных колебаний осесимметричной конструкции распадаются на группы с заданным числом волн по окружности m , то для каждого целого значения m матрицы K и M вычисляются отдельно. Соответственно, для каждого значения m образуется отдельный ряд возрастающих собственных значений Xk( m ), к = 1,2,.... В результате выбор наименьшего собственного

значения 2min можно осуществить только путем перебора m в диапазоне от 0 до определенной величины, начиная с которой наименьшие (для данного m ) собственные значения начинают монотонно возрастать. При этом необходимо для каждого значения m вычислять соответствующие матрицы масс и жесткостей и с их использованием вычислять наименьшее собственное значение \ (m) .

Это заметно усложняет задачу вычисления 2min . Однако это усложнение несущественно по сравнению с задачей трехмерного моделирования таких конструкций.

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

Результаты расчета. Рассмотрим оболочечную конструкцию с полусферическими днищами и толщиной h , заполненную жидкостью и погруженную в воду (рис. 2). Радиус оболочки R = 1 м, высота H = 1 м, характеристики материала: E = 2-1011 Н/м2, v = 0,3 , р = 7800 кг/м3. Плотность жидкости внутри оболочки р1 = 800 кг/м3 и глубина жидкости H = 2 м. Цилиндр заполнен водой вне оболочки высотой H2 = 15 м.

Расчеты были выполнены при различных величинах внутреннего давления оболочки от 0 до 1 атм с шагом 0,1 атм. Результаты расчетов для исследованных значений давления представлены в виде зависимостей собственных частот колебаний от числа волн m при различных толщинах оболочки, включая критическое значение толщины.

/ (m ) 80 70 60 50 40 30 20 10 0 -10 -20 -30

10

-т--т-—«-

Рис. 2. 3D и конечно-элементная модель конструкции

к- = 3,5 5 мм

п = 3,1542 мм

0 12 3 4 1 5 5 ( 5 7 8 9 10 11 12 13 1 4 15

к = = 3 мм

Рис. 3. График зависимости собственной частоты f (ж) от количества волн ж при внутреннем давлении р = 0

На рис. 3 показаны значения собственных частот f (ж) при внутреннем давлении р = 0, где f (ж) — низшая собственная частота для собственных форм с ж волнами. На рис. 4 показаны соответствующие им собственные значения \ (ж) уравнения (26). При этом, если

некоторое собственное значение \ < 0, то на графике рис. 3 откладывается величина а = . В этом случае однородное уравнение собственных колебаний имеет два независимых решения Е,в~аа и Е,в+а, из которых первое реализует форму потери устойчивости.

Аналогично результаты получены для всех исследованных величин внутреннего давления до 1 атм и построен график критических значений толщины оболочки (рис. 5). Этот график представляет собой границу области устойчивости на плоскости р - Н .

Л1 (/и) 1-Ю5

1,5 105

1,0 105

0,5105

-0,5 105

Л, мм 3,16 3,15 3,14 3,13 3,12 3,11 3,10

/1

\\ \\ V V Н = Н = : : 3,5 мм 3,1542 мм

0 1 2 3 5 4 5 ; 6 1 8 3 9 10 1 Н = 3 11 мм 1. 3 14

Рис. 4. График зависимости собственного значения \ (т) от количества волн т при внутреннем давлении р = 0

О бласть устойчивости

Область неустойчивости

0 0,1 0,2 0,3 0,4 0,5 0,6 0,7 0,8 0,9 1,0 р, атм

Рис. 5. График зависимости критической толщины оболочки Н от внутреннего давления р

0

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

Построена область устойчивости на координатной плоскости двух параметров: p — внутреннее давление и h — толщина оболочки. Часть плоскости ниже границы соответствует область неустойчивости.

ЛИТЕРАТУРА

[1] Григорьев В.Г. Применение метода конечных элементов к расчету колебаний упругих оболочечных конструкций, содержащих жидкость. Труды научного семинара ТГУ. Динамика упругих и твердых тел, взаимодействующих с жидкостью. Томск, ТГУ, 1978, с. 55-60.

[2] Григорьев В.Г. Расчет динамических характеристик сложных оболочечных конструкций с жидкостью. Колебания упругих конструкций с жидкостью. Сборник научных докладов IV симпозиума. Москва, ЦНТИ "Волна", 1980, с. 102-107.

[3] Григорьев В.Г., Григорьева Е.В. Контактное взаимодействие ограниченного объема жидкости с деформируемым твердым телом под влиянием гравитационных сил. Известия Российской академии наук. Механика твердого тела, 2011, № 2, с. 147-159.

[4] Пак С., Григорьев В.Г. Устойчивость тонкостенных осесимметричных соос-ных конструкций, содержащих жидкость, при многофакторных нагрузках. Труды МАИ, 2021, № 119. DOI: 10.34759/trd-2021-119-08

[5] Новожилов В.В. Основы нелинейной теории упругости. Ленинград; Москва, Гостехиздат, 1948, 211 с.

[6] Кармишин А.В., Лясковец В.А., Мяченков В.И., Фролов А.Н. Статика и динамика тонкостенных оболочечных конструкций. Москва, Машиностроение, 1975, 375 с.

[7] Димитриенко Ю.И. Механика сплошной среды. Т. 4. Основы механики твердых сред. Москва, Изд-во МГТУ им. Н.Э. Баумана, 2013, 624 с.

[8] Димитриенко Ю.И. Обобщенная трехмерная теория устойчивости упругих тел. Часть 1: конечные деформации. Вестник МГТУ им. Н.Э. Баумана. Сер. Естественные науки, 2013, № 4, с. 79-95.

[9] Григорьев В.Г. Методология исследования динамических свойств сложных упругих и гидроупругих систем. Дисс. д-ра техн. наук. Москва, 2000, 328 с.

[10] Голованов А. П., Тюленева О. Н., Шигабутдинов А. Ф. Метод конечных элементов в статике и динамике тонкостенных конструкций. Москва, Физматлит, 2006, 392 с.

[11] Агапов В.П. Метод конечных элементов в статике, динамике и устойчивости. Москва, АСВ, 2000, 248 с.

[12] Кирилин А.Н., Ахметов Р.Н., Соллогуб А.В. Проектирование, динамика и устойчивость движения ракет-носителей: методы, модели, алгоритмы, программы в среде MathCad. Москва, Машиностроение, Машиностроение-Полет, 2013, 294 с.

[13] Фиалко С.Ю. Прямые методы решения систем линейных уравнений в современныхМКЭ-комплексах. Москва, Изд-во СКАД СОФТ, 2009, 159 с.

[14] Zienkiewicz O.C., Taylor R.L. The finite element method: its basis and fundamentals. Butterworth-Heinemann, Elsevier, 2013, 756 p.

[15] Zienkiewicz O.C., Taylor R.L., Fox D.D. The finite element method for solid and structural mechanics. Butterworth-Heinemann, Elsevier, 2014, 624 p.

[16] Zienkiewicz O.C., Taylor R.L., Nithiarasu P. The finite element method for fluid dynamics. Butterworth-Heinemann, Elsevier, 2014, 544 p.

Статья поступила в редакцию 31.05.2022

Ссылку на эту статью просим оформлять следующим образом: Пак Сонги, Григорьев В.Г. Моделирование динамической устойчивости тонкостенных конструкций, частично заполненных жидкостью, при гидростатическом воздействии. Математическое моделирование и численные методы, 2022, № 3, с. 3-17.

Пак Сонги — аспирантка кафедры «Проектирование и прочность авиационно-ра-кетных и космических изделий» Московского авиационного института. e-mail: haaz7790@naver.com

Григорьев Валерий Георгиевич — д-р. техн. наук, профессор кафедры «Проектирование и прочность авиационно-ракетных и космических изделий» Московского авиационного института. e-mail: valgrigg@rambler.ru

Modeling of dynamic stability of thin-walled structures partially filled with liquid under hydrostatic action

© SongYi Park, V.G. Grigor'ev Moscow Aviation Institute (National Research University), Moscow, 125993, Russia

In this paper, we consider the problem of stability of a thin-walled shell structure with two hemispherical bottoms of the same thickness, partially filled with liquid, which is immersed in an external liquid medium and is under hydrostatic pressure. The dynamic characteristics of such a structure containing a limited volume of liquid under internal pressure and hydrostatic pressure are obtained. The developed program for calculating the dynamic characteristics of axisymmetric shell structures containing liquid is based on the finite element method. The finite elements have an annular shape when rotated around the axis of symmetry. The program is implemented in Excel spreadsheet using Visual Basic for Applications (VBA). It allows to calculate the natural frequencies of thin-walled elastic structures interacting with an arbitrary number of liquids, considering the influence of the static stress-strain state caused by hydrostatic and internal pressure and other external forces that do not violate the axial symmetry of the problem.

At a fixed value of the internal pressure, the calculation of the lowest natural frequencies of vibrations with different numbers of waves along the circumference is performed. By successive refinement, the critical thickness of the shell is determined, at which at least one of the natural frequencies reaches zero. The internal pressure p varies from 0 to 1 atm. in increments of 0,1 atm. and the calculations are repeated to obtain each critical value. At each pressure value, curves are plotted on the graph "number of waves — natural frequency". On the coordinate plane "pressure — shell thickness" the boundary of the instability region is constructed.

Keywords: dynamic characteristics, natural frequency, thin-walled axial symmetric structure, finite element method, submersible vehicle, elastic shell with liquid

REFERENCES

[1] Grigoriev V. G. Primenenie metoda konechnyh elementov k raschetu kolebanij uprugih obolochechnyh konstrukcij, soderzhashchih zhidkost' [Application of the finite element method to the calculation of vibrations of elastic shell structures containing liquid]. Trudy nauchnogo seminara TGU. Dinamika uprugih i tverdyh tel, vzaimodejstvuyushchih s zhidkost'yu [Proceedings of the TSU scientific seminar. Dynamics of elastic and solid bodies interacting with a liquid]. Tomsk, TSU, 1978, pp. 55-60.

[2] Grigoriev V.G. Raschet dinamicheskih harakteristik slozhnyh obolochechnyh konstrukcij s zhidkost'yu [Calculation of dynamic characteristics of complex shell structures with liquid]. Kolebaniya uprugih konstrukcij s zhidkost'yu. Sbornik nauchnyh dokladov IV simpoziuma [Vibrations of elastic structures with liquid. Collection of scientific reports of the IV Symposium]. Moscow, Central Research Institute "Volna", 1980, pp. 102-107.

[3] Grigor'ev V.G., Grigor'eva E.V. Contact interaction of a bounded liquid volume with a solid under the action of gravity forces. Mechanics of Solids, 2011, vol. 46, no. 2, pp. 280-290.

[4] Park S., Grigoriev V. Stability of thin-walled axisymmetric coaxial structures containing liquid under multifactor loads. Trudy MAI, 2021, no. 119. DOI: 10.34759/trd-2021-119-08

[5] Novozhilov V.V. Osnovy nelinejnoj teorii uprugosti [Fundamentals of the nonlinear theory of elasticity]. Leningrad; Moscow, Gostekhizdat Publ., 1948, 211 p.

[6] Karmishin A.V., Lyaskovets V.A., Myachenkov V.I., Frolov A.N. Statika i dinamika tonkostennyh obolochechnyh konstrukcij [Statics and dynamics of thin-walled shell structures]. Moscow, Mashinostroenie Publ., 1975, 375 p.

[7] Dimitrienko Yu.I. Mekhanika sploshnoy sredy. Tom 4. Osnovy mekhaniki tverdogo tela [Continuum Mechanics. Vol. 4. Fundamentals of solid mechanics]. Moscow, BMSTU Publ., 2013, 624 p.

[8] Dimitrienko Yu.I. Generalized three-dimensional theory of elastic bodies stability. Part 1: finite deformations. Herald of the Bauman Moscow State Technical University. Series Natural Sciences, 2013, no. 4, pp. 79-95.

[9] Grigoriev V.G. Metodologiya issledovaniya dinamicheskih svojstv slozhnyh uprugih i gidrouprugih sistem [Methodology for the study of dynamic properties of complex elastic and hydroelastic systems]. Diss. Dr. Sci. (Eng.), Moscow, 2019, 328 p.

[10] Golovanov A.P., Tyuleneva O.N., Shigabutdinov A.F. Metod konechnyh ele-mentov v statike i dinamike tonkostennyh konstrukcij [Finite element method in the statics and dynamics of thin-walled structures]. Moscow, Fizmatlit Publ., 2006, 392 p.

[11] Agapov V.P. Metod konechnyh elementov v statike, dinamike i ustojchivosti [Finite element method in statics, dynamics and stability]. Moscow, ASV Publ., 2000, 248 p.

[12] Kirilin A.N., Akhmetov R.N., Sollogub A.V. Proektirovanie, dinamika i ustojchivost' dvizheniya raket-nositelej: metody, modeli, algoritmy, programmy v srede MathCad [Design, dynamics and stability of the motion of launch vehicles: methods, models, algorithms, programs in the MathCad environment]. Moscow, Mashinostroenie Publ., Mashinostroenie-Polet Publ., 2013, 294 p.

[13] Fialko S.Yu. Pryamye metody resheniya sistem linejnyh uravnenij v sovremennyh MKE-kompleksah [Direct methods for solving systems of linear equations in modern FEM complexes]. Moscow, SCUD SOFT Publ., 2009, 159 p.

[14] Zienkiewicz O.C., Taylor R.L. Zhu J.Z. The finite element method: its basis and fundamentals. Butterworth-Heinemann, Elsevier, 2013, 756 p.

[15] Zienkiewicz O.C., Taylor R.L., Fox D.D. The finite element method for solid and structural mechanics. Butterworth-Heinemann, Elsevier, 2014, 624 p.

[16] Zienkiewicz O.C., Taylor R.L., Nithiarasu P. The finite element method for fluid dynamics. Butterworth-Heinemann, Elsevier, 2014, 544 p.

Park S., Postgraduate student, Department of Design and Strength of Aeronautical Missile and Space Products, Moscow Aviation Institute (National Research University). e-mail: haaz7790@naver.com

Grigor'ev V.G., Dr. Sc. (Eng.), Professor, Department of Design and Strength of Aeronautical Missile and Space Products, Moscow Aviation Institute (National Research University). e-mail: valgrigg@rambler.ru