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

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

Том 7, № 5, 2002

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

ЖИДКОСТИ

М.Ю. Балаганский, Ю.Н. Захаров, В. А. Ханефт Кемеровский государственный университет, Россия e-mail: stzun@ic.kemsu.ru

Nonstationary problems on ideal incompressible stratified fluid flow in Boussinesq approximation in endless channels are given. Implicit difference approximations are suggested. Iterative schemes with component-wise variational optimization are used.

Введение

Одной из простейших задач, описывающих движение идеальной несжимаемой экспоненциально стратифицированной жидкости в приближении Буссинеска в бесконечном канале, является следующая [1]:

t > 0, (1) Хз е (0,1), (2)

(3)

где x = (xi,x3); u = -0(x,t)e-ex3; в > 0, ^(x,t) — функция тока; ^ = 2в# — частота Вайсяля — Брента; g — ускорение свободного падения, направленное вниз по оси x3; П = {x1 е (0, то), x3 е (0,1)} — область решения; Г1, Г2 — стенки канала для x3 = 0, x3 = 1 соответственно; ^(x3,t) — заданная функция. Для существования и единственности решения этой задачи не требуется задания условия на бесконечности (при x1 ^ то) [1].

При решении разностными методами таких задач возникают несколько проблем, затрудняющих их использование. Во-первых, любая разумная разностная схема, аппроксимирующая уравнение (1), имеет пятиточечную часть шаблона на верхнем слое по времени, что затрудняет применение прогонки для реализации этой схемы; во-вторых, необходимость решения задачи (1)-(3) в конечной области (x1 < X < то) и использование разностных методов требуют задания на этой границе (x1 = X) краевого условия. Но, как показано в [1], даже при стационарных краевых условиях и финитной по времени функции ^(x3,t) может возникать так называемое незатухающее "реликтовое течение", которое не позволяет получать легко реализуемые краевые условия.

© М.Ю. Балаганский, Ю.Н. Захаров, В.А. Ханефт, 2002.

д2

dt2 [uxixi + ux3x3 в u] + ^0uxi XI — x е

u(x, 0) = ut(x, 0) = 0, x е П, u(x,t)| 0 = ^(x3; t), u(x,t)| „ = u(x,t)1 „ =0, t > 0,

4 ' 'Ix€i1 4 'Ix€i 2 —

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

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

Введем в области П прямоугольную неравномерную по х1 и хз согласованную с границей Г и Г2 сетку

П = (хн = хН-1 + Нн, хз? = хз^_1 + Нз?, % е I, з е 3} ,

где 1, 3 — множество индексов. Аппроксимируем на сетке Пл задачу (1)-(3) разностной задачей

(Аь - в2Е) ( -Д2-^ ) + А^1 = 0, (4)

и/1 • — ?/°

<+1 = ^(хз; ¿ш), «0 = = 0, = = 0, (5)

% = 1,... ,т1,..., з = 1, 2,...,т3 — 1, к = 1, 2,...

Здесь Ал = Аш + А

зл ,

1

„Л+1 — „Л+1

"¿+1;? г?

„Л+1 — ик+1

г? и'г_1?

1г

Н

1г+1

Н

1г

к+1.

1

Нз?

ик+1 — ик+1 "¿^'+1 "г?

ик+1 — ик+1 г? г;?_1

Н

3?+1

Н

3?

где Е — единичный оператор; Ни = 2 (Нн+1 + Нн), Нз? = 1 (Нз?+1 + Нз?); 4 = кАЬ, к = 0,1, 2,..., А£ — положительная постоянная.

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

д1_) — в 2е

л+1

М1? , ,2 д(_1). к+1

I + ^0 А1Л ит1 ?

0, з = 1, 2,..., тз — 1, (6)

аЛ_) =

д<_-)

,к+1

' - — 2ик+1 + иЛ д(_1)ик+1 = "тц' 2 ит1 _ 1 ? + ит1_2?

Н1т1

к+1

+ А

зЛ:

При построении сетки Пл мы полагаем Н1т1_1 = Н1т1.

