Санкт-Петербургский государственный университет
Механика и математическое моделирование
Прикладная аэродинамика
Селиванов Никита Андреевич
Определение аэродинамических характеристик прямоугольной призмы
с помощью пакета программ SU2.
Бакалаврская работа
Научный руководитель:
Гл. науч. сотр., д. ф.-м. н., с. н. с. Рябинин А.Н.
Рецензент:
Доц. каф. физики и химии, к. т. н., Ежов О. Н.
Санкт-Петербург
2016
SAINT-PETERSBURG STATE UNIVERSITY
Mechanics and mathematical modeling
Applied aerodynamics
Selivanov Nikita
Evaluation of aerodynamic rectangular plate characteristics with the aid of
SU2 software package.
Bachelor’s Thesis
Scientific supervisor:
Principal Researcher, D.Sc. Ryabinin Anatoly
Reviewer:
Associate prof. of the department of Ph. and Chem., PhD, Ezhov Oleg
Saint-Petersburg
2016
Содержание
1. Введение. О необходимости изучения аэродинамических
характеристик призмы……………………………………………………..4
2. Математическая постановка задачи. Уравнения НавьеСтокса……….................................................................................................6
3. Осреднение уравнений по Рейнольдсу…………………………………...6
4. Модель турбулентности Спаларта-Аллмараса…………………………..8
5. Построение геометрической модели призмы и сетки с помощью
программы Gmsh………………………………………………………….10
6. Моделирование обтекания воздухом с помощью программы SU2 и
построение конфигурационного файла………………………………….12
7. Визуализация полученных результатов с помощью программы
Paraview……………………………………………………………………14
8. Сравнение полученных данных с экспериментальными данными……18
9. Изучение возможности галопирования для призмы данной формы….19
10. Заключение………………………………………………………………..20
11. Список используемой литературы………………………………………21
Введение.
Изучение аэродинамических характеристик призмы сегодня имеет большое
значение: призмы прямоугольной или квадратной формы очень часто
используются в машиностроении, при строительстве мостов, зданий и
многих других конструкций. Известны случаи, когда из-за сильного ветра
происходили разрушения различных объектов, в частности, обрушение
мостов.
Для примера, рассмотрим самый известный случай: крушение Такомского
висячего моста в США 7 ноября 1940 года. Крушение произошло через 3
месяца после запуска моста при относительно небольшой скорости ветра 19
м/с (до этого мост выдерживал и более сильные ветровые нагрузки).
Мост начал раскачиваться под действием горизонтального ветра, и
постепенно, изгибные колебания перешли в изгибно-крутильные:
4
Когда закрутка достигла максимума, угол проезжей части к горизонту
достигал 45 градусов. После этого самые слабые элементы моста – подвески
не выдержали, и мост рухнул. Таким образом, причиной разрушения стало
одновременное явление аэроупругого флаттера (сочетание незатухающих
изгибающих и крутящих автоколебаний) и галопирования. Процесс
разрушения был заснят на 16-миллиметровую цветную киноплёнку, что
позволило в дальнейшем детально изучить причины разрушения.
Этот случай послужил важным толчком для изучения аэродинамических
характеристик призмы. Но аналогичные проблемы возникают до сих пор: в
2010 году амплитуда колебаний волгоградского моста составляли 1 метр.
Галопирование – это автоколебания упругой системы в ветровом потоке
(аэроупругие колебания). Оно характерно для гибких сооружений с особыми
формами поперечного сечения, например такими, как прямоугольные или Dобразные. При галопировании в таких сооружениях возможны колебания с
большими амплитудами в перпендикулярном потоку направлении (в 10 или
даже в значительно большее число раз превышающими размеры самого
сечения в этом направлении) при частотах, которые значительно ниже частот
срыва вихрей, характерных для того же самого сечения.
Многие учёные изучали галопирование плохообтекаемых тел. Г. В.
Паркинсон и Н. П. Брукс изучали обтекание призмы с квадратным сечением,
а М. Новак исследовал призму с прямоугольным сечением. Значительный
вклад в исследование обтекания цилиндров внёс Алонсо. Указанные выше
авторы рассматривали тела с бесконечным удлинением.
Изучать галопирование для тел различной формы можно в аэродинамической
трубе. Но создание таких объектов – трудоёмкий процесс. Оптимальный с
5
экономической точки зрения вариант - моделирование обдувания призмы на
компьютере.
Я буду проводить исследование обтекания призмы численным методом с
помощью пакета программ SU2. «Stanford University 2» - бесплатная
программа с открытым исходным кодом, разрабатываемая кафедрой
аэронавтики и астронавтики стэнфордского университета.
Но программа SU2 имеет ограниченные возможности, в частности в ней нет
возможности строить модель объекта (сетку) и визуализировать результаты.
Поэтому сетку я строил с помощью программы Gmsh, а визуализировать
результаты буду с помощью программы ParaView.
Математическая формулировка задачи. Уравнения Навье-Стокса.
В своей работе я буду искать решение с помощью уравнения Навье-Стокса,
осреднённых по Рейнольдсу, при расчёте буду использовать модель
турбулентности Спаларта-Аллмараса.
Уравнения Рейнольдса (в англоязычной литературе: RANS (Reynoldsaveraged Navier–Stokes)) мы можем вывести из уравнений Навье-Стокса[1]
(где u i – компонента скорости, p – давление, ρ – плотность, μ – коэффициент
динамической вязкости, τ ij= ρui’uj’ – тензор турбулентных (Рейнольдсовых)
напряжений) .
Осреднение по Рейнольдсу
6
Уравнения Навье Стокса преобразуются [1] с помощью процедуры
осреднения:
Где f – осредняемая функция, t – время, 2T– период осреднения, который
берём намного большим, чем временные масштабы всех турбулентных
неоднородностей, присутствующих в рассматриваемом течении, и
достаточно малым по сравнению с характерным временным масштабом
осредненного течения.
После применения процедуры осреднения для сжимаемого совершенного
газа уравнения Навье-Стокса могут быть представлены в следующем виде
(знаки осреднения для простоты опущены):
Здесь u – вектор скорости осреднённого течения с компонентами u, v и w,
τm и τt – молекулярная и турбулентная составляющие тензора вязких
напряжений,
E=CvT + ½*(u2+v2+w2) – полная энергия газа,
H=E + ρ/p =CpT + ½*(u2+v2+w2) – его полная энтальпия,
qm и qt – молекулярная и турбулентная составляющие вектора плотности
теплового потока,
T – температура,
Cv=(Cp – R/m) – удельная теплоёмкость газа при постоянном объёме,
Cp – удельная теплоёмкость газа при постоянном давлении,
R = 8,314 Дж/(моль*К) – универсальная газовая постоянная,
7
m – молярная масса газа.
Величины молекулярных составляющих тензора напряжений и вектора
плотности теплового потока определяются соответственно с помощью
реологического закона Ньютона и закона Фурье:
Где S – тензор скоростей деформации, I – единичный тензор, μ(T) и λ(T) –
коэффициенты молекулярной динамической вязкости и теплопроводности.
Модель турбулентности Спаларта-Аллмараса.
В полученной системе уравнений число неизвестных больше количества
уравнений. Так как соотношение между турбулентными составляющими
тензора напряжений tτ и вектора плотности теплового потока t с
параметрами осредненного течения, являющимися основными
переменными этой системы, неизвестно. Мы можем определить это
соотношение с помощью дополнительных соотношений, которые собственно
и составляют модель турбулентности Спаларта-Аллмараса (SA модель) [2].
Эта модель была опубликована в 1992 году. Она достаточна проста в
использовании и идаёт результаты с высокой степенью точности, поэтому SA
модель стала очень популярной среди прочих моделей турбулентности [3]. В
данной модели решается одно дифференциальное уравнение:
Модифицированная кинематическая турбулентная вязкость ν,
относительно которой записано это уравнение, и “истинная”
турбулентная вязкость v t=μt/ρ, входящая в выражение для Рейнольдсовых
напряжений, связаны соотношением: v t=fv1 ν, где
а ν – молекулярная вязкость газа.
8
Генерационный и деструктивный члены в правой части P ν и Dν определяются
следующим образом:
Где d w – ближайшее расстояние до твёрдой стенки, k – постоянная Кармана.
Величина Ω представляет собой модуль тензора завихрённости:
а функция ft2 , обеспечивающая подавление так называемого
“спонтанного” или “численного” перехода от ламинарного режима
течения в пограничном слое к турбулентному, определяется выражением
Наконец, последний член в правой части, предназначенный для
инициирования ламинарно-турбулентного перехода в заданной точке (в
случае двумерных течений) или на заданной линии (в случае трехмерных
течений) на обтекаемой поверхности (“tripterm”), рассчитывается с помощью
формул:
9
Здесь нижний индекс “trip” относится к величинам, определяемым в той
точке на линии перехода, которая находится на минимальном расстоянии
от рассматриваемой точки потока, Δl t = (Δx2 + Δz 2)1/2 - диагональ ячейки
сетки на обтекаемой поверхности в этой точке, а величина ΔU – это
модуль разности векторов скорости в рассматриваемой точке и в точке
перехода.
Эмпирические константы модели равны:
Построение геометрической модели призмы и сетки с помощью
программы Gmsh.
Для изучения обтекания призмы необходимо построить её модель. Но задать
границы призмы будет недостаточно, нужно построить так называемую
сетку, которая будет описывать каждую точку на поверхности с
определённой степенью приближения. Для этого я использовал программу
Gmsh.
Полученный на выходе файл будет иметь расширение .geo
Вот образец команд, которые используются для построения модели в данной
программе:
b=1; /* задаём размер переменной * /
delta=0.025;
/* задаём размер переменной * /
cl1=S rt(b^2+h^2)/100; /* задаём характеристическую длину около
точек * /
Point(1) = {0,-b1,0,cl1}; /* задаём положение первой точки, всего 8
(для призмы). В скобках пишется id, первые 3 значения – координаты
точек, а последнее значение –размер сетки в данной точке. * /
Line(1) = {1,2};
/* задаём линию между первой и второй точкой * /
LineLoop(1) = {-1,14,5,-15};/* по заданным линиям создаём контур * /
PlaneSurface (1) = {1};/* по заданному контуру создаём поверхность * /
Аналогичные команды используются для создания внешней области [4].
10
Я построил объект, длина, ширина и высота которого измеряются в условных
единицах программы Gmsh (высота равна 0,025, а длина и ширина равны 1
единице):
В свою очередь призма расположена в области 50*50 единиц, в которой я
буду моделировать обдувание:
Процесс построения сетки называется триангуляцией. Элементами разбиения
в двумерном случае являются треугольники. Разбиение должно быть более
мелким вблизи нашей призмы, так как именно возле её границ возможно
появление больших градиентов искомых функций. Программа Gmsh
11
автоматически генерирует сетку. После построения сетки получается
следующая картина:
Данная сетка является неструктурированной, то есть узлы сетки
располагаются в произвольном порядке. Но при генерировании сетки
программа использует алгоритм, основанный на критерии Делоне, согласно
которому треугольники стремятся к равносторонним, и соседние
треугольники должны быть приблизительно равны по площади.
После триангуляции программа Gmsh сохраняет сетку во многих форматах
известных конечно-элементных пакетов, мы будем сохранять её в файле с
расширением .su2. В полученном файле записывается положение каждого
узла сетки. Количество полученных узлов: 22822.
Моделирование процесса обдувания с помощью программы SU2.
Кроме файла сетки в формате .su2 и файла с геометрической моделью в
формате .geo для моделирования процесса обдувания призмы с помощью
программы SU2 нам нужно создать ещё конфигурационный файл в формате
.cfg, в котором будут указаны все основные параметры, такие как скорость
потока, температура, плотность воздуха, числа Маха и Рейнольдса и так
далее.
Образец конфигурационного файла можно скачать на сайте разработчиков
программы.
Я буду пользоваться программой SU2 версии 2.0.1.0. Обратите внимание на
то, что сейчас уже выпущены более новые версии программы: третья и
12
четвёртая, но я использовал именно вторую версию, так как она наиболее
изучена. Мой конфигурационный файл подходит только для второй версии
программы, при запуске его на других версиях возникают ошибки.
Конфигурационный файл разбит на 16 разделов:
1) Постановка задачи. Здесь мы говорим, что будем работать с
уравнением Навье-Стокса, и добавляем, что будем использовать
модель Турбулентности Спаларта-Аллмараса.
2) Задаём параметры набегающего сжимаемого потока: число Маха =
0,059, угол атаки, температура = 293 К, число Рейнольдса Re = ρVD/η =
1,3*106 (где D –диаметр, η – коэффициент вязкости).
3) Задаём константы сжимаемого и несжимаемого потока. Например,
здесь задаём газовую постоянную равной 287,87.
4) Определение эталонных значений. Например, плотность: 1.2886 Кг/м3,
давление: 101325.0 Н/м2, и так далее.
5) Задаём граничные условия.
6) Общие параметры, определяющие численный метод. Например, задаём
условие Куранта-Фридрихса-Леви или условия для метода РунгеКутты.
7) Определение линейного расчёта.
8) Параметры сетки.
9) Определение численного метода потока.
10)
Определение численного метода турбулентности.
11)
Стратегия разделения на секции.
12)
Стратегия адаптации сетки.
13)
Параметры деформации сетки.
14)
Параметры сходимости.
15)
Входные / выходные даны. Здесь мы задаём путь к файлу сетки,
определяем формат сетки, задаём названия файлов на выходе и
уточняем, что файлы на выходе должны быть построены для
программы Paraview.
16)
Определение оптимальной формы объекта.
Я взял следующие начальные данные:
Скорость набегающего потока: 19,5 метров в секунду.
Температура: 20 градусов по Цельсию.
И провёл исследование для 11 разных углов атаки от -20 до 20 градусов.
13
В процессе вычислений приближенные значения интересующих нас
коэффициентов устанавливаются не сразу, а примерно после 15000 – го шага
программы. Вычисление для каждого угла атаки происходит относительно
долго, около 6 часов. Для удобства программа создаёт файл restart_flow.dat, с
помощью которого можно продолжить вычисление в случае
непредвиденного завершения работы программы. После 20000 шага
интерации я прерывал расчёты, так как значения коэффициентов подъёмной
силы и сопротивления становятся практически постоянными. На выходе
получаются файлы с детальным отображением всех величин для каждого
шага программы с расширением .csv (его открыть можно как с помощью
Paraview, так и с помощью Microsoft Excel) и файлы для визуализации
результатов с расширением .vtk. Хотелось бы отметить, что результаты
вычислений можно визуализировать также с помощью программ TECPLOT и
STL. Но я буду использовать для визуализации программу Paraview.
Все программы, которые я использовал бесплатны. Таким образом всё это
делает изучение гидроаэромеханики доступным для широкого круга
общественности. Кроме этого, они не требуют мощных и современных
компьютеров, для любой из этих программ достаточно обычного ПК (чего
нельзя сказать про ANSYS).
Визуализация результатов с помощью программы Paraview.
Итак, приступим к визуализации результатов.
Вот так будет выглядеть распределение плотности около призмы при угле
атаки 0 градусов:
14
А вот так будет выглядеть распределение числа Маха вблизи призмы при
угле атаки 20 градусов.
Больше всего нас интересуют коэффициент подъёмной силы C y и
коэффициент лобового сопротивления C x и их зависимость от угла атаки.
Cx =
Cy =
[5]
.
ЗдесьY – подъёмная сила, X–сила лобового сопротивления, ρ–плотность
воздуха, v–скорость, S – характерная площадь.
Для примера, график зависимости коэффициента лобового сопротивления от
шага программы для угла атаки 0 градусов:
15
А вот так будет выглядеть зависимость коэффициента подъёмной силы от
шага интерации для угла атаки -10 градусов:
Как Вы видите, значение коэффициентов перестаёт колебаться после 17000го шага программы, поэтому после 20000-го шага я прекращал вычисления.
Я получил следующие значения этих коэффициентов для различных углов
атаки:
-20 -15 -10 -5
-2
0
2
5
10
15
20
Cx
1,27 1,20 1,21 1,27 1,30 1,29 1,30 1,27 1,22 1,20 1,26
5
9
4
6
3
9
7
5
9
Cy
0,37 0,48 0,39 0,21 0,10 0,04
6
4
2
5
3 0,02 0,15 0,36 0,43 0,34
5
4
5
3
8
Cycosα+Cxs
0,15 0,17 0,1 0,05 0,04 0,02
0,10
inα
0,08 2
6
9
3
0,04 0,14 0,10 7
1
1
6
6
Cn=Cycosα + Cxsinα мы считаем для получения коэффициента нормальной
силы (аэродинамические коэффициенты программа высчитывает для
скоростной системы координат набегающего потока, поэтому мы определяем
проекции на оси связанной системы координат).
16
В итоге мы получаем следующую зависимость коэффициента нормальной
силы Cn от угла атаки:
Сравнение результатов, полученных с помощью программы SU2 и
практических результатов.
Естественно, интересно сравнить значения Cn, полученные программой SU2
со значениями, получаемыми экспериментально. Для примера приведу
17
сравнение с результатами, полученными В.Д. Люсиным и А.Н. Рябининым
[6]
в аэродинамической трубе АТ-12 с открытой рабочей частью.
Параметры моих расчётов в программе SU2: скорость 19,5 м/с, плоское
обтекание.
Рабочие параметры экспериментальных вычислений: скорость потока 27,8
м/с, удлинение сечения призмы λ = 20.
Результаты:
α,
град
Cn
-28,6
-22,9
-17,2
-5,7
0
5,7
14,9
22,9
28,6
-0,1
-0,05
0,35
0,15
0,05
-0,05
-0,37
-0,15
0,12
Как мы видим, значения коэффициента нормальной силы приблизительно
совпадают с данными, полученными экспериментальным путём.
Изучение возможности галопирования для призмы данной формы.
Одна из основных целей моей работы – определение возможности
галопирования призмы. Как уже упоминалось, галопирование – процесс
18
колебания различных тел плохообтекаемой формы в потоке газа или
жидкости перпендикулярно направлению потока. Если предположить, что
призма будет двигаться только относительно оси y, то после возникновения
колебаний Vr будет скоростью потока относительно призмы и будет
складываться из скорости призмы относительно оси y и скорости потока.
Колебания будут увеличиваться до определённого предела, поэтому угол
атаки тоже будет иметь гармоническую зависимость, хотя тело совершает
только поступательные движения и не поворачивается. При отсутствии же
колебаний угол атаки будет постоянно равен нулю.
Вследствие этого, аэродинамические силы могут обеспечивать приток
энергии в систему, отток же в свою очередь будет получаться за счёт вязкого
демпфирования (гашения колебаний). Таким образом, когда получается
баланс между притоком и оттоком энергии, колебания получаются
установившимися.
Как мы видим из графика, около нуля C n, и, соответственно, сама нормальная
аэродинамическая сила будут иметь противоположный знак с углом атаки.
Наличие диапазона углов, при которых будут наблюдаться такие
зависимости и является необходимым условие возможности появления
галопирования [7], так как для тела, совершающего колебания это означает,
что при некотором отношении собственной скорости тела к скорости потока
аэродинамическая сила будет действовать в том же направлении, что и
собственная скорость тела, а это в свою очередь влечёт к увеличению
амплитуды.
Следовательно, в нашем случае будет наблюдаться галопирование.
19
Заключение.
Итак, в рамках своей дипломной работы я изучал аэродинамические силы
для прямоугольной призмы. Модель призмы и сетка были построены с
помощью программы Gmsh. Была построена сетка вблизи призмы с
размерами 1*1*0,025 .
Затем, я моделировал процесс обдувания призмы ветровым потоком со
скоростью 19,5 м/с с помощью программы SU2. Уравнение Навье-Стокса
описывает рассматриваемую мной модель, которая с помощью осреднения
по Рейнольдсу приводится к изменённой системе уравнений. Но данная
система уравнений имеет большее количество неизвестных, чем уравнений,
поэтому мы дополняем её дополнительным соотношением: модель
турбулентности Спаларта-Аллмараса.
В программе SU2 я произвёл моделирование и расчёт коэффициентов для
углов атаки от – 20 до 20 градусов. В результате мы получаем зависимость
коэффициента нормальной силы C n от угла атаки, и исходя из этой
зависимости заключаем, что для данной призмы галопирование возможно.
Результаты моей работы достаточно точно совпадают с результатами других
авторов, полученными в аэродинамической трубе.
В итоге, я увидел огромные возможности в программе SU2, кроме того, она
обладает важными преимуществами, поэтому я предполагаю, что в будущем
большинство задач аэродинамики и гидродинамики будут решаться с
помощью таких программ.
Список используемой литературы.
[1] Рейнольдс О. Динамическая теория движения несжимаемой вязкой
жидкости и определение критерия. Проблемы турбулентности. – М.; Л.:
ОНТИ, 1936. – стр. 135- 227.
[2] SpalartP. R., AllmarasS. R. A one-equation turbulence model for
aerodynamic flows // AIAA Paper 1992-0439.
[3] Гарбарук И. В., Стрелец М. Х., Шур М. Л. Моделирование
турбулентности в расчётах сложных течений. СПб.:Издательство
политехнического университета, 2012, страницы 42-49.
20
[4] ChristopheGeuzaine. «GMSH»
[5] ВалландерС. В..Лекции по гидроаэромеханике. Л.: Издательство
ленинградского университета, 1978.
[6] В. Д. Люсин, А. Н. Рябинин. Исследование влияния удлинения призмы на
её аэродинамические характеристики и амплитуду колебаний при
галопировании. // Вестник СПБГУ, Серия 1, 2011,вып. 2, с. 1-7.
[7] Дж М. Т. Томпсон. Неустойчивости и катастрофы в науке и технике.
Москва: МИР, 1985, стр 38-39, 181-189.
21
Отзывы:
Авторизуйтесь, чтобы оставить отзыв