Применение итерационных схем при решении задач о движении стратифицированной жидкости Текст научной статьи по специальности «Математика»

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

Том 11, Специальный выпуск, 2006

ПРИМЕНЕНИЕ ИТЕРАЦИОННЫХ СХЕМ ПРИ РЕШЕНИИ ЗАДАЧ О ДВИЖЕНИИ СТРАТИФИЦИРОВАННОЙ ЖИДКОСТИ

А. С. Васильев, Ю. Н. Захаров Кемеровский государственный университет, Россия e-mail: vasves@rambler.ru, junz@kemsu.ru

A multi-step iterative scheme is used to solve evolutional problems of an ideal stratified fluid dynamics. The proposed scheme has the convergence rate that is similar compared with the Chebyshev's method and it is stable with regard to variations of the initial assignment of the input information which is required for its construction. Results of application of the proposed method to a number of model problems are given.

Введение

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

При численной реализации моделей движения жидкости со стратификацией по плотности возникают системы линейных алгебраических уравнений с одинаковой на каждом слое по времени матрицей и различными правыми частями. Свойства матриц этих систем таковы, что решение их с помощью методов расщепления [2-4] сопряжено с рядом трудностей, использование же итерационных методов на каждом временном слое требует знания свойств матрицы решаемой системы. Например, применение "чебышевского" метода затруднено, так как могут быть неизвестны минимальное и максимальное собственные числа матрицы или она может быть несамосопряженной. Использование же градиентных методов нецелесообразно из-за большого количества арифметических операций при их реализации. В связи с этим возникает необходимость в методах, которые имеют скорость сходимости, как и у "чебышевского" итерационного метода, и при этом используют минимальную информацию о свойствах матрицы решаемой системы линейных алгебраических уравнений.

© Институт вычислительных технологий Сибирского отделения Российской академии наук, 2006.

1. Постановка задачи

Рассмотрим двумерные движения идеальной стратифицированной жидкости, В декартовой системе координат (х, х3) такие движения описываются уравнением

д2

Ор [их 1*1 + - Р2и\ + ^0^1*1 = 0, (1)

где и = ф (х1; х3, £) е—вхз, ф — функция тока; = 2вд — частота Вейсяля — Брента; в — некоторая положительная постоянная, задающая параметры стратификации жидкости; д — ускорение свободного падения [1]. К уравнению (1) необходимо добавить следующие начальные условия:

и (х, 0) = иг (х, 0) = 0. (2)

Исходное уравнение (1) решается в области

О = (х е (х1, х3) е Я2 : 0 < х1 < 0 < х3 < Ь3) ,

представляющей собой полубесконечный прямоугольный канал. На ее границе ставятся условия

и (х,^)|хз=0,хз=Ьз =0; (3)

и (х,^)|х1=0 = ^ (хз, (4)

где ^ (х3, £) — известная функция.

Задача (1)-(4) решается в бесконечной области, и, как показано в [1], для существования и единственности решения не требуется задавать условия на бесконечности.

2. Аппроксимация

Дифференциальное уравнение (1) аппроксимируем на неравномерной прямоугольной сетке узлов О^ = {хн = х1г-1 + Ни, х3] = х3^-1 + Н3]} , г е I, ] е 3 (I, 3 — множество индексов), неявной разностной схемой

(дл - в2Е)

и

П 2иП] + иП]

] 1

д^2

+ ^02Д1ьиП]+1 = 0, г е I, ] е з.

(5)

Здесь Дн = Д^ + Дзь;

Д^Г = т-

«1г

1 I и^+1 ] и] и] и»

Ни

1»+1

Нт»

п+1

иг]+1 иг]

3]

Н

3]+1

— 1

Н

3]

(6)

Е - единичный оператор; по времени, п = 1, 2, 3,...

= пД£; Д£ — шаг

1

Таким образом, для решения поставленной дискретной задачи необходимо на каждом (п + 1)-м слое по времени решать систему линейных алгебраических уравнений

Ли = /, (7)

