Numerical solution of the Navier - Stokes equations in the modeling of two-dimensional viscous incompressible fluid flows.pdf Задача о течении вязкой несжимаемой жидкости в квадратной каверне с подвижной крышкой является своеобразным «полигоном», на котором вот уже несколько десятков лет многочисленные исследователи демонстрируют свои достижения в области численного решения уравнений Навье - Стокса. Привлекательность этой задачи объясняется, с одной стороны, простотой её формулировки: решение зависит фактически от одного параметра - числа Рейнольдса. А с другой - наличием разрывов и сингулярностей полей решения в верхних углах каверны и фактически бесконечные цепочки вихрей в каждом из нижних углов гарантируют возможность постоянного улучшения решения и, как следствие, постоянный интерес исследователей к этой задаче. Накопленные за последнее время обширные материалы позволяют адекватно оценивать степень правильности новых результатов, обеспечивая тем самым их высокую достоверность. Среди работ на эту тему до сих пор большой популярностью заслуженно пользуется работа [1], сформировавшая современный подход к изложению материалов исследования. Несмотря на использование умеренных по современным меркам равномерных сеток типа 129x129 и 257x257, авторам удалось построить решение задачи по числам Рейнольдса вплоть до Re = 10 000 и оценить параметры вихрей вплоть до правого нижнего вихря третьего уровня. В абсолютном большинстве работ исследуются характеристики течения при Re < 10 000. Более того, ряд авторов высказали предположение, что не существует стационарного решения задачи при больших значениях числа Рейнольдса. Однако некоторые публикации последнего времени опровергают это предположение. В частности, в работе [2] приведены параметры вихрей при Re = 12 500 вплоть до четвёртого уровня (так называемые BL4 и BR4) в нижних углах каверны, но при дальнейшем увеличении числа Re отмечается неустойчивость решения. А в работе [3] (и ряде других публикаций тех же авторов) удалось поднять число Рей-нольдса до 21 000, но при этом, в силу равномерности разностных сеток, не удалось выявить вихри четвёртого уровня. При этом там же утверждается, что эффективным средством борьбы против неустойчивости являются более подробные сетки. В частности авторами были использованы сетки 257x257 вплоть до Re = 10 000, при 10 000 < Re < 15 000 - сетки 513x513, при Re > 15 000 - сетки 601x601 и подробнее вплоть до 1025x1025. Следует отметить, что приблизительно половина исследователей предпочитают постановку задачи в переменных у-ю. Возможно, это связано с тем, что основное сравнение результатов осуществляется по значениям экстремумов полей у-ю (центров вихрей), в то время как использования естественных переменных (и, v, p), несмотря на удобство и абсолютную точность формулировки граничных условий, требует численного восстановления полей ю и у, что не может не сказаться на точности представления сравниваемых параметров течения. Среди последних работ в переменных (и, v, p) следует отметить [4] и другие публикации тех же авторов. Использование высокого порядка аппроксимации и кусочно-равномерной сетки - свыше 65 тысяч расчетных узлов (так называемой сетки М3) - позволило получить высокоточные результаты при Re = 1000 по параметрам вихрей вплоть до третьего уровня. На базе литературных данных можно составить обобщенную схему стационарного течения (см. рис. 1). В ее основу положено решение конкретного варианта задачи при Re > 10 000. При этом необходимо заметить, что вихри третьего и четвертого уровней в нижних углах изображены вне масштаба - в расчетах они получаются значительно меньших размеров. BL4 BR4 X Рис. 1. Система координат и обобщенная схема течения В настоящем исследовании математическая постановка задачи в безразмерном виде представляет собой систему нестационарных двумерных уравнений Навье -Стокса: ut + (V-V)и = -Vxp + Re-1 V2и ; (1) V^ + (V -V) v = -V p + Re-1 V2v ; (2) (3) VV = 0, а также начальных u = v = 0 и краевых условий: на верхней границе области 0 < x < 1, y = 1: u = 1, v = 0, (4) и на остальных границах u = v = 0. (5) Здесь t - время, x, у - пространственные координаты, u, v - x, y компоненты вектора скорости потока V, p - давление, Re = UL/u - число Рейнольдса, и - кинематический коэффициент вязкости жидкости. В начальные моменты времени 0 < t < 1 скорость движения крышки плавно увеличивается до 1 и далее не меняется. Для построения линий тока ищется решение вспомогательной задачи относительно функции тока у. Внутри области решения Ду = -ю, (6) на границах области решения у = 0. (7) Здесь ю = dv/dx - du/dy - завихренность потока. Задача (1)-(5) решается численно. Область решения покрывается так называемой МАС-сеткой [5], в которой узлы определения давления располагаются внутри прямоугольной сеточной ячейки, а узлы определения компонент скорости - на её гранях. В общем случае сетка слабонеравномерная, то есть \к, - hi-1 | ~ O(h2), где h - сеточный шаг. Степень неравномерности определяется параметром сжатия kG, равным отношению среднего шага сетки к минимальному. На каждом шаге по времени реализуется классическая трёхшаговая схема расщепления [6], на первом шаге которой определяются предварительные значения компонент скорости: u * = un-1 -x[(V"-1 -V) u * -Re-1 V 2u * +Vxpn-1 ] ; (8) v * = vn-1 -t[(Vn-1 -V) v * -Re-1 V2v * +Vypn-1 ] . (9) На втором шаге определяется поле поправки давления p' путём решения уравнения Пуассона Дp'=x-1V V*. (10) И наконец, на третьем шаге поля скорости и давления корректируются по формулам Vn = V* -xVp', pn = pn-1 +ap' . (11) Здесь т - шаг по времени, n - номер слоя по времени, ст - итерационный параметр. Граничные условия для (8), (9) берутся по аналогии с (4), (5) с учетом естественной замены u ^ u*, v ^ v*; а для (10) - условия Неймана, которые следуют из (11) в силу стремления V* ^ Vn при n ^ да. Во всех расчётных вариантах методом установления ищется стационарное решение задачи (1)-(5). Путём первоначального численного эксперимента было определено оптимальное значение ст = 1,8, которое в дальнейшем не менялось. Проверка полной аппроксимации уравнения движения Vn = Vn-1 - т [(Vn-1 - V) (Vn + TVp') - Re-1 V2 (Vn + TVp') - CTVp' + Vpn ] = = Vn-1 -t[(Vn-1 -V)Vn -Re-1V2Vn +Vpn -ctVp']-t2 [(Vn-1 -V)Vp'-Re-1V2(Vp')] = = Vn-1 -t[(Vn-1 -V)Vn -Re-1V2Vn +Vpn] + O(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], а также настоящего исследования говорят о том, что при удвоении сеточного разрешения вдоль каждого координатного направления (учетвере-ния общего количества сеточных узлов) или увеличении числа Re на порядок затраты машинного времени возрастают приблизительно на порядок. Подытоживая вышесказанное, можно сделать вывод о том, что подобным трёхшаговым схемам численного решения задач о динамике несжимаемого течения вязкой жидкости присущи две проблемы: 1) поточечная сходимость общего алгоритма и 2) трудности решения СЛАУ, порожденных задачей Неймана, для поправки давления. И если первая проблема снимается переходом на другие методы (например, метод искусственной сжимаемости), что впрочем, тут же порождает проблемы, свойственные этим методам, то вторая проблема может быть преодолена путём использования мощного итерационного метода, слабо зависящего от разрешения и градиентности разностных сеток. Как следствие, целью настоящей работы является демонстрация возможностей неявного итерационного поли-нейного рекуррентного метода решения систем разностных эллиптических уравнений, способного с успехом преодолеть вторую проблему трёхшаговой схемы расщепления. Нетрудно заметить, что рассматриваемый алгоритм решения задачи является «многоитерационным». Действительно, во-первых, имеют место итерации по слоям времени, поскольку методом установления ищется стационарное решение. При этом применяется следующее условие сходимости:

Скачать электронную версию публикации