Содержание к диссертации
Введение
Глава 1. Аналитический обзор литературы
1.1. Актуальность 7
12. Цели и задачи 10
Глава 2. Физико-математическая постановка
2.1. Основные уравнения сохранения 28
2.2. Граничные и начальные условия 29
Глава 3. Численный метод
3.1. Получение дискретного аналога уравнения конвекции-диффузии..31
3.2. Алгоритм связки «скорость-давление» 36
3.3. Реализация модели турбулентности Spalart-Allraaras 40
3.4. Особенности программной реализации численного метода 42
Глава 4. Численный метод
4.1. Тестовые задачи 44
4.2. Тестирование схемы дискретизации конвективных потоков 45
4.3. Тестирование алгоритма «связки» скорость-давление 51
4.4. Тестирование численного метода и программного комплекса 57
Глава 5. Опыт использования UNSS в прикладных аэродинамических исследованиях
5.1. Исследование компоновки самолета-амфибии Бе-200 74
5.2. Анализ результатов 89
Заключение 90
Литература 92
- Граничные и начальные условия
- Алгоритм связки «скорость-давление»
- Тестирование схемы дискретизации конвективных потоков
- Исследование компоновки самолета-амфибии Бе-200
Введение к работе
Определяющие уравнения ньютоновской гидроаэродинамики -нестационарные уравнения Навье-Стокса - известны уже 150 лет и даже более. Однако разработка укороченных форм этих уравнений по-прежнему остается областью активных исследований, также как и проблема замыкания осредненньтх по Рейнольдсу уравнений Навье-Стокса, применяемых в теории турбулентности. В области неньютоновской гидроаэродинамики уровень описания химически реагирующих потоков и двухфазных течений, пока что не является столь высоким.
Экспериментальная гидроаэродинамика сыграла важную роль при проверке справедливости выведенных зависимостей и установлений пределов пригодности различных приближений в определяющих уравнениях. Аэродинамическая труба (АДТ), применяемая как оборудование для проведения физического эксперимента, является эффективным средством моделирования реальных течений. Традиционно АДТ сыграла роль альтернативы натурным измерениям, позволившей снизить временные и материальные затраты. В процессе проектирования конструкций, геометрия которых напрямую связана с характером поведения потока жидкости или газа, например, при проектировании самолета, натурные измерения не оправдывают себя экономически.
Устойчивое повышение скорости существующих ЭВМ и объема их памяти, начавшиеся в 1950-х годах, привело к возникновению вычислительной гидроаэродинамики (ВГАД). Эта ветвь гидро аэродинамики дополняет ее экспериментальную и теоретическую ветви, представляя собой альтернативное экономически эффективное средство моделирования реальных течений. ВГАД предоставляет возможности численного моделирования таких явлений, экспериментальное моделирование которых крайне сложно и дорогостояще, а зачастую вообще технически невозможно. Например, эксперименты в аэродинамических трубах ограничиваются
определенным диапазоном чисел Рейнольдса, которые оказываются обычно на один или два порядка меньше величины натурных чисел.
Еще одно преимущество вычислительной гидроаэродинамики состоит в том, что при желании можно отбросить те или иные члены определяющих уравнений. Тем самым, открывается путь к опробации новых теоретичесішх моделей.
Появление более эффективных ЭВМ стимулировало интерес к ВГАД. В результате ВГАД стало к настоящему времени наиболее предпочтительным средством проверки качества разработок в авиационной промышленности, в промышленности турбодвигателей и, в несколько меньшей степени автомобильной промышленности.
Граничные и начальные условия
При решении задач большой размерности, если не предпринимаются специальные приемы в решении систем дискретных аналогов уравнения для поправок давления, возникает другая проблема, связанная с наличием остаточных низкочастотных составляющих невязки в поле давления, которые подавляются в течение слишком большого числа глобальных итераций, иначе говоря, решение сходится медленно. Специальными приемами могут быть, либо применение многосеточных методов [2] [21], либо решение систем дискретных аналогов уравнения для поправок давления до очень «жестких» критериев сходимости. Первый вариант сегодня достаточно широко применяется в различных коммерческих пакетах. Однако у многосеточных методов существует ряд недостатков. Основными недостатками являются: проблема сходимости решения системы линейных алгебраических уравнений, если в решении появляются кризисные зоны и сложность реализации, особенно на неструктурированных сетках. Второй вариант: решение систем дискретных аналогов уравнения для поправок давления до очень «жестких» критериев сходимости (10"8 и ниже) неэффективен, даже если удалось создать метод решения с суперлинейной сходимостью, так как рассчитываются поправки к давлению, и минимизация нормализованной невязки ниже значений 10 3 становится бессмысленной.
Дальнейшим развитием алгоритма SIMPLE является его модификация SIMPLER, в которой вместо получения поля давления путем добавления поправки решается уравнение Пуассона для давления. В этом случае удается сэкономить до 40-50% итераций, но приходится решать два эллиптических уравнения вместо одного. К недостаткам метода опять же относится проблема шахматного распределения давления, хотя она стоит не так «остро», как в случае SIMPLE. Алгоритм PISO по своим свойствам, с позиции сходимости, близок к SIMPLER, но также неэкономичен с точки зрения количества вычислений за одну итерацию. В нем опять же приходится решать два уравнения эллиптического типа: уравнения для первых и вторых поправок к давлению.
Наиболее естественным по своей природе является алгоритм искусственной сжимаемости, разработанный Chorin в 1967. Это первый известный алгоритм, разработанный на основе физических представлений. Его идея проста. Нарушение закона сохранения массы в рассматриваемом контрольном объеме выражается в виде источника массы, который при заданном шаге интегрирования по времени пересчитывается в приращение плотности в этом объеме. Далее, полученное приращение плотности выражается в приращении давления, которое связывается с плотностью через коэффициент искусственной сжимаемости, так как рассматривается несжимаемая среда. По полученному полю давления, используя уравнения переноса компонент количества движения, рассчитывается новое поле скоростей, которое лучше удовлетворяет уравнению неразрывности. Так итерации выполняются до полной сходимости решения. Немного позже была разработана модификация метода для сжимаемых течений, в которой плотность с давлением связывалась через уравнение состояния. Практика показала достаточную универсальность этого метода, за что он стал популярным [2], [13]. Метод является локальным при расчете поля давления, то есть явным, и хорошо работает на структурированных качественных сетках. В случае неструктурированных сеток появляются проблемы, связанные с получением устойчивой сходимости решения, особенно при достаточно больших числах CFL.
Для решения системы дискретных аналогов уравнения конвекции-диффузии и уравнения для давления были реализованы метод Якоби и метод сопряженных градиентов с предобусловителем соответственно.
Матрица коэффициентов системы дискретных аналогов уравнения конвекции-диффузии является несимметричной в связи с применением транспортивных схем дискретизации. В таком случае решение этой системы возможно с помощью либо «поточечных» методов Якоби, Гаусса-Зейделя, [20] и их модификаций, либо методов крыловского подпространства таких, как GMRES, CGSTAB, [2],[21]. С целью тестирования эти методы были реализованы в солвере. Наиболее эффективным методом в данном случае оказался метод Якоби благодаря тому, что в матрице коэффициентов выполняется условие жесткого диагонального преобладания на любой расчетной сетке, и потребные критерии сходимости (величина нормализованной невязки) составляют значения не ниже 10"3. Метод Гаусса-Зейделя становится оправданным в случае специальной шахматной нумерации (метод Red-Black) [21], которая в солвере не была использована. Конечно, при использовании GMRES, CGSTAB требовалось значительно меньшее число итераций на решение той же системы до требуемых критериев, но это преимущество сводилось «на нет» громадным количеством арифметических операций за одну итерацию.
Алгоритм связки «скорость-давление»
. Итерации повторяются до достижения заданных критериев сходимости. Алгоритм хорошо работает на широком классе задач и структурированных сетках. При попытке перенести его на рассматриваемые неструктурированные сетки с большой величиной деформации элементов (тонкие призмы в пограничном слое) возникают сложности с получением устойчивого итерационного процесса решения при достаточно больших числах CFL. Это связано с явным подходом в расчете поля давления, т. е. новые значения давления в узлах не связаны с новыми значениями в соседних узлах на текущей итерации.
Получаем общий вид дискретного аналога уравнения для давления. Значения величин в центре массы грани контрольного объема вычисляются, используя формулу (3.3). Запись нестационарного, источникового члена, и диффузионного потока выполняются по методике описанной в разделе 3.1. Таким образом, выполнив дискретизацию всех членов уравнения давления, получаем его дискретный аналог в виде линейного алгебраического уравнения:
В результате модификации классический вариант метода искусственной сжимаемости принимает вид: 1. Задается начальное приближение, т.е. ноля скоростей и давления. 2. Записывается система дискретных аналогов уравнения переноса компонент импульса (3.26), (3.27), (3.28) не включая градиенты давления в источники, а именно рассчитываются компоненты вектора псевдоскорости й, v, ш и коэффициент для давленияd. 3. Формируется и решается система дискретных аналогов уравнения для давления (2.30). 4. Рассчитывается новое поле скоростей, путем решения системы дискретных аналогов уравнения переноса компонент импульса (3.26), (3.27), (3.28) используя в источниках новое поле давления. 5. Возвращаемся к шагу 2. Итерации повторяются до достижения заданных критериев сходимости.
Наличие коэффициента сжимаемости /7 в нестационарном члене позволяет ограничить скорость распространения возмущений в поле давления. Изменяя соответствующим образом коэффициент сжимаемости /3 в течение итерационного процесса решения можно достигнуть быстрой сходимости при высокой численной устойчивости процесса. 3.3. Реализация модели турбулентности Spalart-Alimaras
Так как модель турбулентности Spalart-Alimaras является частным случаем уравнения конвекции-диффузии метод получения дискретного аналога, описанный выше, полностью применим и к ней. Особенностью реализации модели является линеаризация источника [3].
В связи с тем, что диструктивный член У. присутствует в источнике с отрицательным знаком, то возможны ситуации, когда в поле модифицированной турбулентном вязкости v появляются отрицательные значения. В этом случае диструктивный член необходимо сделать зависимым от переменной v.
Как уже было сказано ранее, модель турбулентности Spalart-Alimaras является низкореинольдсовои моделью. Поэтому расчетная сетка в пограничном слое должна разрешать не только турбулентный слой, но и ламинарный подслой до значений Y+ 1. Причем количество слоев призм в ламинарном подслое должно быть не менее 10. Только в этом случае появляется возможность достаточно точного расчета потери импульса и касательных напряжений в пограничном слое. Использование же функции стенки не дает требуемой точности. Разработанный расчетный код является объектным и позволяет реализовать другую необходимую модель турбулентности. В частности, предполагается реализовать модель SST, которая дает несколько более точные результаты, чем модель Spalart-Allmaras, но требующая больше вычислительных ресурсов и более капризная в части сходимости численного решения задачи. 3.4. Особенности программной реализации численного метода
Программный комплекс написан на языке программирования C++ с использованием элементов объектно-ориентированного программирования и является кроссплатформенным, т.е. может быть откомпилирован на платформах Windows, Unix и др. платформах, поддерживающих стандарт C++. В связи с этим никакие дополнительные библиотеки не использовались.
Одними из главных требований, предъявляемых к программному комплексу, являлись: - экономичность использования оперативной памяти; - производительность; - надежность;
С целью достижения этих целей была разработана архитектура программного комплекса, состоящего из двух модулей, препроцессора и солвера (решателя). Препроцессор и солвер представляют собой консольные, независимые приложения. Изначально расчетная сетка создается в каком-либо программном комплексе ICEMCFD, Gambit, Ansys и выгружается в текстовом формате. Затем запускается препроцессор с последующей загрузкой расчетной сетки. Структура и особенности сетки анализируются, формируется текстовый файл статистики, содержащий всю необходимую информацию о количестве и размерах векторов, которые будут выделены в памяти на начальном этапе расчета. Затем запускается сам солвер, и загружаются файлы расчетной сетки, созданный препроцессором файл статистики, файл содержащий граничные, начальные условия и физические параметры среды и файлы выбранного профиля. Выполняется расчет, и расчетные данные выгружаются в текстовый файл для последующей обработки в системе визуализации.
Тестирование схемы дискретизации конвективных потоков
Одними из характерных задач, демонстрирующих недостатки схем, являются задачи конвективного переноса скачка скаляра и П-импульса скаляра. Так как схема ориентирована на стационарные течения, исследуем ее только на схемную диффузию от пространственной дискретизации. Для устранения вклада в схемную диффузию от дискретизации по времени используем схему Эйлера с чрезвычайно малым шагом CFL - 0.01.
Распределения массовой концентрации скаляра. Граф.1 - идеальное ступенчатое распределение. Граф.2 - расчетное распределение по оси у = 0.1, z = 0.1. Граф.3 - расчетное распределение по оси у = 0.5, z = 0.5. Граф.4 - расчетное распределение по оси у = 0.9, z = 0.9. Как видно из результатов расчета переноса скачка скаляра (рис. 4.2.), фронт ступеньки «размазан» иа 3 контрольных объема. Величина дисперсии не превышает 1.2% от амплитуды скачка. Несмотря на неоднородную структуру расчетной сетки в кубе распределения массовой концентрации скаляра по трем осям хорошо совпадают.
Расчет выполнялся на неструктурированной сетке из тетраэдров, заполняющей расчетную область, имеющую форму параллелепипеда размером 3x1x0.075 рис.4.3. Средний шаг сетки по всем координатам составил 0.05.
По всей расчетной области задано конвективное поле/ - {1,0,0}. Начальное условие ф - 0. Граничное условие на входе ф = 1 при 0 t 1.334 и = 0 при t 1.334. Граничное условие на выходе — = 0. dx Нестационарный расчет проведен до времени t = 2.24. При решении этой задачи в расчетном пакете Fluent в качестве уравнения конвективного переноса скаляра использовано уравнение переноса энергии с нулевым источниковым членом и нулевым коэффициентом диффузии.
Рассмотрим результаты расчетов. Схема 1-го порядка аппроксимации (рис 4.4. Граф. 2) дает монотонное решение, но обладает характерной для нее большой схемной диффузией. В результате использования рассматриваемой схемы 3-го порядка аппроксимации (рис 2.2.2. Граф. 3) передний фронт ступеньки «размазан» на 5 контрольных объемов после прохождения 45 контрольных объемов, а задний - на 4 после 18. Величина дисперсии не превышает 2.2% от амплитуды импульса. При сравнении распределений массовой концентрации скаляра, полученных с помощью программы UNSS (рис 4.4. Граф. 3) и расчетного пакета Fluent (рис. 4.4. Граф. 4) выявлено, что для достижения идентичных распределений скаляра по оси X, при использовании UNSS потребуется в 1.8 раза меньше контрольных объемов вдоль той же оси, чем при использовании Fluent. Задача переноса синус-импульса скаляра.
Расчет выполнялся на неструктурированной сетке из тетраэдров, заполняющей расчетную область, имеющую форму параллелепипеда размером 3x1x0.075 рис.4.3. Средний шаг сетки по всем координатам составил 0.05.
С целью тестирования алгоритма решалась задача ламинарного обтекания крыла конечного размаха при больших числах Рейнольдса. При такой постановке задачи контроль над конвективными потоками со стороны диффузионных достаточно мал. В природе в этом случае благодаря потери устойчивости слоистого течения ламинарный поток переходит в турбулентный. В данном расчете этот эффект может проявляться виде неустойчивости решения.
Расчетная модель представляет собой Z - симметричную консоль крыла малого удлинения. Неструктурированная сетка строилась в пакете ICEMCFD. Размерность сетки составила 49049 узлов, 60523 тетраэдров и 72790 призм в 10 слоях. Число узлов на поверхности крыла - 3665.
Характеристики среды; плотность р = 1.225 кг/мЗ, динамическая вязкость JLI = 0.00356 Па с. Скорость на бесконечности Vjnf = 50м/с. Число Re = 12000. Угол атаки крыла составил 4е.
Численная формулировка в программе UNSS следующая. Метод связки «скорость - давление» - последовательный алгоритм на основе метода искусственной сжимаемости. Стационарное решение было получено методом установления.
В расчетном пакете Fluent была использована следующая численная формулировка. Метод связки «скорость - давление» - последовательный алгоритм SIMPLE.. Распределение коэффициента давления на поверхности крыла. Как видно из истории сходимости итерационного процесса решения (рис. 4.6.), тестируемая схема в сочетании с последовательным алгоритмом на основе метода искусственной сжимаемости демонстрирует хорошую устойчивость и достаточно быструю сходимость. Относительные нормализованные невязки для каждой степени свободы достигли своего минимума за 10 итераций. Длительный процесс установления поля течения (особенно, это видно по невязке Vz) обусловлен формированием пелены и отрывного течения на задней кромке. История сходимости решения в пакете Fluent не приведена, так как без введения большой начальной вязкости расчет расходился.
Как видно из рис. 4.8-4.10 распределение коэффициента давления на поверхности крыла, полученное в результате расчета поля течения с помощью программы UNSS, хорошо совпадает с распределением, полученным в пакете Fluent. Следует отметить, что число контрольных объемов в случае постановки в программе UNSS в 2.7 раз меньше, чем в случае постановки в расчетном пакете Fluent, хотя расчетная сетка использовалась одна и та же. Это объясняется разными методами формирования контрольного объема.
Исследование компоновки самолета-амфибии Бе-200
В целях поиска мер по улучшению аэродинамических характеристик самолета Бе-200РР, проведены численные расчеты полей течений вокруг модели компоновки самолета. Исследовалось поле течения в районе бортотсека крыла (центроплан) и перед воздухозаборником. Так же исследовалось поле течения в районе концевой части крыла с целью повышения его несущих свойств. А именно интересовали горизонтальные скосы потока в районе пилона поплавка.
Расчетная модель представляет собой Z-симметричную часть компоновки самолета. С целью упрощения модели было убрано хвостовое оперение исходного самолета, не оказывающие влияние на исследуемые зоны. Форма расчетной области на бесконечности - параллелепипед. Граница расчетной области удалена от компоновки самолета, более, чем на 10 длин фюзеляжа.
Исследования выполнены при натурных числах Рейнольдса на крейсерском режиме полета. Расчетная сетка состояла из 300000 тетраэдров и 1200000 призм, заполняющих пограничный слой. Закон изменения толщин призм логарифмический. Параметр пограничного слоя Y+ не превышал 1. Ниже приведены фрагменты результатов расчетов, которые включают: - поверхностные линии тока (аналог «масляных пленок» в физическом эксперименте); - пространственные линии тока; - распределения коэффициента давления по поверхности модели и в Z -сечениях; - изобары в плоскости симметрии двигателя и на поверхности модели. Рис.5.1. Поверхностные элементы пространственной расчетной сетки і
В результате расчетных исследований были определены мероприятия по повышению летно-технических характеристик Бе-200, с минимальными конструктивными изменениями. Во-первых, это создание несимметричного профиля пилона поплавка, тем самым, подстраивая его под горизонтальный скос потока, благодаря чему устраняется разрыв циркуляции на крейсерском режиме полета и повышается аэродинамическое качество. Во-вторых, это перепрофилировка пилона двигателя для уменьшения разгона потока в районе носка, что видно из распределения коэффициента давления на рис. 5.12, с целью затягивания волнового кризиса по числу Маха. Необходимость этой модификации так же выявлена в экспериментальных исследованиях компоновки в АДТ на трансзвуковом режиме.
При выполнении расчета объем используемой оперативной памяти составил 1.7 Гбайт, Время расчета на одном процессоре Pentium 4 3.06 ГГц около 50 часов. Заключение
Как показала опытная эксплуатация, разработанный программный комплекс UNSS, являющийся основой диссертационной работы, позволил повысить эффективность аэродинамического проектирования. Проведенные расчетные исследования продемонстрировали те качества, которые необходимы для решения задач проектирования летательных аппаратов. Прежде всего, это высокая численная устойчивость численного метода реализованного в UNSS, практически не зависящая от качества расчетной сетки и «умеренные аппетиты» к вычислительным ресурсам.
Конечно же, узлы расчетной сетки должны быть распределены исследователем по расчетной области таким образом, чтобы разрешить все особенности исследуемой геометрии и физические особенности течения. Это, сегодня самый сложный этап в численном эксперименте, требующий от исследователя большого опыта, от уровня выполнения которого, напрямую зависит точность расчета.
Что касается требований к объему оперативной памяти, то в среднем на каждый расчетный узел расчетной задачи требуется около 1.4 Кбайт, что в 2.5 раза меньше чем требует Fluent и в 3.5 раза меньше, чем требует StarCD. Поэтому, приведенный расчет симметричной части компоновки в крейсерской конфигурации можно выполнять на компьютере Pentium 4 с 2 Гбайт оперативной памяти, а в случае использования расчетного паїсета Fluent необходимо использовать, как минимум серверную конфигурацию с двумя процессорами в режиме HyperThreading или четырьмя процессорами и 4 Гбайт оперативной памяти. Другой вариант собрать сетевой кластер из четырех Pentium 4. В данном случае распараллеливание необходимо для возможности проведения расчета во Fluent такого масштаба, так как платформа является 32-разрядной, а адресное пространство ограничено 2 Гбайтами для одного процесса. Использование пакетом Fluent динамического перевыделения понижает размер доступной оперативной памяти до 1.4 Гб, благодаря ее дефрагментации. Таким образом, для решения задачи с потребной оперативной памятью 4.5 Гбайт во Fluent потребуется четыре процесса.
В настоящее время одной из главных нерешенных проблем программного комплекса UNSS является отсутствие возможности распараллеливания расчета. В таком виде UNSS сильно проигрывает на больших задачах другим расчетным пакетам, обладающих распараллеливаемым солвером, во времени счета.
