CC BY
46
Вычислительные технологии
Том 12, № 3, 2007
ИССЛЕДОВАНИЕ СВОЙСТВ МЕТОДА КОЛЛОКАЦИИ И НАИМЕНЬШИХ КВАДРАТОВ РЕШЕНИЯ КРАЕВЫХ ЗАДАЧ ДЛЯ УРАВНЕНИЯ ПУАССОНА И УРАВНЕНИЙ НАВЬЕ-СТОКСА*
В. И. Исаев, В. П. Шапеев, С. А. Еремин Институт теоретической и прикладной механики СО РАН,
Новосибирск, Россия e-mail: vadimx@gorodok.net, shapeev@itam.nsc.ru
A new version of the collocation and least squares (CLS) method is proposed and implemented. The dependence of method’s properties on its parameters (coordinates of collocation points, coordinates of matching conditions points on cells’ boundaries, weight parameters of equations) is investigated. Essential dependence of smoothness and accuracy of the approximate solution and of other properties of the proposed method on the parameters is observed. It is shown that the research of that kind is necessary.
Введение
В данной работе проведено исследование свойств метода КНК численного решения краевых задач. С целью выявления общих свойств метода рассматриваются две существенно различные задачи: первая из них — для уравнения Пуассона, вторая — для системы уравнений Навье—Стокса. В формулы метода КНК вводятся параметры: координаты точек коллокации внутри ячеек, координаты точек записи условий согласования численного решения на границе ячеек, весовые множители при отдельных слагаемых в условиях согласования и весовые множители при уравнениях в переопределенной системе. Условно будем называть их управляющими параметрами. При решении конкретных краевых задач необходимо задавать значения этих параметров. В случае, когда метод КНК применяется для решения обыкновенного дифференциального уравнения (ОДУ), можно указать выбор значений отдельных параметров при фиксированной системе базисных функций, дающий наименьшую погрешность численного решения [1]. При применении метода КНК для решения задач для уравнений с частными производными проблема выбора управляющих параметров ранее практически не исследовалась. Однако в численных экспериментах можно наблюдать, что за счет выбора значений параметров удается добиться существенного улучшения точности приближенного решения и обусловленности задач линейной алгебры, возникающих при использовании метода.
* Работа выполнена при финансовой поддержке Российского фонда фундаментальных исследований (грант № 06-01-00080-а).
© Институт вычислительных технологий Сибирского отделения Российской академии наук, 2007.
Реализация метода КНК сводится к решению переопределенной системы линейных алгебраических уравнений (СЛАУ). Отдельный интерес вызывает вопрос о том, насколько при этом должна быть переопределена СЛАУ. Здесь исследовано влияние количества уравнений в СЛАУ на свойства метода. Рассмотрены несколько способов реализации метода, в которых переопределенная СЛАУ имеет различное число уравнений.
В данном исследовании не преследовалась цель указать оптимальные значения параметров. Дело в том, что, с одной стороны, выбор их оптимальных значений зависит от конкретной задачи. С другой стороны, наблюдаются целые области значений параметров, при которых свойства метода существенно не меняются и для всех этих значений практически одинаково хороши.
С целью экономии места в данной работе приводятся лишь некоторые результаты из имеющихся в наличии. Остальные результаты аналогичны приведенным здесь. При изучении отдельных свойств метода в данной работе возникали ограничения на количество ячеек используемой сетки, поскольку для проведения исследований требовались более мощные вычислительные ресурсы, чем те, что имеются в наличии. Однако полученные результаты выборочно проверялись и подтверждались и на более мелких по сравнению с приведенными здесь сетках, т. е. сетках с большим числом ячеек.
1. Постановка задач
В качестве первой рассматривается краевая задача для уравнения Пуассона в прямоуголь ной двумерной области П = [х1,х1] х [х0,х1>]:
Дп = f, (х\,х2) € П, п|ап = U,
(1)
где Ди
д2и д2и
+ дж2'
дх\
В качестве второй рассматривается краевая задача для системы стационарных уравнений Навье—Стокса в области П:
диз
джі
+ М2+ -дХ - 7Г“Ди? = /, (жі,Ж2) Є П, і = 1, 2, дж2 дхз Не
,і ^^2
div(V) = 0, (х1,х2) Є П.
, VIдп = V, р(ж1,ж0) = Р
(2)
где функции р(х1,х2), п(х1,х2), ] = 1, 2 — давление и компоненты вектора скорости V соответственно; константа Ке — число Рейнольдса.
Расчетная область П покрывается равномерной сеткой П^2 с прямоугольными ячейками размера 2^4 х 2Л,2. В рассматриваемых здесь вариантах метода приближенное решение отыскивается в виде кусочно-полиномиальной функции на этой сетке. Для удобства реализации метода внутри каждой ячейки вводятся локальные переменные
Уі
ж1 — ж1т
У2
ж2 — ж2т
ІІ2 ІІ2
где точка (ж1т, ж2т) — центр т-й ячейки. Уравнение Пуассона после перехода к локальной системе координат имеет вид
1 д2м 1 д2 и
+ - — = /• (3)
^2 дУ?
^2 дУІ
Уравнения Навье—Стокса в локальных переменных записываются следующим образом:
1 др 1
1 дuj 1 дuj
—2 « ду1 + —2 «2 ду2 + Л-2 дyj
1 д 2и?-
Б,е \ —2 ду2
+
1 д 2и?-
Ь2 дУ2
Л, = 1, 2,
(4)
div(v) = 0.
Внутри одной ячейки сетки локальная координата у1 £
а у2 £ [— 11]-
Л1 Л1 —2 ’ —2_
В каждой ячейке приближенное решение ищется в виде линейной комбинации базисных элементов
и1
и
«2
р
12
'У ] ^3 т У] •, j=1
(5)
где т — номер ячейки. Решение задачи для уравнения Пуассона здесь в каждой ячейке аппроксимируется полиномом второго порядка. В этом случае базис состоит из шести элементов. При построении приближенного решения для уравнений Навье—Стокса компоненты вектора скорости в каждой ячейке будем аппроксимировать полиномами второго порядка, а давление — полиномами первого порядка. Кроме того, потребуем, чтобы каждый базисный элемент тождественно удовлетворял уравнению неразрывности. Нетрудно подсчитать, что в этом случае в базисе остается 12 элементов. Следовательно, получаемое здесь приближенное решение будет тождественно удовлетворять уравнению неразрывности. В табл. 1 и 2 приведены используемые здесь базисные элементы в вариантах метода КНК для уравнения Пуассона и уравнений Навье—Стокса соответственно.
Для определения неизвестных коэффициентов {^т} (5) в методе КНК для задачи (1) используются коллокации в четырех внутренних точках ячейки уравнения (3), умноженного на Ь|, условия согласования на границах между ячейками, краевые условия (для ячеек, примыкающих к границе дП). Для задачи (2) также используются коллокации в четырех внутренних точках ячейки первых двух уравнений системы (4), умноженных на Ь|, условия согласования и краевые условия.
Условия согласования записываются в точках с локальными координатами
± —1, ±С ) и (±£—1, ±1 ), граничные условия — в точках (± —1, ±£ ) и ( ±£—1, ±1
Ь2
—2
—2
—2
Таблица 1. Базисные функции для уравнения Пуассона
3 1 2 3 4 5 6
1 У1 У2 У2 22 У1 У2
Таблица 2. Базисные функции для уравнений Навье—Стокса
3 1 2 3 4 5 6
1 0 У1 У2 0 0
Vj 0 1 У2 0 У1 0
0 0 0 0 0 1
3 7 8 9 10 11 12
У2 22 У 0 -2 У1 У2 У2 + У2 -2 У1 У2
-2 У1 У2 0 У2 22 У -2 У1 У2 у2 + У2
2 У1/И,е 2 У1/Ие 2 У2/Ке 2 У2/Ке 0 0
ближенного решения задачи (1) в каждой ячейке выписывается переопределенная СЛАУ из 12 уравнений на шесть неизвестных, а для задачи (2) — СЛАУ из 24 уравнений на 12 неизвестных. В обоих случаях для удобства изложения введем для этой системы обозначение Ьм = Ь, где Ь — прямоугольная матрица, элементами которой являются коэффициенты при неизвестных {ajm} в указанных уравнениях. В качестве приближенного решения системы Ьм = Ь в рассматриваемом методе принимается решение, получаемое методом наименьших квадратов.
Как условие согласования решения на границе между ячейками в методе для уравнения Пуассона используется требование непрерывности в указанных на рис. 1 точках линейной комбинации функции м и ее производной по направлению внешней к ячейке нормали с положительными весовыми параметрами к1 и к2 соответственно:
где п — вектор единичной внешней нормали к границе ячейки, м+ и и- — значения функции и по разные стороны от нее.
Если сторона ячейки примыкает к границе области П, то на ней в двух точках, которые схематично изображены на рис. 1 черными кружками, записываются краевые условия. В методе для уравнения Пуассона краевые условия выглядят следующим образом:
жительный весовой параметр к3.
В качестве условий согласования в методе КНК для уравнений Навье—Стокса берутся требования непрерывности в указанных на рис. 1 точках следующих выражений:
&1М+ + &2^— = к1М + &2^Т,
дп дп
дм+ - дм
j
При записи переопределенной СЛАУ для задачи (1) уравнения (6) умножаются на поло-
+ П2^т=г - ПзР, П^т + П2^ дп дп
Точки записи уравнений коллокации
Точки записи граничных услови
У///?///////////
Точки записи
условий согласования
Рис. 1. Схема расположения точек записи граничных условий, условий согласования, уравнений коллокации в ячейке сетки
где п — вектор единичной внешней нормали к границе ячейки; ип, ит — нормальная и касательная компоненты скорости; Цj, 3 = 1, 2, 3 — положительные управляющие весовые параметры.
Граничные условия для уравнений Навье — Стокса реализуются по формулам, аналогичным формуле (6) для задачи (1). При записи переопределенной системы для задачи (2) каждое из них умножается на положительный весовой параметр п4. Уравнения коллокации в задаче (2) при записи переопределенной СЛАУ умножаются на положительный весовой параметр п5.
Вычислительные эксперименты показали, что от выбора значений введенных параметров зависит обусловленность возникающих при реализации метода КНК задач линейной алгебры, а также гладкость и точность получаемого приближенного решения.
Здесь для нахождения приближенного решения задач (1) и (2) используется метод итераций по подобластям (альтернирующий метод Шварца), в котором каждая ячейка сетки является подобластью. При решении задачи (2) линеаризация уравнений производится по Ньютону. Итерации по подобластям и по нелинейности совмещаются. Будем называть глобальной итерацией однократное последовательное выполнение итераций по всем ячейкам области П. Остановка итерационного процесса здесь происходила при выполнении условия
тах Кт1 - «Пт| < е, (7)
j т J 3
где {апт} — значения коэффициентов {ajm} на п-й глобальной итерации; е — малая положительная константа.
В реализованных ранее вариантах метода КНК [2-5] для решения переопределенной СЛАУ Ьм = Ь в каждой ячейке использовался метод наименьших квадратов. В этом методе от исходной СЛАУ переходят к определенной системе
ЬТ Ьм = ЬТЬ, (8)
где ЬТ — матрица, полученная транспонированием матрицы Ь. Число обусловленности V(ЬТЬ) матрицы ЬТЬ равно квадрату числа V(Ь). Здесь число обусловленности V(Ь) определяется, как отношение атах/атш, где атах — максимальное, а ат;п — минимальное син-
гулярные числа матрицы Ь. Уравнения системы (8) следуют из требования того, чтобы решение системы Ьм = Ь, получаемое в методе наименьших квадратов, являлось точкой минимума функционала невязки
ф(м) = ^ ((Ьм)г - ^
г
где (Ьм)г, Ьг — г-е компоненты векторов Ьм и Ь соответственно. Для того чтобы избежать большой аналитической работы при выводе формул метода, в [3] были введены два функционала Ф1, Ф2. Первый из них состоял только из суммы квадратов невязок уравнений, соответствующих граничным условиям и условиям согласования, второй — только из суммы квадратов невязок уравнений коллокации. Определенная система в работе [3] выведена из условий минимума функционалов Ф1 и Ф2, причем минимум Ф1 взят по коэффициентам ajm, 3 = 1, • • • , 10, при фиксированных а11т и а12т, а минимум Ф2 — по ajm, 3 = 11, 12, при фиксированных ajm, 3 = 1, ... , 10.
В данной работе для решения переопределенной СЛАУ Ьм = Ь предлагается использовать ортогональный метод с выбором главного элемента в столбце. Построенное таким
образом решение совпадает при отсутствии ошибок округления с тем, что получается методом наименьших квадратов при минимизации одного функционала, составленного из квадратов невязок всех уравнений системы Ьм = Ь. Алгоритм метода сводится к следующим действиям.
1. На 5-м этапе среди элементов /г5, где г > в (в =1, 2, ... , п2), матрицы Ь = {lij, 1 < г < п1,1 < 3 < п2,п2 < п1} ищется максимальный по модулю. Пусть к — это номер строки, в которой находится такой элемент. Производится перестановка строк с номерами в и к, а также перестановка соответствующих компонент вектора правых частей Ь.
2. Последовательно исключаются элементы, находящиеся под главной диагональю матрицы Ь в столбце с номером в, т. е. элементы {/г5, г > в}. Чтобы исключить элемент lis, матрица Ь умножается слева на матрицу вращения Qгs = {^т^}, которая получается из единичной путем замены ее элементов ^, С, С, на сов(<^), в1п(<^), - в1п(<^), сов^^) соответственно. Величина угла (pis выбирается таким образом, чтобы элемент ^^Ь)^ матрицы QгsL = {^гб1Ь)тп} был равен нулю, т. е. определяется равенствами
/ \ / \ lss
81Щ^) = —. , СОв^Л = —. .
Правая часть СЛАУ Ьм = Ь на каждом шаге ортогонального исключения тоже умножается на матрицу Qгs.
После приведения матрицы Ь к верхнетреугольному виду, т. е. к такому виду, что элементы lij = 0, при г > 3, из первых п2-х уравнений системы Ьм = Ь составляется определенная СЛАУ
Мм = (9)
Нетрудно показать, что решение полученной системы при отсутствии ошибок округления совпадает с решением по методу наименьших квадратов для СЛАУ Ьм = Ь.
В результате применения метода наименьших квадратов ухудшается обусловленность задачи линейной алгебры по отношению к исходной СЛАУ. Преимущество предлагаемого здесь подхода, основанного на применении ортогональных преобразований, заключается в том, что числа обусловленности матриц М и Ь совпадают при отсутствии ошибок округления.
2. Исследование влияния управляющих параметров на свойства метода КНК для уравнения Пуассона
При исследовании свойств метода КНК для уравнения Пуассона особое внимание здесь было уделено анализу обусловленности вспомогательных задач линейной алгебры, к которым сводится реализация метода. Были построены графики зависимости числа обусловленности матриц этих систем от различных пар управляющих параметров. В каждом численном эксперименте из шести параметров к1, к2, к3, и, £, £ выбирались два конкретных параметра. Их величины менялись в некоторой прямоугольной области на плоскости значений. В данном случае к1 = 1, к2 = 1. Величины остальных параметров фиксировались и принимали следующие значения:
к3 = 1, и = 0.5, £ = 0.5, ( = 0.5.
(10)
Если, например, выбраны параметры &1 и к2, то их величины меняются в ходе эксперимента, а значения остальных (к3, ^, £, £) задаются в соответствии с тем, как указано в (10).
Будем называть далее приграничными те ячейки, что имеют общие точки только с одной стороной прямоугольника, являющегося границей области П, угловыми — те, что имеют общие точки с двумя соседними сторонами того же прямоугольника, внутренними — остальные ячейки. Все переопределенные СЛАУ, которые необходимо решать на итерациях по ячейкам, будем называть локальными системами.
Исследована зависимость числа обусловленности матрицы Ь локальной системы для внутренних ячеек от параметров и к2 (рис. 2). Были также построены графики зависимости величины V(Ь) от параметров к1, к2 для приграничной (рис. 3, а) и угловой (рис. 3, б) ячеек сетки. Видно, что при малых значениях параметра к1 обусловленность систем во многих рассмотренных случаях в той или иной степени ухудшается. Исследована зависимость величины V(Ь) от параметров к1, к3 для локальных систем, соответствующих приграничной (рис. 4, а) и угловой (рис. 4, б) ячейкам. При малых значениях параметра к3 наблюдается существенное ухудшение обусловленности СЛАУ для угловой ячейки. Система при этом становится близкой к вырожденной из-за того, что уравнения, реализующие
Рис. 2. Зависимость величины V(£) от параметров &1 и ^2 для внутренней ячейки
Рис. 3. Зависимость величины V(£) от параметров &1 и ^2 для приграничной (а) и угловой (б) ячеек
N Щ
5
к3 7.5
N Щ
5
к3 7.5
Рис. 4. Зависимость величины V(£) от параметров &1 и кз для приграничной (а) и угловой (б) ячеек
а
граничные условия, имеют малый вес по отношению к остальным. По этой же причине наблюдается ухудшение обусловленности при малых значениях параметра к1.
В численных экспериментах было установлено, что величины параметров к1 и к2 влияют на непрерывность и гладкость получаемого здесь приближенного решения. Чем больше значение параметра к1 по сравнению с остальными двумя, тем меньше становятся разрывы приближенного решения. Аналогично, чем больше значение параметра к2 по сравнению с к1 и к3, тем меньше становятся разрывы производных приближенного решения. Это подтверждают рис. 5 и 6. В качестве тестового решения здесь была использована функция 8т(10ж + 10у). Расчет проведен на сетке 10 х 10. От выбора значения к3 зависит точность выполнения граничных условий. Чем больше значение параметра к3 по сравнению с к1 и к2, тем точнее выполняются граничные условия. Влияние параметра к3 на свойства получаемого здесь приближенного решения в приграничных и угловых ячейках более существенно, чем во внутренних. Во всех трех описанных здесь случаях погрешность выполнения одного из условий уменьшается одновременно с ростом величины соответствующего управляющего параметра. Однако при этом у остальных условий погрешность увеличивается. В предложенном здесь варианте метода КНК для уравнения Пуассона значение весового параметра уравнений коллокации было зафиксировано и положено равным единице. Если выбирать значения параметров к1, к2, к3 намного большими единицы, то в результате ухудшается точность, с которой удовлетворяются уравнения коллокации. Погрешность получаемого приближенного решения в этом случае будет велика. Таким образом, имеется возможность за счет ухудшения точности выполнения некоторых менее важных уравнений (условий) добиться уменьшения погрешности выполнения других путем выбора значений параметров к1, к2, к3. Это может быть полезно, например, когда точное решение задачи имеет особенность, в окрестности которой из-за чрезмерных требований гладкости в методе КНК будут с большой погрешностью выполнены уравнения коллокации. Устранить этот недостаток можно при помощи уменьшения параметров, управляющих гладкостью решения. Возможность использования такого приема продемонстрирована в конце разд. 3 на примере решения задачи о течении вязкой жидкости в каверне с движущейся верхней крышкой. Здесь погрешностью выполнения уравнения называется невязка, которая возникает при подстановке в него приближенного решения.
Значение параметра ш в данной реализации метода КНК решения краевой задачи для
Рис. 5. Фрагменты графика приближенного решения для девяти соседних ячеек сетки при следующих значениях параметров: а — &1 = 0.01, к2 = 1, &3 = 1; б — &1 = 1, к2 = 1, &3 = 1
Рис. 6. Фрагменты графика приближенного решения для девяти соседних ячеек сетки при следующих значениях параметров: а — ^ = 1, = 10-5, &3 = 1; б — ^ = 1, ^2 = 0.1, &3 = 1
уравнения Пуассона (1) не влияет на обусловленность возникающих задач линейной алге-
Я
0.0003-
0.0002-
0.0001
0 0.25 0.5 0.75 1
Рис. 7. График зависимости величины К от параметра и
бры, однако от него зависит величина погрешности приближенного решения Л, которая здесь вычислялась как максимум модуля разности значений точного и приближенного решений задачи (1) по всем узлам сетки П^1^2• Шаг Л,і = Л,і/10, шаг Л,2 = Л.2/10, где &н1н2 — сетка, на которой построено решение. График зависимости величины Л от параметра ш изображен на рис. 7. Значения остальных параметров выбирались так, как указано в формуле (10). В качестве тестового решения здесь использовалась функция ехр(х1 + х2). Эксперименты с другими тестами (например, с функцией віп(ж1 + х2)), а также при иных значениях остальных управляющих параметров привели к подобному результату. Расчет проведен на сетке 20 х 20. Исследования погрешности приближенного решения выборочно проводились и на более мелких сетках, например, на сетке 160 х 160 (табл. 3).
Также исследована зависимость от управляющих параметров числа обусловленности матрицы А глобальной СЛАУ, являющейся совокупностью систем вида (8), взятых из каждой ячейки сетки. При построении графика, изображенного на рис. 8, расчеты проводились на сетке 10 х 10. Матрица В другой глобальной СЛАУ, являющейся совокупностью систем вида (9), взятых из каждой ячейки сетки, отличается от матрицы А, но, как показали вычислительные эксперименты, области значений управляющих параметров, при которых число обусловленности у них существенно возрастает, в значительной мере пересекаются (эти области близки по размерам и по расположению на плоскости значений параметров). Результаты исследований приведены на рис. 8. Видно, что при к1 Є [2, 5], к2 Є [1, 5] число обусловленности мало меняется, и оно существенно меньше, чем при некоторых значениях параметров вне указанной области.
Численные эксперименты показали, что области значений параметров, при которых матрицы А, В, Ь хорошо обусловлены, в значительной мере пересекаются с областями, при которых наблюдается наилучшая точность приближенного решения, и с областями, при которых требуется меньшее количество итераций для выполнения условия (7). С целью экономии места соответствующие графики здесь не приводятся.
Рис. 8. Зависимость числа обусловленности V(А) от выбора значений параметров &1, к2 (а) и ( (б)
3. Исследование влияния управляющих параметров на свойства метода КНК для уравнений Навье—Стокса
При исследовании свойств метода КНК для уравнений Навье—Стокса здесь, как и в предыдущем разделе для уравнения Пуассона, особое внимание уделено анализу обусловленности задач линейной алгебры, к которым сводится реализация метода. Построены графики зависимости числа обусловленности глобальной матрицы В от различных пар управляющих параметров. В каждом численном эксперименте выбирались два конкретных параметра. Их величины менялись в некоторой прямоугольной области на плоскости значений. Величины остальных параметров фиксировались. При построении графиков, изображенных на рис. 9, 11-14 значения параметров выбирались следующим образом:
П1 = 1, П2 = 1, Пз = 1, П4 = 1, П5 = 1, ^ = 0.5, С = 0.5, ( = 0.5. (11)
Если, например, выбраны параметры п1 и п2, то их величины меняются в ходе эксперимента, а значения остальных задаются в соответствии с тем, как указано в (11). Для рис. 10 значения параметров п2, П3 равны 0.5.
Исследование обусловленности глобальной матрицы В требует существенно больших затрат ресурсов ЭВМ, чем вычисление приближенного решения на той же сетке. Поэтому сетки, для которых здесь приведены результаты исследований обусловленности матрицы В (самая мелкая сетка — 40 х 40), существенно грубее (имеют меньшее количество ячеек) сеток, для которых приведены значения приближенного решения (самая мелкая сетка — 160 х 160). При построении графиков, изображенных на рис. 9-13, использовалось тестовое решение
м1(ж1,ж2) = — 008(2пх1) 8т(2п(х2 — 0.5)), м2(ж1,ж2) = 8т(2пх1) 008(2п(х2 — 0.5)), р(х1,х2)= 0.5( 008(0.5п х1) + 008(0.5п х2)). (12)
Из графиков, изображенных на рис. 9, видно, что число обусловленности V(В) существенно возрастает при близких к нулю значениях параметров пз, П5. Причиной этого
Рис. 9. Зависимость числа обусловленности V(В) от параметров Пз, П5 при Ие = 1 на сетке 10 х 10 (а) и при Ие = 400 на сетке 20 х 20 (б)
Рис. 10. Зависимость числа обусловленности V(В) от числа Рейнольдса Ие на сетке 40 х 40 (а) и от параметров £, ( на сетке 20 х 20 (б)
Е Кр
0.018- 0.11
0.015- 0.09
0.012- 0.07
0.009- 0 0.25 ‘ . 0.5 0.75 ' 1 ш 1 0.05 0 ' 0.25 0.5 0.75 1
0.006- 0.03
аб
Рис. 11. Зависимость погрешности компонент скоростей Н (а) и давления Нр (б) приближенного решения от параметра и
Рис. 12. Зависимость числа N от выбора величин параметров Пз и П5 при значениях константы е = 2 ■ 10-5 (а) и е = 5 ■ 10-8 (б)
здесь, как и в методе КНК для уравнения Пуассона, является то, что при таких значениях параметров соответствующие уравнения в локальных переопределенных СЛАУ имеют малый вес по отношению к остальным. Существенный рост числа V(В) наблюдался также при уменьшении значений параметров пъ П2, П4. Исследовано влияние числа Рейнольдса И,е на обусловленность глобальной СЛАУ. Из графика, приведенного на рис. 10, а, видно, что с ростом числа И,е обусловленность существенно ухудшается.
Предпринимались попытки уменьшения количества уравнений в переопределенных локальных системах за счет отбрасывания некоторых условий согласования. Это приводило к росту числа обусловленности глобальной СЛАУ, к возникновению разрывов и ухудшению гладкости приближенного решения, к увеличению погрешности и к расходимости итерационного процесса. Также были предприняты попытки уменьшения количества уравнений за счет совмещения всех четырех точек коллокации в центре ячейки, что соответствует значению параметра ш = 0. Результаты численных экспериментов с некоторыми тестами показали, что при таком совмещении может возрасти погрешность приближенного решения по сравнению со случаем, когда ш = 0. При этом, однако, обусловленность глобальной СЛАУ мало зависит от параметра ш и не ухудшается при ш = 0. При построении графиков, изображенных на рис. 11, расчеты проводились с тестом (12) на сетке 20 х 20. Величина К" — это максимальная из погрешностей компонент скорости приближенного решения, Кр — погрешность давления. Для каждой из компонент скорости и для давления погрешность вычислялась здесь таким же образом, как величина К в разд. 2. Графики в экспериментах с тестом
и (Ж1,Ж2) = ехр(ж1 + Х2), V (Ж1,Ж2 ) = — ехр(х + Х2), Р (Ж1,Ж2) = ехр(х + х2)
получились подобными тому, что приведен на рис. 7. Таким образом, на точность численного результата параметр ш влияет по-разному в зависимости от рассматриваемого решения. Аналогичный результат был получен для параметров £, £.
На рис. 10, б приведен график зависимости числа обусловленности V(В) от параметров £, £. Видно, что глобальная СЛАУ плохо обусловлена, когда на всех сторонах каждой ячейки сетки две точки записи условий согласования совмещены в одну, что соответствует
Рис. 13. Зависимость погрешности компонент скоростей Н (а) и давления Нр (б) приближенного решения от выбора значения параметров Пз и П5
г
0.05- _ .........
0.04 - _ ■ • ‘ '
0.03- •
0.020.01—1---------1------1-----1-------|_ Нз
5-10-4 0.25 0.5 0.75 1
Рис. 14. Зависимость максимума величины г от параметра Пз при Ие = 1
значению параметра £ = 0. В численных экспериментах при £ = 0 итерационный процесс расходится. Для всех испытанных здесь тестов при значениях параметров ш = 0.5, £ = 0.5, £ = 0.5 получаемое приближенное решение имело хорошую точность. При другом выборе величин параметров погрешность не была существенно меньше. Видимо, указанные значения можно взять в качестве универсальных.
Были реализованы варианты метода КНК для уравнений Навье—Стокса, в которых использовалось большее количество условий согласования, чем в том, что описан в разд. 1. В экспериментах с гладкими тестами получаемое с помощью этих вариантов приближенное решение имело хорошую точность. Однако в случае, когда тестовое решение имеет особенность, такие варианты по точности оказались хуже того, что был описан ранее. Причина этого заключается в том, что из-за усиления требований гладкости происходит относительное уменьшение веса уравнений коллокации в переопределенной СЛАУ, и они в итоге выполняются с боольшей погрешностью.
На рис. 12 приведены графики зависимости числа N от параметров пз, П5 при двух различных значениях константы е, где N — это количество глобальных итераций, потребовавшихся для выполнения условия (7) при счете на сетке 10 х 10. На графике (рис. 12, а) при значениях параметров пз, близких к 0.8, и П5, близких к 0.05, есть существенный “провал”. При тех же величинах параметров пз, П5 на графике числа обусловленности на рис. 9, а “провал” не наблюдается. Причиной этого является нелинейность уравнений Навье—Стокса. Матрица В глобальной системы меняется в течение итерационного процесса. При недостаточно маленьких е матрица В из-за небольшого числа итераций существенно отличается от той, что была использована при построении графика числа обусловленности. Последний строился в результате проведения достаточно большого числа итераций. Поэтому возможны такие несоответствия упомянутых графиков. Однако при уменьшении числа е возрастает количество необходимых итераций, и, начиная с некоторой итерации По, величины коэффициентов разложения решения (5) становятся близкими к тем значениям, для которых выполнено условие (7). При этом элементы матрицы В меняются уже не так существенно, как на начальном этапе итерационного процесса, а скорость сходимости уже во многом определяется величиной числа обусловленности матрицы В, полученной после остановки итераций. И указанные несоответствия графиков должны быть моеньшими. Это подтверждает сравнение графиков, изображенных на рис. 9, а и 12, б.
На рис. 9 параметры пз £ [0,10] и п5 £ [0,10], а на рис. 12 и 13 пз £ [0,1],
п5 £ [0,1]. Построение в больших областях значений параметров графиков, аналогичных изображенным на рис. 12, 13, является трудоемкой работой, результаты которой не дают информации о свойствах метода КНК, отличной от уже имеющейся. Величины, графики
Таблица 3. Определение порядка сходимости приближенного решения метода КНК для уравнения Пуассона на последовательности сеток
г N1 х N2 Кг Лг-1 /Лг
1 20 х 20 7.2 ■ 10-4 —
2 40 х 40 9.9 ■ 10-5 7.3
3 80 х 80 1.4 ■ 10-5 7.1
4 160 х 160 1.7 ■ 10-6 8.2
которых приведены здесь, выборочно были вычислены и на более мелких сетках. При этом расположение областей наилучших значений параметров практически не менялось по отношению к тем, что имеются на указанных здесь рисунках.
Численные эксперименты для задачи (2), как и для задачи (1), показали, что области значений параметров, при которых матрица В хорошо обусловлена, в значительной мере пересекаются с областями, при которых наблюдается наилучшая точность приближенного решения, и с областями, при которых требуется моеньшее количество итераций для выполнения условия (7). Например, это видно при сравнении графиков, изображенных на рис. 9 (а), 12 (б), 13.
Численные эксперименты для задачи (2) показали, что параметры п1 и п2 влияют соответственно на непрерывность и гладкость компонент скорости в приближенном решении, а пз — на непрерывность давления, аналогично тому как это было относительно параметров к1, к2 в задаче (1). При решении задачи (2), как и в случае задачи (1), можно добиться уменьшения погрешности, с которой выполняются уравнения коллокации, за счет выбора значений параметров п1, п2, пз. В частности, это существенно проявляется, когда решение задачи имеет особенность. В качестве одного из тестов рассматривалось решение задачи о течении вязкой жидкости в каверне с движущейся верхней крышкой. На рис. 14 изображен график зависимости максимума величины г по всей расчетной области от параметра пз при И,е = 1. Здесь г — это невязка, возникающая после подстановки приближенного решения в исходные дифференциальные уравнения в точках коллокации. Использовалась сетка 10 х 10. Видно, что при ослаблении требования непрерывности давления путем уменьшения параметра пз повышается точность, с которой выполняются уравнения коллокации.
4. Численные эксперименты на сходимость приближенного решения
Проведены численные эксперименты на последовательности сеток для определения порядка сходимости приближенного решения, получаемого здесь при помощи метода КНК для задач (1) и (2). В этих экспериментах для уравнения Пуассона использовалось тестовое решение ехр(х1 + х2). Были взяты следующие значения параметров: к1 = 1, к2 = 10, кз = 7, ш = 0.5, £ = 0.5, £ = 0.5. В табл. 3 приведены значения погрешности приближенного решения Лг, которые здесь вычислялись так же, как величина Л в разд. 2.
Таблица 4. Определение порядка сходимости приближенного решения метода КНК для уравнений Навье—Стокса на последовательности сеток, Ие=1
Сетка Скорость Давление
г N1 X N2 п» гС^ -р» / -р» гг_1/ г г? г?_1/гр
Эта работа
1 10 X 10 2.0 10_2 — 2.1 10_ -1 —
2 20 X 20 3.1 10_з 6.5 3.3 10_ -2 6.4
3 40 X 4 О 5.0 10_4 6.2 6.1 10_ -3 5.4
4 80 X 8 О 9.2 10_5 5.4 1.3 10_ -3 4.7
Работа [3]
1 10 X 10 3.0 10_2 — 3.3 10_ -1 —
2 20 X 20 6.2 10_з 4.8 9.2 10_ -2 3.6
3 40 X 4 О 1.4 10_з 4.4 2.5 10_ -2 3.7
4 80 X 8 О 3.1 10_4 4.5 7.4 10_ -3 3.4
Таблица 5. Характерные значения компонент скорости на вертикальной (г1тщ) и горизонтальной (^2тт, ^2тах) средних линиях каверны
Метод Сетка г1 тт г2 тт г 2 тах
Ие = 1
Работа [3] 80 X 80 -0.2074 -0.1846 0.1839
Эта работа 80 X 80 -0.2077 -0.1849 0.1841
Ие = 100
Работа [6] 128 X 128 -0.2106 -0.2521 0.1786
Работа [7] 129 X 129 -0.21090 -0.24533 0.17527
Работа [8] 160 X 160 -0.2140423 — —
Работа [9] 40 X 40 -0.2203 -0.2614 0.1846
Работа [3] 80 X 80 -0.2125 -0.2497 0.1756
Эта работа 80 X 80 -0.2135 -0.2540 0.1798
Ие = 400
Работа [7] 129 X 129 -0.32726 -0.44993 0.30203
Работа [10] 128 X 128 -0.32751 -0.45274 0.30271
Работа [11] 49 X 49 -0.2960 -0.4051 0.2662
Работа [3] 80 X 80 -0.3222 -0.4443 0.2970
Эта работа 80 X 80 -0.3280 -0.4526 0.3022
Ие = 1000
Работа [6] 256 X 256 -0.3764 -0.5208 0.3665
Работа [7] 129 X 129 -0.38289 -0.51550 0.37095
Работа [10] 128 X 128 -0.38511 -0.52280 0.37369
Работа [11] 129 X 129 -0.3689 -0.5037 0.3553
Работа [8] 160 X 160 -0.388569 — —
Работа [3] 80 X 80 -0.3590 -0.4901 0.3463
Эта работа 160 X 160 -0.3850 -0.5228 0.3721
Для уравнений Навье—Стокса эксперименты проводились с тестовым решением (12). Результаты, приведенные в табл. 4, получены при И,е = 1. Использовались следующие значения параметров: ^1 = 1, П2 = 0.5, Пз = 0.5, п4 = 1, П5 = 1, ^ = 0.5, £ = 0.5, ( = 0.5.
Погрешности Я», Я? здесь вычислялись так же, как величины Я», Я? в разд. 3. Для сравнения приводятся численные результаты с тестовым решением (12), полученные в работе
[3] с применением другого варианта метода КНК, в котором давление также аппроксимировалось полиномами первого порядка. Из табл. 4 видно, что рассматриваемый здесь вариант метода КНК позволяет получать более точное приближенное решение. Видно, что с уменьшением шага сетки вдвое погрешность приближенного решения задачи (2) уменьшается не менее чем в 4 раза для компонент скорости и для давления. Следовательно, сходимость метода для гладких решений и умеренных чисел Рейнольдса не хуже чем второго порядка для скорости и для давления. Для задачи (1) наблюдается третий порядок сходимости.
Проведены также вычислительные эксперименты с решением задачи о течении вязкой жидкости в прямоугольной каверне с движущейся верхней крышкой при различных числах Рейнольдса: Ке =1, Ке = 100, Ке = 400, Ке = 1000 (табл. 5). В этих расчетах слабые вихри в левом и правом нижних углах каверны наблюдаются уже при И,е = 1. Это свидетельствует о том, что предложенный здесь вариант метода КНК обладает хорошей точностью. При вычислении величин, приведенных в табл. 5, использовались следующие значения параметров: п1 = 10, п2 = 1, П3 = 0.001 • Ке, п4 = 10, П5 = Ке, и = 0.5, £ = 0.5, £ = 0.5. Полученные здесь результаты хорошо согласуются с известными результатами, указанными в [3, 6, 7, 8, 9, 10, 11].
Заключение
В данной работе на примере решения краевых задач для уравнения Пуассона и уравнений Навье—Стокса рассмотрен новый вариант метода КНК. Он отличается от предыдущих вариантов [3, 4, 5, 12, 6] количеством управляющих параметров в формулах метода, способом решения переопределенных СЛАУ, возникающих при его реализации. При рассмотрении двух существенно различных краевых задач выявлено несколько общих свойств метода. В обоих случаях в численных экспериментах наблюдалось, что области значений параметров, при которых глобальная СЛАУ хорошо обусловлена, в значительной мере пересекаются с областями, при которых наблюдается наилучшая точность приближенного решения, и с областями, при которых требуется меньшее количество итераций для выполнения условия (7). Показано, что управляющие параметры, которые являются весовыми множителями при слагаемых в условиях согласования, одинаковым образом при решении обеих рассмотренных здесь задач влияют на гладкость приближенного решения во всей расчетной области и на величину невязки, возникающей после подстановки приближенного решения в исходные дифференциальные уравнения в точках коллокации. Также показано, что минимизация только одного функционала, составленного из квадратов невязок, при решении переопределенной СЛАУ в методе КНК в сочетании с подходящим выбором управляющих параметров может существенно улучшить точность приближенного решения, обусловленность задачи, скорость сходимости итераций в альтернирующем методе Шварца. Хорошее согласование результатов для задачи о течении в каверне с движущейся верхней крышкой, полученных здесь, с известными результатами других работ [3, 6, 7, 8, 9, 10, 11], полученными при помощи иных методов, свидетельствует о том, что новый вариант метода КНК позволяет получать решение с достаточно хорошей точностью. Таким образом, показано, что свойства метода КНК существенно зависят от управляющих параметров. Поэтому исследования, подобные проделанному здесь, необходимы для
успешного использования этого метода.
Список литературы
[1] de Boor C., Swartz B. Collocation at Gaussian points // SIAM J. Numer. Anal. 1973. Vol. 10, N 4. P. 582-606.
[2] Karamyshey V., Koyenya V., Sleptsoy A. Adaptive methods for Navier—Stokes equations // Computational Fluid Dynamics. Proc. of the Third ECCOMAS Conf. on Computational Fluid Dynamics. 1996. Vol. 1. (Paperback) P. 301-307.
[3] Семин Л.Г., Шапеев В.П. Метод коллокаций и наименьших квадратов для уравнений Навье—Стокса // Вычисл. технологии. 1998. Т. 3, № 3. С. 72-84.
[4] Семин Л.Г., Слепцов А.Г., Шапеев В.П. Метод коллокаций — наименьших квадратов для уравнений Стокса // Вычисл. технологии. 1996. Т. 1, № 2. С. 90-98.
[5] Shapeey V.P., Semin L.G., Belyaev V.V. The collocation and least squares method for numerical solution of Navier—Stokes equations // Proc. of the Steklov Institute of Mathematics, Suppl.2. 2003. P. S115-S137.
[6] Bruneau C.H., Jouron C. An efficient scheme for solving steady incompressible Navier—Stokes equations // J. Comput. Phys. 1990. Vol. 89, N 2. P. 389-413.
[7] 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, N 4. P. 387-411.
[8] Гаранжа В.А., Коньшин В.Н. Численные алгоритмы для течений вязкой жидкости, основанные на консервативных компактных схемах высокого порядка аппроксимации // Журн. вычисл. математики и мат. физики. 1999. Т. 39, № 8. С. 1378-1392.
[9] Шапеев А.В. Разностная схема четвертого порядка для уравнений Навье—Стокса на пятиточечном шаблоне // Динамика сплошной среды. Новосибирск, 2000. Вып. 116. С. 119-122.
[10] Deng G.B., Piquet J., Queutey P., Visonneau M. A new fully coupled solution of the Navier—Stokes equations // Intern. J. for Numer. Methods in Fluids. 1994. Vol. 19, N 7. P. 605639.
[11] Chen C.J., Chen H.J. Finite analytic numerical method for unsteady two-dimensional Navier— Stokes equations // J. Comput. Phys. 1984. Vol. 53, N 2. P. 209-226.
[12] Bochev P., Cai Z., Manteuffel T.A., McCormick S.F. Analysis of velocity-flux first-order system least-squares principles for the Navier—Stokes equations. Part I // SIAM J. Numer. Anal. 1998. Vol. 35, N 3. P. 990-1009.
Поступила в редакцию 21 ноября 2006 г.
