CC BY
306
2014 Математика и механика № 3(29)
УДК 519.632.4:532.516.4
А.А. Фомин, Л.Н. Фомина
ЧИСЛЕННОЕ РЕШЕНИЕ УРАВНЕНИЙ НАВЬЕ - СТОКСА ПРИ МОДЕЛИРОВАНИИ ДВУМЕРНЫХ ТЕЧЕНИЙ ВЯЗКОЙ НЕСЖИМАЕМОЙ ЖИДКОСТИ
Анализируется эффективность применения неявного итерационного полиней-ного рекуррентного метода решения систем разностных эллиптических уравнений, возникающих при численном моделировании двумерных течений вязкой несжимаемой жидкости. Исследование проводится на примере задачи о стационарном двумерном течении в квадратной каверне с подвижной крышкой, сформулированной в естественных переменных (и, V, р). Показано, что применение полинейного рекуррентного метода позволяет заметно сократить общее время решения задачи. В качестве иллюстрации к достигнутым результатам приведена структура течения при Яе = 15 000.
Ключевые слова: течение в каверне, уравнения Навье - Стокса, неявный полинейный рекуррентный метод.
Задача о течении вязкой несжимаемой жидкости в квадратной каверне с подвижной крышкой является своеобразным «полигоном», на котором вот уже несколько десятков лет многочисленные исследователи демонстрируют свои достижения в области численного решения уравнений Навье - Стокса. Привлекательность этой задачи объясняется, с одной стороны, простотой её формулировки: решение зависит фактически от одного параметра - числа Рейнольдса. А с другой - наличием разрывов и сингулярностей полей решения в верхних углах каверны и фактически бесконечные цепочки вихрей в каждом из нижних углов гарантируют возможность постоянного улучшения решения и, как следствие, постоянный интерес исследователей к этой задаче. Накопленные за последнее время обширные материалы позволяют адекватно оценивать степень правильности новых результатов, обеспечивая тем самым их высокую достоверность.
Среди работ на эту тему до сих пор большой популярностью заслуженно пользуется работа [1], сформировавшая современный подход к изложению материалов исследования. Несмотря на использование умеренных по современным меркам равномерных сеток типа 129x129 и 257x257, авторам удалось построить решение задачи по числам Рейнольдса вплоть до Яе = 10 000 и оценить параметры вихрей вплоть до правого нижнего вихря третьего уровня.
В абсолютном большинстве работ исследуются характеристики течения при Яе < 10 000. Более того, ряд авторов высказали предположение, что не существует стационарного решения задачи при больших значениях числа Рейнольдса. Однако некоторые публикации последнего времени опровергают это предположение. В частности, в работе [2] приведены параметры вихрей при Яе = 12 500 вплоть до четвёртого уровня (так называемые БЬ4 и БЯ4) в нижних углах каверны, но при дальнейшем увеличении числа Яе отмечается неустойчивость решения. А в работе [3] (и ряде других публикаций тех же авторов) удалось поднять число Рей-нольдса до 21 000, но при этом, в силу равномерности разностных сеток, не удалось выявить вихри четвёртого уровня. При этом там же утверждается, что эф-
фективным средством борьбы против неустойчивости являются более подробные сетки. В частности авторами были использованы сетки 257x257 вплоть до Яе = 10 000, при 10 000 < Яе < 15 000 - сетки 513x513, при Яе > 15 000 - сетки 601x601 и подробнее вплоть до 1025x1025.
Следует отметить, что приблизительно половина исследователей предпочитают постановку задачи в переменных у-ю. Возможно, это связано с тем, что основное сравнение результатов осуществляется по значениям экстремумов полей у-ю (центров вихрей), в то время как использования естественных переменных (и, V, р), несмотря на удобство и абсолютную точность формулировки граничных условий, требует численного восстановления полей ю и у, что не может не сказаться на точности представления сравниваемых параметров течения. Среди последних работ в переменных (и, V, р) следует отметить [4] и другие публикации тех же авторов. Использование высокого порядка аппроксимации и кусочно-равномерной сетки - свыше 65 тысяч расчетных узлов (так называемой сетки М3) - позволило получить высокоточные результаты при Яе = 1000 по параметрам вихрей вплоть до третьего уровня.
На базе литературных данных можно составить обобщенную схему стационарного течения (см. рис. 1). В ее основу положено решение конкретного варианта задачи при Яе > 10 000. При этом необходимо заметить, что вихри третьего и четвертого уровней в нижних углах изображены вне масштаба - в расчетах они получаются значительно меньших размеров.
В1-4 ВР4 X
Рис. 1. Система координат и обобщенная схема течения
В настоящем исследовании математическая постановка задачи в безразмерном виде представляет собой систему нестационарных двумерных уравнений Навье -Стокса:
И + (V-V)и =-Ухр + Яе-1 V2и ; (1)
V + (V -V) V = -У р + Яе-1 V^ ; (2)
VV = 0,
(3)
а также начальных и = V = 0 и краевых условий: на верхней границе области
0 < х < 1, у = 1: и = 1, V = 0, (4)
и на остальных границах
и = V = 0. (5)
Здесь t - время, х, у - пространственные координаты, и, V - х, у компоненты вектора скорости потока V, р - давление, Яе = иЫи - число Рейнольдса, и - кинематический коэффициент вязкости жидкости. В начальные моменты времени 0 < t < 1 скорость движения крышки плавно увеличивается до 1 и далее не меняется.
Для построения линий тока ищется решение вспомогательной задачи относительно функции тока у. Внутри области решения
Ду = -ю, (6)
на границах области решения
у = 0. (7)
Здесь ю = ду/дх - ди/ду - завихренность потока.
Задача (1)-(5) решается численно. Область решения покрывается так называемой МАС-сеткой [5], в которой узлы определения давления располагаются внутри прямоугольной сеточной ячейки, а узлы определения компонент скорости - на её гранях. В общем случае сетка слабонеравномерная, то есть \к, - к,-11 ~ 0(к2), где к - сеточный шаг. Степень неравномерности определяется параметром сжатия к0, равным отношению среднего шага сетки к минимальному.
На каждом шаге по времени реализуется классическая трёхшаговая схема расщепления [6], на первом шаге которой определяются предварительные значения компонент скорости:
и* = и"-1 -т[^"-1 -V)и *-Яе-1 V2и * +Ухр"-1 ] ; (8)
V* = V"-1 -т[^"-1 -V) V * -Яе-1 V2у * +Vyp"-1 ] . (9)
На втором шаге определяется поле поправки давления р' путём решения уравнения Пуассона
Др'=т^ V*. (10)
И наконец, на третьем шаге поля скорости и давления корректируются по формулам
V" = V* -^р\ р" = р"-1 +стр'. (11)
Здесь т - шаг по времени, " - номер слоя по времени, ст - итерационный параметр. Граничные условия для (8), (9) берутся по аналогии с (4), (5) с учетом естественной замены и ^ и*, V ^ V*; а для (10) - условия Неймана, которые следуют из (11) в силу стремления V* ^ V" при " ^ да. Во всех расчётных вариантах методом установления ищется стационарное решение задачи (1)-(5). Путём первоначального численного эксперимента было определено оптимальное значение ст = 1,8, которое в дальнейшем не менялось.
Проверка полной аппроксимации уравнения движения
V" = V "-1 - т [(V "-1 - V) (V" + ^р') - Яе-1 V2 (V" + ^р') - стVp' + Vp" ] =
= V"-1 -т[(У"-1 " - 2V" +Vр" -стVp']-т2 [(V"-1 -V)Vp ' - Re-1V2^р ')] =
= V"-1 -т[(У"-1 " -Re-1V2V" +Vp"] + О^р') (12)
говорит о том, что при стремлении Vp' ^ 0 используемая трёхшаговая схема аппроксимирует уравнение движения неявным образом, что отличает данную схему от явной схемы в [6]. Подобный подход позволяет выбирать числа Куранта значительно больше 1, ускоряя тем самым процесс установления решения.
Разностная аппроксимация (8) - (10) производится методом контрольного объёма [7] со вторым порядком по пространству и первым по времени. Здесь необходимо обратить внимание на особенности получаемых систем разностных эллиптических уравнений. Для уравнений движения (8), (9) эти системы легко разрешимы: оценка чисел обусловленности по первой норме матриц данных систем в диапазоне сеточного разрешения от 65x65 до 1025x1025 (сетки как равномерные, так и неравномерные) и в диапазоне чисел Рейнольдса от 100 до 15 000 не превысила величины 103. Следовательно, для нахождения решения таких СЛАУ не обязательно использовать мощные современные методы, вполне подойдёт, например, метод последовательной верхней релаксации [5].
Иная ситуация складывается со СЛАУ, которые возникают при разностной дискретизации задачи Неймана для поправки давления на базе уравнения Пуассона (10). В силу того, что, строго говоря, такая задача не имеет единственного решения, числа обусловленности матриц соответствующих СЛАУ равны бесконечности. В литературе описаны приёмы выделения единственного решения, например путем фиксации давления в какой-либо точке [7,8] или использования интегрального балансового соотношения с последующей корректировкой искомой величины на границах области и/или источникового члена в самом уравнении [8, 9]. Однако они слабо меняют ситуацию с точки зрения практических расчётов - даже после их применения процедура оценки чисел обусловленности всё равно заканчивается переполнением порядка. В настоящей работе в силу замкнутости области течения использован прием фиксации давления (обнуления поправки давления) в точке на границе области.
С другой стороны естественной характерной чертой трёхшаговой схемы расщепления является её поточечная сходимость, что не внушает оптимизма по части затрат машинного времени на достижение заданной точности решения. Результаты работ [10-12], а также настоящего исследования говорят о том, что при удвоении сеточного разрешения вдоль каждого координатного направления (учетвере-ния общего количества сеточных узлов) или увеличении числа Яе на порядок затраты машинного времени возрастают приблизительно на порядок.
Подытоживая вышесказанное, можно сделать вывод о том, что подобным трёхшаговым схемам численного решения задач о динамике несжимаемого течения вязкой жидкости присущи две проблемы: 1) поточечная сходимость общего алгоритма и 2) трудности решения СЛАУ, порожденных задачей Неймана, для поправки давления. И если первая проблема снимается переходом на другие методы (например, метод искусственной сжимаемости), что впрочем, тут же порождает проблемы, свойственные этим методам, то вторая проблема может быть преодолена путём использования мощного итерационного метода, слабо зависящего от разрешения и градиентности разностных сеток. Как следствие, целью настоящей работы является демонстрация возможностей неявного итерационного поли-нейного рекуррентного метода решения систем разностных эллиптических уравнений, способного с успехом преодолеть вторую проблему трёхшаговой схемы расщепления.
Нетрудно заметить, что рассматриваемый алгоритм решения задачи является «многоитерационным». Действительно, во-первых, имеют место итерации по слоям времени, поскольку методом установления ищется стационарное решение. При этом применяется следующее условие сходимости:
7*
шу V*
<е. (13)
Во-вторых, на каждом слое по времени решаются три системы линейных уравнений, получаемые путём разностной дискретизации уравнений движения (8), (9) и уравнения Пуассона (10). Системы эти решаются, как правило, неявным итерационным полинейным рекуррентным методом со вторым порядком экстраполяции, ускоренным в подпространствах Крылова, так называемым ЬЯ^К [13]. Понятно, что также может использоваться какой-либо иной метод решения СЛАУ, что при необходимости и делается в настоящем исследовании. В качестве критерия сходимости решения системы линейных уравнений используется соотношение
||ф* -фк-1|[ <фк|[, (14)
где фк - приближение вектора решения (и*, V* или р'), - наперед заданная точность. В силу выполнения соотношения ||У||1 = 1 и требования не менее 6-и верных знаков для элементов векторов компонент скорости (в этом случае при сеточном шаге к ~ 10-3 завихренность ю определяется как минимум до трех верных знаков, а функция тока у - до шести) точность с запасом полагается | = 10-7 - 10-8 при решении СЛАУ для и* и V*. Соответственно при решении системы уравнений для р' принимается = 10-8 - 10-9 в силу того, что 5р ~ (к/х)5^. Здесь 5р - характерная поэлементная ошибка вектора поправки давления, 8у - характерная поэлементная ошибка вектора компоненты скорости, при этом, как правило, число Куранта выбирается порядка ~ 10. При этом, как показали последующие расчеты, величину для р' можно повысить до 10-5 - 10-6 без ущерба к точности конечного результата.
Значение е в (13) следует из цепочки оценок:
10-5-10-6 (15)
к к2
при ожидаемых значениях у ~ 10-11 - 10-13 (литературные данные) и к ~ 10-3.
Следует подчеркнуть, что по причине использования метода установления промежуточные поля решения не удовлетворяют исходным уравнениям задачи. И только на момент установления полей решения и, V, р с заданной точностью величина р' и, как следствие, Vp' (в силу естественного допущения о гладкости поля р') становятся малыми, откуда следует соленоидальность поля V*. Условие (13) использовано по причине простоты оценки величины е, хотя в качестве критерия установления решения можно также воспользоваться и условием для |р'||: с требованием по точности в виде комплекса (ек/т). Дополнительно в расчетах необходимо еще следить за величинами ||и" - ии-1||:/т и - Vй-1 ||1/т, которые на момент установления в рамках условия (15) не должны превышать значений ~ 5х10-5. Данное значение получено в ходе вычислительных экспериментов, которые показали, что во всем рассматриваемом диапазоне чисел Рейнольдса и сеточного разбиения расчётной области выполнение вышеперечисленных условий является на-
дежным критерием установления решения. При их выполнении дальнейшее итерирование не приводит к сколько-нибудь заметному изменению решения.
Все расчеты выполнены на компьютере PC Intel Core i5-750 (четыре ядра), 2,66 GHz, RAM 12Gb.
Проверка правильности алгоритма решения задачи, как обычно, проводилась путём сравнения с результатами других исследователей. На рис. 2 представлены профили горизонтальной и вертикальной компонент скорости вдоль линий, параллельных стенкам каверны и проходящих через её центр, в широком диапазоне чисел Re. А в табл. 1 приведены параметры вихрей для Re = 1000. В обоих случаях видно хорошее согласование результатов, что говорит о достоверности получаемых в настоящем исследовании решений.
u
0.8
Re = 100
• Ghia, Ghia, Shin 129x129 + Каштанова, Окулова 150x150 -Настоящая работа 65x65
♦ Ghia, Ghia, Shin 129x129 ^ Каштанова, Окулова 150x150 -Настоящая работа 257x257
- Ghia, Ghia, Shin 257x257
- Erturk, Gokgol 257x257 -Настоящая работа 1025x 1025
-0.4
/ ^Ч2 Ъ
к
-0.2
-0.4
-0.6
Ghia, Ghia, Shin 129x129 Каштанова, Окулова 150x150
- Настоящая работа 65x65 Ghia, Ghia, Shin 129x129 Каштанова, Окулова 150x150
- Настоящая работа 257x257 Ghia, Ghia, Shin 257x257 Erturk, Gokpol 257x257
- Настоящая работа 1025x 1025
Re =100 '
Re =1000
Re = 10 000
0.0
0.2
0.4
0.8
Рис. 2. Профили горизонтальной (а) и вертикальной (Ь) компонент скорости вдоль линий, проходящих через центр каверны: 1 - Яе = 100, 2 - Яе = 1000, 3 - Яе = 10 000; сетки равномерные
Таблица 1
Параметры вихрей при Ке = 1000
Вихрь Параметры ОЫа и., ОЫа ЕЖ, БЫл С.Т. 257x257 И Во1е1к О., Реуге11 Я. 160x160 О Исаев В.И., Шапеев В.П. М3 Ег1игк Е., Оок5б1 С. 256x256 И Настоящая работа 1025x1025
СТ X 0,5313 0,5308 0,5308 0,5300 0,5315
У 0,5625 0,5652 0,5652 0,5650 0,5645
V -0,1179 -0,1189 -0,1189 -0,1189 -0,1186
О( -2,0497 -2,0678 - -2,0678 -2,0590
БЬ1 X 0,0859 0,0833 0,0833 0,0833 0,0827
У 0,0781 0,0781 0,0781 0,0783 0,0782
V 2,31-10-4 2,33-10-4 2,33-10-4 2,33-10-4 2,33-10-4
О( 0,3618 0,3523 - 0,3543 0,3494
ВЬ2 X - 0,0048 0,0048 0,0050 0,0049
У - 0,0048 0,0048 0,0050 0,0049
V - -6,40-10-9 -6,40-10-9 -6,52-10-9 -7,29-10-9
О( - - - -2,96-10-3 -2,79-10"3
БЬ3 X - - 2,71-10-4 - 2,55-10-4
У - - 3,18-10-4 - 2,52-10-4
V - - 2,1910-13 - 6,55-10"14
О( - - - - 1,22-10-5
БЯ1 X 0,8594 0,8640 0,8640 0,8633 0,8644
У 0,1094 0,1118 0,1118 0,1117 0,1118
V 1,75-10-3 1,73-10"3 1,73-10-3 1,73-10-3 1,73-10-3
О( 1,1547 1,1098 - 1,1182 1,1057
ВЯ2 X 0,9922 0,9923 0,9923 0,9917 0,9923
У 0,0078 0,0077 0,0077 0,0067 0,0076
V -9,32-10-8 -5,04-10-8 -5,04-10-8 -4,92-10-8 -5,39-10-8
О( -8,53-10-3 - - -7,69-10"3 -8,41-10-3
БЯ3 X - - 0,99950 - 0,99965
У - - 4,56-10-4 - 3,54-10-4
V - - 1,44-10-12 - 1,2110-13
О( - - - - 3,75-10-5
При повышении числа Яе до 10 000 с использованием неравномерной сетки 1025x1025 с параметром сжатия к0 = 10 в нижних углах каверны удается выявить вихри третьего и даже четвертого уровней. С другой стороны, существует аналитическое решение задачи о медленном (ползущем) движении вязкой жидкости в углу, описываемом уравнением Стокса [14]. Схема такого течения, представляющая собой бесконечную цепочку так называемых вихрей Моффатта, приведена на рис. 3. В этой ситуации интересно сравнить результаты численного и аналитического решений. Аналитическое решение Моффатта позволяет рассчитать относительные величины размеров вихрей г(5) = р(5)/р(0) и так называемых «интенсивно-стей» вихрей = Р^-1/^0-1. Здесь р - характерный размер вихря, V - величина линейной скорости циркуляции в точке стыка вихрей, 5 - номер вихря, нарастающий по мере «заглубления» в угол (5 = 0 - базовый вихрь).
В качестве базовых были выбраны вихри второго уровня, поскольку их центры, по сравнению с вихрями первого уровня, оказались слабо сдвинутыми относительно биссектрис углов из-за незначительного влияния конвективного переноса, что соответствует выраженному ползущему характеру течения. Параметры нижних вихрей второго - четвертого уровней, измеренные непосредственно по рассчитанным полям компонент скорости потока, имеют следующие значения (табл. 2).
Таблица 2
Параметры bl2 BL3 BL4 br2 BR3 BR4
s 0 1 2 0 1 2
x 0,0212 0,0014 1,5-10-4 0,9353 0,9957 0,9997
y 0,0260 0,0014 1,5-10-4 0,0757 0,0044 3,0-10-4
dx 0,0679 0,0032 1,0-10-4 0,1600 0,0104 6,0-10-4
dy 0,0493 0,0032 1,0-10-4 0,1525 0,0105 6,0-10-4
p(, 0,0586 0,0032 1,0-10-4 0,1575 0,0105 6,0-10-4
yis) 6,26-10-4 2,74-10-7 9,20-10-10 1,83-10-2 8,99-10-6 6,20-10-9
Здесь х, у координаты центров вихрей, ёх и ёу соответственно горизонтальный и вертикальный размеры вихрей, определяется в точке А (рис. 3), 1*Г) - в точке В, У(2) - в точке С.
На основании данных параметров вихрей нетрудно вычислить значения г(х) и и сравнить их с решением из [14] (см. табл. 3).
Таблица 3
s r(s)
[141 BL BR [141 BL BR
1 0,06140 0,05460 0,06600 4,86-10-4 4,38-10-4 4,93-10-4
2 0,00377 0,00170 0,00380 2,36-10-7 1,47-10-6 3,40-10-7
Учитывая приближенность сеточного решения и особенно оценочный характер значений параметров вихрей четвертого уровня, полученное соответствие численного и аналитического решений (особенно для правых вихрей) можно считать очень хорошим.
Эффективность применения полинейного рекуррентного метода при численном моделировании динамики несжимаемой вязкой жидкости демонстрируется на рис. 4, а, где изображены зависимости ||&уУ*||: и кр от номера слоя по времени п.
Рис. 4. Сравнение использования методов ЬЯ28К (а) и БьСО81аЪ РБ (Ь) для решения СЛАУ на базе задачи Неймана для поправки давления р'.
Здесь-||(Иу У*||ь---кр; 1 - Яе = 1000, сетка 257x257 равномерная,
2 - Яе = 10 000, сетка 1025x1025 неравномерная (а) или равномерная (Ь)
Здесь кр - количество итераций, необходимых для сходимости метода решения СЛАУ для р'. Решены простая задача: Яе = 1000, сетка 257x257 равномерная, и сложная: Яе = 10 000, сетка 1025x1025 неравномерная, ка = 10. Во всех расчетах е = 10-5. Для сравнения на рис. 4, Ь приведены аналогичные зависимости при использовании метода бисопряженных градиентов со стабилизацией [15] с предобу-славливателем на базе явного метода Булеева. В этом случае пришлось упростить сложную задачу применением равномерной сетки, в противном случае поиск решения грозил затянуться в лучшем случае на многие сотни часов. Затраты машинного времени составили на решение простой задачи с применением ЬК2бК. -149 с, с применением В1-Св81аЬ РВ - 421 с; на решение сложной задачи с применением ЬЯ^К - 47 162 с, с применением В1-Св81аЬ РВ - 680 687 с.
Хорошо видно, что независимо от используемого метода и сложности задачи имеет место слабая корреляция кривых ||&у V ||1 и кр между собой. Также следует отметить, что значение кр не сильно откликается на сложность задачи, особенно на начальном этапе процесса установления.
Основной результат сравнения заключается в том, что ЬК2бК. тратит на решение СЛАУ 10-20 итераций на каждом слое по времени, в то время как В1-Св81аЬ РВ - 100-400 итераций. И хотя одна итерация ЬК2бК по длительности исполнения соответствует 2-6 итерациям В1-Св81аЬ РВ (в зависимости от конфигурации вычислительной системы), преимущество применения ЬК2бК бесспорно. Этот же вывод подтверждается и данными табл. 4, где приведены времена (сек.), затраченные на решение задачи при Яе = 1000 для различных сеточных разрешений с постоянным шагом.
Таблица 4
Время (с), затраченное на решение задачи при Ке = 1000 (е = 10-5, шаг постоянный)
Bi-CGStab PLU Bi-CGStab PB LR1 LR2 LR1sK LR2sK
65x65 6,0 5,0 3,0 4,0 3,0 3,0
129x129 42,0 32,0 18,0 21,0 21,0 16,0
257x257 507,0 421,0 157,0 218,0 157,0 149,0
513x513 9 766,0 8 898,0 1 952,0 3 087,0 1 705,0 1 903,8
1025x1025 153 144,0 143 076,0 19 021,0 35 678,0 15 732,0 23 396,0
Здесь Bi-CGStab PLU - метод бисопряженных градиентов со стабилизацией с предобуславливателем на базе неполной LU-факторизации исходной матрицы системы линейных уравнений, соответственно Bi-CGStab PB - с предобуславливателем на базе явного метода Булеева, LR1 - полинейный рекуррентный метод с первым порядком экстраполяции приращения решения, LR2 - со вторым порядком экстраполяции приращения решения, соответственно LR1sK и LR2sK - это LR1 и LR2, ускоренные в подпространствах Крылова [13]. Как уже отмечалось ранее, имеет место катастрофическое нарастание затрат машинного времени при каждом учетверении количества расчётных узлов.
Кроме сеточного разрешения дополнительными усложняющими процесс вычислений факторами являются степень неравномерности сетки и величина числа Рейнольдса. Их действие выражено не так явственно, как в случае увеличения сеточного разрешения, в силу того, что невозможно для каждого расчетного вари-
анта подобрать оптимальные значения числа Куранта и параметра компенсации для ьк2бк, тем самым поставив все расчеты в одинаковые сравнимые условия. Однако однозначно можно говорить о тенденции замедления установления решения при переходе на неравномерную сетку и увеличении числа Яе. Например, при увеличении Яе с 1000 до 10 000 (сетка 257x257 равномерная, е = 10-5, ЬЯ^К) количество шагов по времени увеличилось с 505 до 2 077, а общие затраты машинного времени - со 153 до 502 сек. Аналогично, при переходе с равномерной сетки на неравномерную с параметром сжатия ко = 10 (Яе = 10 000, сетка 513x513, е = 10-5, ьк2бк) количество шагов по времени увеличилось с 8 521 до 9 442, а общие затраты машинного времени - с 13 497 до 16 764 сек.
В качестве иллюстрации возможностей применения метода была решена задача при Яе = 15 000 на сетке 513x513 в двух вариантах: равномерной и неравномерной с параметром сжатия ко = 10. На первом этапе алгоритма расщепления для решения СЛАУ относительно и* и V* использован метод ЬЯ2. При точности Е, = 10-8 метод сходится за 2-3 итерации. Для решения СЛАУ относительно р' использован ЬЯ^К. При точности = 10-9 метод сходится за 15-25 итераций в среднем.
Рис. 5. Общая картина течения при Яе = 15 000, сетка 513x513, неравномерная (ко = 10)
На рис. 5 представлена общая картина течения при Яе = 15 000. Кроме центрального вихря хорошо видны верхние левые вихри первого (ТЬ1) и второго (ТЬ2) уровней, а также нижние левые и правые вихри первых двух уровней БЬ1, БЬ2 и БЯ1, БЯ2 соответственно. Увеличение масштаба представления приблизи-
тельно в 100 раз позволяет выявить в нижних левом и правом углах вихри третьего и четвёртого уровней БЬ3, БЬ4 и БЯ3, БЯ4 соответственно (рис. 6).
Рис. 6. Вихри третьего и четвертого уровней в нижних углах каверны при Яе = 15 000, сетка 513x513, неравномерная (кв = 10)
В табл. 5 приведены характеристики вихрей.
Таблица 5
Параметры вихрей при Re = 15 000 (е = 5-10 6, 513x513 G)
Вихри dx dy ю У Qx Qy
CT 0,5177 0,5247 1,0 1,0 -0,9420 -0,0791 -0,0791 -0,0791
TLi 0,0677 0,9145 0,1497 0,3247 2,2018 2,35-10-3 2,36-10"3 2,36-10"3
BLi 0,0512 0,1902 0,3727 0,3051 4,0107 2,52-10-3 2,52-10"3 2,52-10"3
BR1 0,8098 0,0540 0,2946 0,3983 7,6080 4,93-10-3 4,90-10-3 4,90-10-3
tl2 0,0055 0,8391 0,0081 0,0311 -0,1696 -7,40-10-7 -6,90-10-7 -6,80-10-7
bl2 0,0736 0,0600 0,1468 0,1400 -0,3886 -1,78-10-4 -1,78-10-4 -1,78-10-4
br2 0,9424 0,1134 0,1534 0,1966 -1,0988 -6,68-10-4 -6,68-10"4 -6,68-10"4
BL3 0,0042 0,0042 0,0108 0,0107 3,19-10-3 5,42-10"9 5,18-10-9 5,18-10-9
BR3 0,9925 0,0078 0,0188 0,0206 9,08-10-3 5,44-10-8 5,37-10-8 9,95-10-8
BL4 2,27-10-4 2,25-10-4 0,0006 0,0005 -1,89-10-5 -3,22-10-14 -3,0-10"13 -3,0-10"13
BR4 0,9996 0,0004 0,0006 0,0006 -4,33-10-5 -2,50-10"13 -1,80-10-12 3,60-10"10
Как и ранее, когда сравнивалось решение при Яе = 10 000 с аналитическим решением [14], здесь следует обратить внимание на положение центров нижних вихрей. Для первого уровня они значительно удалены от биссектрис углов - сказывается сильное влияние конвективного переноса. По мере увеличения номера уровня вихря это влияние резко сходит на нет, движение жидкости становится фактически Стоксовым, то есть определяется только силами вязкости и интенсивностью внешнего течения. Соответственно частично со второго и уже полностью с третьего уровней центры вихрей располагаются на биссектрисах углов. Об этом же говорят и оценки локальных чисел Рейнольдса Яел = 15 000-¥-ё (табл. 6), построенных по характерным безразмерным параметрам вихрей: размеру ё и скорости циркуляции V.
Таблица 6
Оценка локальных чисел Рейнольдса, построенных по характерным параметрам вихрей (Ке = 15 000, б = 5-10-6, 513x513 С)
Вихрь BLi BR1 bl2 br2 BL3 BR3 BL4 BR4
ReA « 250 520 5 20 10-4 10-3 10-8 10-7
Хорошей внутренней проверкой полученных результатов служит определение экстремумов функции тока в центрах вихрей через расходные характеристики потока. Действительно, пусть xc, yc - координаты центра вихря, тогда из определения функции тока у и граничного условия на стенке каверны у = 0 следует, что
yc xc
У( xc, Ус ) = Qx = { u(xc, y)^y = -Qy =-{ v(x yc )dx
0 0
Хорошее совпадение у (xc, yc) с Qx и Qy в табл. 5 говорит о точности решения задачи (6), (7), включая определения поля завихренности ю. При этом полный расход через произвольное сечение каверны составляет величину порядка 10-8 -10-6 в зависимости от местоположения сечения, что объясняет несовпадение У (xc, yc) с Qx и Qy для вихрей четвертого уровня, характерные размеры которых сопоставимы с шагом сетки.
На основании изложенного материала можно сделать вывод о высокой эффективности неявного итерационного полинейного рекуррентного метода при численном моделировании динамики несжимаемой вязкой жидкости. Его использование позволило уменьшить остроту одной из двух проблем трёхшаговой схемы расщепления - решение задачи Неймана на базе уравнения Пуассона для поправки давления p'.
ЛИТЕРАТУРА
1. Ghia U., Ghia K.N., Shin C.T. High-Re solution for incompressible flow using the Navier-Stokes equations and a multigrid method // J. Computational Physics. 1982. V. 48. P. 387411.
2. Barragy E., Carey G.F. Stream function-vorticity driven cavity solution using p finite elements // Computers & Fluids. 1997. V. 26. No. 5. P. 453-468.
3. ErturkE., Corke T. C., Gôkçôl C. Numerical solutions of 2-D steady incompressible driven cavity flow at high Reynolds numbers // Int. J. Numer. Meth. Fluids. 2005. V. 48. P. 747-774.
4. Исаев В.И., Шапеев В.П. Варианты метода коллакации и наименьших квадратов повышенной точности для численного решения уравнения Навье - Стокса // ЖВМ и МФ. 2010. Т. 50. № 10. C. 1758-1770.
5. Роуч П. Вычислительная гидродинамика. М.: Мир. 1980. 616 c.
6. Белоцерковский О.М., Гущин В.А., Щенников В.В. Метод расщепления в применении к решению задач динамики вязкой несжимаемой жидкости // ЖВМ и МФ. 1975. Т. 15. № 1. C. 197-207.
7. Патанкар С. Численные методы решения задач теплообмена и динамики жидкости: пер. с англ. М.: Энергоатомиздат, 1984. 152 c.
8. КузнецовА.Е., Стрелец М.Х. Численное моделирование существенно дозвуковых стационарных неизотермических течений однородного вязкого газа в каналах // Численные методы механики сплошной среды. Новосибирск: ВЦ СО АН СССР, 1983. Т. 14. № 6. C. 97-114.
9. Фомин А.А. Численное исследование влияния граничных условий на решение задач термогравитационной конвекции в открытых областях. Томск: ТГУ, 1985. 51 с. Деп. в ВИНИТИ 22.11.85 № 8069-В.
10. Деги Д.В, Старченко А.В. Численное решение уравнений Навье - Стокса на компьютерах с параллельной архитектурой // Вестник Томского государственного университета. Математика и механика. 2012. № 2. C. 88-98.
11. Wright N.G. Multigrid solutions of elliptic fluid flow problems. Ph. D. thesis. University of Leeds, 1988. 185 p.
12. Каштанова С.В., Окулова Н.Н. Математическое моделирование течения вязкой теплопроводной жидкости с использованием метода Ls-Stag // Вестник МГТУ им. Н. Э. Баумана. Сер. Естественные науки. 2012. Спец. вып. 2: Математическое моделирование в технике. С. 86-97.
13. Фомин А.А., Фомина Л.Н. Неявный итерационный полинейный рекуррентный метод решения разностных эллиптических уравнений. Кемерово: КемГУ, 2012. 314 с.
14. Moffatt H.K. Viscous and resistive eddies near a sharp corner // J. Fluid Mechanics. 1964. V. 18, part 1. P. 1-18.
15. Van der VorstH.A. Bi-CGStab: a fast and smoothly converging variant of BI-CG for solution of nonsymmetric linear systems // SIAM J. Sci. Stat. Comput. 1992. V. 13. No. 2. P. 631-644.
Статья поступила 23.12.2013 г.
Fomin A.A., Fomina L.N. NUMERICAL SOLUTION OF THE NAVIER - STOKES EQUATIONS IN THE MODELING OF TWO-DIMENSIONAL VISCOUS INCOMPRESSIBLE FLUID FLOWS. In this paper, the effectiveness of the implicit iterative line-by-line recurrence method for solving difference elliptical equations arising in numerical simulations of two-dimensional viscous incompressible fluid flows is analyzed. The research is carried out by an example of the problem about a steady two-dimensional lid-driven cavity flow formulated in primitive variables (u, v, p). It is shown that applying the line-by-line recurrence method allows one to reduce the total time for solving the problem in comparison with the use of the present-day effective bi-conjugate gradients method with stabilization. As an illustration of the achieved results, the structure of the flow at Re = 15000 is shown. Here, in terms of the use of a non-uniform grid, it became possible to obtain a sequence of bottom-corner vortices up to the fourth level. As a validation of the received solution, the comparison of basic parameters of all vortices with results of other authors was carried out at Re = 1000. In addition, the mass imbalance was estimated; it did not exceed 10-8^10-6 depending on the location of the cross section in the cavity, and a comparison of the relative size and 'intensity' of bottom-corner vortices of the third and fourth levels with the Moffatt analytical solution of the problem of a viscous fluid flow near a sharp corner was carried out.
Keywords: lid-driven cavity flow, Navier-Stokes equations, implicit iterative line-by-line recurrence method.
FOMIN Alexander Arkad'evich (Candidate of Physics and Mathematics, T. F. Gorbachev Kuzbass State Technical University, Kemerovo, Russian Federation)
E-mail: fomin_aa@mail.ru
FOMINA Lubov Nikolaevna (Candidate of Physics and Mathematics, Assoc. Prof., Kemerovo State University, Kemerovo, Russian Federation)
E-mail: lubafomina@mail.ru
REFERENCES
1. Ghia U., Ghia K. N., Shin C. V. High-Re solution for incompressible flow using the Navier-Stokes equations and a multigrid method (1982) Journal of Computational Physics. V. 48, pp. 387-411.
2. Barragy E., Carey G.F. Stream function-vorticity driven cavity solution using p finite elements (1997) Computers & Fluids. V. 26. No. 5, pp. 453-468.
3. Erturk E., Corke V.C., Gôkçôl C. Numerical solutions of 2-D steady incompressible driven cavity flow at high Reynolds numbers (2005) Int. J. Numer. Meth. Fluids. V. 48, pp. 747-774.
4. Isaev V.I., Shapeev V.P. Varianty metoda kollakatsii i naimen'shikh kvadratov povyshennoy tochnosti dlya chislennogo resheniya uravneniya Nav'e - Stoksa (2010.) ZhVM i MF. V. 50. No. 10, pp. 1758-1770. (in Russian)
5. Rouch P. Vychislitel'naya gidrodinamika. Moskow, Mir Publ., 1980. 616 p.
6. Belotserkovskiy O.M., Gushchin V.A., Shchennikov V.V. Metod rasshchepleniya v primenenii k resheniyu zadach dinamiki vyazkoy neszhimaemoy zhidkosti (1975) ZhVM i MF. V. 15. No 1, pp. 197-207. (in Russian)
7. Patankar S. Chislennye metody resheniya zadach teploobmena i dinamiki zhidkosti. Per. s angl. Moskow, Energoatomizdat Publ., 1984. 152 p. (in Russian)
8. Kuznetsov A.E., Strelets M.Kh. Chislennoe modelirovanie sushchestvenno dozvukovykh statsionarnykh neizotermicheskikh techeniy odnorodnogo vyazkogo gaza v kanalakh (1983) Chislennye metody mekhaniki sploshnoy sredy. Novosibirsk, VTs SO AN SSSR Publ. V. 14. No. 6, pp. 97-114. (in Russian)
9. Fomin A.A. Chislennoe issledovanie vliyaniya granichnykh usloviy na reshenie zadach termogravitatsionnoy konvektsii v otkrytykh oblastyakh. Tomsk, TGU, 1985. 51 s. Dep. v VINITI 22.11.85 No 8069-B. (in Russian)
10. Degi D.V., Starchenko A.V. Chislennoe reshenie uravneniy Nav'e-Stoksa na komp'yuterakh s parallel'noy arkhitekturoy (2012) Vestnik Tomskogo gosudarstvennogo universiteta. Matematika i mekhanika. No. 2, pp. 88-98. (in Russian)
11. Wright N.G. Multigrid solutions of elliptic fluid flow problems. Ph. D. thesis. University of Leeds, 1988. 185 p.
12. Kashtanova S.V., Okulova N.N. Matematicheskoe modelirovanie techeniya vyazkoy teploprovodnoy zhidkosti s ispol'zovaniem metoda Ls-Stag (2012) Vestnik MGTU im. N.E. Baumana. Ser. Estestvennye nauki. Spets. vyp. 2: Matematicheskoe modelirovanie v tekhnike, pp. 86-97. (in Russian)
13. Fomin A.A., Fomina L.N. Neyavnyy iteratsionnyy polineynyy rekurrentnyy metod resheniya raznost-nykh ellipticheskikh uravneniy. Kemerovo: KemGU Publ., 2012. 314 p. (in Russian)
14. Moffatt H.K. Viscous and resistive eddies near a sharp corner (1964) Journal of Fluid Mechanics. V. 18, part 1, pp. 1-18.
15. Van der Vorst H.A. Bi-CGStab: a fast and smoothly converging variant of BI-CG for solution of nonsymmetric linear systems (1992) SIAM J. Sci. Stat. Comput. V. 13. No. 2, pp. 631644.
