Методика расчета осесимметричных колебаний оболочек вращения с жидкостью по дифференциальной модели Текст научной статьи по специальности «Физика»

УДК 539.3

DOI: 10.18698/2308-6033-2024-2-2338

Методика расчета осесимметричных колебаний оболочек вращения с жидкостью по дифференциальной модели

© К.М. Нгуен1, Д Р. Шелевая1,2, Д.А. Красноруцкий1,3

1НГТУ, Новосибирск, 630073, Российская Федерация 2ИГиЛ СО РАН, Новосибирск, 630090, Российская Федерация 3СибНИА им. С. А. Чаплыгина, Новосибирск, 630051, Российская Федерация

Предложена методика расчета динамического деформирования осесимметричных ортотропных оболочек вращения на основе решения дифференциальных уравнений методом конечных разностей. В частности, для решения сопряженной задачи о малых колебаниях бака с идеальной несжимаемой жидкостью для генерирования локальных весовых коэффициентов конечных разностей на произвольном трафарете узловых точек расчетной области используются полигармонические радиальные базисные функции. Деформирование меридиана осесимметричной оболочки описывается системой шести дифференциальных уравнений движения, полученных на основе общих пространственных уравнений для разрешающих функций в глобальной системе координат. Для описания поворотов используется вектор конечного поворота (вектор Эйлера). Уравнения не содержат начальной кривизны, описывают большие повороты, перемещения и деформации, учитывают утонение/утолщение и влияние поперечного сдвига при деформировании для толстой оболочки. Предложенная методика реализована в собственном программном комплексе DARSYS. Рассчитанные по предложенной методике частоты колебаний баков с жидкостью сравниваются с частотами, полученными другими методами. Анализ сходимости выполняется с помощью комплекса ANSYS, а также программы, реализующей метод конечных и граничных элементов.

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

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

Для определения динамических характеристик тонкостенных конструкций с жидкостью проводят экспериментальные исследования [2-4]. При выполнении расчетов используются разные подходы: точные решения для частных случаев геометрии и граничных условий [1, 5]; приближенные методы оценки на основе усеченного разложения решения по некоторому базису [3, 6]; интегральные подходы на основе минимизации функционалов энергий, в том числе варианты метода конечных элементов (МКЭ) [7-12], метода граничных элементов (МГЭ) [13, 14], и др. Наиболее популярный и универсальный метод анализа — МКЭ, реализованный во множестве программных продуктов и постоянно развиваемый учеными в узких направлениях.

Цель данной работы — программная реализация и совершенствование альтернативного МКЭ и МГЭ подхода к расчету тонкостенных оболочечных конструкций, взаимодействующих с жидкостью. В основе предлагаемой методики лежит решение дифференциальных уравнений с помощью метода конечных разностей (МКР). В нем весовые коэффициенты разложения производных на произвольном трафарете узлов вычисляются с помощью сплайн-интерполяций на основе радиальных базисных функций. К идеологически наиболее близким исследованиям в данном направлении можно отнести представленные в работах [15-17], в которых применяется метод дифференциальных квадратур.

Радиальные базисные функции (РБФ) [18-24] используются для построения интерполирующих многочленов, которые в свою очередь применяются при решении краевых задач для дифференциальных уравнений в разных областях математической физики. Можно выделить два основных направления применения РБФ:

бессеточные методы — линейные комбинации РБФ используются для аппроксимации разрешающих функций задачи. Неизвестными параметрами являются коэффициенты разложения разрешающих функций в ряды по РБФ. Первые определяются из решения СЛАУ, полученной после подстановки в дифференциальные уравнения и краевые/граничные условия соответствующих многочленов на основе РБФ;

метод конечных разностей — интерполяция по РБФ используется для вычисления локальных весовых коэффициентов для разностных трафаретов, аппроксимирующих дифференциальный оператор и краевые/граничные условия по соседним узловым точкам. При этом подлежащими определению неизвестными являются узловые значения разрешающих функций, как и в традиционном МКР. Такая разновидность МКР в англоязычной литературе имеет название Radial Basis Function Finite Difference method (RBF-FD). Этот подход и принят в данной работе.

Методика расчета осесимметричных колебаний оболочек вращения с жидкостью...

Общие уравнения ортотропной оболочки. Рассмотрим бесконечно малый элемент dsl х ^ криволинейной срединной поверхности оболочки, уравнение которой задается радиусом-вектором (рис. 1):

Г ( ¿2) = Х ( ¿2)) + х2 (^ ^)/2 + Хз (^ ^ )7з = Хк (^ ^)Iк, (1)

где ¡к — базисные векторы глобальной системы координат; ¿1, ¿2 — криволинейные координаты, по повторяющимся индексам ведется суммирование от 1 до з.

Рис. 1. Перемещение и деформирование малого элемента оболочки: и — вектор перемещения точки О в О*

В точке О (см. рис. 1) находится локальная система координат, базисные векторы которой выражаются следующим образом:

е1 =

д г

д

е2 =

дг

д

е3 = е1 х е2.

(2)

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

е1 "Р11 Р12 Р13 " ¿1

е2 > = Р21 Р22 Р23 < е (3)

Л, _Р31 Р32 Р33 _ ¿3

Деформация малого элемента (см. рис. 1) состоит:

- из перемещения параллельно самому себе как жесткого целого

*

и поворота относительно точки О ;

- изменения длин его сторон;

- изменения угла между сторонами.

Поворот малого элемента представим таким образом:

h ^11 ^12 ^13

h* h = ^ 21 ^ 22 ^ 23

h* s3 ^32 ^33

1 - cos ш

К =1 Tw (ш? +ш1

ш

)Ч

° sn = ^nksk;

sin ш 1 - cos ш

-шк +-

(4)

ш

ш

2 ш1ш j,

=-

Sin ш

-шк +-

1 - cos ш "2

(5)

ш1 ш j, i = 1, 2, 3; j = 2, 3,1; к = 3,1, 2,

ш ш"

где X — компоненты матрицы поворота; шк — проекции вектора конечного поворота (вектора Эйлера) на оси глобальной системы координат.

В результате перемещения элемента параллельно самому себе

*

как жесткого целого точка О займет новое положение — О , длины отрезков dsl = ОЦ = |ОА2| станут равными О А = Ж*,

O £

*

= ds2 соответственно. Относительные удлинения сторон мало-

*

_ ds, - ds,

го элемента вдоль направлении S1 и S2 обозначим 81 = ■

ds.

и 82 =

_ ds2 - ds2

ds2

соответственно.

При деформировании векторы локального базиса еъ е2 переходят

* *

в е, . С учетом (3) и (4) повернутые орты будут иметь вид

К = = Р,к Хки4. (6)

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

" ' (7)

ek =akjej-

С cos х sin х^

где а =

sin х cos xj Из (2) и (3) можно получить

dr = ¿d! = P^i^A, dr = e2 ds2 = p2k4ds2, а также учитывая (6) и (7),

(8)

