CC BY
56
Вычислительные технологии
Том 17, № 4, 2012
О некоторых проблемах конструирования разностных
>к
схем на двумерных подвижных сетках*
А. Ф. СоммЕР1, Н. Ю. ШокинА2 1 Новосибирский государственный университет, Россия 2Институт вычислительных технологий СО РАН, Новосибирск, Россия e-mail: nisei@sibmail.ru, nina.shokina@ict.nsc.ru
На примере схемы предиктор-корректор для линейного уравнения переноса с переменными коэффициентами, являющегося модельным для тестирования алгоритмов расчёта двумерных течений несжимаемой жидкости, обсуждаются общие проблемы, возникающие при построении адаптивных сеток и конструировании разностных схем на подвижных сетках.
Ключевые слова: численное моделирование, конечно-разностная схема, уравнение переноса, предиктор-корректор, геометрический закон сохранения, монотонная схема, адаптивная сетка, метод эквираспределения, результаты расчётов.
Введение
Для расчёта распространения длинных волн в морских акваториях и последующего наката таких волн на берег широкое распространение получили численные методы, основанные на аппроксимации уравнений модели мелкой воды. При приближении волны к берегу криволинейная линия уреза становится подвижной, поэтому решение задачи приходится искать в области не только сложной формы, но и с подвижной границей. Вследствие этих обстоятельств, процесс нелинейного взаимодействия длинных волн с реальной береговой линией считается трудным для численного моделирования.
В последние годы большое распространение получили алгоритмы решения задач наводнения-осушения на основе метода конечных элементов [1-6]. Тем не менее для таких задач продолжает оставаться актуальной проблема разработки новых и совершенствования существующих конечно-разностных методов, в частности, разностных схем на адаптивных сетках, так как на их основе тоже можно получать вполне адекватное представление о картине волновых процессов в областях c криволинейной формой подвижных границ [7-13].
Использование подвижных сеток, адаптирующихся к решению, может приводить к заметному повышению точности расчёта [14, 15]. Однако решать задачу на адаптивной сетке значительно сложнее, чем на равномерной, поскольку при использовании адаптивных сеток приходится кроме основных уравнений задачи дополнительно решать нелинейные уравнения для определения координат узлов, причём для нестационарных задач уравнения для сетки необходимо решать на каждом временном шаге. При конструировании схем на адаптивных сетках возникают и другие трудности, на которые
* Работа выполнена при финансовой поддержке РФФИ (гранты № 10-05-91052-НЦНИ и 12-01-00721), Совета по грантам Президента РФ по государственной поддержке ведущих научных школ РФ (грант № НШ-6293.2012.9) и Проекта IV.31.2.1. программы фундаментальных исследований СО РАН.
в публикациях, как правило, не обращается внимания. В настоящей работе описываются практические рецепты преодоления такого рода трудностей. Основное внимание уделено тем вопросам, которые не были затронуты в предыдущих публикациях или были недостаточно подробно освещены. Приводится, в частности, конечно-разностная схема второго порядка аппроксимации на адаптивной сетке, выводится геометрический закон сохранения в разностной форме, рассматриваются особенности построения сеток, адаптирующихся к разрывным решениям, обсуждаются реализации краевых условий на подвижной сетке.
Обсуждаемые проблемы демонстрируются на примере разностных схем для простейшего двумерного линейного уравнения переноса
I+dur+tr=°, <> ° (D
с переменными коэффициентами u = u(x,y) и i = i(x,y). Тем не менее сформулированные рекомендации и выводы имеют общий характер и могут быть полезными, например, при численном моделировании течений несжимаемой жидкости в рамках плановой модели мелкой воды. Отметим, что уравнение переноса (1) очень часто используется для тестирования новых конечно-разностных схем (см., например, [14, 16]), а также для сравнения эффективности различных численных методов.
Настоящую статью можно рассматривать как вторую часть работы [17], посвящён-ной анализу проблем построения адаптивных сеток для решения одномерных нестационарных задач и конструированию дивергентных разностных схем, сохраняющих монотонность численного решения.
1. Постановка задачи
Уравнение (1) является модельным для тестирования алгоритмов расчёта двумерных течений несжимаемой жидкости, а его коэффициенты u и i представляют собой аналоги компонент вектора скорости. Поэтому далее будем использовать термины гидродинамики, называя заданный вектор u = (u,i) вектором скорости и предполагая, что его компоненты удовлетворяют уравнению неразрывности модели несжимаемой жидкости
du di ^
ш + di = °, x = (x,y) * «, (2)
где П — односвязная область, представляющая собой криволинейный четырёхугольник с неподвижной границей Г, состоящей из четырёх частей (левой, нижней, правой, верхней):
Г = Глев и Гн U Гпр U Гв. Равенство (2) позволяет переписать уравнение (1) в недивергентной форме
dr dr dr
Ж + udX + vdy = °. (3)
Предполагается, что через участок границы Глев "жидкость" втекает в область П, через Гпр — вытекает, а участки Гн и Гв являются непроницаемыми стенками. Таким образом, вектор скорости удовлетворяет соотношениям
u ■ n \ < u ■ n \ > u ■ n I = (4)
1хеГлев ' 1хеГпР ' 1хеГниГв ' v '
где п — единичная внешняя нормаль к границе Г. Из соотношений (4) следует [18], что краевое условие для уравнения (1) требуется задавать только на участке втока Глев:
^(М) г = ^лев (х,г), г > 0. (5)
1 лев
Кроме (5) предполагается заданным начальное условие
р(х, 0) = ^о(х), х е П = П и Г. (6)
Отметим, что для выполнения уравнения неразрывности (2) необходимо, чтобы на границе Г вектор скорости удовлетворял условию
и ■ п ¿в = 0, (7)
г
где в — длина дуги Г, отсчитываемая от некоторой точки при обходе границы против часовой стрелки. С учетом последнего из соотношений (4) равенство (7) принимает следующий вид:
— J и ■ п ¿в = J и ■ п ¿в. (8)
Глев Гпо
2. Схема предиктор-корректор на равномерной сетке
В работе [17] была исследована схема предиктор-корректор, аппроксимирующая уравнение переноса с одной пространственной переменной. Обобщим схему [17] на двумерный случай, причём вначале рассмотрим двумерную схему на прямоугольной равномерной сетке с узлами и шагами Н\, к2. Фрагмент этой сетки показан на рис. 1. Далее для сокращения записи будут использоваться обозначения, указанные на данном рисунке, согласно которым узлу присвоен локальный номер 0, узлу Хj1-1,j2 — 1, середина отрезка, соединяющего узлы и х^1+1,^2, обозначена буквой Е, центр ячейки хл+1/2,^2+1/2 — буквой С и т.д.
8 I WN
—о—
81
nwl п
SW)
5Т
4 EN
-о-
N
С
г—<?—?
Щ 0 е
6----Ь----^
—о-WS
-Л
и
-о-
^ ES
\т
В TSE
ь-
Рис. 1. Локальная нумерация узлов и обозначения для центров и середин сторон ячеек прямоугольной равномерной сетки
3
1
2
S
Уравнение переноса, записанное в двух формах (3) и (1), аппроксимируется на шагах предиктор и корректор соответственно. При этом, как и в одномерном случае [17], используется сетка с разнесенными узлами: предикторные величины вычисляются в центрах ячеек, а корректорные — в их вершинах х^-. Например, для ячейки с центром С (см. рис. 1) естественным обобщением шага предиктор из [17] будет уравнение
УС УС + ((1 + 01 + ((1 + 02Х)С = 0, (9)
т/2 V ^ V 2' rVo
где
yg + yg + yg + yg (in) Уа =-4-> (1П)
yg + y? - yg - yg y? + y? - yg - yg
y-a =-2h-' yg'C =-2h-• (11)
Предполагается, что сеточная функция u (разностный аналог первой компоненты вектора скорости) определена в серединах вертикальных, а v — горизонтальных сторон ячеек, поэтому в центрах ячеек компоненты скорости находятся по формулам
un + uNE ve + vew
ua = —2— ' va = —2— •
Корректорные величины yg+1 во внутренних узлах сетки вычисляются на основе разностного уравнения
yg+1 - yg + (uy*)E - (uy*)w + (vy*)N - (vy*)g = 0 (12)
т hi h2 '
где
(uy*)B + (uy* )c (uy * )a + (uy * )d
(uy )e =-2-' (uy =-2-'
, (vy * )d + (vy * )a , *, (vy * )a + (vy * )B
(vy )n =-2-' (vy =-2-•
Нетрудно показать, что если 0a = 0(т) (a = 1, 2), то на гладких решениях схема (9), (12) имеет второй порядок аппроксимации. При 0a = 0 она переходит в двухшаговую схему Лакса—Вендроффа.
Схемные параметры 0a (a = 1, 2) выбираются так, чтобы при нулевом значении одной из компонент скорости и постоянном значении другой (u = 0 и v = const или v = 0 и u = const) схема (9), (12) совпала со схемой [17] для одномерного уравнения переноса. Такое совпадение будет иметь место, если схемные параметры вычислять, например, по следующим формулам:
0 при |ga,a1 < |ga,a|' ga,a ■ ga,a > 0'
7«,a = \ I 1 - g^ ) при |ga,a1 > ^a.C|' ga,a ■ ga;a > 0' (13)
y«,a /
^,C при ga,a ■ ga,a < 0'
где
1т
0a,a = - 1' Ka = — |ua| ' ga = |ua| (1 - Ka) Ух
a
Ка < 1, и1 = и, и2 = V, х1 = х, х2 = у, С — центр ячейки, соседней к ячейке с центром С:
_ Г хл+1/2-«1 ,.72 + 1/2 пРи а = 1 С = < = sgn .
[ х.1+1/2,.2 + 1/2-«2 ПРИ а = 2
Схему (9), (12) можно рассматривать как каноническую форму явных двухслойных схем для двумерного уравнения переноса (1) в том смысле, что все они могут быть записаны в виде (9), (12) при подходящем выборе формул для вычисления величин /п в центрах ячеек и параметров Например, при ва = 0 схема (9), (12) совпадает со схемой Лакса —Вендроффа, а при
^ п п _ п0 га _ +
и > 0, иа = , = —^—
получается следующий вариант противопоточной схемы:
( — V МгК — 7)о(0п!
^п+1 /гП 1
^0 — Ро +1
т 2
ив/2 — и А+ ис—
1
+ 2
зд — гА/5 + "с/0 — ^в/2
Л-2 Л-2
Выбирая ту или иную формулу для параметров , можно получить явные схемы с различными свойствами: первого или второго порядка аппроксимации, абсолютно или условно аппроксимирующие, условно устойчивые или абсолютно неустойчивые. Далее при использовании схемы на подвижных сетках будут использоваться аналоги формул (10) и (13).
3. Уравнение переноса в новых независимых переменных
Как и в одномерных нестационарных задачах [17], в двумерном случае для построения разностной схемы на подвижной сетке нужно сначала записать исходную задачу в вычислительной области. Построение подвижных сеток основано на преобразовании координат
х = х^1,д2 ,£), у = у(д\д2,*), (14)
которое в каждый момент времени £ устанавливает взаимно-однозначное непрерывно дифференцируемое соответствие между физической областью П с криволинейными границами и вычислительной областью Q простой формы — единичным квадратом в плоскости д1Од2. Будем предполагать, что части Глев, Гн, Гпр и Гв границы Г являются образами при отображении (14) соответственно левой 7лев, нижней 7н, правой 7пр и верхней 7в сторон квадрата Q, граница которого обозначена через 7 =
В настоящем разделе мы ограничимся рассмотрением в новых координатах д1, д2, £ только уравнения переноса. В силу равенств
/(х,у,£) = /(x(q1,q2,t), у(?1 ,д2= = /(У^у^ 52(х,у,£)
уравнение (3), записанное относительно функции /(д1 ,д2,£), примет вид
Ж + ^ + ^ = °, (15)
где V" (а =1, 2)— контравариантные компоненты вектора скорости (д1 дд1 дд1 дд1 2
V1 =
(И
дд1 + дд1 + дд1 д£ дх ду
а
Учитывая равенства дд1 1 ду
дх 3 дд2
дд2 дх
1 дУ
Здд1
дд1 дУ
дд2 + дд2 + дд2 д£ дх ду
1 дх дд2 1 дх 3 дд2' ду 3 дд1;
:1б)
получаем
V1-
дд1
1 3
ду дх дд2 дд2
V2-
дд2
1 3
ду + дх дд1 дд1
:17)
или
и = 3
и
дх \ ду / ду \ дх
- д^Уд^2 - V - Ж Уд^2
V = —
/ дх \ ду + / ду \ дх
д£ / дд1
д£ / дд1
где 3 — якобиан преобразования (14)
3
дх ду дх ду
дд1 дд2 дд2 дд1
> 0,
и использованы равенства
дд1
~Ж
1 /ду дх дх ду 3 \д£ дд2 д£ дд2
дд2 1 / дх ду
"Ж = 3 V д^1
ду дх д£ дд1
:18)
Уравнение неразрывности (2) можно записать в новых координатах в следующем
виде:
дд1
т1 1 дд1
+
дд2
Т1 2 дд2
или с учетом тождества
д3 д / дд1
Ж + д^1 V
А
дд2
+ 13
как
д3 д3г>1 д3^2 д£ дд1 дд2
дд2
0.
19)
(20)
(21)
Следовательно, (15) можно записать в дивергентной форме
+
дд1
+
дд2
0.
(22)
Замечание 1. Тождество (20) выполняется для произвольного отображения (14) и называется геометрическим законом сохранения. Его можно переписать в эквивалентной форме
д 3 д / дх ду ду дх \ д / дх ду ду дх д£ дд1 \ д£ дд2 д£ дд2/ дд2 \ д£ дд1 д£ дд1
0.
(23)
1
0
0
При построении разностных схем на подвижных сетках необходимо добиваться [17], чтобы разностный аналог этого закона также выполнялся (см. ниже раздел 6).
Замечание 2. На границе 7 = д^ контравариантные компоненты (16) вектора скорости удовлетворяют условиям
V1 > 0, V1 | > 0,
V2 | ,, = 0,
17н и 7в 7
' Тлев ' ОЛр
а равенство (8) принимает в новых координатах следующий вид:
(24)
1
т( 1
1
«*2=/' (V1 - £
о
д!=0
(25)
«1=1
Замечание 3. В силу предположения о неподвижности границы Г ливы равенства | | | |
дд2
дд1
7ле
дд1
~Ж
0,
7пр
дд2
~Ж
0.
дП справед-(26)
7
7
н
в
4. Схема предиктор-корректор на подвижной сетке
Аналогом схемы предиктор-корректор (9), (12) в случае подвижной сетки будет схема, на шаге предиктор которой аппроксимируется недивергентное уравнение (15), а на шаге корректор — дивергентное (22). Первое из уравнений аппроксимируется в центрах ячеек равномерной прямоугольной сетки покрывающей вычислительную область а второе — во внутренних узлах qj этой сетки. Совокупность узлов криволинейной сетки, покрывающей область П, обозначена через П^. Узлы хп сетки являются образами узлов q^• сетки при отображении (14), при этом ] = (^'1,^2), qj = , = .7«^«,
^ = 0,...,Ж„, = 1/Ж„, а =1, 2.
Поскольку положение узлов сетки П^ меняется со временем, то с каждым узлом связана скорость его перемещения
х = хГ1 - хП
' т
а с каждой ячейкой сетки ^ — сеточный аналог якобиана отображения (14)
тга ( \п
'.7+1/2 = Iх«1 у«2 — Х«2 У«^ ^+1/2,
где разностные производные х«1, х«2, у«1, у«2 в центрах ячеек вычисляются с использованием центральных разностей
х«1,с =-т-, х«2,с =---, (27)
Д2
при этом в серединах сторон применяется осреднение
хо + х4 хз + х7 хо + хз х4 + х7
хм = —2—, = —2—, хе = —2—, = —2—
(здесь использованы обозначения рисунка 1 для узлов, середин сторон и центров ячеек сетки Qh).
Прежде чем выписать аппроксимацию уравнения (15), поясним, как следует вычислять в центрах ячеек коэффициенты V1, V2, чтобы уравнение неразрывности (21) (или (19)) было справедливым и для сеточных функций. Если на данном слое по времени ввести функцию тока — (д1, д2)
5д1\
3 V1 — ^ = 3 V2 — = — ^
5д
5д2
5д*
(28)
то уравнение (19) выполнится автоматически, а для определения функции — получится уравнение
5 /#22 _ #12 + + #11 А = _
5д1 \ 3 5д1 3 5д2/ 5д2 \ 3 5д1 3 5д2/
(29)
с известной правой частью
1 / 5^2
3 \ 5д2 5д1
где V« — ковариантные компоненты вектора скорости
V« = #«0 + #а1 V1 + 9а2У2, а = 1, 2,
#«0, , — компоненты метрического тензора
#а0 = х4хда + у4уда , #11 = х^1 + у^1, #22 = х^2 + у^2 , #12 = #21 = х^1 хд2 + у^1 уд2.
Для уравнения (29) ставится задача Дирихле. С учетом (26) и последнего из соотношений (24) можно положить
—(д1, 0) = 0, — (д1,1) = —в = ео^, 0 < д1 < 1.
(30)
Используя формулы (17), (26), (28), получаем граничные значения для — на двух других участках границы 7:
^п, — (1,д2)
д1=0
5у 5х \
5д2 5д2 У
91 = 1
(31)
Заметим, что выполнение условия (3) позволяет вычислить постоянную фв
/( 5у 5х \ № —
5д2 5д2
^д2
1=0
5у 5х \
--V-)
5д2 5д2у
^д2 = —(1,1).
1=1
Решив численно задачу (29)-(31), найдем функцию тока —п на п-м слое по времени и определим в серединах сторон ячеек величины из левых частей равенств (28). Например,
т{ 1 5д1
3Г — "Ж
N
ФП — —п _ -п
3 V2-
2 5д2
5£
Е
^^ = ——п1,е . (32)
2
2
1
1
Полученные выражения используются при вычислении контравариантных компонент скорости также в серединах сторон ячеек, например,
= (^ - + , = ( - ^ + - У^х^!)е, (33)
при этом производные хда в серединах сторон ячеек вычисляются по формулам вида (32), например,
п _ х4 х0 п х3 х0
х„2 N — -:-, х„1
а скорости узлов — следующим образом:
_ х4,0 + х4,4 _ х4,0 + хМ
х4,М — -^-, х*,е = -2-'
И наконец, определяются контравариантные компоненты вектора скорости в центрах ячеек. Например, для центра ячейки, обозначенного на рис. 1 буквой С, полагаем
)С — + , — . (34)
С учетом принятых соглашений уравнение шага предиктор запишется в виде
+ ((1 + ^41) П + ((1 + 02)^2) П — 0, (35)
где величина ^С определяется как среднее арифметическое (10), а разностные производные в центрах ячеек вычисляются по формулам вида (11).
Величины ^П+1 вычисляются на шаге корректор из конечно-разностного уравнения, аппроксимирующего дивергентное уравнение (22):
(^)п+1 - (^)п + (^У)Е - + ^У^ - — 0 (36)
т Кг К2
где
тп I тп \ тп I тга
^ — + 'б + 'с + 'в , (37)
— (^)Б + (^)С, ^^ — + , (38)
(^2— ('7^)в +(, — (+(, (39)
а величины (3^°)* в центрах ячеек полагаются равными среднему арифметическому величин (34) на временных слоях п и п +1, например,
— ЗЗЦЗ!^. (40)
При описании схемы (9), (12) на равномерной сетке было указано, что схемные параметры 0а выбираются аналогично схеме [17] для одномерного уравнения переноса. В работе [17] была предложена также формула для параметра 0, при применении которой схема предиктор-корректор сохраняет монотонность численного решения на одномерной подвижной сетке. В настоящей работе аналог этой формулы используется в схеме предиктор-корректор (35), (36) на двумерной подвижной сетке, а именно, параметры 0а (а — 1, 2) вычисляются по формуле (13), в которой полагается да — 3 (1 - Ка) ,
^ Т. а| - / Ял+1/2-*и2+1/2 при а — 1 „
К — — , С — < — Sgn .
Ка I 4^1+1/2,^2+1/2-82 пРи а — 2,
5. Вычисление граничных значений
При выполнении шага предиктор в приграничных ячейках требуются значения сеточной функции / в граничных узлах я. Е 7^, поэтому необходимо указать способ определения этих граничных значений на каждом шаге по времени. На прообразе входа 7лев значения /п+1 определяются из граничного условия (5). В других граничных узлах я. Е 7^,н и 7^,пр и 7^,в краевые условия для функции / не заданы, поэтому для вычисления значений <^п+1 необходимо иметь какое-то разностное уравнение. Воспользоваться уравнением (36) нельзя, поскольку оно справедливо только для внутренних узлов сетки. Разностное уравнение для граничных узлов я. Е 7^,н и 7^,пр и 7^,в можно получить, используя следующий приём.
Возьмем разложение
/ (я.,£п+0 = / ^ ,*п) + т/ (я.,*п) + у(я.,*п) + °(т3)
(41)
и с помощью уравнения (15) исключим из него производные по времени. Для гладкого решения данного уравнения справедливы следующие выражения:
V1 /д! — V2 ,
= — VlVg1 — V1/q1í — V2/д2 — V2/q2í, = — + V2 ^ д1 , = — (V1 + V2 ^ д
Поэтому для решения уравнения (15) выполняется равенство
/ (я. ,£п+0 — / (я. ,£п)
+
т
1 2 т / 1 2
V + V + 2 + V* —
— + ^^2) д1 — + ^^2) ) (я. , £п) = 0(т 2).
Данное равенство используется для получения разностных уравнений в граничных узлах я. Е 7^,пр, при этом производные по переменной д1 заменяются односторонними разностями второго порядка. В силу (24), V2 = 0 на участках 7н и 7в, поэтому разностные уравнения в узлах я. Е 7^,н и 7^,в выводятся на основе более простого равенства
/ (я. ,£п+0 — / (я. ,£п)
+
т
^^ + ^^— V) ,£п)=0(т2)
6. Геометрический закон сохранения на подвижной сетке
Как было указано в [17], выполнение геометрического закона сохранения (далее г.з.с.) является необходимым условием согласованности разностных формул, связанных с вычислением характеристик подвижной неравномерной сетки, и означает не только согласованность разностных формул для якобиана и скоростей движения узлов, но и гарантию того, что разностная схема, как и аппроксимируемое этой схемой дифференциальное уравнение, имеет в качестве своего точного решения постоянную функцию, что в свою очередь отвечает общему требованию о сохранении разностной схемой как можно большего числа свойств аппроксимируемого дифференциального уравнения. Г.з.с. означает, что ячейка на новом временном слое полностью определяется ее положением на
2
2
текущем слое и скоростями движения границ. Таким образом, выполнение этого закона — гарантия того, что движение сетки не вносит искусственные изменения в область решения задачи.
Г.з.с. был впервые сформулирован в работе [20] (где был назван пространственным законом сохранения (space conservation law (SCL)) как дополнительное уравнение к законам сохранения массы, момента и энергии. В [21] закон был переоткрыт и записан в виде ограничения на численные аппроксимации законов сохранения. Обзоры работ по подходам к реализации г.з.с. даны, например, в [22, 23].
Так как постоянная функция является решением дифференциального уравнения (1), то естественно потребовать, чтобы построенная схема предиктор-корректор также была точна на постоянных функциях. Покажем, что для схемы (35), (36) выполняется разностный г.з.с., из чего будет следовать сохранение данной схемой постоянной функции.
Лемма 1. В каждой ячейке сетки выполняется разностный аналог геометрического закона сохранения (23):
jn+1 _ Jn
—C + (- +y^) qi + (x*y« 1 - ^Xl1) q2,c =0' (42)
где использованы обозначения рис. 1, разностные производные в центре ячейки вычисляются по формулам вида (27) и величины с символом "*" получены осреднением по времени:
xn | xn+1
x* = q , q . (43)
q
2
Доказательство. На примере разностной производной (ж«У?2 ),1 выведем формулу для производной от произведения двух функций:
(xty,
q we v q / n
q1,C
hi
_ + (ж«)^ (у?2 - (у?2 )N + (у?2 + (у?2 )N (ж^Д - (ж^)^
= 2 2 ' В соответствии с формулами (27) полученное равенство означает, что
(ж^ )91>с _ (ж4У?2,1 + У?2 (ж V )с ' (44)
Рассмотрим вторые производные в правой части (44):
? _ _ У7 - У? У? - У?
у,2^ _ и _ иД ^ и
_ 1 (У? - У? У? - У? ^ _ (^Оем—^О-Д - 7
_ ^ йГ) _ и _ , ^'
т. е. имеет место равенство смешанных производных
У?2,1,С _ У?1,2,С' (45)
Нетрудно доказать, что выполняется и следующее равенство смешанных производных:
(46)
В самом деле,
N _ 1 I
91,С
К1
К1
п+1 п
'X дг 'X
N
т
1 / 'п+1 _ хп+1 1 / ^
К1
К1
хп+1 _ хп
— К1)
4,С
Учтём равенства (45), (46) и перепишем формулу в виде (44)
('^О д1,с — ('^д2 + У^2 (х91
-г1-2 + Уп2 (х91
(47)
Очевидно, что аналогичная формула имеет место и на (п+1)-м слое по времени, поэтому с учётом обозначения (43) получим
('¿У*2)д1,С — ('У*1д2 + У* 2 (х91
(48)
Используя последнее равенство и его аналог для производной по переменной д2 от произведения двух функций
'Уд^ 92,С — У*1д2 + У*1 ('д2 )0С
(49)
можно выполнить следующие преобразования:
( - '¿У* 2 + 91,с + (х4у* 1 - д2 ,с — - 'У*1д2 - У*2 (хд1 )4 +
+у4ж!1„2 + '«2 (Уд1 )4 + 'У* 1-2 + У* 1 ('д2 )4 - У4х91„2 - <1 (Уд2 )4
с
- У*2 (х91 )4 - х-1 (Уд2 ^ + х-2 (У-1 )4 + У* 1 ('д2 )4
V ^ Уд
д2 V уд1
V ^ д2
с
Уп + Уп+1 хп+1 хп хп + хп+1 Уп+1 Уп Уд2 + Уд2 'д1 'д1 'д1 + 'д1 Уд2 Уд2
т
тп + 'п+1 Уп+1 Уп Уп + Уп+1 'п+1 -п + 'д2 + 'д2 Уд1 Уд1 + Уд1 + Уд1 'д2 'д2
2
т
2
т
с
т
- ('д1 Уд2 )п+1 + ('д1 Уд2 )п + ('д2Уд1 )п+1 - ('д2 Уд1 )Г
тп+1 _ тп
3С 3С
с
т
Таким образом, разностный геометрический закон сохранения (42) действительно выполняется.
2
2
1
Лемма 2. В каждой ячейке сетки ^ выполняется разностный аналог уравнения неразрывности (21):
З^1 - 3
п
с Зс + (л1):.,0+с»2):.,с—о, (50)
т
где
(З»а)
_ _ (3»а)п + (3»а)п+1
Доказательство. Согласно определению (27) производной в центре ячейки и формуле (33) для контравариантных компонент скорости, имеем цепочку равенств
('»Ч д1,С + (3»2) д2,С
•**2- ж4У*2 + У^д^ д1,с + ( - 1 + ж4У* 1- у^) д2 ,с —
— ^*2д1,С - 1д2,С + ( - 'У*2 + У'дОд1,С + ('¿У* 1 - У^д0д2,С.
Учитывая теперь геометрический закон сохранения (42) и равенство "*2д1 с — 1д2 с, убеждаемся в справедливости разностного уравнения неразрывности (50).
Лемма 3. В каждом внутреннем узле сетки ^ выполняется следующий разностный аналог уравнения неразрывности (21):
зп+1 - 3оп + + — 0, (51)
т К К2
где использованы обозначения (37)-(40).
Доказательство. Запишем уравнение неразрывности (50) для четырёх ячеек с центрами С, В, Д и А:
3С+1 - ЗС + , (З»2)^ - (3»2)Е — 0,
Т К1 К2
3Б+1 - зб , - , (3»2)Е - — 0,
т К1 К2
3В+1 - ЗВ + — 0,
т К1 К2 '
3А+1 - ЗА , (З»1 & - З»1)^ , (З»2- (З»2)^ — 0.
т К1 К2
Сложив эти уравнения, поделив на четыре и приняв во внимание обозначения (34) (37)-(40), получим разностное уравнение (51).
2
Теорема 1. Разностная схема (35), (36) сохраняет постоянную функцию.
Доказательство. Пусть = M = const. Тогда из уравнения (35) получим, что в центре любой ячейки ^*+1/2 = M. Следовательно, на шаге корректор (36) получается уравнение
Jcn+1 - JonM + M (Jv% - (Jv1)^ + M (Jv2)N - (Jv2)S =0.
T h h-2
Перепишем это уравнение в эквивалентной форме
Jn+1(^n+1 - M) +
T
, M J1 - Jn , (Jv1)^ - (Jv1 )W , (Jv2)N - (Jv2 )s , =0.
V T h h '
Учитывая уравнение неразрывности (51), получим, что = М. Таким образом,
разностная схема (35), (36) точна на постоянных функциях.
7. Метод эквираспределения для построения двумерных неравномерных подвижных сеток
В нестационарных задачах адаптивная сетка должна перестраиваться на каждом шаге по времени, поскольку решение задачи зависит от времени. Укажем используемый нами способ вычисления координат узлов х™+1 на (п + 1)-м слое в предположении, что на п-м слое по времени сетка уже имеется и на ней вычислено решение .
Для определения местоположения узлов х™+1 можно воспользоваться методом эквираспределения [24]. Однако при решении одномерных задач было обнаружено [17], что использование классического метода эквираспределения для построения подвижных адаптивных сеток может приводить к осцилляциям траекторий узлов адаптивной сетки и немонотонности изменения отношений длин соседних ячеек сетки. Кроме того, решение классического уравнения очень чувствительно к возмущениям управляющей функции и>, и построенная сетка может сильно осциллировать по времени даже в том случае, когда решение меняется слабо. Поэтому в работе [17] было рекомендовано добавить в уравнение классического метода эквираспределения правую часть, отвечающую за гладкость траекторий узлов сетки. Данной полученной эмпирическим путём рекомендацией воспользуемся и в двумерном случае: для построения подвижных сеток будем применять двумерный аналог одномерного эволюционного уравнения из [17]:
д / дх\ д ( дх\ дх
дд! Г^22д?) + д? дф) = вЖ> (52)
где в — положительный параметр, подбираемый экспериментальным путем в целях уменьшения осцилляций траекторий узлов сетки. При малом значении параметра в влияние члена вх незначительно, но при больших в величины смещений узлов уменьшаются и сетка становится "малоподвижной".
Далее дифференциальное уравнение (52) заменяется конечно-разностным
х™+1 _ х™
Лхп+1 = в * т * , Я, е Ян, (53)
где т — шаг по времени, а оператор Л отличается от оператора стационарного (классического) уравнения [19] только тем, что теперь в его коэффициентах ковариантные компоненты метрического тензора да«1 берутся с (п+1)-го временного слоя, а управляющая функция адп вычисляется по решению с п-го слоя. При в > 0 наличие в разностных уравнениях (53) правой части не только способствует уменьшению осцилляций траекторий узлов, но и увеличивает диагональное преобладание системы этих уравнений, что положительно влияет на скорость сходимости итерационного метода переменных направлений. В качестве начального приближения для итерационного процесса берётся сетка хп. Заметим, что начальная сетка х0 должна быть адаптирована к начальной функции (6). Для построения этой сетки используется классический метод эквираспре-деления (уравнение (52) с в — 0).
Положение граничных узлов также меняется при переходе от одного временноого шага к другому, поэтому перед тем, как решать уравнения (53), необходимо построить на границе области новую сетку хп+1. Пусть некоторый участок границы Г задан в параметрической форме ' — '(р), У — у(р), где р — параметр. Тогда координаты узлов хп+1 на границе определяются по формулам 'п+1 — '(рп+1), Уп+1 — у(рп+1 ), в которых
п+1
для поиска значений р^ используется нелинейное разностное уравнение
1 ( рп+1 _ рп+1 рп+1 _ рп+1\ рп+1 _ рп
К I (^п3п+1)^+1/2р^+1 К - (^п3п+1)^-1/2^ ДМ — в
с функцией 3 — (('р)2 + (уР)2)1/2 и управляющей функцией ад, вычисленной в граничных узлах по известному на старой сетке решению
Как было указано в [17], при расчёте разрывных решений на подвижных сетках возникают проблемы, связанные с разрешимостью уравнений для сетки и с адаптацией сетки к разрывному решению. К таким проблемам относится, например, невозможность сосредоточить достаточное число узлов сетки в зонах скачкообразного изменения решения. Другой пример - возникновение осцилляций управляющей функции, приводящих к немонотонному изменению длин соседних ячеек: происходит чередование длинных и коротких ячеек сетки, в результате чего понижается точность расчёта. Для устранения указанных проблем может быть использовано сглаживание управляющей функции ад, при котором в область разрыва попадает боольшее количество узлов и происходит плавное изменение длин соседних ячеек, что позволяет лучше воспроизводить численное решение.
Здесь использована неявная процедура сглаживания управляющей функции [19] на основе неявной схемы для двумерного уравнения теплопроводности, которая для узла, обозначенного на рис. 1 цифрой 0, выглядит следующим образом:
ад0 - ад0 ад3 - 2ад0 + ад4 - 2ад0 + ад2
— а К2 +а2 К2 '
где w — управляющая функция, зависящая от решения задачи, ад — сглаженная управляющая функция,
т — + К2) а а-1 2
т — 4 ' а — К\ + Л|' а '2'
о > 0 — параметр сглаживания. Для решения этого уравнения применяется метод расщепления по направлениям, каждый этап которого реализуется скалярными прогонками.
8. Результаты решения тестовых задач
В настоящем разделе приводятся некоторые результаты решения тестовых задач, являющихся двумерными аналогами задач из работы [17]. Область решения П изображена на рис. 2, а. Ее граница Г состоит из прямолинейных отрезков
Гл Г
пр
- 1 < х < -0.5, у = 0} х = 0, 0.5 < у < 1}
и криволинейных участков, заданных в параметрической форме:
Гн Гв
{(х,у) {(х,у)
х х
— сое р, у = — втр} , —0.5 сое р, у = —0.5вт р}
где параметр р пробегает значения из отрезка [0, 3п/2]. Эти части границы являются дугами окружностей с центром в начале координат и радиусами 1 и 0. 5.
Компоненты вектора скорости не зависят от времени и задаются по формулам
и = —у, V = х.
(54)
Очевидно, что уравнение неразрывности (2) выполняется, а на границе Г справедливы соотношения (4). Для заданного вектора скорости (54) задача (29)-(31) имеет точное решение
Ф = —
х2 + у2 2 ;
(55)
поэтому указанная задача относительно функции ф численно не решается и в формулах (33) для контравариантных компонент скорости используются точные значения (55) в узлах криволинейной сетки хп.
У
1.0 -I 0.5 0.0 -0.5 -| -1.0
и 1.0
0.5
0.0 -1.0
-1.0 -0.5 0.0 0.5 1.0
1.0
0.5 У
1.0
-1.0
а
Рис. 2. Сетка, построенная с использованием управляющей функции (58). N1 = 100, N2 = 20, а = 0 (а); график начальной функции (56) (б)
8.1. Задача с гладким решением
Начальная функция (6) задаётся в следующем виде:
Po(x) = e- (a |x - x0|)2, (56)
где a = 5, x0 = (—0.75л/2/2; —0.75л/2/2)• График начальной функции показан на рис. 2, б. Точное решение задачи (1), (5), (6), (54), (56) описывается формулой
<^(x,y,t) = <£0(x cos t + y sin t, y cos t — x sin t), (57)
при этом на входе в область полагается ялев (x, t) I „ = ю (x, t) I „ . Графиком
I xfci лев I xfci лев
решения является экспоненциальная "шапочка" единичной высоты, которая без изменения формы и с постоянной угловой скоростью вращается как твердое тело вокруг начала координат.
На рис. 3, а изображено точное решение (57) в момент времени t = 2.87, на рис. 3, б-г приведены графики численного решения в тот же момент времени на "квазиравномер-
Рис. 3. Графики точного решения (а) и численного решения, полученного по схемам предиктор-корректор (б), Лакса — Вендроффа (в) и противопоточной (г)
У 1.00.50.0-0.5-1.0-
-1.0 -0.5 0.0 0.5
1.0
X
и
1.0 У
б
1.0
0.5 У
1.0
-1.0
Рис. 4. Адаптивная сетка в момент времени Ь = п (а) и график численного решения, полученного по схеме предиктор-корректор на динамически адаптивной сетке (б). N1 = 100, N2 = 20, а = 50
а
ной" сетке (см. рис. 2, а), полученной при использовании управляющей функции
'1+1/2 = 1 + а|^"+1/2! (58)
с параметром а = 0. Эта сетка адаптируется лишь к криволинейной границе области и не учитывает поведение решения. Видно, что амплитуда решения, полученного по схеме предиктор-корректор, стала меньше амплитуды точного решения и решения по схеме Лакса — Вендроффа. Вместе с тем использование схемы предиктор-корректор позволяет избежать появления численных осцилляций, возникающих в схеме Лакса — Вендроффа, и сильного размазывания решения, происходящего при применении про-тивопоточной схемы.
При использовании адаптивной сетки обеспечивается сильное сгущение узлов под "шапочкой" (рис. 4, а), при этом область максимальной концентрации узлов вращается вместе с "шапочкой" против часовой стрелки. Из рисунка 4, б видно, что форма и амплитуда решения на подвижной адаптивной сетке сохраняются с течением времени значительно лучше, чем на неподвижной (сравни с рис. 3, б).
8.2. Задача с разрывным решением
Рассмотрим задачу о движении ступеньки в области П. Начальная функция (6) задаётся здесь в следующем виде:
(Х) = Ыр) = { 1 при Р - Р0' (59)
7 ^ 0 при р>р0,
где ро = п/4,
п — агссов(ж/уж2 + у2) при у — 0,
п + агссо8(ж/\/ж2 + у2) при у> 0.
-1.0
1.0 1.0 а
1.0
-1.Г
-1.0
0.5 \ ^ 0.5
1.0 1.0
б
1.0 1.0
Рис. 5. 1 — Графики численного решения, полученного по схеме предиктор-корректор (а — в); 2 — графики точного решения (а) и численного решения, полученного по схеме Лак-са — Вендроффа (б) и по противопоточной схеме (в)
и
0.5
-1.0
в
Точное решение задачи (1), (5), (6), (54), (59) описывается формулой
^(ж,y,t) — ^0(р - ^ (60)
и является в каждый момент времени кусочно-постоянной функцией, прямолинейная линия -' 8т(р0 + ¿) + у ео8(р0 + ¿) — 0 разрыва которой вращается с единичной угловой скоростью вокруг начала координат.
На рис. 5 представлено сравнение результатов расчётов по трём схемам. Схема предиктор-корректор для этой задачи обеспечивает монотонное решение и приемлемое качество описания движения фронта разрыва. Схема Лакса — Вендроффа, как и ожидалось, дает нефизичные осцилляции, а противопоточная схема сильнее размазывает разрыв в сравнении со схемой предиктор-корректор.
Заключение
В данной работе схема предиктор-корректор, исследованная в [17] для одномерного линейного уравнения переноса с постоянным коэффициентом, обобщена для двумерного линейного уравнения переноса с переменными коэффициентами. Описан способ аппроксимации контравариантных компонент скорости, гарантирующий выполнение уравнения неразрывности для сеточных функций на подвижных криволинейных сетках. Указан выбор схемных параметров, при котором сохраняется монотонность численного решения. Доказано выполнение разностного аналога геометрического закона сохранения, что гарантирует сохранение схемой предиктор-корректор постоянной функции. Предложена модификация классического метода эквираспределения построения подвижных сеток, позволяющая избежать возникновения осцилляций траекторий узлов и резкого изменения площадей соседних ячеек сетки. Указан метод построения сетки на границе подвижной области. Для решения проблем, связанных с разрешимостью уравнений для сетки и с адаптацией сетки к разрывному решению, предложено использование неявной процедуры двумерного сглаживания управляющей функции. Приведено сравнение результатов решения тестовых задач с гладким и разрывным решениями на неподвижных и динамически адаптивных сетках.
Результаты, полученные для двумерного линейного уравнения переноса с переменными коэффициентами, имеют общий характер и могут быть использованы для численного решения задач наводнения-осушения в рамках модели мелкой воды.
Список литературы
[1] Bradford S.F., Sanders B.F. Finite-volume model for shallow-water flooding of arbitrary topography //J. Hydraul. Eng. 2002. Vol. 128. P. 289-298.
[2] Bokhove O. Flooding and drying in discontinuous Galerkin finite-element discretizations of shallow-water equations. Part 1: One dimension //J. Sci. Comput. 2005. Vol. 22-23. P. 47-82.
[3] Westerink J.J., Feyen J.C., Atkinson J.H. et al. A basin to channel scale unstructured grid hurricane storm surge model applied to Southern Louisiana // Monthly Weather Rev. 2008. Vol. 136, No. 3. P. 833-864.
[4] Ern A., Piperno S., Djadel K. A well-balanced Runge-Kutta discontinuous Galerkin method for the shallow-water equations with flooding and drying // Intern. J. Numer. Meth. Fluids. 2008. Vol. 58. P. 1-25.
[5] Bunya S., Kubatko E.J., Westerink J.J., Dawson C. A wetting and drying treatment for the Runge-Kutta discontinuous Galerkin solution to the shallow water equations // Comp. Meth. Appl. Mech. Eng. 2009. Vol. 198, No. 17-20. P. 1548-1562.
[6] KArnA T., de Brye B., Gourgue O. et al. A fully implicit wetting-drying method for DG-FEM shallow water models, with an application to the Scheldt Estuary // Ibid. 2011. Vol. 200, No. 5-8. P. 509-524.
[7] ВольцингЕР Н.Е., Клеванный К.А., Пелиновский Е.Н. Длинноволновая динамика прибрежной зоны. Л.: Гидрометеоиздат, 1989.
[8] LoNG-Wave Runup Models / Eds. H. Yeh, P. Liu and C. Synolakis. Singapore: World Sci., 1996.
[9] Федотова З.И. Обоснование численного метода для моделирования наката длинных волн на берег // Вычисл. технологии. 2002. Т. 7, № 5. С. 58-76.
[10] Федотова 3.И. О применении разностной схемы Мак-Кормака для задач длинноволновой гидродинамики // Там же. 2006. Т. 11. Спец. выпуск, посвящённый 85-летию со дня рождения академика Н.Н. Яненко. Ч. 2. С. 53-63.
[11] CuSYAKOv V.K., Fedotova Z.I., Khakimzyanov G.S. et al. Some approaches to local modelling of tsunami wave runup on a coast // Rus. J. Numer. Anal. Math. Model. 2008. Vol. 23, No. 6. P. 551-565.
[12] Synolakis C.E., Bernard E.N., Titov V.V. et al. Validation and verification of tsunami numerical models // Pure and Appl. Geophisics. 2008. Vol. 165, No. 11-12. P. 2197-2228.
[13] ХАкимзянов Г.С., ШокинА Н.Ю. Определение зоны затопления при накате длинных волн на берег // Труды шестого Совещания российско-казахстанской рабочей группы по вычислительным и информационным технологиям. Алматы: Казак университета, 2009. С. 305-314.
[14] Arney S.D., Flaherty J.E. A two-dimensional mesh moving technique for time-dependent partial differential equations // J. Comput. Phys. 1986. Vol. 67, No. 1. P. 124-144.
[15] Bautin S.P., Deryabin S.L., Sommer A.F. et al. Use of analytic solutions in the statement of difference boundary conditions on a movable shoreline // Rus. J. Numer. Anal. Math. Model. 2011. Vol. 26, No. 4. P. 353-377.
[16] Александрикова Т.А., Галанин М.П., Еленина Т.Г. Нелинейная монотонизация схемы К.И. Бабенко для численного решения уравнения переноса // Матем. моделирование. 2004. Т. 16, № 6. С. 44-47.
[17] ХАкимзянов Г.С., ШокинА Н.Ю. Некоторые замечания о схемах, сохраняющих монотонность численного решения // Вычисл. технологии. 2012. Т. 17, № 2. С. 78-98.
[18] ОлЕЙник О.А., Радкевич Е.Б. Уравнения второго порядка с неотрицательной характеристической формой. Математический анализ, 1969 (Итоги науки). М.: ВИНИТИ, 1971.
[19] ХАкимзянов Г.С., Шокин Ю.И., Барахнин В.Б., ШокинА Н.Ю. Численное моделирование течений жидкости с поверхностными волнами. Новосибирск: Изд-во СО РАН, 2001.
[20] Trulio J.G., Trigger K.R. Numerical Solution of the One-Dimensional Hydrodynamic Equations in an Arbitrary Time-Dependent Coordinate System. Tech. Rep. UCLR-6522. Univ. of California, Lawrence Radiation Laboratory, 1961.
[21] Thomas P.D., Lombart C.K. Geometric conservation law and its application to flow computations on moving grids // AIAA J. 1979. Vol. 17, No. 10. P. 1030-1037.
[22] Etienne S., Garon A., Pelletier D. Perspective on the geometric conservation law and finite element methods for ALE simulations of incompressible flow //J. Comput. Phys. 2009. Vol. 228. P. 2313-2333.
[23] Baines M.J., Hubbard M.E., Jimack P.K. Velocity-based moving mesh methods for nonlinear partial differential equations // Commun. Comput. Phys. 2011. Vol. 10, No. 3. P. 509-576.
[24] ХАкимзянов Г.С., ШокинА Н.Ю. Метод эквираспределения для построения адаптивных сеток // Вычисл. технологии. 1998. Т. 3, № 6. С. 63-81.
Поступила в 'редакцию 29 июня 2012 г.
