Численный метод расчета течений сжимаемого вязкого газа в широком диапазоне чисел Маха Чирков Денис Владимирович

Вычисление конвективных разностных потоков

Неявные методы решения пространственных задач можно дифференцировать также по способу обращения стабилизирующего оператора. Метод попеременных направлений (ADI), использованный в работах [3, 56], сводится к последовательности векторных прогонок по каждому из пространственных направлений. К недостатку метода попеременных направлений следует отнести его условную устойчивость в трехмерном случае. В работах [45, 46, 47] использованы варианты метода LU-факторизации неявного оператора, В [4] используется релаксация Гаусса — Зейделя, В работе [61] на примере задачи плоского обтекания головной части затупленного тела сверхзвуковым потоком проведено сравнение различных методов факторизации для неявных схем и показано, что LU-факторизация обеспечивает наиболее оптимальную скорость сходимости.

В работах [46, 47, б, 4, 9] использованы в качестве зависимых консервативные переменные Q = (р, ри, pv, pw, е)т. Такой выбор оказывается удобен при решении уравнений Эйлера. Однако, в случае решения уравнений Навье-Стокса более плодотворным является использование примитивных переменных, например Q = [р, и, v, w, Т)т [49,22, 30]. Во-первых, такой подход позволяет наиболее полно учесть якобианы вязких членов в неявном операторе, что увеличивает запас устойчивости метода. Во-вторых, вектор примитивных переменных легко допускает переход к относительному давлению рд = р—ро при малых числах Маха, что позволяет устранить ошибки машинного округления при аппроксимации градиента давления [49, 39].

Актуальной проблемой при разработке численных методов решения уравнений Навье — Стокса является верификация метода, в частности сравнение с экспериментальными данными. Особенно это относится к трехмерным программам. Возникает потребность собрать банк достоверных экспериментальных данных и результатов численных расчетов задач вязкого ламинарного обтекания тел простой формы, для которых, с одной стороны, легко строится сетка, а с другой стороны проявляются основные особенности вязких сжимаемых течений. К этим особенностям можно отнести ударные волны, тонкие пограничные слои, отрыв и присоединение потока, зоны возвратного течения.

