close

Вход

Забыли?

вход по аккаунту

?

Численная реализация метода решения уравнений Навье

код для вставкиСкачать
Федеральное агентство по образованию
Государственное образовательное учреждение высшего
профессионального образования
Санкт-Петербургский государственный морской технический университет
Факультет кораблестроения и океанотехники
Кафедра прикладной математики и математического моделирования
Разработка универсального программного
кода для решения трехмерных нестационарных
задач механики сплошных сред
Автор
Медведев А. А.
Научный
Рыжов В.А.
работы
руководитель
д.т.н.,
Санкт-Петербург
2010
профессор,
Реферат
Настоящая
дипломная
работа
посвящена
математическому
моделированию нестационарных течений вязкой несжимаемой жидкости в
двумерных областях сложной геометрии с подвижными границами.
Предложенный метод основывается на решения
Рейнольдсу уравнений Навье-Стокса, замкнутых при
турбулентности Спаларта-Аллмараса.
осредненных по
помощи модели
Расчетный алгоритм построен с использованием метода искусственной
сжимаемости, что позволяет избежать возникновения неустойчивости решения
при наложении условия несжимаемости. Получение монотонного решения при
сохранении точности обеспечивается путем применения противопоточной схемы
высокого порядка для расчета конвективных слагаемых. Для пространственной
дискретизации определяющих уравнений применяется метод конечных объемов
на неструктурированных треугольных сетках. Нахождение конвективного и
диффузионного потоков через границу контрольной ячейки осуществляется при
помощи полиномиальной аппроксимации второго порядка точности на
неструктурированных сетках.
Для дискретизации по времени в общем случае используется неявная
схема второго порядка. При расчете нестационарного течения применяется метод
двойных шагов по времени. Результирующая система линейных уравнений
решается методом бисопряженных градиентов.
Проанализировано влияние порядка аппроксимации слагаемых уравнений
Навье-Стокса на точность вычислений и устойчивость решения. Верификация
расчетного метода показала удовлетворительное соответствие полученных
результатов расчета для различных режимов течения в двумерных областях с
известными и расчетными данными.
Оглавление
ОСНОВНЫЕ ОБОЗНАЧЕНИЯ.......................................................................................... 5
ВВЕДЕНИЕ ........................................................................................................................ 6
КРАТКИЙ ОБЗОР МЕТОДОВ ВЫЧИСЛИТЕЛЬНОЙ ГИДРОДИНАМИКИ .................... 8
1.1
Математическая модель движения жидкости ................................................................. 8
1.1.1 Уравнения движения сплошной среды ................................................................................ 8
1.1.2 Уравнения Навье-Стокса ....................................................................................................10
1.1.3 Динамическое подобие ......................................................................................................10
1.1.4 Приближения уравнений Навье-Стокса. Моделирование турбулентности ...................11
1.1.5 Выводы ................................................................................................................................19
1.2 Численные методы решения задач вычислительной гидродинамики ..........................20
1.2.1 Дискретизация определяющих уравнений .......................................................................20
1.2.2 Разбиение расчетной области на элементы .....................................................................22
1.2.3 Уравнения Навье-Стокса для несжимаемой жидкости ....................................................25
1.2.4 Методы решения СЛАУ .......................................................................................................28
1.2.5 Выводы .................................................................................................................................30
ЧИСЛЕННАЯ РЕАЛИЗАЦИЯ МЕТОДА РЕШЕНИЯ УРАВНЕНИЙ НАВЬЕ-СТОКСА. 31
2.1 Общая постановка задачи и основные допущения ............................................................31
2.2 Особенности метода расчета ..................................................................................................31
2.3 Дискретизация по пространству. Метод конечных объемов ...........................................31
2.3.1 Аппроксимация невязкого потока .......................................................................................31
2.3.2 Полиномиальная аппроксимация.......................................................................................31
2.3.3 Аппроксимация вязкого потока ...........................................................................................31
2.3.4 Граничные условия ..............................................................................................................31
2.4 Неявная дискретизация по времени ......................................................................................31
2.4.1 Двойной шаг по времени .....................................................................................................31
2.4.2 Одинарный шаг по времени ...............................................................................................31
2.4.3 Построение матрицы СЛАУ ................................................................................................31
2.5 Расчет турбулентной вязкости ...............................................................................................31
2.5.1 Модель Спаларта-Аллмараса ............................................................................................31
2.5.2 Дискретизация по пространству и времени ......................................................................31
2.6 Расчет гидродинамических характеристик ..........................................................................31
2.7 Программная реализация. Расчетные алгоритмы ..............................................................31
2.7.1 Алгоритм расчета нестационарного течения ....................................................................31
2.7.2 Решение СЛАУ .....................................................................................................................31
2.7.3 Краткая характеристика вычислительных средств ..........................................................31
3 ВЕРИФИКАЦИЯ РАСЧЕТНОГО МЕТОДА ................................................................. 31
3.1 Оценка сходимости метода .....................................................................................................32
3.2 Результаты тестовых расчетов ..............................................................................................32
3.2.1 Задача турбулентного обтекания плоской пластины .......................................................32
3.2.2 Задача обтекания двумерного кругового цилиндра в канале .........................................32
3.2.3 Задача обтекания двумерного канала с обратным уступом ...........................................32
3.2.4 Задача обтекания двумерного квадратного цилиндра в неограниченной жидкости ....32
3.2.5 Задача обтекания двумерного кругового цилиндра, совершающего заданные
колебания поперек потока ....................................................................................................................32
4 ИССЛЕДОВАНИЕ ПРОПУЛЬСИВНЫХ ХАРАКТЕРИСТИК МАШУЩЕГО КРЫЛА. 32
4.1 Постановка задачи обтекания машущего крыла на упругих связях ...............................32
4.2 Пропульсивные характеристики машущего крыла ............................................................32
4.3 Анализ результатов расчета пропульсивных характеристик движителя типа
машущее крыло .......................................................................................................................................32
4.3.1 Угловые заданные колебания профиля ............................................................................32
4.3.2 Вертикально-угловые заданные колебания профиля ......................................................32
4.3.3 Колебания жесткого профиля на упругих связях ..............................................................32
ЗАКЛЮЧЕНИЕ ................................................................................................................. 32
ЛИТЕРАТУРА .................................................................................................................. 32
ПРИЛОЖЕНИЯ ................................................................................................................ 34
ИЛЛЮСТРАЦИИ .............................................................................................................. 34
Основные обозначения

v
– вектор скорости
u , v , w , vi
– компоненты вектора скорости
x, y , z , xi –
декартовы координаты точки
t
– физическое время
p
– гидродинамическое давление
– плотность жидкости

 ij – тензор напряжений
 ij – тензор вязких напряжений

f – вектор напряженности массовых сил

n – единичный вектор нормали к поверхности
– число Рейнольдса
Re
– число Струхаля
St
 ,  t – молекулярная и турбулентная динамические вязкости
 ,  t – молекулярная и турбулентная кинематические вязкости
– удельная кинетическая энергия турбулентности
k
– скорость диссипации турбулентности

 – удельная диссипация турбулентности
 – параметр искусственной сжимаемости
 – искусственное время
S
– площадь контрольной ячейки
– вектор физических переменных в ячейке
q
 q
R
– приращение вектора физических переменных в ячейке
, F – вязкий и невязкий потоки через границу ячейки
 – ограничивающий множитель
