Содержание к диссертации
Введение
ГЛАВА 1. Постановка задачи течения вязкого сжимаемого газа в плоском канале с внезапным расширением 20
1.1. Система уравнений гидромеханики 20
1.2. Граничные и начальные условия 25
ГЛАВА 2. Разностный метод высокого порядка точности решения уравнений гидромеханики 31
2.1. Метод интегрирования уравнений гидромеханики по времени 31
2.2. Метод интегрирования уравнений гидромеханики по пространственным координатам 39
2.3. Алгоритм расчета пространственных производных в окрестности стенки 47
ГЛАВА 3. Анализ результатов методических расчетов 56
3.1. Обоснование сходимости процесса вычислений при интегрировании уравнений гидромеханики 56
3.2. Точки отрыва и присоединения ламинарных и турбулентных течений за обратным уступом 62
3.3. Пульсационные характеристики дозвукового турбулентного течения за обратным уступом 70
3.4. Расчет одномерного энергетического спектра 75
ГЛАВА 4. Параметрическое исследование сжимаемых течений газа в плоском канале с внезапным расширением 81
4.1. Условия численного моделирования 81
4.2. Симметричные ламинарные стационарные течения 82
4.3. Несимметричные ламинарные стационарные течения 90
4.4. Нестационарные ламинарные течения 97
4.5. Переходные течения в канале с боковыми стенками 102
Заключение 115
Список литературы 117
- Граничные и начальные условия
- Метод интегрирования уравнений гидромеханики по пространственным координатам
- Точки отрыва и присоединения ламинарных и турбулентных течений за обратным уступом
- Несимметричные ламинарные стационарные течения
Введение к работе
Существование резко различающихся ламинарных и турбулентных режимов течения было замечено еще в первой половине XIX в., но начало теории турбулентности положено лишь в конце того же столетия в работах Осборна Рейнольдса [1—5]. Именем именно этого ученого названо впоследствии число, являющееся общим критерием перехода ламинарного течения в турбулентное.
Наиболее распространенной является интерпретация числа Рейнольдса как меры относительной значимости сил инерции и сил вязкости, действующих внутри жидкости или газа. Силы инерции, если они существенно превосходят силы вязкости, что соответствует большим числам Re, вызывают перемешивание конечных объемов газа, движущихся с разными скоростями. В результате осуществляется передача энергии от крупномасштабных структур к менее крупным, образующимся за счет потери устойчивости более крупных вихрей [6 - 11]. Поглощая энергию основного потока, эти структуры оказываются сильно анизотропными, завихренными и существенно отличаются от течения к течению. Поэтому возникает необходимость детального изучения потоков газа в технических системах, так как характер течения может сильно повлиять на работоспособность устройства и иные его характеристики. Процесс потери устойчивости и переход к турбулентному течению происходит практически скачкообразно, следовательно, важно знание параметров, при которых наступает этот переход, и где та граница, при превышении которой происходит разрушение существующего течения.
Ламинарные потоки, по сравнению с турбулентными, наверное, наиболее изучены и экспериментально и теоретически. Тем не менее, подробных параметрических исследований особенно для сжимаемых сред проведено не много. Решить уравнения гидромеханики аналитически не представляется возможным из-за их нелинейности. Поэтому все
5 теоретические исследования уравнений основываются на предположениях, значительно упрощающих систему уравнений гидромеханики. Точные решения для течений Пуазейля, между двумя параллельными пластинами, и Куэтта, когда одна из пластин неподвижна, а другая движется с постоянной скоростью, получены при существенном допущении о несжимаемости среды [4, 8, 12]. Кроме того, при всех других более сложных точных решениях предполагается, что потоки среды стационарны. Однако преобладающее большинство течений, с которыми приходится иметь дело на практике, не являются идеализированными. Поэтому целесообразно дальнейшее исследование более сложных вариантов течений, в том числе и ламинарных, а особенно - турбулентных. Тем более, что проблема расчетного предсказания характеристик движения имеющего реальный практический интерес далека от решения и чрезвычайно актуальна.
В настоящее время существует три основных подхода' к численному описанию турбулентной конвекции [13]. Традиционный подход, основанный на решении уравнений, возникающих вследствие применения реинольдсова осреднения уравнений Навье-Стокса (Reynolds-Averaged Navier-Stokes, RANS), менее трудоемкий и не требует огромных вычислительных затрат [14- 19]. До сих пор метод RANS остается наиболее распространенным подходом к моделированию турбулентных течений. Вместе с тем, результаты расчетов по методу RANS очень чувствительны к выбору той или иной замыкающей полуэмпирической модели турбулентности, а иногда и просто не способны отразить характерные особенности, присущие термоконвективным течениям. Опыт трехмерных расчетов турбулентной температурной конвекции на основе ряда моделей RANS показал, что свойственная этим моделям генерация высокого уровня турбулентной вязкости препятствует развитию крупномасштабных трехмерных пульсаций, которые в действительности определяют структуру осредненного движения и конвекции в целом [20]. Этот подход не в состоянии обеспечить приемлемую для практики точность описания турбулентных течений при наличии в
потоке обширных отрывных зон. Более того, хотя возможности усовершенствования полуэмпирических моделей, в принципе, еще не исчерпаны, существенный прогресс в этой области едва ли возможен. Это объясняется специфическими физическими особенностями отрывных течений, в частности, наличием в них так называемых организованных (когерентных) вихревых нестационарных структур, параметры которых определяются конкретными геометрическими характеристиками рассматриваемого течения и граничными условиями. Ясно, что это делает построение универсальной полуэмпирической модели турбулентности для расчета отрывных течений исключительно сложной, если вообще разрешимой задачей [21].
Метод моделирования крупных вихрей (Large Eddy Simulation, LES) предполагает аккуратный расчет переноса импульса и тепла лишь крупными, энергетически важными структурами, что позволяет рассчитывать термоконвективные течения при высоких значениях числа Рейнольдса с привлечением сравнительно простых замыкающих моделей [22 - 26]. Особая привлекательность метода LES применительно к расчетам термоконвективных задач объясняется способностью адекватно воспроизводить эволюцию во времени определяющих конвекцию крупномасштабных вихревых структур. Однако моделирование турбулентных течений в присутствии твердых границ на основе метода LES в чистом виде сопровождается требованиями по сеточному разрешению пристенных областей, в которых присутствуют относительно мелкие вихри [27]. Тем не менее, стоит отметить, что стремление преодолеть ограничения RANS и LES привело к появлению гибридного подхода в 1997г. В работе [28] был сформулирован новый подход к моделированию отрывных течений, получившего название метода Моделирования Отсоединенных Вихрей (Detached Eddy Simulation, DES). Этот метод успешно применяется в расчетах сложных отрывных течений в задачах внешней аэродинамики.
7 Среди трех подходов к численному описанию турбулентности все возрастающей привлекательностью обладает метод прямого численного моделирования (Direct Numerical Simulation, DNS) [13]. Однако метод DNS обеспечивает надежность результатов расчетов только при полном разрешении всех составляющих движения. Выполнение данного условия налагает жесткие требования к вычислительным ресурсам, быстро возрастающие при желании продвинуться вверх по числу Рейнольдса. Этот подход наиболее точен и универсален, однако полноценное использование прямого численного моделирования в задачах с геофизическими масштабами даже по самым оптимистичным прогнозам будет возможно лишь по прошествии нескольких десятилетий. Поэтому характерной особенностью течений, исследованных до настоящего времени в рамках DNS, является их пространственная ограниченность (течения в канале, пограничный слой) и сравнительно небольшое число Рейнольдса. Хотя появление многопроцессорных вычислительных систем с возможностью распараллеливания вычислительного процесса уже дает то быстродействие, которого достаточно для проектирования технических устройств, функционирование которых связано с турбулентными потоками.
Кроме высокопроизводительной вычислительной системы при
проведении расчетов с помощью прямого численного моделирования
необходим эффективный численный метод. Существует множество способов
численного решения уравнений Навье-Стокса [29]. Но для решения задач с
очень маленькими масштабами необходимы высокоточные
пространственные методы. С одной стороны, количество узлов, требуемое методами высокого порядка точности обычно в несколько раз меньше, чем для методов первого или второго порядков [30]. При моделировании реальных потоков с мелкими пространственными структурами методами низкого порядка требования к машинной памяти будут огромны, и возникает необходимость в высокоточных схемах. Решение высокой точности достигается с использованием спектральных методов. Эти методы, однако,
8 более трудны в использовании, особенно при наличии сложной геометрии. Основная проблема, которая возникает при применении высокоточных схем, заключается в постановке и обработке граничных условий. Для схем высокого порядка точности вычислительный шаблон становится шире, и выбор численных граничных условий с сохранением точности и устойчивости - не тривиальная задача. С другой стороны, использование схем высокого порядка ведет к увеличению времени обработки каждой точки (узла), особенно если применяется неявная схема по времени. Однако высокий порядок точности по пространству позволяет снизить количество временных шагов, то есть увеличить шаг по времени. Поэтому основными направлениями в рамках DNS являются разработка схем высокого порядка точности и проведение численных экспериментов с использованием качественных математических моделей. Стоит сказать, что методы высокого порядка привлекают в последнее время все больше внимания, так как их теоретические разработки получают широкое применение в разнообразных научных направлениях, включая моделирование глобальных атмосферных явлений, аэродинамику, океанографию, термодинамику, теоретическую химию [31]. Более того, существует специальная конференция International Conference on Spectral Application and High-Order Methods (ICOSAHOM) посвященная проблемам высокоточного моделирования.
Теория турбулентности далека от своего завершения. Появляются все новые подходы к ее изучению. Растет число моделей, предлагаемых для понимания отдельных ее свойств. В работах [32, 33] делается обзор именно таких подходов, которые еще не стали хрестоматийными. Анализ Фурье, теория динамических систем, теория фракталов, вейвлет-анализ - вот далеко не полный перечень областей науки, которые дают новые идеи в исследовании турбулентности.
Далее, проведем обзор работ, ставящих во главу угла высокий порядок дискретизации, так как именно на это делается акцент и в данной работе при прямом численном моделировании гидромеханических потоков.
9 В работе [34] численно исследуется течение около криволинейного цилиндра с числами Рейнольдса 100 и 500. Вычисления проводятся с использованием спектрального /г/?-метода. Данный метод обладает высоким порядком дискретизации, причем одновременно может повышаться число элементов (/г-точность) и порядок полинома разложения в области элемента (р-точность). Метод имеет экспоненциальную сходимость с показателем степени полинома, р. Алгоритм для решения несжимаемых уравнений Навье-Стокса основывается на абсолютно устойчивой схеме расщепления. Нелинейные члены рассчитываются по явной временной схеме, линейные -по неявной. Вычисления на начальной стадии выполняются с использованием полиномов второго порядка, затем порядок повышается до четвертого и, наконец, до шестого. При числе Рейнольдса 100, когда имеет место стационарное течение, разница моделируемых параметров при переходе со второго порядка на четвертый составила приблизительно 4%, а при переходе с четвертого на шестой - менее 0.1%. При числе Рейнольдса 500 авторы получают нестационарное трехмерное течение и методы ниже четвертого порядка не используют.
Похожий алгоритм интегрирования несжимаемых уравнений Навье— Стокса используется в работе [35], в которой моделируется течение за круговым цилиндром. Применяется абсолютно устойчивая схема расщепления по времени высокого порядка, когда линейные и нелинейные члены рассчитываются на разных промежуточных шагах. Дискретизация же по пространству основывается на методе спектральных элементов. Интерполирование внутри элемента осуществляется полиномами Лежандра-Лагранжа восьмого порядка. При интерполировании с десятым порядком разница чисел Струхаля получается в пределах 1%, и авторы делают вывод, что полиномы восьмого порядка дают приемлемое пространственное разрешение.
Спектральный метод коллокации Чебышева используется в работе [36] при численном исследовании двумерной локальной нестационарности,
10 возникающей в течении Пуазейля. В работе [37] несжимаемые уравнения Навье-Стокса решаются с помощью псевдоспектрального метода: в горизонтальных направлениях вычислительная область дискретизируется с помощью рядов Фурье, а в нормальном к стенке направлении используются полиномы Чебышева. В работе решается задача управления с обратной связью возрастающим погранслоем.
Однако спектральные методы имеют ряд недостатков. Как указывалось выше, они достаточно сложны в использовании, особенно в задачах со сложной геометрией и не применимы, когда нестационарность вызвана движением стенки. Поэтому в последнее время значительное внимание уделяется разработке и применению схем, использующих в своей основе методы конечных разностей и конечных объемов.
Например, в работе [38] используется пространственная схема центральных разностей четвертого порядка точности при моделировании сжимаемого вязкого потока около движущейся каверны в совокупности со схемой Рунге-Кутта второго порядка по времени. Показано, что аналогичная схема второго порядка применима только для низких чисел Рейнольдса (до 300). Кроме того, данные схемы четвертого и второго порядка апробируются при решении уравнения Бюргерса. И, если при числе Рейнольдса, равном 100, ошибки дискретизации сопоставимы, то при числе Рейнольдса 200 и более точность схемы четвертого порядка на порядок выше.
Компактная схема конечных разностей четвертого порядка так же используется в работе [39] при решении уравнений Навье-Стокса около движущейся каверны. Авторы приходят к выводу, что использование схем с более низким порядком точности допустимо только в задачах при характерных числах Рейнольдса до 100 даже при достаточном насыщении сетки. С возрастанием числа Рейнольдса ошибка схем низкого порядка оказывается около 10%, тогда как метод четвертого порядка имеет погрешность меньше 0,32%, а в сравнении со спектральным методом шестого порядка [40] - 0,45%).
Сравнение схем второго и четвертого порядка точности при моделировании потока около движущейся каверны в паре с явным методом Рунге-Кутта интегрирования по времени делается и в работе [41]. Авторы делают вывод, что схема четвертого порядка более эффективна.
Стационарно - нестационарный переход двумерного отделившегося пузыря исследуется в работе [42] с помощью прямого численного моделирования. Несжимаемые уравнения Навье-Стокса дискретизируются с помощью компактной конечно-разностной схемы четвертого порядка по пространству и явной схемой Рунге-Кутта третьего порядка по времени. Прямое численное моделирование ламинарно- турбулентного перехода течения в плоском канале осуществляется в работе [43]. Авторы используют полуспектральный метод, когда в продольном и поперечном направлениях применяется разложение Фурье, а в нормальном к стенке направлении -компактная конечно- разностная схема четвертого порядка точности. В работе [44] при численном исследовании пограничного слоя на вогнутой поверхности используется классический метод Рунге-Кутта четвертого порядка при интегрировании по времени и компактная конечно-разностная схема 6 порядка при дискретизации по пространству, за исключением пристенной области, где интегрирование осуществляется с пятым порядком точности.
В работах [45-47] моделируется сжимаемое дозвуковое и сверхзвуковое течение в канале. Уравнения Навье-Стокса записываются в формулировке давление-скорость-энтропия. Компактная противопоточная схема Адамса и Шерифа пятого порядка используется для дискретизации членов уравнения Эйлера, компактная схема Лиля шестого порядка - для молекулярных членов. По времени уравнения интегрируются методом Рунге-Кутта третьего порядка. Кроме того, в работе [45] результаты прямого численного моделирования сравниваются с результатами LES моделирования.
Сравнение результатов DNS и LES моделирования проведено и в работе [48]. Исследовался несжимаемый поток в канале с обратной ступенькой. Для пространственной дискретизации конвективных и диссипативных членов была выбрана компактная схема Эрмита четвертого порядка для неравномерной сдвинутой сетки.
Говоря о схемах высокого порядка точности, возникает вопрос об их эффективности, особенно с точки зрения вычислительных затрат. Именно этой проблеме уделяют внимание авторы работы [49]. Понятия высокая точность и высокий порядок не всегда синонимы. По мнению авторов, главный критерий точности пространственной схемы - возможность минимизировать ошибку для максимальных волновых чисел. Имея в виду данное условие, схемы высокого порядка точности позволяют уменьшить необходимое число расчетных точек, а соответственно и сэкономить машинное время. При решении двумерных задач высокоточные методы позволяют уменьшить объем вычислений на два порядка, для трехмерных расчетов это соотношение еще больше. При моделировании двумерного невязкого потока за цилиндром с использованием метода четвертого порядка точности при интегрировании по пространству, основанном на разложении в ряд Тейлора, авторами показана экономия машинного времени в 98% по сравнению с расчетами, выполненными без использования схем высокого порядка.
Существуют работы, в которых совмещено использование различных методов пространственной дискретизации, например [50]. Авторы используют метод конечных разностей шестого порядка точности для дискретизации в двух направлениях (по потоку и нормальному к стенке), а в третьем направлении (поперечном) используют спектральный метод. Для интегрирования по времени используется схема Рунге- Кутта четвертого порядка точности.
Для эффективной реализации методов прямого численного моделирования турбулентности возникает необходимость в использовании
13 высокопроизводительных вычислительных систем. Наиболее подходящими системами в современных условиях являются многопроцессорные комплексы. Использование таких вычислительных машин требует разработки экономичных алгоритмов, что является самостоятельной достаточно сложной задачей. В работе [51] осуществлена реализация компактных схем высокого порядка (четвертого и выше) на многопроцессорном кластере с использованием технологии МРІ для DNS и LES моделирования турбулентных потоков. Авторы используют фильтры высокого порядка для подавления ложных численных осцилляции, возникающих при использовании компактных разностных схем высокого порядка. В работе проводится анализ эффективности численных методов разного порядка точности при расчете на разном числе процессоров.
Одна из проблем, возникающих при использовании разностных методов высокого порядка, - получение монотонных решений, особенно при наличии разрывов. Для решения данной проблемы существуют методы реконструкции решения. Наиболее известные из них ENO (essentially non-oscillatory) и WENO (weighted essentially non-oscillatory) методы с автоматическим анализом гладкости решения. Они позволяют автоматически достигнуть высокого порядка, не приводя к появлению случайных колебаний в окрестностях разрывов, при решении задач, которые содержат ударные волны и другие сложные структуры течения. Подробный анализ методов и примеры применения содержатся в работе [52]. В работах [53, 54] также показана возможность использования ENO-схемы в сочетании с методом Рунге-Кутта третьего порядка интегрирования по времени.
Прямое численное моделирование турбулентных потоков в канале с внезапным расширением на входе выполнено в работе [55]. Приведен алгоритм расчета с высоким требуемым порядком точности, как по времени, так и по пространству. Проведен анализ сходимости решения при аппроксимации производных и устойчивости схем. Там же в работе содержится обзор методов моделирования гидромеханических процессов,
14 имеющих в своей основе высокий порядок точности дискретизации частных производных.
Научная новизна работы заключается в следующем:
-Реализован метод для прямого численного моделирования гидромеханических процессов, основанный на решении полных уравнений гидромеханики, описывающих трехмерные течения вязкого сжимаемого газа, с помощью устойчивых разностных схем высокого порядка точности; показана высокая точность, хорошая работоспособность предложенных схем в сравнении с имеющимися экспериментальными и расчетными данными.
-Проведено численное моделирование трехмерных ламинарных, переходных и турбулентных течений в канале с обратной ступенькой; исследованы структура и параметры течений в зонах отрыва и присоединения потока; получены осредненные и мгновенные картины течения; проведено сравнение с экспериментальными и теоретическими данными.
-Впервые проведены детальные численные исследования структуры и параметров ламинарных, трехмерных переходных и турбулентных течений в канале с резким расширением на входе; исследованы все основные стадии эволюционирования потока: отрыв и связанное с ним образование рециркуляционных зон, установление течения; вихреобразование и зарождение нестационарности; диссипация и переход к развитому турбулентному течению.
-Впервые исследовано влияние линейных размеров прямоугольного канала с внезапным расширением и различных граничных условий на характеристики потока.
-Получены статистические характеристики крупномасштабной турбулентности в ядре потока, найдены распределения пульсационных характеристик скорости по спектрам.
Достоверность результатов подтверждается следующим:
-Построенная математическая модель основывается на системе полных уравнений гидромеханики и базируется на фундаментальных законах механики сплошной среды.
-Разработанные численные алгоритмы апробированы при решении тестовых задач и показывают высокую точность в широком диапазоне варьируемых параметров.
-Полученные численные результаты согласуются с известными аналитическими, экспериментальными и расчетными данными.
Личный вклад автора состоит в разработке математической модели и создании алгоритма для проведения численного исследования гидромеханических процессов. Проведено сравнение полученных численно результатов с известными расчетными и экспериментальными данными. Автором исследованы симметричные [88 - 90] и асимметричные [91, 92] ламинарные стационарные потокрі газа, ламинарные нестационарные и переходные течения [93 — 95]. Все указанные исследования выполнены на основе анализа численных результатов, полученных лично автором. Анализ выполнен совместно с академиком A.M. Липановым.
На защиту выносится:
-Алгоритм метода высокого порядка точности по времени и пространству для расчета течений вязкого сжимаемого газа в плоском канале с внезапным расширением на входе.
-Результаты тестовых расчетов с использованием метода высокого порядка точности интегрирования уравнений гидромеханики при решении задачи о течении газа в канале с обратной ступенькой
-Результаты спектрального анализа колебаний компоненты вектора скорости во времени.
-Результаты прямого численного моделирования течения в плоском канале с расширением. Влияние характерного числа Рейнольдса на
закономерности течения. Исследование влияния линейных размеров канала на характер течения и изменение гидромеханических параметров.
По главам содержание работы распределено следующим образом. В первой главе описывается математическая модель течения вязкого сжимаемого газа в плоском канале с расширением на входе, приводится система уравнений гидромеханики, определяются начальные и граничные условия.
Во второй главе описывается численный метод решения системы уравнений гидромеханики, алгоритм интегрирования по времени, конечно-разностная схема применительно к внутренним и граничным точкам, реализация граничных условий.
В третьей главе анализируются выполненные методические расчеты, приводятся результаты численного моделирования течения газа в канале с обратной ступенькой. Проводится сравнение результатов с известными экспериментальными и расчетными данными. Приводятся расчеты турбулентного течения газа в канале с внезапным расширением на входе. Спектральное распределение энергии продольной компоненты вектора скорости сравнивается с теоретической кривой.
Четвертая глава посвящена детальному исследованию течения в канале с внезапным расширением. Исследуется зависимость от числа Рейнольдса режима течения. Показывается влияние линейных размеров расчетной области (канала) на параметры потока. Строятся распределения пульсаций скорости по спектрам.
Результаты докладывались на IX Всероссийском съезде по теоретической и прикладной механике, г. Нижний Новгород, 2006 г., на III научно-практической конференции «Проблемы механики и материаловедения», г. Ижевск, 2006 г., на конференции молодых ученых «Численные методы в математике и механике», г. Ижевск, 2007 г., на международной конференции «XVIII сессия Международной школы по
17 моделям механики сплошной среды», г. Саратов, 2007 г., на II Всероссийской научно-технической конференции «Безопасность критичных инфраструктур и территорий», г. Екатеринбург, 2008 г.
Основные результаты опубликованы в работах [88-95].
Работа выполнена при финансовой поддержке в рамках гранта молодых ученых и аспирантов УрО РАН и стипендии президента Удмуртской Республики 2006-2007 гг.
Граничные и начальные условия
Однако спектральные методы имеют ряд недостатков. Как указывалось выше, они достаточно сложны в использовании, особенно в задачах со сложной геометрией и не применимы, когда нестационарность вызвана движением стенки. Поэтому в последнее время значительное внимание уделяется разработке и применению схем, использующих в своей основе методы конечных разностей и конечных объемов.
Например, в работе [38] используется пространственная схема центральных разностей четвертого порядка точности при моделировании сжимаемого вязкого потока около движущейся каверны в совокупности со схемой Рунге-Кутта второго порядка по времени. Показано, что аналогичная схема второго порядка применима только для низких чисел Рейнольдса (до 300). Кроме того, данные схемы четвертого и второго порядка апробируются при решении уравнения Бюргерса. И, если при числе Рейнольдса, равном 100, ошибки дискретизации сопоставимы, то при числе Рейнольдса 200 и более точность схемы четвертого порядка на порядок выше.
Компактная схема конечных разностей четвертого порядка так же используется в работе [39] при решении уравнений Навье-Стокса около движущейся каверны. Авторы приходят к выводу, что использование схем с более низким порядком точности допустимо только в задачах при характерных числах Рейнольдса до 100 даже при достаточном насыщении сетки. С возрастанием числа Рейнольдса ошибка схем низкого порядка оказывается около 10%, тогда как метод четвертого порядка имеет погрешность меньше 0,32%, а в сравнении со спектральным методом шестого порядка [40] - 0,45%). Сравнение схем второго и четвертого порядка точности при моделировании потока около движущейся каверны в паре с явным методом Рунге-Кутта интегрирования по времени делается и в работе [41]. Авторы делают вывод, что схема четвертого порядка более эффективна.
Стационарно - нестационарный переход двумерного отделившегося пузыря исследуется в работе [42] с помощью прямого численного моделирования. Несжимаемые уравнения Навье-Стокса дискретизируются с помощью компактной конечно-разностной схемы четвертого порядка по пространству и явной схемой Рунге-Кутта третьего порядка по времени. Прямое численное моделирование ламинарно- турбулентного перехода течения в плоском канале осуществляется в работе [43]. Авторы используют полуспектральный метод, когда в продольном и поперечном направлениях применяется разложение Фурье, а в нормальном к стенке направлении -компактная конечно- разностная схема четвертого порядка точности. В работе [44] при численном исследовании пограничного слоя на вогнутой поверхности используется классический метод Рунге-Кутта четвертого порядка при интегрировании по времени и компактная конечно-разностная схема 6 порядка при дискретизации по пространству, за исключением пристенной области, где интегрирование осуществляется с пятым порядком точности.
В работах [45-47] моделируется сжимаемое дозвуковое и сверхзвуковое течение в канале. Уравнения Навье-Стокса записываются в формулировке давление-скорость-энтропия. Компактная противопоточная схема Адамса и Шерифа пятого порядка используется для дискретизации членов уравнения Эйлера, компактная схема Лиля шестого порядка - для молекулярных членов. По времени уравнения интегрируются методом Рунге-Кутта третьего порядка. Кроме того, в работе [45] результаты прямого численного моделирования сравниваются с результатами LES моделирования. Сравнение результатов DNS и LES моделирования проведено и в работе [48]. Исследовался несжимаемый поток в канале с обратной ступенькой. Для пространственной дискретизации конвективных и диссипативных членов была выбрана компактная схема Эрмита четвертого порядка для неравномерной сдвинутой сетки.
Говоря о схемах высокого порядка точности, возникает вопрос об их эффективности, особенно с точки зрения вычислительных затрат. Именно этой проблеме уделяют внимание авторы работы [49]. Понятия высокая точность и высокий порядок не всегда синонимы. По мнению авторов, главный критерий точности пространственной схемы - возможность минимизировать ошибку для максимальных волновых чисел. Имея в виду данное условие, схемы высокого порядка точности позволяют уменьшить необходимое число расчетных точек, а соответственно и сэкономить машинное время. При решении двумерных задач высокоточные методы позволяют уменьшить объем вычислений на два порядка, для трехмерных расчетов это соотношение еще больше. При моделировании двумерного невязкого потока за цилиндром с использованием метода четвертого порядка точности при интегрировании по пространству, основанном на разложении в ряд Тейлора, авторами показана экономия машинного времени в 98% по сравнению с расчетами, выполненными без использования схем высокого порядка.
Существуют работы, в которых совмещено использование различных методов пространственной дискретизации, например [50]. Авторы используют метод конечных разностей шестого порядка точности для дискретизации в двух направлениях (по потоку и нормальному к стенке), а в третьем направлении (поперечном) используют спектральный метод. Для интегрирования по времени используется схема Рунге- Кутта четвертого порядка точности.
Метод интегрирования уравнений гидромеханики по пространственным координатам
Нестационарность и нелинейность гидромеханических процессов, а также интенсивное изменение гидромеханических параметров в турбулентных потоках делают необходимым их рассмотрение в пространственно-временном континууме.
При численном решении систем дифференциальных уравнений с частными производными вся область интегрирования покрывается достаточно густой ортогональной пространственной сеткой. В общем случае шаг Ах может быть функцией номера ячейки, так же как и шаги Ау и Az. Однако в дальнейшем в работе и при проведении численных расчетов считаем, что Ах = const, Ay = const, Az = const, хотя шаги в разных направлениях могут быть и не равны. Шаг интегрирования по времени обозначим At и определим его разностью At = tn+1n, где п - номер временного слоя.
В условиях интенсивного изменения ГМП использование разностных схем, применяемых для расчета плавно изменяющихся параметров, когда достаточно первого и второго порядков точности при аппроксимации частных производных по пространству, будет грубым. Как указывалось выше, такие методы обладают значительной схемной вязкостью, что целиком искажает реальную картину течения, особенно когда дело касается турбулентных потоков с мелкими пространственно-временными масштабами. Отсюда возникает необходимость в разработке метода высокого порядка аппроксимации как по времени, так и по пространству. Выразим величину ГМП W на (и + 1)-м временном слое в точке (i,j,k) через величины параметров на п -м слое в той же точке [56] Аґ - малый параметр.
Ряд Тейлора (2.1) позволяет отделить расчет параметров во времени от их расчета вдоль пространственных координат при переходе с одного временного слоя на другой. Причем этот переход возможен с любой требуемой точностью. Также стоит отметить, что с помощью ряда (2.1) достигается разделение изменения параметра W во времени от его изменения в пространстве лишь на малом временном отрезке At. Такой подход удобнее ввиду очень сложной зависимости ГМП от времени и от пространственных координат. Разделение пространственной и временной зависимостей на маленьком временном отрезке позволяет достичь сходимости по времени, и даже использование всего двух членов ряда (2.1) делает переход с одного временного слоя на другой по сути нелинейным, в отличие от методов Рунге-Кутта, когда на промежуточных итерациях (относящихся к дробным частям шага At) переход происходит по прямой (линейно). Это особенно важно ввиду того, что нерегулярное и хаотическое поведение потока, особенно при наличии развитой турбулентности связано со свойством нелинейных систем экспоненциально быстро разводить решения в ограниченной области пространства [33].
С помощью ряда (2.1) реализуется асимптотическая сходимость при интегрировании по времени уравнений (1.8)-(1.12). При этом с заданной точностью минимизируется влияние схемной вязкости. Кроме того, необходимо учитывать, что удовлетворяющий условию устойчивости шаг интегрирования по явной временной схеме At столь мал, что даже четырех членов ряда (2.1) достаточно, чтобы ошибка между численным результатом и решением уравнений была за пределами последней значащей цифры числа в машинном представлении.
Точки отрыва и присоединения ламинарных и турбулентных течений за обратным уступом
В соответствии с изложенным во второй главе содержанием приведем некоторые результаты в обоснование процессов сходимости при интегрировании уравнений гидромеханики. Рассмотрим в данном разделе течение в канале с резким расширением на входе, которое имеет ряд особенностей. Наличие двух уступов приводит к отрыву и развитию за их угловыми точками двух сдвиговых слоев смешения, сложным образом взаимодействующих между собой и со стенками канала. Такие сдвиговые слои при определенных режимах течения содержат крупномасштабные вихревые структуры, эволюционирующие по потоку. Такие эффекты течений в рассматриваемом диапазоне чисел Рейнольдса и Маха (Re 103-104, М = 0.6) объясняются «резонансными» условиями вихреобразования [55].
В работах [55, 56] показана сходимость процесса интегрировании по времени. Учитывая временной шаг интегрирования, который во всех расчетах не превышал 0.002, такая сходимость становится достаточно очевидной, так как значение At2 лежит на границе последних значащих цифр в машинном представлении числа с плавающей точкой. При более мелком шаге значение At3 вообще выходит за эту границу.
В тех же работах [55, 56] и в работах [52, 57-59] показано, что восьмого порядка точности при интегрировании по пространственным переменным достаточно для получения достоверных результатов на приемлемой для расчета сетке (имеется в виду общее количество узлов интегрирования в сопоставлении с машинными мощностями).
Не смотря на высокий порядок аппроксимирующего выражения при интегрировании уравнений гидромеханики по пространственным переменным, если число элементарных объемов (узлов расчетной сетки) недостаточное, пространственное сглаживание решения будет большим. Особенно важно то, что не будет достигнута требуемая точность получаемого решения. Далее, проведем анализ выбора оптимального значения числа элементарных объемов (ЭО) в области интегрирования (ОИ). Проведем расчет течения с характерными числами Re = 5000 и М = 0.6 на сетках с разным числом узлов. На рис. 3.1 приведены поля векторов скорости движения газа в плоскости симметрии для различных вариантов измельчения ОИ на ЭО. Объем интегрирования в данном случае представлял собой трехмерный канал длиной = 12 с квадратным поперечным сечением со сторонами Н = 2 Hz = 3 (все величины указаны в безразмерных единицах). На рис. 3.2 показано распределение продольной компоненты вектора скорости на оси симметрии X (Y = 1.5, Z = 0 ).
Видно, что с ростом числа ЭО картина течения детализируется, становится более подробной. При 89725 ЭО намечаются лишь грубые контуры траекторий движения газа. Такую же оценку можно сделать и при 208593 ЭО. Кривые U(x) в случаях а и б на рис. 3.2 более сглажены и не воспроизводят многих колебаний, которые появляются при переходе к 690361 и 5415025 элементарным объемам. В сущности, только начиная со случая в можно анализировать процесс асимптотической сходимости картины течения для плоскости симметрии ОИ и в целом. А, рассмотрев рис. 3.2, можно говорить о достаточной близости полученных кривых U(x) только в случаях виг. Для количественной оценки результатов рассмотрим вектор завихренности Q = rot V, V = (f/,F, ), определяемый проекциями [3, 64]:
Как видно из данных таблицы переход с 690361 ЭО к 5415025 ЭО позволил уменьшить разницу Латах до 1%, таким образом, с приемлемой точностью обеспечить сходимость приближенного решения к точному при квантовании области интегрирования на элементарные объемы. Векторы скорости в поперечном сечении (YZ) и векторы вихря - в продольном (XY) для двух последних вариантов сетки показаны соответственно на рис. 3.3 и рис.3.4.
Таким образом, можно сделать вывод, что при интегрировании по пространственным переменным уравнений гидромеханики с восьмым порядком точности при Re 103-104 сетки, пропорциональной третьему варианту, достаточно для получения достоверных результатов, хотя, в некоторых задачах (требующих получения мгновенных характеристик) может потребоваться увеличение мощности сетки. Очевидно, что при расчете с более низкими числами Re ламинарных и переходных течений, сходимость будет, по крайней мере, не хуже, а даже на порядок лучше. Поэтому, третий вариант расчетной сетки будем считать удовлетворительным.
Пригодность математической модели (метода моделирования) движения сплошной среды сильно зависит от точности предсказания явлений отрыва и присоединения. Контроль и регулирование отрывных течений -одни из основных задач механики сплошных сред. В данном разделе рассматривается классический случай отрывного течения в канале с обратной ступенькой (backward-facing step). Результаты численного моделирования сравниваются с известными результатами физических экспериментов. Геометрия канала показана на рис. 3.5.
На рис. 3.5,а показано пространственное расположение области интегрирования и координатных осей. Рис. 3.5,6 (показана проекция канала на плоскость XY) соответствует ламинарному режиму течения с низкими числами Re, когда на входе в канал задается профиль скорости Пуазейля с учетом формулы (1.16). Число Re определяется по формуле
Несимметричные ламинарные стационарные течения
Для моделирования течений в плоском канале со скачком площади поперечного сечения на входе была составлена параллельная программа с использованием технологии MPI. С помощью программы были проведены расчеты в канале в широком диапазоне чисел Рейнольдса. Расчеты проводились на многопроцессорных вычислительных системах с использованием ресурсов МСЦ РАН. Область интегрирования однонаправлено делилась на равные подобласти вдоль оси X.
В ламинарной области изменения чисел Рейнольдса интегрирование уравнений гидромеханики по времени осуществлялось методом Рунге-Кутта со вторым порядком точности. Если же число Рейнольдса соответствовало турбулентному или переходному диапазону изменения, то интегрирование по времени осуществлялось методом, описанным в главе 2, с использованием формулы (2.8).
Интегрирование по пространственным переменным выполнялось с восьмым порядком точности. Внутри области интегрирования производные аппроксимировались центральными разностями восьмого порядка на девятиточечном шаблоне. В окрестности границ производные вычислялись из внутренних точек методом конечных разностей, однако, по причине несимметричности шаблона в этом случае он был расширен до десяти точек для сохранения восьмого порядка точности при интегрировании по пространству. Диффузионные члены, содержащие производные второго порядка, рассчитывались с применением формулы для вторых производных.
Геометрия канала показана на рис. 1.1. Рассматривались два вида каналов: с боковыми стенками и без них. Величина высоты входного отверстия и величины ступенек равны между собой и принимались равными h, в безразмерном виде - 1. Данная величина h использовалась как масштаб для линейных величин. Ширину рассматриваемого канала, как и его длину с высотой, будем выражать в долях h. В зависимости от решаемой задачи длину канала варьировали (от 10 до 60), для течений с более низкими числами Рейнольдса она была меньше. Число Маха в расчетах принималось равным М = 0.6, число Прандтля Рг = 0.7, отношение теплоємкостей = 1.4. Стационарные ламинарные потоки по сравнению с нестационарными ламинарными и турбулентными потоками, наверное, наиболее изучены и экспериментально, и теоретически [4]. Тем не менее, подробных параметрических исследований особенно для сжимаемых сред не проводилось. Симметричные потоки реализуются при малых величинах чисел Рейнольдса Re. В работах [88-90] симметричный поток получен при числе Re = 50, но уже при Re = 150 поток, хотя и остается стационарным, но оказывается несимметричным. Проследим за процессами установления и изменения по длине канала величин параметров ламинарного потока, когда реализуется симметричная струя в двумерной области интегрирования. Расчеты начнем выполнять с числа Рейнольдса Re = 40, постепенно его повышая. Длину канала в данном случае возьмем равной 10. Внесем в поток трехмерное возмущение в течение 0.2 единиц безразмерного времени (100 итераций) согласно формуле (1.20), как описано в главе 1. Здесь и далее будем следить за изменением третьей компоненты W вектора скорости потока. Об эволюции W можно судить, анализируя рис.4.1, на котором показана зависимость W от х и у в плоскости симметрии (Z = 0) в различные моменты времени. На рис. 4.1,а переменная W соответствует моменту, когда на входе в канал ее значение оказывается равным нулю не только на границе, но и в ближайшей внутренней точке (реализуется условие dW/dx = 0 при х = 0). Это условие реализуется в момент времени / = 3.4 и соответствует числу Рейнольдса, равному 50. С ростом времени третья компонента W вектора скорости эволюционирует, и, хотя её пространственная конфигурация сохраняется, крутизна поверхности W{x,y) и ее максимальная величина только уменьшаются. Возмущение распространяется вниз по потоку. К моменту времени / = 8 передний фронт возмущения достигает примерно 5.5 (рис.4.1,6). При этом максимум величины W уменьшается более чем в 3 раза, по сравнению с её величиной при / = 3.4 (рис. 4.1,а). При дальнейшем увеличении времени поверхность W[x,y) становится более пологой с максимумом около её переднего фронта (рис. 4.1,6, / = 14.2). Она по-прежнему имеет две возвышенности в продольном направлении канала, заканчивающиеся общим фронтом, но в целом её максимальная величина уже на порядок меньше. И, наконец, для t 20 третья компонента вектора скорости исчезает. Уже при выходе на правую границу ОИ (рис. 4.1,в) величина W не превышает нескольких единиц в 4-ой значащей цифре и ею можно пренебречь. В вязком симметричном потоке при числе Re = 50 3-я - компонента со временем постоянно уменьшается и для / 14.2 её влиянием. можно пренебречь. А при t 20 поток становится двумерным. С ростом числа Re время, при котором исчезает 3-я компонента W вектора скорости на входе и во всем канале, растет, причем в обоих случаях почти линейно. Так, если для числа Re = 50, как было сказано выше, условие dWJdx = 0 при х = 0 начинает выполняться при / 3.4, то для Re = 40 начиная с t - 2.8, для Re = 60 - с t = 3.9, для Re = 80 - с t = 4.8 и для Re = 90 -с ( = 5.6.