В последнее время одним из основных инструментов проектирования турбинных установок является численное моделирование. При моделировании течения воды в элементах проточного тракта гидротурбин для решения уравнений Эйлера и Навье-Стокса несжимаемой жидкости зачастую используется метод искусственной сжимаемости [62, 63]. При этом для аппроксимации уравнений в основном применяются схемы с направленными разностями с расщеплением матрицы Якоби по знакам ее собственных значений, как, например, в [64, 65]. Численные эксперименты показывают, что точность подобных схем существенным образом зависит от параметра псевдосжимаемости /?ПС который, как правило, берется постоянным для всей расчетной области. С точки зрения сходимости и точности расчета основного потока оптимальное значение параметра предобуславливания Д,с = ifv ,p, где v p — характерная скорость потока, а К = 5 -г 10. Однако в случае, когда расчетная область содержит обширные подобласти стагнации потока, где jv g; [vjxap, точность метода в этих подобластях падает. Поэтому до сих пор актуальной остается проблема исследования точности алгоритмов типа [64, 65], основанных на методе искусственной сжимаемости. Один из подходов, позволяющий оценить точность подобных алгоритмов расчета течений капельной жидкости состоит в сравнения полученных с их помощью решений с результатами расчетов, проведенных в рамках модели сжимаемой жидкости при М-С1.

Таким образом, сохраняется потребность в экономичном и универсальном численном методе, позволяющем рассчитывать движения сжимаемого вязкого газа в реальных аэрогидродинамических установках в широком диапазоне характерного числа Маха потока. Идея создания подобного алгоритма состоит в построении эффективного метода решения уравнений Навье—Стокса при умеренных и больших числах Маха, и дальнейшем распространении области его применимости на течения сМ 1 при помощи методики предобуславливания. Диссертационная работа посвящена: 1) Разработке эффективного численного метода решения уравнений Навье-Стокса сжимаемого газа, обладающего высокой точностью разрешения основных особенностей газодинамических течений, таких как ударные волны, тонкие пограничные слои, явления отрыва и присоединения потока, а также имеющего абсолютную устойчивость и приемлемую сходимость. Этот алгоритм берется в качестве базового и в дальнейшем распространяется на случай расчета течений с малыми числами Маха. 2) Построению и исследованию алгоритма предобуславливания, который позволит существенно улучшить сходимость метода установления при М g; 1 и обеспечит сохранение точности базового алгоритма при малых числах Маха. 3) Созданию комплекса программ для моделирования течений вязкого газа, основанного на разработанном численном алгоритме. 4) Приложению комплекса к ряду задач, связанных с проектированием рабочих колес турбомашин. Настоящая работа продолжает цикл работ [64, 66, 67, 68], посвященных моделированию течений вязкой несжимаемой жидкости, и распространяет предложенные в них удачные идеи на случай расчета сжимаемых течений. Диссертационная работа состоит из введения, четырех глав, заключения, приложения и списка литературы. В главе 1 строится базовый численный метод CoNS3D решения системы уравнений Навье — Стокса/Эйлера динамики сжимаемого газа на структурированных сетках с ячейками в виде криволинейных параллелепипедов. Построение численного алгоритма ведется в рамках идеологии метода конечных объемов. Для аппроксимации конвективных членов в уравнениях движения используется противопотоковая TVD-схема Чакраварти-Ошера [44] с ограничителями Ван-Лира [57]. На гладких решениях метод обладает вторым порядком аппроксимации по пространству. Реализация метода сводится к вычислению невязких и вязких разностных потоков через грани ячейки на n-м слое, и обращению бегущим счетом LU-факторизованного стабилизирующего оператора. В стабилизирующем операторе учтены как якобианы невязких потоков, так и якобианы вязких потоков. Дано обоснование выбора в качестве зависимых переменных вектора примитивных переменных Q = (рд,іл, v, w,T)T. Подробно изложены линеаризация вязких потоков и способы реализации краевых условий на твердых стенках, входных, выходных и удаленных границах.

Взаимодействие ударной волны с погранслоем на пластине

Эта тестовая задача позволяет оценить точность расчета особенностей течения, связанных с отрывом и присоединением ламинарного пограничного слоя.

Косой скачок уплотнения, взаимодействующий с пограничным слоем, приводит к отрыву последнего и формированию в месте взаимодействия зоны возвратного течения (рис. 7). Угол наклона косого скачка уплотнения в = 32.6. Расчетная область представляет собой прямоугольник размера 1.6L х 1.1L. Здесь L = 0.049 м — расстояние по оси х от переднего края пластины до точки падения скачка уплотнения, взятое в качестве характерной длины для определения числа Рейнольдеа Re = 2.96 10s. Число Маха набегающего потока М = 2.0. Расчет проводился на сетке 90 х 168, которая сгущалась к поверхности пластины. Высота первой ячейки сетки составляла 0.0002L. На пластине ставилось условие адиабатичности дТ/дп = 0 для температуры и условие прилипания для скорости, реализованное по формуле (1.34). Значение Т для этой задачи бралось равным 167 К, что соответствует температуре торможения потока Го = 300 К.

Расчеты задачи на последовательности сеток 60 х112, 90 х168,120 х 224 и 180 х 336 показали, что решение, полученное на сетке 90 х 168, является предельным и не меняется при дальнейшем сгущении сетки. Результаты численных расчетов сравнивались с экспериментальными данными из работы Хаккинена и соавторов [78] и расчетами Томаса — Уолтерса [4]. Распределение давления р (нормированного на р ) и коэффициента трения су на поверхности пластины представлено на рис. 8, а, б. Коэффициент трения вычислялся по формуле (2.2).

Видно, что давление с хорошей точностью совпадает как с расчетом по схеме Томаса — Уолтерса, так и с экспериментальными данными. Черными кружками на рис. 8, б обозначены точки, в которых экспериментально зафиксировано отрывное течение, но значения коэффициента трения с/ получены не были. Для предложенного метода можно отметить увеличение длины отрывной рециркуляционной зоны (рис. 8, б, область, где с/ 0) по сравнению с данными Томаса — Уолтерса и экспериментом, а также заниженную величину с/ в задней части рециркуляционной зоны (x/L а 1.1). Очевидно общее для двух численных методов занижение с/ за отраженным скачком [x/L 1.2). В целом результаты расчета удовлетворительно согласуются с экспериментом.

