CC BY
48
Вычислительные технологии
Том 3, № 3, 1998
МЕТОД КОЛЛОКАЦИЙ И НАИМЕНЬШИХ КВАДРАТОВ ДЛЯ УРАВНЕНИЙ НАВЬЕ-СТОКСА*
Л. Г. Семин, В. П. Шапеев Институт теоретической и прикладной механики СО РАН Новосибирск, Россия
Светлой памяти Анатолия Георгиевича Слепцова посвящается
A grid method of collocation and least squares for solving the boundary problem for a set of Navier — Stokes equations is created. The method is implemented on a regular square grid. Formulas of the method are tested on a range of test problems. Convergence of numerical solution is established in computational experiments. The computation of the flow in a rectangular lid-driven cavity is performed. The results of computations by the present method for various Reynolds numbers are compared with the results of other researchers, as well as with the physical experiment. All results are shown to be in good agreement. Formulas of the method are derived using the REDUCE computer algebra system.
Введение
В последнее время благодаря простоте физической интерпретации, математической формы и гибкости численного алгоритма все более широкое распространение получают методы конечных элементов. В настоящей работе представлен сеточный метод коллокаций и наименьших квадратов (КНК) решения краевых задач для уравнений Навье — Стокса, который является одним из вариантов метода конечных элементов. Метод КНК успешно применялся при решении уравнений в частных производных эллиптического и параболического типов [1], хорошо зарекомендовал себя при решении эллиптических задач на адаптивных сетках с малым параметром при старшей производной, где могут возникать пограничные и/или внутренние слои [2]. Для систем уравнений метод впервые использовался нами при решении краевых задач для уравнений Стокса [3].
При конечно-элементной аппроксимации краевой задачи для системы уравнений На-вье — Стокса одна из трудностей состоит в обеспечении хорошей аппроксимации уравнения неразрывности divv = 0 [4] и реализации граничных условий. В предлагаемом здесь методе уравнение неразрывности выполняется точно в каждой ячейке сетки за счет выбора базиса, легко передаются условия на границах. Для определения коэффициентов применяется метод наименьших квадратов, который дает лучше обусловленную по сравнению с другими методами матрицу системы для определения коэффициентов разложения решения по базису “результирующей системы".
* Работа выполнена при финансовой поддержке Российского фонда фундаментальных исследований, грант №96-01-01888.
© Л. Г. Семин, В. П. Шапеев, 1998.
Опыт решения эллиптических задач на анизотропных сетках [2] показал, что при использовании лагранжевых интерполянтов в каждой ячейке сетки и классических согласованных конечных элементов получаются плохо обусловленные системы линейных алгебраических уравнений. Поэтому для аппроксимации были применены кусочно-полиномиальные функции, согласованные в слабом смысле (см. [2]). Использование таких аппроксимаций при решении краевой задачи для уравнений Навье — Стокса несколько затрудняет применение хорошо известных итерационных методов, разработанных для этих задач [5]. Поэтому здесь использован метод итераций по подобластям, с которым можно познакомиться, например, в [6]. Для ускорения сходимости линейных итераций применялся метод, основанный на проекции погрешности на подпространство невязок [7].
С помощью метода КНК проведен расчет течения в прямоугольной каверне с движущейся верхней границей. Задача о течении в каверне интересна тем, что она представляет собой относительно простую модель, на которой можно проверять и сравнивать численные методы. Поэтому ряд исследователей [11-15] и использовали эту задачу для сравнения своих результатов.
1. Постановка задачи
Рассматривается краевая задача для уравнений Навье — Стокса:
- V,д3 - г,д3 - дХ
дх , дх 2 дХз
= , (Х,,Х2) Є П, і = 1, 2,
(1)
= 0, фп = V.
Приближенное решение отыскивается в виде кусочно-полиномиальной функции на регулярной сетке. Вводятся следующие обозначения: к — половина ширины ячейки, (х1т, х2т) — координаты центра т-й ячейки, у1 = (х1 — х1т)/к, у2 = (х2 — х2т)/к — локальные координаты в ячейке, й(у1,у2) = у(х1,х2), д(у1,у2) = кр(х1,х2). Тогда уравнения (1) в локальных переменных примут следующий вид:
д—
+
дд
дУ2) ' дУ] — = 0,
—|ап = и.
ан2!3, і = 1,2,
(2)
(3)
(4)
Пусть известно некоторое приближенное решение (й1,й2,д). Искомое уточненное решение представим в виде и1 = й1 + м1, й2 = й2 + й2, д = д + й Подставив это представление в уравнения (2) и пренебрегая членами второго порядка малости, получим уравнения Озеена:
Аи3 — Кк(й1й^;у1 + й2йj,y2 + й1 щ,у 1 + й2йj,y2) — Rgyj — , ] — 1, 2
где Fj = R(к2fj + к(йlйj,yl + й22) + ) — Айj.
В каждой ячейке уточнение решения отыскивается в виде
(5)
где ф] — базисные функции (векторы с тремя компонентами), т — номер ячейки.
Коэффициенты а]т будут определяться из уравнений коллокаций и условий согласования на границах соседних ячеек или краевых условий на дП.
2. Приближенные уравнения
Базисные функции выбираем так, чтобы компоненты вектора скорости представлялись в виде полиномов второго порядка, а давление выражалось в виде линейной функции, причем требуем, чтобы решение (6) удовлетворяло уравнению неразрывности и = 0. Учитывая это, базисные функции берем в следующем виде:
к 1 2 3 4 5 6 7 8 9 10 11 12
1 0 VI У2 0 0 V2 V22 0 v2 + V2
фк 0 1 1 0 VI 0 2УіУ2 0 V2 V2 -^ 1V2 V2 + V2
0 0 0 0 0 1 2Я-1уі 2Я-1уі 2 1- 2 1- 0 0
Для получения единственного кусочно-полиномиального решения необходимо на границах между ячейками задать условия согласования. В качестве последних рассматриваются условия непрерывности выражений:
дип ТТ ди
ЦП - р + пип, ш + (7)
Здесь ип и Иг — нормальная и касательная к границе ячейки компоненты скорости, п — вектор единичной внешней к ячейке нормали, п — некоторый положительный параметр, с помощью которого можно управлять обусловленностью получающейся системы линейных алгебраических уравнений и влиять на скорость сходимости итераций.
Для определения коэффициентов а]т в представлении (6) был использован метод итераций по подобластям, в котором подобластями являлись каждая из ячеек сетки. При построении уточнения решения в т-й ячейке применялись условия (7). Вообще говоря, эти условия не могут быть удовлетворены на всей границе между ячейками. Удовлетворим им, требуя, чтобы они выполнялись в смысле наименьших квадратов в восьми точках границы ячейки.
Если граница ячейки лежит на границе области, то в двух точках этой границы задается вектор скорости. Кроме того, в левом нижнем углу области задается давление.
К условиям согласования (краевым условиям) добавляются условия коллокации уравнений (5) в четырех точках внутри ячейки. В результате получается следующая система линейных алгебраических уравнений:
12
^2/В1к акт = Рг, I = 1,..., 24. (8)
к=1
В ней первые четыре уравнения получены из условий согласования (7) или краевых условий на нижней границе, следующие четыре — из соответствующих условий на правой границе, уравнения при I = 9, ..., 12 — из условий на верхней, а уравнения при I = 13, ... , 16 — из условий на левой границе. Уравнения при I = 17, ... , 24 получены из условий коллокации. Условия согласования в локальных координатах задаются в точках с координатами (±1, ±£), (±С, ±1), краевые условия — в точках (±1, ±£), (±£, ±1), условия коллокации — в точках (±ш, ±ш). Величины £, ^, ш положительны, не превосходят
единицы, их можно выбирать различным образом для получения хорошо обусловленной матрицы системы линейных алгебраических уравнений.
При построении вариантов описанного метода могут быть полезны найденные формулы. Часть из них приведена в приложении. Элементы матрицы В и вектора правых частей Г, полученные из условий согласования, приведены в приложении 1, а полученные из краевых условий — в приложении 2. В этих формулах а™, а™, ав, аЛ — коэффициенты аг из нижней, правой, верхней и левой ячеек соответственно, Ьг — коэффициенты аг из данной ячейки, взятые с предыдущей итерации.
Элементы матрицы В и вектора правых частей Г, найденные из условий коллокации уравнений (2), громоздки, и поэтому в работе не приводятся.
Система (8) переопределена, в ней 24 уравнения и 12 неизвестных. Для того, чтобы определить, что понимается под решением этой системы, рассмотрим два функционала:
Первый функционал соответствует сумме квадратов невязок уравнений, полученных из условий согласования или краевых условий, второй — сумме квадратов невязок уравнений коллокаций. Решение системы (8) находится из условия минимума этих функционалов, причем минимум 1 берется по коэффициентам а]т, ] = 1, ... , 10 при фиксированных а11т и а12т, а минимум 2 — по а]т, ] = 11,12 при фиксированных а]т, ] = 1, ... , 10.
Таким образом, для определения коэффициентов а]т, ] = 1, ... , 12, имеем систему уравнений вида:
решение которой находим методом прямого исключения.
Для определения решения задачи применяется итерационный метод, в котором решение уточняется для каждой ячейки области в отдельности. На каждой итерации при уточнении решения в m-й ячейке решение в соседних ячейках либо уже уточнялось, либо берется с предыдущей итерации.
Различались несколько типов ячеек: все стороны ячейки находятся внутри области (“внутренние"ячейки); одна сторона ячейки лежит на границе области; две соседние стороны ячейки лежат на границе области. Аналитические выражения для элементов матрицы D и вектора правых частей F системы (10) для всех указанных типов ячеек находились программой, написанной в системе REDUCE [8]. Исходные уравнения, базисные функции, вид решения записаны в ней в качестве входных данных. Координаты точек постановки условий согласования, граничных условий, уравнений коллокаций, параметр п входят в полученные формулы в символьном виде, что позволяет при организации численного счета рассматривать различные варианты, связанные с изменением этих величин, не выводя формулы метода заново. Здесь очевидна польза от применения системы компьютерной алгебры.
Правильность полученных формул, сходимость и работоспособность метода проверялись на ряде тестовых задач с известным точным решением. Среди них следующие.
1. Течение между двумя параллельными пластинами. Разность давлений на входе и выходе равна нулю, верхняя пластина движется, внешние силы отсутствуют.
12
l = 1, ..., 12,
(10)
j=i
2. Течение между двумя параллельными покоящимися пластинами. Давление на входе и выходе разное, внешние силы отсутствуют.
3. Течение между двумя параллельными пластинами. Верхняя пластина движется, давление на входе и выходе разное, внешние силы отсутствуют.
4. Расчетная область — единичный квадрат. В качестве точного решения взяты выражения для (у1,у2,р) в виде тригонометрических функций. Подстановка этого решения в уравнения (1) вносит в правые части члены типа внешних сил и источников, которые обеспечивают движение жидкости, заданное этим решением. Для полученной таким образом системы уравнений с дифференциальными операторами (1) применяем метод КНК для нахождения численного решения.
Во всех перечисленных тестах в процессе итераций численное решение сходилось к точному.
3. Расчет течения в прямоугольной каверне с движущейся верхней границей
0.0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1.0 х 0.0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1.0 х
0.0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1.0 х 0.0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1.0 х
Рис. 1. А = 1, сетка 80x80. Я= 1 (а), 100 (б), 400 (в), 1000 (г).
С помощью метода КНК проведен расчет течения в прямоугольной каверне с движущейся верхней границей. Расчеты проводились для чисел Рейнольдса И=1, 100, 400, 1000. Картина линий тока для различных чисел Рейнольдса в квадратной каверне представлена на рис. 1. Знаком ’+’ здесь обозначен центр вихря. Видно, что при увеличении числа Рейнольдса происходит рост угловых придонных вихрей и смещение центра главного вихря сначала в направлении движения крышки, а затем к центру каверны. Этот эффект подтверждается физическими экспериментами, описанными в работе [10], и расчетами, проведенными в [9, 12 - 15].
0.0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1.0 ^ 0.0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1.0 х
в г
0.001 0.003
0.0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1.0 х 0.0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1.0 х
Рис. 2. А = 0.5, сетка 80x80. Я=1(а), 100 (б), 400 (в), 1000 (г).
Сравнение максимальных и минимальных значений скоростей вдоль вертикальной и горизонтальной средней линии каверны с результатами, опубликованными в [12 - 15], приведено в табл. 1. Все данные приведены для А =1. Видно, что характерные величины профилей скорости, найденные нами, находятся в хорошем соответствии с результатами расчетов других исследователей. На рис. 4 приведены профили компонент скорости У\ и для различных чисел Рейнольдса в квадратной каверне.
Координаты центров вихрей и значения функции тока ф в них при различных числах Рейнольдса для квадратной каверны приведены в табл. 2. По значениям функции тока видно, что вторичные вихри на два порядка слабее главного. Тем не менее предлагаемый в данной работе метод хорошо их передает.
Обозначим через А отношение глубины каверны к ее ширине. На рис. 2, 3 приведены картины линий тока в случаях А = 0.5 и 2, а на рис. 5 показаны соответствующие профили горизонтальной компоненты скорости вдоль вертикальной средней линии каверны при А = 2.
Видно, что в отличие от случаев А = 1 (квадратная каверна) и 0.5 при А = 2 образуется второй вихрь, занимающий всю нижнюю часть каверны даже при малых числах Рейнольдса. Это подтверждается в экспериментах [10], а также расчетами, представленными в [10 - 12]. При увеличении числа Рейнольдса размер второго вихря возрастает, хотя последний остается весьма слабым.
На рисунках буквами и и V обозначены соответственно горизонтальная и вертикальная компоненты скорости и X и У — оси Х\ и х2.
Л. Г. Семин, В. П. Шапеев Таблица 1.
Метод (ссылка на работу) Сетка тт ^2 тт ^2 тах
И= 1
Настоящая работа 40 х 40 -0.2071 —0.1844 0.1837
То же 80 х 80 -0.2074 —0.1846 0.1839
И=100
[12] 128 х 128 -0.2106 -0.2521 0.1786
[13] 129 х 129 -0.21090 -0.24533 0.17527
Настоящая работа 20 х 20 -0.2100 -0.2462 0.1761
То же 40 х 40 -0.2133 -0.2513 0.1782
То же 80 х 80 -0.2125 -0.2497 0.1756
И=400
[13] 129 х 129 -0.32726 -0.44993 0.30203
[14] 64 х 64 -0.32368 -0.44862 0.29925
То же 128 х 128 -0.32751 -0.45274 0.30271
[15] 49 х 49 -0.2960 -0.4051 0.2662
Настоящая работа 20 х 20 -0.2780 -0.3899 0.2477
То же 40 х 40 -0.3041 -0.4174 0.2768
То же 80 х 80 -0.3222 -0.4443 0.2970
И=1000
[12] 256 х 256 -0.3764 -0.5208 0.3665
[13] 129 х 129 -0.38289 -0.51550 0.37095
[14] 64 х 64 -0.37436 -0.51015 0.36364
То же 128 х 128 -0.38511 -0.52280 0.37369
[15] 129 х 129 -0.3689 -0.5037 0.3553
Настоящая работа 40 х 40 -0.3360 -0.4752 0.3182
То же 80 х 80 -0.3590 -0.4901 0.3463
Таблица 2. Координаты центров вихрей и значения функции тока в них в квадратной каверне.
И, Главный вихрь Правый п эидонный вихрь
Значение ф Координаты (х,у) Значение ф Координаты (х,у)
1 100 400 1000 -0.100 -0.103 -0.112 -0.110 (0.5; 0.7625) (0.6125; 0.7375) (0.55; 0.6125) (0.5375; 0.575) 2е-6 1е-5 0.0006 0.0016 (0.9625; 0.0375) (0.95; 0.0625) (0.8875; 0.125) (0.8625; 0.1125)
0.0 0.2 0.4 0.6 0.8 1.0 х 0.0 0.2 0.4 0.6 0.8 1.0 х
0.0 0.2 0.4 0.6 0.8 1.0 X О.о 0.2 0.4 0.6 0.8 1.0 х
Рис. 3. А = 2, сетка 40x80. И,= 1 (а), 100 (б), 400 (в), 1000 (г).
4. Заключение
Для системы стационарных уравнений Навье — Стокса построен сеточный метод КНК на регулярной квадратной сетке. С его помощью проведен расчет течения в прямоугольной каверне с движущейся верхней границей при различных числах Рейнольдса. Сравнение значений горизонтальной компоненты скорости у1 вдоль вертикальной средней линии каверны и вертикальной компоненты скорости у2 вдоль горизонтальной средней линии каверны с результатами, приведенными в работах [9 - 15], показывает хорошее соответствие расчетов по методу КНК с расчетами по методам, предлагаемым в этих работах. Наиболее близкое совпадение значений скорости наблюдается при малых числах Рейнольдса. При больших числах Рейнольдса наилучшее совпадение получено с результатами работы [15]. Расхождение данных составляет менее 3%, причем наши расчеты по методу КНК проводились на сетке 80 х 80, а в работе [15] — на сетке 129 х 129.
Картины линий тока, полученные по методу КНК при различных числах Рейнольдса, и аналогичные картины, приведенные в указанных работах, идентичны. С увеличением числа Рейнольдса при всех значениях параметра А происходит смещение центра главного
Рис. 4. Профили горизонтальной компоненты скорости и вдоль вертикальной средней линии каверны (слева) и вертикальной компоненты скорости V вдоль горизонтальной средней линии каверны (справа), А = 1, И,= 1(1), 100(2), 400(3), 1000(4).
Рис. 5. Профили горизонтальной компоненты скорости и вдоль вертикальной средней линии каверны, А = 2, сетка 40x80. И,= 1(1), 100(2), 400(3), 1000(4).
вихря сначала в направлении движения крышки, а потом к центру каверны. Картины течения с указанными выше деталями наблюдались в физическом эксперименте [10].
Авторы благодарят А. Г. Слепцова за внимание к работе
Приложение 1. Элементы матрицы B и вектора правых частей F переопределенной системы (8), полученных из условий согласования.
Элементы матрицы B
l = 1: B„ = - , B,2 = B,5 = BK = Bl9 = 0, B,3 = n+iС, Bu = 1, Bn = -^(2,
Bis = - f+r, Bi 10 = B, 12 = 2Z, B, i i = -(<? ++)n+2.
l = 2: По сравнению со случаем l = 1 изменились три элемента: B,3 = -n+jZ, B, 10 =
Bl 12 = -2C.
1 = 3: B,1 = B,4 = B12 = B,10 = -n+15 B,3 = -1 B,5 = n+1Z, B,6 = - n+1 ’ B,7 = 2((П+12) Z,
B,8 = П+1 Z ’ Bl9 = - П^+12 ’ B,11 = 2Z ’ B,12 = - (Z +^j)1?+2 •
4: По сравнению со случаем l = 3 изменились четыре элемента: B,5 = - -+1 Z,
B,7 = - 2((П+12) Z; Bis = - n+IZ; B,11 = - 2C •
5: B,1 = B,4 = B,6 = B,8 = 0j B,2 = П+1 > B,3 = П+1Z’ B,5 = 1 B,7 = B,11 = 2C; B,9 = -+I’
p = П z 2 p = (С2+1)П+2
Bl10 = П+1Z , B,1r = П+1 •
6: По сравнению со случаем l = 5 изменились три элемента: B,3 = - -+1 Z, B,7 = B,11 = -2Z. 2
7: B,1 = П+1 ’ B,2 = B,5 = 0j B,3 = 1J B,4 = -П+1Z; B,6 = -П+1; B,7 = -+1 > BIS = -П+Г’
B,9 = П+1 z, B,10 = 2(П+кГ‘) Z • Bi11 = 1£!++1П±2 • B,12 = 2Z.
/ T/T^A/T^TTT/T ТТТ/ГГ1!^ TJ^T^T-TT»^ QTT^TV/T^TT'T'Pl • Hi, -
n+1
- П+1Z ’ B,10 = - 2(П?+12) Z > B,12 = -2Z•
= 8: По сравнению со случаем l = 7 изменились четыре элемента: B,4 = -+1Z, B,g =
= 9: B,1 = П+1, B,2 = B,5 = B,6 = B,9 = 0, B,3 = --+1Z, B,4 = 1, B,7 = -+1Z2, B,s = -+f,
B110 = B,12 = 2Z, B,11 = (c2±±)-±2 •
= 10: По сравнению со случаем l = 9 изменились три элемента: B,3 = -+1Z, B,10 =
B,12 = -2Z •
= 11: B,1 = B,4 = 0, B,2 = B,10 = -+1; B,3 = - 1 B,5 = --+1Z; B,6 = --+1; B,7 = Z;
Bis = -+1Z, B,9 = -J+r, B111 = 2Z, B,12 = I^!±I)lS+2 •
= 12: По сравнению со случаем l = 11 изменились четыре элемента: B,5 = -+1Z,
B,7 = - rП±±Г) Z, BIS = - П+1Z, B,11 = -2Z •
= 13: B,1 = B,4 = B,6 = B,s = 0J B,2 = - -+1 ; B,3 = - -+1Z; B,5 = 1J B,7 = B,11 = 2Z;
p П+2 p П Z 2 b (C 2±1)П±2
B,9 = - -+1, B,10 = - -+1Z , B,12-------------------
+1 i10 +1 i12 +1
14: По сравнению со случаем l = 13 изменились три элемента: B,3 = -+1Z, B,7 =
Bi11 = -2Z •
= 15: B,1 = - П+1, B,2 = B,5 = 0j B,3 = 1J B,4 = -+1Z; B,6 = - -+1; B,7 = - -+
Bis = --j+r, B,9 = -+1Z, b,10 = ^ Z, B,11 = - , B,12 = 2Z •
16: По сравнению со случаем l = 15 изменились четыре элемента: B,4 = - -+1Z,
2 Z Р _____ 2(rl±2,
П+1Z, Pi10 = П+Г
Компоненты вектора F
F1 = {^Z2(Ъ7 + Ъ11 - a7 - al1) + nZ[a3 - Ъ3 - 2(Ъ10 + Ъ12 + al0 + al2)] + П(Ъ1 - Ъ4 + Ъ8 + Ъ11 -
a” - a7 - a7 - a^) + 2Z(a^ + a^ - &10 - Ъ^) - &4 + a7 + 2(bs + Ъц + a7 + a71)}/(n + 1),
F2 = {nZ2(Ъ7 + Ъ11 - a7 - al1) + nZ[Ъ3 - a3 + 2(Ъ10 + Ъ12 + al0 + al2)] + П(Ъ1 - Ъ4 + Ъ8 + Ъ11 -
a7 - a7 - a7 - a71) - 2Z(a7^ + a72 - Ъю - Ъ^) - Ъ4 + a7 + 2(Ъ8 + Ъц + a7 + a71)}/(n + 1),
F3 = -{[(Ъ5 - a5 + 2(Ъ7 + Ъ11 + a7 + al1))Z - (Ъ9 + Ъ12 - a9 - al2)Z2 - Ъ2 - Ъ3 - Ъ10 - Ъ12 + a2 - a3 + a10 + al2]n - (Ъ6 - a6)R +2Z (Ъ8 + Ъ11 - a8 - al1 +2(Ъ7 - a7)) - Ъ3 + a3 +2(Ъ9 - Ъ12 + ag - al2)}/(П + 1)J
F4 — {[(Ъ5 - a7 + 2(Ъ7 + Ъц + a7 + a71))Z + (Ъ9 + Ъ^ - a7 - a72)Z + Ъ2 + Ъз + Ъю + Ъ^ - a7 + a7 - a7o -a 12]n + (Ъб - a7)R + 2Z(Ъ8 + Ъц - a7 - a71 + 2(Ъ7 - a7)) + Ъз - a7 + 2(-Ъ9 + Ъ^ - a7 + a72)}/(n + 1), F5 = -{nZ2(Ъ10 + Ъ12 - al0 - al2) + nZ [Ъ3 - a3 + 2(Ъ7 + Ъ11 + a7 + al1)] + П(Ъ2 + Ъ5 + Ъ9 +
Ъ12 - a33 + a!3 - a7 - a^) + 2Z(Ъ7 + Ъц - a33 - a71) + Ъ5 - a!3 + 2(Ъ9 + Ъ12 + a7 + a72)}/(n + 1),
F6 = -{nZ2(Ъ10 + Ъ12 - a!0 - a!2) - nZ[Ъ3 - a3 + 2(Ъ7 + Ъ11 + a7 + al1)] + П(Ъ2 + Ъ5 + Ъ9 +
Ъ12 - a7 + a!7 - a7 - a^) - 2Z(Ъ7 + Ъц - a77 - a^) + Ъ5 - a!7 + 2(Ъ9 + Ъ12 + a7 + a72)}/(n + 1),
F7 = {[(Ъ4 - a4 - 2(Ъ10 + Ъ12 + a!0 + al2))Z - (Ъ8 + Ъ11 - a8 - al1)Z2 - Ъ1 - Ъ3 - Ъ7 - Ъ11 + al - a3 + a7 + a11 ]П + (Ъ6 - a6)R - 2Z(Ъ9 + Ъ12 - ag - al2 + 2(Ъ10 - al0)) - Ъ3 + a3 + 2(Ъ8 - Ъ11 + a8 - al1)}/(П + 1)J
FS = {[-(Ъ4 - a4 - 2(Ъ10 + Ъ12 + al0 + al2))Z - (Ъ8 + Ъ11 - a8 - al1)Z2 - Ъ1 - Ъ3 - Ъ7 - Ъ11 + al - a3 +
a7 + a^n + (Ъб - a7)R+2Z (Ъ9 + Ъ12 - a7 - a^ + 2(Ъю - a^)) - Ъз + a7 + 2(Ъ8 - Ъц + a7 - a71 )}/(n + 1),
F9 = -{nZ2(Ъ7 + Ъ11 - a7 - a11) + nZ[a3 - Ъ3 + 2(Ъ10 + Ъ12 + al0 + a12)] + П(Ъ1 + Ъ4 + Ъ8 + Ъ11 -
a1 + a7 - a7 - a^) + 2Z(Ъю + Ъ12 - a1o - a^) + Ъ4 - a7 + 2(Ъв + Ъц + aS + a^)}/(n + 1),
F10 = -{nZ2 (Ъ7 + Ъ11 - a7 - al1) + nZ [Ъ3 - a3 - 2(Ъ10 + Ъ12 + al0 + al2)] + П(Ъ1 + Ъ4 + Ъ8 +
Ъц - a1 + a7 - a7 - a^) - 2Z(Ъю + Ъ^ - a1o - a12) + Ъ4 - a4 + 2(Ъв + Ъц + aS + a^)}/(n + 1),
F11 = {[(Ъ5 - a5 - 2(Ъ7 + Ъ11 + a7 + al1))Z - (Ъ9 + Ъ12 - a9 - al2)Z2 - Ъ2 + Ъ3 - Ъ10 - Ъ12 + a2 + a3 + a10 + al2]n +(Ъ6 - ag)R- 2Z (Ъ8 + Ъ11 - a8 - a11 + 2(Ъ7 - a7))+ Ъ3 - a3 + 2(Ъ9 - Ъ12 + a9 - a12)}/(n +1),
F12 = { [(Ъ5 - a5 - 2(Ъ7 + Ъ11 +a7 + a11 ))Z + (Ъ9 + Ъ12 - a9 - al2)Z2 +Ъ2 - Ъ3 +Ъ10 + Ъ12 - a2 - a3 - a10 -
a12]П - (Ъ6 - a6)R - 2Z(Ъ8 + Ъ11 - a8 - al1 + 2(Ъ7 - a7)) - Ъ3 + a3 + 2(-Ъ9 + Ъ12 - a9 + a12)}/(П + 1),
F13 = {nZ2(Ъю + Ъ12 - a?0 - a?2) + nZ[Ъз - ag - 2(Ъ7 + Ъц + a? + a^)] + n^2 - Ъ5 + Ъ9 + Ъ12 - a? - a? - a? - a?2) - 2Z(Ъ7 + Ъц - a? - a?1) - Ъ5 + a? + 2(Ъ9 + Ъ12 + a? + a?2)}/(n + 1),
F14 = {nZ2(Ъю + Ъ12 - a?o - a?2) - nZ[Ъз - a? - 2(Ъ7 + Ъц + a? + a^)] + n^2 - Ъ5 + Ъ9 +
Ъ12 - a? - a? - a? - a?2) + 2Z(Ъ7 + Ъц - a? - a?1) - Ъ5 + a? + 2(Ъ9 + Ъ12 + a? + a?2)}/(n + 1),
F15 = {[-(Ъ4 - a4 + 2(Ъ10 + Ъ12 + al0 + al2))Z + (Ъ8 + Ъ11 - a8 - a11 )Z2 + Ъ1 - Ъз + Ъ7 + Ъ11 - al - a3 -a? - ^Jn + (Ъ6 - a6)R- 2Z(Ъ9 + Ъ12 - a9 - al2 + 2(Ъ10 - al0)) - Ъз + a3 - 2(Ъ8 - Ъ11 + a8 - al1)}/(n +1),
F16 = {[(Ъ4 - a? + 2(Ъю + Ъ12 + a'Jo + a^X + (Ъв + Ъц - a" - a^Z2 + Ъ1 - Ъз + Ъ7 + Ъц - al'1 - a" - a? -a11 ]n + (Ъ6 - a6)R + 2Z(Ъ9 + Ъ12 - a9 - a!2 + 2(Ъ10 - al0)) - Ъз + a3 - 2(Ъ8 - Ъ11 + a8 - al1)}/(n + 1)
Приложение 2. Элементы матрицы B и вектора правых частей F переопределенной системы (8), полученных из краевых условий.
Элементы матрицы B
l = 1: Bi1 = Bi8 = 1 Bi2 = Bi5 = Bi6 = Bi9 = ° Р,з = - С Bi4 = -1, Bi7 = ^
Bi10 = Bi12 = -2^ Bi11 = С2 + 1
l = 2: По сравнению со случаем l = 1 изменились три элемента: В,з = С, B,10 = B,12 = 2С•
l = 3: Bi1 = Bi4 = Bi6 = BiS = ° Bi2 = В,з = Bi10 = 1, Bi5 = -С Bi7 = Bi11 = -2C,
Bi9 = С ^ Bi12 = С 2 + 1
l = 4: По сравнению со случаем l = 3 изменились три элемента: В,5 = С, В,7 = В,11 = 2£-
l = 5: В,1 = В,з = В,7 = 1, Bi2 = Bi5 = Bi6 = Bi9 = ° Bi4 = -С, BiS = С2, Bi10 = Bi12 = 2C,
Bi11 = C2 + 1
l = 6: По сравнению со случаем l = 5 изменились три элемента: В,4 = С, В,10 = В,12 = -2С-
l = 7: В,1 = Bi4 = Bi6 = BiS = ° Bi2 = Bi5 = Bi9 = 1, В,з = С, Bi7 = Bi11 = 2C, Bi10 = С2,
Bi12 = C2 + 1
l = 8: По сравнению со случаем l = 7 изменились три элемента: В,з = -С, В,7 = В,11 =
- 2С-
l = 9: Bi1 = Bi4 = BiS = 1, Bi2 = Bi5 = Bi6 = Bi9 = 0, В,з = -С, Bi7 = ^ Bi10 = Bi12 = 2С,
Bi11 = С2 + 1
l = 10: По сравнению со случаем l = 9 изменились три элемента: В,з = С, В,10 = В,12 =
-2С-
l = 11: Bi1 = Bi4 = Bi6 = BiS = ° Bi2 = Bi10 = 1 В,з = -1, Bi5 = -С, Bi7 = Bi11 = 2С,
Bi9 = С2, Bi12 = С2 + 1
l = 12: По сравнению со случаем l = 11 изменились три элемента: В,5 = С, В,7 = В,11 = -2С-
l = 13: В,1 = В,7 = 1, В,2 = В,5 = В,6 = В,9 = 0, В,з = -1, В,4 = -С, B,s = С2,
В,10 = В,12 = -2С, В,11 = С2 + 1
l = 14: По сравнению со случаем l = 13 изменились три элемента: В,4 = С, В,10 = В,12 = 2С
l = 15: В,1 = В,4 = В,6 = В,8 = 0, В,2 = В,9 = 1, В,з = С, В,5 = -1 В,7 = В,11 = -2С,
В,10 = С2, В,12 = С2 + 1
l = 16: По сравнению со случаем l = 15 изменились три элемента: В,з = -С, В,7 = В,11 = 2С
Компоненты вектора F
F1 = (Ъз + 2(Ъю + Ъ12))С - (С2 + 1)Ъц - Ъ7С2 - Ъ1 + Ъ4 - Ъв + С/1(-С, -1),
F2 = - [(Ъз + 2(Ъю + Ъ12))С + (С2 + 1)Ъц + Ъ7С2 + Ъ1 - Ъ4 + Ъв - [Л(С, -1)],
Fз = (Ъ5 + 2(Ъ7 + ЪИ))С - (С2 + 1)Ъ12 - Ъ9С2 - Ъ2 - Ъз - Ъю + Сг(-С, -1),
F4 = -[(Ъ5 + 2(Ъ7 + Ъц))С + (С2 + 1)Ъ12 + Ъ9С2 + Ъ2 + Ъз + Ъю - U^2(С, -1)],
F5 = - [(2(610 + ъ 12) - Ъ4)С + (С2 + 1)Ъц + ЪвС2 + Ъ1 + Ъз + Ъ7 - [/1(1, -С)],
F6 = (2(610 + 612) - 64)С - (С2 + 1)611 - ЪвС2 - 61 - Ъз - 67 + [1 (1,С),
F7 = -[(2(67 + Ъп) + Ъз)С + (С2 +1)612 + ЪюС2 + 62 + 65 + 69 - [/2(1, -С)],
Fs = (2(67 + 611) + Ъз)С - (С2 + 1)612 - ЪюС2 - 62 - 65 - 69 + [/2(1,С),
F9 = (Ъз - 2(Ъю + ъ 12))С - (С2 +1)611 - 67с2 - 61 - 64 - bs + [/1 (-С, 1),
F10 = - [(Ъз - 2(Ъю + Ъ12))С + (С2 + 1)611 + 67С2 + 61 + 64 + bs - [71(С, 1)],
F11 = (65 - 2(67 + Ъп))С - (С2 + 1)612 - 69С2 - 62 + Ъз - Ъю + [^2( С, 1),
F12 = -[(65 - 2(67 + ъп))с + (С2 +1)612 + 69с2 + 62 - Ъз + Ъю - [/2(С, 1)],
Fw = (2(Ъю + 612) + 64)С - (С2 + 1)611 - bsC2 - 61 + Ъз - 67 + U1 (-1, -С),
F14 = -[(2(610 + 612) + 64)С + (С2 +1)611 + bsC2 + 61 - Ъз + 67 - [/1(-1,С)],
F15 = (2(Ъ7 + Ъ11) - Ъз)С - (С2 + 1)Ъ12 - Ъ10С2 - Ъ2 + Ъ5 - Ъ9 + f^-1 -С),
F16 = -[(2(67 + Ъп) - Ъз)с + (С2 +1)612 + ЪюС2 + 62 - 65 + 69 - [^2( 1,С)]•
Список литературы
[1] Plyasunova A. V., Sleptsov A. G. Collocation — grid method for solving nonlinear parabolic equations^ Rus. J. Theoret. and Appl. Mech., 1, No 1, 1991, 15—26^
[2] СлЕпцов А. Г., Шокин Ю.И. Адаптивный проекционно-сеточный метод для эллиптических задач^ Журн. вычислит. математики и математ. физики, 37, №5, 1997, 572-586-
[3] Семин Л. Г., Слепцов А. Г., Шапеев В. П. Метод коллокаций — наименьших квадратов для уравнений Стокса^ Вычислит. технологии, 1, №2, 1996, 90-98^
[4] СьярлЕ Ф. Метод конечных элементов для эллиптических задач. М^ 1980^
[5] Дьяконов Е. Г. Минимизация вычислительной работы. Асимптотически оптимальные алгоритмы для эллиптических задач. Наука, М., 1989.
[6] ЛЕБЕДЕВ В. И., Агошков В. И. Операторы Пуанкаре — Стеклова и их приложения в анализе. Изд-во ОВМ АН СССР, М^ 1983^
[7] Слепцов А. Г. Об ускорении сходимости линейных итераций IL Моделирование в механике• Новосибирск, 3 (20), №5, 1989, 118-125^
[8] Hearn A. C. REDUCE. User’s ManuaL Version 3^3^ RAND Publ., Santa Monica, 1987^
[9] BURGGRAF O. R. Analytical and numerical studies of the structure of steady separated flows^ J. Fluid Mech., 24, pt 1, 1966, 113-151
[10] Pan F., AcrivOS A. Steady flows in rectangular cavities^ Ibid., 28, pt 4, 1967, 643-655^
[11] BOZEMAN J.D., Dalton C. Numerical study of viscous flow in a cavity^ J. Comput. Phys., 12, No^ 3, 1973, 348-363^
[12] BRUNEAU C.H., JOURON C. An efficient scheme for solving steady incompressible Navier — Stokes equations^ Ibid., 89, No^ 2, 1990, 389-413^
[13] Ghia U., Ghia K.N., Shin C.T. High-Re solutions for incompressible flow using the Navier — Stokes equations and a multigrid method^ Ibid., 48, No^ 4, 1982, 387-411 •
[14] Deng G.B., Piquet J., Queutey p., Visonneau M. A new fully coupled solution of the Navier — Stokes equations^ Int. J. Numerical Methods in Fluids, 19, No^ 7, 1994, 605-639^
[15] Chen C.J., Chen H.J. Finite analytic numerical method for unsteady two-dimensional Navier— Stokes equations^ J. Comput. Phys., 53, No^ 2, 1984, 209-226^
Поступила в редакцию 5 января 1998 г.
