Взаимодействие упругих конструкций с потоком жидкости

В первой главе представлен обзор современного состояния проблемы взаимодействия упругих конструкций с потоком жидкости.
Экспериментальному определению вибрационных характеристик ВУ посвящены работы С.И. Девнина, В.Н. Зайцева, В.В. Тюшкевича, М.Ю. Григорова.
Расчетные методы определения гидродинамических характеристик (ГДХ) обтекателей ВУ основаны на численном решении уравнений Навье–Стокса, чему посвящены работы О.М. Белоцерковского, А.А. Самарского, И.А. Белова, С.А. Исаева, Ю.В. Гурьева, И.В. Ткаченко, К.Н. Волкова, В.Н. Емельянова и многих других отечественных и зарубежных авторов.
Вопросы расчетного определения характеристик вибрации упругих конструкций рассматриваются в работах С.П. Тимошенко, Ю.А. Шиманского, М.В. Келдыша, К.С. Колесникова, Р. Бишопа и У. Прайса.
Подход, основанный на совместном численном решении уравнений гидродинамики и механики деформируемого твердого тела (FSI) используется в работах Г.Б. Крыжевича, А.А. Родионова, В.А. Коршунова, Д.А. Пономарева, D.J. Piro и K.J. Maki.
Теоретические оценки тяги и КПД машущего крыла представлены в работах В.В. Голубева. Двумерное обтекание крыла потоком идеальной жидкости рассматривается Л.И. Седовым, М.В. Келдышем, М.А. Лаврентьевым, Я.Ц. Фыном, Р.Л. Бисплингхоффом, Х. Эшли и Р.Л. Халфмэном.
Различные варианты конструкции плавникового движителя рассматриваются в работах В.А. Рыжова, Е.П. Носова и С.В. Тарасова. Нестационарному обтеканию тонких крыловых профилей посвящены работы К.В. Рождественского.
Во второй главе рассматриваются вопросы расчетного определения ГДХ ВУ.
Представлен обзор основных подходов к разрешению характеристик турбулентного течения. Для достижения целей диссертационного исследования был выбран RANS–подход, используемый совместно с моделью турбулентности SAS. Он позволяет разрешать отрыв крупномасштабных вихрей и при этом характеризуется приемлемой вычислительной трудоемкостью.
Валидация данного численного метода проводилась на основе сопоставления результатов расчета и экспериментальных данных по обтеканию стержня в виде кругового цилиндра большого удлинения.
Полученные изоповерхности завихренности представлены на рисунке 1. Структура потока неоднородна по длине стержня и является принципиально трехмерной. Поле завихренности показано на рисунке 2. Расчетное значение коэффициента стационарной силы cx st = 0,66 хорошо согласуется с экспериментальным cx st = 0,65.
Рисунок 1 – Изоповерхности завихренности (1/с)
Рисунок 2 – Поле завихренности (1/с) для обтекателя цилиндрической формы
Форма сечений современных выдвижных устройств значительно отличается от круговой. На рисунке 3 показано поле завихренности течения около обтекателя измененной формы, типичной для современных ВУ. Сравнение ГДХ обтекателей цилиндрической и измененной формы представлено в таблице 1. Обтекатель измененной формы отличается более низкими значениями действующих
также малые
представлен аналитический
метод расчета колебаний ВУ.
В первом приближении ВУ может
быть представлено в виде
прямого однородного стержня
переменного сечения. Расчетная
схема показана на рисунке 4.
Примем, что все инерционные
характеристики стержня включают в себя значения присоединенных масс жидкости. Уравнение поперечных колебаний стержня имеет вид:
∂ ∂ ∂ y ∂ y
1 h ∂t ∂x EJ x ∂x m x ∂t q x, t , (1)
Рисунок 3 – Поле завихренности (1/с) для обтекателя измененной формы
В третьей главе
Таблица 1 – ГДХ обтекателей
ГДХ
гидродинамических сил, а способностью демпфировать колебания.
Цилиндрический Измененный обтекатель обтекатель
cx st cy vibr Sh сyα
0,66 0,44 0,43 0,27 0,2 0,22 – 3,92

