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

УДК 519.635.1

О ЧИСЛЕННОМ ИНТЕГРИРОВАНИИ НЕСТАЦИОНАРНОГО НЕОДНОРОДНОГО БИГАРМОНИЧЕСКОГО УРАВНЕНИЯ В ЗАДАЧАХ КОНДУКТИВНОЙ СВОБОДНОЙ

КОНВЕКЦИИ

В. И. Ряжских, М. И. Попов

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

Ключевые слова: бигармоническое уравнение, кондуктивно-ламинарная свободная конвекция, конечно-разностная схема

Введение.

Режим "ползущего" течения наиболее часто используется для описания поведения слабодефор-мируемых жидких сред и, в частности, ньютоновских жидкостей с помощью уравнений Стокса [1], в том числе, в задачах кондуктивно-ламинарной свободной конвекции [2]:

dv 1 , _ _

— = — grad p + vV v -ptg, дг p

div v = 0,

* = «V ^t, дг

(1) (2) (3)

где v, p, t - вектор скорости, давление, температура; p, v, a, p - плотность, коэффициенты кинематической вязкости, температуропроводности и объемного расширения среды; г - текущее время, g - вектор ускорения силы тяжести. Для исключения давления в системе (1)-(3), ее записывают в переменных Гельмгольца [3]:

дю _ , _

— = vV ю-pgradt х g, дг

V2y = -ю,

* = aV 2t,

дг

где ю = rot v, v = rot y - вихрь и функция тока, или

исключением ю приводят к бигармоническому виду:

d(V2y) _

дг

- = vV y — p gradt х g,

^ = «V2t.

дг

(4)

(5)

В стационарном случае система (4)-(5) может даже допускать аналитическое решение с соответствующими граничными условиями, например [4,5] или численно проинтегрирована по общеизвестной

Ряжских Виктор Иванович - ВГТУ, д-р техн. наук,

профессор, тел. (473) 246-42-22

Попов Михаил Иванович - ВГУИТ, аспирант, тел.

8-904-210-19-98

неявной [6], или в рамках метода установления, явной схемами [7]. В общей постановке конструктивные методы решения системы (4)-(5) пока отсутствуют.

В данной работе предлагается восполнить этот пробел путем разработки эффективной конечно-разностной схемы для численного интегрирования нестационарного квазибигармонического уравнения для квадратной области с постановкой граничных тепловых условий, соответствующих первой тестовой задаче [8].

Постановка задачи и построение численной схемы.

В безразмерном виде система (4)-(5) декомпозируется на две последовательно решаемые задачи:

тепловую

дТ (X, 9) 1 д2Т (X, 9)

----1—, (6)

д9 Рг дХ2

Т (X ,0) = 0, Т (0,9) = 1, Т (1,9) = 0; (7)

и гидродинамическую

д (д2Ф(Х,У,9) д2Ф(X,У,9)| = д4Ф(X,У,9) 59^ дX2 + 5У2 ) дX4 +

д4Ф( X ,У, 9) д4Ф( X ,У, 9) дТ (X, 9)

+2-;-;--1-----, (8)

дX дУ2 дУ4 дX

Ф(X ,У ,0) = 0, (9)

Ф(0, У, 9) = Ф(1, У, 9) = Ф^ ,0,9) =

= Ф^ ,1,9) = 0, (10)

дФ(0, У, 9) _ дФ(1, У, 9) _ дФ(X,0,9)

дХ

дХ дУ

дФ( X ,1,9)

дУ

= 0,

(11)

где Т (X, 9) - температура, Ф^, У, 9) - функция тока, Рг - число Прандтля, 9, X, У - текущие время и декартовы координаты, заданную в замкнутой области Б = [0,1]х [0,1] функций.

Решение тепловой задачи (6)-(7) известно [9]

2 ¥ (_1) Р

Т (X, 9) = 1 _ X + - У • р р=1 Р

2 2

• sin[(1 -X)pp]exp| -^p-e

(12)

и поэтому сосредоточим основное внимание на системе (8)-(11), имея ввиду (12).

