Метод коллокаций - наименьших квадратов для уравнений Стокса Текст научной статьи по специальности «Математика»

Вычислительные технологии

Том 1, № 2, 1996

МЕТОД КОЛЛОКАЦИЙ - НАИМЕНЬШИХ КВАДРАТОВ ДЛЯ УРАВНЕНИЙ СТОКСА*'*'

Л. Г. СЁмин

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

Новосибирск, Россия

A. Г. Слепцов

Институт вычислительных технологий СО РАН Новосибирск, Россия

B. П. Шапёёв

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

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

1. Введение

Широкое распространение метода конечных элементов объясняется простотой его физической интерпретации, математической формы и гибкостью численного алгоритма. При конечно-элементной аппроксимации краевой задачи для системы уравнений Стокса основная трудность состоит в правильном учете уравнения неразрывности ^у V = 0. Один из подходов состоит в использовании пространств конечных элементов, в которых предполагается, что условие ^у ^ = 0 выполняется точно. Однако этот процесс часто приводит к усложнению элементов [1]. При другом подходе уравнение неразрывности выполняется приближенно.

*© Л. Г. Сёмин, А. Г. Слепцов, В. П. Шапеев, 1996.

^ Работа выполнена при финансовой поддержке Российского фонда фундаментальных исследований, гранты №96-01-01888, №95-01-01357а

Одним из вариантов метода конечных элементов является сеточный метод коллока-ций — наименьших квадратов. Впервые примененный для решения обыкновенных дифференциальных уравнений, в дальнейшем он использовался при решении уравнений в частных производных эллиптического и параболического типа [2-6], он хорошо зарекомендовал себя при решении эллиптических задач с малым параметром при старшей производной, где могут возникать пограничные и/или внутренние слои. В этом случае он использовался на адаптивных сетках [5]. Однако для систем уравнений метод не применялся.

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

Наш опыт решения эллиптических задач на анизотропных сетках [5] показал, что при использовании лагранжевых интерполянтов в каждой ячейке сетки и при применении классических согласованных конечных элементов получаются плохо обусловленные системы линейных алгебраических уравнений. В результате этого сетки сгущались в окрестности особенностей, но увеличения точности при этом не происходило. Поэтому для аппроксимации были применены кусочно-полиномиальные функции, согласованные в слабом смысле (см. [5]). Использование таких аппроксимаций для решения задачи Стокса затрудняет применение хорошо известных итерационных методов, разработанных для этих задач [7]. Поэтому мы и в данной работе применяем метод итераций по подобластям, с которым можно познакомиться, например, в [8]. Для ускорения сходимости линейных итераций применялся метод, основанный на проекции погрешности на подпространство невязок [9, 10].

2. Постановка задачи

Рассматривается краевая задача для уравнений Стокса [11, 12]:

Я-1 Avj — рх. = ¡3, (хъх2) е П, 3 = 1, 2,

&у V =0, (1)

VI ап = V.

Приближенное решение отыскивается в виде кусочно-полиномиальной функции на регулярной сетке. Вводятся следующие обозначения: к — половина ширины ячейки, (х1т,х2т) — координаты центра т-й ячейки, у1 = (х1 — х1т)/к, у2 = (х2 — х2т)/к —

локальные координаты в ячейке, И(у1,у2) = V(x1,x2), д(у1,у2) = кр(х1,х2). Тогда уравне-

ния (1) в локальных переменных примут вид

Ащ — Яду. = Як21з, 3 = 1, 2, (2)

&у и = 0, (3)

и|дп = и. (4)

Пусть известно некоторое приближенное решение (а1,112,д). Искомое уточненное решение представим в виде: и1 = и1 + и1, и2 = и2 + и2, д = д + д. Подставляя эти выражения в уравнения (2), получим:

Ли, - Еду. = , з = 1, 2, (5)

где ^ = Е{к21, + ду.) - Ли,.

В каждой ячейке приближенное решение отыскивается в виде

\д/ = азтфз , (6)

' ' 3

