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

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

Том 2, № 5, 1997

О НЕКОТОРЫХ СПОСОБАХ АППРОКСИМАЦИИ УРАВНЕНИЯ ДЛЯ ПОТЕНЦИАЛА *

Г. С. Хлкимзянов, Э. В. ЧУБАРОВА Институт вычислительных технологий СО РАН Новосибирск, Россия e-mail: khak@adm.ict.nsc.ru

We present some approximations of mixed problem in a complex domain for the equation with respect to potential of velocity vector and Dirichlet or Neuman conditions on different parts of boundary based on integro-interpolational approach. Coefficients of 9-point difference equations are also presented. Several properties of these equations are shown, the results of test computations are discussed.

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

= 0, (1.1)

решение которого ищется в криволинейной области П плоскости х1Ох2.

При использовании адаптивных сеток подвижная область П(£), занятая жидкостью, отображается в каждый момент времени £ на неподвижную расчетную область Q с помощью преобразования

ха = ха(дв,£), а, в = 1,2. (1.2)

Здесь для простоты изложения в качестве Q рассматривается единичный квадрат в плоскости д1Од2. В новой координатной системе д1,д2,1 уравнение (1.1) имеет вид

д (h д^ + h 9tp\ + д (h д^ + h 9tp \ -п (1 3)

dq1 Г11 д«1 + hl2W2) + W2 V s«1 + h22 5?) -n' (1'3)

где к11 = д22/3, к12 = к21 = -д12/3, к22 = д11/3, да@ — ковариантные компоненты метрического тензора, 3 — якобиан преобразования (1.2).

Область Q покрывается равномерной прямоугольной сеткой с шагом ка и количеством узлов Ыа в направлении оси Ода. Совокупность внутренних узлов ^ сетки (г — мультииндекс, г = (г1 ,г2)) будем обозначать Qh, граничных — 7^. Будем считать, что в узлах этой сетки определены сеточные функции (р и х“, в то время как коэффициенты кав могут быть определены как в узлах с целыми индексами, так и с полуцелыми. На сетке Qh уравнение

* Работа выполнена при финансовой поддержке Российского фонда фундаментальных исследований, грант №97-01-00819.

© Г. С. Хакимзянов, Э.В. Чубарова, 1997.

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

Согласно алгоритму решения задачи о потенциальных течениях [1], после определения значений потенциала на свободной границе Гf (рис. 1) вычисляются значения р внутри области О, на нижней и боковых ее границах. В работе предполагается, что на границах Го и Гь задано условие непротекания

тг = »■

дп

(1.4)

которое в координатах да принимает вид

к дг + к дГ

МІТТТ + к12ТГ^ дд1 од2

70

П к дг + к дг

П> к21ТТТ + к2 2ТГ2 дд1 од2

П.

(1.5)

1ъ

На границе 7/ значение р считается известным:

гк = р(дг)> П < д1 < 1

2. Получение разностных уравнений основывается на интегроинтерполяционном методе.

Рис. 1. Физическая П и вычислительная Q области. Для этого уравнение (1.3) записывается в виде контурного интеграла

\ дГ + к дГ \ л~2

к11ТТТ + к12ТТ^ дд1 дд2

)йд2 - (к21 Цр + к22Цр) ¿д1 = 0

(2.1)

где С — произвольный контур, гомеоморфный окружности.

В качестве С возьмем прямоугольник ЛЕСО (рис. 2, а), стороны которого параллельны координатным осям и делят расстояния до соседних узлов пополам. Для этого контура интегральное соотношение (2.1) примет вид