Таким образом, чтобы найти решение и^1 разностной задачи (4), (6), необходимо для каждого к = 1, 2,... решать систему линейных алгебраических уравнений (СЛАУ)

Аи = /, (7)

где и = {и^1}, % е (1,т1), з е 3; f, — известный вектор правых частей, зависящий от краевых условий и от и?. Размерность вектора и равна количеству точек разбиения т области П1 = {х1 е (0,Х1 ], хз е (0,1)}. Очевидно, что даже в случае равномерной

по х1 сетки матрица А системы является несамосопряженной. Но если сетка является адаптивной или форма границы канала зависит от времени, то матрица А, оставаясь шестидиагональной, также будет меняться. Таким образом, решение задачи (4)-(6) может сводиться к многократному решению системы (7) с различными матрицами А. При этом свойства этих матриц, вообще говоря, трудно прогнозируемы. Следовательно, для решения задачи (4) - (6) необходимо использовать итерационные схемы, быстросходящиеся в случае произвольной матрицы А.

2. Итерационные схемы

Для решения (7) рассмотрим явную итерационную схему

ип+1/2 = ип _ Гп+1Гп (8)

ип+1 = ип+1/2 _ ап+^п, п = 0,1,2,..., (9)

и0 — начальные данные, которые мы задавали как решение задачи (4)-(6) на момент времени ¿к; гп = Аип _ ^ — невязка; ъп — произвольный вектор с ненулевыми компонентами; тп+1 выбирается из условия минимума нормы невязки гп+1/2 = Аип+1/2 _ f; ап+1 — диагональная матрица итерационных параметров |аР+1}, р =1, 2, ...,т, и, ^ £ Нт, Нт — т-мерное евклидово пространство.

Далее будем считать, что матрица А может быть незнакоопределенной и почти особенной, т. е. det А ^ 10-20. В этом случае система (7) формально имеет единственное решение, но если А будет незнакоопределенной матрицей, для решения (7) нельзя использовать 1-ю трансформацию Гаусса (т.е. умножить (7) на А*).

Пусть Д — линейный самосопряженный положительно определенный оператор в Нт. Относительно vn = Д1/2(ип _ и) схема (8), (9) имеет вид

vn+1/2 = Бп+1Уп, (10)

уп+1 = vn+1/2 _ Д1/2ап+12п, п = 0,1, 2,..., (11)

5п+1 = Е _ Тп+1Д1/2АД-1/2.

В [4] были предложены варианты выбора тп+1, ап+1 в виде некоторых констант. В продолжение работы [5], в которой излагаются некоторые алгоритмы выбора итерационного параметра ап+1 в виде матрицы, рассмотрим несколько способов выбора элементов матрицы ап+1, обеспечивающих полную или частичную оптимизацию перехода от vn+1/2 к

Перепишем (11) в виде

т

vn+1 = vn+1/2 _ ^ ап+1Д1/2ъп,

г=1

где «п+1 — элемент главной диагонали матрицы ап+1; Ър — вектор с одной ненулевой ¿-й

т

компонентой, такой что ^ Ър = ъп.

г=1

Будем находить ага+1 = diag{an+1} из условия минимума ||уга+11|2. Имеем

т

|уп+1м2 = ||уп+1/2|

г=1

тт

+ 2 ^(Я1^?, ) + К+^Р^, + 2 ^(Я1^, Я1/2^) + ■ ■ ■ +

г=2

г=3

+ (ат-1)2(^т-ь^т^) + 2 ^р1^-^1^) + к+^р^т,^2о (12)

Полагая частные производные от ||уп+1 и2 по

алгебраических уравнений

по аП+1 нулю, получаем систему линейных

) (Я1^ .Я1^») ... (Я1^, Я1/2*т)/

(Я1/2 zn, ^+1/2)

X

а.

п+1 2

п+1

\ат+1/

(13)

В силу самосопряженности и положительной определенности матрицы Я система (13) дает условие минимума |^га+1|| . Если система векторов ... , zm} линейно независима, то матрица системы (13) является матрицей Грама и ее определитель не равен нулю. Тем самым (13) имеет единственное решение («П+1,... , а»г+1)Т, и оно доставляет минимум ||уга+1||2 итерационного шага (11).