где m (x) – погонная масса; y – смещение в направлении Y; q (x, t) – интенсивность внешней распределенной нагрузки; E – модуль Юнга; J (x) – момент инерции сечения с абсциссой x; h – коэффициент внутреннего трения.
Рисунок 4 – Расчетная схема
Для определения собственных форм fn и частот колебаний ωn рассматривается случай отсутствия внешней распределенной нагрузки q (x, t) = 0. Общее решение (1) можно представить в виде:

y x, t f x q t .
(2) (3)
(4)
(5)
Подставим (2) в (1):
EJ x f ω m x f 0,
Граничные условия имеют вид:
f 0 0 EJ x f | 0
Введем замену:
u f x , u f x , u EJ x f x , u EJ x f x ,
с учетом которой (3) сводится к системе:
u 0100 du u ⎛ 0 0 1 0⎞
, f L 0 . f L 0
A ω,x u , u u , A ω,x ⎜ EJ x ⎟. (6) dx u 0001
n 1, 2, … .
⎝ω m x 0 0 0⎠
Общее решение имеет вид:
u C u C u C u C u . (8)
С учетом (4) получаем:
C u L C u L 0 . (9)
C u L C u L 0
Система (9) имеет нетривиальное решение в случае
D ω u L u L 0. (10)
u L u L
Вычисляя значения D при различных ω, можно найти корни этой функции – ωn. Это позволяет найти fn:
f x C u x u ω , L u x . (11) u ω , L
При наличии опоры в сечении Ln стержень разбивается на участки «a» и «b»: u C u C u C u C u C u C u C u C u , (12)
u C u C u C u C u C u C u C u C u . (13)
Примем, что реакция опоры пропорциональна перемещениям с коэффициентом k. В сечении x = Ln выполняются условия:
u L u L , u L u L , u L u L , u L u L ku L .
(14)
Учитывая (14) и наличие наделки массой mnd на свободном конце, находим выражения для fn:
f x C u x m ω u x C u x , (15) f x C u x C u x C u x C u x . (16)
Упругая конфигурация ВУ под действием стационарной нагрузки описывается уравнением:
d d z ρv
dx EJ x dx q x c x 2 d . (17)
где qst (x) – погонная стационарная нагрузка; c – коэффициент гидродинамического сопротивления сечения; v – скорость потока; dH – характерный размер сечения.
Система, аналогичная (6), имеет вид:
0100 0 du ⎛0010⎞ 0
d x A u B , A ω , x ⎜ E J x ⎟ , B 0 . ( 1 8 )
0001 q x ⎝0 0 0 0⎠

10 Если на наделку действует сила Pndst, то:
u L 0, u L 0, u 0 0, u 0 P . (19) Решение (18) с учетом (19) имеет вид:
u C u C u C u C u u , (20)
u C u C u C u C u u , (21)
где u 0 0 0 0 0 .
Определив Cn, получим упругую конфигурацию ВУ и распределение силовых
факторов.
Нестационарные гидродинамические силы могут быть представлены в виде
q (x, t) = q (t)ꞏP (x):
q hq EJ x f m x f q q t P x . (22)
(23)

m q 2 ε q ω q Q ,
k
m m x f dx, 2ε hm ,
k
ω m , k EJ x f dx , Q q t P x f dx.
Если гидродинамическая нагрузка действует по гармоническому закону, то:
q 2ε q ω q Решение (24) имеет вид:
m
cosωt cosωt, m
n 1,2,….
P x f dx Q
(24) (25)
q A cos ωt φ , A Q ,
n 1,2,….
где
Разработанный метод расчета колебаний ВУ в потоке жидкости был реализован в виде комплекса программ, позволяющего определять такие характеристики вибрации
φ arctg 2ε ω . (26) ω ω
y x, t f x A cos ωt φ . (27)
m ω ω 4ε ω
Окончательное выражение для вынужденных колебаний можно записать как:

ВУ, как распределения максимальных значений перемещений и напряжений по длине ВУ. Главное окно программного комплекса представлено на рисунке 5.
В четвертой главе описан расчет колебаний упругого стержня в потоке жидкости.
Рассматривались малые колебания прямого однородного стержня в потоке жидкости. Данная задача состоит из гидродинамической части, в результате решения которой были определены ГДХ заделанного стержня, а также динамической части, в которой на основе данных о гидродинамических воздействиях были получены реакции конструкции – перемещения и напряжения. Для решения динамической части задачи использовались разработанный аналитический метод, а также численное интегрирование уравнений механики деформируемого твердого тела с помощью метода конечных элементов (МКЭ). Результаты, полученные с помощью обоих подходов, хорошо согласуются между собой и показаны на рисунке 6.
Рисунок 5 – Главное окно программы для расчета характеристик вибрации ВУ
Одним из способов борьбы с вихревой вибрацией ВУ является применение гасителя гидродинамического типа. При внесении струи в вихревой след за конструкцией наблюдается уменьшение амплитуд гидродинамических сил вследствие стабилизации следа. Поля завихренности при различных значениях отношений скорости струи v1 к скорости потока v0 представлены на рисунках 7 и 8. На рисунке 9 показан вариант реализации обтекателя с гидродинамическим гасителем.

Рисунок 6 – Максимальные значения перемещений, полученных при помощи разработанного метода и МКЭ
Результаты расчета колебаний стержня с гасителем подтверждают его эффективность в борьбе с вихревой вибрацией. Рассматривался цилиндрический стержень, взаимодействующий с потоком жидкости при неработающем (v1 = 0) и работающем (v1/v0 = 2,5) гасителе. На рисунке 10 представлены распределения максимальных значений перемещений по длине стержня при различных состояниях гасителя. Использование данного способа борьбы с вихревой вибрацией позволяет уменьшить максимальные значения перемещений в 3 раза.
Рисунок 8 – Поле завихренности (1/с) при v1/v0 = 5
Рисунок 7 – Поле завихренности (1/с) при v1/v0 = 2,5
Рисунок 9 – Обтекатель, оснащенный гидродинамическим гасителем вибрации
Рисунок 10 –Максимальные значения перемещений стержня
В пятой главе рассматривается численный метод расчета ГДХ колеблющегося в потоке жидкости тела.
В качестве объекта рассматривалось прямоугольное крыло малого удлинения. Экспериментальная модель, показанная на рисунке 11, совершала вращательные колебания в воздушном потоке.
Расчетная сетка, построенная для данной задачи, показана на рисунке 12. Ее внутренняя сферическая подобласть совершала колебания вместе с крылом относительно внешней неподвижной подобласти. Экспериментальное амплитудное значение коэффициента подъемной силы cyA = 0,27 хорошо согласуется с расчетным значением cyA = 0,26.
Расчеты для водной среды проводились для частот n = 1,14 Гц и n = 4,56 Гц. Скорость потока составляла 1,57 м/с. Полученное значение производной коэффициента подъемной силы по углу атаки составило c = 3,12.
Рисунок 11 – Экспериментальная модель Рисунок 12 – Сечение расчетной сетки
ГДХ колеблющегося в потоке жидкости крыла могут быть получены при помощи теории квазиустановившегося обтекания (ТКО) крыла. На рисунках 13 и 14 приведены расчетные и теоретические зависимости коэффициента cy(t). При n = 1,14 Гц наблюдается хорошее согласование теории и численного расчета. При n = 4,56 Гц заметна роль нестационарных эффектов. Это приводит к расхождениям в численных и теоретических результатах. Поля скорости при n = 1,14 Гц и n = 4,56 Гц представлены на рисунках 15 и 16.
Рисунок 13 – Сравнение теоретических и расчетных зависимостей cy (t)
(n = 1,14 Гц)
Рисунок 14 – Сравнение теоретических и расчетных зависимостей cy (t)
(n = 4,56 Гц)