Структура течения в угле сжатия 9 — 24 схематично представлена на рис. 9. Угол сжатия состоит из горизонтальной пластины и следующего за ней клина. В случае достаточно больших углов клина (9 18е) при его обтекании образуется зона возвратного течения в окрестности вершины угла. Входные данные для этой задачи (см. таблицу 1) соответствуют эксперименту Холдена и Мозеллы [79], с которым и проводилось сравнение. В большинстве работ, посвященных численному исследованию этой задачи, она решается в двумерной постановке. При этом в численных расчетах [80] даже на очень подробных сетках размеры рециркуляционной зоны по крайней мере в полтора раза превосходят полученные в эксперименте. В [81] показано, что конечная ширина модели, использованной в эксперименте, играет существенную роль. Действительно, ширина установки (ширина пластины и клина) в эксперименте d= 0.6096 м сравнима с длиной горизонтальной части хс = 0.439 м. Поэтому в настоящей работе наряду с двумерными расчетами на сетке 100 х 100 представлены также результаты расчета этой задачи в трехмерной постановке, проведенные на сетке 100 х 100 х 25. Расчетная область для моделирования трехмерного обтекания показана на рис. 10 штриховой линией. На рисунке твердая поверхность угла сжатия ограничена сплошной линией. Разбиение области в поперечном направлении выбиралось так, чтобы на поверхность модели между плоскостью симметрии и боковой стороной (соответствующей точке F) приходились 19 ячеек сетки и б ячеек — на область течения между боковой стороной модели и внешним потоком (отрезок FA). Сетка сгущалась к твердой поверхности в нормальном направлении, и к переднему краю пластины и вершине угла в направлении оси х. Для реализации условий на границах с внешним потоком, соответствующих граням ABODE, AFGE, EGHD, использовались одномерные инварианты Римана (1.41).

Рассмотрим для начала результаты расчета задачи в плоском (двумерном) приближении. На рис. 11 показано распределение коэффициента давления Ср = 2р/р0Ои о вдоль поверхности, полученное экспериментально, настоящим методом и с помощью программы CFL3D [3] (данные взяты из [81]). Видно, что при решении задачи в плоском приближении точка отрыва погранслоя хв, которой соответствует скачок давления в окрестности х/хе = 0.5, для обоих численных методов предсказывается неудовлетворительно, а именно, смещен влево на 0.25гс. Далее, максимум в распределении коэффициента давления (соответствующий точке взаимодействия ударной волны от передней кромки с поверхностью клина) сильно смещен вправо относительно эксперимента. На рис. 12 представлено соответствующее распределение коэффициента теплопередачи с , определяемого по формуле в сечении z = const, проходящем по центру ячеек, примыкающих к плоскости симметрии. Результаты трехмерных расчетов по настоящему методу и CFL3D гораздо лучше совпадают с экспериментальными данными, о чем свидетельствуют кривые, представленные на рис. 11 и 12. Изолинии давления и линии тока на поверхности угла сжатия приведены на рис. 10. Видна тенденция к сужению рециркуляционной зоны с приближением к краю пластины.

Таким образом, еще раз подтверждается вывод работы [81] о том, что для проведения корректного сравнения с экспериментом Холдена — Мозеллы необходимо учитывать пространственную структуру течения в угле сжатия. При адекватной постановке задачи предложенный численный метод, равно как и CFL3D, с хорошей точностью передает особенности течения.

Обзор работ и методология предобуславливания

