CC BY
37
УДК 519.632.4:532.516.5
О стационарном решении задачи течения несжимаемой вязкой жидкости при больших числах Рейнольдса
1 2 © А.А. Фомин , Л.Н. Фомина
1 Кузбасский государственный университет им. Т.Ф. Горбачева, г. Кемерово, 650000, Россия
2Кемеровский государственный университет, г. Кемерово, 650043, Россия
Проанализированы вопросы сходимости итерационного процесса и достоверности решений, получаемых методом установления, на примере численного решения задачи стационарного течения несжимаемой вязкой жидкости в плоской квадратной каверне с подвижной верхней крышкой. Задача решается при числах Рейнольдса 15 000 <Re <20 000 и шагах сеточного разбиения 1/128 >h > 1/2048. Показано, что не при всех соотношениях Re и h итерационный процесс установления решения сходится, а полученные стационарные решения достоверны хотя бы на качественном уровне. В системе координат (Re, 1/h) проведен качественный анализ результатов решения задачи с точки зрения сходимости итераций, достоверности получаемых решений и затрат машинного времени.
Ключевые слова: уравнения Навье — Стокса, течение в каверне, сходимость итерационного процесса.
и -►
г ть2 N ст 1
\ № \рвь 2 (Я) 1.
О
Введение. Численному моделированию течения несжимаемой вязкой жидкости в плоской квадратной каверне с подвижной верхней крышкой посвящено множество работ [1-12] в силу того, что эта задача отражает основные особенности, присущие данному классу течений. Схема течения и обозначение возникающих вихрей приведены на рис. 1. Обширный исследовательский материал, накопленный по данной задаче, позволяет авторам работ в режиме сравнительного анализа демонстрировать новые способы разностной аппроксимации исходных дифференциальных уравнений, методы решения получаемых систем линейных алгебраических уравнений (СЛАУ), способы разрешения сингулярностей в верхних углах каверны и т. д. При этом исследователи периодически сталкиваются с проблемой сходимости итерационного процесса нахождения стационарного решения при тех или иных значениях параметра задачи — числа Рейнольдса Re, а также степени подробности сеточного разрешения
Рис. 1. Схема течения и обозначение возможных вихрей
области расчета. Одно время — с первой половины 1990-х годов и до второй половины 2000-х — даже существовало убеждение, что при Re > 104 стационарного решения задачи в принципе не существует (см., например, [3-7]). Однако вышедшие затем работы [8-11] убедительно опровергли эту точку зрения. На сегодняшний день можно с уверенностью говорить о том, что число Рейнольдса в задаче о стационарном течении несжимаемой вязкой жидкости в каверне поднято, как минимум, до уровня 30 000.. .35 000 [10, 11].
Тем не менее вопросы сходимости итераций при построении стационарного решения задачи продолжают оставаться актуальными. Впервые они были, по-видимому, наглядно продемонстрированы в [1], где при Re = 400 на грубой сетке с шагом к = 0,1 возникали осцилляции итерационных приближений решения; на более подробной сетке (к = 0,05) удавалось получить решение с заданным уровнем точности, а при дальнейшем наращивании сеточного разрешения до к = 0,033 процесс сходимости прерывался по исчерпании лимита итераций. Увеличение Re до 1 000 только обострило проблему: с заданной точностью решение не было построено ни для одной из используемых сеток за разумное по времени число итераций. Сорок лет спустя вопрос сходимости итераций при построении стационарного решения по-прежнему остается актуальным, но уже при существенно других параметрах расчетов. Например, в работах [8, 9] для обеспечения сходимости итерационного процесса при росте числа Re авторам приходилось мельчить разностную сетку, поскольку, по их мнению, решающую роль в подавлении осцилляций приближений решения играет величина сеточного числа Рейнольдса Rek. Так, при Re < 10 000 в этих работах использовался сеточный шаг к = 1/256; при 10 000 < Re < 15 000 — к = 1/512; а при Re > 15 000 — к = 1/600 и менее вплоть до к = 1/1024. Следует отметить, что авторам настоящей работы тоже приходилось использовать данный прием в целях обеспечения решения. При этом платой за достижение сходимости итераций были чрезвычайно большие затраты машинного времени [13].
Из вышеприведенного следует, что при численном моделировании течения жидкости практически всегда возникает следующая проблема оптимизации расчетов: как при заданном параметре (параметрах) задачи — в данном случае это число Re — применить минимально подробную сетку и, соответственно, минимизировать расчетное время, получив при этом физически корректное решение. Понятно, что при этом используемый алгоритм нахождения решения должен обеспечить сходимость итераций с заданной точностью.
Настоящее исследование посвящено анализу возникающих при различных сочетаниях Re и к качественно различающихся ситуаций в процессе численного решения задачи о стационарном течении несжимаемой вязкой жидкости в плоской каверне с подвижной верхней крышкой.
Постановка задачи и метод решения. Математическая постановка задачи представляет собой систему безразмерных нестационарных уравнений Навье — Стокса
ды ды ды др 1
— + ы — + V— = —— + — & дх ду дх Яе
дv ду ду др 1
-+ ы-+ V-=--— +-
д1 дх ду ду Яе
ды ду ^
— + — = 0, дх ду
д ы
д ы
дх2 ду2
д ^ ^
дх2 ду2
а также начальных ы = V = 0 и краевых условий: на подвижной крышке каверны ы = 1, V = 0 при у = 1, 0 < х < 1; на остальных стенках каверны ы = V = 0. В начальные моменты времени при 0 < I < 1 крышка по закону синуса плавно разгоняется от состояния покоя до максимального значения 1. Здесь I — время; х, у — декартовы координаты; ы, V — компоненты вектора скорости потока V вдоль координат х и у соответственно; р — давление; Яе — число Рейнольдса, Яе = иЬЫ, где и — максимальная скорость движения крышки; Ь — размер стороны квадратной каверны; V — кинематический коэффициент вязкости жидкости.
Задача решается численно классическим трехшаговым методом расщепления [14] с учетом небольшой модификации [13]: 1) на первом шаге расщепления учитывается давление с предыдущего слоя по времени, а разностные схемы для ы и V записываются в неявном виде; 2) соответственно на втором шаге расщепления задача Неймана формулируется для поправки давления р', которая равна разности давлений на текущем и предыдущем слоях по времени. Разностная аппроксимация производится методом контрольного объема со вторым по пространству и первым по времени порядком аппроксимации с применением экспоненциального профиля, приближенного степенной зависимостью пятой степени [15]. Во всех расчетах использована равномерная сетка. Решение задачи считается найденным (установившимся) при выполнении условия
Оу =
п п-1 п п-1
и - и + V - V
I I
т V п
I
<в.
где п — номер слоя по времени; и, V — векторы сеточных функций компонент скорости ы и V; т — шаг по времени; в — заданная точность установления решения.
Системы разностных эллиптических уравнений решают неявным итерационным полинейным рекуррентным методом, ускоренным в подпространствах Крылова [16]. В качестве критерия сходимости решения СЛАУ используют соотношение
фк _ фк 1
< ^ф Ф\ (1)
с учетом контроля относительного уменьшения первой нормы начальной невязки системы Яф . Здесь Ф — один из векторов сеточных функций: и, V или р', (р' — вектор сеточной функции поправки давления р'); к — номер итерации при решении СЛАУ; — точность сходимости решения СЛАУ для переменной Ф. Выбор данного критерия сходимости связан с удобством контроля количества верных знаков решения системы линейных уравнений. Расчеты показали, что всегда при выполнении условия (1) выполнялось неравенство
Rk
Ф
/
R 0 R Ф
< 0,10^Ф.
Шаг по времени определяют из соотношения т = C min (h, Re h2), где С — число Куранта; h — шаг сеточного разбиения области решения. В силу неявной формы записи разностных схем число Куранта принимали С >> 1, однако при определенных соотношениях Re и h из соображений устойчивости С понижали вплоть до 1. Путем многочисленных вычислительных экспериментов в диапазоне сеточных разбиений от 129x129 до 2049x2049 и числа Рейнольдса 1 000 < < Re < 20 000 было получено, что значения s = 10-5 с избытком хватает для надежного установления численного решения задачи (безотносительно к его физической достоверности хотя бы на качественном уровне). Что касается значения то из соображений его взаимосвязи через С, h и Re с точностью s [13] оно полагалось 10-7 при Ф = u или Ф = v и 10-8 при Ф = р'.
Все расчеты проведены на ПЭВМ Intel Core i5-750, 2,66 ГГц, RAM 12 Гбит.
Примеры решения задачи при Re > 15 000. В силу того что в [4] признаки потери устойчивости расчетного алгоритма имели место для чисел Рейнольдса в диапазоне 12 500 < Re < 16 000, первоначальный расчет настоящего исследования был выполнен при Re = 15 000. На рис. 2 приведены линии тока течения, полученные для различных сеточных разбиений, причем последний фрагмент (рис. 2, г) взят из работы [9], в которой использована разностная аппроксимация четвертого порядка по пространству. Видно, что варианты решения на рис. 2, б—г очень похожи между собой (изображения на рис. 2, в, г практически совпадают), что говорит о насыщении решения в диапазоне сеток от 513x513 до 1025x1025 и его достоверности в силу совпадения
Рис. 2. Картина течения при Re = 15 000 для сеточных разбиений расчетной
области:
а — 257x257; б — 513x513; в — 1025x1025; г — 601x601
с решением из [9]. С другой стороны, картина течения на рис. 2, а принципиально отличается от остальных вариантов решения задачи вследствие наличия артефакта — небольшого пристенного вихря в нижней части правой стенки каверны. Отсюда следует естественный вывод: несмотря на сходимость итераций, существует ограничение снизу по степени подробности сеточного разрешения области, при достижении которого установившееся решение не может считаться правильным.
Повышение числа Рейнольдса до Re = 20 000 (рис. 3) привело к следующим результатам: на грубой сетке 257x257 (рис. 3, а) имеет место установление решения, которое, тем не менее, нельзя признать
Рис. 3. Картина течения при Яе = 20 000 для сеточных разбиений расчетной
области:
а — 257x257; б — 1025x1025; в — 1793x1793; г — 601x601
верным; на сетке 1025x1025 (рис 3, б) решение не установилось с требуемой точностью — имели место постоянные изменения параметров течения без уменьшения значения Оу, поэтому приведенные на данном фрагменте «линии тока» имеют условный смысл. Следует заметить, что изображение на рис. 3, б очень похоже на соответствующие картины «течения» в работах [6, 7], демонстрирующих бифуркации Хопфа решения задачи в случаях отсутствия сходимости итераций. Изображение на рис. 3, в: стационарное решение задачи, полученное на подробной сетке 1793x1793 и совпадающее с соответствующим решением работы [9] (рис. 3, г).
С точки зрения поведения динамических систем процесс установления решения представлен на рис. 4, на котором параметриче-
ские кривые в координатах (1§ Бу, 1§ || р' Ц1) демонстрируют взаимосвязь степени установления решения Бу и первой нормы поправки давления || р' ||ь При этом в качестве параметра выступает время
а б
Рис. 4. Параметрические кривые установления решения при Яе = 20 000 в зависимости от сеточного разрешения расчетной области:
а — 1025x1025; б — 1793x1793
В случае установления решения обе величины, Бу и || р' Ц1, становятся заведомо малыми (рис. 4, б). В противном случае (отсутствие сходимости итераций) их значения непрерывно и, вообще говоря, хаотично изменяются, оставаясь все это время ограниченными снизу некоторыми предельными значениями, не удовлетворяющими условию установления. На рис. 4, а хорошо видно, что за все время расчета варианта задачи на сетке 1025x1025 величина Бу не опускается ниже стартового значения 5,20-10 3. Стрелками на рисунке показаны направления перемещения вдоль кривых по мере увеличения параметра Из приведенных графиков и результатов других расчетов следует, что процесс установления решения при больших Яе, как правило, сопряжен со стадией хаотического изменения параметров Бу и || р' ||ь В дальнейшем он либо остается на этой стадии навсегда, либо переходит в стадию монотонной сходимости итераций вплоть до достижения наперед заданной точности установления решения.
Для понимания причин отсутствия сходимости итерационного процесса на сетке 1025x1025 и, наоборот, сходимости итераций на более мелкой сетке 1793x1793 следует обратиться к выводам работ [9, 17], согласно которым сходимость итерационного процесса напрямую связана с величиной сеточного числа Рейнольдса: чем оно меньше, тем больше вероятность того, что решение в конечном счете установится.
На рис. 5 приведены поля величины ЯеИ = Яе |У| И, соответствующие картинам течения, приведенным на рис. 3, а-в. В качестве критического значения выбрана величина ЯеИ = 2, поскольку ее превышение может приводить к так называемым пилообразным осцил-ляциям решения, не совместимым с физически осмысленными результатами расчетов [18]. На рис. 5, а область ЯеИ > 2 занимает почти всю каверну и, тем не менее, итерации сходятся. Следовательно, здесь основную роль играют другие факторы, такие, например, как
0,4 0,6 0,8 1,0 б
Рис. 5. Поля сеточного числа Рей-нольдса Res для сеточных разбиений расчетной области:
а — 257x257, max Res = 78,2; б — 1025x1025, max Re* = 19,5; в — 1793x1793, max Res = 11,2
в
стабилизирующие влияние схемной вязкости (при ЯеИ > 2 схема Па-танкара схожа с противопотоковой [15]) и (или) понижающее число обусловленности матрицы СЛАУ из-за относительно небольшого количества уравнений — порядка 65 000. Что касается более подробных сеток (рис. 5, б, в), то тут уже превалирующую роль в смысле сходимости итерационного процесса играет величина ЯеИ. И хотя площадь областей ЯеИ > 2 на этих фрагментах приблизительно одинакова, конфигурация проходящих вблизи ЯеИ = 2 изолиний, а также величины максимумов ЯеИ и наличие областей возмущений решения в нижних
углах каверны на фрагментах рис. 5, б объясняют, почему на сетке 1793x1793 итерации сошлись, а на сетке 1025x1025 — нет.
Качественный анализ результатов решения задачи. Приведенные результаты расчетов говорят о том, что при различных сочетаниях Яе и к процесс установления решения не всегда позволяет получить стационарное решение, а в случае успешной сходимости итераций сами решения могут быть как физически достоверными, так и наоборот. Сюда же следует отнести вопрос принудительного ограничения времени расчета варианта, поскольку вычисления не могут длиться бесконечно: их продолжительность всегда ограничена, как правило, предельным количеством итераций или астрономическим временем исполнения программы. Соответственно, из самых общих соображений в системе координат входных параметров расчетов (Яе, 1/к) можно провести линии, которые качественно разграничивают координатную плоскость на зоны, в которых реализуются или не реализуются следующие факторы: сходимость итераций, достоверность решения, ограничение по итерациям. Такие линии представлены на рис. 6. Приведем обоснование поведения каждой из них на качественном уровне.
Рис. 6. Принципиальная схема результативной области П0 и переходов между вариантами решений задачи, полученных различными
исследователями
Линия АВ — достоверность решения. Увеличение числа Рей-нольдса приводит к уменьшению характерных размеров структурных элементов течения (например, толщины пограничного слоя вдоль стенок), что, в свою очередь, при фиксированном к приводит к огрублению картины течения. Следовательно, в целях обеспечения достоверности решения для больших значений Яе требуется большее разрешение, иными словами, между Яе и 1/к имеет место прямая пропорциональная зависимость.
Линия FE — ограничение по итерациям. Как показывает практика расчетов [1, 13], не только увеличение сеточного разрешения (уменьшение h), но и увеличение числа Re приводит к резкому росту времени вычислительного процесса при прочих равных условиях (в первую очередь — при условиях точности сходимости решения). Поскольку должно существовать разумное ограничение такого времени, то связь между Re и 1/h через это ограничение должна носить обратно пропорциональный характер: чем больше значение Re задачи, тем менее подробным должно быть сеточное разрешение, чтобы за разумное время вычислений построить решение, и наоборот.
Линия BCE — сходимость итераций. Как упоминалось выше, в многочисленных публикациях 1990-х и начала 2000-х годов утверждается, что неограниченное увеличение числа Re в конечном счете приводит к прекращению сходимости итераций при построении стационарного решения задачи. Соответственно, отсутствие сходимости итерационного процесса или, наоборот, ее наличие можно объяснить двумя причинами. Приведем их.
Линия BC — левая ветка. Увеличение разрешения приводит к увеличению размерности системы разностных уравнений, откуда следует увеличение числа обусловленности матрицы разностной СЛАУ, что, в свою очередь, увеличивает ошибку вектора решения системы алгебраических уравнений. Как следствие, в результате накопления от итерации к итерации ошибки вектора решения возникает потеря устойчивости итерационного процесса, которая, в общем случае, не позволяет достичь установления решения всей задачи. Соответственно, между величинами Re и 1/h имеет место обратно пропорциональная зависимость, поскольку увеличение одной из них для сохранения сходимости итераций необходимо компенсировать уменьшением другой.
Линия CE — правая ветка. Согласно рекомендациям [9, 17], подтвержденным авторскими расчетами, уменьшение сеточного Reh позволяет добиться сходимости итерационного процесса построения стационарного решения задачи, откуда следует, что увеличение числа Re должно компенсироваться уменьшением сеточного шага h, т. е. между Re и 1/h имеет место прямая пропорциональная зависимость.
В итоге часть координатной плоскости, ограниченная замкнутой линией ABCEF, образует так называемую результативную область Q0 — множество таких координатных пар (Re, 1/h), использование которых в качестве входных параметров позволяет получить хотя бы качественно правильное стационарное решение задачи. На этом же рисунке приведены переходы между вариантами решений по параметрам Re и 1/h, построенные по материалам работ О.Р. Бургграфа (Burggraf) [1] и Э. Эртурка (Erturk) с соавторами [10, 17]. Обоснование второй линии устойчивости рассмотрим далее.
Нетрудно видеть, что в случае справедливости приведенной конфигурации области Q0 должны существовать такие пары входных параметров (Re, 1/h), которые обеспечили бы цепочку переходов между ва-
риантами решений задачи, обозначенной на рис. 6 узлами а-ж. Такие параметры действительно были подобраны, и полученные на их основе картины течения представлены на рис. 7. При этом, строго говоря,
Рис. 7. Картины течений, соответствующие входным параметрам задачи:
а — Яе =15 000, И = 1/256; б — Яе =15 000, И = 1/320; в — Яе =16 250, И = 1/320; г — Яе =16 250, И = 1/352; д — Яе =16 250, И = 1/512; е — Яе =16 250, И = 1/1024
рис. 7, д не является решением с заданной точностью, поскольку в этом расчете итерации не сошлись, а минимально достигнутое значение Пу составляло 1,4-10 3. Также понятно, что пару входных параметров (Яе, 1/к), соответствующую узлу ж на рис. 6, для любого фиксированного Яе нетрудно определить подходящим выбором малого сеточного шага к. Дело в том, что в основе установления решения лежит механизм поточечной сходимости [14], приводящий при увеличении количества расчетных узлов (уменьшении к) к катастрофическому нарастанию количества итерационных шагов по времени [13]. В итоге можно заключить, что построенная из общих качественных рассуждений схема области О0 получила частное количественное подтверждение на примере вариантов решений задачи, приведенных на рис. 7.
Для большего обоснования справедливости схемы на рис. 6 были проведены систематические расчеты в диапазоне чисел Рейнольдса 15 000 < Яе < 20 000 и сеточных шагов 1/128 > к > 1/2048. Результаты расчетов с точки зрения сходимости итераций, физической корректности решения и превышения допустимого количества итераций по времени представлены в таблице в виде чисел Куранта, использованных при расчете различных вариантов.
Результаты решений задачи в зависимости от числа Яе и шага сеточного разбиения к
Яе 1/к
128 256 320 352 384 448 512 768 896 1024 1280 1536 1792 2048
20 000 3 16 3 8 3 1 3 (1) (3) (8) (16) (16) (16) (32) (16) (32) (32) (48) (64) (96) (32) 64 (96) 80
18 750 3 16 3 8 3 1 (1) (3) (8) (16) (16) (16) (32) (16) 32 (64) 64 64 80
18 250 3 16 3 8 3 (1) (1) (3) (8) (16) (16) (16) (32) (64) 48 64 64 80
17 500 3 16 3 8 3 (1) (1) (3) (3) (3) (8) (32) (3) (16) (32) (3) (16) (32) 32 32 (16) 32 64 64 80
16 250 3 16 3 8 3 1 (1) (2) (3) (1) (3) (3) (8) (3) 16 (32) 32 32 32 48 48 72
15 000 3 16 3 8 1 3 1 3 4 8 32 32 32 32 32 48 70
В левой части таблицы для диапазона шагов 1/128 > к > 1/352 курсивом обозначены те расчеты (числа Куранта), в которых итерационный процесс сошелся, но сами решения получились качественно
неверными, аналогичными тем, что приведены на рис. 7, а, в. В правой части таблицы для к = 1/2048 полужирным выделены результаты, каждый из которых был получен со значительными затратами машинного времени (более 250 ч), за счет чего они отнесены к случаю превышения ограничения по итерациям (для сравнения: при построении решения при Яе = 15 000 и к = 1/2048 было затрачено приблизительно 180 ч). Понятно, что предельно допустимое количество итераций (время вычислений) — величина субъективная, однако в данном случае речь идет не о конкретном значении этого параметра, а о его существовании в принципе, и для определенности такой предел выбран в размере 250 ч.
Данные в скобках в средней части таблицы отвечают за расчеты, в процессе которых не удалось достигнуть стационарного решения. Для этих вариантов минимальные значения Пу изменяются в пределах от 0,33 (число Куранта С = 3 при Яе = 17 500 и к = 1/896) до 4,0 10-5 (число Куранта С = 48 при Яе = 20 000 и к = 1/1536). В случае некоторых пар входных параметров (Яе, 1/к) было выполнено несколько расчетов при различных значениях числа Куранта, для того чтобы убедиться в отсутствии установления решения либо все-таки отыскать его, как это было сделано для следующих пар входных параметров: Яе = 16 250 и к = 1/768, Яе = 18 750 и к = 1/1280, Яе = 20 000 и к = 1/1792. Обычным начертанием в таблице представлены варианты расчетов, в которых удалось получить стационарные физически корректные решения.
Конечно, нельзя утверждать со стопроцентной уверенностью, что найденные разграничения между расчетами, в которых итерации сошлись и не сошлись, абсолютно верны, тем более что они зависят и от других факторов, таких как заданная точность сходимости и (или) точность и метод решения СЛАУ и т. д. Главное здесь другое — систематические расчеты показывают, что при больших числах Рей-нольдса по параметру 1/к между зонами стационарных физически корректных и некорректных решений существует зона, в которой решение не устанавливается в принципе. При этом конфигурация группы ячеек таблицы, соответствующих стационарным и достоверным хотя бы на качественном уровне решениям, очень похожа на конфигурацию области 00 на рис. 6.
Еще одной особенностью полученных результатов является понижение числа Куранта от левой и правой границ таблицы к центру — прежде всего здесь следует обратить внимание на строку Яе = 15 000. Любопытно, что уменьшение шага по времени не всегда приводит к сходимости итерационного процесса, о чем свидетельствуют, например, результаты расчетов при Яе = 16 250, к = 1/768 и другие.
Еще одной интересной особенностью являются результаты для Яе = 18 750 и Яе = 20 000 при к = 1/352. Вопреки первоначальным обсуждениям схемы на рис. 6, итерационный процесс здесь сошел-
ся, хотя картина течения получилась такой же физически недостоверной, как и решения в соседних слева ячейках таблицы. Объяснение здесь кроется в следующем: первоначальные общие рассуждения не касались способа аппроксимации исходных дифференциальных уравнений. Однако известно, что в случае использования в конвективных членах противопотоковых разностей процесс установления ведет себя устойчиво, поскольку схемная вязкость гасит возмущения приближений решения, препятствующих его сходимости [18]. С другой стороны, использование центрально-разностной схемы аппроксимации конвективных членов на грубых сетках (большие к) приводит к превышению Яек значения 2 и, следовательно, к возможности возникновения пилообразных осцилляций, а значит, в данной ситуации — к возможному отсутствию сходимости итераций. Поскольку применяемая в настоящей работе схема Патанкара при Яек < 2 схожа с центрально-разностной аппроксимацией, а при Яек > 2 — с противопотоковой, то существует вторая линия устойчивости, которая, в отличие от остальных линий схемы, «действует» только вверх. Это означает, что в расчетах при параметрах Яе и 1/к выше этой линии итерации будут сходиться, в противном случае — в зависимости от других вышеперечисленных факторов. Понятно, что корреляция Яе и 1/к при этом должна быть прямо пропорциональной, поскольку коэффициент схемной вязкости пропорционален сеточному шагу к, откуда следует, что чем меньше к, тем большим значениям Яе соответствует область преобладания схемной вязкости над физической. В подтверждение данного вывода для параметров Яе = 20 000, к = 1/352 выполнен расчет с использованием центрально-разностной аппроксимации конвективных членов и числа Куранта С = 1. В результате решение не установилось, а минимальное значение Пу составило 0,66.
Затраты машинного времени для вариантов, в которых было достигнуто установление решения, значительно различаются по мере увеличения степени подробности сеточного разбиения области. Так, на сетке 129x129 для установления решения потребовалось от 1,5 до 10 мин, а на сетке 2049x2049 — от 180 до 468 ч в зависимости от значений чисел Куранта и Рейнольдса. При этом увеличение первого числа приводило к ускорению сходимости решения, а увеличение второго — к его замедлению.
Заключение. Проанализированы результаты численного решения задачи стационарного течения несжимаемой вязкой жидкости в плоской каверне с подвижной верхней крышкой при больших числах Рейнольдса (15 000 < Яе < 20 000). Показано, что в зависимости от соотношений входных параметров вычислительного алгоритма числа Яе и сеточного шага к могут возникать различные ситуации как по сходимости итерационного процесса решения задачи, так и по физической корректности получаемого решения. В системе координат (Яе, 1/к) на основании самых общих рассуждений показано наличие так
называемой результативной области Q0 — совокупности точек, координаты которых, взятые в качестве входных параметров, обеспечивают сходимость итерационного процесса установления решения и физическую корректность решения задачи хотя бы на качественном уровне.
По опубликованным данным исследований и результатам проведенных систематических расчетов для больших чисел Рейнольдса в диапазоне сеточных шагов 1/128 > h > 1/2048 можно сделать следующие выводы:
• сходимость итерационного процесса установления решения на грубых сетках зависит от степени влияния схемной вязкости и (или) числа обусловленности матриц разностных эллиптических СЛАУ, а на подробных сетках — от величины сеточного числа Рейнольдса Reh;
• при больших числах Re сходимость итераций достигается либо на очень грубых сетках, либо на очень подробных, причем ширина зоны отсутствия сходимости итераций по параметру 1/h увеличивается по мере увеличения числа Re;
• поскольку на грубых сетках решение физически некорректно, а использование подробных сеток приводит к большим затратам машинного времени, дальнейший рост числа Рейнольдса в задаче сопряжен с повышением порядка аппроксимации исходных дифференциальных уравнений.
ЛИТЕРАТУРА
[1] Burggraf O.R. Analytical and numerical studies of the structure of steady separated flows. J. of Fluid Mechanics, 1966, vol. 24, pp. 113-151.
[2] Ghia U., Ghia K.N., Shin C.T. High-Re solution for incompressible flow using the Navier — Stokes equations and a multigrid method. J. of Computational Physics, 1982, vol. 48, pp. 387-411. Available at: http://dx.doi.org/10.1016/ 0021-9991(82)90058-4
[3] Bruneau C.-H., Jouron C. An efficient scheme for solving steady incompressible Navier Stokes equations. J. of Computational Physics, 1990, vol. 89, pp. 389-413. URL: http://dx.doi.org/10.1016/0021-9991(90)90149-U
[4] Barragy E., Carey G.F. Stream function-vorticity driven cavity solution using p finite elements. Computers & Fluids, 1997, vol. 26, no. 5, pp. 453-468. Available at: http://dx.doi.org/10.1016/S0045-7930(97)00004-2
[5] Marinova R.S., Christov C.I., Marinov T.T. A fully coupled solver for incompressible Navier-Stokes equations using operator splitting. Int. J. of Computational Fluid Dynamics, 2003, vol. 17, issue 5, pp. 371-385. Available at: http://dx.doi.org/10.1080/1061856031000114300
[6] Bruneau C.-H., Saad M. The 2D lid-driven cavity problem revisited. Computers & Fluids, 2006, vol. 35, pp. 326-348. Available at: http://dx.doi.org/10.1016/ j.compfluid.2004.12.004
[7] Kumar D.S., Kumar K.S., Kumar M.D. A fine grid solution for a lid-driven cavity flow using multigrid method. Engineering Applications of Computational Fluid Mechanics, 2009, vol. 3, no. 3, pp. 336-354. Available at: http://dx.doi.org/ 10.1080/19942060.2009.11015275
[8] Erturk E., Corke T.C., Gokgol C. Numerical solutions of 2-D steady incompressible driven cavity flow at high Reynolds numbers. Int. J. for Numerical
Methods in Fluids, 2005, vol. 48, pp. 747-774. Available at: http://dx.doi.org/10.1002/fld.953
[9] Cardoso N., Bicudo P. Time dependent simulation of the Driven Lid Cavity at High Reynolds Number. ArXiv: D809.3098v2[physics.fly-dyn], 2009, pp. 1-20. Available at: http://arxiv.org/pdf/0809.3098.pdf
[10] Erturk E., Gokgol C. Fourth-order compact formulation of Navier — Stokes equations and driven cavity flow at high Reynolds numbers. Int. J. for Numerical Methods in Fluids, 2006, no. 50, pp. 421-436. Available at: http://dx.doi.org/ 10.1002/fld. 1061
[11] Wahba E.M. Steady flow simulation inside a driven cavity up to Reynolds number 35000. Computers & Fluids, 2012, vol. 66, pp. 85-97. Available at: http://dx.doi.org/10.1016/j.compfluid.2012.06.012
[12] Басараб М.А. Численно-аналитический метод решения двумерных задач естественной конвекции в замкнутых полостях. Математическое моделирование и численные методы, 2014, № 1, с. 18-35.
[13] Фомин А. А., Фомина Л.Н. Неявный итерационный полинейный рекуррентный метод в применении к решению задач динамики несжимаемой вязкой жидкости. Компьютерные исследования и моделирование, 2015, т. 7, № 1, с. 35-50.
[14] Белоцерковский О.М., Гущин В. А., Щенников В.В. Метод расщепления в применении к решению задач динамики вязкой несжимаемой жидкости. Журнал вычислительной математики и математической физики, 1975, т. 15, № 1, с. 197-207.
[15] Патанкар С. Численные методы решения задач теплообмена и динамики жидкости. Пер. с англ. Москва, Энергоатомиздат, 1984, 152 c.
[16] Фомин А.А., Фомина Л.Н. Ускорение полинейного рекуррентного метода в подпространствах Крылова. Вестник Томского государственного университета. Математика и механика, 2011, № 2 (14), c. 45-54.
[17] Erturk E. Discussions on driven cavity flow. Int. J. for Numerical Methods in Fluids, 2009, vol. 60, pp. 275-294. Available at: http://dx.doi.org/ 10.1002/fld.1887
[18] Роуч П. Вычислительная гидродинамика. Москва, Мир, 1980, 616 c.
Статья поступила в редакцию 03.09.2015
Ссылку на эту статью просим оформлять следующим образом:
Фомин А.А., Фомина Л.Н. О стационарном решении задачи течения несжимаемой вязкой жидкости при больших числах Рейнольдса.
Математическое моделирование и численные методы, 2015, № 4 (8), с. 92-102.
Фомин Александр Аркадьевич родился в 1958 г., окончил Томский государственный университет в 1980 г. Канд. физ.-мат. наук, ведущий эксперт отдела развития и международного сотрудничества Кузбасского государственного технического университета им. Т.Ф. Горбачева. Автор монографии и 15 статей в научных рецензируемых журналах. e-mail: [email protected].
Фомина Любовь Николаевна родилась в 1956 г., окончила Томский государственный университет в 1978 г. Канд. физ.-мат. наук, доцент кафедры вычислительной математики Кемеровского государственного университета. Автор монографии и 9 статей в научных рецензируемых журналах. e-mail: [email protected]
On stationary solution of the problem of an incompressible viscous fluid at high Reynolds numbers
© A.A. Fomin1, L.N. Fomina2
:T. F. Gorbachev Kuzbass State Technical University, Kemerovo, 650000, Russia
2Kemerovo State University, Kemerovo, 650043, Russia
The research explored questions of the convergence of iterative processes and correctness of the solutions on the example of the problem about a steady-state flat square lid-driven cavity flow of incompressible viscous liquid. The problem is solved for Reynolds numbers of 15000 < Re <20000 and steps of grid 1/128 >h >1/2048. The findings of the research illustrate that not for all relationships between Re and h the convergence of iterative processes is stable and the resulting steady-state solutions are qualitatively correct. We conducted a qualitative analysis of the solutions of the problem in the coordinate system (Re, 1/h) in terms of the convergence of iterative process, solution correctness and the required computing time. According to the literature and the results of systematic calculations we conclude that the stability of the convergence of iterative process on the coarse grid depends on the degree of influence of the artificial viscosity and/or the condition number of the matrix of difference elliptical linear algebraic equations, and on the detailed grid it depends on the grid Reynolds number. At high Reynolds numbers steady calculations can be carried out either on very coarse grids, or on very detailed ones. The width of the zone of instability in terms of parameter 1/h increases with increasing Reynolds number. Since the coarse grid solution is incorrect, and the use of detailed grid leads to very high costs of computer time, the further increase of the Reynolds number in the problem is associated with increasing the order of approximation of the differential equations.
Keywords: Navier—Stokes equations, a lid-driven cavity flow, a convergence of iterative process.
REFERENCES
[1] Burggraf O.R. Analytical and numerical studies of the structure of steady separated flows. Journal of Fluid Mechanics, 1966, vol. 24, pp. 113-151.
[2] Ghia U., Ghia K.N., Shin C.T. High-Re solution for incompressible flow using the Navier-Stokes equations and a multigrid method. Journal of Computational Physics, 1982, vol. 48, pp. 387-411. Available at: http://dx.doi.org/10.1016/ 0021-9991(82)90058-4
[3] Bruneau C.-H., Jouron C. An efficient scheme for solving steady incompressible Navier Stokes equations. Journal of Computational Physics, 1990, vol. 89, pp. 389-413. Available at: http://dx.doi.org/10.1016/0021-9991(90)90149-U
[4] Barragy E., Carey G.F. Stream function-vorticity driven cavity solution using p finite elements. Computers & Fluids, 1997, vol. 26, no. 5, pp. 453-468. Available at: http://dx.doi.org/10.1016/S0045-7930(97)00004-2
[5] Marinova R.S., Christov C.I., Marinov T.T. A fully coupled solver for incompressible Navier — Stokes equations using operator splitting. Int. J. of Computational Fluid Dynamics, 2003, vol. 17, iss. 5, pp. 371-385. Available at: http://dx.doi.org/10.1080/1061856031000114300
[6] Bruneau C-H., Saad M. The 2D lid-driven cavity problem revisited. Computers & Fluids, 2006, vol. 35, pp. 326-348. Available at: http://dx.doi.org/10.1016/ j.compfluid.2004.12.004
[7] Kumar D.S., Kumar K.S., Kumar M.D. A fine grid solution for a lid-driven cavity flow using multigrid method. Engineering Applications of Computa-
tional Fluid Mechanics, 2009, vol. 3, no. 3, pp. 336-354. Available at: http://dx.doi.org/10.1080/19942060.2009.11015275
[8] Erturk E., Corke T.C., Gokgol C. Numerical solutions of 2-D steady incompressible driven cavity flow at high Reynolds numbers. Int. J. for Numerical Methods in Fluids, 2005, vol. 48, pp. 747-774. Available at: http://dx.doi.org/ 10.1002/fld.953
[9] Cardoso N., Bicudo P. Time dependent simulation of the Driven Lid Cavity at High Reynolds Number, arXiv: D809.3098v2[physics.fly-dyn], 20 Nov 2009, pp. 1-20. Available at: http://arxiv.org/pdf/0809.3098.pdf
[10] Erturk E., Gokgol C. Fourth-order compact formulation of Navier-Stokes equations and driven cavity flow at high Reynolds numbers. Int. J. for Numerical Methods in Fluids, 2006, no. 50, pp. 421-436. Available at: http://dx.doi.org/ 10.1002/fld.1061
[11] Wahba E.M. Steady flow simulation inside a driven cavity up to Reynolds number 35000. Computers & Fluids, 2012, vol. 66, pp. 85-97. Available at: http://dx.doi.org/10.1016/j.compfluid.2012.06.012
[12] Basarab M.A. Matematicheskoe modelirovanie i chislennye metody — Mathematical Modeling and Computational Methods, 2014, no. 1, pp. 18-35.
[13] Fomin A.A., Fomina L.N. Kompyuternye issledovaniya i modelirovanie — Computer Research and Modeling, 2015, vol. 7, no. 1, pp. 35-50.
[14] Belotserkovskiy O.M., Gushchin V.A., Shchennikov V.V. Zhurnal vychislit-elnoy matematiki i matematicheskoy fiziki — USSR Comp. Math. Math. Phys.,
1975, vol. 15, no. 1, pp. 190-200. Available at: http://dx.doi.org/10.1016/ 0041-5553(75)90146-9
[15] Patankar S.V. Numerical Heat Transfer and Fluid Flow. Hemisphere Publishing Co, 1980, 197 p. [In Russ.: Patankar S. Chislennye metody resheniya zadach teploobmena i dinamiki zhidkosti. Moscow, Energoatomizdat Publ., 1984, 152 p.].
[16] Fomin A.A., Fomina L.N. Vestnik Tomskogo gosudarstvennogo universiteta, Matematika i mekhanika — Tomsk State University J. of Mathematics and Mechanics, 2011, no. 2 (14), pp. 45-54.
[17] Erturk E. Discussions on driven cavity flow. Int. J. for Numerical Methods in Fluids, 2009, vol. 60, pp. 275-294. Available at: http://dx.doi.org/ 10.1002/fld.1887
[18] Roache P.J. Computational Fluid Dynamics. Albuquerque, Hermosa Publs,
1976, 446 p. [In Russ.: Rouch P. Vychislitelnaya gidrodinamika. Moscow, Mir Publ., 1980, 616 p.].
Fomin A.A. (b.1958) graduated from Tomsk State University in 1980. Cand. Sci. (Phys. & Math.), the leading expert of the Department of Development and International Cooperation at T.F. Gorbachev Kuzbass State Technical University. Author of a monograph and 15 papers published in scientific peer-reviewed journals. SPIN-code: 3326-4190. e-mail: [email protected].
Fomina L.N. (b.1956) graduated from Tomsk State University in 1978. Cand. Sci. (Phys. & Math.), Associated Professor, Department of Computational Mathematics at Kemerovo State University. The author of a monograph and 9 papers published in scientific peer-reviewed journals. SPIN-code: 5658-7206. e-mail: [email protected]