где Л — пятидиагональная матрица системы; и = {и"+1, г € I, ] € — неизвестный

/п (п — 1)-м слоях по времени.

Так как область решения уравнения (1) не ограничена, матрица и векторы системы (7) будут содержать бесконечное число элементов. Поэтому при численном решении задачи (1)-(4) необходимо ограничивать область О прямой х1 = Ь1, на которой должно выполняться само уравнение (1), вследствие чего в схеме (5) индекс г ограничен постоянной т1.

В этом случае система (7) будет недоопределенной и возникнет необходимость ее доопределения па прямой х1т1 = Ь1. Если считать, что па этой прямой течение слабо изменяется со временем, то можно поставить какое-либо краевое условие, не искажающее решения [5]. Однако в нашем случае характер течения нестационарной жидкости на рассматриваемом участке границы определить трудно, так как оно нетривиально [1, 6], Поскольку на прямой х1т1 = Ь1 должно быть выполнено само уравнение (1), предлагается замыкать систему (7) путем аппроксимации этого уравнения внутрь области решения. Такая аппроксимация, например, может иметь вид

(и*+1— оип + ип-1\ (А, - Р2Е) ( ^ + ) + и,02Д1Л«К = О, з = .1; (8)

Д1Ь и

ип+1 — ип+1 ип+1 — ип+1 п+1 _ 1 I итц' т1 — 1у _ ит1-1з ит1-2]

д „ п+1 _ 1 т 1 3+1 "т^ "т^ "т 13-1 I /пч

Ш ня 1' (9)

^1т1 — 1 \ h1ml hlml — 1 1 / ип+1 _ ип+1 ип+1 _ ип+1

Ч? \ "33+1 "33

Для решения подобных систем существует достаточно богатый арсенал итерационных методов [3], Например, пользуясь тем, что матрица системы (7) не меняется при переходе от слоя к слою по времени, выгодно было бы использовать итерационный метод Чебышева, Но для его реализации требуется найти значения спектральных границ 7^ 72 операто-Л

использовании приближенных значений 71, 72 итерационный метод Чебышева может быть неустойчивым [3],

Вышесказанное верно, если матрица системы самосопряжена, но в нашем случае Л = Л*. В связи с этим возникает необходимость в использовании метода, который имеет оценку сходимости "чебышевского" итерационного метода в самосопряженном случае, сходится при неточном задании входных параметров и может быть применен в случае

Л

говая итерационная схема [6, 7].

3. Метод решения

Пусть В — самосопряженный положительно определенный оператор, действующий в гильбертовом пространстве И. Введем энергетическое подпространство Нсостоящее из эле-

ментов Н со скалярным произведением (х,у)Г) = (Их, у) и нормой ЦжЦд = л/(Их, х). Пусть уравнение (7) задано в этом пространстве И, * Рассмотрим многошаговую итерационную схему [7]

где

ик+2р _

вр-+ Аик = /, (10)

т*

Вк — —-Вк_ 1 (Е + Ск-\) 1;

1 + ^к

Ск-1 = Е _ г*в-11А;

ик= Шк \—, к = 1, 2,... ,р ; В0 = В(т*)] ш0 = р*; (11)

2 _ шк-1

т* — значение оптимального итерационного параметра; £0 — оператор шага и р* = ||£0|| — норма оператора шага одношаговой итерационной схемы при оптимальном значении итерационного параметра.

Для схемы (10), (11) в Ив в случае, когда С = С * = Д1/2В-1АД-1/2 > 0, справедлива следующая оценка [7|:

где (р*) = 2

К+2Р < (р*) |И|Д

2Р , I-ч 2Р

1 + 31-1 +1-31-1

р* р* р* р*

-1

т* р*

"чебышевский" итерационный процесс, т.е. схема (10), (11) есть другая форма записи 2р-шагового "чебышевского" итерационного метода.

Предположим, что вместо т* известно некоторое его приближение т0, выбранное, например, экспериментально, при котором ||$0 (т0)|| = р0 > р*. В этом случае точное значение р0 обычно неизвестно и вместо него используется некоторая оценка рс < 1, которую можно взять в (11) в качестве о>0 = рс. Можно показать [7], что если

0<Р1<-^, (12)

1 + р0

то для оператора Ск справедлива оценка

||Ск|| = рк < рк-1 < ... < р0к = 1, 2,..., р,

где

2рк-1 Шек-1 п . 2 ^ 2 —---, 0 < рс < Рол = / 2 _ Шск-1 рк = ^ , ,2 0 2 Шск-1 2 ^ 2 ^ 2р0

,2 , РО < Рс <

2 _ Г0 гс 1 + р2

Когда Ро ^ Рс ^ -2' значение Рк = (Рс) И при этом схема сходится со скоростью

1 + р20

;'чебышевского" итерационного процесса в соответствии со значением рс. Если значение рс

удовлетворяет неравенству

1 + Р§

< р2 < 1,

то схема (10), (11) так же является сходящейся. Таким образом, она будет сходиться для всех р2 € (0,1).

Используя соотношения Ск = (1 + ок) — окЕ (к = 1, 2,..., г), С0 = Е _ тВ-1А, можно построить следующую реализацию схемы (10), (11):

ик+1 _

В--— + Аик = /, к = 0,1,..., 2', (13)

Тк

где тк = то + а^) > к =1, 2,..., 2Р. Коэффициенты а^ вычисляются по рекуррентным формулам