где — базисные функции (векторы с тремя компонентами), т — номер ячейки.

Коэффициенты а,т будут определяться из уравнений коллокаций и условий согласования на границах соседних ячеек или краевых условий на дП.

3. Приближенные уравнения

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

k 1 2 3 4 5 6 7 8 9 10 11 12

1 0 Vi V2 0 0 У2! V2 0 -2ViV2 V2 + v2 -2ViV2

ук 0 1 -У 2 0 Vi 0 -2ViV2 0 V2 v2 -2ViV2 V2 + v2

0 0 0 0 0 1 i Vi 1 CN 2R-iVi 2 iV i- R 2R-iV2 0 0

Для получения единственного кусочно-полиномиального решения необходимо на границах между ячейками задать условия согласования. В качестве последних рассматриваются условия непрерывности выражений:

д]]<‘ + и ди, + и

ЦП - р + пи"’ ап + (7)

ТП - p - Wn’ тп- nUt■ (8)

Здесь Un и Ut — нормальная и касательная к границе ячейки компоненты скорости, n — вектор единичной внешней к ячейке нормали, а ц — некоторый положительный параметр, с помощью которого можно управлять обусловленностью получающейся системы линейных алгебраических уравнений и влиять на скорость сходимости итераций.

Для определения коэффициентов ajm в представлении (6) был использован метод итераций по подобластям, в котором подобластями являлись каждая из ячеек сетки. При уточнении решения в m-й ячейке применялись условия (7), которые принимают вид (8) при переходе к соседним ячейкам. В общем случае эти условия не могут быть выполнены в каждой точке границы между ячейками. Поэтому мы требуем, чтобы они выполнялись в смысле наименьших квадратов в восьми точках границы ячейки.

Если граница ячейки лежит на границе области, то в двух точках этой границы задается вектор скорости. Кроме того, в левом нижнем углу области задается давление.

К условиям согласования (краевым условиям) добавляются условия коллокации уравнений (2) в четырех точках внутри ячейки. В результате получается следующая система линейных алгебраических уравнений:

12

Bik a km = Fi, l = 1, , 24. (9)

k=1

В этой системе записаны уравнения коллокации условий согласования при l = 1, ... , 4 — на нижней, при l = 5, ... , 8 — на правой, при l = 9, ... , 12 — на верхней и при l = 3, ..., 16 — на левой границах. l = 17, ... , 24 соответствуют условиям коллокации дифференциальных уравнений (2). Условия согласования в локальных координатах задаются в точках с координатами (±1, ±Z), (±С, ±1), краевые условия — в точках (±1, ±£), (±£, ±1), условия коллокации дифференциальных уравнений (2) — в точках (±ш, ±ш). Величины Z, £, & положительны, не превышают единицы, их можно выбирать разными способами для получения хорошо обусловленной матрицы системы линейных алгебраических уравнений.

Здесь ввиду громоздкости выражений приведем выборочно только часть строк матрицы B и вектора F.

Для случая условий согласования на границах ячейки в принятых обозначениях элементы матрицы B имеют следующий вид:

Если на границах ячейки задавались условия согласования, то в этих обозначениях элементы матрицы B имеют следующий вид

l = 1: Bn = -^, Bi2 = Bi5 = Bi6 = Bn = 0, Віз = ^Z, Bi4 = 1, Bl7 = -Z2,

Bis = - n+2, Biio = Bii2 = 2Z, Bin = -(Z2++)n+2.

l = 5: Bii = Bi4 = Bi6 = Bis = ° Bi2 = , Bi3 = Z, Bi5 = 1, Bi7 = Bi11 = 2Z, Bi9 = П+2,

B = JÜL- Z 2 B = (Z 2+1)n+2 i10 n+1Z , i12 n+1 .

l = 9: Bi1 = n+-1, Bi2 = Bi5 = Bi6 = Bi9 = 0, Bi3 = -^Z, Bu = 1, B„ = n+lZ2, Bl8 = П+-,