Для предотвращения проблем со сходимостью при М 1 в середине 80-х годов был предложен метод предобуславливания (preconditioning). Идея предобуславливания при М —У 0 состоит в модификации члена с производной по времени в исходных уравнениях движения путем домножения его на матрицу Р-1: которая подбирается так, чтобы скорости распространения возмущений (суть — собственные числа матриц РА и РВ, где A = 9E/9Q и В = 9F/0Q — матрицы Яко-би) имели один порядок при М — 0. При этом на дифференциальном уровне при установлении (9Q/9i = 0) система (3.1) описывает стационарное решение исходных уравнений Эйлера. Этот прием аналогичен методу искусственной сжимаемости, предложенному в [62, 63] для решения уравнений динамики несжимаемой жидкости. Метод предобуславливания применительно к уравнениям газовой динамики при малых числах Маха впервые использован в работе Чоя и Меркла [20]. Предложенная ими матрица предобуславливания для системы, записанной в переменных q = (р, щ v,p)T имеет простой диагональный вид Р"1 = diag(l, 1,1, М-2), а при переходе к консервативным переменным превращается в имеют один порядок при М 1, то есть модифицированные уравнения (3.1) уже не являются жесткими. В [20] использование такого предобуславливания позволило получить для задачи о течении в решетке профилей не зависящую от числа Маха Моо сходимость метода установления при М 0.05.

В 1987 г. Туркел опубликовал теоретическую работу [28], посвященную методам предобуславливания для уравнений несжимаемой и сжимаемой жидкости. В ней автор распространил метод искусственной сжимаемости Чорина [63] на случай сжимаемого газа, получив двухпараметрическое семейство матриц предобуславливания, отличных от найденной в [20], но так же обеспечивающих одинаковый порядок собственных значений модифицированной матрицы Якоби РА при всех М 1.

В это же время Стрелец и Шур [49] успешно применили «метод масштабирования сжимаемости» (по сути — метод предобуславливания) для расчета вязких течений сжимаемого газа при М ; 1. Предложенный здесь способ модификации нестационарного члена уравнений выглядит достаточно искусственным, не нашедшим должного обоснования и исследования. Однако, к особенностям этой работы следует отнести использование в качестве одной из зависимых переменных относительного давления j/ = р — Ро, где ро — характерное давление в потоке), что весьма плодотворно при расчете течений с Мхар- ; 1, а также гибридный, соединяющий в себе направленные и центральные разности, способ аппроксимации производных конвективных членов и членов с давлением. При этом для направленной аппроксимации производных конвективных членов в уравнениях сохранения импульса и члена типа div v в уравнении неразрывности выбиралось направление против потока, а для аппроксимации членов типа gradj/ — по потоку. То есть для построения противопотоковой схемы использовалось расщепление вектора потока по давлению [89, 48]: Использование подобной аппроксимации позволило в [49] построить численный алгоритм, работающий как при М 1, так и при М ; 1, включая предельный случай М = 0.

Чой и Мёркл в [22] показали, что предложенная ими в [20] матрица предобу-славливания не пригодна для расчета вязких течений, а именно, дает неустойчивый даже в линейном приближении алгоритм. В этой же работе авторы разработали матрицу предобуславливания для уравнений Навье — Стокса и продемонстрировали ее эффективность на задачах обтекания цилиндра, течения в каверне с движущейся крышкой и течения в каверне, индуцированного горизонтальным перепадом температур при числах Маха вплоть до М = Ю-5. Необходимо отметить, что в работе [22], так же как и в [20, 49] введение предобуславливания в численный алгоритм требовало модификации только члена с djdt.

В первых работах, посвященных методу предобуславливания [20, 49, 28, 22], он рассматривался исключительно как способ ускорения сходимости к стационарному решению схем установления. Как правило, эти схемы были центральноразностными, или, как в случае со схемой Стрельца—Шура [49], схемами со специальной аппроксимацией. Вероятно, поэтому, а также потому, что расчеты проводились в основном для течений с М 0.01, серьезных проблем с точностью полученных решений авторы не встретили.

Идея модификации пространственной аппроксимации противопотоковых схем, позволяющая устранить обнаруженное в [23] ухудшение точности при М впервые была предложена в работе Ван-Лира, Ли и Роу [29] применительно к распространенной схеме Роу [50]. Для восстановления точности противопотоковых схем при малых числах Маха авторы работы [29] предложили вычислять член AAj-ny2Q, отвечающий в схеме Роу за схемную вязкость, в соответствии с собственными значениями АІ и диагональным разложением матрицы Якоби РА предобусловленной системы