В общем случае, когда матрица полностью заполнена, решение полученной СЛАУ достаточно сложная задача, но при решении разностных задач вида (4) - (6) можно использовать различные способы или точного, или приближенного решений задачи минимизации

V

п+11

Рассмотрим несколько алгоритмов выбора «П+1.

Алгоритм 1. Пусть I = (г : 1 < г < т} — множество всех индексов. Разобьем это множество на непересекающиеся подмножества /1г .. ,1р. Тогда

^+1 = ^+1/2

г

ге/1

ге/2

Е

ге/„

<+1^4

(14)

параметры «П+1 при г € 11 будем находить из условия минимума

К+1/2 «П+1Я1^

ге/1

«П+1 при г € 12 — из условия минимума

п+1/2

ге/1

ге/2

и т. д.

т

Для нахождения а'+1, г £ 11,... ,1Р необходимо каждый раз решать систему вида (13), только меньшей размерности. Таким образом, можно подобрать такое число р, что находить итерационные параметры будет несложно, и получить схему, в которой оптимизируется переход к ип+1 по группам компонент ип+1/2.

Далее рассмотрим случай, когда матрица Б1/2 имеет шесть ненулевых диагоналей ("ше-стидиагональна"). Тогда матрица Б будет "одиннадцатидиагональна" и у вектора Бх' будут только 11 ненулевых компонент. Вследствие этого матрица ZD системы (13) также является "одиннадцатидиагональной".

Алгоритм 2. Разобьем множество всех индексов I = {г : 1 < г < т} на непересекающиеся подмножества 11,... ,1& по следующему принципу: индексы к, I принадлежат одному множеству , если