[ (к дг + к дг ' ^

і і дд1 +к12 дд2

(ВС) (АВ)

Уд2 - / (к11 ^ + к12Ц) ід2+

+ / (к'21 0+к22 Шід1 - / (к21 ^+к22 Шід1

(ВС) (АВ)

(2.2)

Рис. 2. Шаблон для внутреннего (а)и граничного Рис. 3. Нумерация типов узлов

(б) узлов. в вычислительной области.

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

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

/ (к" ъд +к12)Лд2 ~

(ВС)

н (к11(П) + к11(3) р3 - р0 , 1 (к (3) р7 - р6 , к (П) р4 - р2

- л, (------2-----------нт + 2 + МП)_ЖГ.

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

(3.1)

1

а=1

1^2

+ 4 [(ка/3рда )дв + (кав)дв + (ка/3рда )дв + (ка/3)дв] = П.

Здесь и / означают [2] разности вперед и назад соответственно.

(3.2)

Разностные уравнения (3.2) выписываются во всех узлах, не принадлежащих верхней границе 'у/ расчетной области. Для внутреннего узла используется девятиточечный шаблон (см. рис. 2), для граничных узлов шаблон состоит из меньшего количества точек. На рис. 2, б показан шаблон и контур интегрирования для узла дг, расположенного на нижней границе области Q. В силу граничного условия (1.5) интеграл по стороне ЛЕ в равенстве

2

(2.2) равен нулю, интеграл по ПС аппроксимируется аналогично случаю внутреннего узла, а для интеграла по ВС вместо (3.1) получим

(вс)

кц + кі21^2 )

од1 од2

к2 Г кіі(0) + кіі(3) ^3 — ^0 + 1

ТІ 2 к- 2

к12 (3)

к

^3 ,7 /п\^4

------+ к-2(0) —

^0

к2

Полученные разностные уравнения можно записать в виде

8

^=0

(3.3)

(3.4)

В табл. 1 приведены выражения для коэффициентов в зависимости от типа узла, определяемого его расположением на сетке. Внутренним узлам приписан тип 0, узлам на левой границе — тип 1 и т. д. (рис. 3). На верхней границе расчетных узлов нет, так как там значение функции известно.

Таблица 1

Тип узла аі а2 аз а4 а5 ав а7 а8

0 ві в2 вз в4 в5 вб в7 в8

1 0 % — 72 вз Ах 2 — 74 0 вв в7 0

2 1Г — 7і 0 # — 73 в4 0 0 в7 в8

3 ві + 72 0 в4 + 74 в5 0 0 в8

5 0 0 — 73 <2 — 74 0 0 в7 0

6 § — 7і 0 0 «і + 74 0 0 0 в8

В таблице использованы следующие обозначения:

Ь Ь

в3 = (М°) + ад)), 3 = 1, 3; в3 = 2ьЬ2(к22(0) + k22(j)), 3 = 2, 4;

в1 — 4(к-2(І — 4) + кі2(і - 3)), 3 — 5, 7; Рі — — 4(к-2(і - 4) + кі2(9 — ^)), і — 6,8;

13 = 4(к12(0) - к12(3)), 3 = 1,2; 73 = 4(к12(3) - к12(0)), 3 = з, 4.

Коэффициент ао вычисляется по формуле

8

ао = — ^ ' а^. (3.5)

3=1

Коэффициенты уравнения (3.4) удовлетворяют равенствам

а3(^гх,£2) а1(^гх + 1,£2), а4(^гх,£2) а2(Qiх,І2 + l'),

а7(^гх,г2) а5 (^гх + 1,г2 + 1), а8(0_11,12) а6 (^гх-1,£2+1) , (3.6)

что означает симметричность матрицы системы уравнений (3.4).

Разностное уравнение (3.2) аппроксимирует дифференциальное уравнение (1.3) на его гладких решениях р со вторым порядком. Так как линейная функция р(х,у) = ах + Ьу + с удовлетворяет уравнению Лапласа (1.1), было бы желательно, чтобы она являлась решением и сеточных уравнений (3.2). Справедливы следующие утверждения.

Утверждение 1. Функции р = х“(^1, д2), а =1, 2, д = (д1, д2) € < удовлетворяют уравнению (1.3).

Утверждение 2. Пусть рг = ха(дг), дг € где ха — функции отображения (1.2). Тогда

А(1)рг = 0, д € Ян. (3.7)

4. При аппроксимации уравнения (1.3) вторым способом коэффициенты кар вычисляются на сторонах ячеек: к11 и к12 — в точках дгх+1/2,г2, к21 и к22 — в точках дгх,г2+1/2. Производные х^з аппроксимируются центральными разностями Одвха. Формулы, по которым вычисляются разностные производные, зависят от того, в какой точке (в узле, в середине стороны или в центре ячейки) требуется определить производную. Так как функции р и ха определены только в узлах сетки, то для вычисления Одв используются при необходимости усреднения. Например,

Одхр(Е) Ь , Одхр(С) 2 [Одхр(дгх+1/2,г2) + Одхр(дгх+1/2,г2+1)] ,

П,хр(0) = - [П,хр(Е) + П,хр(Ж)] , П,хр(М) = - [Одхр(С) + П,хр(О)]. (4.1)

2 [Пд. р(Е) + О,, р(И-’)], О,, р(М) = 1

Теперь вместо (3.1) будем иметь следующее выражение для интеграла по стороне ВС:

к11 Ър + к12 ^д2 ~ Ь2 [к11(Е)Одх р(Е) + к12 (Е)Од2 р(Е)]

(ВС)

а разностное уравнение во внутренних узлах сетки примет вид

2

Л(2)рг = ^ О,а (ка/3Од13р) (дг) = 0. (4.2)

а,в=1

Справедливо

Утверждение 3. Пусть рг = ха(дг), где ха — функции отображения (1.2) . Тогда

Л(2)рг = 0, дг € Яъ.

Доказательство следует из тождеств

П (1 П 17 П Л / ч у5 — у6 + у7 — у8

Одх (к11Одхх + к12Од2х) (дг) т ,

4Ъ1Ъ2

П П П 17 П ( \ у5 — у6 + у7 — у8

Од2 (к21Одхх + к22Од2х) (дг) =-------77-т-------,

4Ь1Ь2

и аналогичных тождеств для функции р = у.

Таким образом, разностное уравнение (4.2) обладает тем же свойством, что и дифференциальное уравнение (1.3): на линейных функциях оно выполняется тождественно.

Коэффициенты а^ для разностного уравнения (4.2), записанные в виде (3.4), приведены в табл. 2. Однако для этого метода матрица системы (3.4) не является симметричной, равенства (3.6) в этом случае не выполняются.

Таблица 2

Тип узла «1 а2 аз а4 «5 «6 «7 «8

0 в1 - 74+ +72 в2 - 73 + +71 вз + 74-72 в4 + 7з--71 71 + 72 -72 - 7з 7з + 74 1 4 1 1

1 0 в2/2- -7з + 72 вз + 74-72 в4/2 + +7з - 74 0 -72 - 7з 7з + 74 0

2 вх/2- -74 + 71 0 вз/2+ +74 - 7з в4 + 7з--71 0 0 7з + 74 -74 - 71

3 в1 - 74 + +72 в2/2+ +71 - 72 0 в4/2+ +74 - 71 71 + 72 0 0 -74 - 71

5 0 0 вз/2+ +74 - 7з в4/2 + +7з - 74 0 0 7з + 74 0

6 в1/2- -74 + 71 0 0 в4/2 + +74 - 71 0 0 0 1 ? 1

Здесь

А = Ъ2 кц(Ш), в2 = ^ М£), вз = ^ кп(Е), в4 = ^ к22(н),

Ъ П2 Ъ П2

71 = ^ к12(Ж), 72 = 4 к21(Б), 73 = 4 кл2 (Е), 74 = 1 к2л(^).

5. При третьем способе аппроксимации уравнения (1.3) коэффициенты ка@ вычисляются в центрах ячеек сетки — в точках дг1+1/2,г2+1/2. Для интегралов в (2.2) по сторонам прямоугольника ЛБСО применяется формула трапеций. Тогда вместо (3.1) для аппроксимации интеграла по БС имеем следующее выражение:

кц + к121 ^

\ од1 од2

(ВС)

Ъ

- -2 [ки(Б)Дд¡¿>(Б) + к12(Б)Д^р(Б) + кп(С)Ддр(С) + МС)^р(С)] (5.1)

и разностное уравнение во внутреннем узле д. примет вид

2

Л(3) ^ Ода (кавОдв <р) (д.) = 0, (5.2)

а,в=1

где

Утверждение 4. Пусть <^г = ха(дг), где ха — функции отображения (1.2) . Тогда

Л(з)^ = 0, д € Ян.

В табл. 3 приведены коэффициенты разностного уравнения (5.2), записанного в виде

(3.4). Матрица системы уравнений (3.4) симметрична. Так же как и для первой аппроксимации, разностный оператор получается самосопряженным и положительно определенным. Интересно отметить, что в декартовых координатах на квадратной сетке эта схема переходит в известную пятиточечную схему “косой крест"[3].

Таблица 3

Тип узла «1 а- аз а4 а5 ав аг а8

0 в1 + ^4 - в1 - в2 в2 + вз - вз - в4 в5 вб в7 в8

1 0 -в2 в2 + вз -вз 0 вв вг 0

2 в4 0 вз - вз - в4 0 0 вг в8

3 в1 + в4 -в1 0 - в4 в5 0 0 в8

5 0 0 вз -вз 0 0 вг 0

6 в4 0 0 - в4 0 0 0 в8

Здесь

А = 4^к11<А>- -Ъ-МА), А = ^МВ)- -Ъ-*'22<В). А = к“(с) - 4Й2• А = 4^ М^) - -^12

Ъ Ъ 1 Ъ Ъ 1

в5 = ^кц(А) + ^к--(А) + 2МА вб = ^ки(В) + ^МВ) - 2к!2(В),

Ъ Ъ 1 Ъ Ъ 1

вг = ^кц(С) + ^МС) + 2ЫС), в8 = ^кц(Б) + 4ъ-М^) - 2^Р).