C f – коэффициент местного трения
C d , C t – коэффициенты сопротивления и силы тяги
C l – коэффициент подъемной силы
C m – коэффициент момента
C p – коэффициент затраченной мощности
Введение
Среди разнообразия современных движительных комплексов особое место
занимают пропульсивные системы с машущим крылом. Это объясняется в первую
очередь тем, что принципы, лежащие в основе подобных устройств, являются
практическим продолжением идей, заложенных природой в живых объектах,
прошедших огромный путь эволюционного развития.
Существующие на сегодняшний день достижения в области разработки
малых подводных
необитаемых аппаратов (AUV – Autonomous Underwater
Vehicle) и микролетательных аппаратов (MAV — Micro Aerial Vehicle) с
пропульсивно-несущей системой типа машущее крыло [44] говорят о
перспективах их использования на практике и возможностях дальнейшего
совершенствования в целях получения уникальных гидродинамических,
гидроакустических и экологических качеств, присущих гидробионтам.
В отличие от традиционных типов движителей, данный тип является
экологически чистым, так как его воздействие на окружающую среду (акустическое
излучение, силовое воздействие и др.) является минимальным. С использованием
рассматриваемого устройства, функционирующего в режиме генератора, может
решается проблема освоения возобновляемых природных источников энергии
(морских волн и воздушных потоков).
Современные достижения в области новых конструкционных материалов,
микроэлектромеханических систем, источников и преобразователей энергии,
открывают новые горизонты в практическом развитии данного направления
исследований и требуют разработки математических моделей, адекватно
отражающих реальные физические процессы, и, следовательно, позволяющих
прогнозировать оптимальные режимы функционирования пропульсивных
устройств.
Безусловно, важный вклад в процесс создания высокоэффективного
машущего движителя дают теоретические исследования, определяющие на
начальных этапах разработки применимость тех или иных технических решений
для практических целей, выявляющие оптимальные значения проектных
параметров, требуемых для достижения необходимых пропульсивных,
маневренных и энергетических характеристик для различных режимов движения и
условий эксплуатации.
Подобные исследования проводятся в настоящее время различными
исследовательскими центрами, среди которых необходимо упомянуть:
Massachusetts Institute of Technology, Naval Postgraduated School, СанктПетербургский государственный морской технический университет, Tokyo Institute
of Technology, Necton Research LLC, UC Berkeley, GTRI, George Washington
University и пр. [44]. Можно констатировать, что к настоящему времени этими и
другими
школами
достигнуты
определенные
результаты
в
области
математического моделирования пропульсивных характеристик машущего крыла.
Очевидно, что наиболее достоверную информацию можно получить в
рамках общей модели нестационарного обтекания потоком вязкой жидкости
упругого телесного крыла конечного размаха произвольной формы в плане,
совершающего колебания по сложному закону с большими амплитудами с учетом
взаимодействия с телом (корпусом). Существующие на сегодня решения
получены с допущениями, упрощающими общую модель. Наибольшее количество
исследований выполненных ранее относится к модели обтекания крыла потоком
невязкой несжимаемой жидкости. В рамках подобной модели в линейной и
нелинейной постановках проведен обширный анализ влияния проектных
параметров на аэрогидродинамические характеристики машущего крыла.
Многочисленные результаты получены как для двумерного, так и для трехмерного
случаев.
Что касается модели обтекания машущего крыла потоком вязкой жидкости,
то количество результатов здесь существенно меньшее и относятся они в
основном к моделированию обтекания жесткого машущего профиля. Работ по
моделированию трехмерного вязкого обтекания машущего крыла - единицы. При
этом исследование гидродинамических характеристик проведено для небольшого
количества проектных параметров, численный анализ не носит систематического
характера.
Настоящая работа имеет своей целью продолжение проводимых в СПбГМТУ
исследований пропульсивных характеристик машущего движителя с разработкой
универсального программного кода для решения трехмерных нестационарных
задач механики сплошных сред.
Для достижения указанной цели в работе решены следующие задачи:

построена математическая модель нестационарного обтекания тела потоком
вязкой несжимаемой жидкости,

разработан универсальный программный код для решения трехмерных
нестационарных задач механики сплошных сред (”SmartFlow”),

проведена верификация расчетного метода для различных модельных задач,
выполнено сравнение полученных расчетных результатов с известными
экспериментальными и расчетными данными,

выполнено
численное
моделирование
нестационарного
обтекания
колеблющегося крыльевого профиля NACA0012 для случаев угловых и
совместных вертикально угловых колебаний, а также вынужденных
гидроупругих колебаний профиля, закрепленного на упругих связях,

проведен анализ полученных расчетных зависимостей пропульсивных
характеристик машущего движителя для различных проектных параметров.
Результаты
конференциях:
данной
работы
представлены
на
следующих
научных
Краткий обзор методов вычислительной гидродинамики
1.1 Математическая модель движения жидкости
Необходимым условием успешного решения задачи механики является
правильный выбор математической модели, адекватно отражающей исследуемые
динамические процессы.
При рассмотрении безвихревого обтекания тел потоком идеальной
жидкости возникает парадокс Эйлера-Даламбера - отсутствие каких-либо сил,
действующих на тело. Причиной его появления является пренебрежение
процессами вихреобразования в данной модели течения. Вихревые структуры в
потоке жидкости возникают под действием вязкости, поэтому для их исследования
необходимо использовать модели вязкой среды. Известны различные режимы
течения потоков вязкой жидкости - ламинарный, турбулентный и смешанный
ламинарно-турбулентный. Появляющиеся в таких потоках вихри имеют некоторые
отличия в своей динамике, обусловленные различием механизмов, действующих
в потоке.
1.1.1 Уравнения движения сплошной среды
Прежде, чем перейти к непосредственному численному моделированию
течений
вязкой
жидкости,
рассмотрим
основополагающие
уравнения,
описывающие это течение.
В гидромеханике применяют два основных метода изучения движения
жидкости - подходы Лагранжа и Эйлера. Любой жидкий объем можно представить
состоящим из большого числа жидких частиц. В соответствие с этим, к
исследованию движения жидкой частицы возможен такой же подход, как и в
теоретической механике. То есть, для каждой точки, однозначно определяемой
начальными координатами, в любой момент времени известны ее скорость и
ускорение. Такой метод хорош при рассмотрении задач диффузии, при описании
одномерных потоков. В более сложных случаях он приводит к громоздким
вычислениям. Метод предложен Лагранжем и носит его имя.
Метод характеристики движения, при котором в каждой точке задаются
функции зависимости характеристик течения от времени, но частицы теряют свою
индивидуальность, называется методом Эйлера. Например, для случая
нестационарного течения жидкости, поле скоростей задается в виде
⃗ = ⃗ (⃗, )
Большинство
приборов
измеряют
характеристики
жидкости
в
фиксированном месте (датчик прибора неподвижен), то есть определяют
Эйлерову характеристику среды.
Все основные уравнения движения сплошной среды представляют собой
фундаментальные законы сохранения [10, 14]. Для вывода уравнений движения
жидкости обычно рассматривается малый контрольный объем и требуется, чтобы
для жидкости, протекающей через этот объем, выполнялись законы сохранения
массы и энергии и количества движения.
Согласно закону сохранения вещества, для произвольного неподвижного
объема Ω скорость изменения массы внутри него равна потоку массы через
поверхность , ограничивающую этот объем. Уравнение закона сохранения массы
(уравнение неразрывности) для некоторого объема в инерциальной системе
координат, записанное в интегральной форме, имеет вид
∂
∫
∂ Ω
 Ω + ∫ (⃗, ⃗⃗)  = 0
(1)
где  -- плотность жидкости. Эквивалентное дифференциальное уравнение
в частных производных
∂
+ ∇ ⋅ (⃗) = 0
∂
(2)
То-же в координатной форме записи
∂
∂
+
∂
∂
=0
(3)
Для несжимаемой жидкости, учитывая, что плотность есть величина
постоянная ( =  )
∂
∂
=0
(4)
В соответствие со вторым законом Ньютона, скорость изменения
количества движения равна сумме действующих сил. Уравнение закона
сохранения количества движения (уравнение динамики в напряжениях) для
некоторого объема, имеет вид
∂
∫
∂ Ω
⃗ Ω + ∫ ⃗(⃗, ⃗⃗) − ⃗⃗  = ∫Ω ⃗ Ω
(5)
где  - тензор напряжений, ⃗ - вектор напряженности массовых сил.
Переходя к уравнениям в частных производных, получаем консервативную
(дивергентную) форму уравнения сохранения количества движения
∂( )
∂
+
∂(  )
∂
=
∂
∂
+ 
(6)
Эквивалентная неконсервативная (конвективная) форма
∂( )
∂
+ 
∂( )
∂
=  +
∂
∂
В механике вязкой несжимаемой жидкости предполагается, что перенос
тепла происходит мгновенно в силу большой скорости передачи тепла в
несжимаемой жидкости, поэтому изменения температуры пренебрежимо малы. В
силу этого нет необходимости определять изменение термодинамического
состояния системы по балансу внутренней и механической энергий, и уравнение
сохранения энергии не используется.
В результате используются только уравнения, выражающие законы
сохранения массы (4) и количества движения (6). Эта система является
незамкнутой в силу неопределенности тензора напряжений. Для ее замыкания
вводятся гипотезы о связи компонентов тензора напряжений со скоростями
потока, то есть используются реологические соотношения.
1.1.2 Уравнения Навье-Стокса
Обобщенный закон Ньютона для вязкой жидкости [10] устанавливает
линейную связь между тензором напряжений и тензором относительных
скоростей деформации
∂
∂
2
∂
 = −  +  [(∂ + ∂ ) − 3  ∂ ] , , ,  = 1,2,3



(7)
где  -- гидродинамическое давление,  -- молекулярная динамическая
вязкость,  -- дельта-функция Кронекера.
В случае несжимаемой жидкости, тензор напряжений выглядит следующим
образом
∂
∂
 (∂ + ∂ ) ,
 =

при  ≠ 