Bi10 = Bi12 = 2Z, Bin = (z2++n+2.

l = 13: Bi1 = Bi4 = Bi6 = Bi8 = ° Bi2 = - n+1, Bi3 = - n+1Z, Bi5 = 1, Bi7 = Bi11 = 2Z,

B n+2 B n Z 2 b (Z,2 + 1)rl+2

Bi9 = - n+1, Bi10 = - n+1Z , Bi12 = n+1 .

Соответствующие компоненты вектора F выглядят следующим образом:

F1 = {nZ2(bj + bn - ay - a^) + nZ№ - Ьз - 2(bw + bu + a1¡0 + )] + n(h - b4 + bs + bn -

al - a4 - a8 - al1) + 2Z (aW + al2 - b10 - b12) - b4 + a4 + 2(bS + b11 + a8 + an)} / (П + ^

F5 = -{nZ2(b10 + b12 - a1o - a12) + nZ[Ьз - a? + 2(b7 + bn + a? + ah)] + n(b2 + b5 + b9 + bu -a% + a? - a? - a12) + 2Z(b7 + bn - a? - a^) + b5 - a? + 2(b9 + bu + a? + a12)}/(n + 1),

F9 = -{nZ 2 (b7 + b11 - a7 - a11) + nZ [a3 - b3 + 2(b10 + b12 + a10 + a<!2)] + n(b1 + b4 + bS + b11 -

al + a4 - as - al1) + 2Z (b10 + b12 - a\o - a12) + b4 - a4 + 2(bS + b11 + as + a11)}/(n + ^

F13 = {nZ 2(b10 + b12 - aio - a?2) + nZ [Ьз - - 2(b7 + bn + aЛ + a^)} + n(b2 - b5 + bg + bu -

Здесь аН, о,П, а®, аЛ — коэффициенты аг из нижней, правой, верхней и левой ячеек соответственно, Ьг — коэффициенты аг из данной ячейки, взятые с предыдущей итерации.

Если же на границах ячейки задавались краевые условия, то соответствующие строки матрицы В записываются иначе:

1 — 1: В11 — В18 — 1, В12 — В15 — В16 — В19 — О, В13 — —£, В14 — -1, В17 — £2,

В110 — В112 — —2C, ВИ1 — С2 + 1.

I — 5: В11 — В13 — В17 — 1, В12 — В15 — В16 — В19 — 0, В14 — —^ В18 — C2, В110 — В112 — 2C,

в И1 — С2 + 1.

I — 9: Вц — В14 — В18 — ~1, В12 — В15 — В16 — В19 — 0, В13 — — С, В17 — В110 — В112 — 2С,

Вт — С2 + 1.

I — 13: Вц — Ви — 1, В12 — В15 — В16 — В19 — О, В13 — -1, Вц — -С, В18 — С2,

Вао — ВИ2 — —2С, Вт — С2 + 1

Система (9) переопределена, в ней 24 уравнения и 12 неизвестных. Для того, чтобы сформулировать, что понимается под решением этой системы, рассмотрим два функционала:

Первый функционал соответствует сумме квадратов невязок уравнений, полученных из условий согласования или краевых условий, второй — сумме квадратов невязок уравнений коллокаций. Решение системы (9) определяется из условия минимума этих функционалов, причем минимум 1 берется по коэффициентам От, ] — 1, ..., 10, при фиксированных а11т и а12т, а минимум 2 — по а3-т, ^ — 11, 12, при фиксированных а3-т, ] — 1, ... , 10.

Таким образом, для определения коэффициентов а3-т, ] — 1, ... , 12, получаем систему уравнений вида

3 = 1

решение которой находим прямым методом исключения.

Для определения решения задачи применяется итерационный метод, в котором решение уточняется для каждой ячейки области в отдельности. На каждой итерации при

— -[{2(Ь10 + Ь12) — Ь4)С + (С,2 + 1)Ь11 + Ь8С 2 + Ь1 + Ь3 + Ь7 — и1(1, —

Р9 — (Ь3 — 2(Ь10 + Ь12))С — (С,2 + 1)Ь11 — Ь7С,‘2 — Ь1 — Ь4 — Ь8 + и1 (-С, ^

Р13 — (2(Ь10 + Ь12) + Ь4)С — (С,2 + 1)Ь11 — Ь8С 2 — Ь1 + Ь3 — Ь7 + и1 ( —1, —С).

16 12 24 12

1 — ^(^ В3 азт — Р1 )2, 2 — ^(^ В3 а3т — Р^2.

(10)

1=1 3=1

1=17 3=1

12

(11)

уточнении решения в m-й ячейке решение в соседних ячейках либо уже уточнялось, либо берется с предыдущей итерации. Очевидно, что здесь можно достичь высокой степени распараллеливания вычислений.

Различались несколько типов ячеек: все грани ячейки находятся внутри области (“внутренние” ячейки); одна грань ячейки лежит на границе области; две соседние грани ячейки лежат на границе области.

Аналитические выражения для элементов матрицы D и вектора правых частей F системы (11) для всех указанных типов ячеек находились программой, написанной в системе REDUCE [13]. Исходные уравнения, базисные функции, вид решения записаны в ней в качестве входных данных. Точки постановки условий согласования, граничных условий, уравнений коллокаций, параметр ц входят в полученные формулы в символьном виде, что позволяет при организации численного счета рассматривать различные варианты, не изменяя формул.

4. Тестирование

Правильность полученных формул, сходимость и работоспособность метода проверялись на ряде тестовых задач.

1. Область квадратная, на границах скорость равна нулю, давление в левом нижнем углу области равно единице, внешние силы отсутствуют. Точное решение: v1 = v2 = 0, p = 1.

2. Течение между двумя параллельными пластинами. Разность давлений на входе и выходе равна нулю, верхняя пластина движется, внешние силы отсутствуют. Точное решение: линейный профиль скорости, давление постоянно.

3. Течение между двумя параллельными покоящимися пластинами. Давление на входе и выходе разное, внешние силы отсутствуют. Точное решение: параболический профиль скорости, линейный закон распределения давления.

4. Комбинация двух предыдущих тестов. Течение между двумя параллельными пластинами. Верхняя пластина движется, давление на входе и выходе разное, внешние силы отсутствуют.

В этих тестах в качестве начального приближения выбиралось возмущенное точное решение. Расчеты проводились на сетках 10 х 3, 20 х 6, 30 х 9. Во всех тестах 1 - 4 в процессе итераций численное решение сходилось к точному.

5. Область — единичный квадрат, внешние силы f1 = п sin(nxi) — R-1n2 sin(nx2), f2 = п sin(nx2) — R-1n2 sin(nxi). Краевые условия: на левой и правой границах (x1 = 0 и x1 = 1) v1 = sin(nx2), v2 = 0; на нижней и верхней границах (x2 = 0 и х2 = 1) v1 = 0, v2 = sin(nx1); p(0,0) = 1. Точное решение: v1 = sin(nx2), v2 = sin(nx1), p = (cos(nx1) + cos(nx2))/2. В качестве начального приближения брались нулевые значения скорости и давления.

№ n1 п2 Погрешность Величина шага сетки

1 4 4 0.03568730 0.25

2 8 8 0.00518029 0.125

3 16 16 0.00051893 0.0625

4 32 32 0.00006088 0.03125

На этом тесте исследовалась сходимость метода на последовательности сеток. Значения точек постановки условий согласования, граничных условий, уравнений коллокаций

Рис. 1. Рис. 2.

в локальных координатах £ — С — ш — 0.5. Результаты приведены в таблице. Здесь п1 и п2 — количество ячеек по осям Х\ и Х2 соответственно. Погрешность вычислялась как максимум модуля разности между точным и численным решением. В программе в каждой ячейке этот максимум вычислялся по 121 точке, равномерно распределенной по ячейке.

Видно, что погрешность при переходе от сетки 4 х 4 к сетке 8 х 8 уменьшается примерно в 6.89 раза, при переходе от сетки 8 х 8 к сетке 16 х 16 — в 9.98 раза, при переходе от сетки 16 х 16 к сетке 32 х 32 — в 8.52 раза. Можно сделать вывод, что порядок сходимости метода на гладких решениях не хуже второго.

Рис. 3.

5. Расчет течения в прямоугольной каверне с движущейся верхней границей

С помощью коллокационно-сеточного метода проведен расчет течения в прямоугольной каверне с движущейся верхней границей. Высота каверны принималась равной единице. Расчеты проводились для небольших чисел Рейнольдса. Обозначим через A отношение высоты каверны к ее длине. На рис. 1-3 приведены картины линий тока в случаях A = 0.5, 1 и 2, соответственно. Видно, что при A = 0.5 и 1 в каверне существует один вихрь, а при A = 2 происходит образование второго (медленного) вихря. Этот результат находится в полном соответствии с данными физических экспериментов, полученными Паном и Акривосом [14], и расчетами, проведенными Бургграфом [15], Паном и Акриво-сом [14] методом релаксации для уравнений Навье—Стокса при малых числах Рейнольдса. В расчетах при A = 2 нижний вихрь по интенсивности на два порядка слабее верхнего. Тем не менее рассмотренный нами метод хорошо его передает. Это также косвенно подтверждает, что предложенный численный метод имеет высокий порядок аппроксимации.

Авторы благодарят С. В. Мелешко за обсуждение работы.

Список литературы

[1] СьярлЕ Ф. Метод конечных элементов для эллиптических задач. М., 1980.

[2] Слепцов А. Г. Коллокационно-сеточное решение эллиптических краевых задач. Моделирование в механике, 5(22), №2, 1991, 101-126.

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

[4] PlyasüNüva A. V., Sleptsüv A. G. Collocation-grid method for solving nonlinear

parabolic equations. Rus. J. of Theoretical and Applied Mechanics, 1, №1, 1991, 15-26.

[5] СЛЕПЦОВ А. Г., Шокин Ю. И. Адаптивный проекционно-сеточный метод для эллиптических задач. Журн. вычисл. матем. и матем. физ., 37, №5, М., 1997, 164-167.

[6] Sleptsüv A. G. Grid-projection solution of elliptic problem for a irregular grid. Russ. J. Numer. Analys. and Math. Modelling, 8, №6, 1993, 501-525.

[7] Дьяконов Е. Г. Минимизация вычислительной работы. Асимптотически оптимальные алгоритмы для эллиптических задач. Наука, М., 1989.

[8] Лебедев В. И., Агошков В. И. Операторы Пуанкаре-Стеклова и их приложения в анализе. Изд-во ОВМ АН СССР, М., 1983.

[9] СЛЕПЦов А. Г. Об ускорении сходимости линейных итераций. Моделирование в механике, 3(20), ч. 1, №3, 132-147, ч. 2, №5, 1989, 118-125.

[10] Karamyshev V. B., Küvenya V. M., Sleptsüv A. G., Cherny S. G. Variational

method of accelerating linear iterations and its applications. Computers Fluids, 25, №5, 1996, 467-484.

[11] Кочин Н. Е., КИБЕЛЪ И. А., Розе Н. В. Теоретическая гидромеханика. Ч. 2., М., 1963,

[12] ТЕМАМ Р. Уравнения Навье-Стокса. Теория и численный анализ. Мир, М., 1981.

[13] HEARN A. C. REDUCE-3 User’s Manual, version 3.3. Rand Corporation Publication CP78 (7/78), 1987.

[14] PAN F., AcrivOS A. Steady flows in rectangular cavities. J. Fluid Mechanics, 28, pt. 4, 1967, 643-655.

[15] BüRGGRAF O. R. Analytical and numerical studies of the structure of steady separated flows. J. Fluid Mechanics, 24, pt. 1, 1966, 113-151.

Поступила в редакцию 23 декабря 1996 г.