dr = dsi =(! + SiИJpJk4ni„dsi,

dr = e2 ds2 = ( + )o2JpJkXk¿nds2> J = 1> 2

Поскольку вектор перемещения U = r * - r, получим

(9)

OU ^

— = [( + Si )JpJkXkn -pin ] iП,

ou ^

— = [(i + S2 ) OJpJkXkn - p2n ] 4.

(i0)

Уравнения (10) описывают кинематику деформирования срединной поверхности оболочки. Они связывают поворот элемента оболочки, растяжение его сторон и изменение угла между ними с перемещениями элемента.

В основе рассматриваемой математической модели лежит гипотеза Кирхгофа — Лява, подразумевающая малую толщину оболочки и сохранение ортогональности нормалей к срединной поверхности при деформировании. Однако кинематические соотношения (10) легко модифицируются для приближенного учета поперечных сдвиговых деформаций, которые вносят ощутимый вклад в изменение кривизны при деформировании толстых оболочек и пластин.

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

= Д^,* = —«= =

1 11 у у * ^ 1 у у * .У 1

кцзИ кЦзИ

= 'у'уРиЧи1п Р . г ,* = Т1/РзкXк р „ ^ * = (11)

- „ ,* Рз рКрп1паъ1~ „ ,* Рз рКрп1паъ1~ (11)

к013п к013п

Т ■ ^ * ^ *

= „ , * ХфзкРзрХрЛп^Ч = ,

КЦзИ

где — вектор перемещения за счет сдвига; — перерезывающее усилие; к — корректирующий коэффициент сдвига прямоугольного сечения, к^л:2/12; Т1 = Т1 j■ij■, ' = Т2 — векторы внутренних погонных сил; — модуль сдвига ортотропного материала; И* — толщина деформированной оболочки; ¿1п — проекция вектора перемещения за счет сдвига.

Аналогично в другом направлении перемещение за счет сдвига будет иметь следующее выражение:

Т2

= ¿2^2 =

'2 Л

к023И

фък Р3 р Х рп*п^2 = ¿2п*п^2,

(12)

где ¿2 — вектор перемещения в глобальной системе координат;

Т2 л — векторы внутренних погонных сил.

В итоге кинематические соотношения (10) с учетом поперечного сдвига (11) и (12) можно записать так:

дО

дs1

дО д^2

(1 + 81 )(а1 лр лк X кп + ¿1п )-р1, (1 + 82 )(а 2 лР Лк Х кп + ¿2п )-р

гп;

1п.

(1э)

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

к е к*=д1 к1~ л е1, к1 - _ *

э -?**

„1 , п>2

к2 =

деэ

дs0

- ь* — e2, к2 =

де

э

д£*

(14)

Рассмотрим динамическое равновесие элемента оболочки под действием распределенной нагрузки и моментов от внешних факторов Ц, т и инерции самой оболочки ци, ти. По граням элемента прикладываются внутренние погонные усилия и моменты, полученные в результате осреднения напряжений по толщине оболочки (рис. 2). Дифференциалы параметров поверхности представим в виде ds1 = А^а1з ds2 = А2dа2, где а1з а2 — безразмерные параметры на поверхности, А1, А2 — параметры Ламэ.

дМ2 т,ти

+ ——с1а

9, 9и

г1+ ао^"2'^1

дТ-,

Т2 + (1а2 * да2

(к2 = А2(1а2 «¿у, =А\йа

Т, +-г—с1а, 1 дои 1

дМл

Мл +——с1а.л 1 5а! 1

-М,

Ч -Т2

Рис. 2. Схема динамического равновесия элемента оболочки

Методика расчета осесимметричных колебаний оболочек вращения с жидкостью...

Для деформированной оболочки справедливы равенства: ds* = А^а^ = А^а^. Уравнения динамического равновесия сил, записанные для деформированного состояния, имеют вид

1 а а) ! а( а*т* )

А* А* да1 А* А* да*

1 а(*М1) 1 д( А*м*)

+ я = — Яи; (15)

А* А* да1 А* А* да*

1 дт - 1 дт - - --XТ +—-— XТ* + да = -ти,

(16)

А1 да1 А* да * где я = я^ — вектор внешней распределенной нагрузки; яи — сила инерции; М1 = М1 ^, М* = М* ^ — векторы внутренних моментов; даи — инерционный момент; да = дау1у — вектор внешнего распределенного момента.

Здесь необходимо отметить, что инерционная сила и момент вследствие сохранения массы элементарного объема зависят только от ускорения. Поэтому в формулу (15) включена сила инерции

я = рьи Яи = " (1 + в1 )(1 + 8* ), а в (16) введен инерционный момент

рк (О^ + )

да =-----

и 1*(1 + 81 )(1 + 8*) '

где £1 * — угловые ускорения вокруг главных осей (точки обозначают дифференцирование по времени ¿).

Для замыкания системы уравнений необходимо добавить связь между параметрами деформации и напряжениями и граничные условия. Запишем физические соотношения для ортотропного материала в плоском напряженном состоянии [*5, *6]. Принимаем, что оси ортотропии совпадают с направлениями е1, е*. Для плоского случая независимыми являются четыре постоянные: два модуля упругости Е^ Е*, модуль сдвига Оп и коэффициент Пуассона у1* или V*1 (у*1Е1 =у1*Е*) [*5]: растяжение:

Е1к (8 , V 8 \ лг = Е*к

N = 1-1-(81 +V*18*), Ы* = --*-(8* +Vl28l); (17)

1 — V12V21 1 — V12V21

изгиб:

E1h M, = 1

^7271-)+ v21к2)' M2 =

1211 — v12v21)

E2 h3

12 ( — v12v21)

(K2

v12 K

i); (18)

сдвиг:

кручение:

S = G12X;

H =

12

т.

(19)

Здесь N1 = 71^1, N2 = T2e2 — продольные силы; h — толщина оболочки; S1, s2 — параметры продольной деформации; M1 = M1e2, M2 = Me — изгибающие моменты; S = Te* = T^e* — сдвиговое усилие; K2 — параметры изменения кривизны срединной поверхности; х — параметр деформации сдвига; H = M1e1 = M2 e2 — крутящий момент; т - параметр деформации кручения.

Осесимметричная ортотропная оболочка. При осесимметрич-ном деформировании оболочки вращения (рис. 3) координатная сетка на оболочке остается ортогональной, т. е. параметр деформации сдвига х = 0. Введем в рассмотрение цилиндрическую систему координат с ортами

ir = cosфi + sinфi2, /ф =— sinфi1 + cosфi2, = i3. (20)

Радиус-вектор срединной поверхности

r (ф) = r (s) Tr + z( s)Iz, (21)

где r(s), z(s) — параметрические уравнения меридиана оболочки.

"С

-Г»

Ir

Рис. 3. Геометрия и деформация оболочки вращения

3

В качестве параметров поверхности возьмем sx = Гф, = s. Запишем векторы локального базиса:

e1 = - sin ф/j + cos ф12 = 1ф,

e2 = r,s cos Ф h1 + r,s sin Ф h + z,s h = rA + z,s hz, (22)

e3 = z, s c0s ф h + z, s sin ф 4 - rs h = ZJr - rA,

где символ в нижнем индексе после запятой означает дифференцирование по этой переменной.

Таким образом, матрица (3) для осесимметричной оболочки будет иметь вид

Р =

- sin ф cos ф rs cos ф rs sin ф

0

(23)

zs cos ф zs sin ф -rs Запишем кривизны недеформированной срединной поверхности:

k = = — (-zs sinф/ + zs cosфi2)(-sinф/ + cosфi2) = —; (24)

k2 =

двз

ds2

(

-e =

z,sscos ф h + + z,ss sin фh2 - r,ssh

л f

rs cos ф l1 +

+ rs sinф4 + z,sh3

Л

= z r - r z (25)

,ss ,s ,ss ,s \ J

В случае осесимметричного деформирования вектор поворота (5) определяется так:

Ю = -Yei = -ypjk4 = Ysin Щ -Y cos ф4.

(26)

Таким образом, ш1 = у sin ф, ш2 =-ycos ф, ш3 = 0, |ш| = у.

Элементы матрицы поворота для преобразования (4) примут вид

2 2 Xn = 1 - (1 - cos y)cos ф, X22 = 1 - (1 - cos у)sin ф, X33 = cos у,

X12 = -(1 - cos у) sin ф cos ф, X23 = sin у sin ф, X31 = - sin у cos ф, (27)

X13 = sin у cos ф, X21 =-(1 - cos у)sin ф cos ф, X32 =- sin у sin ф.

В силу осевой симметрии деформирования векторы погонных внутренних усилий и векторы моментов имеют вид

T1 = Т1ф1ф, T2 = T2rlr + T2ziz, M1 = M1se2 = M1srJr + M1szJz,

M2 = M2ф1ф.

(28)

Кривизны координатных линий на деформированной поверхности с учетом (6), (23) и (27) можно записать следующим образом:

, * de3 1 de

k1 = —г ei = ds1

i—ei*= i—, ^cosу + sinу— (29)

r (1 + s1 j d9 r (1 + s1 — 7

t * de3 1 de3 1 / \

k2 = 3Гe2 =—-\ — e2 = —-dz,ssr,s " r,ssz,s + Уj. (30)

ds* (1+ S2) ds (1 + 82)

Изменение кривизн при деформировании вычислим по формулам: * 1 í \ Z s

Ak1 = k1 - k1 =-т;—cos r+rssin r)—-; (31)

r (1+ 81) 4 ' r

M2 = k2 - k2 = 77+1-—,ssr,s - r,ssZ,s +r,s ) - ((s - r,ssZ,s ) (32)

I1 + 82 J

Следует отметить, что в выражения для изменения кривизн (31) и (32) входят начальные кривизны. Начальная кривизна меридиана k2 = zssrs - rssz s содержит вторые производные функций r (s) и z (s),

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

При малых продольных деформациях (удлинениях) выражения для приращения кривизн (31) и (32) примут вид

М1 =1 _ z, s (cos Г-1) + rs sin у], (33)

^2 =r,s>

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

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

U = Urir + Uziz = Ur cos ф/1 + Ur sin ф/2 + Uzi3. (34)

Преобразуем выражения для производных перемещений:

dU = 1 dU = 1 -г = г

~ = ~ = Ur гф=81гф,

ds1 r dф r

Методика расчета осесимметричных колебаний оболочек вращения с жидкостью...

откуда следует выражение для окружной деформации (удлинения)

61 = ~иг. (35)

г

Из соотношения (13) с учетом (27) и (12) можно записать

дб Сб

д*2 С1

:[(! + 62)( ] 1Г +[(1 + 82)(-dlQ)-ZsS]I, (36)

где

-, сг = ((* У-^ ^пу), dz =(г* вту + z,s СОБу). (37)

dzT2r - СгТ2z

к023И

Таким образом, получены выражения для изменения кривизн (31), (32) и для удлинения в окружном направлении (35), а также дифференциальные уравнения для перемещений (36).

Теперь получим уравнения равновесия и физические соотношения с учетом изменения объема при больших перемещениях. В приведенном выше материале использовалась естественная параметризация функций — *1, *2, *, за исключением уравнений равновесия (15) и (16), где введены безразмерные параметры на поверхности — а17 а2 и параметры Ламэ А1, А2. Однако на практике стоит предпочесть произвольную безразмерную параметризацию. Для случая осе-симметричной оболочки введем безразмерный параметр длины меридиана £, т. е. * = * (£). Определим параметры Ламэ с помощью следующих соотношений:

с*

С*1 = А1Са1 = гСф, С*2 = А2Са2 = — С£,, (38)

С £

т. е. параметры Ламэ А1 = г, А2 = *£ и С а1 = С ф, С а 2 = С £.

Вычислим производные, входящие в уравнения равновесия (15), (16) по безразмерным параметрам поверхности:

д[А2 (1 + 82 )Т ] = С [А2(1 + ^Уф] А (1 + 8 ) = А (1 + 8 )т .

I = = А2(1 + 62) = -А2(1 + 82)т1ф'г,

да1 С ф С ф

д[ А1 (1 + 81) Т2 ] = С [ г (1 + 81Т ] = С [г (1 + 81)Т2 г ] 1 + С [г (1 + 81 )Т2 z ] г да 2 С £ С £ г С £ ^

д[ А2 (1 + 82 )М1 ] . ( )С (М1*гЛ +M1sZ,s4 ) , ( .

= А2 (1 + 82) —--'- = А2 (1 + 82)М1*гЛ,

да1 С ф

5[ А (1 + £1 )М2 ]_ й [г (1 + 8!) М2фГф]_ й [г (1 + 81) Ы^] ^

5а 2

й %

й%

1ф ■

Таким образом, уравнения динамического равновесия сил в проекциях на оси цилиндрической системы координат 1Г, ¡г имеют вид

й [г (1 + 81 )Г2г й %

й [ Г (1 + 81) г й%

- А2(1 + 82)Т1ф + гА2 (1 + 81) (1 + 82 ) _ АрИО,, + А (1 + 81) (1 + 82 )дг _ тА2рИИг,

(39)

где Ц _ цггг + — внешняя нагрузка.

Уравнение равновесия в проекции на ось /ф удовлетворяется тождественно.

Уравнение равновесия моментов в проекции на ось /ф можно представить так:

й [Г (1 + 81 )М2ф

й %

+ А2 (1 + 82 )М1,Г, +

+ А (1 + 81 )(1 + 82)

Т2г ( ^ Ч + с°8 у)- Т2 г ( С°8 У- ^ У)

1 рИ3..

_-гА2-—у.

2 12

(40)

Физические соотношения (17)-(19) для ортотропной осесиммет-ричной оболочки примут вид

_ Т1ф _—ЕИ-(81 + V2182);

1 -^21

(41)

N2 _ Т2г

Сг5 собу-V- ^п уу

+ Т

2 г

С г5 у + ^ + г 5 соб у

(( +Vl28l); (42)

1 -^21

Ы1 _ Мь _

Е1И

*3

12 (1 -V12V2l)

*3

М2 _-М2ф _

Е2И

12 (1 -V12V21)

(( + V 21М2);

(Ак2 + Vl2 А),

(43)

(44)

где 81, 82 — меры продольных деформаций; И — толщина деформированной оболочки.

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

81 = 1п (1 + 81), 82 = 1п (1 + 82 ). (45)

Необходимо отметить, что при малых деформациях или при геометрически линейной постановке задачи удлинения равны и меры деформации: 81 =81, 82 = 82. Это следует в том числе из формул (45).

Определим изменение толщины при деформировании через изменение объема. До деформации элементарный объем

V = р^С*1С*2, после деформирования его можно записать так:

V * = р/г*С*1*С*2 = ри/И (1+81 )(1+82) С*1С*2,

где новую толщину представим в виде произведения старой толщины

и некоторой функции, описывающей изменение толщины, подлежа*

щей определению: И = И/и.

Имеет место следующее отношение:

V = М„ (1 + е,)(1+Е2 ) = И ( +82 ),

но, с другой стороны, известно, что

*

V = (1 + 81 )(1 + 82 )(1 + 83 )

следовательно, / = 1 + 83.

Запишем закон Гука для ортотропного материала:

81 = т — -у31 тт,

Т1 Т2 Т3

82 =-У12 Т + Т Уз2 Т,

М 2 ^3

83 — VI3 V23 I .

Г Т7 Т7

Поскольку рассматриваются тонкие оболочки, положим = 0, тогда

8 = У^^+У^8 , V21V13 + V23 8 ь3~ 1 Ь1 + 1 2,

V12V 21 -1 V12V21 -1

и с учетом (45) в итоге для функции изменения тол-

щины получим выражение

/и _ ехР

^23 +^3" , V21V13 +V23

V V12V 21 - 1 1 ^ 21 -1 2,

(46)

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

5 _5, И_И, г_Г, г_г, и_и, Т= Т,

£ £ £ £ £ ЕИ £

- 12(1 -V2 )

М _ ^ 2 ! М,

ЕИ £

(47)

где £ — характерный геометрический размер оболочки (например, длина меридиана, радиус), V _ ^^21, Е _У]Е11Е2.

Кроме того, перейдем к параметризации функций с помощью %. Физические соотношения (42)-(44) с учетом (47) можно переписать следующим образом:

Т1ф _ Л Е(81 +V2182 );

(48)

соб у-г,% вШ у) +

+ Т2г ((% 81П У + С08 у)

Е

Е

(49)

М15 _ЛЕИ£11 (( +v21Ак2);

\Е.

(50)

М2ф_ ЕИ (Ак2 +^2Ак1).

(51)

Из (49) определим меру деформации в меридиональном направлении:

Е 1

8 _

Е2 5,%Л

Т2Г (г% соб у - г, % ЯП у) + Тг ((% Б1П у + г, % соб у)] - Vl28^^ (52)

Здесь необходимо отметить, что в выражение (52) входит (46), в результате чего образуется нелинейное уравнение относительно меры деформации 82, которое в программной реализации решается методом простых итераций.

Упростим и перепишем уравнения равновесия (39) с учетом (47) и (38):

1

dT2

2r

d £

( (1 + б!)

- + J£ С1+82) T _

T2r + - ¡Л . _ \ -1ф

(1 + 81)

_ JA(1 + 82 q + чЙ^ (1 _V12V21) TUГ

hJEE r 4ёЁ2 (1+81)

dT2

2 г

d £

íi+_8u_

( (1+81)

-2 г _

- (1 + )(1 _V12V21) + J,£P^2 (1 _ V12V21) TT _(1+82) - ^ ^ qz + ■ ^-^u,;

г jee (1+81)

dM.

2ф

d £

_ h(1 + 82)

i1,£

(1 + 81)

M2ф_

(1 + 82 ) M

2ф —

(1 + 81)

Mu _

T2r ((

£ sin y + z.£ COS y)_ _T2г (( COS У_ г,£ Sin У)

J.£Ph l2 (1 _ V12V21)

■JEE (1 + 81 )

y-

Из физического соотношения (51) выразим у5, при больших

продольных деформациях из выражений для изменения кривизн (31), (32) получим

У - =■

где А/ =■

"(1 + 82 )/Ет ^f _ V12 (1 + 82 )Ak1 + 82 j _ (JJZJ ). (53) VE2 lhJh

, [cos у _ (1 + 81)] + rs sin У

(1 + 81 )

При малых продольных деформациях из выражений для кривизн (33) получим

E M 2ф А/

У. J = _J _ V12 ^1.

(54)

\е2 ihf3

где А/1 = — [z-s (cos у _ 1) + rs sin у].

Таким образом. вычислив у s по формуле (53) или (54). изменение кривизны А/* можно рассчитать по (32) или (33) соответственно. и тогда M- возможно определить из физического соотношения (50).

С учетом (47) и переходом к параметризации по £ окончательно в итоге имеем разрешающую систему шести геометрически нелинейных

дифференциальных уравнений первого порядка относительно шести функций — иг, иг, у, Т2г, Т2г, М2ф, описывающих динамическое

напряженно-деформированное состояние осесимметричной ортотроп-ной оболочки при больших продольных деформациях, перемещениях и поворотах:

кинематические соотношения:

1) _(1 + 82 )(( + 2) йит _(1 + 82 )(( - 3,0)-7%; (55)

уравнение упругости при больших продольных деформациях:

3)4% _ ЕМ2ф - Vl2 (1 + 82)Ак1 + 52((%-г, %), (56)

J С

где Akj =

W

S cY

(j + )

z? [cos y-(j + 8j )]-

Y? Sin у

—2

Ak2 _

s,i s2(^ ^

—3

(j + S2 )

уравнение упругости при малых продольных деформациях:

3а)

d у

-M2, - v,2s cAk,

dc hf3Vф" 124 j'

__j __у

где Akj z> (cos у-i) Y? sin у], Ak2 ;

S?Y L 'S 's J S ?

уравнения равновесия:

dT2 - - s>(1 + s2 ) ^

"df _ -^y BT, - C ' h Ъ + DUr;

dT2 - S-£(1 + s2) ^ 5) _ -AT2z - C ^ ^ + DUZ;

d С

h

(57)

(58)

(59)

6)

dM-

2,

d С

AM2, - Y BMjs -12 (1 + S2 ) (Td - T2zdr ) - Dh у; (60)

дополнительные выражения и обозначения: dr _( cos у-Z"c sin у), dz _(Y? sin у + Z^ cos у), 9

A _

(j + S1)

fYL и - .1 dU Л

_2 r —

Y2 -

B _

dzT2r drT2z kG23\?fhC

C _-

(j + S2 ) C _ 1 -v12v21

r ) (i+Sj)

С

j

Методика расчета осесимметричных колебаний оболочек вращения с жидкостью. " 1 и 8 = ехР(8 ) 1 8 = 1Т1 Т2гаг + 12zaz V

12ь1,

stр£ 1 _

Э = С / „ 8 "

(1+81)

Л (82 ) = ехР

1 — \ 1 - /Т1 + ТЛ,

1 == иг, 82 = ехр (82)-1, 82 = -^8,

V12V23 + ^3" , V21V13 +V23

V12V 21 -1 1 V12V 21 - 1 2

Т1Ф = Л Т-/и ( + V2182 ), Мь = IТг-И/И (( + V2,Ак2 ).

V Т2 V т2

Линеаризуя систему уравнений (55)-(60) относительно деформаций, поворотов, перемещений, усилий и моментов, получим геометрически линейную систему уравнений движения осесимметричной оболочки:

1) = Г'^82 - 2>У + 2) = ^82 + - ^

Т2 \%/и (82 )

81 = 1п (1 + 81),

д У = Ъ Е1 М

3) ^ = -Т/Г* ЕТМ2ф-^2^У а ? И/, у т2 "

4) ^ = (Vl2 - 1)Т2г + Vl2 Т2, + (1 - Vl2V2l))

а ?

^ ет

7 \ Е2

81 -

- 1 ^^ Чг+

(61)

И4Е\Е2 7 , л/Е1Е2

5)

ат2

а ? 7

7,? ^ _ 1 -V1^ V

12У21 Ч2 + ^р*2,12*21 и,;

1 -V, 2V-,, -ч-

И ^/Е,Е2 . , ^Е1Е

—2

ам2ф . . г? _ . хг? и ¡е,

6) ^Г = (Vl2 - 1)м2ф - (1 - Vl2V2l ))И. Е У -

а ? г ^ г2 . ? V Е2

где 8, =

0 =

- И (Г,-*

81 = I1 и г, 82 = ^ |-Ет (г + Т 2 г^?)-^ 8l,

7 V Е2

Г,?Т2г - г,?Т2, \1Е1Е2

,? V 2 у,

(Т - V12V21)

Для задач статики деформирования тонких и толстых оболочек под давлением уравнения (55)-(60) (при иг =у = 0) были успешно протестированы вне рамок данной статьи. Расчет геометрически

линейного и нелинейного деформированного состояния сопоставлялся с расчетами в среде АКБУБ. Рассчитанные контрольные величины перемещений совпали в пределах разброса вычислений по разным конечным элементам АКБУБ (8ЫБЬЬ181, 43, 93). Вычисление изменений кривизны по формуле (57) дали лучшее приближение к решениям, получаемым в АКБУБ.

Сопряженная задача гидроупругости. Рассмотрим расчет малых колебаний осесимметричного бака с идеальной несжимаемой жидкостью (рис. 4). Для такой модели жидкости ее малые перемещения определяются потенциалом перемещений ф [14], а перемещения стенки бака — системой дифференциальных уравнений движения осесимметричной оболочки вращения (55)-(60), которые в сокращенном виде можно записать так:

Р и, X, X),

(62)

где X /)_{иг иг у Т2г Т2г М2ф) — вектор разрешающих функций.

<5ср дп

У

= О

Дф = 0

Ф = о

и-п = -т— у дп

Рис. 4. Расчетная схема бака с жидкостью

В случае малых колебаний решение разыскивается в виде X(%, I) _ АХ(%) ехр (Ш).

(63)

Теперь, подставляя решение (63) в (62), линеаризуя и вычитая исходное равновесное состояние, получим уравнение для амплитуд малых колебаний

й АХ

ар - 2 др_ аХ ю аХ

АХ.

(64)

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

тате решения задачи о взаимодействии. Потенциал перемещений на смоченной стенке определяет давление жидкости на эту стенку:

Р _-Ржф® , (65)

где рж — плотность жидкости; ш — частота собственных колебаний.

В рамках допущений для потенциала малых перемещений жидкости ф получим краевую задачу [14]: Аф _ 0 — в объеме жидкости; (66)

ф_ 0 — на свободной поверхности жидкости; (67)

дф / дп _ 0 — на оси симметрии и/или неподвижной части бака; (68)

ЦУп _дф/ дп, Т2п _ -ржф®2 — на смоченной части оболочки, (69)

где п — внешняя нормаль.

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

Применение КБЕ-ЕВ. Для решения краевой задачи (66)-(69) применим метод конечных разностей, а для генерирования весовых коэффициентов разностных схем на произвольном трафарете будем использовать радиальные базисные функции.

Радиальная функция — это любая вещественная функция, значение которой зависит только от расстояния до начала координат

ф(г)_ф(|Щ) или от расстояния между некоторой другой точкой,

называемой центром с : ф(г,с)_ф(||г -с||). В качестве нормы обычно выступает евклидово расстояние. Существует много видов РБФ, наиболее часто используемыми являются: мультиквадратичная [20]:

ф(г)_,! 1 + (вг )2, (70)

и полигармоническая [18, 19, 24]:

ф( г )_ г 2т 1п (г ), (71)

где в — параметр формы; т — степень четной полигармонической РБФ.

Использование мультиквадратичной РБФ сопряжено с необходимостью выбора параметра формы £, зависящего от расположения узлов интерполяции. Выбор оптимального значения, минимизирующего ошибку аппроксимации, — дополнительная задача, для решения которой предложены разные подходы [20, 27]. Полигармонические РБФ (71) не имеют параметра формы, однако есть некоторые умолчания

относительно практического применения, не освещенные или недостаточно освещенные в литературе. В данной работе используются полигармонические РБФ.

Тестовые расчеты задачи о колебаниях бака с жидкостью показали, что для использования (71) в качестве РБФ интерполирующего сплайна необходимо приводить нормирование размера трафарета, для того чтобы избежать плохой обусловленности разрешающей матрицы, приводящей из-за ошибок округления чисел в ЭВМ к критическим ошибкам при определении весовых коэффициентов конечных разностей. Авторы опытным путем установили, что «рабочей» областью полигармонической РБФ вида (71) является интервал для радиусов больше единицы, т. е. необходимо нормировать координаты точек в трафарете таким образом, чтобы все расстояния между любыми точками были больше 1 (на практике можно установить число, близкое к 1, например 1,01).

Классический метод конечных разностей предполагает разбиение геометрической области узлами, лежащими на координатных линиях (прямоугольная сетка). Для такого расположения узлов записываются простые и достаточно точные разностные схемы аппроксимации различных производных по координатам через узловые значения разыскиваемых функций. Использование РБФ позволяет вычислять весовые коэффициенты таких разложений на произвольной (нерегулярной) сетке, имеющей сгущение там, где это требуется (например, в области взаимодействия жидкости и оболочки), а в других местах иметь разреженную структуру для экономии ресурсов ЭВМ. Можно добиться аналогичных нерегулярностей с ортогональной сеткой, используя деревья, но при их реализации возникают дополнительные трудности.

Быстрый алгоритм генерирования нерегулярной сетки узлов в прямоугольной области предложен в [28]. Этот алгоритм был слегка модифицирован для генерирования сетки из предварительно дискрети-зированной выпуклой граничной области с произвольной дискретизацией. Пример результата работы этого алгоритма приведен на рис. 5, где наряду с граничными и внутренними точками области показаны законтурные (вспомогательные) узлы (англ. ghost nodes). Они вводятся для уравнивания числа уравнений и неизвестных. Таким образом, уравнение Лф = 0 должно выполняться во всех узлах границы и внутренней области, а законтурные узлы участвуют в формировании трафаретов конечных разностей. Число законтурных узлов строго равно числу граничных, для которых заданы граничные условия (на зеркале жидкости, на оси симметрии, на смоченной части). Трафарет, составленный из 22 точек, показан на рис. 5. Подобные трафареты строятся

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

1Д 1,0 0,9 0,8 0,7 0,6 0,5 0,4 0,3 0,2 ОД 0

-0,1 0 0,1 0,2 0,3 0,4 0,5 0,6 Рис. 5. Разбиение расчетной области узлами

Рассмотрим алгоритм вычисления весовых коэффициентов для разностных схем. Интерполяционный сплайн на основе РБФ записывается в стандартном виде с дополнительным полиномом первой степени [29], гарантирующим уникальность решения [18], а также с условиями ортогональности весовых коэффициентов:

п

8 (х у ) = Ф(х - хк, у - Ук )+У1+У2х+Узу; (72) к=1

п п п

к кхк = кУк = ° (73)

к=1 к=1 к=1

где (хк, Ук) — узлы интерполяции; п — число узлов интерполяции; Хк — весовые коэффициенты РБФ; У123 — коэффициенты дополнительного многочлена.

Сплайн (72) может быть расширен до пространственного случая, а также до произвольного количества размерностей.

Пусть интерполируемая функция имеет узловые значения /,

тогда из (72), (73) СЛАУ для определения весовых коэффициентов будет иметь вид

ф(( - г) ф(( - Г) . . ф( Г1 - гп) 1 Х1 У1 V /1

ф((- Г) ф (( - Г) . . ф(Г2 - П,) 1 х2 У2 Ч2 /2

ф(Гп - Г) ф(Гп - ?2) . . ф(( - Гп ) 1 Хп Уп Ч п » = < /п

1 1 .. 1 0 0 0 У1 0

Х1 х2 . хп 0 0 0 У2 0

У1 У2 . Уп 0 0 0 .У3. 0

(74)

[А]

СЛАУ (74) имеет уникальное решение, если узлы интерполяции не повторяются (это условие обеспечивается генератором распределения узлов), а ее матрица [А] не меняет свой вид при аппроксимации произвольного дифференциального оператора Ь [29], что легко показать. Пусть необходимо найти весовые коэффициенты Жу разложения по узловым значениям функции /у дифференциального оператора Ь в центральной точке хс, ус трафарета из п узлов:

I У

у =1

п п п

п+1I4 у + Жп+2I4 уХу + Жп+3I4 уУу

У=1 У =1 У =1

ЬБ (х, у)|х=Хс . (75)

У=Ус

Здесь в силу условия ортогональности (73) слагаемые в квадратных скобках равны нулю, т. е. значимыми в разложении будут только п весов ж у. С учетом (72) и (74) имеем

(

+ Ж

I жу IЧкф (ху - хк, У у -Ук) + У1 + У2ху + УзУ у у=1 Ц=1

п п п

IЧУ + Жп+21ЧУХУ + Жп+31Ч уУу = ЬБ (x, У |

У=1

У=1

У=1

х=хс ■< У=Ус

п п

1ЧкIжуф(ху -хк, Уу -Ук) + У11"у +У21У +У31жуУу

к=1 У=1

У=1 У=1

п

У=1

+ ж

п+1

I Чк + Жп+21 ЧкХк + Жп+31 ЧкУк = Ь Б ( У)|х=Хс ;

к=1

к=1

к=1

У=Ус

п ( п

Ё ■*к Ё ^Ф(( " Хк, У/ " Ук ) + ^п+1 + ^п+2хк + ^п+3Ук к=1

+ У1Ё +У2 Ё + Уз Ё = 7=1 /=1 /=1

( п Л

Ё* кЬф(х " Хк, У - Ук ) + У1Ь 1 + У 2Ьх + У 3 ЬУ V к=1

х= хс У=Ус

В итоге получим СЛАУ для определения весовых коэффициентов разложения дифференциального оператора по узловым значениям функции //, в которой матрица совпадает с (74):

[А]

м>1 Ьф( х - хь У - У1Г

ь -©- - ^ У - У2 )

wn =■ = ■= ь -©- - хп, У - Уп )

^п+1 ь 1

^п+2 Ьх

. ^п+3. ЬУ

(76)

х=хс

У=Ус

Расчетная область, занятая жидкостью, разбивалась неравномерной сеткой, имеющей сгущение на границах и разреженность внутри (см. рис. 5). Это дает существенную экономию размерности задачи и незначительно сказывается на точности решения. Для интерполяции характерного размера сетки узлов от границ внутрь области применялся многочлен (72), но в качестве РБФ использовалась функция

ф(х, У) = х2 + У2, более устойчивая с вычислительной точки зрения,

чем полигармоническая РБФ (71), однако она не подходит для аппроксимации дифференциальных операторов рассмотренной задачи.

Алгоритм решения уравнений ортотропной осесимметричной оболочки совместно с жидкостью добавлен в функциональные возможности программы БАЯБУБ [30].

Тестирование методики расчета собственных колебаний ортотропной оболочки без жидкости. Рассмотрим тестовую задачу расчета собственных колебаний ортотропной оболочки вращения, состоящей из защемленной снизу цилиндрической и свободной конической частей (см. рис. 4). Радиус цилиндра составляет 0,5 м, высота цилиндра — 0,8 м, длина образующей конуса — 0,3 м, угол конуса относительно вертикали — ж /6, толщина стенки — 0,02 м. Для

п

подтверждения правильности разработанных уравнений рассмотрим процесс колебания оболочки из трех материалов: 1) изотропного как частный случай ортотропного, т. е. Е1 = Е2; 2) ортотропного с соотношением модулей упругости Е1 = 2Е2; 3) ортотропного с соотношением Е = 5Е2. Примем модуль упругости материала в окружном направлении Е1 = 2 -10 Па, коэффициент Пуассона у12 = 0,3, плотность материала р = 7850 кг/м3. Частоты собственных колебаний оболочки без жидкости, рассчитанные методом конечных элементов в программе В.Е. Левина [14], в пакете АКБУБ и по предложенной методике в БЛЕФУЯ, приведены в табл. 1-3.

Таблица 1

Частоты собственных колебаний, Гц, оболочки из изотропного материала при Е1 = Е2, определяемые с помощью трех способов

Тон МКЭ [14] ЛШУЯ БЛЯЯУЯ

1 10,488 10,464 10,463

2 15,345 15,326 15,324

3 16,538 16,514 16,510

4 17,613 17,578 17,575

5 18,337 18,272 18,258

Таблица 2

Частоты собственных колебаний, Гц, оболочки из ортотропного материала с соотношением Е1 = 2Е2, определяемые тремя способами

Тон МКЭ [14] ЛШУЯ БЛЯЯУЯ

1 7,94367 7,951 7,9475

2 14,3333 14,314 14,340

3 16,1666 16,120 16,206

4 16,5995 16,518 16,585

5 16,8880 16,982 17,087

Таблица 3

Частоты собственных колебаний, Гц, оболочки из ортотропного материала с соотношением Е1 = 5Е2, определяемые тремя способами

Тон МКЭ [14] ЛШУЯ БЛКБУБ

1 5,197 5,200 5,1979

2 12,091 12,062 12,072

3 15,955 16,087 15,967

4 16,205 16,252 16,192

5 16,356 16,676 16,417

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

Тестирование методики расчета гидроупругих колебаний. Пример 1. Рассмотрим приведенную в статье [13] тестовую задачу расчета колебаний заполненного водой полусферического открытого бака, шарнирно закрепленного по радиусу R = 5,08 м. Толщина бака

0,0254 м, модуль упругости E = 7 -1010 Па, коэффициент Пуассона

v = 0,3, плотность материала р = 2770 кг / м3.

Три низшие частоты осесимметричных гидроупругих колебаний, рассчитанные различными методами в разных программах, приведены в табл. 4. Следует отметить, что частоты, взятые из статьи [10], не сопровождаются информацией о сходимости в зависимости от дискретизации модели. Частоты, полученные расчетом в программе В.Е. Левина [14] по МКЭ—МГЭ, практически не меняются при увеличении числа конечных и граничных элементов. Частоты, определенные в ANSYS Workbench 14.5, возникали при ~50 разбиениях вдоль меридиана.

Таблица 4

Низшие частоты колебаний полусферической оболочки с водой, Гц, определенные шестью разными способами

Тон Частота колебаний, Гц

МКЭ-МГЭ [13] МКЭ [10] ANSYS [10] (Shell63) МКЭ-МГЭ [14] ANSYS (Shell181, Fluid80) DARSYS

1 23,59 22,00 22,07 22,04 22,333 20,296

2 35,70 33,38 33,41 33,32 33,714 35,172

3 43,92 42,02 41,30 41,10 39,091 40,637

Из данных табл. 4 следует, что расхождение частоты первого тона, полученной разными методами, достигает 16 %. Для анализа сходимости были проведены расчеты с различной дискретизацией модели. Графики сходимости рассчитанных частот в ANSYS и по разработанной методике представлены на рис. 6, где хорошо видно, что кривые пересекаются при возрастании числа разбиений. Сходимости частот в ANSYS Workbench 14.5 (Shell181 + Fluid80) получить не удалось, а частоты, рассчитанные по предлагаемой методике, имеют выраженную асимптотику при увеличении числа разбиений меридиана. На основании этого можно предположить, что расчет по предлагаемой методике более точен.

100 200 Число разбиений вдоль меридиана

300

20,4

Рис. 6. Сходимость частот гидроупругих колебаний полусферической оболочки с водой

100 150 200

Число разбиений меридиана

300

Рис. 7. Низшая частота колебаний в зависимости от дискретизации меридиана и числа точек в трафарете п, степень РБФ т = 1

Зависимости низшей частоты от дискретизации меридиана при различных размерах трафарета конечных разностей п = 15, 22, 50 для степеней четных полигармонических РБФ т = 1 и т = 2 соответственно показаны на рис. 7 и рис. 8. На рисунках видно, что область разброса вычисляемых частот при изменении дискретизации меньше для т = 2. Примеры распределения узлов, в которых средний размер дискретизации границы увеличивается в 2 раза при приближении

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

20,4

20,2 I 20,0

| 19,8 19,6 19 4

'0 50 100 150 200 250 300

Число разбиений меридиана

Рис. 8. Низшая частота колебаний в зависимости от дискретизации меридиана и числа точек в трафарете n, степень РБФ m = 2

1) ■ (I) 19.B4S Гц t) . (1) 10,1» Гц 1) - (1) :0,307 Гц

Рис. 9. Первая форма колебаний с различной дискретизацией

Пример 2. Частоты собственных колебаний составной оболочки из изотропного материала (из первой тестовой задачи), заполненной водой до уровня 0,5 м, приведены в табл. 5. Были выполнены расчеты частот собственных осесимметричных колебаний в программе [14] по МКЭ—МГЭ (~60 элементов), в ANSYS Workbench 14.5 моделировалась четверть модели с условиями симметрии (—119 000 элементов) и по предложенной методике в программе DARSYS при 300 разбиениях меридиана (см. табл. 5).

Таблица 5

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

Тон ANSYS (Shell181, Fluid80) МКЭ-МГЭ [14] DARSYS

1 8,464 9,233 9,394

2 10,969 11,256 11,307

3 13,801 14,137 14,130

4 16,377 16,474 16,448

5 16,763 16,937 16,869

6 18,487 18,756 18,577

По данным табл. 5 видно, что наибольшее различие, которое приходится на низшую частоту, достигает 11 %, а с ростом номера тона разница падает до 0,5 %. Следует отметить, что уточнить результаты, полученные в ANSYS, путем увеличения числа КЭ не удалось.

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

Пример 3. Рассмотрим задачу из работы [11], в которой происходят колебания цилиндрического бака радиусом 7,25 м и высотой 12,6 м, наполненного водой до уровня 11,49 м. Толщина стенки бака 7 мм,

модуль упругости E = 2-1011 Па, коэффициент Пуассона 0,3, плотность материала р = 7850 кг/м3. Частоты осесимметричных колебаний, определенные в ANSYS Workbench 14.5, в программе из [14] и по разработанной методике, приведены в табл. 6. Низшая частота, полученная в ANSYS, отличается на 32 %, однако с ростом номера тона это различие уменьшается до 8 % на тоне 5. При этом частоты, полученные по МКЭ—МГЭ [14] и по предлагаемой методике, хорошо согласуются между собой.

Таблица 6

Частота колебаний цилиндрического бака с водой, Гц, определенные тремя методами

Тон ANSYS (Shell181, Fluid80) МКЭ-МГЭ [14] DARSYS

1 6,410 8,414 8,441

2 14,696 15,751 16,484

3 19,503 20,700 21,026

4 22,966 24,567 24,547

5 25,680 27,834 27,485

Пример 4. Рассматривается тороидальный бак, наполовину заполненный водой, образующая окружность имеет радиус 0,5 м, внутренний радиус 0,5 м, внешний — 1,5 м. Оболочка толщиной 0,02 м из

изотропного материала с модулем упругости Е = 2 -10 Па, коэффициент Пуассона V = 0,3, плотность

р = 7850 кг / м3. Тор шарнирно закреплен по внешнему радиусу. Частоты и формы колебаний, рассчитанные в программе [14] (140 КЭ, 140 ГЭ вдоль меридиана), и по предлагаемой методике в программе БЛЯБУБ [30] (140 разбиений меридиана, т = 2, п = 22), показаны на рис. 10. Частоты осесимметричных колебаний тороидального бака первых трех низких тонов для сравнения приведены в табл. 7. По ее данным видно, что частоты близки.

Рис. 10. Формы трех низших тонов осесимметричных колебаний тора, полученные при расчете в программах МКЭ—МГЭ [14] (а) и БЛЯЗУБ (б)

Таблица 7

Частоты колебаний тороидального бака с водой, Гц, вычисленные двумя способами, и разница между ними

Тон МКЭ-МГЭ [14] БЛЯЗУБ Разница, %

1 0,5616 0,5480 2,42

2 3,2617 3,2118 1,53

3 3,9653 3,8970 1,72

4 4,4991 4,4930 0,14

5 5,2194 5,1741 0,87

6 6,0651 5,9365 2,12

Зависимости трех низших частот от числа разбиений меридиана, рассчитанные в БЛЯБУБ, приведены на рис. 11, на котором видно, что с ростом числа разбиений частоты стремятся к горизонтальным прямым без заметных осцилляций.

О 50 100 150 200 250 300

Число разбиений меридиана

Рис. 11. Зависимости низших частот тороидального бака от дискретизации

Заключение. Из общей системы уравнений движения оболочки, записанной для разрешающих функций в глобальной системе координат, представлен вывод системы уравнений ортотропной осесиммет-ричной оболочки, описывающих геометрически нелинейное динамическое деформирование с учетом утонения/утолщения, поперечного сдвига и больших продольных деформаций. На основе этих уравнений путем линеаризации и вычитания исходного равновесного состояния формулируется линейная краевая задача на собственные значения, в результате решения которой вычисляются частоты и формы малых колебаний. Для описания малых движений идеальной несжимаемой жидкости используется потенциал перемещений, краевая задача для которого аппроксимируется с помощью конечных разностей на произвольной нерегулярной сетке узлов. Весовые коэффициенты конечных разностей вычисляются с помощью сплайн интерполяции на основе полигармонических радиальных базисных функций. Составляется матричная обобщенная задача собственных значений, построенная на дискретизации МКР уравнений осесимметричной оболочки и краевой задачи для жидкости, которая решается численно. Разработанная методика расчета внедрена в собственный программный комплекс для расчета стержневых систем DARSYS и протестирована на нескольких задачах. Показана сходимость рассчитываемых частот при увеличении числа разбиений меридиана оболочки для разных размеров трафарета МКР и двух степеней четной полигармонической РБФ. Предлагаемый подход к расчету частот гидроупругих колебаний осесимметричных ортотропных оболочек с помощью РБФ-МКР (RBF-FD) будет применен при расчете общего случая пространственных оболочек и пластин.

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

ЛИТЕРАТУРА

[1] Колесников К.С., Рыбак С.А., Самойлов Е.А. Динамика топливных систем ЖРД. Москва, Машиностроение, 1975, 172 с.

[2] Wei Liu, Chang Xiao, Hao Zhou, Chenyan Wang. Experimental investigation of liquid-tank interaction effects on full containment LNG storage tanks through shaking table tests. Thin-Walled Structures (2023). https://doi.org/10.1016/j.tws.2023.111527

[3] Грибков В.А., Адаменко Р.А. Двумерная модель жидкости для расчета собственных частот колебаний осесимметричных гидрооболочечных систем. Инженерный журнал: наука и инновации, 2017, вып. 3 (63), с. 5. http://dx.doi.org/10.18698/2308-6033-2017-3-1593

[4] Шупиков А.Н., Мисюра С.Ю., Ярещенко В.Г. Численное и экспериментальное исследование гидроупругих колебаний оболочек. ВосточноЕвропейский журнал передовых технологий, 2014, т. 6, № 7 (72), с. 8-12. DOI: 10.15587/1729-4061.2014.28861

[5] Jhung M.J., Jo J.C., Jeong K.H. Modal analysis of conical shell filled with fluid. J. Mech. Sci. Technol, 2006, vol. 20, no. 11, pp. 1848-1862.

[6] Гончаров Д.А., Пожалостин А.А., Кокушкин В.В. Моделирование осесим-метричных колебаний упругого бака с жидкостью с учетом сил поверхностного натяжения посредством механического аналога. Наука и образование: научное издание МГТУ им. Н.Э. Баумана, 2015, № 6, с. 372-383. DOI: 10.7463/0615.0779724

[7] Phan, Hoang Nam, Fabrizio Paolacci. Fluid-structure interaction problems: An application to anchored and unanchored steel storage tanks subjected to seismic loadings. arXiv: Numerical Analysis (2018): n. pag.

[8] Шклярчук Ф.Н., Рей Чжунбум. Расчет неосесимметричных колебаний оболочек вращения с жидкостью методом конечных элементов. Вестник МАИ, 2013, т. 20, № 2, с. 49-58.

[9] Шклярчук Ф.Н. Расчет колебаний оболочек вращения с жидкостью методом конечных элементов. Проблемы машиностроения и надежности машин, 2015, № 1, с. 17-29.

[10] Мокеев В.В. Исследование динамики конструкций с жидкостью и газом с помощью метода конечных элементов. Изв. РАН. Механика твердого тела, 1998, № 6, с. 166-174.

[11] Yi Lu, Ji Jing. Modal analysis on anchored tank considering shell and fluid coupling. Advanced materials research, 2012, vol. 549, pp. 903-907. Crossref, DOI: 10.4028/www.scientific.net/amr.549.903

[12] Dubois Т.Т., de Rouvray A.L. An improved fluid superelement for the coupled solid-fluid-surface wave dynamic interaction problem. Earthquake Eng. Struct. Dynam., 1978, vol. 6, no. 3, pp. 235-245.

[13] Гнитько В.И., Дегтярев К.Г., Кононенко Е.С., Тонконоженко А.М. Сравнение методов конечных и граничных элементов в задачах о колебаниях составной оболочки вращения с жидкостью. Вiсник Харювського нацюналь-ного унiверситету iменi В.Н. Каразiна, 2019, вип. 42, c. 38-45.

[14] Левин В.Е. Метод конечных и граничных элементов в динамике конструкций летательных аппаратов: 05.07.03 «Прочность и тепловые режимы

летательных аппаратов»: дис. ... д-ра техн. наук. Новосибирск, Новосибирский государственный технический университет, 2001, 341 с.

[15] Бочкарев С.А., Лекомцев С.В., Матвеенко В.П. Собственные колебания усеченных конических оболочек, содержащих жидкость. Прикладная математика и механика, 2022, т. 86, № 4, с. 505-526.

DOI: 10.31857/S0032823522040038

[16] Nurul Izyan M.D., Viswanathan K.K., Aziz Z.A. et al. Free vibration of layered truncated conical shells filled with quiescent fluid using spline method. Compos. Struct., 2017, vol. 163, pp. 385-398.

[17] Mohammadi N., Aghdam M.M., Asadi H. Instability analysis of conical shells filled with quiescent fluid using generalized differential quadrature method. In: The 26th Annual Int. Conf. of Iranian Society of Mechanical Engineers-ISME2018, 24-26 April, 2018. School of Mechanical Engineering, Semnan Univ., Semnan, Iran, ISME 2018-1216.

[18] Flyer N., Fornberg B., Bayona V., Barnett G.A. On the role of polynomials in RBF-FD approximations: I. Interpolation and accuracy. Journal of Computational Physics, 2016, vol. 321, pp. 21-38. DOI: 10.1016/j.jcp.2016.05.026

[19] Shankar V., Wright G.B., Kirby R.M., Fogelson A.L. A Radial Basis Function (RBF)-Finite Difference (FD) Method for Diffusion and Reaction-Diffusion Equations on Surfaces. J. Sci. Comput., 2016 Jun 1, vol. 63 (3), pp. 745-768. DOI: 10.1007/s10915-014-9914-1

[20] Chein-Shan Liu, Dongjie Liu. Optimal shape parameter in the MQ-RBF by minimizing an energy gap functional. Applied Mathematics Letters, 2018, vol. 86, pp. 157-165. DOI: 10.1016/j.aml.2018.06.031

[21] Rahimi A., Shivanian C., Abbasbandy S. Analysis of new RBF-FD weights, calculated based on inverse quadratic functions. J. Math., 2022, (2022).

[22] Álvarez D., González-Rodríguez P., Kindelan M. A local radial basis function method for the Laplace-Beltrami operator. J. Sci. Comput., 2021, vol. 86, 28.

[23] Tominec I., Larsson E., Heryudono A. A Least Squares Radial Basis Function Finite Difference Method with Improved Stability Properties. SIAM Journal on Scientific Computing, 2021, vol. 43, no. 2, pp. A1441-A1471.

DOI: 10.1137/20M1320079

[24] Kalani Rubasinghe, Guangming Yao, Jing Niu, Gantumur Tsogtgerel. Polyhar-monic splines interpolation on scattered data in 2D and 3D with applications. Engineering Analysis with Boundary Elements, 2023, vol. 156, pp. 240-250. DOI: 10.1016/j.enganabound.2023.08.001

[25] Великанов П.Г., Артюхин Ю.П. Общая теория ортотропных оболочек. Часть I. Вестник Самарского университета. Естественнонаучн. серия, 2022, т. 28, № 1-2, c. 46-54. DOI: 10.18287/2541-7525-2022-28-1-2-46-54

[26] Великанов П.Г., Артюхин Ю.П. Общая теория ортотропных оболочек. Часть II. Вестник Самарского университета. Естественнонаучн. серия,

2022, т. 28, № 3-4, c. 40-52. DOI: 10.18287/2541-7525-2022-28-3-4-40-52

[27] Jie Hou, Ying Li, Shihui Ying, Iterative optimization method for determining optimal shape parameter in RBF-FD method. Applied Mathematics Letters,

2023, vol. 145, art. ID 108736. DOI: 10.1016/j.aml.2023.108736

[28] Fornberg B., Flyer N., Fast generation of 2D node distributions for mesh-free PDE discretizations. Computers & Mathematics with Applications, 2015, vol. 69, iss. 7, pp. 531-544. DOI: 10.1016/j.camwa.2015.01.009

[29] Shankar V. The overlapped radial basis function-finite difference (RBF-FD) method: A generalization of RBF-FD. J. Comput. Phys., 2017, vol. 342, pp. 211-228.

[30] Красноруцкий Д.А., Лакиза П.А., Шелевая Д.Р. Программный комплекс для моделирования механики системы тонких упругих стержней. Краевые задачи и математическое моделирование: темат. сб. науч. ст. Новокузнецк, Изд-во КГПИ КемГУ, 2023, с. 57-60.

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

Ссылку на эту статью просим оформлять следующим образом: Нгуен К.М., Шелевая Д.Р., Красноруцкий Д.А. Методика расчета осесиммет-ричных колебаний оболочек вращения с жидкостью по дифференциальной модели. Инженерный журнал: наука и инновации, 2024, вып. 2. http://dx.doi.org/10.18698/2308-6033-2024-2-2338

Нгуен Кыонг Мань — аспирант кафедры «Прочность летательных аппаратов» НГТУ. Область научных интересов: численное моделирование механики оболочек. e-mail: mckq1985@gmail.com

Шелевая Дарья Руслановна — мл. науч. сотр., Институт гидродинамики им. М.А. Лаврентьева СО РАН; аспирант и ассистент кафедры «Прочность летательных аппаратов» НГТУ. Область научных интересов: численное моделирование механики оболочек. e-mail: shelevaya.d.r@hydro.nsc.ru

Красноруцкий Дмитрий Александрович — канд. техн. наук, доцент кафедры «Прочность летательных аппаратов» НГТУ; старший научный сотрудник СибНИА им. С.А. Чаплыгина. Область научных интересов: численное моделирование механики стержневых систем и оболочек. e-mail: krasnorutskiy@corp.nstu.ru

Method for calculating axisymmetric oscillations of the shells of revolution with liquid according to the differential model

© CM. Nguyen1, D.R. Shelevaya1,2, D.A. Krasnorutskiy1,3

1 Novosibirsk State Technical University, Novosibirsk, 630073, Russian Federation 2 Lavrentyev Institute of Hydrodynamics SB RAS, Novosibirsk, 630090, Russian Federation 3 Siberian Aeronautical Research Institute named after S. A. Chaplygin, Novosibirsk, 630051, Russian Federation

The paper proposes a technique for calculating dynamic deformation of the axisymmetric orthotropic shells of revolution shells based on solving differential equations by the finite difference method. In particular, it uses the polyharmonic radial basis functions to solve the adjoin problem of the tank insignificant oscillations with the ideal incompressible fluid generating the local finite difference weight coefficients on an arbitrary stencil of the computational domain nodal points. The axisymmetric shell meridian deformation is described by a system of six differential motion equations obtained on the basis of general spatial equations for resolving functions in the global coordinate system. To describe rotations, the final rotation vector (Eulerian vector) is used. Equations do not contain the initial curvature; they describe large rotations, displacements and deformations and take into account thinning/thickening and influence of the transverse shear in the thick shell deformation. The proposed technique is implemented in the proprietary DARSYS software package. Oscillation frequencies of the liquid tanks are calculated using the proposed technique and are compared with frequencies obtained by the other methods. To analyze convergence, the ANSYS system is used, as well as a program that implements the finite and boundary element methods.

Keywords: tank hydroelastic oscillations, ideal incompressible fluid, displacement potential, radial basis functions, finite difference method (RBF-FD)

REFERENCES

[1] Kolesnikov K.S., Rybak S.A., Samoylov E.A. Dinamika toplivnykh sistem ZhRD [Dynamics of the rocket engine liquid fuel systems]. Moscow, Mashinostroenie Publ., 1975, 172 p.

[2] Wei Liu, Chang Xiao, Hao Zhou, Chenyan Wang. Experimental investigation of liquid-tank interaction effects on full containment LNG storage tanks through shaking table tests. Thin-Walled Structures, 2023. https://doi.org/10.1016/j.tws.2023.111527

[3] Gribkov V.A., Adamenko R.A. Dvumernaya model zhidkosti dlya rascheta sobstvennykh chastot kolebaniy osesimmetrichnykh gidroobolochechnykh sis-tem [Two-dimensional fluid model for calculating the natural vibration frequencies of axially symmetric hydro-shell systems]. Inzhenerny zhurnal: nauka i in-novatsii — Engineering Journal: Science and Innovation, 2017, iss. 3 (63), p. 5. https://doi.org/10.18698/2308-6033-2017-3-1593

[4] Shupikov A.N., Misyura S.Yu., Yareshchenko V.G. Chislennoe i eksperimen-talnoe issledovanie gidroupruigikh kolebaniy obolochek [A numerical and experimental study of hydroelastic shell vibrations]. Vostochno-Evropeyskiy zhur-nal peredovykh tekhnologiy — Eastern-European Journal of Enterprise Technologies, 2014, vol. 6, no. 7 (72), pp. 8-12. https://doi.org/10.15587/1729-4061.2014.28861

[5] Jhung M.J., Jo J.C., Jeong K.H. Modal analysis of conical shell filled with fluid. J. Mech. Sci. Technol, 2006, vol. 20, no. 11, pp. 1848-1862.

[6] Goncharov D.A., Pozhalostin A.A., Kokushkin V.V. Modelirovanie osesimmetrich-nykh kolebaniy uprugogo baka s zhidkostyu s uchetom sil poverkhnostnogo natya-zheniya posredstvom mekhanicheskogo analoga [The mechanical of the small axisymmetric oscillations of the with the surface tension forces in elastic tank]. Nauka i obrazovanie: nauchnoe izdanie MGTU im. N.E. Baumana — Science and Education: Scientific Edition of Bauman MSTU, 2015, no. 6, pp. 372-383. https://doi.org/10.7463/0615.0779724

[7] Phan Hoang Nam, Fabrizio Paolacci. Fluid-structure interaction problems: An application to anchored and unanchored steel storage tanks subject to seismic loadings. arXiv: Numerical Analysis, 2018, n. pag.

[8] Shklyarchuk F.N., Rei Juhnbum F. Ray Zhongbum. Raschet neosesimmetrich-nykh kolebaniy obolochek vrashcheniya s zhidkostyu metodom konechnykh elementov [Calculation of non-axisymmetric vibrations of shells of revolution with liquid by finite element method]. Vestnik MAI - Aerospace MAI Journal, 2013, vol. 20, no. 2, pp. 49-58.

[9] Shklyarchuk F.N. Raschet kolebaniy obolochek vrashcheniya s zhidkostyu metodom konechnykh elementov [Calculation of oscillations of shells of revolution with liquid using the finite element method]. Problemy mashinostroeniya

i nadezhnosti mashin — Journal of Machinery Manufacture and Reliability, 2015, no. 1, pp. 17-29.

[10] Mokeev V.V. Issledovanie dinamiki konstruktsiy s zhidkostyu i gazom s pomoshchyu metoda konechnykh elementov [Investigation of the dynamics of structures with liquid or gas by the finite element method]. Izv. RAN. Mekhanika tverdogo tela — Mechanics of Solids. A Journal of Russian Academy of Sciences, 1998, no. 6, pp. 166-174.

[11] Yi Lu, Ji Jing. Modal analysis on anchored tank considering shell and fluid coupling. Advanced materials research, 2012, vol. 549, pp. 903-907. Crossref, https://doi.org/10.4028/www.scientific.net/amr.549.903

[12] Dubois T.T., de Rouvray A.L. An improved fluid superelement for the coupled solid-fluid-surface wave dynamic interaction problem. Earthquake Eng. Struct. Dynam, 1978, vol. 6, no. 3, pp. 235-245.

[13] Gnitko V.I., Degtyarev K.G., Kononenko E.S., Tonkonozhenko A.M. Sravnenie metodov konechnykh i granichnykh elementov v zadachyakh o kolebaniyakh sostavnoy obolochka vrashcheniya s zhidkostyu [Comparison of finite and boundary element methods in problems of vibrations of a composite shell of revolution with the fluid]. Visnik Kharkivskogo natsionalnogo universitetu imeni V.N. Karazina — Bulletin of V.N. Karazin Kharkiv National University, 2019, iss. 42, pp. 38-45.

[14] Levin V.E. Metod konechnykh i granichnykh elementov v dinamike konstruktsiy letatelnykh apparatov: 05.07.03 "Prochnost i teplovye rezhimy letatelnykh ap-paratov. Dis. ... d-ra tekhn. nauk [Finite and boundary element method in the dynamics of aircraft structures: 05.07.03 "Strength and thermal conditions of aircraft". Diss. ... Dr. Sc. (Eng.)]. Novosibirsk, Novosibirsk State Technical University, 2001, 341 p.

[15] Bochkarev S.A., Lekomtsev S.V., Matveenko V.P. Sobstvennye kolebaniya usechennykh konicheskikh obolochek, soderzhashchikh zhidkost [Natural vibrations of truncated conical shells containing fluid]. Prikladnaya matematika i mekhanika — Journal of Applied Mathematics and Mechanics, 2022, vol. 86, no. 4, pp. 505-526. https://doi.org/10.31857/S0032823522040038

[16] Nurul Izyan M.D., Viswanathan K.K., Aziz Z.A., et al. Free vibration of layered truncated conical shells filled with quiescent fluid using spline method. Compos. Struct, 2017, vol. 163, pp. 385-398.

[17] Mohammadi N., Aghdam M.M., Asadi H. Instability analysis of conical shells filled with quiescent fluid using generalized differential quadrature method. In: The 26th Annual Int. Conference of Iranian Society of Mechanical Engineers-ISME2018, April 24-26, 2018. School of Mechanical Engineering, Semnan Univ., Semnan, Iran, ISME 2018-1216.

[18] Flyer N., Fornberg B., Bayona V., Barnett G.A. On the role of polynomials in RBF-FD approximations: I. Interpolation and accuracy. Journal of Computational Physics, 2016, vol. 321, pp. 21-38. https://doi.org/10.10167j.jcp.2016.05.026

[19] Shankar V., Wright G.B., Kirby R.M., Fogelson A.L. A Radial Basis Function (RBF)-Finite Difference (FD) Method for Diffusion and Reaction-Diffusion Equations on Surfaces. J Sci Comput., 2016 Jun 1, vol. 63 (3), pp. 745-768. https://doi.org/10.1007/s10915-014-9914-1

[20] Chein-Shan Liu, Dongjie Liu. Optimal shape parameter in the MQ-RBF by minimizing an energy gap functional. Applied Mathematics Letters, 2018, vol. 86, pp. 157-165. https://doi.org/10.1016Zj.aml.2018.06.031

[21] Rahimi A., Shivanian C., Abbasbandy S. Analysis of new RBF-FD weights, calculated based on inverse quadratic functions. J. Math., 2022.

[22] Álvarez D., González-Rodríguez P., Kindelan M. A local radial basis function method for the Laplace-Beltrami operator. J Sci Comput, 2021, no. 86, p. 28.

[23] Tominec I., Larsson E., Heryudono A. A Least Squares Radial Basis Function Finite Difference Method with Improved Stability Properties. SIAM Journal on Scientific Computing, 2021, vol. 43, no. 2, pp. A1441-A1471. https://doi/org/10.1137/20M1320079

[24] Kalani Rubasinghe, Guangming Yao, Jing Niu, Gantumur Tsogtgerel. Polyhar-monic splines interpolation on scattered data in 2D and 3D with applications. Engineering Analysis with Boundary Elements, 2023, vol. 156, pp. 240-250. https://doi/org/10.1016/j.enganabound.2023.08.001

[25] Velikanov P.G., Artyukhin Yu.P. Obshchaya teoriya ortotropnykh obolochek. Chast I [General theory of orthotopic shells. Part I]. Vestnik Samarskogo uni-versiteta. Estestvennonauchn. Seriya — Vestnik of Samara University. Natural Science Series, 2022, vol. 28, no. 1-2, pp. 46-54. https://doi.org/10.18287/2541-7525-2022-28-1-2-46-54

[26] Velikanov P.G., Artyukhin Yu.P. Obshchaya teoriya ortotropnykh obolochek. Chast II [General theory of orthotopic shells. Part II]. Vestnik Samarskogo uni-versiteta. Estestvennonauchn. Seriya — Vestnik of Samara University. Natural Science Series, 2022, vol. 28, no. 3-4, pp. 40-52. https://doi.org/10.18287/2541-7525-2022-28-3-4-40-52.

[27] Jie Hou, Ying Li, Shihui Ying. Iterative optimization method for determining optimal shape parameter in RBF-FD method. Applied Mathematics Letters, 2023, vol. 145, art. ID 108736. https://doi.org/10.1016/j.aml.2023.108736

[28] Fornberg B., Flyer N. Fast generation of 2D node distributions for mesh-free PDE discretizations. Computers & Mathematics with Applications, 2015, vol. 69, iss. 7, pp. 531-544. https://doi.org/10.1016/j.camwa.2015.01.009

[29] Shankar V. The overlapped radial basis function-finite difference (RBF-FD) method: A generalization of RBF-FD. J Comput Phys, 2017, vol. 342, pp. 211-228.

[30] Krasnorutskiy D.A., Lakiza P.A., Shelevaya D.R. Programmnyi kompleks dlya modelirovaniya mekhaniki sistemy tonkikh uprugikh sterzhney [A software package for modeling the mechanics of a system of thin elastic rods]. In:

Kraevye zadachi i matematicheskoe modelirovanie: temat. sb. nauch. st. [Boundary value problems and mathematical modeling: thematic collection of scientific articles]. Novokuznetsk, KGPI KemGU Publ., 2023, pp. 57-60.

Nguyen Cuong Manh, Postgraduate, Department of Aircraft Strength, Novosibirsk State Technical University (NSTU). Scientific interests: numerical simulation of the shell mechanics. e-mail: mckq1985@gmail.com

Shelevaya D.R., Junior Researcher, Lavrentyev Institute of Hydrodynamics SB RAS; Postgraduate and Assistant, Department of Aircraft Strength, Novosibirsk State Technical University (NSTU). Scientific interests: numerical simulation of the shell mechanics. e-mail: shelevaya.d.r@hydro.nsc.ru

Krasnorutskiy D.A., Cand. Sc. (Eng.), Associate Professor, Department of Aircraft Strength, Novosibirsk State Technical University (NSTU); Senior Researcher, Siberian Aeronautical Research Institute named after S.A. Chaplygin. Scientific interests; numerical simulation of rod systems and shells mechanics. e-mail: krasnorutskiy@corp.nstu.ru