− + 2
∂
∂
,
при  = 
(8)
{
Тензор напряжений часто разделяют на две части
∂
∂
 = −  +  (∂ + ∂ ) = −  + 


(9)
где  -- тензор вязких напряжений. Подстановка (8) в уравнения для
напряжений (6) дает известные уравнения Навье-Стокса
∇ ⋅ ⃗
+ ⃗ ⋅ ∇⃗
∂
⃗⃗
∂
=
=
0
1
⃗
 −  ∇ + Δ⃗
(10)
где,  = / -- молекулярная кинематическая вязкость.
Для получения конкретных решений, при интегрировании системы (10)
должны быть использованы граничные, а в случае нестационарного движения --граничные и начальные условия. На твердых границах задаются условия
''непротекания'' и ''прилипания'' ⃗| = 0. Начальные условия ставятся в задачах
нестационарного движения и представляют собой заданные в некоторый
начальный момент времени поля скоростей и давлений.
1.1.3 Динамическое подобие
Для того, чтобы наиболее оптимальным образом (с точки зрения
проведения минимального числа расчетов или экспериментальных наблюдений)
получить картину течений у тел подобной конфигурации, желательно
сгруппировать все параметры (такие как длина тела, скорость набегающего
потока и т. п.) в ряд безразмерных параметров [14]. Два потока динамически
подобны, если безразмерные числа, определяющие течения, равны.
Безразмерные переменные вводятся следующим образом

∗ =   ,  ∗ =
0
0 
0

, ∗ =   , ∗ =
0
−0
 02
Уравнения движения жидкости (10) в безразмерной координатной форме
принимают вид

∂∗
+ ∗
∂ ∗
∂∗
∂∗
∂∗
1 ∂2 ∗
=

∂∗2
∂∗

=
∂∗
1
− ∂ ∗ +  2 
0
В данном выражении использованы следующие безразмерные параметры
 =
0 0

,
0
 = 
0 0
,

 = ( 0)1/2
0
Здесь  -- число Рейнольдса,  -- число Струхаля и  -- число Фруда
соответственно. Число Рейнольдса характеризует отношение инерционных сил к
вязким, и является критерием, определяющим этапы перехода от ламинарных
течений к турбулентным.
1.1.4 Приближения уравнений Навье-Стокса. Моделирование турбулентности
Основная трудность расчета потоков вязкой несжимаемой жидкости
частично связана с широким диапазоном изменения масштаба турбулентности.
Прямой расчет полных уравнений Навье-Стокса для трехмерного турбулентного
потока требует значительных вычислительных ресурсов и не под силу даже
существующим суперкомпьютерам. В этой связи, важную роль играет
турбулентная модель, позволяющая учесть влияние турбулентности в расчетах.
Вторая трудность в расчете вязкого потока связана с необходимостью
использования чрезвычайно мелких сеток при расчете течения в турбулентном
пограничном слое. Поскольку вычислительная устойчивость существующей схемы
решения непосредственно связана с размером минимальной ячейки, то, если не
принимать достаточно малый шаг по времени, возникают проблемы устойчивости
расчета. В результате, и увеличение разрешения и уменьшение шага по времени
при вычислениях влекут за собой резкое увеличение требуемых вычислительных
ресурсов.
Поэтому для решения задач гидродинамики применяются различные
подходы, основной целью которых является уменьшение ''вычислительной
стоимости'' методов решения уравнений Навье-Стокса при минимально
возможной потере точности.
Различают два основных подхода моделирования вязких течений
• Прямое численное моделирование (DNS) --- решение полных уравнений
Навье-Стокса.
• Моделирование с использованием осредненных уравнений НавьеСтокса, а именно: по времени (RANS), по пространству (LES), гибридные
модификации (DES).
Прямое численное моделирование: DNS
Среди известных методов численного моделирования трехмерных
турбулентных течений необходимо выделить прямое численное моделирование
турбулентности (DNS --- Direct Numerical Simulation of turbulent flows).
Метод DNS представляет собой прямое численное решение полной
нестационарной системы уравнений Навье-Стокса, при таком подходе
разрешаются все масштабы движения [22]. В результате возникает
необходимость
строить
чрезвычайно
мелкую
сетку
для
больших
пространственных областей. Известна следующая оценка числа узлов при
прямом моделировании турбулентности
 = ( 9/4 )
Для реальных чисел Рейнольдса порядка 106−8 число расчетных узлов
должно составлять  = 1013−15 . То есть для использования DNS требуются
достаточно мощные вычислительные ресурсы, и на сегодняшний день
возможности применения метода ограничиваются лишь случаями достаточно
простых течений и весьма невысоких чисел Рейнольдса.
Метод моделирования крупных вихрей: LES
В методе моделирования крупных вихрей (LES --- Large Eddy Simulation)
осуществляется решение отфильтрованных по пространству уравнений НавьеСтокса и разрешеатся движение только крупных вихрей [31].
Метод основан на двух предположениях. Первое состоит в возможности
разделения поля скорости на движение крупных и мелких вихрей, причем
движение крупных вихрей может быть рассчитано отдельно, что связано с
достаточной
изотропностью
и
универсальностью
мелких
масштабов
турбулентного
движения.
Второе
предположение
--в
возможности
аппроксимации нелинейных взаимодействий между крупными и мелкими вихрями
только о крупным вихрям с использованием моделей подсеточного масштаба SGS
(SubGrid Scale models).
Для отделения крупноых масштабов от мелких, применяется операция
фильтрации, определяемая следующим образом
(⃗) = ∫ (⃗) (⃗, ⃗, Δ) 
(11)
где  -- фильтрационная функция, Δ -- ширина фильтра, определяющая
наименьший масштаб турбулентности, допустимый фильтром. Наиболее
популярные и часто используемые фильтрационные функции --- Гаусса,
идеальный и ''top-hat'' фильтры.
Фильтр дает формальное определение процесса осреднения и отделяет
способные к разрешению масштабы от подсеточных. Фильтрация используется,
чтобы вывести уравнения для разрешимых масштабов. Для течения несжимаемой
жидкости отфильтрованные уравнения Навье-Стокса принимают следующую
форму
∂
=
∂
∂
∂
+
∂(  )
∂
=
0
1 ∂
−
∂
−
∗
∂
∂
+
∂2 
∂2
(12)
Здесь воздействие мелкомасштабных структур на движение жидкости
представляется через тензор напряжений подсеточного масштаба
∗

=   −  
(13)
Среди применяемых подсеточных моделей можно выделить модель
Смагоринского, двухточечные замыкания, динамические модели, модели одого
уравнения [40].
Популярность метода моделирования крупных вихрей для проведения
расчетов сложных турбулентных течений с достаточно высокими числами
Рейнольдса объясняется тем, что он требует меньших вычислительных затрат по
сравнению с DNS. Общее соотношение количества узлов сетки для этих методов
определяется зависимостью
 ≈ 0,4() 1/4 
Необходимо отметить, что на сегодняшний день опробовано значительное
количество подсеточных моделей, фильтров, граничных условий и расчетных
схем. Несмотря на это, пока не ясны ни оптимальный вариант подсеточной
модели, ни обоснование выбора такого варианта. Тем не менее, LES является
перспективным направлением в развитии методов расчета турбулентных течений
и представляется весомой альтернативой DNS и RANS.
Уравнения осредненного движения: RANS
Как было отмечено, решение полных уравнений Навье-Стокса в
трехмерном пространстве при больших (турбулентных) числах Рейнольдса
остается на сегодняшний день довольно сложной задачей. Поэтому для описания
трехмерных течений часто используют осредненные по времени уравнения
Навье-Стокса. В турбулентном течении локальные давление и составляющие
вектора скорости изменяются во времени случайным образом. Основная идея
осреднения состоит в том, чтобы разделить в потоке стационарные и случайные
части.
Система уравнений Навье-Стокса для описания движения вязкой
несжимаемой жидкости при отсутствии массовых сил, использующая
консервативную форму записи уравнения изменения количества движения, может
быть представлена в скалярно-тензорной форме следующим образом
∂
∂
∂
∂
∂(  )
+
∂
=
=
0
−
1 ∂
 ∂
+
1 ∂
(14)
 ∂
С учетом уравнения неразрывности, компоненты тензора напряжений
записываются так
1 ∂
 ∂
= 
∂2 
∂2
(15)
Согласно
подходу
Рейнольдса,
любые
мгновенные
значения
гидродинамических параметров потока представляются в виде суммы
осредненной по времени величины и ее пульсационной составляющей [4].
Фактически это означает, что гидродинамическая величина является случайной,
осреднение которой во времени дает математическое ожидание, а пульсационная
составляющая - дисперсию случайной величины. Обозначая осредненную во
времени величину ( ) , а пульсационную ( )′ , можно записать для давления,
составляющих скорости, и тензора напряжений следующие выражения
 =  + ′,
 =  + ′ ,
 =  +  ′
Следует отметить, что среднее значение, несмотря на интегрирование по
времени, может изменяться во времени. Это означает, что период
интегрирования  должен быть малым по сравнению с характерным временем
нестационарного изменения величины. Кроме того, период осреднения
выбирается так, чтобы оно приводило к величине, не изменяющейся при
повторном осреднении.

1
 () =  ∫0  () 
(16)
Применяя операцию осреднения по времени (16) к уравнениям системы
(46), с учетом уравнения неразрывности, получим
∂
=
∂
∂( )
∂
+
∂(  )
∂
0
∂
∂
= − ∂ + ∂ ( −  ′ ′ )

(17)

где  ′ ′ -- составляющие тензора напряжений Рейнольдса, или
рейнольдсовых напряжений. Они являются шестью дополнительными
неизвестными к гидродинамическим параметрам осредненного движения ( , ).
Таким образом, система уравнений и (17) является незамкнутой.
Вопрос замыкания полученной системы решается различными способами.
Простейший путь --- использование эмпирической информации о характеристиках
турбулентности, наиболее сложный --- заключается в выводе уравнений
относительно рейнольдсовых напряжений, где широкое применение получили
модели турбулентной вязкости
∂
∂
− ′ ′ =  (∂ + ∂)


(18)
где  -- турбулентная динамическая вязкость.
Используя зависимость (18) и опуская знак осреднения по времени,
уравнения Навье-Стокса в форме Рейнольдса (RANS --- Reynolds Averaged
Navier-Stokes equations) приводятся к виду
∇ ⋅ ⃗
+ ⃗ ⋅ ∇⃗
∂
⃗⃗
∂
=
=
1
−
0
∇ + ∇[( +  )∇⃗ ]
(19)
где  =  / -- турбулентная кинематическая вязкость.
Само по себе уравнение (18) не вводит модели турбулентности, а только
характеризует структуру такой модели. При этом основной задачей является
определение коэффициента турбулентной вязкости  . В отличие от
коэффициента молекулярной кинематической вязкости  , коэффициент 
определяется состоянием турбулентности и не связан со свойствами жидкости.
Он может сильно изменяться от точки к точке пространства и в зависимости от
типа течения. Так, например, в зонах циркуляционного течения  может на
несколько порядков превышать .
Модели
турбулентности
обычно
классифицируются
по
числу
дифференциальных уравнений, вводимых в дополнение к уравнениям НавьеСтокса. Различают модели ''0-уравнений'' (алгебраические), ''1-уравнений''
(модель Спаларта-Аллмараса), ''2-уравнений'' ( −  и  −  модели). Рассмотрим
наиболее часто используемые турбулентные модели.
Алгебраические модели турбулентности
Алгебраические модели принадлежат к простейшим типам моделей
турбулентности, в которых связь между турбулентной вязкостью и параметрами
осредненного потока задается алгебраическими соотношениями. Отсюда следуют
основные достоинства моделей такого типа вычислительная эффективность,
простота калибровки и модификаций с учетом специфики рассматриваемых
течений. Однако очевидна и узкая специализация этих моделей, поскольку они
опираются на априорную (эмпирическую) информацию о структуре конкретного
рассматриваемого течения.
Модель пути смешения Прандтля (Prandtl, 1925)
Модель для описания распределения  впервые была предложена
Прандтлем в 1925 году и известна как модель пути смешения. Доказано, что она
довольно хорошо воспроизводит тонкие вязкие слои. Рассматривая осредненные
сдвиговые течения без градиента давления, Прандтль постулировал, что
характерный масштаб пульсаций скорости равен градиенту осредненной
скорости, умноженному на характерный масштаб длины  , который он назвал
путем смешения
2
 = 
∂|| ∂
(20)
Длина пути смешения определяется эмпирически. При рассмотрении
течения в пограничном слое полагается
 = 
(21)
где  ≈ 0.39 -- универсальный коэффициент пропорциональности, не
зависящий от числа Рейнольдса. Таким образом, путь перемешивания
пропорционален расстоянию от стенки .
Успех
модели,
предложенной
Прандтлем,
определяется
тем
обстоятельством, что для многих типов сдвиговых течений  может быть
выражен относительно несложными формулами.
Модель Болдуина-Ломакса (Baldwin-Lomax, 1978) [16]
Модель была сформулирована для расчета потока в тех случаях, когда
параметры пограничного слоя (толщина и скорость на границе) трудно
определить. Такая ситуация часто возникает при численном моделировании
отрывных течений, в особенности, течений со скачками уплотнения.
Турбулентная вязкость определена следующим соотношением
 ,  ≤ 
 = { ,  > 
(22)
где  и  являются значениями  во внутреннем и внешнем слоях, а 
есть наименьшее значение , при котором  =  . Вязкости во внутреннем и
внешнем слоях принимаются так
внутренний слой
2
 = 
||
 = (1 −  −
(23)
+ /+

)
внешний слой
 =    
(24)
2
 = min (  ;   
/ )
1
 =  [max( ||)]
где || -- величина вектора завихренности,  -- величина , при котором
 || достигает максимальной величины,  -- максимальная величина скорости
в пограничном слое. Функция  является функцией перемежаемости
Клебанова
Коэффициенты замыкания
 = 0.40,
{ = 1.6,
 = 0.0168,
 = 0.3,
+
 = 26
 = 1
Алгебраические модели, безусловно, являются наиболее простыми из всех
турбулентных моделей. Они концептуально очень просты и редко вызывают
неожиданные вычислительные трудности. Однако следует всегда помнить о
проблеме неполноты информации, получаемой с их помощью. Эти модели
хорошо работают только при анализе тех потоков, на которые они были
предварительно настроены.
Ограниченность моделей такого типа заключается в их природе --- в
локальном равновесии моделируемой турбулентности. Это означает, что в каждой
точке пространства наблюдается баланс генерации и диссипации турбулентной
энергии, на который не влияют ни перенос из соседних точек, ни предыдущее
развитие процесса. Трудности для сложных типов течений представляет задание
пути смешения. Однако для простых ситуаций, в частности при описании
сдвиговых слоев, модель вполне пригодна.
Модели одного уравнения
Такие модели дают описание турбулентности с помощью одной переменной
величины, для которой строится дифференциальное уравнение переноса. Другие
турбулентные характеристики связываются с ней при помощи алгебраических или
иных соотношений. К данному классу относится модель Спаларта-Аллмараса.
Модель Спаларта-Аллмараса (Spalart-Allmaras, 1992) [37]
Данная модель относится к классу низкорейнольдсовых. Первоначально
она была развита для получения разумных расчетных оценок для двумерных
смешанных течений, следов и пограничного слоя на плоской пластине. Испытания
показали достоинство этой модели при расчете потоков с неблагоприятными
градиентами давления по сравнению с  −  и  −  моделями [4]. Определяющие
уравнения модели таковы
Кинематическая турбулентная вязкость
 = ̃1
(25)
Уравнение турбулентной вязкости
̃
∂
̃
∂
̃
∂
1 ∂
+  ∂ = 1 ̃̃ − 1  ̃() 2 +  ∂ [( + ̃) ∂ ] +
∂



2 ∂
̃ ∂
̃
 ∂ ∂
(26)
Коэффициенты замыкания и вспомогательные соотношения
1 = 0.1355, 2 = 0.622, 1 = 7.1,  = 2/3
{1 = 3.239, 2 = 0.3,
3 = 2,
 = 0.41
̃

 =  , 1 =  3  3 + 31 , 2 = 1 − 1 − 1
1/6
6
6
 = [1 + 
6 + 
]
3
3
̃

,  = ̃ 22 ,  =  + 2 ( 6 − )
̃

̃ =  + 2 2 2 ,  = √2 Ω Ω
{
где  -- расстояние до ближайшей стенки, Ω -- тензор вращения.
Опыт использования модели Спаларта-Аллмареса показал, что ее
реальные возможности значительно шире, чем предполагалось при ее создании.
Более того, после введения в нее поправок на кривизну линий тока и вращение,
границы применимости модели заметно расширились. Модель является
удовлетворительной для многих инженерных приложений. В особенности она
применима для расчета обтекания профилей и крыльев, для которых она была
калибрована.
Резюмируя,
следует
отметить,
что
класс
моделей
с одним
дифференциальным уравнением обладает большой приемлемостью к описанию
турбулентных течений с учетом сжимаемости, кривизны линий тока и отрыва
потока. Однако, объектами их приложения, как правило, являются простые
конфигурации потоков с минимальным набором структурных элементов.
Модели с двумя уравнениями
Модели турбулентности с двумя дифференциальными уравнениями
являются наиболее представительной группой дифференциальных моделей.
Первая модель такого типа была предложена в классической работе Колмогорова
[27]. В качестве одного из уравнений все развитые модели, также как и Модель
Колмогорова, используют уравнение переноса  -- кинетической энергии
турбулентных пульсаций. Причиной применения этого уравнения является то, что
оно строго следует из уравнений Навье-Стокса, а также то, что для его замыкания
необходимо промоделировать только два члена: диффузионный и диссипативный
[4]. В качестве примера рассмотрим  −  и  −  модели.
 −  модель Лаундера-Шармы (Launder-Sharma, 1974)
Моделирование турбулентности осуществляется на основе стандартного
или низкорейнольдсового вариантов  −  модели. Уравнения модели включают
формулу Колмогорова-Прандтля для турбулентной вязкости и уравнения
переноса кинетической энергии турбулентных пульсаций и скорости ее
диссипации.
Кинематическая турбулентная вязкость
 =    2 
(27)
Кинетическая энергия турбулентности
∂
∂
∂
∂
∂
∂
+  ∂ =  ∂  −  + ∂ [( +  / ) ∂ ]



(28)

Степень диссипации турбулентности
∂
∂
 ∂
+  ∂ = 1 1   ∂ − 2 2
∂


2

∂
∂
+ ∂ [( +  / ) ∂ ]


(29)
Коэффициенты замыкания
1 = 1.44, 2 = 1.92,  = 0.09,  = 1.0,  = 1.3
При постановке граничных условий на твердых поверхностях в рамках
стандартной  −  модели предполагается, что в турбулентном пограничном слое
имеет место универсальный логарифмический профиль скорости. При этом
полагается  = 1 = 2 = 1.
При использовании низкорейнольдсового варианта  −  модели, вблизи
твердой поверхности принимается
∂
= 0,
∂
3/4  3/2
 Δ
 = 
(30)
а также
 =  −2.5/(1+0.02) ,
1 = 1.0,
где  =  2 / -- турбулентное число Рейнольдса.
 −  модель Вилкокса (Wilcox, 1998) [39]
2
2 = 1 − 0.3  
Данная модель может применяться как в потоках с твердыми стенками, так
и в течениях без касательных напряжений. Использование  −  модели
оказывается численно более устойчивым, чем упомянутых  −  моделей.
Соответствующие определяющие уравнения следующие
Кинематическая турбулентная вязкость
 = 
(31)
где  -- удельная диссипация турбулентности.
Кинетическая энергия турбулентности
∂
∂
∂
∂
∂
∂
+  ∂ =  ∂  −  ∗   + ∂ [( +  ∗  ) ∂ ]



(32)

Скорость диссипации
∂
∂
∂

∂
∂
∂
+  ∂ =    ∂ −  2 + ∂ [( +  ) ∂ ]



(33)

Коэффициенты замыкания и вспомогательные соотношения
 = 0.52, 0 = 0.072, 0∗ = 0.09,  = 0.5,  ∗ = 0.5
Ω Ω 
 = 0  ,  = |
(0 )3
1 ∂ ∂
 ∗ = 0∗ ∗ ,  = 3 ∂

1+70
|,  = 1+80
∂

1,
, ∗ = {
1+6802
1+4002
 ≤ 0
,  > 0
{
где составляющие осредненных
деформации определяются так
1 ∂
∂
Ω ≡ 2 (∂ − ∂ ),


тензоров
вращения
1 ∂
и
скоростей
∂
 ≡ 2 (∂ + ∂ )


1.1.5 Выводы
Исходя из особенностей описанных подходов, и учитывая имеющиеся
вычислительные возможности, в данной работе в качестве модели движения
вязкой несжимаемой жидкости выбраны осредненные по Рейнольдсу уравнения
Навье-Стокса.
В качестве замыкающей модели турбулентности, для представляющих
интерес режимов обтекания нестационарных движущихся объектов, выбрана
модель Спаларта-Аллмараса, относящаяся к классу низкорейнольдсовых. В
результате не возникает необходимости применения пристеночных функций, что
позволяет избежать связанных с ними сложностей при моделировании отрывных
течений.
1.2
Численные
гидродинамики
методы
решения
задач
вычислительной
1.2.1 Дискретизация определяющих уравнений
Выбор численной схемы дискретизации закона сохранения зависит от
формы в которой этот закон представлен [17].
В качестве примера рассмотрим модельное уравнение переноса величины
. Интегральная форма записи этого закона сохранения будет иметь вид
∂
∫
∂ Ω
 Ω + ∫ ⃗ () ⋅ ⃗⃗  = 0
(34)
Применяя теорему о дивергенции, и устремляя объем области Ω к нулю,
можно перейти к дифференциальной форме записи
∂
∂
 + ∇ ⋅ ⃗ () = 0
(35)
Домножая это выражение на произвольную функцию  и интегрируя по
частям, можно получить вариационную форму закона сохранения
∂
∫
∂ Ω
 Ω − ∫Ω ∇ ⋅ ⃗ () Ω + ∫ ⃗ () ⋅ ⃗⃗  = 0
(36)
Дискретизация дифференциальной формы записи осуществляется при
помощи
конечно-разностных
методов.
Интегральная
форма
записи
дискретизируется при помощи метода конечных объемов, а вариационная форма
--- при помощи метода конечных элементов.
2.1.1 Метод конечных разностей: FDM
Метод конечных разностей (FDM --- Finite Differences Method) является
одним из самых первых и самых простых методов решения дифференциальных
уравнений, особенно в случае его использования в задачах с простой геометрией.
Наиболее часто метод применяется на регулярных сетках, когда линии
координатной сетки служат локальными координатными линиями. В методе
конечных разностей исходное дифференциальное уравнение аппроксимируются
системой линейных алгебраических уравнений, где неизвестными являются
значения переменных решения в узлах сетки. При этом Каждое слагаемое
исходного дифференциального уравнения представляется соответствующим
конечно-разностным
отношением.
В
результате
получается
система
алгебраических уравнений относительно неизвестных узловых значений.
В принципе, метод конечных разностей может быть применен к любому
типу сетки. На регулярных сетках метод конечных разностей очень прост и
эффективен. Особенно просто в этом случае получить схемы более высокого
порядка точности. Недостатком конечно-разностных методов является то, что
законы сохранения не учитываются без специальной заботы об этом. Кроме того,
имеет место ограничение по сложности геометрии расчетной области.
Метод конечных объемов: FVM
Метод конечных объемов (FVM --- Finite Volume Method) представляет
собой главный способ решения связанных уравнений переноса импульса и
турбулентности [22].
В отличие от метода конечных разностей, данный метод использует
формулировку уравнений в интегральной форме. Расчетная область разбивается
на определенное количество контрольных объемов (ячеек), каждому из которых
сопоставляется неизвестная величина, представляющая собой среднее значение
переменной по этому объему. Для того, чтобы получить алгебраическое
уравнение, соответствующее интегральному, записанному для некоторой
контрольной ячейки, необходимо необходимо осуществить два этапа
аппроксимации
• Приближенные значения интегралов, входящих в уравнение, при
помощи квадратурных формул, выражаются через значения подынтегральных
выражений в точках границы.
• Значения переменных в точках границы ячейки интерполируются по их
значениям, заданным в узловых точках.
Интегральное уравнение выполняется как для каждого отдельного
контрольного объема в отдельности, так и для расчетной области в целом. Таким
образом метод конечных объемов обладает свойством глобального сохранения,
что является важным преимуществом этого метода.
Метод конечных объемов может применяться с любым типом сетки, так что
он применим для сложных геометрий. Сетка определяет только границы
контрольного объема и не нуждается в привязке к системе координат. По
сравнению с методом конечных элементов, метод конечных объемов более
приемлем для большинства программистов, менее сложен с математической
точки зрения и требует меньшей памяти компьютера при том же числе расчетных
узлов.
Метод обладает преимуществами несложного программирования,
математической простоты и физической адекватности. Вследствие этих
достоинств, большинство разработанных коммерческих программ численного
решения задач гидродинамики используют метод конечных объемов.
Метод конечных элементов. FEM
Метод конечных элементов (FEM --- Finite Element Method) во многом
подобен методу конечных объемов [20]. Область разбивается на раздельные
объемы, или конечные элементы, которые в общем случае являются
неструктурированными. В двумерном пространстве это обычно треугольники или
четырехугольники, а в трехмерном наиболее часто используются тетраэдры или
шестигранники.
В методе конечных элементов искомая функция аппроксимируется
линейной комбинацией координатных функций. Для получения уравнений метода,
исходные уравнения интегрируются по всей расчетной области с некоторым
весом, в качестве весовых функций принимаются координатные функции. В самых
простых методах используются функции формы, линейные в пределах каждого
элемента, что гарантирует непрерывность решения на границах элементов.
Важным преимуществом методов конечных методов является их
применимость для задач со сложными пространственными конфигурациями.
Данные методы относительно просто анализировать математически.
1.2.2 Разбиение расчетной области на элементы
Для численной реализации того или иного метода необходимо
сгенерировать сетку для дискретизации определяющих уравнений. Процесс
построения сетки относится к ключевым моментам численного моделирования,
так как рациональным выбором сетки можно значительно упростить решение
уравнений [20].
Регулярные сетки
Традиционно, при решении задач вычислительной гидродинамики
применялись регулярные сетки (четырехсторонние на поверхности и
шестигранные в пространстве). Регулярность заключается в том, что сетка
представляет собой упорядоченную по определенным правилам структуру данных
с явно выраженными сеточными направлениями, которые, в общем случае,
образуют
криволинейную
систему
координат.
В
преобразованном
(вычислительном) пространстве ячейки сетки являются топологическими
прямоугольниками (двумерные задачи) или параллелепипедами (трехмерные
задачи). Для дискретизации уравнений Навье-Стокса на регулярных сетках
используют, как правило, метод конечных разностей.
При построении регулярных сеток в геометрически сложных областях
применяют преобразование координат общего вида, связанных с поверхностью
тела. Основная цель данного преобразования заключается в построении
равномерной расчетной сетки в преобразованном пространстве таким образом,
чтобы координатные линии совпали с границами физической области. Введение
криволинейной системы координат общего вида повышает эффективность
расчетов в вычислительной области, так как система координат является
прямоугольной. Вместе с тем, при записи уравнений в криволинейных
координатах, появляются дополнительные члены --- параметры преобразования,
определяющие отображение физической области на пространство обобщенных
координат.
Следует особо отметить случаи ортогональных и конформных сеток. В
первом случае при дискретизации модельных уравнений обнуляются некоторые
параметры преобразования --- компоненты метрического тензора преобразования
(матрицы Якоби), находящиеся не на главной диагонали данного тензора.
Следствием является уменьшение погрешности, и, следовательно, повышение
точности решения. Использование конформных преобразований позволяет
сохранить такую же структуру модельных уравнений, записанных в
вычислительной системе координат, как и в декартовом пространстве.
Для построения регулярных сеток применяются следующие методы
• Алгебраические методы
• Дифференциальные методы
• Методы теории функций комплексного переменного
Из анализа литературы последних лет следует, что наиболее широко
распространенными являются методы построения сеток на основе решения
эллиптических уравнений и алгебраические методы --- трансфинитная
интерполяция и метод многих поверхностей. Однако, ни один из
вышеперечисленных подходов не является универсальным по отношению ко
всему спектру задач вычислительной гидродинамики.
Блочные сетки
Для
структурированных
сеток
сравнительно
легко
реализуются
вычислительные алгоритмы на основе современных методов высокого порядка
точности.
Однако,
диапазон
геометрических
объектов,
описываемых
структурированными сетками, ограничен. Как правило, невозможно построить
единую сетку для всей расчетной области, в связи с чем производится
разделение поля течения на подобласти, в каждой из которых генерируется своя
сетка регулярной структуры. Блочный подход предоставляет широкие
возможности для использования эффективных численных методов внутри
отдельных блоков. Основной недостаток блочного подхода состоит в достаточно
сложной процедуре сшивки решений, полученных в различных подобластях.
Можно выделить следующие подходы
• Многоблочные структуры. Физическая область разбивается на несколько
зон или блоков, причем границы блоков могут не соответствовать границам
физической области. Затем для каждого блока отдельно строится сетка в
соответствии с граничными условиями для каждой подобласти. Сетки из разных
блоков могут как стыковаться точно по поверхностям раздела физической области
на зоны, так и пересекаться между собой.
• Иерархические блочные структуры. Данный метод подразумевает
иерархическую вложенность блоков сетки друг в друга. Нижестоящие по иерархии
сетки ''погружены'' в вышестоящие.
Неструктурированные сетки
Характерной
особенностью
неструктурированных
сеток
является
произвольное расположение узлов сетки в физической области. Произвольность
расположения узлов понимается в том смысле, что отсутствуют выраженные
сеточные направления и нет структуры сетки, подобной регулярным сеткам.
Число ячеек, содержащих каждый конкретный узел, может изменяться от узла к
узлу. Узлы сетки объединяются в многоугольники (двумерный случай) или в
многогранники (трехмерный случай). Как правило на плоскости используются
треугольные и четырехугольные ячейки, а в пространстве --- тетраэдры и призмы.
Основное преимущество неструктурированных сеток перед регулярными состоит
в большей гибкости при дискретизации физической области сложной формы, а
также
в
возможности
полной
автоматизации
их
построения.
Для
неструктурированных сеток легче реализуются локальные сгущения и адаптация
сетки в зависимости от поведения решения. Современные программы генерации
сеток позволяют за приемлемое время строить сетки для сколь угодно сложных
геометрических объектов. Для дискретизации уравнений Навье-Стокса на
неструктурированных сетках применяются методы конечных элементов и
конечных объемов. Метод конечных разностей к таким сеткам неприменим.
Недостатком неструктурированных сеток по сравнению с регулярными
сетками является необходимость хранить информацию о структуре сетки, что
повышает требования к вычислительным ресурсам.
Тем не менее, данный подход
распространение по следующим причинам
получает
все
более
широкое
• Процесс генерации сетки сравнительно легко формализуется, и
автоматизируется.
Генераторы
неструктурированных
сеток
обладают
универсальностью по отношению к широкому диапазону прикладных задач.
• В случае очень сложной геометрии расчетной области, время генерации
неструктурированной сетки на порядок меньше времени генерации регулярной
сетки.
• Измельчение произвольных участков сетки осуществляется
естественным образом, что дает преимущество при использовании адаптивных и
многосеточных методов.
Как показывает практика, выбор структуры для представления
триангуляции оказывает существенное влияние на трудоемкость алгоритмов,
использующих данную структуру, а также на скорость конкретной реализации [13].
Рассмотрим структуры данных, применяемые для двумерных треугольных сеток.
• ''Узлы с соседями''. Для каждого узла триангуляции хранятся его
координаты на плоскости и список указателей на соседние узлы, с которыми есть
общие ребра. Треугольники при этом не представляются вообще, что является
обычно существенным препятствием для дальнейшего применения триангуляции.
• ''Узлы, ребра и треугольники''. В явном виде задаются все виды
объектов триангуляции узлы, ребра и треугольники. Недостатком данной
структуры является большой расход памяти.
• ''Узлы и треугольники''. Для каждого треугольника хранятся три
указателя на образующие его узлы и три указателя на смежные треугольники.
Данная структура триангуляции наиболее часто применяется на практике в силу
своей компактности и относительного удобства в работе.
Гибридные сетки
Гибридные сетки комбинируют регулярные и неструктурированные области
сетки. Данный подход позволяет сочетать достоинства и снизить влияние
недостатков, присущих каждому типу сеток. Простейший пример гибридной сетки -- течение около системы профилей. Область около профилей покрывается
ортогональной регулярной сеткой, а области между профилями и далекое поле --неструктурированной.
Вычислительный
алгоритм
содержит
процедуру
переключения вычислительных схем на различных сетках, и, при необходимости,
переноса информации с одного типа сетки на другой.
1.2.3 Уравнения Навье-Стокса для несжимаемой жидкости
Для несжимаемой жидкости выполняются свойства
 = 0,
 = inf
где  -- число Маха,  -- местная скорость звука. Предположение о
несжимаемости является также хорошим приближением, например для воздуха
при скорости  < 100 м/с , или числе Маха  < 0,3.
При рассмотрении задач о несжимаемой жидкости используется частный
случай уравнений для сжимаемой жидкости. Эти уравнения, в случае отсутствия
массовых сил и подвода тепла извне, записываются следующим образом
∇ ⋅ ⃗
+ ⃗ ⋅ ∇⃗
∂
⃗⃗
∂
=
=
1
−
0
∇ + Δ⃗
(37)
Эти уравнения образуют смешанную эллиптически-параболическую
систему относительно неизвестных (, ⃗). Моделирование вихревых течений на
основе численного решения уравнений Навье-Стокса для несжимаемой жидкости
сопровождается рядом трудностей математического характера.
Уравнение неразрывности для несжимаемой жидкости содержит лишь
составляющие скорости, в связи с чем нет прямой связи с давлением, которая для
сжимаемых течений осуществляется через плотность. Вследствие этого одной из
главных проблем является наложении условия несжимаемости, и возможное
возникновение неустойчивости решения при наложении условия несжимаемости.
Другой проблемой при решении системы уравнений Навье-Стокса является
нелинейность, связанная с конвективными слагаемыми в уравнениях движения,
которая может приводить к появлению осцилляций решения в областях с
большими градиентами. В случае превалирования конвекции над диффузией
происходит ухудшение решения из-за жесткости, которую вносят в систему
уравнений конвективные члены, и несимметричности матрицы системы линейных
алгебраических уравнений. Наиболее эффективными методами учета
конвективных слагаемых являются противопотоковые схемы, реализованные
методом конечных объемов.
Метод завихренности и функции тока
В методе завихренности и функции тока [14], выполняется замена
переменных для перехода от компонент скорости (, ) к завихренности  и
функции тока , которые в двумерных декартовых координатах определяются так
∂
∂
 = ∂ − ∂ ,
∂
∂
∂
= ,
∂
= −
(38)
Используя новые независимые переменные, можно скомбинировать
уравнения сохранения импульса, исключая из них давление, что дает
∂
∂
∂
∂
∂2 
∂2 
+  ∂ +  ∂ =  ( ∂ 2 + ∂ 2 )
(39)
Это параболическое уравнение в частных производных называется
уравнением переноса завихренности.
Зависимости (38) позволяют получить уравнение для  и 
∂2 
∂ 2
∂2 
+ ∂ 2 = −
(40)
Это выражение представляет собой уравнение Пуассона, являющееся
эллиптическим.
Итак, в результате замены переменных, смешанная эллиптическипараболическая система уравнений разделилась на одно параболическое и одно
эллиптическое уравнение. Обычно эти уравнения решают методом установления
по времени.
Чтобы определить давление, необходимо решить уравнение Пуассона для
давления
∂2 
∂2 
+ ∂ 2 = 2[∂2 () ∂ 2 ∂2 () ∂ 2 ]
∂ 2
(41)
Распространение подхода с использованием завихренности и функции тока
в качестве независимых переменных на трехмерные задачи осложнено тем, что
для действительного трехмерного течения нельзя ввести функцию тока.
Проекционные методы и метод коррекции давления
Конечно-разностные и конечно-объемные методы решения делятся на
методы, использующие процедуру коррекции давления (pressure-based algorithm),
и методы, основанные на принципе расщепления неизвестных (pressure-velocity
coupling) [26].
Метод расщепления по физическим факторам подразумевает разделение
системы уравнений Навье-Стокса на последовательность более простых
уравнений, таких как уравнения диффузии, переноса и уравнение Пуассона.
Разработка численных методов для этих уравнений оказывается значительно
проще, чем непосредственного для уравнений Навье-Стокса.
Фактически, производится расчет поля скоростей в два этапа. На первом
этапе вычисляется промежуточное поле скоростей без учета уравнения
неразрывности. На втором этапе осуществляется коррекция поля скоростей,
чтобы обеспечить выполнение уравнения неразрывности. Тем самым
осуществляется ''проекция'' поля скоростей на пространство векторов с нулевой
дивергенцией, поэтому методы называются проекционными.
Схема расщепления по физическим факторам (метод проекции),
применяемая для дискретизации уравнений Навье-Стокса, записанных в
физических переменных [7].
Пусть в момент времени   известны поле скорости ⃗ и поле давления .
Тогда для расчета неизвестных функций в момент времени  +1 используется
схема расщепления [6], состоящая из трех этапов.
На этапе 1 предполагается, что перенос
осуществляется только за счет конвекции и диффузии
⃗⃗ ∗ −
⃗⃗ 

Δ
= −⃗  ⋅ ∇⃗  + Δ⃗  + ⃗ 
количества
движения
(42)
Несмотря на то, что промежуточное поле скорости ⃗ ∗ не удовлетворяет
уравнению неразрывности, оно имеет физический смысл, так как сохраняет
вихревые характеристики  ⃗ ∗ =  ⃗ +1 .
На этапе 2 по найденному промежуточному полю скорости ⃗ ∗ , с учетом
соленоидальности вектора скорости ⃗ +1 , рассчитывается поле давления
1
Δ+1 = Δ ∇⃗ ∗
(43)
Для решения уравнения Пуассона (43) на каждом шаге по времени могут
использоваться как итерационные, так и прямые методы.
На этапе 3 предполагается, что перенос количества движения
осуществляется только за счет градиента давления (конвекция и диффузия
отсутствуют)
⃗⃗ +1 −
⃗⃗ ∗

Δ
= −∇+1
(44)
Следует отметить, что уравнение Пуассона (43) получено путем взятия
дивергенции от обеих частей равенства (44) с учетом уравнения неразрывности
∇⃗ +1 = 0.
Метод искусственной сжимаемости
Рассматривается формулировка метода искусственной сжимаемости, в
которой в уравнение неразрывности включено слагаемое, содержащее
производную давления по искусственному времени [26, 22]. Полученное
уравнение вместе с уравнениями сохранения количества движения образуют
гиперболическую систему уравнений, для которой скорость распространения
возмущений давления конечна.
Преимущество метода при решении задач стационарного течения,
заключается в том, что он не требует получения поля скоростей с нулевой
дивергенцией на каждой итерации. Дивергенция поля скоростей автоматически
обращается в ноль, когда решение устанавливается во времени. Кроме того,
данный метод может быть обобщен на случай нестационарного движения.
За счет введения искусственного
неразрывности приводится к виду
∂
1 2 ∂ +
уравнения
∂( )
∂
=0
состояния,
уравнение
(45)
Параметр , представляющий собой искусственную скорость звука должен
быть выбран с учетом двух факторов. С одной стороны, искусственная скорость
звука должна быть близка к скорости движения жидкости, обеспечивая лучшую
обусловленность системы. С другой стороны, искусственная скорость звука
должна быть достаточно велика, чтобы возмущения поля давлений успевали
распространяться достаточно далеко до того, как произойдет их затухание. Тем
самым осуществляется переход к установившемуся состоянию, при котором
выполняется условие несжимаемости. Для обезразмеренных уравнений этим
условиям удовлетворяет диапазон значений  от 0.1 до 10 , предложенный в
работе [29].
С точки зрения линейной алгебры, система линейных уравнений,
возникающая при дискретизации уравнений метода искусственной сжимаемости,
оказывается хорошо обусловленной (при правильном выборе параметра  ) по
сравнению с исходными уравнениями. Любые другие методы решения должны
использовать специальные приемы для получения корректных распределений
давления.
В случае применения метода искусственной сжимаемости для расчета
нестационарных течений, применяется метод двойных шагов по времени (dual
time step) [33].
В этом случае интегрирование по физическому времени  осуществляется
только для уравнения сохранения количества движения. На каждом шаге
физического времени осуществляется итерационный процесс интегрирования по
искусственному времени  следующей системы
∂
∂( )
∂
∂
∂(  )
1 2 ∂ +
∂
∂
+
∂
+
∂
=
=
0
1 ∂
−
∂
+
∂2 
(46)
∂2
Таким образом, на каждом шаге физического времени обеспечивается
выполнение условия несжимаемости.
1.2.4 Методы решения СЛАУ
Одним из наиболее емких, с точки зрения затрат машинных ресурсов,
этапов вычислительной процедуры является решение системы линейных
алгебраических уравнений (СЛАУ) большой размерности. Необходимость
решения такой системы возникает в следующих двух случаях
• в результате применения
определяющих уравнений по времени,
неявных
схем
для
дискретизации
• в результате дискретизации уравнения Пуассона при использовании
проекционных методов наложения условия несжимаемости.
Рассмотрим систему линейных алгебраических уравнений, записанную в
виде
 = 
(53)
где  = { } - матрица системы, имеющая размерность  ×  , вектор
 = (1 , 2 , … ,  ) - вектор решения,  = (1 , 2 , … ,  ) - вектор правых частей.
Численные методы решения систем данного вида [3, 8] принято разделять
на два класса: прямые методы и итерационные.
Методы, позволяющие получить решение системы уравнений за конечное
число арифметических операций, называются прямыми. К ним относятся метод
Крамера, метод исключения Гаусса, метод -разложения и ряд других методов.
Основным недостатком прямых методов являются жесткие требования к
быстродействию и памяти. Например, метод Гаусса требует выполнения порядка
3 арифметических операций и хранения порядка 2 переменных. Кроме того,
если  велико, то машинная погрешность вычислений будет оказывать
значительное влияние на конечный результат.
Итерационный метод общего вида [2, 8] основан на последовательном
улучшении начального приближения решения. Построение последовательности
приближений
осуществляется
посредством
единообразного
процесса,
называемого процессом итераций. Вычислительный процесс заканчивается, когда
изменение решения при переходе к следующей итерации становится достаточно
малым, или когда невязка уменьшается до заданного значения. Итерационные
методы требуют хранения порядка  переменных, а время решения зависит от
обусловленности матрицы и качества начального приближения. Например для
методов Якоби и Гаусса-Зейделя количество арифметических операций
составляет порядка  2 , где  -- количество проведенных итераций. Особенностью
итерационных методов является необходимость исследования сходимости
каждого метода.
Общая структура итерационных методов связана с представлением
матрицы в виде  =  −  и видоизмененной формой исходного уравнения
(  −  )  =  . Различные итерационные методы отличаются друг от друга
способом выбора матрицы  . Примером обычных итерационных методов могут
служить метод Якоби (метод простых итераций), метод Зейделя, метод
последовательной релаксации. К отдельному классу следует отнести
вариационные итерационные методы: метод наискорейшего спуска, минимальных
невязок, минимальных поправок, минимальных погрешностей и другие.
Классические итерационные методы
Рассмотрим для системы (53) итерационный процесс общего вида,
заданный следующим образом
 +1 =   −   (    −  )
где   -- невырожденная матрица. Перепишем выражение в виде
 +1 =     +   
(54)
где   =  −    -- оператор -го шага итерационного процесса,  -единичная матрица.
Пусть  ∗ -- точное решение системы. Введем обозначения
  =   −  ∗,
где   -- вектор ошибки,
место соотношения
 =  − 
  -- вектор невязки. В таком случае имеют
   =  ,
 +1 =    
Итерационный
процесс
(54)
называется
сходящимся,
если
последовательность {   } сходится к  ∗ при любом начальном приближении
 0 . Итерационный процесс называется стационарным, если   не зависит от
номера итерации . В противном случае процесс называется нестационарным.
Методы вариационного типа
В методах вариационного типа [2, 8] решение линейной алгебраической
системы заменяется эквивалентной экстремальной задачей. Пользуясь
обозначениями (53), образуем квадратичный функционал следующего вида
Φ(  ) = (   ,  ) − 2(  ,  )
(55)
где символ (⋅,⋅) обозначает скалярное произведение двух векторов. Можно
показать, что если матрица  симметрична и положительно определена, то
задача решения системы линейных уравнений и задача минимизации
функционала (55) эквивалентны. Для случая произвольной невырожденной
матрицы  , не являющейся симметричной, можно сформулировать аналогичную
задачу минимизации функционала.
Конструирование
применения к задаче
функционала.
итерационного
процесса
(55) различных методов
осуществляется
путем
численной минимизации
Итерационные методы Крылова
[8] Наиболее эффективными и устойчивыми среди итерационных методов
являются
проекционные
методы,
связанные
с
проектированием
на
подпространства Крылова (Krylov subspace methods). По сравнению с
классическими итерационными методами, они не содержат эмпирически
подбираемых параметров и позволяют получить более высокую скорость
сходимости, несмотря на увеличение числа операций на каждой итерации.
В семейство итерационных методов Крылова входят, в частности
• обобщенный метод минимальных невязок (Generalized Minimum
Residual, GMRES),
• метод сопряженных градиентов (Conjugate Gradients, CG),
• метод квадратичных сопряженных градиентов (Conjugate Gradients
Squared, CGS),
• метод бисопряженных градиентов (Bi-Conjugate Gradients, BiCG),
• метод бисопряженных градиентов со стабилизацией (BiCGStab).
В программном коде «SmartFlow» используются обобщенный метод
минимальных невязок (GMRES) и метод бисопряженных градиентов (BiCG).
1.2.5 Выводы
Исходя из поставленной задачи расчета течений в сложных областях с
подвижными границами, для разделения расчетной области на элементы
наиболее
перспективным,
с
точки
зрения
автора,
оказывается
неструктурированный тип сетки.
Для пространственной дискретизации определяющих уравнений выбран
метод конечных объемов с высоким порядком аппроксимации, который, наряду с
другими достоинствами, хорошо зарекомендовал себя при расчетах на
неструктурированных сетках.
Для наложения условия несжимаемости в уравнениях Навье-Стокса выбран
подход
с
искусственной
сжимаемостью.
Результирующие
уравнения
дискретизируются по времени при помощи неявной схемы второго порядка
точности. Возможность применения метода искусственной сжимаемости для
решения нестационарных задач обеспечивается за счет применения метода
двойных шагов по времени.
Учитывая, что система линейных алгебраических уравнений, возникающая
в результате применения неявной схемы дискретизации по времени, имеет
большую размерность и обладает разреженной несимметричной матрицей, для
решения этой системы выбран итерационный метод бисопряженных градиентов.
Численная реализация метода решения уравнений Навье-Стокса
2.1 Общая постановка задачи и основные допущения
2.2 Особенности метода расчета
2.3 Дискретизация по пространству. Метод конечных объемов
2.3.1 Аппроксимация невязкого потока
2.3.2 Полиномиальная аппроксимация
2.3.3 Аппроксимация вязкого потока
2.3.4 Граничные условия
2.4 Неявная дискретизация по времени
2.4.1 Двойной шаг по времени
2.4.2 Одинарный шаг по времени
2.4.3 Построение матрицы СЛАУ
2.5 Расчет турбулентной вязкости
2.5.1 Модель Спаларта-Аллмараса
2.5.2 Дискретизация по пространству и времени
2.6 Расчет гидродинамических характеристик
2.7 Программная реализация. Расчетные алгоритмы
2.7.1 Алгоритм расчета нестационарного течения
2.7.2 Решение СЛАУ
2.7.3 Краткая характеристика вычислительных средств
3 Верификация расчетного метода
3.1 Оценка сходимости метода
3.2 Результаты тестовых расчетов
3.2.1 Задача турбулентного обтекания плоской пластины
3.2.2 Задача обтекания двумерного кругового цилиндра в канале
3.2.3 Задача обтекания двумерного канала с обратным уступом
3.2.4 Задача обтекания двумерного квадратного цилиндра в неограниченной
жидкости
3.2.5 Задача обтекания двумерного кругового цилиндра, совершающего
заданные колебания поперек потока
4 Исследование пропульсивных характеристик машущего крыла
4.1 Постановка задачи обтекания машущего крыла на упругих связях
4.2 Пропульсивные характеристики машущего крыла
4.3 Анализ результатов расчета
движителя типа машущее крыло
пропульсивных
характеристик
4.3.1 Угловые заданные колебания профиля
4.3.2 Вертикально-угловые заданные колебания профиля
4.3.3 Колебания жесткого профиля на упругих связях
Заключение
Литература
[1]
Андерсон Д., Таннехил Дж.,
гидромеханика и теплообмен. М.: Мир, 1990.
Плетчер
Р.
Вычислительная
[2]
Баландин М.Ю., Шурина Э.П. Методы решения СЛАУ большой
размерности. // Новосибирск: Изд­во НГТУ, 2000.
[3] Бахвалов Н.С., Жидков Н.П., Кобельков Г.М. Численные методы. // М.:
Наука, 1987.
[4] Белов И.А., Исаев С.А. Моделирование турбулентных течений: Учебное
пособие. // СПб: БГТУ, 2001.
[5]
Белов И.А., Исаев С.А., Коробков В.А. Задачи и методы расчета
отрывных течений несжимаемой жидкости // Л.: Судостроение, 1989.
[6] Белоцерковский О.М. Численное моделирование в механике сплошных
сред. М.: Физматлит, 1994.
[7] Волков К.Н. Реализация схемы расщерления на разнесенной сетке для
расчета
нестационарных
течений
вязкой
несжимаемой
жидкости.
//
Вычислительные методы и программирование. 2006. 269-282.
[8] Голуб Дж., Ван Лоун Ч. Матричные вычисления. // М.: Мир, 1999.
[9]
Елизарова Т.Г., Калачинская И.С., Шеретов Ю.В., Шильников Е.В.
Численное моделирование отрывных течений за обратным уступом.
[10] Лойцянский Л.Г. Механика жидкости и газа. // М.: ГИТТЛ, 1970.
[11] Пейре Р., Тейлор Т.Д. Вычислительные методы в задачах механики
жидкости. // Л.: Гидрометеоиздат, 1986.
[12] В.А., Тарасов С.В. Метод численного расчета течений вязкой жидкости
с использованием осредненных по Рейнодльдсу уравнений Навье-Стокса //
Тезисы докладов научно-технической конференции ''Проблемы мореходных
качеств судов и корабельной гидромеханики'' (XLII Крыловские чтения). СПб,
2006. 17-19.
[13] Скворцов А.В. Обзор алгоритмов построения триангуляции Делоне. //
Вычислительные методы и программирование. 2002. 14-39.
[14] Флетчер К. Вычислительные методы в динамике жидкостей. М.: Мир,
1991.
[15]
Anderson J.D. Computational Fluid Dynamics. The Basics with
Applications. 1995.
[16] Baldwin B.S., Lomax H. Thin-Layer Approximation and Algebraic Model for
Separated turbulent Flows. // AIAA Paper, 78-257. 1978.
[17] Barth T.J. Aspects of Unstructured Grids and Finite-Volume Solvers for the
Euler and Navier-Stokes Equations. // Unstructured Grid Methods for AdvectionDominated Flows. AGARD, 1992.
[18] Barth T.J., Frederickson P.O. Recent Developments in High Order K-Exact
Reconstruction on Unstructured Meshes. // AIAA paper 93-0668, 1993.
[19] Biedron R.T., Vatsa V.N., Atkins H.L. Simulation of Unsteady Flows Using
an Unstructured Navier-Stokes Solver on Moving and Stationary Grids. // AIAA paper
2005-5093, 2005.
[20] Chung T.J. Computational fluid dynamics. // CUP, 2002.
[21] Fekken G. Numerical Simulation of Free-Surface Flow with Moving Rigid
Bodies. 2004.
[22] Ferziger J.H., Peric M. Computational methods fluid dynamics. // Springer,
2001.
[23] Guilmineau, E. and Queutey, P. A Numerical simulation of vortex shedding
from an oscillating circular cylinder // J. Fluids Struct. Vol. 16. 2002. 773–794.
[24] Hino T. Navier-Stokes Computations of Ship Flows on Unstructured Grids
// Twenty-Second Symposium on Naval Hydrodynamics, 1998. 463-475.
[25]
Hino T. Unsteady Flow Simulation Around a Moving Body by an
Unstructured Navier-Stokes Solver.
[26] Langtangen H.P., Mardal K. Numerical Methods for Incompressible Viscous
Flow.
[27] Kolmogorov A.N. Equations of Turbulent Motion of an Incompressible Fluid.
1942.
[28] Kravchencko A.G., Moin P. Numerical Studies of Flow Over a Circular
Cylinder at ReD=3900 // Physics of Fluids. 2000. Vol. 12. 403-417.
[29]
Kwak D.C., Chang J.L., Shanks S.P., Chakravarthy S.K. A ThreeDimensional Incompressible Navier-Stokes Solver Using Primitive Variables.
[30] Ollivier-Gooch C., Van Altena M. A High Order Accurate Unstructured
Mesh Finite-Volume Scheme for the Advection-Diffusion Equation. 2002.
[31]
Piomelli U., Scotti A., Balaras E. Large-Eddy Simulations of Turbulent
Flows, from Desktop to Supercomputer (Invited Talk). 2002.
[32] Roe P.L. Characteristic-Based Scheme for the Euler Equations // Annual
Review on Fluid Mechanics. 1986. Vol. 18. 337-365.
[33] Rogers S.E., Kwak D. Upwind Differencing Scheme for the Time-Accurate
Incompressible Navier-Stokes Equations // AIAA Journal. 1990. Vol. 28. 253-262.
[34] Rogers S.E., Kwak D. Upwind Differencing Scheme for the Incompressible
Navier-Stokes Equations // Applied Numerical Mathematics. 1991. Vol. 8. 43-64.
[35] Schafer M., Turek S. Benchmark Computations of Laminar Flow Around a
Cylinder // Notes on Numerical Fluid Mechanics. 1996. 856-869.
[36] Sohankar A., Davidson L., Norberg C. Numerical Simulation of Unsteady
Flow Around a Square Two-Dimensional Cylinder // In Proc. 12-th Australasian Fluid
Mechanics Conference, 1995. 517-520.
[37]
Spalart P.R., Allmaras S.R. A One Equation Turbulence Model for
Aerodynamic Flows. AIAA Paper, 92-439. 1992
[38]
Van Altena M. High-Order Finite-Volume Discretisation for Solving a
Modified Advection-Diffusion Problem on Unstructured Triangular Meshes. 1999.
[39] Wilcox D.C. Multiscale Model for Turbulent Flows. // AIAA Journal. Vol. 26.
1988.
[40] Wilcox D.C. Turbulence modelling for CFD. // DCW Industries, 1993.
[41] Yang J., Balaras E. An embedded-boundary formulation for large-eddy
simulation of turbulent flows interacting with moving boundaries. 2005.
Приложения
Иллюстрации
Документ
Категория
Теория систем управления
Просмотров
745
Размер файла
147 Кб
Теги
1/--страниц
Пожаловаться на содержимое документа