= ^ «S = ^1-1 + (1 + ^i-i) ± Ja^ (i + , l = 2,..., p, (14)

= «22 ± ^(1 + а£), Г = 1, 2,..., 21, I = 1, 2,... (15)

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

Следует заметить, что самосопряженность оператора ДВ-1А используется только при оценке скорости сходимости алгоритма (13)—(15) и не требуется при его реализации. Это делает возможным применение многошаговой итерационной схемы и в случае несамосопряженного оператора А.

4. Модельные расчеты

Для практической проверки свойств рассмотренного метода проведена серия численных экспериментов на модельной задаче, в качестве которой была использована задача Дирихле для уравнения Пуассона, заданного в единичном квадрате с границей Г:

—AU = 2n2sin (—nx1) sin (nx2); (16)

U |г = 0. (17)

Уравнение (16) аппроксимировалось па квадратной сетке узлов с шагом h = 1/M разностной схемой второго порядка. Полученная система линейных уранений имеет самосопряженную положительно определенную матрицу, границы спектра которой известны точно.

Решение этой системы проводилось явным 2р-циклическим методом Чебышева с упорядочиванием параметоров по Самарскому [3] и многошаговой схемой (13)—(15). Начальный вектор во всех расчетах брался единичным. Итерации останавливались при выполнении условия ||rk||/||r0|| ^ 10-4.

Таблица 1. Результаты вычислений для самосопряженного оператора

V М = = 50 М = 100 М = 150

Л^сЬ Nm Л^сЬ Nm Л^сЬ Nm

2 884 723 3420 2867 7864 6671

3 528 415 1864 1279 4112 3127

4 304 206 1040 702 2160 1438

5 223 151 608 350 1184 670

6 176 134 447 253 768 382

Таблица 2. Количество итераций, в зависимости от рс

Рс Nm

0.001 2253

0.400 2065

0.990 308

р*= 0.998026 134

0.999 158

0.9999 446

0.99999 1791

0.999999 7167

Таблица 3. Количество итераций для решения задач с использованием многошаговой схемы (для несамосопряженного оператора)

V М = 25 М = 150

т0 = 0.0001 т0 = 0.0003 т0 = 0.0004 т0 = 0.000005 т0 = 0.00001 т0 = 0.000011

2 1704 564 414 25216 12622 11478

3 1074 352 241 12962 6512 5928

4 896 288 176 7168 3648 3328

5 872 276 160 4832 2512 2304

В табл. 1 приведены результаты расчетов при точно заданных границах спектра 71 и 72 оператора А, Здесь р — величина, определяющая длину цикла итерационных параметров; и — количество итераций при использовании "чебышевского" метода и итерационной схемы,

В табл. 2 приведено число итераций, получаемое при решении с помощью схемы (10), (11) при р = 6 и М = 50 для различных значений рс. Итерационный параметр т* при этом

2

задавался точно — т*

71 + 72

Поскольку при построении и реализации многошаговой схемы (10), (11) не используется свойство самосопряженности оператора А, мы посчитали целесообразным проверить возможность ее применения в несамосопряженном случае. Для получения системы с несамосопряженным оператором в задаче (16), (17) па стороне квадрата х1 = 1 граничное

условие и\ = 0 заменили на (тг~ — 1000С/ ) 1 \дх1 )

= — пет (х2) и эту дифференци-

Х1 = 1