6. Приведем результаты численных расчетов при использовании различных аппроксимаций уравнения (1.3). В качестве области П в тестовой задаче брался криволинейный

четырехугольник, ограниченный слева и справа прямыми х = 0 и х = 1, снизу — прямой у = -1, сверху — кривой

у = п(х) = -а сов(2пх). (6.1)

В качестве точного решения задачи (1.1), (1.4) была выбрана функция

рвх(х, у) = — * ' совЬ[2п(у + 1)] ■ сов(2пх). (6.2)

С08П(2п)

Для отображения области П на единичный квадрат Я использовалось преобразование

х = д1, у = -1 + д2(п(х) + 1). (6.3)

Рис. 4. Отклонение численного решения от точного для методов

1(и), 2 (□), 3(Д).

11 41 161

N

Рис. 5. Оптимальные значения параметра релаксации для методов

!(■), 2 (□), 3(Д).

Система разностных уравнений, полученная с помощью каждого из рассмотренных здесь способов аппроксимации, решалась методом последовательной верхней релаксации. Расчеты показали, что все три схемы, несмотря на несимметричность матрицы для метода 2 и ненулевую невязку в (3.7) для линейной функции в первом методе, имеют порядок точности, близкий ко второму. В табл. 4 приведены оптимальные значения параметра релаксации т, количество необходимых для сходимости итераций к1а8і и норма разности между точным и численным решениями ||р — рех||с на различных сетках при первом способе аппроксимации и а = 0.5. На рис. 4 представлена зависимость величины ||р — рех||с от числа узлов N = N = квадратной сетки для всех трех разностных схем. Можно заметить, что схема 3 дает несколько более высокую точность на мелких сетках. На рис. 5 показано влияние размера сетки Qh на оптимальные значения параметра в итерационном методе последовательной верхней релаксации.

