БЕЗДИВЕРГЕНТНЫЙ МЕТОД КОЛЛОКАЦИЙ И НАИМЕНЬШИХ КВАДРАТОВ ДЛЯ РАСЧЕТА ТЕЧЕНИЙ НЕСЖИМАЕМОЙ ЖИДКОСТИ И ЕГО ЭФФЕКТИВНАЯ РЕАЛИЗАЦИЯ Текст научной статьи по специальности «Математика»

Вестн. Сам. гос. техн. ун-та. Сер. Физ.-мат. науки. 2020. Т. 24, № 3. С. 542-573

ISSN: 2310-7081 (online), 1991-8615 (print) d https://doi.org/10.14498/vsgtu1758

УДК 519:63.4:532.51.5

Бездивергентный метод коллокаций и наименьших квадратов для расчета течений несжимаемой жидкости и его эффективная реализация

© Е. В. Ворожцов1, В. П. Шапеев1'2

1 Институт теоретической и прикладной механики им. С.А. Христиановича СО РАН,

Россия, 630090, Новосибирск, ул. Институтская, 4/1.

2 Новосибирский национальный исследовательский университет,

Россия, 630090, Новосибирск, ул. Пирогова, 2.

Аннотация

Рассматривается проблема ускорения итерационного процесса численного решения методом коллокаций и наименьших квадратов (КНК) краевых задач для уравнений с частными производными. Для ее решения предложено применять одновременно три способа ускорения итерационного процесса: предобуславливатель, многосеточный алгоритм и метод Крылова. Предложен метод нахождения оптимальных значений параметров двухпараметрического предобуславливателя. Использование найденного предобуславливателя существенно ускоряет итерационный процесс. Исследовано влияние на итерационный процесс всех трех способов его ускорения: каждого по отдельности, а также при их комбинированном применении. Наибольший вклад дает применение алгоритма, использующего подпространства Крылова. Комбинированное применение одновременно всех трех способов ускорения итерационного процесса решения краевых задач для двумерных уравнений Навье-Стокса уменьшило время их решения на компьютере до 362 раз по сравнению со случаем, когда применялся только один из них — предобуславливатель.

Ключевые слова: предобуславливание, подпространства Крылова, многосеточные алгоритмы, уравнения Навье-Стокса, метод коллокаций и наименьших квадратов.

Получение: 25 ноября 2019 г. / Исправление: 29 июля 2020 г. / Принятие: 24 августа 2020 г. / Публикация онлайн: 21 сентября 2020 г.

Научная статья

3 ©® Контент публикуется на условиях лицензии Creative Commons Attribution 4.0 International (https://creativecommons.org/licenses/by/4.0/deed.ru) Образец для цитирования

Ворожцов Е. В., Шапеев В. П. Бездивергентный метод коллокаций и наименьших квадратов для расчета течений несжимаемой жидкости

и его эффективная реализация // Вестн. Сам. гос. техн. ун-та. Сер. Физ.-мат. науки, 2020. Т. 24, № 3. С. 542-573. https://doi.org/10.14498/vsgtu1758. Сведения об авторах

Евгений Васильевич Ворожцов А https://orcid.org/0000- 0003- 2753- 8399 доктор физико-математических наук, профессор; ведущий научный сотрудник; лаб. физики быстропротекающих процессов; e-mail: vorozh@itam.nsc.ru Василий Павлович Шапеев © https://orcid.org/0000-0001-6761-7273 доктор физико-математических наук, профессор; главный научный сотрудник; лаб. термомеханики и прочности новых материалов; e-mail: shapeev@itam.nsc.ru

542 © Самарский государственный технический университет

Введение. Для численного решения уравнений Навье—Стокса вязкой несжимаемой жидкости к настоящему времени получили широкое распространение конечно-разностные методы [1], методы конечных элементов [2] и методы конечного объёма [3].

В течение последних трёх десятилетий за рубежом получает всё большее распространение класс численных методов, в которых дискретизация краевой или начально-краевой задачи для уравнений Навье-Стокса приводит к переопределенной системе линейных алгебраических уравнений (СЛАУ). Одним из таких методов яляется коллокационный LSFEM — «коллокационный метод конечных элементов и наименьших квадратов» [2,4], в котором метод конечных элементов (МКЭ) комбинируется с методом наименьших квадратов. Методу LSFEM присущи некоторые недостатки, в частности, нарушается закон сохранения массы, который для несжимаемой жидкости выражается уравнением div V = 0, где V — вектор скорости жидкости [5,6]. Другой недостаток метода LSFEM — трудности, возникающие при вычислении давления с использованием вектора скорости.

В методе граничных элементов (BEM) [7] также возникает переопределенная СЛАУ. Причиной этого является то, что в этом методе используются условия непрерывности решения и его производной в направлении нормали к каждому граничному элементу. Включение этих условий в СЛАУ как раз и приводит к её переопределенности.

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

Известно, что СЛАУ, полученное аппроксимацией дифференциальной задачи методом коллокации без комбинирования его с методом наименьших квадратов, даже при умеренных числах Рейнольдса имеет плохую обусловленность [8]. Поэтому актуальна задача расширения возможностей существующих методов за счет улучшения их свойств.

Метод коллокаций и наименьших квадратов (КНК) численного решения краевых задач для дифференциальных уравнений возник недавно. В настоящее время предложены и опубликованы несколько различных вариантов этого метода, которые имеют заметные преимущества перед первоначальными вариантами [9-13]. В частности, в [9] было предложено в методе КНК дополнить систему уравнений коллокаций линейными условиями согласования локального решения в каждой ячейке с локальными решениями, взятыми во всех соседних с ней ячейках.

В работе [14] было предложено вводить пять управляющих (регулируемых) параметров в условия согласования с целью уменьшения числа обусловленности. Было показано с помощью численных экспериментов, что область значений параметров, при которых глобальная СЛАУ, определяющая глобальное решение задачи в методе КНК, хорошо обусловлена, в значительной мере пересекается с областью, где наблюдается наилучшая точность численного решения задачи. В работах [13, 15] показано, что включение условий согласования в локальные СЛАУ позволяет существенно уменьшить число

обусловленности этих СЛАУ от значений 105 до значений в пределах от 3 до 10 по сравнению со случаем, когда они не используются.

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

(i) в этом методе, в отличие от LSFEM, обеспечивается точное выполнение закона сохранения массы благодаря использованию соленоидально-го базиса;

(ii) отсутствует проблема увязывания давления с вектором скорости, потому что давление вычисляется в методе КНК одновременно с составляющими вектора скорости.

Отметим также, что метод КНК позволяет эффективно решать задачи для эллиптических, параболических и гиперболических уравнений в частных производных (УЧП) на различных адаптивных сетках с прямоугольными и треугольными ячейками.

Для достижения большего ускорения решения приближенной задачи в данной работе рассматривается комбинированное применение трех способов ускорения итерационного процесса: предобуславливателя, операции продолжения на многосеточном комплексе, являющейся составной частью метода Федоренко [16], и метода Крылова [17,18].

В работе [13] для аппроксимации составляющих скорости использованы многочлены второй степени по пространственным переменным, а для аппроксимации давления — многочлены первой степени, так что общее количество базисных вектор-функций составляло 12 в представлении приближенного решения в пространстве полиномов. Здесь с целью повышения точности численного решения применялся многочлен второй степени и для аппроксимации давления. В этом случае в общей сложности имеются 15 независимых базисных функций — полиномов в выбранном пространстве. Их совокупность можно назвать соленоидальным базисом, так как каждый базисный вектор является бездивергентным. Поэтому уравнение неразрывности и, следовательно, закон сохранения массы удовлетворяется численным решением задачи во всей расчётной области. В дальнейшем метод КНК, в котором использовались 12 базисных векторов [13,15], будем называть методом КНК12, а метод КНК с 15 базисными векторами — методом КНК15.

При использовании конечно-разностных методов и методов конечного объема для численного решения задач с открытой границей, через которую жидкость может свободно вытекать из канала, существует проблема устойчивых граничных условий на открытой границе. Некорректная формулировка и/или плохая аппроксимация (реализация) принятой формулы граничного условия задачи на такой границе приводит к неустранимой погрешности, которая распространяется на всю область решения задачи. В литературе описано несколько видов записи граничных условий на открытой границе. В частности, в [19,20] применялось граничное условие Неймана J^ = 0, которое выполняется точно только в пределе для канала бесконечной длины. Здесь xi — координата, отсчитываемая вдоль стенки канала, ф — зависимая переменная. В работе [21] использовалось граничное условие вида = 0, = 0, где

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

Целью исследований, описываемых в предлагаемой статье, является повышение точности и скорости сходимости нового варианта метода КНК для численных расчетов двумерных стационарных ламинарных течений вязких несжимаемых жидкостей. Осуществление поставленной цели достигается реализацией следующих новых элементов вычислительного алгоритма:

1) новый вариант метода КНК — метод КНК15;

2) предлагается использовать двухпараметрический предобуславливатель в варианте метода КНК 15;

3) показано, как в отличие от конечно-разностных методов можно реализовать в различных вариантах метода КНК граничное условие в выходном сечении канала без использования рядов законтурных ячеек и без искусственного увеличения в расчетах длины канала.

Это существенное преимущество метода КНК перед конечно-разностными методами и методами конечного объема.

1. Описание метода КНК. Рассмотрим краевую задачу для системы уравнений Навье—Стокса

(V -V)V + Vp = ReAV - f, div V = 0, (xi,x2) e Q, (1)

VL = g (2)

в области Q с границей dQ. В уравнениях (1) Х1, Х2 —декартовы пространственные координаты; V = (v1(x1, х2), v2(x1, х2)) —вектор скорости; р = = р(х1,х2) —давление; f = (f\, f2) —заданная вектор-функция, Re — число

о 2 о 2 __о о

Рейнольдса, А = + щ, (V -V) = V1 ^ + V2 -¡^. Система (1) решается

при граничных условиях Дирихле (2), где g = g(x1 ,х2) = (g1,g2) — заданная вектор-функция. Давление определяется из (1), (2) с точностью до постоянной. В дальнейшем будем подбирать эту постоянную так, чтобы выполнялось условие

JJ pdx1dx2 = 0. (3)

п

В качестве области решения в настоящей работе рассматривается прямоугольник

Q = {(Х1,Х2)1 0 ^ Ж1 < L, 0 < Х2 < Н}, (4)

где L > 0 и Н > 0 — заданные длины сторон области (4) вдоль осей х1 и Х2, соответственно. Величина Н использовалась в конкретных расчетах в качестве характерной длины при обезразмеривании переменных, и она же входит естественным образом в определение числа Рейнольдса Re в (1). Далее краевую задачу для УЧП будем называть дифференциальной задачей.

В данной задаче (1)—(4) область (4) покрывается сеткой из квадратных ячеек Qij, i = 1,... ,1, j = 1,... ,J, I ^ 1, J ^ 1. Для поиска кусочно-аналитического решения задачи удобно ввести локальные координаты у1, у2

в каждой ячейке П^. Зависимость локальных координат от глобальных координат х\, х2 задается формулами ут = (хт — )/к, т = 1,2, где хт— значение координаты хт в центре ячейки П ^, а к — половина длины стороны квадратной ячейки, к = Ь/(21) = Н/(23). Пусть и(у\,у2) = (и1,и2) = = У(%1 + Х1,ъ,з, ку2 + х2,ъ,з), д(у1,у2) = р(ку! + ,ку2 + ). В локальных переменных уравнения Навье—Стокса принимают следующий вид:

/ ди*^^ дит дд \ Аит — Ре Мы—— + + ^— = Ре к2/*, т = 1,2; (5)

^ ду1 ду2 дут)

1 ( _1 + ^ = о к V дух ду2 ' ,

где А = 7^2- + . Линеаризация по Ньютону уравнений (5) приводит к формуле

е[А<+х—(РеК){и\изт+11+и{+1ит,У1+и2ит+х2+и2+1итУ2+с1)] = (6)

Здесь 8 = 0,1, 2,... — номер итерации по нелинейности, и1, и2, дя — известное приближение решения на з-ой итерации, начиная с выбранного начального приближения с индексом в = 0, Рт = Ре[к2/т — к(и1ит,,У1 + ^п*,^)], где итМ = дит/ду1, дУт = дд/дут, I = 1, 2. Здесь, как и в [13,15], введён задаваемый пользователем параметр £ с целью управления величиной числа обусловленности переопределённой системы линейных алгебраических уравнений (СЛАУ), которая должна решаться в каждой ячейке П^.

Приближенное решение в каждой ячейке П^ ищется в виде линейной комбинации базисных вектор-функций ^:

(и{,и2,д*)Т = (7)

где верхний индекс Т обозначает операцию транспонирования, а ть — задаваемое пользователем количество базисных вектор-функций. В рассматриваемом варианте метода ^ являются многочленами. Таким образом, искомое приближенное глобальное решение задачи является кусочно-полиномиальным. А часть решения (7) в каждой ячейке представляет собой локальное решение в окрестности начала локальной системы координат в ячейке.

Базисные вектор-функции имеют следующий вид:

V1 = (1, 0, 0)Т, ^2 = (У1, У2, 0)т, ^з = (У2, 0, 0)Т,

^ = (у2, —2У1У2, 0)т, ^ = (—2У1У2,У2, 0)т, <р6 = (у1, 0, 0)т,

= (0,1, 0)т, = (0, У1, 0)т, (^9 = (0,у1, 0)т, (8)

<Р10 = (0, 0,1)т, <Р11 = (0, 0,У1)т, <Р12 = (0, 0,У2)т,

<р13 = (0, 0, у2)т, Ч>14 = (0, 0, У1У2)т, <Р15 = (0, 0, у|)т

Их совокупность можно назвать соленоидальным базисом, так как ёгу ^ = 0 У1. Таким образом, достоинством предложенного варианта метода является то, что уравнение неразрывности и, следовательно, закон сохранения массы удовлетворяется численным решением задачи во всей расчётной области.

Набор базисных функций, применявшийся в [13,15], получается из набора (8), если положить ть = 12 в (7), оставив в (8) только первые 12 базисных вектор-функций.

Локальной СЛАУ будем называть далее систему уравнений, определяющую приближённое решение в каждой ячейке, а объединение всех локальных СЛАУ — глобальной СЛАУ. Основными уравнениями, определяющими решение приближённой задачи, являются уравнения коллокаций. Они получаются в результате подстановки в уравнения (6) выражений (7) и координат точек коллокаций. Количество этих точек и их расположение внутри ячейки может варьироваться в различных вариантах метода. В данной работе были реализованы три варианта задания координат точек коллокаций. Обозначим через КС число точек коллокации внутри каждой ячейки. При = 2 локальные координаты точек коллокаций таковы: (ш, ш), (-ш, ш), где ш — задаваемое пользователем значение в интервале 0 < ш < 1. При ЫС = 4 локальные координаты точек коллокаций имеют вид (±ш, ±ш). В случае ЫС = 8 координаты точек коллокаций задавались следующим образом: расположение первых четырех точек было взято таким же, как при ЫС = 4, а координаты следующих четырех точек задавались по формулам (±ш, 0), (0, ±ш). Подставляя (7), а также численные значения координат каждой точки коллокации в (6), получим 2ЫС линейных алгебраических уравнений относительно искомых Ь

<1 • ЪТ = К, V = 1,..., 2МС.

Следуя [13], дополним систему уравнений приближённой задачи в ячейке П^ линейными условиями согласования локального решения в каждой ячейке с локальными решениями, взятыми во всех соседних с ней ячейках. Эти условия записываются в отдельных точках (называемых точками согласования) на сторонах ячейки П^, которые являются общими с соседними ячейками. Условия согласования берутся в виде

л Щ51 + п(п+Т = л ^ + п(и-Т, (9)

л ^ + («+Г = л ^ + (и-У, (10)

д+ = д-. (11)

Здесь л^ = Н(п1 -щ-+П2= Пхщ--Щ2; п = (пх ,п2) — внешняя нормаль к стороне ячейки П^; ( • )га, ( • )т — нормальная и касательная составляющие вектора скорости по отношению к стороне ячейки; и+, и- —пределы функции и и д при стремлении аргументов к их значениям в точке согласования изнутри и снаружи ячейки П^. Здесь, как и в [13,15], введён параметр ^ с целью управления величиной числа обусловленности матрицы СЛАУ, которая должна решаться в каждой ячейке П^.

Для однозначного определения давления в решении задаём его значение в одной точке области либо аппроксимируем условие (3) по формуле

л (Ц д йухйу^ = л (-I * + Я д*Лук1у2^. (12)

Здесь I* —интеграл по всей области П, рассчитанный как сумма интегралов по каждой ячейке на предыдущей итерации, д* —давление в ячейке на предыдущей итерации.

Обозначим через Nm число точек согласования для составляющих вектора скорости на сторонах каждой ячейки. При Nm = 4 координаты этих точек согласования задаются формулами (±1, 0), (0, ±1). При Nm = 8 координаты точек согласования таковы: (±1, —(), (±1, (), (—(, ±1), ((, ±1), где 0 < ( < 1. В расчетах, результаты которых представлены ниже, использовалось значение ( = 1/2. Условия согласования для давления (11) задаются в четырех точках с координатами (±1, 0), (0, ±1). Используя (7), подставим координаты этих точек в каждое из трех условий согласования (9)—(11). Из первых двух условий получим 2Nm линейных алгебраических уравнений для составляющих скорости. Подстановка представления (7) в (11) дает ещё четыре линейных алгебраических уравнения согласования.

В настоящей работе давление задавалось в нижней левой вершине ячейки Qi,i или же использовалось условие (12). Если сторона ячейки лежит на границе области Q, то в соответствующих точках на стороне ячейки вместо условий согласования в локальной СЛАУ выписываются граничные условия: ит = дт, т = 1,2.

Объединяя уравнения коллокаций, согласования и уравнения, полученные из краевых условий, если ячейка Qi,j граничная, в каждой ячейке получим СЛАУ вида

^ • x°+1 = f°f\ (13)

где 1 = (Ъ1+\,... ,Ц+1тъ)т. В варианте метода КНК, описываемом в настоящей работе, система (13) переопределенная. — матрица полного ранга с системой независимых вектор-столбцов. Глобальную СЛАУ, полученную объединением всех локальных СЛАУ, решаем итерационно. При этом на каждой глобальной (s + 1)-й итерации последовательно перебираем все ячейки в области Q и решение локальных СЛАУ (13) отыскиваем, применяя QR— декомпозицию матриц с ортогональными матрицами Гивенса либо Хаусхол-дера. В правой части уравнений (9)-(11) в качестве и-, q- берем либо значения решения на (s + 1)-й итерации, если они уже сосчитаны на этой итерации, либо их значения на предыдущей итерации. Как известно, решение СЛАУ, отыскиваемое применением Q^-декомпозиции матриц, можно получить применением метода нормальных уравнений. Оба способа при отсутствии ошибок округлений дают один и тот же вектор Xjj, называемый псевдорешением переопределенной СЛАУ, матрица которой имеет независимые вектор-столбцы [23,24]. На этом решении достигается минимум евклидовой нормы невязки Afj • 1 — Известно также, что первый из указанных

способов решения СЛАУ предпочтительнее второго, так как в методе нормальных уравнений при построении решения необходимо обратить матрицу АтА, число обусловленности которой cond(ATA) = (cond(A))2. Поэтому при решении СЛАУ методом нормальных уравнений происходит более быстрое накопление ошибок округлений в процессе арифметических вычислений, чем при применении метода Q^-декомпозиции.

2. Предобуславливатели для метода КНК. Опустим далее в (13) для краткости верхние и нижние индексы:

AX = f. (14)

При численном решении уравнений Навье—Стокса различными численными методами при больших числах Рейнольдса получаются плохо обусловленные СЛАУ. Поэтому из-за ошибок округлений в процессе решения получается неприемлемо большая погрешность или расходящееся решение [24]. Чтобы уменьшить негативное влияние плохой обусловленности СЛАУ на величину погрешности решения, на скорость сходимости итераций при применении итерационных методов традиционные методы решения модифицируют, комбинируя их с дополнительными численными алгоритмами. Один из способов модификации известных методов решения СЛАУ заключается в предобуслав-ливании ее матрицы так, что решение исходной системы сводится к решению системы, матрица которой обусловлена лучше, чем матрица исходной системы. Здесь числа обусловленности матриц A и R вычислялись в спектральной норме [23,24] с использованием программных пакетов Mathematica и MATLAB. В них необходимые для этого вычисления осуществляются путем применения операторов действий с матрицами. Для исключения каких-либо ошибок проводилось сравнение результатов, полученных с помощью этих программных пакетов.

При применении описанной выше процедуры аппроксимации краевой задачи для уравнений Навье—Стокса матрица A получается с независимыми столбцами и матрица Ai = ATA — несингулярная. Спектральное число обусловленности матрицы A вычисляется по формуле

«2(A) = ^||Ai||2 ■ ||A—1 П2 = ^max/^min,

где || ■ 112 — спектральная матричная норма; amax и amin — соответственно максимальное и минимальное сингулярное число матрицы A, A—1 —матрица, обратная по отношению к Ai.

Классический диагональный предобуславливатель строится следующим образом. Сделаем в (14) замену [7]: X = CY, где C — квадратная матрица порядка ть. Тогда система (14) принимает вид ACY = f. Нужно подобрать матрицу C так, чтобы матрица B = AC была близкой к единичной матрице или, говоря в более общем смысле, чтобы число обусловленности матрицы B было меньше в сравнении с обусловленностью матрицы A. В качестве C, как правого предобуславливателя, в [7] использована матрица C = R-1, где R — диагональная матрица, на главной диагонали которой стоят квадратные корни из диагональных элементов матрицы Ai = ATA. После того как решение системы BY = f найдено, вычисляем искомый вектор X по формуле X = CY.

Введение параметров в предобуславливатель увеличивает его возможности для дальнейшего понижения числа обусловленности, так как эти параметры можно затем подбирать из требования минимизации числа обусловленности. В нашем случае мы построили предобуславливатель, зависящий от параметров £ и входящих в уравнения (6) и (9) соответственно.

Обозначим через Acoi матрицу размера 2Nc х ть, получаемую при подстановке в (6) решения (7) и координат точек коллокации; ее элементами являются коэффициенты при ть искомых коэффициентах в представлении решения (7). Заметим, что матрицу Acoi можно представить в виде Acoi = Acoi ■ D, где D = diag({, ...,£) — диагональная матрица порядка ть, а матрица Acoi

получается из СЛАУ (6) при £ = 1. То есть матрицу Б можно рассматривать как правый диагональный однопараметрический предобуславливатель матрицы Асо1.

Давление входит в уравнение количества движения (6) только в виде производных дд/ду1 и дд/ду2, поэтому 10-й столбец матрицы Асо1 нулевой. Вследствие этого матрица Асо1 неполного ранга, в матрице А^А^ при любом Ис > 1 десятый столбец также состоит только из нулей и для этой матрицы не существует обратной.

Чтобы с помощью матрицы Асо1 потенциально можно было определить решение, необходимо включить в нее строку, соответствующую уравнению (12). Эта строка при ть = 15 имеет следующий вид: Ь2^с+1 = {0, 0, 0, 0, 0, 0, 0, 0, 0, к, 0, 0, к/3, 0, к/3}. Пусть

В этой матрице десятый столбец ненулевой, и Асо\ становится при Ис ^ 7 матрицей полного ранга. При заданном численном значении полушага к элементы матрицы А зависят от £ и

Характеристические уравнения, отвечающие методам КНК12 и КНК15,— алгебраические уравнения, соответственно, 12-й и 15-й степени. Поэтому невозможно получить в замкнутом аналитическом виде выражения для собственных чисел соответствующих матриц А1. Тем не менее мы показываем ниже, что можно получить информацию о некоторых свойствах матрицы А1 и, следовательно, о свойствах рассматриваемого двухпараметрического предобуславливателя, исследуя аналитические выражения для ее элементов. Выпишем в аналитическом виде эти выражения для случая внутренней ячейки (г, ]). Так как матрица А1 симметричная, достаточно привести выражения для элементов, входящих в верхнюю треугольную часть этой матрицы. Порядок матрицы А1 равен ть, где ть = 15 для рассматриваемого здесь метода КНК15. Размер верхней треугольной части матрицы А1 равен ть • (ть + 1)/2, так что при ть = 15 в этой части матрицы А1 содержатся 120 элементов. Обозначим через ^^ элементы этой матрицы, г,] = 1,...,ть. Предполагается, что число точек коллокаций в ячейке Ыс = 8, поэтому общее число уравнений коллокаций 2Ыс = 16. Кроме того, в матрицу А включены 4 строки, соответствующие четырем условиям согласования для давления (по одному в середине каждой стороны ячейки), и 16 строк, соответствующих условиям согласования для составляющих вектора скорости (в двух точках на каждой из четырех сторон). И, наконец, одна строка учитывает интегральное условие для давления (12).

Сначала приведем выражения для тех элементов ^^^, в которые входит параметр ц и/или полушаг сетки к:

(/л, и) = (1,1), (1, 2), (1, 4), (1, 6), (2, 2), (2, 4), (2, 5), (2, 6), (2, 7), (2, 9), (3, 3),

где

(3, 5), (4, 4), (4, 6), (4, 8), (5, 5), (5, 7), (5, 9), (6, 6), (7, 7), (7, 9), (8, 8), (9, 9), (11,11), (12,12), (13,13), (15,15);

5jS2 ■ 4 ц + Sj5i ■ (1 + 4 гf ) + Sjöl ■ (12 + + öjÖ2Vx

x (10 + 8 V2) + öjSt ■ 12ц - öj¿5 ■ 12 77 + öjö6„ ■ ц - öjöl ■ 4ц - Sjö9 ■ Ц+ + öj■ (16 + Ц2) - öj¿5 ■ 2ц + öpt ■ (^ + 8ц2) + öj■ (3 + Ц2)-

öffi ■ 2ц + öjö5 ■ (^ + 8 ц2) + ■ (^ + 36) + öjöl ■ (4 + 4r?2)5j59x

2

x (12 + ц2) + öjöl ■ (16 + ц2) + öjö9 ■ (36 + +

+ 2^ + 4 öj2 öl2 + öj3 513{ f9 +2) + öj5 ¿15 ( у + 4).

Здесь öj — символ Кронекера. Остальные элементы матрицы A 1 можно записать в виде

2Nc

_ 2N _c

ßi,j = %2У2 \У2аш,гаР,3, i,j = 1,...,mb. (16)

P= l

Величины am,i, m = 1,..., 2Nc, 1 = 1,..., ть, входящие в (15), (16), являются коэффициентами уравнений коллокаций, получаемых из уравнения (6) при подстановке в него разложений (7), значения £ = 1 и локальных координат Ш, У2 точек коллокаций. Согласно (6), эти коэффициенты зависят от решения на предыдущей итерации и от числа Рейнольдса. Анализ выражений (15), (16) приводит к следующим выводам.

1. Параметр £ входит в величины ßi,j только как множитель вида £2. Отсюда следует, что поверхность К2 = К2(£, ц) симметрична относительно оси ц. Это позволяет ограничиться поиском оптимального значения параметра £ только в полуплоскости £ > 0.

2. Выражения, содержащие параметр ц, входят в элементы ßi,j как аддитивные слагаемые, при этом имеются как первая, так и вторая степени этого параметра. Поэтому ясно, что поверхность К2 = К2(£, ц) не является ни четной, ни нечетной функцией параметра ц.

3. Полушаг h входит в формулы для ßi,j только аддитивно и только во второй степени. Во многих задачах динамики жидкостей размеры расчетной области в плоскости обезразмеренных пространственных координат обычно являются величинами порядка 0(1). Кроме того, для обеспечения приемлемой точности в рассматриваемом методе КНК нужно применять, как правило, сетку, в которой не менее 10 ячеек в каждом пространственном направлении. Поэтому полушаг h < 1, и тогда h2 ^ 1. В то же время остальные выражения, входящие в ßi,j, являются величинами порядка 0(1). Отсюда следует, что число обусловленности слабо зависит от величины полушага h при h < 1. Подчеркнем, что этот вывод является общим и не зависит от специфики конкретной прикладной задачи. Это обстоятельство можно эффективно использовать при поиске оптимальных значений параметров , , обеспечивающих минимум

числа обусловленности матрицы A1: достаточно найти эти оптимальные значения с использованием численного решения, найденного методом КНК на сравнительно грубой сетке. Затем оптимальные ц можно будет использовать в расчетах на сетках, имеющих намного меньшие шаги.

Так как элементы матрицы Acoi зависят от решения на предыдущей итерации, дальнейшее исследование свойств обусловленности необходимо осуществлять на заданной сетке при решении конкретной задачи. Соответствующие вычислительные эксперименты описываются ниже в п. 5.

В [13] для нахождения оптимальных значений £opt, 'Qopt в любой ячейке пространственной сетки из требования минимизации числа обусловленности K2(C,v) применялся метод равномерного поиска с переменным шагом. Поиск минимума функции K2(^opt, ^opt) был осуществлен в [13] для случая метода КНК12. Оказалось, что в точке минимума число K>2(Copt, Wopt) удовлетворяет неравенствам 3 < K2(£opt, rjopt) < 10. Кроме того, было установлено, что оптимальные значения £opt, ^opt слабо зависят от положения конкретной ячейки в расчетной сетке, по меньшей мере, в случаях тех тестовых и эталонных задач, которые рассматривались в [13].

3. Вариант алгоритма Крылова с редукцией базиса подпространства. Для ускорения сходимости итераций во всех новых вариантах метода КНК, описываемых в настоящей работе, здесь применялся новый вариант известного метода [17], основанный на подпространствах Крылова и подробно изложенный ранее в [15, 22]. Важным отличием от [18] варианта метода Крылова, предложенного в [22], является автоматизация редукции базиса подпространства Крылова в области малых невязок, что позволяет избежать возможные АВОСТы программы в области малых невязок исходного СЛАУ при стремлении получить ее решение с достаточной точностью.

4. Ускорение сходимости итераций с помощью многосеточного алгоритма. Основная идея многосеточных алгоритмов состоит в селективном демпфировании гармоник погрешности решения задачи [16,25]. В методе КНК, как и в других методах, количество итераций, необходимых для достижения заданной точности приближения к предельному решению, зависит от начального приближения. В работе применялась операция продолжения вдоль восходящей ветви V-цикла — расчеты на последовательности измельчающихся сеток — в качестве способа получения хорошего начального приближения для итераций на самой мелкой сетке среди сеток, используемых в многосеточном комплексе. Переход от грубой сетки к более мелкой делается с помощью операторов продолжения. Проиллюстрируем алгоритм операции продолжения на примере составляющей скорости и1(у1, у2,Ь1,..., Ь15). Пусть h1 = h, где h — полушаг грубой сетки, и пусть h2 = h1/2 — полушаг мелкой сетки, на которой нужно найти разложение функции U1 по базису.

Ш!аг 1. Пусть Х1, Х2 — глобальные координаты центра ячейки грубой сетки. Сделаем в полиномиальное выражение для U1 следующие подстановки: yi = (xi — Xi)/h1, l = 1, 2. В результате получаем многочлен

U1(X1,X2 ,Ъ1,...,Ъ1ь) = X1 — X1, Х2 — Х2 ,b1,...,b15^. (17)

Ш!аг 2. Пусть (.Х-!,^) —глобальные координаты центра любой из четырех ячеек мелкой сетки, содержащихся в ячейке грубой сетки. Сделаем в (17) замену Х[ = Х1 + у ■ 1ъ2, I = 1, 2. В результате получим многочлен второй степени 11\ = Р (У1,У2,Ь 1,..., Ь 15) от переменных уц, у2 с коэффициентами Ь1,... ,Ъ 15.

Приведем выражения для коэффициентов у (] = 1,..., 15) представления решения в ячейке мелкой сетки с полушагом Ь,2 в терминах коэффициентов Ь1,..., Ь15 представления решения в ячейке с полушагом Ь,1 = 2^:

Ъ\ = Ь1 - 5X1 (Ь2 - М^) - 5Х2(ЬЗ +265^^1 -Ье8x2), Ь2 = О^СА + М^), Ъз = а1[Ь3 + 2(Ъ5>6Х1 - Ь65х2)}, ЪА =а2Ь4, Ь5 = а2Ь5, Ь6 = <Г2Ь6, У = Ь7 - 5х1 (Ь8 - Ьд5х1) + 6х2Т1, Ь8 = а1(Ь8 - 2Ьд5х1 + 2Ь45х2),

Ьд = (Г2Ъд, Ъю = Ъю - бх^ - бх^Ьы - ^58X2)^ Ь 11 = (71(Т2 - У36X1), У12 = °1(Ъ12 - Ьи§Х1 - 2Ь15§Х2), Ь13 = (Г2Ъ13, Ь14 = 02Ьи, Ь15 = а2^5,

где 5X1 = -(XI - Х{)/Ь,1, 01 = Л.2/Л4, 02 = Т1 = Ь2 - 2Ь45х1 + Т2 = Ьц - Ь13§Х1 - Ьи5х2. Заметим, что приведенные выше выражения для Ь1,... ,Ьд совпадают с выражениями, представленными в [13,15] для случая, когда ть = 12 в (7).

5. Результаты численных экспериментов. Тестирование. Рассмотрим следующее точное решение уравнений Навье—Стокса (1) [26]:

-2(1+ Ж2) 2(1+ Ж1) и1 = Т^--чо . /-, .-, и2 =

(1+ Х1)2 + (1+ Ж2)2' (1+ Х1)2 + (1+ Х2)2'

2

Р = - (1+ ,1)2 + (1+ ,2)2 , 0 ^ ^ ^ 1.

(18)

Заметим, что функции и1(х1 ,х2) и и2(х1,х2) описывают поле скорости с нулевой дивергенцией. Далее,

У У р (1x1(1x2 = 4С - ж 1п 2 - 2г Ы2^-2) - ^^2)

2

и -0.462 613 146 772 815 498 72,

где г = \/-1, С — постоянная Каталана, С и 0.915 965 594177219 015 05, ^2(2) — полилогарифмическая функция. Чтобы обеспечить выполнение равенства (3) с погрешностью, не превышающей погрешность машинных вычислений, давление р в (3) заменялось на величину р = р + 0.462 613 146 772 815 5.

Для вычисления среднеквадратичных величин погрешности решения использовались следующие формулы:

Бгг(и(Л,)) =

Егг(р(Ь)) =

1 13 2

277 ^ ^ - )

г=1]=1 и=1

1 2

I 3

1 Е - р% )2

и

г=1 3 = 1

(19) 553

где uex и — вектор скорости и давление, вычисленные из точного решения (18). Величины u¿j и pij обозначают численное решение, вычисленное в центре ячейки Qij по описанному выше методу КНК. Величины vu и vp — порядки сходимости погрешностей uij и pij соответственно, которые вычислялись по известным формулам [12,27]. Пусть Ъ^г, s = 0,1,... — значения коэффициентов bi j i в (7)) на s-ой итерации. Использовалось следующее условие окончания итераций:

Sbn+1 < е, (20)

где öbn+1 = maxi max 16™+/ — , а e < h2 — малая положительная ве-

i, j г ' г

личина. В дальнейшем будем называть величину öbn+1 псевдопогрешностью приближенного решения. Наряду с условием (20) для окончания итераций по нелинейности также применялся следующий критерий:

¿ura+1 = || ura+1 — ura ||<е2, (21)

где || ■ || —евклидова норма вектора, £2 —заданная малая положительная величина.

Была проведена серия расчетов с целью изучения влияния конкретного вида предобуславливателя на сходимость итераций по методу КНК15. В этой серии расчетов критерием останова расчета было выполнение неравенства öbs < 10-9. Использовалась восходящая ветвь многосеточного алгоритма с переходами от грубой сетки из 5 х 5 ячеек к сеткам из 10 х 10, 20 х 20, 40 х 40 и 80 х 80 ячеек. Результаты представлены на рис. 1. Крестиком на оси п показан тот номер итерации, начиная с которого расчет в рамках многосеточного алгоритма осуществляется на сетке из 80 х 80 ячеек. Видно, что отсутствует рост погрешности Err(ura) на всех рассмотренных сетках с увеличением числа итераций также при отсутствии предобуславливателя, то есть когда £ = ц = 1. Наоборот, в случае метода КНК 12 указанная погрешность начинает расти после перехода к расчету на сетке из 80 х 80 ячеек, см. [15, Fig. 2]. Диагональный предобуславливатель имеет удобство в использовании, так как он не содержит регулируемых параметров и обеспечивает устойчивость счета по методу КНК15.

Изучалось также влияние включения условия (12) в локальную СЛАУ на скорость сходимости метода КНК 15 на сетке из 40 х 40 ячеек при числе Рейнольдса Re = 1000, L = Н = 1 в (4), Nc = 8. Критерием останова счета было выполнение неравенства öbn < 10-12. Для сравнения был проведен расчет, когда условие (12) не использовалось, а вместо него задавалось значение

Рис. 1. Погрешность Err(un), полученная при расчете по методу КНК15 при использовании различных предобуславли-

вателей: (-) £ = £opt = 0.05, -q = ^opt =

= 1.75; (---) £ = v = 1; (......) —диагональный предобуславливатель [Figure 1. Error Err(un) obtained at the computation by the method CLS1b at

the use of different preconditioners: (-)

£ = Copt, v = ^opt; (---) £ = v = i;

(......) — the diagonal preconditioner]

давления в точке с координатами х\ =0, = 0. В процессе итераций по методу КНК15 абсолютная величина интеграла (3) падала со значения порядка 10-3 до величины порядка 10-13-10-14, то есть до величины порядка машинных ошибок округления, накопленных в процессе решения задачи при расчетах по Фортран-программе с двойной точностью. Это может служить одним из критериев правильности программной реализации представленного выше метода КНК15. Число невязок в алгоритме подпространств Крылова равно 10. Оказалось, что использование условия (12) вместо задания давления в одной точке позволяет уменьшить количество итераций по методу КНК15 в 12 раз.

На рис. 2, 3 приводятся результаты численных экспериментов, в которых при решении уравнений Навье—Стокса по методу КНК тестовой задачи с точным решением (18) использовались только два из описанных в предыдущих разделах трех способов ускорения сходимости: двухпараметрический предобуславливатель и метод подпространств Крылова. Число Рейнольдса Ре = 1000, Ь = Н = 1 в (4), Ыс = 8. Расчеты, результаты которых представлены на рис. 2 и 3, были выполнены на сетке из 40 х 40 ячеек. Критерием останова счета было выполнение неравенства 5Ьп < 10-12.

Расчету без применения алгоритма Крылова соответствует на рис. 2 случай к = 0. В этой серии расчетов в переопределенной СЛАУ (13) использовалось условие (12). На рис. 2 видно, что с увеличением числа невязок к, используемых в методе Крылова, скорость сходимости численного решения по методу КНК 15 растет. Количество итераций Иц, необходимых для обеспечения выполнения неравенства 5Ьп < 10-12, составляло 56392, 7753, 5784 и 4936 соответственно при к = 0, 2, 10 и 15. Таким образом, применение алгоритма Крылова при Ре = 1000 с к = 15 позволило уменьшить количество итераций, требуемых для сходимости приближенного решения, в 11.4 раза по сравнению со случаем к = 0. При дальнейшем увеличении числа невязок

Рис. 2. Расчеты по методу КНК1в. Влияние числа невязок к, используемых в алгоритме Крылова, на скорость сходимости величины log10 Err(un), где п — число итераций

[Figure 2. Computations by the CLS1g method. The influence of the number of residuals к employed in the Krylov's algorithm on the convergence rate of the quantity log10 Err(un), where n is the number of iterations]

Рис. 3. Сравнение профилей приближенного и точного решений при х2 = 1/2.

Сетка 40 х 40 ячеек [Figure 3. Comparison of the profiles of the approximate and exact solutions at x2 = 1/2. The 40 x 40 grid]

-- fc = 0

_ fc = o, £ = ?7 = 1

■ • к = 2

- к = 10

о к = 15

10000 20000

30000 п

40000 50000

0.5

_ -е__о— в —о—& —о—в---е—о— -в——о- — 0

к, используемых в алгоритме Крылова, увеличивается пропорционально к3 машинное время, требуемое для решения ортогональным методом переопределенной системы для поправочных коэффициентов при невязках. Это приводит к увеличению суммарного машинного времени, требуемого для численного решения задачи. Поэтому для каждого конкретного случая применения алгоритма Крылова для ускорения итерационного процесса решения задачи методом КНК существует свое оптимальное значение количества невязок к в алгоритме Крылова.

Также был проведен расчет без использования предобуславливателя, то есть когда £ = ^ = 1; кроме того, алгоритм на основе подпространств Крылова не применялся, см. штрихпунктирную линию на рис. 2. Псевдопогрешность öbn в процессе итераций уменьшалась до заданной величины е = 10-12, но, несмотря на это, не было сходимости численного решения для скорости. Этот же эффект наблюдался и в случае метода КНК12 [15].

На рис. 3 дано сравнение профилей компонент точного решения и решения, полученного методом КНК15. Компоненты vi, V2 и р приближенного решения изображены символами △, о и v, те же компоненты точного решения — сплошными, штриховыми и штрих-пунктирными линиями соответственно. Здесь видно хорошее согласие между численными результатами и аналитическим решением.

В [15] были представлены результаты серии тестовых расчетов по методу КНК12 с целью определения порядков сходимости ии, vp при Re = 1000. Аналогичные расчеты были проведены также в случае метода КНК15. Оказалось, что порядок сходимости vu численного решения по этому методу намного выше, чем в случае метода КНК12, и превышает значение 2.15. Кроме того, увеличение количества базисных вектор-функций с 12 до 15 уменьшает погрешность скорости на три десятичных порядка и давления — на два десятичных порядка по сравнению с методом КНК12.

Были проведены расчеты с применением восходящей ветви V-цикла с целью выяснения, как влияет применение только многосеточного алгоритма на ускорение сходимости метода КНК15. Были также проведены расчеты, в которых движение по восходящей ветви многосеточного V-цикла применялось совместно с алгоритмом ускорения, основанным на подпространствах Крылова. В этой серии расчетов, выполненных при Re = 1000, применялись либо все сетки из последовательности 5 ■ 2m х 5 ■ 2т (т = 0,..., 4), либо только сетка из 80 х 80 ячеек. Пусть Kmgr — количество последовательно используемых сеток в многосеточном комплексе. Если Kmgr = 1, то это означает, что в расчете используется только одна сетка, и это самая мелкая сетка с числом ячеек 80 х 80. Расчету без применения ускорения сходимости с помощью подпространств Крылова соответствует случай к = 0, где к — число невязок, используемых в поправке по Крылову. Фактор ускорения итерационного процесса AF в результате применения того или иного способа его ускорения вычисляется как отношение времени счета при Kmgr = 1, к = 0, ко времени счета при применении последовательности из нескольких сеток (Kmgr > 1, к = 0) или же последовательности сеток в сочетании с применением алгоритма Крылова на каждой сетке (Kmgr > 1, к > 1). Во всех расчетах этой серии использовались оптимальные значения £opt = 0.02, ,qopt = 1.08 в двухпа-раметрическом предобуславливателе; кроме того, в локальные матрицы Aij

в (13) включалась аппроксимация (12) интегрального условия для давления. Критерием сходимости было выполнение неравенства 5Ьп < 0.5 ■ 10-9. Были рассмотрены следующие комбинации величин (Kmgr,k): (1,0), (5,0), (1,5), (1,15), (5,5), (5,10), (5,15), (5,20). Оказалось, что величина AF достигает наибольшего значения AF = 362.38 при Kmgr — 5, к — 10.

Общий вывод из данной серии расчетов следующий: наиболее существенный вклад в ускорение итерационного процесса решения уравнений Навье— Стокса по методу КНК12 дает совместное применение всех трех способов ускорения сходимости итераций.

Течение Пуазейля. Эта задача является примером течения в канале с открытой границей свободного истечения жидкости. Как отмечалось во введении, неправильная формулировка граничных условий на такой границе может привести к развитию неустойчивости расчета или к большим ошибкам решения внутри расчетной области. Те формы граничных условий, которые использовались различными авторами, носят искусственный характер, и для уменьшения ошибок, вызываемых применением таких условий, приходится брать очень большую длину канала, в котором рассматривается течение [21].

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

гН гН

/ ^1(ж1,ж2)^ж2 — v1(0,x2)dx2 — const, 0 ^ x1 ^ L. (22) Jo Jo

Расчет течения Пуазейля будем осуществлять в прямоугольной области (4). Следуя [28], введем безразмерные переменные x1,x2,v1,v2,p в (1), (2) по следующим формулам:

_ Х1 _ Х2 _ V1 _ V2 _ р ,поЛ

X1 = н, Х2 = н, V1 = ZV V2 = ZV p = -рй%. (23)

Здесь p — заданная плотность жидкости, р = const > 0, Ucp —средняя величина составляющей скорости V1 (0,Х2) на входе в канал. Зададим функцию v1(0,x2) в виде параболы: г^(0, ж2) = 4UH-2х(Н — х), где U — максимальное

значение составляющей скорости ^(0, ,2). Тогда размерный объемный расход

¡'Н

жидкости Q в сечении х1 = 0 выражается формулой Q = / v1(0,x2)dx2 =

o

= 2UH/3. Отсюда следует, что Ucp = 2U/3. Поэтому максимальная безразмерная скорость U во входном сечении Х1 = 0 равна U = U/(3U) = 1.5. Пусть Q — безразмерный объемный расход в сечении х1 = 0. Тогда Q = = Q/(UcpH) = 1. Число Рейнольдса вводится по формуле [28] Re = UcpH/u, где v — кинематическая вязкость жидкости, v = р/р. С учетом (23) легко найти, что ^1(0, Ж2) = 6,2(1 — %2). Ниже черточки в обозначениях безраз-

мерных величин опущены для краткости. Осуществляя обезразмеривание в формуле (22), легко получить равенство

/ v1(x1,x2)dx2 = / v1(0,x2)dx2 = Q = 1, 0 ^ x\ ^ L/H. Jo J 0

(24)

Начальное приближение для составляющей скорости ^1(^1,^2) задавалось при х1 > 0 из требования выполнения равенства (24): v1 = 4(2 — |ж2 — 11). Из этой формулы следует, что max г>1 = 2. Заметим, что начальное приближение v1(x1,x2) = 6ж2(1 — ж2) при ж1 = 0 и v1(x1,x2) = 2 — |4ж2 — 2| при 0 < Ж1 ^ L/H является разрывным в сечении Ж1 = 0. Это, однако, не вызывало каких-либо проблем при расчетах течения Пуазейля по методу КНК.

При численных расчетах рассматриваемой задачи уравнение (24) включалось в переопределенную систему (13) наряду с условием для давления (12). При этом оно включалось на грани у1 = 1 в каждой ячейке и имело следующий вид:

2h[bi ,jt 1 + bi ,jt 2 + bi 4 + (1/3)6, 6] = fi ,j,

где bi,j,m, m = 1, 2, 4, 6 — искомые коэффициенты разложения решения по

базису (см. (7)), fi,j = 1.0 — Q* + Q*,, Q* = v^x^ij + h, Х2)dx2 — интеграл

0

по всему сечению 0 ^ Х2 ^ 1, вычисленный как сумма интегралов по каждой ячейке с использованием решения, полученного на предыдущей итерации. Здесь

Q* j = V1( X1, i ,j + h, X2) dx2 = 2h[ai ,jt 1 + ai ,jt 2 + ai ,jt 4 + ^i^6 ),

Jx2,i,j-h v 37

ai,j,m, ш = 1, 2, 4, 6 — известные коэффициенты разложения решения по базису, полученные на предыдущей итерации.

Точное решение рассматриваемой тестовой задачи имеет следующий вид:

Ыж1, ж2) = 6Ж2(1 — Х2), V2(X1,X2) = 0,

dv 12 (25)

= — К", 0 < Ж1 < L/H, 0 < Ж2 < 1. ox1 Re

Величина др/дх1 постоянна внутри канала в соответствии с точным решением. Для давления задавалось нулевое начальное приближение во всей расчетной области, за исключением точки с координатами Ж1 = Х2 = 0: в этой точке задавалось значение р = 1. Ненулевой постоянный градиент давления формировался в расчетной области в процессе итераций по методу КНК.

В вертикальном слое ячеек, примыкающем к выходному сечению, задавались условия согласования (9)—(11) в точках стыка правых сторон соседних ячеек, расположенных вдоль линии Ж1 = L/H. Эти условия отражают лишь факт непрерывности решения уравнений Навье-Стокса и могут применяться, наряду с условием постоянства расхода, при численном решении любых стационарных и нестационарных задач, в которых есть открытая граница истечения жидкости.

В табл. 1 приводятся результаты расчетов течения Пуазейля по методу КНК, в которых применялись два из описанных выше трех способов ускорения сходимости: двухпараметрический предобуславливатель и метод подпространств Крылова. Число Рейнольдса Re = 500, L/H = 3, число точек коллокации Nc = 8. Во всех расчетах, результаты которых приведены в табл. 1, в предобуславливателе использовались значения параметров £ = 0.04, ^ = 1.0. Критерием останова счета было выполнение неравенства (20) с е = 10-15. Число невязок в алгоритме подпространств Крылова равно 12.

Вместо погрешности для давления (19) в этой серии расчетов использовалась погрешность вычисления производной dpn/dxi:

Err(

/ \ 1 ^ ^ / \dx\J IJ ЕЕ

1 ^ ^ ( bi,j,u + 12 \ 2

h ReJ

1 1

Аппроксимация интеграла от давления, входящего в (3), подсчитывалась по формуле

I J

МЛ = EEJJ Qn(yi,y2)dyidy2.

i=1 j =1 Qij

Если итерационный процесс по предлагаемому методу КНК сходится, то должно иметь место следующее свойство: Int(gn) ^ 0 при п ^ те. Подсчитывалась также абсолютная погрешность вычисления объемного расхода 5Qn в выходном сечении: 5Qn = IQn — 1|. Наконец, N-lt — количество итераций по методу КНК, требуемое для выполнения неравенства (20). В третьем столбце в скобках приводятся значения Nu, полученные без использования в расчете условий согласования (9)—(11) в выходном сечении.

Оказалось, что применение этих условий в выходном сечении канала позволяет существенно снизить требуемое количество итераций в тех случаях, когда для достижения сходимости требуется задавать в (20) весьма малое значение величины е. Например, при е = 10-15 для сходимости на пространственной сетке из 30 х 10 ячеек понадобилось 7774 итераций при включении условий (9)—(11) в переопределенную систему, а без их применения потребовалось 13359 итераций. Таким образом, счет ускорился в 1.65 раза. При этом погрешность Err(u(h)) = 4.466 ■ 10-13, абсолютная погрешность в безразмерном расходе жидкости, который равен 1 согласно (24), была меньше, чем 10-13. Как видно из табл. 1, эффект от применения условий согласования в выходном сечении канала усиливается при увеличении числа ячеек пространствен-

Таблица 1

Погрешности Err(u"), Err(, 5Qn, Int(g") на последовательности сеток [The errors Err(un), Err( , SQ", Int(qn) on a sequence of grids]

I J Nit Err(u") Err( dpn dx-\ SQn Int(gn)

15 5 1231(1866) 1.234 • 10-13 4.884 10- -14 1.177 • 10- 14 1.911 • 10- 16

30 10 7774(13359) 4.466 • 10-13 1.561 10- -13 7.105 • 10- -1b 1.002 • 10- 17

60 20 36038(141607) 1.013 • 10-12 1.016 • 10- 12 1.665 • 10- 14 2.196 • 10- 17

ной сетки. Например, за счет этого при счете на сетке из 60 х 20 ячеек удается уменьшить требуемое количество итераций в 141607/36038 ~ 3.9 раза.

Проводились также расчеты течения Пуазейля по методу КНК без применения алгоритма подпространств Крылова; на сетке из 30 х 10 ячеек это приводило к увеличению требуемого количества итераций в 34680/7774 ~ 4.5 раза.

Необходимо отметить, что сетка из 30 х 10 квадратных ячеек достаточно грубая, шаг сетки равен 0.1. Поэтому при использовании численного метода второго порядка точности разумно ожидать на такой сетке абсолютную погрешность порядка 0(10-2). Но реальный расчет по методу КНК (см. табл. 1) обеспечивал в норме пространства ¿2 погрешность 0(10-13), которая на 11 десятичных порядков меньше ожидаемой погрешности. Этот эффект можно объяснить тем, что полиномиальная функция

представляющая приближенное решение для т в методе КНК, является собственной функцией решаемой задачи, и она способна обеспечить с компьютерной точностью совпадение численного решения по методу КНК с точным аналитическим решением (25).

Точное решение для VI зависит только от координаты Х2. Поэтому в локальных координатах у1, у2 точное решение должно зависеть только от координаты у2. Это можно использовать для проверки правильности работы компьютерной программы, реализующей метод КНК. Действительно, если в программе нет ошибок, то значения коэффициентов 62, &4, в (26) должны быть либо равны нулю, либо иметь очень малые абсолютные значения. Для проверки рассмотрим ячейку с индексами г = 30 и ] = 5 (расчет течения Пуазейля проводился на сетке из 30 х 10 ячеек). В локальных координатах мы получили, что У,1 = 1.485 + 0.03^2 - 0.015у|. Подставляя в это выражение формулу для у2: у2 = (х2 — х2с)/к = (х2 — 0.45)/0.05, получим VI = 6.0000000000000000^2 — 6.0000000000000000^2, то есть совпадение с точным решением (25).

На рис. 4 представлены результаты расчета течения Пуазейля по методу КНК15 на сетке из 30 х 10 ячеек. Видно, что при всех х1 € [0,Ь/Н] сохраня-

Рис. 4. Поверхности решения, полученные по методу КНК15: слева — поверхность vi = и1(ж1,ж2); справа — поверхность р = р(ж1,ж2)

[Figure 4. Solution surfaces obtained by the CLS16 method: on the left is the surface v1 = и1(ж1,ж2); on the right is the surface p = р(ж1,ж2)|

U1 = 61 + &2У1 + Ь3у2 + Ь4У2 — 2Ь5У1У2 + Ъбу1,

(26)

1.5 «il.O 0.Е

3 0

ется одинаковая параболическая форма профиля составляющей скорости v\. Далее из рис. 4 (справа) видно, что профили давления имеют одинаковый наклон во всех сечениях = const £ [0,1], что согласуется с точным решением.

Обтекание ступеньки. Эта задача при Re = 800 часто применяется для тестирования численных алгоритмов решения уравнений Навье—Стокса. Рассматривается стационарное ламинарное течение вязкой несжимаемой жидкости в канале длиной L и высотой Н. На входе в канал расположена прямоугольная ступенька с высотой Н/2 и длиной Lst (см. рис. 5). Поэтому область, в которой рассчитывается течение жидкости, в отличие от течения Пуазей-ля является не прямоугольником (4), а объединением двух прямоугольников: Q = Q U Q2, где

Hi = {(Ж1, х2) | 0 < XX < Let, Н/2 < Ж2 < н},

(27)

Q = {(х\, Х2) | Let ^ Х\ ^ L, 0 ^ Х2 ^ н}.

Краевые условия для составляющих вектора скорости на входе в канал имеют следующий размерный вид:

16 / Н^

--—

Щ (0, Х2) = ^ и(х2 - ~2)(Н - х2), М0, Х2) = 0,

где и > 0 — заданное максимальное значение составляющей скорости щ в сечении х\ = 0. Размерный объемный расход жидкости в поперечном сечении Х\ = 0 дается формулой

(н 2 Н

Qo = щ {0,х2)(1х2 = - и • —. (28)

/Н/2

3

2

Отсюда следует, что средняя величина составляющей скорости щ (0, Ж2) в интервале [Н/2, Н] есть Ucp = 2U/3. В результате обезразмеривания по формулам (23) краевое условие для безразмерной составляющей скорости щ (0,ж2) и безразмерный расход Q0 принимают следующий вид:

щ(0, Ж2) = 24(^2 - 0.5)(1 - Ж2), Qo = 0.5.

Аналогично случаю течения Пуазейля легко показать, что в каждом сечении х\ = const, 0 ^ Х\ ^ L/H, должно иметь место равенство

/ V\(x\,x2)dx2 = Qo = 0.5,

(29)

где х2ь = 1/2 при 0 ^ X\ ^ Lst/Н и х2ь = 0 при Lst/H < x\ ^ L/H.

X2 x3 >

X2 ,

....... ....... Xi

Lst/H Xi

L/H

Рис. 5. Вид расчетной области в задаче об обтекании ступеньки [Figure 5. Form of the computational region in the backward-facing step problem]

Начальное приближение для составляющей скорости vi(xi, х2) задавалось аналогично случаю течения Пуазейля с применением треугольных профилей:

( 24(^2 - 0.5)(1 - х2), xi = 0,1/2 < х2 < 1, Vi(xi ,Х2) = { 2(1 - 4 |Ж2 - 0.751), 0 < xi < Lst/H, 1/2 < ,2 < 1, (30) { 2 (0.5 - |Ж2 - 0.5|), xi > Lst/H, 0 < ,2 < 1.

Легко проверить, что функция (30) удовлетворяет соотношению (29).

Невязки, входящие в поправку по Крылову, сохранялись в числовых массивах отдельно для каждой из двух подобластей (27). При этом в алгоритме Крылова использовалось 8 невязок. Численное решение переопределенной системы (14) находилось в каждой ячейке с помощью метода отражений Ха-усхолдера [29]. Этот метод обеспечивал ускорение счета по сравнению с методом вращений Гивенса [29] примерно в 1.5 раза. В двухпараметрической матрице A использовались следующие значения параметров £, £ = 0.04, ^ = 1. В качестве примера были выполнены расчеты обтекания ступеньки при Re = 800, L = 14, Lst = 2.5. Расчеты проводились на равномерной сетке из 560 х 40 квадратных ячеек. В качестве критерия окончания итераций по нелинейности использовалось выполнение неравенства (21) с £2 = 5 ■ 10-8.

В случае применения алгоритма Крылова для сходимости потребовалось 109500 итераций, а без его применения — 857000 итераций. Таким образом, применение алгоритма Крылова позволило уменьшить потребное для сходимости число итераций (следовательно, и машинное время) в 7.83 раза.

При отсутствии алгоритма Крылова в методе КНК были получены следующие значения величин Xi, Х2, A3: Xi = 6.10, Х2 = 4.91, Х3 = 10.32.

При вышеуказанном значении числа Рейнольдса образуются два вихря — непосредственно за ступенькой и на верхней стенке канала (см. рис. 6). Координата xi = Xi точки повторного присоединения потока на нижней стенке, координата Xi = Х2 точки отрыва потока от верхней стенки и координата Xi = Х3 точки повторного присоединения потока к верхней стенке вычислялись по алгоритму, описанному в [33]: эти точки находились как точки, в которых вязкое касательное напряжение обращается в ноль на стенке. То есть в этих точках выполняется равенство dvi/dx2 = 0. Из этого равенства и из определения локальной координаты у2 следует равенство dui/dy2 = 0. Так как в каждой ячейке Q^ функция Ui(yi,y2) имеет аналитическое представление (см. табл. 1), производная dui/dy2 вычислялась в ячейке Q^ по формуле dui/ду2 = bi,j,3 -2(bi,j,5yi -bi,j,6y2). Уравнение dui/ду2 = 0 решалось методом половинного деления с заданной погрешностью ö = 10-i5. Для этого сначала задавался интервал для поиска корня уравнения.

В табл. 2 результаты вычисления величин Xi, Х2, Х3 по методу КНК сравниваются с ранее опубликованными данными других авторов. При этом

1.0 0.8 0.6 0.4 0.2 0.0

Рис. 6. Изолинии функции тока при обтекании ступеньки (Re = 800) [Figure 6. Streamlines in the backward-facing step problem (Re = 800)]

мы ограничились в основном работами, опубликованными в 2008 году и в последующие годы. Данные из работы Гартлинга [28] приведены потому, что именно с результатами этой работы авторы более поздних работ часто сравнивали свои результаты. Видно, что имеет место совпадение значения Х\ с данными Гартлинга [28], которые на сегодняшний день считаются эталонными. Что касается данных для величин и Хз, то имеются некоторые расхождения с данными других авторов. Однако данные по Х2 и Хз, полученные методом КНК, находятся внутри интервалов разброса этих данных в работах других авторов. При этом величина Х2 отличается от значения, приведенного в [28], на 4%.

На рис. 7 графики кривых щ(х\,х2) и г>2(ж1,ж2) в сечении х\ = Lst/H + 7 сравниваются с результатами [33]. Данное сечение пересекает верхнюю вихревую зону почти в ее середине, см. рис. 6. Сплошные линии — результат работы [33], темные кружки — расчет по методу КНКi 5 настоящей работы. Видно хорошее соответствие с результатами [33].

Необходимо отметить, что в приведенных расчетах по методу КНК 5 использовалась довольно грубая сетка из 560 х 40 ячеек. Для сравнения: в рабо-

Таблица 2

Сравнение результатов моделирования при обтекании ступеньки; Re = 800 [Comparison of modelling results obtained in the problem of the backward-facing step flow]

Scientific researches X3 — X2 X2 X3

by Gartling [28] 6.10 5.63 4.85 10.48

by Erturk [21] 5.92 5.54 4.74 10.28

by Martynenko [30] 6.10 5.63 4.84 10.47

by Rouizi et al. [31] 5.88 5.57 4.71 10.28

by Parsani et al. [32] 5.84 5.73 4.61 10.34

by Roberts et al. [33] 6.10 5.63 4.85 10.48

by Bustamante et al. [34] 5.99 5.17 4.82 10.00

Current research 6.10 5.40 4.91 10.31

l/l 1)2

Рис. 7. Профили составляющих скорости vi, v2 в сечении xi = Lst/H + 7 при Re = 800 [Figure 7. The profiles of velocity components vi and v2 in the cross section xi = Lst/H + 7 at Re = 800]

те [30] использовалась сетка из 1400 х 100 ячеек. Важно также отметить, что в полученном численном решении закон сохранения массы (28) выполнялся с абсолютной погрешностью, меньшей, чем 5 ■ 10-6.

Течение вязкой жидкости в каверне с движущейся крышкой.

На высокоточных решениях этой эталонной задачи проверяются возможности известных или вновь создаваемых численных методов.

В рассматриваемой задаче расчетная область — каверна — квадрат (4) со стороной Ь = 1, начало координат находится в ее левом нижнем углу. Верхняя крышка каверны движется в безразмерных величинах с единичной скоростью в положительном направлении оси Ох\. Остальные стороны каверны (4) покоятся. На всех сторонах заданы условия прилипания: ы = 1, = 0 при х2 = 1 и ут = 0, т = 1, 2 на остальных сторонах.

Течение в каверне с движущейся крышкой имеет сингулярности в верхних углах области. Их влияние на точность численного решения усиливается с увеличением числа Рейнольдса. Поэтому при больших числах Рейнольд-са для получения более точного решения необходимо применять адаптивные сетки с более мелкими ячейками в окрестности сингулярностей [10,35]. Здесь использованы только равномерные сетки, размер которых не превосходил 256 х 256 ячеек.

Были проведены расчеты рассматриваемой задачи по методу КНК15 при числах Рейнольдса Ре = 100, 500, 1000, 1500, 2000, 2500. Перед каждым из этих расчетов находились оптимальные значения цорь параметров £, ц предобуславливателя, описанного в п. 2, методом равномерного поиска с переменным шагом. Оказалось, что при всех рассмотренных числах Рейнольдса поверхность к2 = симметрична относительно что подтверждает

результаты аналитического исследования, приведенного в п. 2.

С целью уменьшения расхода машинного времени при переходе к расчету течения при более высоких значениях числа Рейнольдса численное решение по методу КНК15, полученное при меньшем числе Рейнольдса, использовалось в качестве начального приближения для решения с более высоким числом Рейнольдса. Это позволяло уменьшить потребное машинное время не менее чем в два раза. Указанный прием применялся ранее в [38].

По аналогии с тестовыми расчетами по методу КНК15, где численное решение сравнивалось с точным решением (18), были проведены расчеты рассматриваемой задачи при Ре = 100 с целью получения данных об эффективности совместного применения двухпараметрического предобуславливателя и алгоритмов Крылова и Федоренко. В этой серии расчетов применялись либо все сетки из последовательности 5 ■ 2т х 5 ■ 2т (т = 0,..., 3), либо только сетка из 40 х 40 ячеек. В критерии сходимости по псевдопогрешности (20) использовалось значение е = 10-9. Во всех этих расчетах полагали = 0.17, Цорь = 1.75 в двухпараметрическом предобуславливателе. Кроме того, в локальные матрицы А^ в (13) включалась аппроксимация (12) интегрального условия для давления (3). Были рассмотрены следующие комбинации величин (Ктё1,к): (1,0), (4,0), (1,9), (4,8), (4,9), (4,10). Оказалось, что величина ЛР достигает наибольшего значения АР = 26.86 при Ктех = 4, к = 9.

В работе [15] описана аналогичная серия расчетов по методу КНК12 с теми же шестью комбинациями величин (КтёХ,к). При этом не применялась аппроксимация (12) интегрального условия для давления, а задавалось дав-

ление р = 1 в точке с координатами х\ =0, = 0. Были также проведены расчеты при том же числе Рейнольдса с включением уравнения (12) в локальные СЛАУ, в результате машинное время счета каждого варианта уменьшилось примерно в 3.6 раза. Наибольшее ускорение сходимости AF = 162.28 было получено при Kmgr = 4, к = 9. Из описанных двух серий расчетов рассматриваемой задачи гидродинамики по методам КНК15 и КНК12 следует, что наилучшим для ускорения сходимости итерационного процесса при умеренных числах Рейнольдса является число невязок к = 9 в описанном выше варианте алгоритма Крылова.

Для случая Re = 100 было проведено сравнение точности методов КНК12 и КНК15 путем вычисления величины

Svi,mb = max|v1,KHKm, (0.5, X2m) - ^1,Ghia(0.5, X2m)|,

" m' b 1

где t>1,Ghia(0.5, X2m) —значения ы, полученные в [36], а ^1,кнкть (0.5, Ж2т) — значения ^1, полученные методом КНКть при = 12 или 15. Было получено, что 5v1}12 = 0.01722, 8v1t15 = 0.01402. То есть применение метода КНК15 позволило уменьшить ошибку 5v1,mb на 18.58 % по сравнению с методом КНК12.

Как отмечалось в [40], в литературе имеется много очень близких друг к другу численных результатов в случае Re ^ 1000, но численные решения начинают заметно отличаться друг от друга при Re > 1000. Ниже мы приводим пример расчета эталонной задачи по методу КНК при числе Рейнольдса Re = 2500. Перед выполнением этого расчета были найдены оптимальные значения £opt = 0.08 и ^opt = 1.75 для двухпараметричекого пре-добуславливателя. На рис. 8 представлены некоторые результаты численных расчетов течения вязкой несжимаемой жидкости в квадратной каверне для рассматриваемого числа Рейнольдса. Расчеты выполнены на равномерной сетке 256 х 256 ячеек. Сплошная линия на рис. 8 (справа) —результат расче-

Рис. 8. Решение эталонной задачи по методу КНК15 при Re = 2500: слева —картина линий тока, справа — профиль составляющей скорости v1 вдоль линии х1 = 0.5 (значки "х" и малые серые точки — результаты работ, соответственно, [38] и [39])

[Figure 8. The benchmark problem solution by the CLS1b method at Re = 2500: on the left is the streamline pattern; on the right is the profile of the velocity component v1 along the line x1 = 0.5 (the symbols "х" and small grey circles are the results of the works [38] and [39],

respectively)]

та по вышеописанному методу КНК15. Полученные результаты сравнивались с результатами работы [38] и более поздней работы [39]. Отметим, что в [38] отсутствуют табличные данные о ^(0.5, Ж2) в интервалах 0.2 < Х2 < 0.5 и 0.5 < х2 < 0.9. Из рис. 8 (справа) видно, что результат расчета по методу КНК15 хорошо согласуется с результатами работ [38,39]. При Ре = 1000 вблизи левого верхнего угла каверны вихрь весьма слабый и плохо выявляется в расчетах [15,33,37,38,40], а при Ре = 2500 он существенно интенсивнее, чем при Ре = 1000, и хорошо передан в расчетах, приведенных на рис. 8 (слева) и в работах [35,40].

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

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

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

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

Такая комбинация способов ускорения итерационных процессов может быть реализована и при применении других численных итерационных методов решения УЧП.

Предложенный выше способ расчета параметров течения в выходном сечении канала показал свою высокую эффективность при его применении в рамках метода КНК. Поэтому представляет значительный интерес разработка аналогичных алгоритмов для их применения в конечно-разностных методах и методах конечного объема.

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

Конкурирующие интересы. Заявляем, что в отношении авторства и публикации этой статьи конфликта интересов не имеем.

Авторский вклад и ответственность. Оба автора принимали участие в разработке концепции статьи и в написании рукописи. Авторы несут полную ответственность за предоставление окончательной рукописи в печать. Окончательная версия рукописи была одобрена обоими авторами.

Финансирование. Результаты исследования получены в рамках выполнения государственных заданий Минобрнауки России № AAAA-A17-117030610134-9 и № AAAA-A17-117030610136-3.

Библиографический список

1. Ferziger J. H., Peric M., Street R. L. Computational Methods for Fluid Dynamics. Cham: Springer, 2020. xviii+596 pp. https://doi.org/10.1007/978-3-319-99693-6.

2. Reddy J. N., Gartling D. K. The Finite Element Method in Heat Transfer and Fluid Dynamics. Boca Raton: CRC Press, 2010. xxiii+501 pp. https://doi.org/10.1201/ 9781439882573.

3. Moukalled F., Mangani L., Darwish M. The Finite Volume Method in Computational Fluid Dynamics. Heidelberg: Springer, 2016. xxiii+791 pp. https://doi.org/10.1007/ 978-3-319-16874-6.

4. Jiang B. N. The Least-Squares Finite Element Method: Theory and Applications in Computational Fluid Dynamics and Electromagnetics. Berlin: Springer, 1998. xvi+418 pp. https://doi.org/10.1007/978-3-662-03740-9.

5. Kim N., Reddy J. N. A spectral/hp least-squares finite element analysis of the Carreau-Yasuda fluids// Int. J. Numer. Meth. Fluids, 2016. vol.82, no. 9. pp. 541-566. https:// doi.org/10.1002/fld.4230.

6. Ranjan R., Chronopoulos A. T., Feng Y. Computational algorithms for solving spectral/hp stabilized incompressible flow problems// J. Math. Res., 2016. vol.8, no. 4. pp. 21-39. https://doi.org/10.5539/jmr.v8n4p21.

7. Ramsak M., Skerget L. A subdomain boundary element method for high-Reynolds laminar flow using stream function-vorticity formulation // Int. J. Numer. Meth. Fluids, 2004. vol. 46, no. 8. pp. 815-847. https://doi.org/10.1002/fld.776.

8. Zhang X., An X., Chen C. S. Local RBFs based collocation methods for unsteady Navier-Stokes equations// Adv. Appl. Math. Mech., 2015. vol.7, no. 4. pp. 430-440. https://doi. org/10.4208/aamm.2013.m337.

9. Плясунова А. В., Слепцов А. Г. Коллокационно-сеточный метод решения нелинейных параболических уравнений на подвижных сетках // Моделирование в механике, 1987. Т. 18, №4. С. 116-137.

10. Исаев В. И., Шапеев В. П. Варианты метода коллокаций и наименьших квадратов повышенной точности для численного решения уравнений Навье-Стокса // Ж. вычисл. матем. и матем. физ, 2010. Т. 50, №10. С. 1758-1770.

11. Исаев В. И., Шапеев В. П. Метод коллокаций и наименьших квадратов повышенной точности для решения уравнений Навье-Стокса // Докл. РАН, 2012. Т. 442, №4. С. 442445.

12. Shapeev V.P., Vorozhtsov E.V. Symbolic-numeric implementation of the method of collocations and least squares for 3D Navier-Stokes equations / Computer Algebra in Scientific Computing. CASC 2012/ Lecture Notes in Computer Science, 7442. Heidelberg: Springer, 2012. pp. 321-333. https://doi.org/10.1007/978-3-642-32973-9_27.

13. Shapeev V. P., Vorozhtsov E. V. Symbolic-numerical optimization and realization of the method of collocations and least residuals for solving the Navier-Stokes equations / Computer Algebra in Scientific Computing. CASC 2016/ Lecture Notes in Computer Science, 9890. Cham: Springer, 2016. pp. 473-488. https://doi.org/10.1007/978-3-319-45641-6_30.

14. Исаев В. И., Шапеев В. П., Еремин С. А. Исследование свойств метода коллокации и наименьших квадратов решения краевых задач для уравнения Пуассона и уравнений Навье-Стокса// Вычислит. технологии, 2007. Т. 12, №3. С. 1-19.

15. Vorozhtsov E. V., Shapeev V. P. On the efficiency of combining different methods for acceleration of iterations at the solution of PDEs by the method of collocations and least residuals// Appl. Math. Comput., 2019. vol.363. pp. 1-19. https://doi.org/10.1016/ j.amc.2019.124644.

16. Федоренко Р. П. О скорости сходимости одного итерационного процесса // Ж. вычисл. матем. и матем. физ., 1964. Т. 4, №3. С. 559-564.

17. Крылов А. Н. О численном решении уравнения, которым в технических вопросах определяются частоты малых колебаний материальных систем // Изв. АН СССР, VII сер., Отд. матем. и естеств. наук, 1931. №4. С. 491-539.

18. Saad Y. Numerical Methods for Large Eigenvalue Problems. Philadelphia: Society for Industrial and Applied Mathematics, 2011. xvi+276 pp. https://doi.org/10.1137/ 1.9781611970739.

19. Kirkpatrick M. P., Armfield S. W., Kent J. H. A representation of curved boundaries for the solution of the Navier-Stokes equations on staggered three-dimensional Cartesian grid// J. Comput. Phys., 2003. vol.184, no. 1. pp. 1-36. https://doi.org/10.1016/ s0021-9991(02)00013-x.

20. Abbassi H., Turki S., Nasrallah Ben S. Channel flow past bluff body: outlet boundary condition, vortex shedding and effects of buyoancy // Comput. Mech., 2002. vol. 28, no. 1. pp. 10-16. https://doi.org/10.1007/s004660100261.

21. Erturk E. Numerical solutions of 2-D steady incompressible flow over a backward-facing step, Part I: High Reynolds number solution // Comput. Fluids, 2008. vol. 37, no. 6. pp. 633-655. https://doi.org/10.1016Zj.compfluid.2007.09.003.

22. Шапеев В. П., Ворожцов Е. В., Исаев В. И., Идимешев С. В. Метод коллокаций и наименьших невязок для трехмерных уравнений Навье-Стокса // Выч. мет. программирование, 2013. Т. 14, №3. С. 306-322.

23. Demmel J. W. Applied Numerical Linear Algebra. Philadelphia: Society for Industrial and Applied Mathematics, 1997. viii+418 pp. https://doi.org/10.1137/1.9781611971446.

24. Saad Y. Iterative Methods for Sparse Linear Systems. Philadelphia: Society for Industrial and Applied Mathematics, 2003. xvi+447 pp. https://doi.org/10.1137/ 1.9780898718003.

25. Wesseling P. An Introduction to Multigrid Methods. Chichester: John Wiley & Sons, 1992. vi+284 pp.

26. Chiu P. H., Sheu T. W. H., Lin R. K. An effective explicit pressure gradient scheme implemented in the two-level non-staggered grids for incompressible Navier-Stokes equations // J. Comput. Phys., 2008. vol.227, no. 8. pp. 4018-4037. https://doi.org/10.1016/j.jcp. 2007.12.007.

27. Shapeev V. P., Vorozhtsov E. V. CAS application to the construction of the collocations and least residuals method for the solution of 3D Navier-Stokes equations / Computer Algebra in Scientific Computing. CASC 2014 / Lecture Notes in Computer Science, 8136. Heidelberg: Springer, 2013. pp. 381-392. https://doi.org/10.1007/978-3-319-10515-4_31.

28. Gartling D. K. A test problem for outflow boundary conditions — flow over a backward-facing step// Int. J. Numer. Methods Fluids, 1990. vol.11, no. 7. pp. 953-967. https:// doi.org/10.1002/fld.1650110704.

29. Воеводин В. В. Вычислительные основы линейной алгебры. М.: Наука, 1977. 303 с.

30. Мартыненко С. И. Совершенствование вычислительных алгоритмов для решения уравнений Навье-Стокса на структурированных сетках // Вестн. МГТУ им. Н.Э. Баумана. Сер. Естественные науки, 2008. №2. С. 78-94.

31. Rouizi Y., Favennec Y., Ventura J., Petit D. Numerical model reduction of 2D steady incompressible laminar flows: Application on the flow over a backward-facing step // J. Comput. Phys., 2009. vol. 228, no. 6. pp. 2239-2255. https://doi.org/10.1016/j.jcp.2008.12.001.

32. Parsani M., Ghorbaniasl G., Lacor C. Analysis of the implicit LU-SGS algorithm for 3rd-and 4th-order spectral volume scheme for solving the steady Navier-Stokes equations // J. Comput. Phys., 2011. vol.230, no. 19. pp. 7073-7085. https://doi.org/10.1016/j.jcp. 2011.05.026.

33. Roberts N. V., Demkowicz L., Moser R. A discontinuous Petrov-Galerkin methodology for adaptive solutions to the incompressible Navier-Stokes equations// J. Comput. Phys., 2015. vol. 301. pp. 456-483. https://doi.org/10.1016/jjcp.2015.07.014.

34. Bustamante C. A., Power H., Florez W. F. A global meshless collocation particular solution method for solving the two-dimensional Navier-Stokes system of equations // Comput. Math. Appl., 2013. vol.65, no. 12. pp. 1939-1955. https://doi.org/10.1016/j.camwa.2013.04. 014.

35. Shapeev A. V., Lin P. An asymptotic fitting finite element method with exponential mesh refinement for accurate computation of corner eddies in viscous flows // SIAM J. Sci. Comput., 2009. vol.31, no. 3. pp. 1874-1900. https://doi.org/10.1137/080719145.

36. Ghia U., Ghia K. N., Shin C. T. High-Re solutions for incompressible flow using the Navier-Stokes equations and a multigrid method// J. Comput. Phys., 1982. vol.48, no. 3. pp. 387-411. https://doi.org/10.1016/0021-9991(82)90058-4.

37. Botella O., Peyret R. Benchmark spectral results on the lid-driven cavity flow// Comput. Fluids, 1998. vol.27, no. 4. pp. 421-433. https://doi.org/10.1016/S0045-7930(98) 00002-4.

38. Erturk E., Corke T. C., Gôkçôl C. Numerical solutions of 2-D steady incompressible driven cavity flow at high Reynolds numbers // Int. J. Numer. Methods Fluids, 2005. vol. 48, no. 7. pp. 747-774. https://doi.org/10.1002/fld.953.

39. Vuorinen V., Larmi M., Schlatter P., Fuchs L., Boersma B. J. A low-dissipative, scale-slective discretization scheme for the Navier-Stokes equations // Comput. Fluids, 2012. vol.70. pp. 195-205. https://doi.org/10.1016/j.compfluid.2012.09.022.

40. Lim R., Sheen D. Nonconforming finite element method applied to the driven cavity flow // Comm. Comput. Phys., 2017. vol.21, no. 4. pp. 1012-1038, arXiv: 1502.04217 [math.NA]. https://doi.org/10.4208/cicp.QA-2016-0039.

Математический институт им. В.А. Стеклова Российской академии наук приступает к работе в рамках Государственного контракта № 13.597.11.0043 по теме «Создание электронного архива выпусков научных журналов по тематическому направлению «Математика, физика, информационные технологии». Архив будет размещен на Общероссийском портале Math-Net.Ru.

Предполагается пополнить коллекцию Math-Net.Ru архивами ряда ведущих журналов по математике, физике и информационным технологиям, а также материалами научных мероприятий.

Проект представлен в социальных сетях: ¥ О н Math-Net.Ru.

Vestn. Samar. Gos. Tekhn. Univ., Ser. Fiz.-Mat. Nauki

[J. Samara State Tech. Univ., Ser. Phys. Math. Sci.], 2020, vol. 24, no. 3, pp. 542-573

d https://doi.org/10.14498/vsgtu1758

ISSN: 2310-7081 (online), 1991-8615 (print)

MSC: 76D05, 76D17, 76G25, 76M25, 76M30

A divergence-free method of collocations and least squares for the computation of incompressible fluid flows and its efficient implementation

© E. V. Vorozhtsov1, V. P. Shapeev1,2

1 Khristianovich Institute of Theoretical and Applied Mechanics,

Siberian Branch of the Russian Academy of Sciences,

4/1, Institutskaya st., Novosibirsk, 630090, Russian Federation.

2 Novosibirsk National Research University,

2, Pirogov st., Novosibirsk, 630090, Russian Federation.

Abstract

The problem of the acceleration of the iterative process of numerical solution by the collocation and least squares (CLS) method of boundary value problems for partial differential equations is considered. For its solution, it is proposed to apply simultaneously three ways to accelerate the iterative process: preconditioner, multigrid algorithm, and Krylov method. A method for finding the optimal values of the parameters of the two-parameter preconditioner is proposed. The use of the found preconditioner significantly accelerates the iterative process. The influence on the iterative process of all three ways of its acceleration is investigated: each separately, and also at their combined application. The application of the algorithm using Krylov subspaces gives the greatest contribution. The combined use of all three ways to speed up the iteration process of solving boundary value problems for two-dimensional Navier-Stokes equations has reduced the CPU time up to 362 times as compared with the case when only one of them, the preconditioner, was applied.

Keywords: preconditioning, Krylov subspaces, multigrid algorithms, Navier-Stokes equations, the method of collocations and least squares.

Received: 25th November, 2019 / Revised: 29th July, 2020 / Accepted: 24th August, 2020 / First online: 21st September, 2020

Research Article

3 ©® The content is published under the terms of the Creative Commons Attribution 4.0 International License (http://creativecommons.org/licenses/by/4.0/) Please cite this article in press as:

Vorozhtsov E. V.,S hapeevV. P. A divergence-free method of collocations and least squares for the computation of incompressible fluid flows and its efficient implementation, Vestn. Samar. Gos. Tekhn. Univ., Ser. Fiz.-Mat. Nauki [J. Samara State Tech. Univ., Ser. Phys. Math. Sci.], 2020, vol. 24, no. 3, pp. 542-573. https://doi.org/10.14498/vsgtu1758 (In Russian). Authors' Details:

Evgenii V. Vorozhtsov https://orcid.org/0000-0003-2753-8399

Dr. Phys. & Math. Sci., Professor; Leading Research Scientist; Lab. of the Physics of Rapid Processes; e-mail: vorozh@itam.nsc.ru

Vasily P. Shapeev © https://orcid.org/0000-0001-6761-7273

Dr. Phys. & Math. Sci., Professor; Chief Research Scientist; Lab. of Thermomechanics and Strength of New Materials; e-mail: shapeev@itam.nsc.ru

570

© Samara State Technical University

Competing interests. We declare that we have no conflict of interest regarding the authorship and publication of this article.

Authors' contributions and responsibilities. Both authors participated in the development of the concept of the article and in the writing of the manuscript. Authors are solely responsible for submitting the final manuscript to the press. The final version of the manuscript was approved by both authors.

Funding. The results have been obtained in the framework of the state tasks of the Ministry of education of Russia no. AAAA-A17-117030610134-9 and no. AAAA-A17-117030610136-3.

References

1. Ferziger J. H., Peric M., Street R. L. Computational Methods for Fluid Dynamics. Cham, Springer, 2020, xviii+596 pp. https://doi.org/10.1007/978-3-319-99693-6.

2. Reddy J. N., Gartling D. K. The Finite Element Method in Heat Transfer and Fluid Dynamics. Boca Raton, CRC Press, 2010, xxiii+501 pp. https://doi.org/10.1201/ 9781439882573.

3. Moukalled F., Mangani L., Darwish M. The Finite Volume Method in Computational Fluid Dynamics. Heidelberg, Springer, 2016, xxiii+791 pp. https://doi.org/10.1007/ 978-3-319-16874-6.

4. Jiang B. N. The Least-Squares Finite Element Method: Theory and Applications in Computational Fluid Dynamics and Electromagnetics. Berlin, Springer, 1998, xvi+418 pp. https://doi.org/10.1007/978-3-662-03740-9.

5. Kim N., Reddy J. N. A spectral/hp least-squares finite element analysis of the Carreau-Yasuda fluids, Int. J. Numer. Meth. Fluids, 2016, vol.82, no. 9, pp. 541-566. https://doi. org/10.1002/fld.4230.

6. Ranjan R., Chronopoulos A. T., Feng Y. Computational algorithms for solving spectral/hp stabilized incompressible flow problems, J. Math. Res., 2016, vol.8, no. 4, 21-39. https://doi.org/10.5539/jmr.v8n4p21.

7. Ramsak M., Skerget L. A subdomain boundary element method for high-Reynolds laminar flow using stream function-vorticity formulation, Int. J. Numer. Meth. Fluids, 2004, vol. 46, no. 8, pp. 815-847. https://doi.org/10.1002/fld.776.

8. Zhang X., An X., Chen C. S. Local RBFs based collocation methods for unsteady Navier-Stokes equations, Adv. Appl. Math. Mech., 2015, vol.7, no. 4, pp. 430-440. https://doi. org/10.4208/aamm.2013.m337.

9. Plyasunova A. V., Sleptsov A. G. Collocation-grid method for solving nonlinear parabolic equations on moving grids, Model. Mekh., 1987, vol. 18, no. 4, pp. 116-137 (In Russian).

10. Isaev V. I., Shapeev V. P. High-accuracy versions of the collocations and least squares method for the numerical solution of the Navier-Stokes equations, Comput. Math. Math. Phys., 2010, vol. 50, no. 10, pp. 1670-1681. https://doi.org/10.1134/S0965542510100040.

11. Isaev V. I., Shapeev V. P. High-order accurate collocations and least squares method for solving the Navier-Stokes equations, Dokl. Math., 2012, vol.85, no. 4, 71-74. https://doi.org/10.1134/S1064562412010255.

12. Shapeev V.P., Vorozhtsov E.V. Symbolic-numeric implementation of the method of collocations and least squares for 3D Navier-Stokes equations, In: Computer Algebra in Scientific Computing. CASC 2012, Lecture Notes in Computer Science, 7442. Heidelberg, Springer, 2012, pp. 321-333. https://doi.org/10.1007/978-3-642-32973-9_27.

13. Shapeev V. P., Vorozhtsov E. V. Symbolic-numerical optimization and realization of the method of collocations and least residuals for solving the Navier-Stokes equations, In: Computer Algebra in Scientific Computing. CASC 2016, Lecture Notes in Computer Science, 9890. Cham, Springer, 2016, pp. 473-488. https://doi.org/10.1007/978-3-319-45641-6_30.

14. Isaev V. I., Shapeev V. P., Eremin S. A. Investigation of the properties of the method of collocations and least squares for solving the boundary-value problems for the Poisson equation and the Navier-Stokes equations, Vychisl. Tekhnol., 2007, vol.12, no. 3, pp. 1-19 (In Russian).

15. Vorozhtsov E. V., Shapeev V. P. On the efficiency of combining different methods for acceleration of iterations at the solution of PDEs by the method of collocations and least residuals, Appl. Math. Comput., 2019, vol.363, pp. 1-19. https://doi.org/10.1016/ j.amc.2019.124644.

16. Fedorenko R. P. The speed of convergence of one iterative process, U.S.S.R. Comput. Math. Math. Phys., 1964, vol.4, no. 3, pp. 227-235. https://doi.org/10.1016/0041-5553C64) 90253-8.

17. Krylov A. N. On the numerical solution of the equation, which determines in technological questions the frequencies of small oscillations of material systems, Izv. Akad. Nauk SSSR, Otd. Mat. Estest. Nauk, VII Ser., 1931, no. 4, pp. 491-539 (In Russian).

18. Saad Y. Numerical Methods for Large Eigenvalue Problems. Philadelphia, Society for Industrial and Applied Mathematics, 2011, xvi+276 pp. https://doi.org/10.1137/ 1.9781611970739.

19. Kirkpatrick M. P., Armfield S. W., Kent J. H. A representation of curved boundaries for the solution of the Navier-Stokes equations on staggered three-dimensional Cartesian grid, J. Comput. Phys., 2003, vol.184, no. 1, pp. 1-36. https://doi.org/10.1016/ s0021-9991(02)00013-x.

20. Abbassi H., Turki S., Nasrallah Ben S. Channel flow past bluff body: outlet boundary condition, vortex shedding and effects of buyoancy, Comput. Mech., 2002, vol.28, no. 1, pp. 10-16. https://doi.org/10.1007/s004660100261.

21. Erturk E. Numerical solutions of 2-D steady incompressible flow over a backward-facing step, Part I: High Reynolds number solution, Comput. Fluids, 2008, vol. 37, no. 6, pp. 633-655. https://doi.org/10.1016/j.compfluid.2007.09.003.

22. Shapeev V. P., Vorozhtsov E. V., Isaev V. I., Idimeshev S. V. The method of collocations and least residuals for three-dimensional Navier-Stokes equations, Num. Meth. Prog., 2013, vol. 14, no. 3, pp. 306-322 (In Russian).

23. Demmel J. W. Applied Numerical Linear Algebra. Philadelphia, Society for Industrial and Applied Mathematics, 1997, viii+418 pp. https://doi.org/10.1137Z1.9781611971446.

24. Saad Y. Iterative Methods for Sparse Linear Systems. Philadelphia, Society for Industrial and Applied Mathematics, 2003, xvi+447 pp. https://doi.org/10.1137/ 1.9780898718003.

25. Wesseling P. An Introduction to Multigrid Methods. Chichester, John Wiley & Sons, 1992, vi+284 pp.

26. Chiu P. H., Sheu T. W. H., Lin R. K. An effective explicit pressure gradient scheme implemented in the two-level non-staggered grids for incompressible Navier-Stokes equations, J. Comput. Phys., 2008, vol.227, no. 8, pp. 4018-4037. https://doi.org/10.1016/jjcp. 2007.12.007.

27. Shapeev V. P., Vorozhtsov E. V. CAS application to the construction of the collocations and least residuals method for the solution of 3D Navier-Stokes equations, In: Computer Algebra in Scientific Computing. CASC 2014, Lecture Notes in Computer Science, 8136. Heidelberg, Springer, 2013, pp. 381-392. https://doi.org/10.1007/978-3-319-10515-4_31.

28. Gartling D. K. A test problem for outflow boundary conditions — flow over a backward-facing step, Int. J. Numer. Methods Fluids, 1990, vol. 11, no. 7, pp. 953-967. https://doi. org/10.1002/fld.1650110704.

29. Voevodin V. V. Vychislitel'nye osnovy lineinoi algebry [Computational Fundamentals of Linear Algebra]. Moscow, Nauka, 1977, 303 pp. (In Russian)

30. Martynenko S. I. Improvement of numerical algorithms for solving the Navier-Stokes equations on structured grids, Herald of the Bauman Moscow State Technical University. Series Natural Sciences, 2008, no. 2, pp. 78-94 (In Russian).

31. Rouizi Y., Favennec Y., Ventura J., Petit D. Numerical model reduction of 2D steady incompressible laminar flows: Application on the flow over a backward-facing step, J. Comput. Phys., 2009, vol. 228, no. 6, pp. 2239-2255. https://doi.Org/10.1016/j.jcp.2008.12.001.

32. Parsani M., Ghorbaniasl G., Lacor C. Analysis of the implicit LU-SGS algorithm for 3rd- and 4th-order spectral volume scheme for solving the steady Navier-Stokes equations, J. Comput. Phys., 2011, vol.230, no. 19, pp. 7073-7085. https://doi.Org/10.1016/j.jcp.2011. 05.026.

33. Roberts N. V., Demkowicz L., Moser R. A discontinuous Petrov-Galerkin methodology for adaptive solutions to the incompressible Navier-Stokes equations, J. Comput. Phys., 2015, vol. 301, pp. 456-483. https://doi.Org/10.1016/j.jcp.2015.07.014.

34. Bustamante C. A., Power H., Florez W. F. A global meshless collocation particular solution method for solving the two-dimensional Navier-Stokes system of equations, Comput. Math. Appl., 2013, vol.65, no. 12, pp. 1939-1955. https://doi.Org/10.1016/j.camwa.2013.04. 014.

35. Shapeev A. V., Lin P. An asymptotic fitting finite element method with exponential mesh refinement for accurate computation of corner eddies in viscous flows, SIAM J. Sci. Comput., 2009, vol.31, no. 3, pp. 1874-1900. https://doi.org/10.1137/080719145.

36. Ghia U., Ghia K. N., Shin C. T. High-Re solutions for incompressible flow using the Navier-Stokes equations and a multigrid method, J. Comput. Phys., 1982, vol. 48, no. 3, pp. 387-411. https://doi.org/10.1016/0021-9991(82)90058-4.

37. Botella O., Peyret R. Benchmark spectral results on the lid-driven cavity flow, Comput. Fluids, 1998, vol.27, no. 4, pp. 421-433. https://doi.org/10.1016/S0045-7930(98)00002-4.

38. Erturk E., Corke T. C., Gokçol C. Numerical solutions of 2-D steady incompressible driven cavity flow at high Reynolds numbers, Int. J. Numer. Methods Fluids, 2005, vol.48, no. 7, pp. 747-774. https://doi.org/10.1002/fld.953.

39. Vuorinen V., Larmi M., Schlatter P., Fuchs L., Boersma B. J. A low-dissipative, scale-slective discretization scheme for the Navier-Stokes equations, Comput. Fluids, 2012, vol. 70, pp. 195-205. https://doi.org/10.1016/j.compfluid.2012.09.022.

40. Lim R., Sheen D. Nonconforming finite element method applied to the driven cavity flow, Comm. Comput. Phys., 2017, vol.21, no. 4, pp. 1012-1038, arXiv: 1502.04217 [math.NA]. https://doi.org/10.4208/cicp.0A-2016-0039.