Содержание к диссертации
Введение
Глава 1. Математическая модель взаимосвязанного течения жидкости и движения системы профилей с кинематическими и упругими связями 14
1.1. Общая постановка задачи 14
1.2. Моделирование течения несжимаемой жидкости методом ВВД
1.2.1. Вычисление конвективной скорости по заданному распределению завихренности 17
1.2.2. Вычисление диффузионной скорости в методе ВВД 19
1.3. Граничные условия 20
1.3.1. Выражение граничного условия прилипания через распределение завихренности 21
1.3.2. Формулировка граничных условий в дискретном виде 22
1.4. Моделирование движения тел с кинематическими и упругими связями 25
1.4.1. Гидродинамические силы, действующие на твёрдое тело 25
1.4.2. Уравнения движения твёрдых тел с упругими связями под действием гидродинамических сил 27
1.4.3. Уравнения связей в общем виде
1.5. Моделирование гибкого профиля 30
1.6. Выводы по главе 1 37
Глава 2. Программный комплекс «Vvflow» 38
2.1. Общая характеристика программного комплекса 38
2.2. Библиотека метода ВВД 38
2.3. Прикладные программы 39
2.4. Выводы по главе 2 43 Стр.
Глава 3. Методические исследования и тестирование вычислительной схемы ВВД 44
3.1. Исследование схемной вязкости метода ВВД 44
3.2. Устойчивость численной схемы метода ВВД 49
3.3. Исследование сходимости на примере обтекания неподвижного кругового профиля 50
3.4. Выводы по главе 3 52
Глава 4. Верификация программного комплекса, примеры решения нестационарных задач взаимодействия тел с жидкостью 53
4.1. Стабилизация следа за круговым профилем, совершающем высокочастотные вращательные колебания 53
4.2. Перестройка дорожки Кармана в дальнем следе за круговым профилем 56
4.3. Исследование переходных процессов при отклонении интерцептора на крыловом профиле 61
4.4. Колебания одиночного и тандема физических маятников в вязкой жидкости 72
4.5. Автономное движение квазибиологического объекта в вязкой среде 77
4.6. Моделирование машущего крылового профиля в потоке вязкой несжимаемой жидкости 80
4.7. Выводы по главе 4 85
Общие выводы и заключение по диссертационной работе 88
Список литературы
- Вычисление конвективной скорости по заданному распределению завихренности
- Прикладные программы
- Исследование сходимости на примере обтекания неподвижного кругового профиля
- Исследование переходных процессов при отклонении интерцептора на крыловом профиле
Вычисление конвективной скорости по заданному распределению завихренности
Ось профиля связана с державкой, характеризуемой координатами Р = {РХ}РУ} и собственным углом поворота Р0 (Рис. 1.2).
Профиль имеет три степени свободы, для каждой из которых определены параметры вязко-упругой связи — коэффициенты жёсткости кх, ку, к0, и демпфирования тх, ту, т0. При упругой связи на профиль действует сила FP, приложенная к оси, и момент сил Мр: FPx = -кхАх - тхих, FPy = -куАу - туиу, МР = -к0А0 - т0ш} (1.1) где кх, ку, к0 — коэффициенты упругости, тх, ту, т0 — коэффициенты демпфирования, А = О — Р. Здесь и далее ввиду двумерности задачи все силы приводятся в смысле погонного значения. Некоторые степени свободы могут отсутствовать, при этом полагается к = оо, и соответствующие уравнения (1.1) заменяются кинематическими соотношениями Рх = Ох, Ру = Оу} Ро = 00. (1.2)
Державка может быть неподвижной, или совершать заданное движение, или быть связанной с другим телом. В случае крепления державки профиля с индексом (ш) к опорному профилю в точке Р(т, на последний, согласно третьему закону Ньютона, действуют те же самые силы и момент с противополож у к к Аж , __ р4 Jo o кх Ту ткуHWWU\h р д„ F LizEJ Л о оТх У////////, О х Рис. 1.2. Принципиальная схема упругого закрепления тела на державке ными знаками. В таком случае уравнения движения опорного тела имеют вид ти + mezuj хйm- muj2R m = FP + FH + Fg - Fpm), Іегш 1mхгаЙ=ег {Мр + Мя + Мд) - ezMpm) (Р{т) б х F Здесь т — масса профиля, / — момент инерции вокруг оси О, Rm — радиус-вектор центра масс, Rm = (Rm — О), F, М — силы и моменты, действующие на тело: Fp, Мр — реакция опоры; F#, МЦ — гидродинамические нагрузки; Fg,Mg — внешние силы и их момент (сила тяжести). Все моменты вычисляются относительно оси вращения профиля. 1.2. Моделирование течения несжимаемой жидкости методом ВВД Метод вязких вихревых доменов (ВВД) [1] предназначен для решения двумерных нестационарных уравнений Навье-Стокса, описывающих течения вязкой несжимаемой жидкости: DV z/W--Vp, VV = 0. (1.3) Р Dt р Здесь = + (V V) — оператор субстанциональной производной, V = V(f,t) — скорость жидкости, p = p(f,t) — давление, v = const — кинематический коэффициент вязкости, р = const — плотность жидкости.
Аналогично методу DVM (diffusion velocity method) [68], в методе ВВД используется понятие диффузионной скорости Vd = -i/Vn/П, где ezQ = V х V (вектор завихренности в плоскопараллельных течениях направлен перпендикулярно плоскости течения). Это позволяет записать уравнение (1.3) в виде V {V + Vd)xezn = --v( p + — }. (1.4) dt " p 2 После применения оператора ротор к обеим частям уравнения (1.4) получается уравнение эволюции завихренности V х Vz х ezQ , Vz = V + Vd, ez— dt котороев скалярной форме в случае плоскопараллельных течений представляет собой уравнение переноса завихренности [78]: = -V (Vsri) . (1.5)
Из этого уравнения видно, что вектор JY, = VY1 имеет смысл потока завихренности. Его составляющие J = VQ и Jj = V fil = — uVQ представляют собой конвективный и диффузионный потоки завихренности соответственно.
Уравнение (1.5) позволяет рассматривать поле Q как плотность некой материальной субстанции и моделировать вихревые области в виде мелких участков (доменов), которые движутся со скоростью VY,. Как следует из уравнения (1.5), при таком движении у каждого домена сохраняется величина 7v = Js Qds, где Sv — площадь домена. Как известно, интеграл завихренности по площади S равен циркуляции скорости V по контуру, ограничивающему эту площадь. Поэтому величину 7v называют циркуляцией домена. Домен также можно рассматривать как частицу, характеризуемую значением интенсивности (циркуляции) и положением в пространстве, которое определяется радиус-вектором контрольной точки Ёv внутри домена.
Общая схема метода ВВД аналогична схеме метода DVM. Вихревые частицы, распределение которых задано в начальной момент, перемещаются со скоростью VY с постоянной циркуляцией 7v При этом на обтекаемых поверхностях образуются новые частицы, циркуляции которых определяются заданными граничными условиями. Существенным отличием метода ВВД от DVM является способ вычисления диффузионной скорости, позволивший адекватно моделировать эволюцию завихренности вблизи поверхностей.
Скорость жидкости V выражается через поле завихренности и скорость на границе области течения с помощью обобщённого разложения Гельмгольца, называемого также формулой Био-Савара: V(R) = 00+ Q(r)ez х K(R )ds + (п х V(f)j х K(R )dl+ S с + ( n- V(r)j K(Rr)dl, K(Rr) = ——2, Rf = R-f. с
Здесь S, С — область течения и её граница соответственно, V(r) — скорость жидкости на поверхности, п — вектор единичной нормали, направленный внутрь области течения. В дискретном представлении интеграл по пространству заменяется суммированием по доменам, а интегралы по контуру — суммами по отрезкам, аппроксимирующим профили: Nv -, Ns V(R) = Voo + У 7vгг/(Лv) + 2,1 j (fij x Vsjj x K(R sJ)+ i=1 j=1 + Y)j (4 Kj ) K(R sj), (16) 3=1 где Ns — количество отрезков, где Nv — количество вихревых доменов, / , Rsj, Vsj — длина j-го отрезка, радиус-вектор его середины и скорость соответственно, v — скорость, индуцируемая единичным вихревым доменом, R!vi = R — Rvi, R sj = R — Rsj. Если выполняется условие прилипания, то на поверхности справедливо равенство V(Rs) = Rs.
Выражения (ftj х Vsj) и {п3 Vsj), стоящие под знаками сумм, удобно рассматривать как поверхностные вихри и источники соответственно: 1-п- х Vsj = ez7sj, Ifij Vsj = qsj. (1.7) При этом выражение (1.6) примет вид V(R) = ї4 + J A + Хл ( ) + Qsjv(R sJ) х ez. (1.8) i=1 j=1 Обычно для нахождения индуцированной скорости v в вихревых методах используется выражение ІЄ(Л ), \й\ є tf(R ) = gz x { 1 , (1.9) , lit I є 2тг є2 что соответствует равномерному распределению завихренности внутри круга радиуса є и отсутствию завихренности вне круга. Такое распределение носит название вихря Рэнкина. В данной работе используется другое выражение: 1 R v(Rf) (1.10) 2тг R 2 + є2
Прикладные программы
Описанный в предыдущей главе алгоритм был разработан в виде комплекса программ. Одной из задач разработки было максимально широко охватить все доступные классы задач. Помимо этого, от программной реализации требовалось предоставить возможность обрабатывать результаты, без которой «сырые» данные не представляют научного интереса. Все эти требования были учтены, и итогом работы явился набор программ, который не только обеспечивает решение различных задач гидродинамики, но и включает в себя ряд полезных утилит пре- и пост-обработки, упрощающих проведение численного эксперимента.
Комплекс предназначен для операционных систем GNU/Linux. Большая часть кода написана на языке C++, в утилитах используется также Python, Awk, и Bash.
Весь программный комплекс можно условно разделить на две части: библиотеку метода ВВД, и набор прикладных программ.
Самой главной компонентой библиотеки является структура TSpace, которая хранит в себе все данные о расчёте. Он же реализует функции сохранения и загрузки результатов. В классе TSpace содержится следующая информация: - служебная информация: название расчёта, время директивного останова, реальное время; - физические параметрыпространства: вязкость жидкости, скорость набегающего потока, текущий момент времени и параметры дискретизации; - массив вихревых доменов: их координаты, циркуляции и скорости; - контуры профилей и их параметры: плотность, координаты державки, деформации вязко-упругих связей и их коэффициенты жёсткости и демпфирования; - бинарное дерево, обеспечивающее быстрое решение задачи N тел [81]. Остальные компоненты библиотеки реализуют формулы, приведённые в главе 1. К ним относятся следующие модули: - заполнение системы линейных уравнений, описанной в разделе 1.4; - решение этой системы; - генерация новых вихревых доменов; - объединение вихревых доменов для уменьшения вычислительной нагрузки; - построение бинарного дерева для быстрого решения задачи N тел; - вычисление размера ядра вихря , описанного в разделе 1.2.1 - вычисление конвективной скорости по формуле (1.6); - вычисление диффузионной скорости по формуле (1.12); - интегрирование уравнений движения доменов и тел.
Среди прикладных программ самой главной является vvflow — она реализует основной цикл алгоритма ВВД [82]. Помимо пошагового решения задачи программа vvflow отвечает за загрузку данных и сохранение результатов расчёта. Программа получает имя входного файла в виде единственного аргумента командной строки и загружает параметры расчёта. На их основе инициализируются необходимые модули из библиотеки метода. После инициализации всех модулей запускается цикл по времени согласно алгоритму метода. На каждом шаге по времени сохраняются интегральные параметры течения, а текущее состояние расчёта сохраняется с периодичностью, указанной в параметрах. По достижении момента окончания расчёта программа завершает свою работу, и можно переходить к обработке результатов.
Для хранения параметров расчёта и полученных результатов используется формат файлов HDF5 [83]. Он предназначен для хранения бинарных данных в иерархическом виде, а также поддерживает встроенное сжатие. Для создания файла с параметрами расчёта в составе комплекса имеется препроцессор vvcompose. Он обладает консольным интерфейсом и позволяет легко редактировать любые параметры класса TSpace. Наиболее трудоёмким является процесс задания контуров: для каждого профиля необходимо задать массив точек, его аппроксимирующих. Для упрощения работы используется вторая утилита препроцессора vvgen, которая генерирует профили различной формы по заданным геометрическим параметрам.Наданный момент поддерживаютсяэл-липс, пластина и другие формы профилей, рассматриваемые в главе 4. Список поддерживаемых форм по необходимости расширяется.
Для обработки результатов применяются программы vvxtract и vvplot. Первая предназначена для чтения формата HDF5 и используется при полуавтоматической обработке результатов скриптами. Назначение программы vvplot — визуализация полученных результатов. Примеры вариантов визуализации приведены на Рис. 2.1: массив вихревых доменов, поле завихренности, поле давления и линии тока. Так как результаты численного моделирования не содержат в себе ничего кроме массива вихревых доменов и информации о профилях, все необходимые вычисления также делает постпроцессор.
Непосредственное создание изображений выполняется свободной программой gnuplot [84]. Назначение программы vvplot заключается в подготовке численных данных в необходимом формате. Массив вихревых доменов и точек обтекаемых профилей в подготовке не нуждаются. Поля завихренности и давления вычисляются непосредственно в процессе визуализации на основе информации о вихревых доменах. Внимания заслуживает визуализация линий тока, для которой предложен новый подход.
Исследование сходимости на примере обтекания неподвижного кругового профиля
На Рис. 4.12 представлены полученные зависимости от времени, начиная с момента начала подъёма интерцептора. Результаты численного моделирования сравниваются с двумя работами: [111] — результат вычислительного эксперимента, и [ 110] — физический эксперимент. В обеих этих работах конфигурация тела и постановка задачи полностью совпадают с рассматриваемыми, за исключением числа Рейнольдса, которое в этих работах было равно 3.5 105, в то время как в расчёте методом ВВД использовалось Re = 105.
На обоих графиках также изображён закон открытия интерцептора 6(t). Как видно из Рис. 4.12, результаты показывают качественное согласие с данными [110; 111].
В начальный момент движения интерцептора аэродинамические характеристики изменяются медленно, затем следует значительное увеличение подъёмной силы и пикирующего момента, после чего устанавливается колебательный режим. Причём в эксперименте амплитуда колебаний затухает, тогда как в расчётных работах этого не наблюдается. Последнее обстоятельство, возможно, связано с двумерностью расчётов. В трёхмерном течении после прекращения движения интерцептора крупномасштабные вихревые структуры начинают ослабляться, передавая энергию мелкомасштабным. В двумерных течениях этот механизм отсутствует.
Рис. 4.12. Эволюция коэффициента подъёмной силы и момента тангажа в процессе открытия интерцептора. Чёрная линия — результат данной работы Из графиков видно, что колебания характеристик, полученныхна основе уравнений Рейнольдса в работе [111], имеют меньшую амплитуду и меньший период, чем в эксперименте, что связано с использованием модели турбулентности. В этом отношении результаты метода ВВД ближе к экспериментальным.
На Рис. 4.13 представлены расчётные зависимости от времени суммарных аэродинамических коэффициентов сопротивления и подъёмной силы на крыле с интерцептором. Текущее положение вихревых доменов позволяет получить мгновенное поле давления и мгновенные линии тока в различные моменты времени. Результаты отражены на Рис. 4.14. В левой колонке изображены мгновенные линии тока на фоне поля коэффициента давления (области повышенного давления более светлые, пониженного — тёмные). Справа представлены поля завихренности (тёмный цвет соответствует областям отрицательной завихренности, светлый — положительной).
Сопоставление полей давления и завихренности позволяет проанализировать динамику вихрей и её взаимосвязь с изменением нагрузок на профиль. Моменты времени, соответствующие рисункам полей давления и завихренности с Рис. 4.14, сопоставлены с мгновенными значениями Cxa(t), Cya(t) на Рис. 4.13.
В начале процесса, когда угол открытия интерцептора мал, он работает аналогично изменению кривизны профиля: увеличивается скорость на верхней
Изменение коэффициентов аэродинамического сопротивления и подъёмной силы в процессе открытия интерцептора. Цифрами отмечены моменты времени, соответствующие Рис. 4.14 Рис. 4.14. Эволюция поля давления (слева) и завихренности (справа) в процессе открытия интерцептора Рис. 4.14. (продолжение) Эволюция поля давления (слева) и завихренности (справа) в процессе открытия интерцептора поверхности, циркуляция и подъёмная сила (момент 1). На коэффициенте сопротивления Сха эти изменения в начальные моменты почти не сказываются.
Далее между интерцептором и поверхностью профиля образуется область пониженного давления, приводящая к искривлению линии тока, огибающей интерцептор, и стеканию с него отрицательной завихренности, так как именно такой знак завихренности имеет пограничный слой на верхней стороне профиля. Накапливаясь за интерцептором, этот вихрь индуцирует на верхней стороне поверхности профиля скорость, направленную от задней кромки к ин-терцептору, что приводит к образованию на ней пограничного слоя с положительной завихренностью (изображено белым цветом). Эта завихренность сносится под интерцептор и накапливается там, отодвигая отрицательный вихрь к задней кромке (момент 2). Этот момент характеризуется максимальной подъёмной силой и максимальным пикирующим моментом из-за пониженного давления. Образовавшаяся пара вихрей засасывает жидкость вначале вдоль поверхности профиля, затем между вихрями и поступает под продолжающий открываться интерцептор между положительным вихрем и новым отрицательным, стекающим с интерцептора. По мере заполнения пространства под интер-цептором первый отрицательный вихрь сносится вниз по течению, частично увлекая с собой положительный вихрь (момент 3). Коэффициенты сопротивления и подъёмной силы уменьшаются.
При прохождении вихря над задней кромкой профиля искривляется линия тока, сходящая с его нижней поверхности, и вблизи задней кромки на верхней поверхности накапливается положительная завихренность (момент (4). Подъёмная сила достигает минимума. Коэффициент сопротивления при этом близок к локальному максимуму, так как новый отрицательно закрученный вихрь создаёт область пониженного давления за кромкой интерцептора.
Далее отрицательный вихрь вытесняет положительный, давление на верхней стороне профиля за интерцептором опять падает и подъёмная сила вновь возрастает (момент (5). Позже можно наблюдать аналогичные процессы поочерёдного образования положительных и отрицательных вихрей, сопровождающиеся колебанием аэродинамических коэффициентов и продолжающихся после остановки интерцептора.
Исследование переходных процессов при отклонении интерцептора на крыловом профиле
При проведении экспериментов [114] по обтеканию эластичного профиля было замечено, что гибкость препятствует нарушению симметрии следа, возникающему в опытах с жёстким профилем. Асимметричные следы в этом случае не наблюдались. При этом область значений числа Струхаля в экспериментах не превышала Sh = 0.45. Численное моделирование трёхзвенной модели показали, что область частот и амплитуд, в которой след симметричен, расширяется по сравнению со случаем жёсткого профиля, однако при достаточно больших значениях частоты и/или амплитуды симметрия всё равно нарушается. Возможно, это связано с недостаточной гибкостью рассмотренной в расчёте модели.
В данной серии численных экспериментов интерес представляет сравнение пропульсивной силы и её отношения к затраченной мощности в случаях жёсткого и гибкого профилей. С этой целью были проведены расчёты однозвен-ного профиля на упругом подвесе (h = К, к2 = h = об) и трёхзвенного профиля с тремя упругими шарнирами (h1 = К,к2 = К/4, к3 = К/16). В качестве основного коэффициента жёсткости использовалась величина К = 800pD2V . В обоих этих случаях задача является сопряжённой, а все секции профиля движутся под действием гидродинамических сил и упругих моментов.
Зависимость пропульсивной силы (тяги) от числа Струхаля для од-нозвенного (жесткого) и трёхзвенного (гибкого) профилей На Рис. 4.34 приведены расчётные зависимости осреднённой по времени погонной пропульсивной силы от числа Струхаля при Re = 250, Ав = 7, р = ръ/р =1,3. Значения силы представлены обезразмеренными на величину pDV. Из графика видно, что во всех случаях при Sh « 0.5 наблюдается локальный максимум тяги, при этом тяга гибкого профиля вновь начинает увеличиваться при Sh 0.9, в то время как в случае жёсткого профиля она стремится к нулю.
На Рис. 4.35 представлены зависимости коэффициента полезного действия г] = FVoo/W от частоты колебаний при тех же параметрах, где W = кіїв — осреднённая мощность, прикладываемая в первом шарнире. Из Рис. 4.35 видно, что к.п.д. гибкого профиля всюду выше, чем у жёсткого, причём при Sh величина г] остаётся приблизительно на неизменном уровне, а у жёсткого профиля КПД резко падает с увеличением Sh.
Выводы по главе Представлены примеры практического применения программного комплекса «Vvflow» для решения исследовательских задач нестационарной гидродинамики. В каждой из рассмотренных шести задач не только протестированы разнообразные вычислительные возможности разработанного комплекса, но и получены новые научные сведения о свойствах нестационарных вихревых течений вязкой жидкости при взаимодействии с твёрдыми телами. Полученные
Зависимость коэффициента полезного действия от числа Струхаля Sh для однозвенного и трёхзвенного профилей расчётные данные сопоставлены с известнымирезультатами физических экспериментов. Тем самым подтверждена способность разработанной вычислительной технологии адекватно отражать сложные нестационарные взаимодействия систем профилей с вязкой жидкостью.
В задаче об обтекании вязкой жидкостью кругового профиля, совершающего высокочастотные вращательные колебания. численно воспроизведён известный из физических экспериментов эффект вырождения дорожки Кармана по мере увеличения частоты колебаний профиля. Идентифицированы тонкие концентрические системы чередующихся вихревых слоёв разноимённой циркуляции, влияющие на количество завихренности, поступающей в ближний след за колеблющимся цилиндром и, в конечном счёте, влияющей на его устойчивость.
В вычислительном эксперименте воспроизведено и объяснено явление перестройки дорожки Кармана в следе за цилиндром. Показано, что область локального понижения скорости на оси следа в неограниченном пространстве связана с неустойчивостью первичной дорожки Кармана к варикозным возмущениям, приводящим к образованию своеобразного вихревого кластера с повышенной дипольной плотностью. Кластер перемещается относительно вихрей Кармана вверх по течению. Это перемещение носит волновой характер, при котором в область сгущения включаются вихри, лежащие выше по потоку, а ниже по потоку остаются деформированные смешавшиеся друг с другом и претерпевшие частичную аннигиляцию вихри. При этом обеспечивается устойчивое положение дипольного кластера относительно цилиндра на неизменном расстоянии, зависящем от числа Рейнольдса.
Решена задача об обтекании аэродинамического профиля NACA-0012 с интерцептором. Воспроизведён наблюдаемый в физических экспериментах парадоксальный положительный «заброс» подъёмной силы и дано его объяснение на основе анализа мгновенных распределений давления в окрестности открывающегося интерцептора. Среднее значение прироста подъёмной силы при любом угле отклонения интерцептора является отрицательным. Величина пика подъёмной силы немонотонно зависит от времени. Угол открытия ин-терцептора, при котором достигается максимум подъёмной силы, монотонно уменьшается с уменьшением угловой скорости интерцептора. Максимальная величина пика наблюдается при времени открытия tA = 4 - 8, критический угол составляет 30 40.
В результате вычислительных экспериментов с одиночным и парой физических маятников, погружённых в вязкую несжимаемую жидкость в однородном поле сил тяжести показано, что энергия одиночного маятника убывает немонотонно; наблюдается рекуперация энергии, связанная с передачей импульса от жидкости маятнику в фазе его подъёма. Декремент затухания энергии зависит в основном только от массы маятника. Также показано, что характер колебаний маятника зависит не только от его физических свойств, но и и от воздействия среды.
Воспроизведено автономное движение квазибиологического объекта «головастик» в вязкой жидкости за счёт периодической деформации его формы. Получены зависимости скорости движения и затрачиваемой мощности от частоты и амплитуды колебаний. Показано, что минимальные энергозатраты при заданной скорости движения достигаются при амплитуде колебаний A opt « 0.227.
В тестовых расчётах обтекания машущих жёстких и гибких многозвенных крыловых профилей в потоке вязкой несжимаемой жидкости воспроизведены наблюдавшиеся в физических экспериментах эффекты перестройки следа при изменении амплитудно-частотных характеристик колебаний. Показано, что гибкость профиля может оказывать существенное влияние на течение жидкости и гидродинамические силы, действующие на профиль. В частности, на режимах существования пропульсивной силы увеличиваются тяга и коэффициент полезного действия движителя.