Таблица 2 – Амплитудные значения cy и mz
Рисунок 15 – Поле скорости (м/с) при n = 1,14 Гц
Рисунок 16 – Поле скорости (м/с) при n = 4,56 Гц
Учесть влияние нестационарных эффектов (силы инерции присоединенных масс жидкости, прирост сил циркуляционной природы, возникающих из-за внезапного увеличения скоса потока) при определении ГДХ колеблющегося в потоке крыла позволяет теория нестационарного движения (ТНД) тонкого крыла.
В таблице 2 представлены амплитудные значения коэффициентов подъемной силы cy и продольного вращательного момента mz. При медленных колебаниях крыла ТКО дает приемлемые результаты. В случае существенно нестационарных движений результаты численного расчета будут находится в диапазоне, нижняя граница которого определяется с помощью ТКО, а верхняя – на основе ТНД. Полученные численно значения cy и mz являются физически обоснованными, что подтверждает адекватность вычислительных методов, применяемых в расчетах с использованием динамических сеток.
Метод расчета
ТКО Численный расчет ТНД
cy
mz
n = 1,14 Гц 0,58 0,58 0,80
n = 4,56 Гц 0,95 1,45 1,75
n = 1,14 Гц 0,07 0,08 0,11
n = 4,56 Гц 0,18 0,30 0,40
Шестая глава посвящена численному определению пропульсивных характеристик СПД.
СПД представляет собой
систему крыловых поверхностей,
последовательно соединенных
жесткими связями. Каждая связь
имеет в составе шарнир c
торсионной пружиной c жесткостью k. Схема простейшего СПД показана на рисунке 17. Закон движения ведущего элемента 0 известен. Элемент 1 совершает колебания, вызванные кинематическим возбуждением элемента 0.
Рисунок 17 – Схема простейшего СПД
В качестве обобщенной координаты выберем угол отклонения хорды элемента 1 от горизонтали φ. Начало системы координат расположено в шарнирной точке. На зависимый элемент действует гидродинамическая сила, которая представима в виде проекций Fx и Fz. Здесь и далее предполагается, что центры масс элементов совпадают с их центрами давления.
Уравнение движения элемента 1 имеет вид:
φ L cosφ mAω cos ωt F F sinφ kφ , (30)
J mL L
где L – расстояние от шарнира до центра масс элемента 1; m – масса элемента 1; J – момент инерции элемента 1; A – амплитуда колебаний элемента 0; ω – частота колебаний элемента 0. Зависимости Fx (t) и Fz (t) определялись в процессе совместного численного интегрирования уравнения движения элемента 1 и уравнений динамики вязкой жидкости.
Эффективность различных режимов движения СПД оценивалась на основе сопоставления значений таких пропульсивных характеристик, как коэффициент силы тяги J и КПД η:
J J 2ρv B∑ b ,
1
J T F dt,
n 0,1, (31) (32)
N J v v F dt η
,
N N
F ∙ Z dt
где bn – хорда элемента n, T – период колебаний, Fx и Fz – проекции гидродинамической силы на оси x и z для системы в целом.
Рассматривались режимы с частотой колебаний f ведущего элемента от 0,5 до 0,8 Гц. На рисунке 18 показаны зависимости φ(t). Максимальная амплитуда колебаний достигается при рабочей частоте 0,7 Гц. Также были рассчитаны значения коэффициента силы тяги J и КПД η для каждого режима. Они представлены в таблице 3. Режим, соответствующий рабочей частоте 0,7 Гц, характеризуется максимальной эффективностью.
Таблица 3 – Пропульсивные характеристики СПД
f,Гц J η

0,5 0,6 0,7 0,8
-0,003 -0,08 0,014 0,34 0,025 0,49
-0,0002 -0,01
Рисунок 18 – Зависимость φ (t)
Седьмая глава посвящена вопросам вывода уравнений движения СПД.
Расчетная схема для рассматриваемого СПД, состоящего из трех крыловых поверхностей, показана на рисунке 19.
элементами имеет вид:
m L J φ J φ T 2 5φ 4φ φ cos φ φ φ 2 2 .
Уравнения динамики
Рисунок 19 – Расчетная схема
СПД могут быть получены с помощью подходов Лагранжа и Виттенбурга. При этом применение аппарата уравнений Лагранжа II рода сопряжено с трудностями при решении практических задач, поскольку требует работы с весьма громоздкими выражениями. Так, уравнение для кинетической энергии СПД с двумя зависимыми
С увеличением количества машущих элементов количество членов в выражении для кинетической энергии растет в геометрической прогрессии.
и матрицу T:
1, S 1, 0,
u выходит из S , u входит в S ,
в других случаях,
u ∈ путимеждуS иS ⎪и направлена к S ,
(34)
⎧ 1,
T ⎨ 1, u ∈путимеждуS иS (35) ⎪и направлена от S ,
⎩ 0, u ∉путимеждуS иS .
Положение центра тяжести тела системы определяет вектор cia. Также вводятся вектора dij:

T S , i,j 1…n.
(36)
Введем понятие дополненного тела (ДТ). В смежную для выбранного тела шарнирную точку присоединяется масса всех тел (кроме тела 0), связанных через нее
(33)
Альтернативный подход Виттенбурга основан на схематизации системы в виде направленного графа, показанного на рисунке 20. Структура взаимосвязей графа определяет единственным образом матрицу с элементами Sia:
Рисунок 20 –Ориентированный граф
с этим телом напрямую или косвенно. Положение центра тяжести ДТ определяется вектором bi0, а выражение для тензора инерции имеет вид:

m , i 1 … n. (37)

Уравнения движения тел системы можно записать в виде:
φ φ φ
A ⋮ B ⋮ SkS φ⋮ R k φ SS , i 1…n. (38)
φ φ
В выражении (38) k – диагональная матрица коэффициентов жесткости пружин
в шарнирах. Элементы матриц A, B и R имеют вид:
⎧K ,
A ⎪M d b cos φ φ α β ,
⎨⎪M d b cos φ φ α β , ⎩0,
M d b sin φ φ α β ,
B M d b sin φ φ α β , 0,
i j, s s ,
s s ,
в других случаях,
s s ,
s s ,
в других случаях,
R M M b r sin φ β r cos φ β
d F sin φ α F cos φ α . :
Углы αij и βi определяют положение векторов dij и bi0 в связанном с телом базисе. Mi – момент внешних сил, Fj1 и Fj2 – проекции внешней силы на оси 1 и 2.
В восьмой главе исследуются колебания СПД в потоке жидкости. Рассматривается СПД, состоящий из четырех машущих элементов. Движение трех зависимых элементов описано с помощью системы уравнений, аналогичной (38), полученной с помощью подхода Виттенбурга. Линеаризованные уравнения динамики СПД были использованы для исследования устойчивости движения системы.
При варьировании различными параметрами задачи относительно базового режима анализировались корни характеристического уравнения. На рисунке 21 показано изменение действительной части одного из корней при варьировании скоростью потока v0. При v0 > 0,7 м/с действительная часть становится положительной, а движение системы – неустойчивым.
Моделирование колебаний СПД в потоке жидкости было проведено в два этапа. Вначале элемент 0 оставался неподвижным, а элементу 3 придавалась некоторая начальная угловая скорость. Режимы колебаний показаны на рисунках 22 и 23.
Затем исследовались режимы движения, когда элемент 0 совершал вертикальные колебания по гармоническому закону, при этом варьировалась частота колебаний. Анализ зависимости максимальных амплитуд зависимых элементов от рабочей частоты, представленный на рисунке 24, показывает, что максимальная амплитуда достигается при частоте, соответствующей основной собственной частоте системы.
Рисунок 21 – Зависимость действительной части корня от скорости потока v0
Рисунок 22 – Устойчивый режим при v0 = 0,3 м/с
Рисунок 23 – Неустойчивый режим при v0 = 0,9 м/с
Рисунок 24 – Зависимость максимальных амплитуд от частоты колебаний элемента 0
ЗАКЛЮЧЕНИЕ
В рамках диссертационного исследования были разработаны математические модели и численные методы решения важных классов задач о вибрации ВУ и колебаниях СПД в потоке жидкости.
Основные результаты диссертационного исследования заключаются в следующем.
1. Разработаны математическая модель и метод численного расчета вибраций ВУ. Для решения гидродинамической части данной задачи были рассмотрены и проанализированы основные подходы к разрешению турбулентного течения. Применение выбранного численного метода позволяет достичь цели исследования с минимальными вычислительными затратами, при этом результаты расчетов хорошо согласуются с экспериментальными данными. Полученные гидродинамические характеристики обтекателей используются в качестве исходных данных для решения динамической части задачи. На основе дифференциальных уравнений колебаний стержня переменного сечения был разработан аналитический метод, позволяющий оперативно определять параметры напряженно-деформируемого состояния ВУ, учитывая при этом влияние промежуточной опорной связи. Разработанный комплекс программ, позволяющий определять вибропрочностные характеристики ВУ при различных скоростях хода подводного объекта, может быть использован в проектно- конструкторской деятельности.
2. Разработанный расчетный метод определения вибрационных характеристик ВУ был применен для моделирования колебаний упругого стержня в потоке жидкости. Достоверность полученных результатов подтверждает хорошее согласование с данными, полученными с помощью численного решения уравнений механики деформируемого твердого тела методом конечных элементов. Разработанный метод позволил произвести оценку эффективности гасителя вибрации гидродинамического типа. Согласно результатам расчетов, применение гидродинамического гасителя предложенной конструкции позволяет значительно снизить параметры вихревой вибрации ВУ и улучшить их эксплуатационные качества.
3. На основе сопоставления с экспериментальными данными и аналитическими оценками была подтверждена адекватность численного метода расчета ГДХ колеблющегося в потоке жидкости крыла. Данный метод использовался в дальнейшем для моделирования работы СПД. С помощью совместного численного интегрирования уравнений динамики вязкой жидкости и уравнений движения машущих элементов были получены гидродинамические и пропульсивные характеристики данного движительного комплекса. Анализ средней тяги и КПД позволил выявить наиболее эффективный режим движения для рассматриваемого варианта конструкции.
4. Рассмотрены основные подходы к выводу уравнений движения элементов СПД. Показано, что классический подход Лагранжа является более трудоемким по сравнению с подходом Виттенбурга, при этом оба подхода позволяют получить одинаковые уравнения. Даже в случае относительно простой конструкции (четыре последовательно соединенных крыла) выражение для кинетической энергии становится настолько громоздким, что дальнейшие операции с ним при использовании подхода Лагранжа вызывают затруднения, также возникает риск появления ошибок при проведении дальнейших преобразований. С помощью подхода Виттенбурга была разработана математическая модель движения элементов СПД, что позволило провести дальнейшее исследование его колебаний в потоке жидкости. Линеаризованные уравнения использовались для проведения анализа устойчивости
движения конструкции. Варьирование таких параметров, как скорости потока, коэффициенты жесткости пружин, коэффициенты демпфирования в шарнирах, позволило выявить устойчивые режимы. Также было проведено моделирование малых колебаний СПД. Анализ различных режимов вынужденных колебаний показал, что максимальная амплитуда машущих элементов достигается при частоте колебаний ведущего элемента, близкой к основной собственной частоте.
Таким образом, в результате выполнения диссертационного исследования была решена научная задача о разработке и обосновании математических моделей и методов расчета вибраций ВУ и колебаний СПД в потоке жидкости. Разработанные численные методы могут быть использованы при проведении расчетов, необходимых для проектирования данных упругих конструкций.
Предлагаемый расчетный метод определения характеристик вибраций ВУ позволит сократить объем экспериментальных исследований и, тем самым, снизить стоимость их создания. В ходе разработки метода была подтверждена адекватность подходов, используемых при решении гидродинамической и динамической частей гидроупругой задачи. Это позволяет моделировать колебания ВУ в сопряженной постановке, что может быть целесообразно при рассмотрении резонансных режимов, характеризующихся большими амплитудами колебаний.
Методы численного определения гидродинамических и пропульсивных характеристик СПД позволяют обоснованно подходить к вопросу выбора его конструктивных параметров. Выведенные с помощью подхода Виттенбурга уравнения движения элементов многосоставных конструкций могут быть использованы для дальнейшего совместного численного интегрирования с уравнениями динамики вязкой жидкости. Это позволит приблизиться к реализации регулируемого гидроупругого эффекта и достичь повышенных пропульсивных характеристик, что приведет к улучшению эксплуатационных характеристик подводного объекта, оснащенного подобным типом движителя.