Для построения конечно-разностной схемы временной отрезок разделен на равные промежутки De . В каждый момент времени kДе, k = 1,... область решения D системы (8)—(11) заменена сеткой с шагами ДХ , DY , wDXDr = {(X,.,Y.) = (iDX, jDY)}, i = 0, . ., n, j = 0, .., m . Вместо функции непрерывного

аргумента на D рассматривается функция дискретного аргумента Ф(Xt, Yj, ek)ДХ DY Де, которая обозначается Fkj. Выбор разбиения осуществлен таким

образом, чтобы в него попадали точки границы области и центр области.

Конечно-разностная аппроксимация бигармо-нического оператора ВД в каждый момент времени представлена в [7].

Разностный аналог вторых производных

д2Ф( X, Y, е) д2Ф( X, Y, е)

-z- и --- записан через цен-

dX dY

трально-разностный оператор второго порядка по соответствующим переменным.

d 2Ф( X, Y,

dX2

д 2Ф( X ,Y,

Ф

i + 1, j

- 2Ф* +Фк-1

i, j,k

(ДГ)2

- +

O [(AX )2 ],

Фк - 2Фк + Фк ,

C _ _,+!. j i, j i-!. j

^—' A Y

dY2

(ДГ)2

Фk . , - 2Фк . + Фк .

4 j+1

i, j,k

C _

ДГ

(^)2

i + O [(ДY )2 ],

Ф- j+1 - 2ф, j +Ф* j-.

(Д7)2

Обозначим СД = -(СДХ + СД7). Конечно-разностная

„ д

аппроксимация производной — представлена опе-

дб

фк+1 - фк

Л 1,1 1,1 ратором первого порядка А = ——-.

Записывая уравнение (8) в дискретном виде и отбрасывая слагаемые более высокого порядка малости, получена конечно-разностная схема

Де

Ф, -2Ф + Ф , Ф ,-2Ф +Ф ,

l+1,j l,j ,-1,j + i,j i, j-1

(ДТ)2

(^)2

^ Ф - 2Ф+Ф Ф - 2Ф+Ф ^

^i+1,j ^ i-1,j + Vi,j+1 ^^i,^i,j-1

(^ ) 2

(AY)2

Ф

i+2, j

- 4Ф+1, + 6Ф- 4фм,: + ф

i-2, j

Ф

(ду)4

- 4Ф,., j+1 +6Ф1, j- 4Ф1, j-1+Ф,

i, j+2

l, j-2

^Y)4

-[(Ф,ч

,2 [V-i+2,j+2 -16Ф,+2,j+1 +

144 (ДУ)2 ^Y)2 +30Ф.+2,j -16Ф

i+2, j-1 + Ф.+2,j-2 ) --16 (Ф,+1 , j+2 16Ф,+1, j+1 + 30Ф ,+1, j 16Ф ,+1, j-1 + Ф+1, j-2 ) +

+30(F,,j+2 - 16Фi,j+1 + 30Фi,j - 16Ф,,j-j +F,.,j-2) --16(F,-1,j+2 - 16Ф,-1,j+1 + 30Ф.-1,j - 16Ф..-1,j-1 +F,.-1,j-2 ) + +(F,-2,j+2 -16Ф,_2,j+1 + 30Ф,.-2,j -

-16F,.-2,j-1 + Ф/-2,j-2 ) ]k - 1 -

¥

-2X(-1)P c°s[(! - (i - l)AX) pP] •

p_1

exp | -kДе

22 p p

Pr

(13)

Начальное условие (9) на сетке таково:

ф<\ =0, 1, ] = 0п. (14)

Дискретные граничные условия, получены сужением (10) на граничные узлы сетки фк = фк = фк = фк =о

0,1 п,1 ^1,0 1,п

к = 1,2,..., 1 = 0п, у = 0п. (15) Соответствующая процедура переводит (11) в следующий вид

ф к - Фк Фк - Ф к, фк - Ф к

1,1 о, 1 _ п,1 п-1,1 _ ^1,1 ^1,0 _