Таблица 4

N1 N2 топт ||^1а81 — ^ех ||с

11 6 1.24 4 0.32 • 10-1

11 11 1.32 4 0.32 • 10-1

21 11 1.48 11 0.13 • 10-1

21 21 1.58 19 0.90 • 10-2

31 31 1.68 39 0.45 • 10-2

41 21 1.72 35 0.36 • 10-2

41 41 1.78 64 0.26 • 10-2

61 61 1.84 133 0.12 • 10-2

81 41 1.82 148 0.96 • 10-3

81 81 1.94 283 0.66 • 10-3

161 81 1.96 442 0.24 • 10-3

161 161 1.96 896 0.17 • 10-3

Таблица 5

N1 N2 г(1) ю г(3)

11 11 0.31 ■ 10-1 0.95 ■ 10-15 0.35 ■ 10-15

21 21 0.22 ■ 10-2 0.11 ■ 10-14 0.55 ■ 10-15

41 41 0.15 ■ 10-3 0.12 ■ 10-14 0.68 ■ 10-15

В табл. 5 приведены значения невязки г(а) = Л(а)р, а = 1, 2, 3, полученные для функции р = у при использовании рассмотренных способов аппроксимации. Эти результаты согласуются с утверждениями 2-4.

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

[1] Шокин Ю. И., Рузиев Р. А., ХАкимзянов Г. С. Численное моделирование плоских потенциальных течений жидкости с поверхностными волнами. Препринт ВЦ СО АН СССР, №12, Красноярск, 1990.

[2] Самарский А. А. Теория разностных схем. Наука, М., 1983.

[3] Самарский А. А., Андреев В. Б. Разностные методы для эллиптических уравнений. Наука, М., 1976.

Поступила в редакцию 14 марта 1997 г., в переработанном виде 31 июля 1997 г.