Вязкое течение в изогнутом канале квадратного сечения

Следующая тестовая задача позволяет выяснить работоспособность метода предобу-славлнвания при расчете вязких течений. Рассмотрим обтекание цилиндра потоком вязкого слабосжимаемого (М К 1) газа с числом Рейнольдса Re = 40. Для этой двумерной задачи имеются данные расчета в рамках уравнений Навье — Стокса несжимаемой жидкости, а также экспериментальные данные о структуре течения (при Re = 41).

Радиус цилиндра Лі =0.5, радиус внешней границы расчетной области Я2 = 20. Сетка содержала 80 ячеек в окружном и 120 ячеек в радиальном направлениях и сгущалась к поверхности цилиндра по экспоненциальному закону На рис. 40 показал фрагмент сетки в окрестности цилиндра. Во всех расчетах шаг по времени г принимался равным 10.

Число Рейнольдса, посчитанное по диаметру цилиндра и параметрам набегающего потока Re = 40, число Прандтля Рг = 0.7. Вязкость вычислялась по степенной формуле (1.8), с параметром w = 0.5. Безразмерные параметры набегающего потока: р = 1, и = 1, v = 0, р = 1/(7 0

В отличие от теста, приведенного разделе 2.6, здесь на внешней границе во всех проведенных расчетах фиксировалось состояние набегающего потока. На твердой стенке ставилось условие прилипания, реализация которого описана в разделе 1.5.2.

В расчетах по алгоритму CoNS3D-P коэффициент предобуславливания ft в матрице Туркела полагался постоянным во всей расчетной области и равным /3 — М .

На рис. 41 представлены изолинии давления и линии тока около цилиндра, полученные по исходному и предобусловленному алгоритму при М, , = 0.1, 0.01, 0.001. Как и для предыдущей тестовой задачи, мы замечаем, что исходный численный метод достаточно хорошо передает картину течения только при М = 0.1. При числе Маха М = 0.01 длина рециркуляционной зоны меньше, изолинии давления искажаются и становятся не перпендикулярны оси ОХ. Уменьшение числа Маха набегающего потока до 0.001 приводит к еще более сильному искажению решения (см. рис. 41, внизу слева). В то же самое время алгоритм, предобусловленный матрицей Туркела при всех числах Маха дает картину течения, которая практически не зависит от числа Маха и близка к экспериментально зафиксированной, изображенной на рис. 20.

На рис. 42 приведены распределения коэффициента давления Ср по поверхности цилиндра в сравнении с данными из работы [83], где решались несжимаемые уравнения Навье — Стокса. Показаны кривые для исходного метода при М = 0.01, 0.001 (кривая для случая М =0.1 не приведена, она практически совпадает с расчетом [83]), а также для предобусловленного метода при MQQ = 0.01, 0.001, 0.0001. Этот рисунок подтверждает сделанный нами вывод о том, что исходный алгоритм при М 0.1 дает неудовлетворительную точность, в то время как предобусловленный метод сохраняет аппроксимапионные свойства даже при существенно малых числах Маха (М = 0.0001).

Проанализируем истории изменения невязки Err с от числа итераций, изображенные на рис. 43. Из рисунка видно, что при уменьшении числа Маха сходимость исходного метода установления падает: если при М = 0.1 понадобилось провести около 8000 итераций, то при М = 0.01 для достижения предельной ошибки 10 6 необходимо провести около 15000 итераций. При М = 0.001 график истории сходимости выходит на режим периодических колебаний. Необходимо отметить, однако, что во всех случаях попе течения визуально устанавливается уже после 2000 итераций.

Теперь посмотрим, что дает введение предобуславливания. Как и в предыдущем тесте, можно отметить, что при уменьшении числа Маха на порядок, предельная ошибка как предобусловленного, так и исходного алгоритма возрастает на 2 порядка. Графики истории сходимости при М = 0.1 и М = 0.01 можно разделить на две части: крутую (от нуля до 1000-2000 итераций), и более пологую (до 5000-8000 итераций). Причем, как крутая, гак и пологая часть имеют одинаковый наклон при разных числах Маха. Более того, наклон пологой части для предобусловленного алгоритма совпадает с наклоном пологой части исходного при М = 0.1. Это позволяет сделать вывод о том, что невязка численного решения складывается из двух составляющих. Одна из них связана с малым числом Маха, и эта составляющая эффективно подавляется предобуславливанием, а вторая (медленная) составляющая является внутренним свойством численного алгоритма и никак не связанная с малым числом Маха. Вероятнее всего, она обусловлена переключениями в TVD-схеме. После приблизительно 2000 итераций эта вторая составляющая становится доминирующей и не изменяется под действием предобуславливания. Остается открытым вопрос: как устранить эту составляющую, чтобы процесс установления шел с такой скоростью, как на первых 1000 итерациях. Судя по всему, решение этой проблемы лежит за пределами возможностей метода предобуславливания.

По результатам этого теста можно сделать вывод: при расчете задач в рамках уравнений Навье — Стокса предобуславливание конвективной части уравнений матрицей Туркела позволяет восстановить точность и хорошую сходимость численного метода при М —» 0.

Цель данного теста — продемонстрировать работоспособность предложенного численного метода CoNS3D-P рассчитывать течения, в которых присутствуют как зоны, где М 1 и течение практически несжимаемо, так и зоны, где число Маха велико и сжимаемость потока существенна. Рассмотрим течение в плоском симметричном относительно оси ОХ сопле, форма стенок которого у(х) определяется зависимостью

Ширина входного сечения Sin = 5, ширина критического сечения S, = ОД, таким образом отношение сечений составляло Sm/S = 50. Задача решалась в невязкой постановке. Краевые условия для задачи подбирались таким образом, чтобы дозвуковой на входе поток достигал в критическом сечении S, скорости звука (М = 1) и далее ускорялся. Соответствующее число Маха на входе в сопло М = 0.011. В силу симметрии задачи расчеты проводились в области, покрывающей только нижнюю половину сопла. Сетка 70x15 в расчетной области сгущалась к критическому сечению сопла и, незначительно, — к твердой стенке (рис. 44). На верхней границе ставилось условие симметрии. Во входном сечении держались температура торможения То = Т -j , энтропия S = —, соответствующие состоянию (р=1, и = 1, М = 0.011), а также нулевая вертикальная компонента скорости v=0. «-компонента вектора скорости на входной границе экстраполировалась изнутри расчетной области. На выходной границе течение сверхзвуковое, здесь все величины экстраполировались изнутри области.

Похожие диссертации на Численный метод расчета течений сжимаемого вязкого газа в широком диапазоне чисел Маха

Развитие адаптивного алгебраического многосеточного метода расчета многогрупповой диффузии нейтроновЗаславский, Леонид Юрьевич

Методы расчета формообразования поверхности при нестационарной электрохимической обработке Федорова Галина Ильясовна

Математические модели и методы расчета магнитных силовых систем Кивокурцев Александр Юрьевич

Математические модели и численные методы расчета характеристик спутных следов и их воздействия на самолетСудаков Георгий Григорьевич

Разработка высокоточных численных методов расчета пространственного поведения плазмы под действием сильных магнитогидродинамических возмущенийХолодов Ярослав Александрович

Методы расчета рассеяния света осесимметричными частицами в сферическом базисеВинокуров Александр Александрович

Методы расчета и моделирование дискретных стохастических систем с парными взаимодействиямиЛанге Андрей Михайлович

Разработка эффективных методов расчета тонкостенных конструкций с учетом пластических и демпфирующих свойств материалаШишкин Виктор Михайлович

Метод расчета конвективного теплообмена, основанный на вихревом подходеАфанасьева Виктория Викторовна

МЕТОДЫ РАСЧЕТА РАССЕЯНИЯ СВЕТА НЕСФЕРИЧЕСКИМИ ЧАСТИЦАМИ С ИСПОЛЬЗОВАНИЕМ РАЗЛОЖЕНИЯ ПОЛЕЙ ПО ВОЛНОВЫМ ФУНКЦИЯМФАРАФОНОВ ЕВГЕНИЙ ВЯЧЕСЛАВОВИЧ