AX

AX

AY

Фk -Fk

_ i,n i,n-1

AY

Учитывая (15), имеем

= 0, i, j = 1, n -1. (16)

Fk ■ = Фk 1 ■ = Фк, = Фк , = 0,

1,j n-1, j l,1 l,n-1 5

к = 1,2,..., 1 = 1,п -1, 1 = 1,п -1 (17)

Таким образом аналогом непрерывной задачи (8)—(11) является конечно-разностная схема (13)-(17).

Схема (13) является явной по времени (индекс к ) и неявной по пространственным переменным, каноническая форма которой

Сдфк+' -фк + Бдфк = фк, к = 0,1,..., ф0 =0,(18)

Д8

где фк -производная температурного поля в момент времени к Д6 .

Для дальнейшего анализа на множестве сеточных функций {и(Х(,У})ДХДГ} вводится Гильбертово

пространство Н ДХ Дг , которое для каждого разбиения представляет собой вещественное пространство (п +1) х (т +1)-мерных векторов [8], с требованием выполнения граничных условий (15)-(17) для каждой функции из НДХДГ и с заданием скалярного произведения в виде [10]

п-2 т-2

(и, у) = V, 1 ДХ Д^

1=2 1=2

и следующих норм

2

+

+

k

+

+

n-2 m-2 1

||u|| != max Iuu |, ||u|| 2 = DXDY(19)

1=2 j=2

Нормы выбраны таким образом, чтобы выполнялось условие согласования норм, при DX, AY ® 0

INI 1 переходит в обычную норму в пространстве C(D), а INI 2 в обычную норму в пространстве L 2( D) [11]. В случае отсутствия индекса у норм в выражениях можно использовать обе.

Вычисление нормы оператора перехода и оптимального шага по времени.

На равномерных сетках AX = AY = h и Ad = t введем конечномерное гильбертово пространство Hh и рассмотрим в нем операторы C = h2Ch и B = h4 Bh не зависящие явно от h.

Для дальнейших рассуждений нам потребуются следующие вспомогательные теоремы о свойствах операторов действующих в конечномерных гильбертовых пространствах[12]

Теорема 1. Если A и B самосопряженные операторы, а a, ß - вещественные числа, то оператор a A + ßB самосопряженный.

Теорема 2. Если A и B положительно-определенные операторы, то оператор c1A + c2B положительно определен при cj > 0, c2 > 0.

Теорема 3. Если оператор A положительно определен, то оператор A-1 тоже положительно определен.

Теорема 4. Если A и B положительно-определенные операторы, то все собственные значения оператора BA положительны.

Теорема 5. Если A и B перестановочные положительно-определенные операторы, то AB положительно определенный оператор.

Докажем следующее свойство оператора B.

Теорема 6. Оператор B положительно-определенный.

Доказательство. Поскольку B = Bj + 2B2 + B3 и согласно теореме 1, линейная комбинация положительно-определенных операторов с положительными коэффициентами положительно определена, достаточно показать лишь положительную определенность операторов B1,B2,B3.

