CC BY
20
УДК 519.63:532.51 DOI: 10.25513/2222-8772.2017.3.70-85
ПРИМЕНЕНИЕ ИНТЕГРАЛЬНОЙ ФОРМЫ УРАВНЕНИЙ КОЛЛОКАЦИЙ И ДИФФЕРЕНЦИАЛЬНЫХ УСЛОВИЙ СОГЛАСОВАНИЯ В МЕТОДЕ КОЛЛОКАЦИЙ И НАИМЕНЬШИХ КВАДРАТОВ
'Институт теоретической и прикладной механики СО РАН им. С.А. Христиановича 2Новосибирский национальный исследовательский университет
Аннотация. С целью повышения точности расчётов по методу коллокаций и наименьших квадратов (КНК) предлагается увеличить число степеней свободы с помощью следующих двух способов: увеличения числа базисных функций и интегрирования линеаризованных уравнений с частными производными (УЧП) по подъячейкам каждой ячейки пространственной расчётной сетки. Показано, что предложенные новые варианты метода КНК обладают более высокой точностью, чем предыдущие версии этого метода. Кроме того, вариант метода КНК, который использует интегральную форму уравнений коллокаций, требует для сходимости меньшего числа итераций, чем «дифференциальный» метод КНК.
Ключевые слова: метод коллокаций и наименьших квадратов, предобу-славливание, подпространства Крылова, многосеточные алгоритмы, колло-кация интегральных соотношений, уравнения Навье-Стокса.
Введение
В настоящее время широко применяется численное моделирование различных процессов в технике и промышленности с помощью численного решения начально-краевых задач для нелинейных уравнений с частными производными (УЧП). В частности, некоторые прикладные задачи, в которых необходимо решать уравнения Навье-Стокса, являются настолько вычислительно трудоёмкими, что требуют для своего решения на суперкомпьютерах от нескольких недель до одного года [1,2]. В этой связи является актуальной, как и прежде, разработка более эффективных методов численного решения систем УЧП, которые позволили бы заметно снижать требуемое машинное время.
Метод КНК, предложенный в [3], получил дальнейшее развитие в последующих работах других авторов, и в настоящее время он является одним
Е.В. Ворожцов1
профессор, д.ф.-м.н., e-mail: [email protected] В.П. Шапеев12
профессор, д.ф.-м.н., e-mail: [email protected]
из методов, позволяющих эффективно решать системы уравнений с частными производными, возникающие в разнообразных приложениях [4-13]. Были построены варианты метода, которые позволили получать численные решения двумерных и трёхмерных эталонных задач гидродинамики, по своей точности не уступающие наиболее точным решениям, известным в настоящее время и полученным другими методами [14,15].
В течение последних трёх десятилетий за рубежом получает всё большее распространение класс численных методов LSFEM (Least-Squares Finite Element Method) — «метод конечных элементов и наименьших квадратов» [16-19]. В нём метод конечных элементов (МКЭ) комбинируется с методом наименьших квадратов. В LSFEM подлежащие решению уравнения с частными производными сначала интегрируются в МКЭ по каждому конечному элементу, который представляет собой подобласть пространственной расчётной области. По аналогии с LSFEM здесь рассмотрим следующие варианты метода КНК: (i) вариант, в котором все уравнения коллокаций, получаемые из системы УЧП, заменяются на их интегральные аналоги, получаемые при интегрировании УЧП по нескольким подъячейкам, на которые разбивается каждая ячейка расчётной области; (ii) вариант, в котором используются и уравнения коллокаций, получаемые из решаемых УЧП, и уравнения, получаемые путём интегрирования УЧП по подъячейкам. Рассматриваемые варианты метода КНК предложены впервые.
Ниже описываются и исходный «дифференциальный» вариант метода КНК, и методы КНК, в которых «дифференциальный» вариант комбинируется с интегральной формой уравнений коллокаций и дифференциальной формой условий согласования. Представлены примеры расчётов, которые показывают, что новые модификации метода КНК позволяют получать более точные результаты, чем в случае «дифференциальных» вариантов метода КНК.
1. «Дифференциальный» метод КНК
1.1. Описание метода
Рассмотрим краевую задачу для системы уравнений Навье-Стокса
в области П с границей дП. В уравнениях (1) х^х2 — декартовы пространственные координаты; V = (^(х!,х2)^2(хьх2)) — вектор скорости; р = р(хьх2) — давление; f = (/1,/2) — заданная вектор-функция; Re — число Рейнольдса, А = щг + щ, (V ■ V) = V!^ + ^2-¿2. Система (1) решается при граничных условиях Дирихле (2), где g = g(x1,x2) = (#ьд2) — заданная вектор-функция. Давление определяется из (1), (2) с точностью до постоянной. В дальнейшем будем подбирать эту постоянную так, чтобы выполнялось условие
(V -V)V + Vp = —AV - f, div V = 0, (xi,x2) e Q,
1
(1) (2)
Jfn pdx1 dx2 = 0.
(3)
В качестве области решения берётся квадрат
П = |(Ж1,Ж2), 0 ^ Хг ^ Ь,1 = 1, 2}, (4)
где Ь > 0 — заданная длина стороны квадрата. Величина Ь использовалась в конкретных расчётах в качестве характерной длины при обезразмеривании переменных, и она входит естественным образом в определение числа Рейнольдса Re в (1). Далее краевую задачу для УЧП будем называть дифференциальной задачей.
В данной задаче (1)-(4) область (4) покрывается сеткой из квадратных ячеек Пу, г,] = 1,... ,1, I ^ 1. Удобно ввести локальные координаты у1, у2 в каждой ячейке Пу. Зависимость локальных координат от глобальных координат х1, х2 задаётся формулами ут = (хт — хт^)/к, т = 1,2, где хт,гу — значение координаты хт в центре ячейки Пу, а к — половина длины стороны квадратной ячейки, к = Ь/(2М). Пусть и(у1,у2) = (п1,п2) = У(ку1 + ,ку2 + х2,гу), я(у1,у2) = р(ку1 + х1}гу,ку2 + х2гу). В локальных переменных уравнения Навье-Стокса принимают следующий вид:
Апт — Re к(п, ^ + «2^ + = Re ■ к21т, т = 1,2; (5) V Оу1 ду2 дут)
1( дп1 + ^ =0, (6)
к\ду1 ду2
где А = дуг + щ. Линеаризация по Ньютону уравнений (5) приводит к формуле
£[Апт — ■ к)КПт+ + п1+ пт,у1 + п2 Пту2 + П2+ Пту2 + )] = т, (7)
где т =1, 2, а в — номер итерации по нелинейности, в = 0,1, 2,..., п\,п2,д— известное приближение решения на в-ой итерации начиная с выбранного начального приближения с индексом в = 0, Ет = Re [к2/т — к (и\п8туу1 + п2пту2)] , пт,У1 = дпт/дуь, дут = дд/дут, 1,т =1,2. Здесь, как и в [12,13], введён задаваемый пользователем параметр £ с целью управления величиной числа обусловленности переопределённой системы линейных алгебраических уравнений (СЛАУ), которая должна решаться в каждой ячейке Пу.
Приближённое решение в каждой ячейке Пу ищется в виде линейной комбинации базисных вектор-функций
(«1,«2,д ^ )Т = £1=1 ЬЬ щ, (8)
где верхний индекс Т обозначает операцию транспонирования, а ть — задаваемое пользователем количество базисных вектор-функций. В рассматриваемом варианте метода ^ являются многочленами. Таким образом, искомое приближенное глобальное решение задачи является кусочно-полиномиальным. А кусок решения (8) в каждой ячейке представляет собой локальное решение в окрестности начала локальной системы координат в ячейке.
Таблица 1. Вид базисных функций tpi
1 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15
1 yi У2 у2 -2У1У2 y2 0 0 0 0 0 0 0 0 0
фг 0 -У2 0 -2У1У2 у2 0 1 У1 у2 0 0 0 0 0 0
0 0 0 0 0 0 0 0 0 1 У1 У2 у2 У1У2 У toto
В работе [12] для аппроксимации составляющих скорости использованы многочлены второй степени по переменным y1;y2, а для аппроксимации давления — многочлены первой степени, так что общее количество базисных вектор-функций в (8) составляло 12.
Ранее в [6] было показано, что можно повысить точность численного решения, получаемого по методу КНК, путём использования многочленов более высоких степеней. В связи с этим мы ниже применяем многочлен второй степени также для аппроксимации давления. В этом случае имеются в общей сложности 18 базисных функций в выбранном пространстве. Поскольку в уравнении неразрывности коэффициенты константы, то его легко удовлетворить за счёт выбора базисных полиномов фг. Нетрудно установить, что для этого требуется удовлетворить ими три линейных соотношения. В итоге из первоначальных 18 базисных полиномов независимыми останутся только 15. Они приведены в табл. 1. Их совокупность можно назвать соленоидальным базисом, так как divфг = 0. Таким образом, достоинством предложенного варианта метода является то, что уравнение неразрывности и, следовательно, закон сохранения массы удовлетворяется численным решением задачи во всей расчётной области. Отметим, что проблема обеспечения закона сохранения массы представляет собой существенную трудность для метода LSFEM [19].
Набор базисных функций, применявшийся в [12], получается из набора, представленного в табл. 1, если положить mb = 12 в (8), то есть если оставить в табл. 1 только первые 12 базисных вектор-функций.
Выпишем систему уравнений, определяющую приближённое решение в каждой ячейке. Назовём её далее локальной СЛАУ, а объединение всех локальных СЛАУ — глобальной СЛАУ. Основными уравнениями, определяющими решение приближённой задачи, являются уравнения коллокаций. Они получаются в результате подстановки в уравнения (7) выражений (8) и координат точек коллокаций. Количество этих точек и их расположение внутри ячейки может варьироваться в различных вариантах метода. В данной работе были реализованы три варианта задания координат точек коллокаций. Обозначим через Nc число точек коллокации внутри каждой ячейки. При Nc = 2 локальные координаты точек коллокаций таковы: (ш,ш), (—ш,ш), где ш — задаваемое пользователем значение в интервале 0 < ш < 1. При Nc = 4 локальные координаты точек коллокаций имеют вид (±ш, ±ш). В случае Nc = 8 координаты точек коллокаций задавались следующим образом: расположение первых четырёх точек было взято таким же, как при Nc = 4, а координаты следующих
четырёх точек задавались по формулам (±ш,0), (0,±ш).
Подставляя (8), а также численные значения координат каждой точки кол-локации в (7), получим 2ЫС линейных алгебраических уравнений:
Ет=14т■ е1 = к, V=1,...,2ыс. (9)
Следуя [12], дополним систему уравнений приближённой задачи в ячейке Пу линейными условиями согласования локального решения в каждой ячейке с локальными решениями, взятыми во всех соседних с ней ячейках. Эти условия записываются в отдельных точках (называемых точками согласования) на сторонах ячейки Пу, которые являются общими с соседними ячейками. Условия согласования берутся здесь в виде
к ^ + п(п+)п = к ^ + П(п-)п, (10)
к ^ + («+Г = к ^ + («-)т, (11)
д + = д-. (12)
Здесь к£ = к (п ^ + П2¿^ = Пщц+п^, п = (п ,п) — внешняя нормаль
д__I п д
Х1 2 дхг
к стороне ячейки Пу, (■)"", (^)т — нормальная и касательная составляющие вектора скорости по отношению к стороне ячейки, «+, п- — пределы функции п и д при стремлении аргументов к точке согласования изнутри и снаружи ячейки . Здесь, как и в [12], введён задаваемый пользователем параметр ц с целью управления величиной числа обусловленности матрицы СЛАУ, которая должна решаться в каждой ячейке Пу.
Для однозначного определения давления в решении задаём его значение в одной точке области либо аппроксимируем условие (3) по формуле
I (Д^ дду^) = Ь (—I* + Д^ д*ду1 ду2). (13)
Здесь I* — интеграл по всей области П, рассчитанный как сумма интегралов по каждой ячейке на предыдущей итерации, д* — давление в ячейке на предыдущей итерации.
Обозначим через Ыт число точек согласования для составляющих вектора скорости на сторонах каждой ячейки. При Ыт = 4 координаты этих точек согласования задаются формулами (±1,0), (0,±1). При Ит = 8 координаты точек согласования таковы: (±1, — (), (±1,(), (—С,±1), (С,±1), где 0 < ( < 1. В расчётах, результаты которых представлены ниже, использовалось значение ( = 1/2. Условия согласования для давления (12) задаются в четырёх точках с координатами (±1,0), (0,±1). Используя (8), подставим координаты этих точек в каждое из трёх условий согласования (10)-(12). Из первых двух условий получим 2Ыт линейных алгебраических уравнений для составляющих скорости. Подстановка представления (8) в (12) даёт ещё четыре линейных алгебраических уравнения согласования.
В настоящей работе давление задавалось в нижней левой вершине ячейки П1,1 или же использовалось условие (13). Если сторона ячейки лежит на границе области П, то в соответствующих точках на стороне ячейки вместо условий
согласования в локальной СЛАУ выписываются граничные условия: ит = дт, т = 1,2.
Объединяя уравнения коллокаций, согласования и уравнения, полученные из краевых условий, если ячейка П^- граничная, в каждой ячейке получим СЛАУ вида
А; ■ 1 = У^1, (14)
где 1 = ,..., Ь+т)Т. В вариантах метода КНК, описываемых в настоящей работе, система (14) является переопределённой.
Из всех возможных псевдорешений полученной переопределённой СЛАУ приближённой задачи здесь наиболее интересное, на котором достигается минимум эвклидовой нормы невязки её уравнений. Возникшую задачу наименьших квадратов линейной алгебры (терминология, использованная в [20]) решим сочетанием способа ^Я-разложения матриц локальных СЛАУ и итерационным методом решения СЛАУ приближенной задачи. Применение для этого способа ^Я-разложения матриц локальных СЛАУ в случае плохо обусловленных задач предпочтительнее, чем способа нормальных уравнений [9,13,20]. Для решения СЛАУ приближённой задачи применялся процесс, который условно можно назвать итерациями Гаусса-Зейделя. Одна глобальная (в + 1)-ая итерация заключалась в том, что в расчётной области П последовательно перебирались все ячейки. При этом в каждой ячейке локальное СЛАУ (14) решалось ортогональным методом (Гивенса или Хаусхолдера), в правой части уравнений (10), (12) в качестве и-, q- берутся либо значения решения на (в + 1)-ой итерации, если они уже сосчитаны на этой итерации, либо их значения на предыдущей итерации.
1.2. Предобуславливатели для метода КНК
В каждой ячейке П+ необходимо решать СЛАУ вида (14). Опустим в (14), для краткости, верхние и нижние индексы:
АХ = /. (15)
Число обусловленности прямоугольной матрицы А вычисляется по формуле
к(А) = А1 |М| А-1 ||, (16)
где предполагается, что квадратная матрица А1 = АтА несингулярная. В нашем случае предобуславливатель зависит от параметров В [12] был описан простой алгоритм для нахождения оптимальных значений £ор!, пр в любой ячейке из требования минимизации числа обусловленности к(£,п). Число к(£ор^пор0 обычно удовлетворяло неравенствам 3 < к(£ор!,пор0 < 10. Оказалось, что оптимальные значения £ор!, пор слабо зависят от положения конкретной ячейки в расчётной сетке, по меньшей мере, в случаях тех тестовых и эталонных задач, которые рассматривались в [12]. Некоторые свойства предобуславливателя были также исследованы в [12]. В частности, было показано, что снижение N сильнее влияет на значение £ор!, чем на значение пор^
1.3. Алгоритм ускорения сходимости на основе подпространств Крылова
Для ускорения сходимости итераций здесь применялся во всех новых вариантах метода КНК, описываемых в настоящей работе, новый вариант известного метода [21], основанный на подпространствах Крылова и подробно изложенный ранее в [9,13]. Ниже кратко описывается соответствующий алгоритм. Пусть СЛАУ имеет вид X = ТХ + /, где вектор X — искомое решение, Т — квадратная матрица, а / — вектор-столбец. Пусть матрица Т имеет полный ранг, и пусть сходится следующий итерационный процесс: Хга+1 = ТХга + /, п = 0, 1, ..., где Хп — приближение решения на п-ой итерации. По определению, гп = ТХп + / - X/п = Xп+1 - Xп — невязка уравнений X = ТХ + /, и нетрудно получить следующее соотношение из вышеприведённых формул: гп+1 = Тгп. Предположим, что сделано к + 1 итераций, начиная с начального приближения X0, то есть вычислены величины X 2,..., Xк+1 и г0, г1,..., г к . Значение Xк+1 затем уточняется по формуле X*к+1 = Xк+1 + Ук+1. Используется поправка вида
/к+1 = Е «г г * (17)
г = 1
с неопределёнными коэффициентами а1,..., ак, которые находятся из условия минимизации функционала невязки Ф(а1,...,ак) =|| X*к+1 — Т^?*к+1 — / ||2, возникающей при подстановке X*к+1 в систему X = ТХ + /. Здесь ||и||2 — евклидова норма вектора и размерности N. Уточнённый вектор к+1-го приближения X *к+1 используется в качестве начального приближения для дальнейшего продолжения последовательности итераций.
1.4. Ускорение сходимости итераций с помощью многосеточного алгоритма
Основная идея многосеточных алгоритмов состоит в селективном демпфировании гармоник ошибки [22,23]. В методе КНК, как и в других методах, количество итераций, необходимых для достижения заданной точности приближения к предельному решению, зависит от начального приближения. В работе применялась операция продолжения вдоль восходящей ветви У-цикла — расчёты на последовательности измельчающихся сеток — в качестве способа получения хорошего начального приближения для итераций на самой мелкой сетке среди сеток, используемых в многосеточном комплексе. Проиллюстрируем алгоритм операции продолжения на примере составляющей скорости и1(у1, у2,Ь1,..., Ь15). Пусть к = к, где к — полушаг грубой сетки, и пусть к2 = ^1/2 — полушаг мелкой сетки, на которой нужно найти разложение функции и1 по базису.
Шаг 1. Пусть X1 и X2 — глобальные координаты центра ячейки грубой сетки. Сделаем следующие подстановки в полиномиальное выражение для и1:
уг = (х — XI)/к1, I = 1, 2.
(18)
В результате получаем многочлен
U1(x1,x2,b1,..., Ь15) = щ (^, ^, bi,..., bi5>). (19)
Шаг 2. Пусть (XbX2) — глобальные координаты центра любой из четырёх ячеек мелкой сетки, содержащихся в ячейке грубой сетки. Сделаем в (19) замену xi = Xi + yi • h2, l = 1, 2. В результате получим многочлен второй степени U1 = P (y1,y2, b1,..., Ь15) от переменных y1,y2 с коэффициентами Ь1,..., Ь15. После приведения подобных оказывается, что координаты X1,X2 и X1,X2 входят в bi (l = 1,..., 15) только в виде комбинаций 5xi = (Xi — Xi)/h1. Согласно (18), величина —8xi = (Xi — Xi)/h1 является локальной координатой в ячейке грубой сетки центра ячейки мелкой сетки.
Приведём выражения для коэффициентов bj (j = 1,..., 15) представления решения в ячейке мелкой сетки с полушагом h2 в терминах коэффициентов ... ,Ь15 представления решения в ячейке с полушагом h1 = 2h2:
Ь1 = b1 — 5x1(b2 — b^5x 1) — 5x2(b3 + 2b56x1 — b&6x2),b2 = 01 (7\ + b56x2),
Ьз = 01[ЬЗ + 2(b56x1 — ^6x2)], b4 = 02b4, b5 = 02h, b6 = 02b6,
b7 = b7 — 6x1(b8 — bg6x1) + 6x2T1, b8 = o\_(b8 — 2bg6x1 + 2b46x2),
bg = 02bg, bw = bw — 6x{T2 — 6x2^12 — h56x2), Ьц = 0^2 — ^36x1),
b12 = О1Ф12 — bU6x1 — 2Ь 156x2), b13 = 02b13, b14 = 02b14, b15 = 02b15,
где 01 = h2/h1, 02 = 02, T1 = b2 — 2b46x1 + b56x2, T2 = bn — bVi6x1 — bu6x2. Заметим, что приведённые выше выражения для Ь]_,...,Ьд совпадают с выражениями, представленными в [12,13] для случая, когда mb =12 в (8).
2. Использование интегральной формы уравнений коллокаций
Выше в разделе 1 был описан «дифференциальный» вариант метода КНК, в котором уравнения коллокаций (7) были получены из дифференциальных уравнений (5). Здесь можно использовать вместо уравнений коллокаций (7) их интегральные аналоги, получаемые путём интегрирования уравнений (5) по нескольким подобластям. В «дифференциальном» варианте метода КНК областью влияния локального решения являются точки коллокации и точки согласования. В интегральном варианте областью влияния является множество всех точек ячейки. В этом варианте можно ожидать, что скорость сходимости итерационного процесса решения задачи будет больше, чем в случае «дифференциального» метода КНК.
Также отметим, что в [13] было показано, что включение аппроксимации (13) интеграла (3) в переопределённую СЛАУ (14) вместо задания давления в одной точке существенно ускоряет сходимость итерационного процесса.
Сначала введём в рассмотрение в каждой ячейке равномерную расчётную сетку, линии которой разбивают каждую сторону ячейки на Nsub отрезков,
где Х8иЬ > 1 — заданное количество отрезков вдоль каждой локальной координаты , к = 1,2. Линии этой сетки разбивают ячейку П? на N = Д,2иЬ подъячеек I = 1,...,ХС. Значение Х8иЬ задаётся так, чтобы величина Д2иЬ была сравнимой с количеством ть неизвестных коэффициентов в (8) или превосходящей ть. Затем представление (8) приближённого решения задачи подставляется в уравнения (7), которые интегрируются по каждой подъячейке
£ I I (¡) (Re ' h)(u1un+yi + u1+ Um,yi + U2Un+y2 + U2+ Um,y2
+ q^1)] dyidy2 = zjj Fndyidy2, m =1, 2; l = 1,... , Nc. (20)
Заметим, что интегрирование в левой части равенства (20) легко осуществляется в аналитическом виде, потому что подынтегральное выражение имеет полиномиальный вид в соответствии с табл. 1.
В вышеописанном варианте метода КНК интегрируются только линеаризованные дифференциальные уравнения с частными производными с целью получения уравнений коллокаций, но при этом условия согласования используются в дифференциальной форме, как и в «дифференциальном» методе КНК. Именно поэтому термин «интегральный», применяемый для краткости в следующем разделе при описании результатов расчётов по варианту метода КНК, описанному в данном разделе, берётся в кавычки.
3. Использование дифференциальной и интегральной форм уравнений коллокаций
Метод КНК является гибким численно-аналитическим методом решения краевых задач для уравнений с частными производными, и эта гибкость позволяет комбинировать в одной компьютерной программе и в одной переопределённой системе (14) уравнения коллокаций, получаемые непосредственно из линеаризованного дифференциального уравнения (7) и получаемые в результате интегрирования уравнения (7) по формуле (20).
В левой части (20) имеются, среди прочих, интегралы вида
Г Ут2
bij,fc / ymdym, m = 1,2, (21)
J Ут1
где ym1, ym2 — заданные пределы интегрирования, k = 2, 3,5,8,11,12,14 согласно табл. 1. Если |ym1| = |ym21, то интеграл (21) равен нулю, и в соответствующей строке матрицы Ajj системы (14) теряется информация об одном или нескольких коэффициентах bjj,fc. Следовательно, нужно подбирать подъячейки в ячейке Qj так, чтобы в каждой подъячейке выполнялось условие |ym1| = |ym2| хотя бы при одном m. Если оставить в матрице Ajj часть или все уравнения коллокаций, получаемые из УЧП (7), то можно добавить, с целью повышения
точности расчёта по методу КНК, несколько уравнений коллокаций вида (20). При этом можно брать в (20) особенно простые пределы интегрирования, например, ym1 = 0, ym2 = 1 или ym1 = -1, ym2 = 0. В этих случаях выражения для интегралов (20) существенно упрощаются, за счёт чего можно добиться некоторого снижения машинного времени, необходимого для решения конкретной задачи.
Пусть Nc,dif — количество уравнений коллокаций, полученных из дифференциальных уравнений (7), и Nc,int — количество уравнений коллокаций в интегральной форме (20). Можно по-разному задавать величины Nc,dif и Nc,int, нужно только следить за тем, чтобы величина Nc = Nc,dif + Nc,int была сравнимой с количеством неизвестных коэффициентов в (8) или превосходила его.
В следующем разделе приводится пример расчёта, когда в системе (14) использовались дифференциальные уравнения коллокаций и к ним ещё добавлялись четыре уравнения коллокаций в интегральной форме (20). При этом интегрирование в (20) при m =1 выполнялось по следующим двум подъячей-кам nj и nj ячейки Пу:
^ = {(У1,У2)|- 1 ^ У1 ^ 0, -1 ^ y2 ^ 1}, = { (yi, У 2) 10 ^ У1 ^ 1, -1 ^ У2 ^ 1}.
При m = 2 интегрирование выполнялось в (20) при другом разбиении ячейки Пу на две подъячейки:
П(3) = {(У1,У2)|-1 ^ У1 ^ 1,-1 ^ У2 ^ 0}, П^ = {(У1,У2)|-1 ^ У1 ^ 1, 0 ^ У2 ^ 1}. Очевидно, Пу = П(1) (J nj = nj U П^.
4. Результаты численных расчётов 4.1. Тестирование
Рассмотрим следующее точное решение уравнений Навье-Стокса (1) [24]:
-2(1 + ал) 2(1+ ал)
U1 = 7—-чо . /-, .-^, U2 ^
(1+ X.)2 + (1+ Ж2)2' 2 (1+ X.)2 + (1+ Х2)2' 2
Р = - (1+ X.)» + (1+ Х2)2 • 0 ^ * ^ « 1 (22)
Заметим, что функции и(х.,х2) и и2(х!,х2) описывают поле скорости с нулевой дивергенцией. Далее,
/ / pdx1dx2 = 4G - п ln2 - 2i '0 Jo
« -0.46261314677281549872
Li2f- 2 - Li2 2
где г = V-!, С - постоянная Каталана, С « 0.91596559417721901505, ЬЬ(г) -полилогарифмическая функция. Чтобы обеспечить выполнение равенства (3) с
погрешностью, не превышающей погрешность машинных вычислений, давление p в (3) заменялось на величину p = p + 0.4626131467728155.
Среднеквадратичные величины погрешности решения вычислялись по формулам:
MM 2 1
Err(u(h)) = [2m£EE(u™ - uVXij)2
i=1 j=1 V=1
1 MM 1
Err(p(h)) = - pex)2
i=1 j=1
где M — количество ячеек вдоль каждого координатного направления, Uj и p^j — вектор скорости и давление, вычисленные из точного решения (22). Величины Ujj и p^j обозначают численное решение, полученное по методу КНК, описанному выше, вычисленное в центре ячейки Qj. Порядки сходимости и vp вычисляются по известным формулам [8,10]. Пусть bj, s = 0,1,... — значения коэффициентов frj в (8) на s-ой итерации. Использовалось следующее условие окончания итераций: £bn+1 < е, где £bn+1 = max^-(max1^l^mb |— bj^l), а e < h2 — малая положительная величина. В дальнейшем будем называть величину #bn+1 псевдопогрешностью приближённого решения. Наряду с условием #bn+1 < е также применялся следующий критерий для окончания итераций по нелинейности:
¿un+1 =| un+1 — un ||<е2,
где || ■ || — евклидова норма вектора, е2 — заданная малая положительная величина.
В работе [13] введено определение степени переопределённости или недо-определённости системы (15) как величины х(А) = mr/mc, где mr и mc — соответственно, число строк и число неизвестных в СЛАУ с матрицей A. Было исследовано влияние величины х(А) на число обусловленности (16) в случае применения «дифференциального» метода КНК с mb =12 в (8). При этом использовались граничные условия Дирихле, соответствующие аналитическому решению (22).
Опустим для краткости полученные табличные данные и лишь перечислим выводы, полученные из проведённого исследования.
1°. Скорость сходимости итерационного процесса в «дифференциальном» методе КНК существенно зависит от чисел обусловленности алгебраических систем, которые необходимо решать в каждой ячейке.
2°. В тех случаях, когда в матрицу А включаются только строки, соответствующие уравнениям коллокаций, и одна строка, соответствующая аппроксимации (13) интеграла от давления, получается очень большое число обусловленности порядка 105, то есть матрица А является плохо обусловленной. При включении в матрицу A строк, соответствующих условиям согласования, происходит уменьшение числа обусловленности на три-пять десятичных порядков в зависимости от количества ячеек сетки, точек коллокаций и точек согласования.
Таблица 2. Погрешности Етт(ип), Егг(рп)
и их порядки сходимости на
последовательности сеток, Ие = 1000, Ь = 1, Мс = 8. Метод КНК12
Таблица 3. Погрешности Етт(ип), Етт(рп)
и их порядки сходимости vu, vp на последовательности сеток, Ие = 1000, Ь = 1, Мс = 8. Метод КНКД15
м Егг(и п) Егг(рп) V« ир
10 1.309е- 02 9.754е-03
20 5.484е- 03 3.664е-03 1.26 1.41
40 1.869е- 03 1.135е-03 1.55 1.69
80 4.938е- -04 2.870е-04 1.92 1.98
м Егг(и п) Егг(рп) V« Vp
10 2.947е- 05 4.731е-05
20 6.971е- 06 2.168е-05 2.08 1.13
40 2.289е- 06 1.046е-05 1.61 1.05
80 9.816е- 07 4.950е-06 1.21 1.09
Таблица 4. Погрешности Етт(ип), Етт(рп) и их порядки сходимости Vи, vp на последовательности сеток, Ие = 1000, Ь =1, Мс = 8. Метод КНКИ15
м Егг(и п) Егг(рп) V« Vp
10 2.910е- -05 5.437е-05
20 6.807е- 06 2.249е-05 2.10 1.27
40 2.325е- -06 1.076е-05 1.55 1.06
80 9.877е- -07 4.893е-06 1.24 1.14
Именно такое значительное снижение числа обусловленности обеспечивает работоспособность метода КНК при решении краевых задач для уравнений с частными производными.
3°. Число обусловленности в граничных ячейках всегда выше, чем в случае внутренних ячеек.
4°. При N = 4 и = 1 существенно замедляется скорость сходимости
метода КНК по сравнению со случаем N = 8, Жта1 = 2.
5°. Полученные данные являются практическим доказательством того, что при использовании переопределённой системы в приближённой задаче (х(А) > 1) соответствующая СЛАУ приближённой задачи получается лучше обусловленной, чем когда СЛАУ при этом не является переопределённой. Это важное преимущество метода коллокации и наименьших квадратов перед методом кол-локации, в котором не используется переопределённое СЛАУ в приближенной задаче.
Для изучения свойств сходимости и точности представленных выше вариантов метода КНК были проведены многочисленные вычислительные эксперименты с использованием аналитического решения (22). Были даны следующие названия различным вариантам метода КНК: КНКД12 — «дифференциальный» метод КНК [12,13], использующий 12 базисных векторов, то есть ть =12 в (8); КНКД15 — «дифференциальный» метод КНК с ть =15 в (8); КНКДИ15,Мс — метод КНК с ть = 15, в котором используются N уравнений коллокаций, полученных из дифференциальных уравнений (7), и четыре уравнения коллокаций в интегральной форме (20); наконец, КНКИ15 — «интегральный» метод КНК с ть =15 в (8).
1.0 -х2
0.8
0.5
> е -о-®- в -о-в- в -о-©- в -о©- -
0.6
0.0 „.ът-Ъ-'Т-Ч
^ Т.4
0.4
0.6
0.8
1.0
-0.5
0.2
0 0.5 1.
Рис. 1. Сравнение профилей приближенного и точного решений при х2 = Ь/2. Сетка 40 х 40 ячеек
Рис. 2. Решение эталонной задачи методом КНКДИ15,4 при Ие = 100 на сетке из 40 х 40 ячеек: профиль составляющей скорости вдоль центральной линии каверны х1 = 0.5 (сплошная линия); (о о о) ОЫа и др. [25]
В табл. 2-4 представлены результаты численных экспериментов, в которых использовались два из описанных выше способов ускорения сходимости: двухпараметрический предобуславливатель и метод подпространств Крылова. В этих расчётах были взяты число Рейнольдса Ие = 1000, Ь =1 в (4) и N = 8 в случае методов КНКД12 и КНКД15; в случае метода КНКИ15 использовали значение N = Д,2иЬ = 16. Критерием окончания расчётов по вариантам метода КНК было выполнение неравенства < 10-10. Из сравнения таблиц 2 и 3 видно, что увеличение количества базисных вектор-функций с 12 до 15 уменьшает ошибки в давлении и скорости в численном решении на два-три десятичных порядка.
На рис. 1 приведены профили компонент приближённого решения, полученного методом КНКД15, и точного решения. Компоненты ц и р приближенного решения нарисованы символами Д, о и V, соответствующие компоненты точного решения — сплошными, штриховыми и штрих-пунктирными линиями. Здесь видно хорошее согласие между численными результатами и аналитическим решением.
4.2. Течение в квадратной каверне с движущейся крышкой
В рассматриваемой задаче расчётная область — каверна — квадрат (4) со стороной Ь = 1, начало координат находится в её левом нижнем углу. Верхняя крышка каверны движется в безразмерных величинах с единичной скоростью в положительном направлении оси Ох1. Остальные стороны каверны (4) покоятся. На всех сторонах заданы условия прилипания: = 1, = 0 при х2 = 1 и = 0, т =1,2 на остальных сторонах.
Течение в каверне с движущейся крышкой имеет сингулярности в верхних углах области. Их влияние на точность численного решения усиливается с увеличением числа Рейнольдса. Поэтому при больших числах Рейнольдса для получения более точного решения необходимо применять приём исключения особенностей решения и/или адаптивные сетки с более мелкими ячейками в окрестности сингулярностей. Здесь использовали только равномерные сетки.
Было проведено сравнение точности различных вариантов метода КНК, опи-
Таблица 5. Погрешности 6п1, полученные при использовании различных вариантов метода
КНК, Ие = 100
Метод £ V K 11 mgr k CPU time, s.
КНКД12 2.0 3.5 4 8 2336 53.59 1.726e-02
КНКД15 2.0 3.5 4 8 2540 99.59 1.150e-02
КНКИ15 0.25 1.75 4 9 1700 110.38 8.521e-03
КНКДИ15,4 2.0 3.5 4 8 3257 84.66 8.443e-03
КНКДИ15,8 2.0 3.5 4 8 2697 78.89 1.005e-02
санных выше, в случае числа Рейнольдса Re = 100. Для сравнения использовалось численное решение из работы [25]. Введём в рассмотрение погрешность iui = maxj |ui,Ghia(0.5,) — mi,clr(0.5,)|, где координаты x2j были взяты из [25]. Во всех расчётах, результаты которых приведены в табл. 5, использовалась комбинация всех трёх способов ускорения сходимости итераций, кратко описанных в разделе 1. Оптимальные значения параметров п, входящих в двухпараметрический предобуславливатель, были найдены с помощью алгоритма [12]. Nit — количество итераций, необходимых для выполнения неравенства < 10-9. Значение k — количество невязок, применявшихся в методе Крылова (см. уравнение (17)). В многосеточном комплексе применялись сетки из 5 х 5, 10 х 10, 20 х 20 и 40 х 40 ячеек; Kmgr — количество последовательно использовавшихся сеток в многосеточном комплексе. Из табл. 5 видно, что метод КНКДИ 15,4 обеспечивает наилучшую точность среди четырёх рассмотренных вариантов метода КНК (см. также рис. 2).
Заключение
Представлены новые варианты метода КНК. В численных экспериментах показано, что предложенные варианты метода КНК позволяют получать более точные численные решения, чем «дифференциальный» вариант КНКД12 метода КНК.
Значительное (на два-три десятичных порядка) увеличение точности результатов, получаемых по методу КНК при небольшом (на 25%) увеличении числа базисных вектор-функций указывает на целесообразность (и осуществимость) обобщения этого приёма на случаи применения вариантов метода КНК для численного решения краевых задач для трёхмерных уравнений Навье-Стокса, так как тогда можно будет использовать для достижения заданной точности сетки с меньшим числом ячеек, чем в случае базисов меньшего размера.
Следует также отметить, что представленные варианты метода КНК очень хорошо подходят для распараллеливания, так как вычисления в каждой ячейке пространственной сетки могут выполняться на каждой итерации независимо от других ячеек. Можно разбить расчётную область на любое количество подобластей, равное числу доступных процессоров.
Литература
1. Li K., Li Q. Three-dimensional gravity-jitter induced melt flow and solidification in magnetic fields // J. Thermophys. Heat Transfer. 2003. V. 17, N. 4. P. 498-508.
2. Bailly O., Buchou C., Floch A., Sainsaulieu L. Simulation of the intake and compression strokes of a motored 4-valve SI engine with a finite element code // Oil & Gas Science and Technol. 1999. V. 54. P. 161-168.
3. Плясунова А.В., Слепцов А.Г. Коллокационно-сеточный метод решения нелинейных параболических уравнений на подвижных сетках // Моделирование в механике. 1987. Т. 18, № 4. С. 116-137.
4. Семин Л.Г., Слепцов А.Г., Шапеев В.П. Метод коллокаций-наименьших квадратов для уравнений Стокса // Вычисл. технологии. 1996. Т. 1, № 2. С. 90-98.
5. Semin L., Shapeev V. Constructing the numerical method for Navier-Stokes equations using computer algebra system // Proc. 8th Internat. Workshop on Computer Algebra in Scientific Computing. LNCS. 2005. V. 3718. P. 367-378.
6. Исаев В.И., Шапеев В.П. Варианты метода коллокаций и наименьших квадратов повышенной точности для численного решения уравнений Навье-Стокса // Ж. вычисл. матем. и матем. физ. 2010. Т. 50, № 10. С. 1758-1770.
7. Shapeev V.P., Isaev V.I., Idimeshev S.V. The collocations and least squares method: application to numerical solution of the Navier-Stokes equations // CD-ROM Proc. 6th ECCOMAS, Sept. 2012, Vienna Univ. of Tech. 2012. ISBN: 978-3-9502481-9-7.
8. Shapeev V.P., Vorozhtsov E.V. Symbolic-numeric implementation of the method of collocations and least squares for 3D Navier-Stokes equations // Proc. 14th Internat. Workshop on Computer Algebra in Scientific Computing. LNCS. 2012. V. 7442. P. 321-333.
9. Шапеев В.П., Ворожцов Е.В., Исаев В.И., Идимешев С.В. Метод коллокаций и наименьших невязок для трёхмерных уравнений Навье-Стокса // Вычислит. методы и программирование. 2013. Т. 14. С. 306-322.
10. 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 // Proc. 15th Internat. Workshop on Computer Algebra in Scientific Computing. LNCS. 2013. V. 8136. P. 381-392.
11. Shapeev V. Collocation and least residuals method and its applications // EPJ Web of Conferences. 2016. V. 108, N. 01009. DOI: 10.1051/epjconf/201610801009.
12. Shapeev V.P., Vorozhtsov E.V. Symbolic-numeric optimization and realization of the method of collocations and least residuals for solving the Navier-Stokes equations // Proc. 18th Internat. Workshop on Computer Algebra in Scientific Computing. LNCS. 2016. V. 9890. P. 473-488.
13. Ворожцов Е.В., Шапеев В.П. О комбинировании способов ускорения сходимости итерационных процессов при численном решении уравнений Навье-Стокса // Вы-числ. методы и программирование. 2017. Т. 18. С. 80-102.
14. 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. V. 31. P. 1874-1900.
15. Botella O., Peyret R. Benchmark spectral results on the lid-driven cavity flow // Comput. Fluids. 1998. V. 27. P. 421-433.
16. Carey G.G., Cheung Y.K., Lau S.L. Mixed operator problems using least squares finite element collocation // Comput. Meth. Appl. Mech. Engng. 1980. V. 22. P. 121-130.
17. Jiang B.N. The Least-Squares Finite Element Method: Theory and Applications in Computational Fluid Dynamics and Electromagnetics. Berlin: Springer, 1998. 418 p.
18. Bochev P.B., Gunzburger M.D. Least-Squares Finite Element Methods. New York: Springer-Verlag, 2009.
19. Ranjan R., Chronopoulos A.T., Feng Y. Computational algorithms for solving spectral/hp stabilized incompressible flow problems // J. Math. Res. 2016. V. 8, N. 4. P. 21-39.
20. Деммель Дж. Вычислительная линейная алгебра. Теория и приложения. М. : Мир, 2001. 430 с.
21. Крылов А.Н. О численном решении уравнения, которым в технических вопросах определяются частоты малых колебаний материальных систем // Изв. АН СССР, Отд. матем. и естеств. наук. 1931. № 4. С. 491-539.
22. Федоренко Р.П. О скорости сходимости одного итерационного процесса // Ж. вы-числ. матем. и матем. физ. 1964. Т. 4, № 3. С. 559-564.
23. Wesseling P. An Introduction to Multigrid Methods.Chichester: John Wiley & Sons, 1992. 284 p.
24. 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. V. 227. P. 4018-4037.
25. 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. V. 48. P. 387-411.
APPLICATION OF THE INTEGRAL FORM OF COLLOCATION EQUATIONS AND DIFFERENTIAL MATCHING CONDITIONS IN THE METHOD OF COLLOCATIONS AND LEAST SQUARES
E.V. Vorozhtsov1
Leading Research Scientist, Dr.Sc. (Phis.-Math.), e-mail: [email protected]
V.P. Shapeev12 Professor, Dr.Sc. (Phis.-Math.), e-mail: [email protected]
1Khristianovich Institute of Theoretical and Applied mechanics SB RAS, Novosibirsk 2Faculty of mechanics and mathematics, Novosibirsk National Research University
Abstract. To increase the accuracy of computations by the method of collocations and least squares (CLS) it is proposed to increase the number of degrees of freedom with the aid of the following two techniques: an increase in the number of basis vectors and the integration of the linearized partial differential equations (PDEs) over the subcells of each cell of a spatial computational grid. It is shown that the proposed new versions of the CLS method possess a higher accuracy than the previous versions of this method. Furthermore, the version of the CLS method, which employs the integral form of collocation equations, needs a lesser number of iterations for its convergence than the "differential" CLS method.
Keywords: method of collocations and least squares, preconditioning, Krylov subspaces, multigrid, collocation of integral relations, Navier-Stokes equations.
Дата поступления в редакцию: 28.04.2017
