ПРИКЛАДНАЯ МАТЕМАТИКА И МЕХАНИКА
Том 68. Вып. 1, 2004
УДК 539.3:534.1
© 2004 г. Р. Е. Лампер, В. Е. Левин
МЕТОДЫ КОНЕЧНЫХ И ГРАНИЧНЫХ ЭЛЕМЕНТОВ В ЗАДАЧАХ ДИНАМИКИ УПРУГИХ СОСУДОВ С ЖИДКОСТЬЮ
Излагается вариант метода расчета собственных продольных колебаний упругого осесимметричного сосуда, частично заполненного жидкостью. В основу метода положено сочетание граничных элементов для жидкости и конечных элементов для упругого сосуда. Приводятся примеры, показывающие эффективность предлагаемого метода.
Основные постановки задач о динамическом взаимодействии тонкостенного упругого сосуда с заполняющей его жидкостью обсуждались ранее [1, 2]. Предложенные алгоритмы расчета динамических характеристик продольных колебаний сосуда, ориентированные на использование ЭВМ, в определенной степени решили часть проблемы для сосудов простых форм [3, 4]. Дальнейшее их развитие связано с прогрессом вычислительной техники [5, 6].
Развитие конструктивных схем ракет-носителей требует совершенствования расчетных алгоритмов. Имеющиеся вычислительные комплексы типа КЛЗТЯЛК не всегда позволяют получить приемлемое решение.
Предлагаемый метод построен на более гибком описании геометрии сосуда с использованием построенных эффективных конечных элементов, что позволяет уточнить результаты расчетов по известным методикам и проводить расчеты сосудов сложной геометрии с учетом их внутренних конструктивных особенностей.
1. Постановка задачи. Рассматривается частично заполненный несжимаемой невязкой жидкостью тонкостенный упругий сосуд, составленный из оболочек вращения, подкрепленных шпангоутами. Движение жидкости считается потенциальным; пренебрегается энергией волнообразования на свободной поверхности жидкости.
Свяжем с сосудом цилиндрическую систему координат г, г, Ф (фиг. 1). Внешнее крепление сосуда реализуется через шпангоуты с помощью кольцевых подвесок, обладающих распределенной жесткостью в продольном и радиальном направлениях, а также угловой жесткостью. В расчетной схеме жесткие связи представляются достаточно большими жесткостями в нужных направлениях. Внутри сосуда могут находиться стержневые конструкции, не влияющие на течение жидкости. К шпангоутам могут быть прикреплены дополнительные массы. Рассматриваются малые гармонические продольные колебания сосуда с круговой частотой ю (все рассматриваемые функции имеют временной множитель ехр(гюО).
Такая конструкция с жидкостью представляет собой систему с бесконечным числом степеней свободы. Для ее приближенного описания системой с конечным числом степеней свободы ниже используются два известных метода: метод конечных элементов (МКЭ) для представления упругих элементов сосуда и метод граничных элементов (МГЭ) для моделирования движений жидкости в сосуде.
При постановке задачи для оболочек в рамках МКЭ имеются разные возможности выбора аппроксимирующих функций формы. Требования к таким функциям известны [7]: эти функции должны восстанавливать перемещение элемента как жесткого целого, а при уменьшении размеров элемента они должны описывать состояние постоянной деформации. Известны проблемы, связанные с построением
Фиг. 1
криволинейного элемента, обладающего такими свойствами [8, 9]. Отдельным вопросом является описание исходной геометрии объекта в виде, удобном для эффективного применения ЭВМ, в частности для упрощения процедуры увеличения количества конечных элементов.
Метод граничных элементов предъявляет особые требования к функциям, аппроксимирующим границу жидкого объема, и к функциям формы, описывающим перемещения этой границы.
2. Описание меридиана осесимметричного сосуда. Параметрические уравнения меридиана - плоской кривой, наиболее просто выглядят в случае, когда параметром является длина дуги кривой s. В рамках такого описания можно приближенно построить участок кривой по информации в соседних узлах m и m + 1 следующим образом [10].
Радиус-вектор участка кривой запишем в виде
г (s) = xk (s) ik, k = 1, 2 где ik - орты глобальной системы координат.
Производная вектора r по дуге имеет единичную длину. Обозначим dx1/ds = cos х, dx2/ds = sinx, где x - угол наклона нормали к кривой. Зададим закон изменения угла х по длине кривой в виде
Х( s)
1- Vх"
s
+ - а«
1+^ sj s
(2.1)
где хт, Хт + 1 - углы в точках т и т + 1 соответственно, В и I - подлежащие определению неизвестный коэффициент и неизвестная длина кривой между точками. Для неизвестных В и I получим систему двух нелинейных уравнений
cos
1- sJX»
+ -Х"
+1
+ BI 1-
s\s
l l.
ds
v2m + 1
= x2 m +1 Sln
Í Sil
1 - + - Xm + 1
+ в 11-
s\ s
ljl
ds
x1 m + 1
= x 1 m +
0
Для решения системы применяется метод Ньютона.
Поскольку в аппроксимации (2.1) заложена возможность описания постоянной кривизны кривой, то по изложенному алгоритму с высокой степенью точности восстанавливается элемент окружности любого угла раствора, и в частности прямолинейный участок.
Таким образом, криволинейный контур сосуда восстанавливается на ЭВМ по координатам отдельных точек и направлениям нормалей к контуру в этих точках. В процессе аппроксимации вычисляются длины криволинейных участков.
3. Описание движения жидкости по МГЭ. Принятые допущения относительно движений жидкости в объеме V являются основанием для введения потенциала перемещений жидкости ф [11]. Потенциал ф удовлетворяет краевой задаче
Дф = 0 в V; ф = 0 на ЦП = и ■ п на
где и - вектор упругого перемещения стенки сосуда, п - вектор внешней нормали к стенке, ^ - свободная поверхность жидкости, - смоченная поверхность сосуда.
Для формулировки задачи о колебаниях сосуда с жидкостью нужно найти связь между граничными значениями ф и Эф/Эя. Такая связь между соответствующими значениями в конечном числе точек смоченного контура сосуда устанавливается в рамках прямого метода граничных элементов (ПМГЭ).
Реализация ПМГЭ основана на применении интегральной теоремы Грина для двух гармонических в области V функций: искомого потенциала ф и фундаментального решения уравнения Лапласа для осесимметричного течения жидкости. Фундаментальное решение задается в таком виде, чтобы получаемое по ПМГЭ решение для потенциала перемещений удовлетворяло условию на свободной поверхности жидкости точно. В рассматриваемой задаче фундаментальное решение имеет вид
Ф(х, х(= О(х, х(- О(х, х(
( 2\ ( 2\
О (х, х(0) = — К 1--2 , О (х, х( °) = -1- К 1 - —2
Ш1 Т2 П Т* Т*2
1 V '1/ 1 V Т1 у
где К(...) - полный эллиптический интеграл первого рода, х - текущая точка границы области, х(1) - точка коллокации на границе, т - расстояние между точками х и х(1), т1 - расстояние между точкой х и точкой х'(1), симметричной точке х(1) относительно оси симметрии области, занятой жидкостью, т* и т* - расстояния между точкой х и
точками, симметричными точкам х(1) и х'(1) относительно поверхности жидкости.
В формуле Грина интегрирование проводится по смоченной поверхности. Схема получения разрешающих соотношений ПМГЭ стандартна при использовании граничных элементов с постоянным распределением неизвестных граничных функций. В этом случае точки коллокации берутся в середине каждого элемента, а для согласования гранично-элементной и конечно-элементной схем на каждом конечном элементе берется по два граничных элемента. Сингулярные интегралы вычисляются по специальной квадратурной формуле, содержащей логарифмы. В результате разрешения системы уравнений ПМГЭ относительно неизвестных узловых значений потенциала ф функция кинетической энергии жидкости (кинетическая энергия без множителя ю2 - квадрата круговой частоты) сводится к квадратичной форме неизвестных узловых значений Эф/Эп, которые, в свою очередь, выражаются через узловые значения обобщенных перемещений конечно-элементной схемы сосуда.
При использовании в качестве Эф/Эп соответствующих функций формы смоченного конечного элемента оболочки проблемы стыковки МКЭ и МГЭ схем не возникает.
Применение этого подхода позволяет записать функцию кинетической энергии жидкости в виде квадратичной формы обобщенных узловых неизвестных конечно-элементной схемы
1 т
Т = ± Ут Т г У
2 1
где Т - матрица присоединенных масс жидкости, У - вектор обобщенных узловых неизвестных конечно-элементной схемы.
4. Конечно-элементное описание оболочки. В рассматриваемой задаче форма оболочек сосуда до и после его деформации определяется плоской кривой - меридианом. Построение конечного элемента такого сосуда эквивалентно построению элемента для плоского криволинейного стержня. Функции формы конечного элемента плоского стержня предлагается получать в результате интегрирования задаваемых изменений углов поворота нормали к стержню и деформации вдоль стержня [12]. Связь деформации е, угла поворота нормали к элементу Ах и производных глобальных компонентов вектора перемещений и1, и2 по дуге ■ контура сосуда имеет вид
йик / = ейхк / + (-1 )кАхйхъ_кШ, к = 1, 2 (4.1)
Эти соотношения сохраняют свой вид и при другой параметризации осевой линии стержня. Их можно положить в основу построения функций формы.
Возможны различные варианты представления деформации и угла поворота конечного элемента. Здесь предположим, что неизвестная деформация е постоянна на элементе. Угол поворота нормали к элементу при его деформировании зададим в виде
Ах( ■) = Ах/1-Ь{з)) + Ах;- +ХЬ{8) + В( 1-Ь{з))Ь{8); Ь{з) = ^^ (4.2)
sj +1 -
где j и j + 1 - узлы конечного элемента, АХj, АХj + 1 - углы поворота нормалей к элементу в этих узлах, Sj + 1, Sj - координаты узлов элемента вдоль дуги, В - постоянная, подлежащая определению вместе с е.
Обозначим: иу, и у - компоненты перемещения узла j, и^ + 1, и у + 1 - компоненты перемещения узла j + 1.
Интегрируя равенства (4.1) с использованием соотношения (4.2), получим
' йХ
ик(■) = Vк + е[Хк(■) - х^ + (-1 )к|Ах(у)-^¿У, к = 1, 2
Постоянные В и е определяются из условий в узле j + 1
ик(+!) = ик} + !, к = 1, 2
и выражаются через обобщенные узловые перемещения конечного элемента.
В частном случае прямого стержня описанный подход приводит к известным функциям формы [7].
Приведем пример эффективности построенного конечного элемента. Рассматривался плоский стержень с осево
Для дальнейшего прочтения статьи необходимо приобрести полный текст. Статьи высылаются в формате PDF на указанную при оплате почту. Время доставки составляет менее 10 минут. Стоимость одной статьи — 150 рублей.