Матрица Bj представляет собой блочную квазидиагональную матрицу, каждый блок которой симметричная пятидиагональная матрица ( 6 -4 1 0 - 00

Bx(n) =

-4 6 -4 1 1 -4 6 -4 0 1 -4 6

0000

6

v J

Вычисление определителя пятидиагональной матрицы

A(n) =

1 aj a2 0 - 0

a1 1 aj a2 - 0

a2 aj 1 aj - • 0

0 a2 aj 1 - • 0

0 0 0 0- • 1

сводится к решению системы разностных уравнений:

Ап _ 4,-1 + а224,-3 _ а244-4 + аА_1 _ а1а2Вп-2 = С1,

В + «2Вп_1 _ «14,-1 =

Здесь Ак = А(к), а Вк-1 = А(1,2 / к) - минор, получаемый из определителя А(к) вычеркиванием первой строки и второго столбца. Исключая Ви-1 и Ви-2 из первого уравнения получим характеристическое уравнение

г5 + (а2 -1)г4 + (а2 - а2)г3 + ^(а2 - а2)71 + а3(1 - а2)7 - а5 = 0.

При а1 = -4/ 6, а2 = 1/6 все корни равны 7к = 16, к = 1..6, поэтому

det[B,(n)] = 6nAn =

= 6й (С1 + С2 п + С3п2 + С, + С5п4) ^ =

= (С1 + С2 п + С3п2 + С4п3 + С5п4), причем постоянные С1,С2,С3,С4,С5 определяются из начальных условий

В1 (0) = 1, В1 (1) = 6, В1 (2) = 20, В1(3) = 50, В1(4) = 105. Вычислив постоянные, получим

7 23 2 1

det[Д(п)] = 1 +—п +--п2 + —п3 +--п4. Поскольку

1 3 12 3 12

det[В; (п)] > 0 "п е N , В1(п) удовлетворяет критерию Сильвестра. Матрица В3 симметрична и в другом ортогональном базисе при замене I на ] совпадает с матрицей В1 , то есть В3 и В1 подобны, поэтому В3 положительно определен. Матрица В2

представляет собой произведение положительно определенных квазидиагональной и квазипятидиа-гональной матриц. Поскольку каждая клетка квазидиагональной матрицы есть пятидиагональная матрица, а каждая клетка и квазипятидиагональной матрицы есть диагональная матрица, не трудно увидеть, что искомое произведение перестановочно. В силу перестановочности этого произведения согласно теореме 4, оно положительно определено. Таким образом, оператор В положительно определенный. ■ Теорема 7. Оператор С положительно-определенный.

Доказательство. Поскольку

С = С, + С2= (-Ф,-+1!, + 2Ф,. 1 -Ф,,, 1) +

+ (-Фу+1 + 2Фу -Фу-. ) ,

то матрица С1 является блочной квазидиагональной матрицей, каждый блок которой имеет вид

n

Г 2 -1 -1 2

0 0

0 ö 0

2

v /

Определитель этой матрицы Dn = 2Dn-1 - Dn_

По-

скольку Д = 2 и Д = 3 , то Бп >0 Vп е N . Таким образом эта матрица, а, следовательно, и матрица С1 удовлетворяет критерию Сильвестра, значит положительно определена.

Матрица С2 представляет собой блочную ква-зитрехдиагональную матрицу Г J -Е

-E J

0 ö 0

J

Ч 0 0

где J - диагональная матрица с диагональными элементами равными 2, Е - единичная матрица. В другом базисе при замене 1 на 1 эта матрица совпадает с матрицей Сх. Матрицы операторов в разных базисах подобны, поэтому матрица С2 положительно определена, следовательно, С - положительно-определенный, как сумма положительно определенных операторов. ■

Задача Коши (18) разрешима, т.к. для любого

разбиения щ , существует обратный оператор С,

существование которого следует из положительности оператора С . Поэтому схему (18) можно представить в виде

®к+1 = ^ + /, /к = Сф,

п = 0,1,... ,

ф0=0,

где

(20) (21)

£ = Е -тС"'Б - оператор перехода.

Норма оператора перехода вычисляется по формуле

М| з = тах Ы, (22)

к

где /лк - собственные значения, тогда из (21) следует

m* =1 -121 *'

(23)

где 1 к - собственные значения оператора С В , что дает возможность сформулировать задачу минимак-са [10]

INI з = min max 11 —^1* |

t>0 * ■ tf" *|. (24)

Лемма 1. Если 0<1 . <1<1 ,t>0, то

min max -1

1 -1

_ max

min max |1 -tl |=

T>0 1e[1min,Amax ] 1

x

2

:PQ>

0 1 +1 ■

max min

Из теорем 3 и 4 следует положительность собственных значений оператора C _1Б, а в силу его

конечномерности и конечность спектра, поэтому согласно лемме 1

212

Ы 3=

-1.

P0, Т =

(25)

1 +1 . 0 1 +1 .

max min max min

Таким образом получено не только значение 1Ы113, но и выражение для оптимального значения шага. Поскольку вычисление спектра оператора C -1Б довольно трудоемко, используется альтернативный способ вычисления по формулам [13]

1max=SUp(C"'Б, И); =.ШС"Б, и)

IUI 2 =1 lul 2 =

с помощью стандартных алгоритмов оптимизации. В таблице 1 приведены расчеты ||Ы|| 3 и оптимально-

го шага т„.

Таблица 1

Расчет входных параметров

Размерность сетки Мелкость разбиения INI

11x11 0.1 0.9 1.8 -10-3

17 x17 1/16 0.96 7 -10-4

19 x19 1/18 0.973 5.5 -10-4

21x 21 0.05 0.979 4.5-10-4

Анализ расчетов показал что при увеличении размерности сетки максимальное собственное значение 1max = ||C "В|| 3 стремится к значению нормы, вычисленному в сеточном аналоге нормы в C (D)

supl |Ви|| 1 ||C-В|| 1 = Ы 1_1

= 100 supl|Би|| 1 9 ,

а минимальное значение стремится к 0.

Вычисление погрешности аппроксимации уравнения разностной схемой.

На равномерной сетке вновь рассматриваем уравнение (13) в операторной форме

АС,и + Б,и = ф. (26)

Введем разность г = и - V, где и - решение задачи (18), а V - решение задачи (8)-(11) [10]. Подставляя и = 7 + V в (26), получим для г задачу во внутренних узлах сетки

АС,и + Би = АС, (г + V) + Б, (г + V) = = (АС, + Бк) г + (АС, + Б, ^ = ф, (АС, + Б,) г = -е и для г выполнены начальное (14) и граничные условия (15)-(17), где е = (АС + Б^ -ф - погрешность аппроксимации задачи (8)-(11) схемой

д 2 (13)-(17). Так как —+ Д^-ф = 0 , то

дв

e = ( ACi + Bi )v -j-

^ (-Av) + A2v -j дв

д

— (-Av) + A2v дв

= (ACh + Bh )v -Погрешность аппроксимации бигармоническо-

ull 1=1

го оператора вычислена в [7]:

| еВ |= М1к2 + М2Н4; там же вычислена погрешность аппроксимации оператора 4 :

I еа |= Мт.

Погрешность аппроксимации оператора С

, . h

I ес I £ — 12

\

Обозначив , (

M 3 = — 3 12

max

( X ,Y ,0)eDx

д V

max

(X ,Y ,0)eDx

дх4

д 4v

дх4

+ max

( X ,Y ,0)eDx

+ max

(X ,Y ,0)eDx

д 4v

dY4

д\

дY4

получим оценку погрешности аппроксимации уравнения

| е | < | еАеС | +1 еВ | < ММ3тИ + М1к2 + М2Н4. Учитывая (3.36) и обозначая

2MM3

M 0 =-3—

0 1 + 1 .

max min

+M„

получим

| e | £ M0 h4 + Mh

(27)

Ф.+J £ M

Устойчивость и сходимость.

Чтобы доказать устойчивость схемы (18) воспользуемся теоремой [10], адаптированной к рассматриваемой задаче.

Теорема 8. Для устойчивости схемы (18) достаточно, чтобы выполнялось условие для разрешающего оператора Тп. = ^Д^...^.:

Т< М " 0 < . < п. при этом верна априорная оценка

( к I

КII+У С-уц. ¡=0 0

Поскольку оператор перехода схемы (18) постоянный с нормой меньшей единицы, очевидно неравенство:

Т¡|| <1Ы1 = Р0 " 0 < . < п.

Таким образом, схема (18) устойчива и выполняется оценка

( к I ||Фк+х\\ < Р0 КИ + УТ1.

V

Далее рассматривается стационарное разностное уравнение, полученное из уравнения (13) отбрасыванием левой части равенства, в операторной форме

Ви = ф. (28)

Пусть и решение уравнения (28) и фк = ф. Поскольку схема (18) точно аппроксимирует уравнение (28) на его решении, разность гк = Фк - и удовлетворяет однородному уравнению

z — z

-L + B =0, k = 0,1,..

Решая его относительно z.

Z0=0.

получим

zk+1 = Szk. Используя рекурсию, получим zk+x = Skz0,

а поскольку ||S|| <1 для любой размерности сетки, то схема (18) сходится к решению стационарной задачи. Действительно

F — U = Fk

= SkzJ £

IM -II zjl <| |S||k||z„

® 0, k ® ¥.

1) 2) 3)

"0|| —М^И Ц^Н

Чтобы показать сходимость решения задачи Коши (18) к решению краевой задачи (8)-(11) воспользуемся теоремой [14], модифицированной под формулировку нашей задачи. Теорема 9. Если:

решение краевой задачи (8)-(11) существует в некотором классе функций, разностная задача (18) аппроксимирует краевую (8)-(11) на классе решения, разностная задача (18) корректна, то: при И ® 0 решение ФИ разностного уравнения стремится к решению Ф дифференциального уравнения, т.е. Ф-ФИН ® 0.

ii ii и

Существование и единственность решения краевой задачи (8)-(11) доказаны. Уравнение (18) аппроксимирует (8) со вторым порядком (27). Начальное (14) и граничное условия (15) аппроксимируются точно, граничное условие (17) аппроксимирует (11) с первым порядком точности. В силу определения пространства сеточных функции НИ разностная задача устойчива по граничному условию (17), поэтому корректность задачи следует из ее устойчивости по начальным условиям и устойчивости по правой части, которые были доказаны в теореме 8 Таким образом, все условия теоремы 9 выполнены, следовательно решение разностной задачи

(18) стремится к решению краевой задачи (8)-(11) при И ® 0 .

Улучшение решения на границе.

На множестве сеточных функций {и^,.,У.)И}

введем Гильбертово пространство НИ, отличающиеся от НИ отсутствием условия (17) для любой функции из Н И. Нормы на Н И введем аналогично

(19). Очевидно вложение НИ с Н И.

Рассмотрим задачу (18) в пространстве НИ и

добавим граничное условие (17). Нетрудно увидеть, что свойства операторов В и С , сохраняются и решение полученной задачи совпадает с решением задачи (18). Граничное условие (17) имеет первый порядок аппроксимации, поэтому интерполяция решения разностной задачи на границе недостаточно гладкая. Для получения гладких приближений решения краевой задачи необходима высокая степень измельчения сетки. Это требует значительных по объему вычислений. Однако, существуют и другие методы решения данной проблемы. Рассмотрим один из них.

Идея основана на пересчете приграничного слоя сеточной функции по разностной схеме аппроксимирующей (11) со вторым (или более высо-

t

ким) порядком и распределения полученной невязки на внутренние узлы сетки. Множество узлов принадлежащих границе определим как нулевой слой, приграничных узлов как первый слой и т. д. в направлении к центральной точке. Разностная аппроксимация производной на левой границе второго порядка имеет вид [15]: дФ( X ,У, 9)

дX

0,1

-3Ф 0,. + 4Ф_ 1 -Ф 2,. 2

-01-11-21 + 0(И2).

2И

Отбрасывая слагаемые более высокого порядка малости и учитывая (11) и (15), получим

Ф1,1 =4 Ф2. . (29)

Аналогично получим выражения на остальных границах.

Введем вспомогательный оператор ¥: НИ ® Нк следующим образом. Значения в приграничных узлах определяются по формуле (29), причем невязка определяется как

N = У Ф.,1 -Ф,,1,

здесь - первый слой. Распределение невязки может быть равномерным, либо пропорциональным по формуле

~ N

<Ф. . = Ф. .(1--),

,, 1 ,, Л ^

где £ - сумма значений во внутренних узлах.

Для определенности обозначим значения в узлах первого слоя до применения J а1, после применения а1; во внутренние узлах с,., среднее значение во внутренних узлах с. Докажем вспомогательную лемму.

Лемма 2.

' а, < 0, <5,. < 0, с , < 0 если а, < а, и

с < 2а1, то для К > 3

К 4К+4 К

2 >,1

4К+4

Ус,21 + У а2 >У (с. 1 - К;)2 + У а2

1,1=1 ,.=1 1,1=1 К ,.=1

где

1 К 4К+4

с =~Р£ Ус,., 1, N = У(а, - а,).

К . , 1 =1 Доказательство:

2

У с,21 + У а,2 - У (с,, 1 - Кг )2 - У а,2 = . , 1 =1 . =1 . , 1 =1 К . =1

N 4К+4

= К 2с2 - К 2(с--^)2 + У (а,2 - а,2) =

К

,.=1

N 2 4К+4

= К2с2 - К2с2 + 2сЖ- — - У а - а2) =

_ N2 4К+4

= 2cN - —2 - У (а, - а, )(а, + а,.)

К ,.=1 Поскольку <5 < 2¿г,, и ¿г,. < а,

N 2 4 К+4

2сЫ - — - У (а, - а, )(а, + а,) > К ,.=1

N2 4К+4 N2 N

2CN - - У (а,. - а,. )с = Ыс - — = N (с - —).

_ N

Так как N < 0 , оценим разность с - ,

_ _ 1_ а. - а. > а. >— с,

' ' ' 2

поэтому N > (2К + 2)с . Таким образом

_ N _ (2К + 2)с 2К + 2Ч

с--2 < с-= с (1--—).

К2 К2 К2

_ N

Для К > 3 с--2 < 0, следовательно

К

4К+4 К

2

У с,21 + У а,2 - У (с,, 1 - |т)2 - У а,2 > 0. ■

4К+4 22

¥ фи11 1 <1N11 =

,,1=1 ,=1 ,,1=1 ,=1 Теорема 10. При равномерном распределении невязки для решений Ф И выполняются оценки

¥Фи|| 2 <1 |Фи|| 2. (30)

Доказательство. Анализ расчетов показал, что во введенных выше обозначениях значения решения

Фк в узлах удовлетворяют условиям леммы 2. Левое неравенство очевидно. Для доказательства правого воспользуемся леммой 2.

N .

¥Фи|| 2 = (УУФ 2,1И2)2 = (У (с,, у + У а2)2 <

4К+4 1

,.=1 1=1

,, 1=1 4К+4 1

= I

К2

< (Ус,21 + У а,2)2=| |Фи|| 2. ■ ,, 1=1 ,=1

Из неравенства (30) и теоремы 8 следует устойчивость разностной задачи с вспомогательным оператором. Нетрудно видеть, что для данной задачи выполнены все условия теоремы 9, что означает сходимость ее решения к решению краевой задачи (8)-(11) при И ® 0.

В табл. 2 приведены результаты вычислительных экспериментов для различных сеток и значений временного шага т в момент времени 9 = 0.1 и Рг = 1.

Таблица 2

Результаты вычислительных экспериментов

Размерность сетки т | итах| Сходимость

11 х 11 2.5 •Ю-3 0.55 расходится

11 х 11 1.8 •Ю-3 2.5 •Ю-4 сходится

11 х 11 10-3 2.45 •Ю-4 сходится

17 х17 8 •Ю-4 9.3 расходится

17 х17 7 •Ю-4 2.96 •Ю-4 сходится

17 х17 5 •Ю-4 2.95 •Ю-4 сходится

21х 21 5 •Ю-4 2.48-104 расходится

21х 21 4.5 •Ю-4 3.12 •Ю-4 сходится

В граф >е | итах | максимальные по модулю зна-

чения решений, графа "Сходимость" показывает сходимость численной схемы. По таблице видно, чем мельче разбиение области, тем при меньшем отклонении от оптимального значения т0 временно-

,=1

,=1

го шага в большую сторону схема расходится. На рис. изображены профили функции тока для сетки 51х 51 при У = 0.5 и Рг = 1 для разного времени, полученные с помощью вспомогательного оператора, улучшающего решение на границе.

Профили функции тока: штриховая 9=0.001, сплошная 9=0.01

Литература

1. Хаппель, Дж. Гидродинамика при малых числах Рейнольдса [Текст] / Дж. Хаппель, Г. Бреннер. - М.: Мир, 1976. - 630 с.

2. Богер, А. А. Расчет кондуктивно-ламинарного режима термоконвекции ньютоновской среды в прямоугольной каверне с вертикальными изотермическими стенками [Текст] / А.А. Богер, С.В. Рябов, В.И. Ряжских, М.И. Слюсарев // Изв. РАН. Механика жидкости и газа. -2010, № 3. - С. 17-21

3. Лойцянский, Л. Г. Механика жидкости и газа

[Текст] / Л. Г. Лойцянский. - М.: Наука, 1987. - 840 с.

4. Selvadurai, A. P. Partial differention equations in mechanics 2 [Text] / A. P. Selvadurai. - New York: Springer, 2004. - 698 p.

5. Слюсарев, М.И. Аналитическое решение первой тестовой задачи свободной конвекции для кондуктивно-ламинарного режима [Текст] / М.И. Слюсарев, Е.Ю. Чертов, В.И. Ряжских // Вестник Воронежского государственного технического университета. - 2010. - Т.6. № 7. - С. 165-167.

6. Сальвадори, М. Дж. Численные методы в технике [Текст] / М.Дж. Сальвадори. - М.: Иностранная литература, 1955. - 244 с.

7. Ряжских, В. И. Численное интегрирование бигар-монического уравнения в квадратной области [Текст] /

B.И. Ряжских, М.И. Слюсарев, М.И. Попов // Вестник СПБГУ. - 2013. - Сер.10, вып. 1. - С. 52-62.

8. Самарский, А.А. Вычислительная теплопередача [Текст] / А.А. Самарский, П.Н. Вабищевич. - М.: Едито-риал УРСС, 2003. - 784 с.

9. Лыков, А. В. Теория теплопроводности [текст] / А.В.Лыков. - М.: Высшая школа, 1967. - 600 с.

10. Самарский, А. А. Введение в теорию разностных схем [Текст] / А.А. Самарский. - М.: Наука, 1971. - 552 с.

11. Годунов, С. К. Разностные схемы [текст] /

C.К.Годунов, В.С.Рябенький. - М. Наука, 1977. - 440 с.

12. Фаддеев, Д. К. Вычислительные методы линейной алгебры [Текст] / Д.К.Фаддеев, В.Н.Фаддеева. - М.: ФИЗМАТЛИТ, 1963. - 735 с.

13. Ахиезер, Н. И. Теория линейных операторов в гильбертовом пространстве [Текст] / Н.И. Ахиезер, И.М. Глазман. - М.: Наука, 1966. - 543 с.

14. Филиппов, А. Ф. Об устойчивости разностных уравнений [Текст] / А.Ф.Филиппов, В.С.Рябенький. - М.: ГИТТЛ, 1956. - 171 с.

15. Андерсон, Д. Вычислительная механика и теплообмен [Текст] / Д. Андерсон, Дж. Таннехил, Р.Плетчер: В 2-х т. - Т. 1. - М.: Мир, 1990. - 384 с.

Воронежский государственный технический университет Воронежский государственный университет инженерных технологий

ABOUT THE NUMERICAL INTEGRATION OF THE NONSTATIONARY NON-HOMOGENEOUS BIHARMONIC EQUATION IN PROBLEMS OF CONDUCTIVE FREE

CONVECTION

V. I. Ryazhskih, M. I. Popov

Synthesis of the finite-difference scheme of a numerical integration of a nonstationary non-homogeneous biharmonic equation is presented at zero boundary conditions of the first and second sort. In the Hilbert space the analysis is carried out it. The approximation error is calculated, stability and convergence of the scheme are proved, the way of identification of an optimum step on time is offered. Approbation of computing procedure on a problem about conductive and laminar free convection in a square cavity confirmed effectiveness of the offered approach

Key words: biharmonic equation, conductive and laminar free convection, finite-difference scheme