альную задачу аппроксимировали на квадратной сетке узлов разностной схемой второго

х1 = 1

мосопряженной и знакоопределенной матрицей А,

В табл. 3 приведены результаты вычислений при использовании многошаговой схемы для М = 25, М = 150 и различных значений т0, выбранных экспериментально. Итерации останавливались при выполнении условия ||гк||/||г0|| ^ 10-7, Для определения значения параметра рс, необходимого для схемы (13)-(15), пр и выбранном т0 проводил ось N итераций с применением одношаговой схемы, значение рс определялось с помощью выражения рс = ||иМс — иМс-1 ||/||иМс-1 — иМс-2Ц, Для М = 25 количество итераций N = 100, а для М = 150 — N = 1000.

5. Решение задачи о течении стратифицированной жидкости

При решении задачи (1)-(4) функция ^ (х1, хз) из (4) задавалась в виде

^ (хз ,*) = 5 (хз ,*) хзе вХ3

5 (хз ,*)

1, (1 _ 2 ■ 10_4 (* _ 10_3)) < хз < 1, * < 6 ■ 10 6 ■ 10_з.

-з

(18)

(19)

Здесь 5 — параметр, влияющий на скорость потока на границе х1 =0 для времени * > 6 ■ 10_з.

Функция (18) задана таким образом, что на входе в канал у вектора скорости отсутствует вертикальная составляющая.

Во время всех вычислений О = {х1 € [0, 6]; хз € [0,1]} число узлов = 300, МХ3 = 50, шаг по времени А* = 10_з, На первом этапе для найденных экспериментально значений т0, рс согласно (14), (15) строился набор итерационных параметров в количестве 2Р (р = 15).

Далее па каждом шаге по времени решение ига+1 разностной схемы (5), (6) и (8), (9) находилось с помощью схемы (10), (11). В качестве начального приближения бралось решение на п-м слое по времени. Вычисления останавливались при выполнении условия ||гк|| /||г°|| ^ 10_6. Выбранная точность вычислений объясняется тем, что при переходе от п-го к (п + 1)-му слою по времени значение ||г0|| может достигать порядка 102,

г = 0.02

г = 0.5

г = 1

1 0.8 0.6 0.4 0.2 о.

1 0.8 0.6 0.4 0.2 0,

0 12 3

г = 0.02

'о

1 0.8 0.6 0.4 0.2 0,

1 2

г = 0.5

1 2

г=1

0

Линии тока при в = 0 (сверху) и в = 0.05 (снизу).

Полученное с помощью описанного алгоритма значение un+1 считалось решением разностной схемы на (n + 1)-м слое по времени. При этом среднее количество итераций на каждом слое составило 178,

На рисунке показаны линии уровня функции Ф для трех значений времени (t = 0.02, t = 0.5, t = 1). Во всех расчетах значение параметра стратификации ß = 5, Показана та область решения, где наблюдается возмущение. При s = 0 для времени t > 6 ■ 10-3 жидкость не проходит через границу области x1 = 0, В случае s = 0.05 для любого значения t ^ 0 через границу xi = 0 поступает жидкость.

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

[1] Габов С.А., Свешников А.Г. Задачи динамики стратифицированных жидкостей. М.: Наука, 1986. 288 с.

[2] Яненко H.H. Метод дробных шагов решения многомерных задач математической физики. М.: Наука, 1967. 196 с.

[3] Самарский A.A., Николаев Е.С. Методы решения сеточных уравнений. М.: Наука, 1978. 592 с.

[4] Марчук Г.И. Методы вычислительной математики. М.: Наука, 1989. 608 с.

[5] Ильгамов М.А. Обзор работ по неотражающим условиям на границах расчетной области // Тр. семинара Казан, физ.-техн. ин-та. 1990. № 4. С. 6-54.

[6] Захаров Ю.Н. Градиентные итерационные методы решения задач гидродинамики. Новосибирск: Наука, 2005. 239 с.

[7] Захаров Ю.Н. Об одном способе построения циклических итерационных схем // Числ. методы механики сплошной среды: Сб. науч.тр. / АН СССР. ВЦ; I I I I IM. 1979. Т. 10, № 4. С. 85-100.

Поступила в редакцию 20 февраля 2006 г. в переработанном виде — 3 апреля 2006 г.