(Б1/2хпк,Б1/2Х') = 0, к = I. (15)

Заметим, что для "одиннадцатидиагональной" матрицы наименьшим числом групп будет 6 (р = 6). Далее, используя алгоритм 1, будем находить итерационные параметры. За счет особой группировки индексов получающиеся СЛАУ, подобные по структуре (13), будут диагональными. Решение таких систем не составляет труда.

Итак, построен алгоритм, с помощью которого без дополнительных вычислительных затрат минимизируется ||уп+11| по большему числу параметров одновременно. Если гп = (1,... , 1), разбиение множества индексов на подмножества нужно произвести только один раз, следовательно, количество затрачиваемых операций возрастет незначительно.

Алгоритм 3. Упростим СЛАУ (13) без понижения ее размерности, занулив элементы вне главной диагонали. Для этого необходимо получить некоторый набор векторов, удовлетворяющих свойству (15) при любых к, I. Это требование приводит нас к необходимости ортогонализации (в пространстве Ир) линейно-независимой системы векторов

¿П хп ¿1 ,. . . ,хт.

Выбор этих векторов осуществим следующим образом:

1. За основу возьмем группу с максимально возможным числом векторов, которые удовлетворяют свойству (15). В случае шестидиагональной матрицы А такая группа имеет размерность порядка т/6.

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

3. Так как матрица системы (13) имеет теперь диагональный вид,

, (уп+1/2,О1/2хп) ап+1 =-:-2—, г = 1,...,т,

ЦБ1/2

п

2„ ,,

где {х'П,... , ¿т} — ортогональная система векторов.

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

Алгоритм 4. Если множества 11,... ,1Р — пересекающиеся или, тем более, среди них есть совпадающие и для выбора {а^+1} используется алгоритм 1, то это означает, что часть компонент вектора ип+1 пересчитывается, а часть нет. Следовательно, этот алгоритм позволяет считать задачи движения жидкости только в тех областях, где решение

не установилось и достаточно быстро меняется. Это позволяет существенно уменьшить время расчета.

Из вышеприведенного можно сделать следующий вывод. Введение в итерационную схему итерационного параметра в виде матрицы позволяет построить сходящийся оптимальным образом (см. также, [5, 6]) алгоритм решения СЛАУ с незнакоопределенной и почти особенной матрицей, не использующий никаких преобразований исходной системы. При этом алгоритм позволяет при переходе от итерации к итерации считать только там, где решение еще не установилось.

3. Численные расчеты

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

Для уравнения (1) рассматривались следующие начальные и краевые условия. Задача 1.

«(ж, £) I

0, и(ж,£)|

1жз€ГЬГ2 1x1=0

«(ж,0) = «¿(ж,0) = 0, х\ е [0,то), х3 е [0,1],

где <р(жз,г)

вт(2пх3 + п)вт(п£), £ = 0

На дне канала расположено препятствие,

0, £ > 0

описываемое функцией вида у (Ж1)

середина препятствия по х1; Г1, Г2 (рис. 3).

Ь

е-(х1-«)2 — е

-0.04

Ь

высота препятствия; а

1 - е-0.04

верхняя и нижняя стенки канала соответственно

Задача 2.

Рис. 1.

Рис. 2.

и(х, £) |хзеГ1 Г2 = 0, и(х, 0) = 0, «¿(ж, 0) = 0,

х1 е (-то, то), ж3 е [0,1],

где Г1 — верхняя граница канала, Г2 — нижняя граница канала (рис. 2). На дне канала расположена тонкая, толщиной Ь = 0.05, колеблющаяся пластинка высотой Ь = 0.3 с

Рис. 3.

Для решения разностных задач на каждом шаге использовалась итерационная схема (8), (9) с различными способами выбора ага+1, как в виде постоянной, так и в виде матрицы (см. более подробно в [5]). В качестве начальных данных для схемы (8), (9) выбиралось решение исходной разностной схемы на время ¿д.. Во всех расчетах итерации на каждом шаге = (к + 1)т останавливались при выполнении условий ||гга+1||/||г0|| < 10-3.

На рис. ?? приведены рассчитанные мгновенные линии тока для задач 1, 2 (рис. 3, а - в и 3, г, д соответственно) в различные моменты времени. Результаты, представленные на

рис. 3, а-3, г, получены при решении задачи в области П = {[0,1] х [0,1]}, £ € [0,1], т1 = 100, т3 = 55, Д£ = 0.00025, Ь = 0.3, а = 0.5, а на рис. 3, д, е — в области П = {[0,1] х [0,1]}, £ € [0,1], т1 = 50, т3 = 50, Д£ = 0.001, высота пластинки Ь = 0.3, а = 0.5. На свободных границах (х1 = 0, х1 = Х1 = 1) аналогично (6) аппроксимировалось уравнение (1) внутрь области П.

Для проверки правильности проведенных расчетов мы увеличили область решения по х1 до Х1 > Х1, оставляя неизменными сетки по пространству и времени. Затем сравнивали решения в общей области. Оказалось, что при увеличении области в два раза евклидова норма разности решений не превышает 2%.

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

[1] ГАБов С. А., Свешников А. Г. Линейные задачи теории нестационарных внутренних волн. М.: Наука, 1990. 342 с.

[2] Захаров Ю.Н. Итерационные схемы неполной аппроксимации / / Численные методы механики сплошной среды: Сб. науч. тр. / СО АН СССР. ИТПМ. 1985. Т. 16, №6. С. 77-83.

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

[4] ЗАхАров Ю. Н. Итерационные схемы неполной аппроксимации решения систем линейных уравнений с незнакоопределенной матрицей // Конструирование алгоритмов и решение задач математической физики: Сб. науч. тр. М.: Наука, 1989. С. 197-201.

[5] Захаров Ю. Н., Терешкова В. В., Шокин Ю. Н. Об одном классе итерационных схем решения систем линейных уравнений с незнакоопределенной матрицей. Красноярск, 1990 (Препр. / АН СССР. Сиб отд-ние. ВЦ; №14). 22 с.

[6] ЗАхАров Ю. Н. Об одном методе последовательных приближений решения линейных операторных уравнений // Теория функций и ее приложения, Кемерово, 1985. С. 86-89.

Поступила в редакцию 13 августа 2002 г.