Быстрые методы вычисления характеристик гидродинамической устойчивости

Демьянко, Кирилл Вячеславович

Быстрые методы вычисления характеристик гидродинамической устойчивости тема автореферата и диссертации по математике, 01.01.07 ВАК РФ

Демьянко, Кирилл Вячеславович АВТОР

кандидата физико-математических наук УЧЕНАЯ СТЕПЕНЬ

Москва МЕСТО ЗАЩИТЫ

2014 ГОД ЗАЩИТЫ

01.01.07 КОД ВАК РФ

Автореферат по математике на тему «Быстрые методы вычисления характеристик гидродинамической устойчивости»

Автореферат диссертации на тему "Быстрые методы вычисления характеристик гидродинамической устойчивости"

На правах рукописи

Демьянко Кирилл Вячеславович

Быстрые методы вычисления характеристик гидродинамической устойчивости

Специальность 01.01.07 — вычислительная математика

Автореферат диссертации на соискание учёной степени кандидата физико-математических наук

15 ЯНВ 2015

005557576

Москва — 2014

005557576

Работа выполнена в федеральном государственном автономном образовательном учреждении высшего профессионального образования «Московский физико-технический институт (государственный университет)»

Научный руководитель: Нечепуренко Юрий Михайлович

доктор физико-математических наук, доцент.

Официальные оппоненты: Амосов Андрей Авенирович,

доктор физико-математических наук, профессор. Федеральное государственное бюджетное образовательное учреждение высшего профессионального образования «Национальный исследовательский университет МЭИ», заведующий кафедрой математического моделирования.

Богачев Кирилл Юрьевич,

доктор физико-математических наук, доцент. Федеральное государственное бюджетное образовательное учреждение высшего профессионального образования «Московский государственный университет имени М.В. Ломоносова», доцент.

Ведущая организация: Федеральное государственное бюджетное учреждение науки Вычислительный центр им. A.A. Дородницына Российской академии наук.

Защита состоится 22 января 2015 г. в 1600 часов на заседании диссертационного совета Д 002.045.01 при Федеральном государственном бюджетном учреждении науки Институте вычислительной математики Российской академии наук (ИВМ РАН), расположенном по адресу: 119333, г. Москва, ул. Губкина, д. 8.

С диссертацией можно ознакомиться в библиотеке и на сайте ИВМ РАН http://www.inm.ras.ru.

Автореферат разослан fß Soj? 2014 года.

V /

Ученый секретарь

диссертационного совета Д 002.045.01 доктор физико-математических наук.

Ml^Z Бочаров Г.А.

Общая характеристика работы

Актуальность темы. Задачи гидродинамической устойчивости возникают, например, при проектировании судов и дозвуковых летательных аппаратов и находятся в центре внимания многих исследователей начиная с конца 19-го века благодаря пионерским работам Гельмгольца, Кельвина, Релеея, Рейнольдса.

Основными характеристиками гидродинамической устойчивости являются энергетическое и линейное критические числа Рейнольдса. Их вычисление сводится к решению частичных обобщенных проблем собственных значений для систем обыкновенных дифференциальных и алгебраических уравнений, возникающих после пространственной аппроксимации линеаризованных уравнений вязкой несжимаемой жидкости. Если течение зависит от одной пространственной переменной или сводится к таковому, то для расчетов достаточно современного ПК и стандартного численного программного обеспечения. Если же скорость основного течения зависит от двух (биглобальная устойчивость) или трех (триглобальная устойчивость) пространственных переменных и/или времени, то, в силу существенно большей алгебраической размерности возникающих вычислительных задач, необходимы специальные алгоритмы, построенные с учетом структуры исследуемых уравнений. Кроме того, актуальна разработка метода вычисления линейного критического числа Рейнольдса с заданной точностью, чего не позволяют делать традиционные подходы.

Одной из задач биглобальной устойчивости является зависимость линейного критического числа Рейнольдса течения Пуазейля в бесконечном канале постоянного прямоугольного сечения от величины Z > 1 отношения длин сторон сечения. Известно, что существует Zc > 1, такое что при Z < Zc линейное критическое число Рейнольдса бесконечно, а при Z > Zc - конечно. В работах Т. Tatsumi, Т. Yoshimura (1990) и V. Theofilis, P.W. Duck, J. Owen (2004) было численно установлено, что Zc ~ 3.2. Однако используемые в этих работах методы не позволяли судить о погрешности результата. Более того, какого-либо теоретического обоснования данной зависимости до недавнего времени получить не удавалось.

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

но простое течение, то для пространственной аппроксимации часто используют метод коллокаций, который приводит к обобщенным проблемам собственных значений с плотными матрицами не очень большого размера. При исследовании сложных течений, метод коллокаций неприменим из-за слишком высокого размера получающихся плотных матриц. В этом случае приходится использовать конечно-разностные или конечно-элементные методы аппроксимации, приводящие к задачам с разреженными матрицами. Наиболее популярными методами решения частичной проблемы собственных значений с большими разреженными матрицами являются метод Арнольди, несимметричный метод Ланцоша и метод Якоби-Дэвидсона. Существенным недостатком этих методов является высокое требование к объему оперативной памяти, так как приближение к собственному вектору ищется в подпространстве, размерность которого растет на каждом шаге, а для нахождения инвариантного или понижающего подпространства приходится использовать более сложные и вычислительно-емкие блочные варианты этих методов. Таким образом, разработка более эффективных алгоритмов решения частичных проблем собственных значений с большими разреженными матрицами, по-прежнему остается актуальной задачей.

Целью диссертационной работы является разработка, обоснование и реализация быстрых методов вычисления характеристик гидродинамической устойчивости. В том числе, развитие и обоснование предложенной A.B. Бойко и Ю.М. Нечепуренко технологии численного анализа систем обыкновенных дифференциальных и алгебраических уравнений, полученных после пространственной аппроксимации линеаризованных уравнений вязкой несжимаемой жидкости. Создание специального варианта технологии для исследования устойчивости течений в бесконечных каналах постоянного прямоугольного сечения. Исследование с его помощью зависимости критического числа Рейнольдса течения Пуазейля от отношения длин сторон сечения канала. Теоретическое обоснование этой зависимости. Разработка, обоснование и реализация методов ньютоновского типа для решения частичных обычных и обобщенных проблем собственных значений с большими разреженными матрицами, возникающих при вычислении критических чисел Рейнольдса.

Основные положения, выносимые на защиту:

1. Развита и обоснована предложенная A.B. Бойко и Ю.М. Нечепуренко технология численного анализа систем обыкновенных дифференциальных и алгебраических уравнений, полученных после пространственной аппроксимации линеаризованных уравнений вязкой несжимаемой жидкости. В частности, обоснованы предварительные преобразования исходной системы, а также предложены и обоснованы модификации стандартных процедур FZERO и FMIN, позволяющие вычислять линейное критическое число Рейнольдса и строить соответствующие нейтральные кривые с заданной относительной точностью.

2. Разработан и обоснован специальный вариант технологии для исследования устойчивости течений в бесконечных каналах постоянного прямоугольного сечения. Численно исследована зависимость линейного критического числа Рейнольдса течения Пуазейля от отношения длин сторон сечения. Выполненные расчеты позволили уточнить известные результаты.

3. Впервые дано теоретическое обоснование зависимости линейного критического числа Рейнольдса течения Пуазейля в бесконечном канале постоянного прямоугольного сечения от отношения длин сторон сечения.

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

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

Научная новизна. Использование предложенных и обоснованных в диссертационной работе модификаций стандартных процедур FZERO и FMIN, позволяет вычислять линейные критические числа Рейнольдса и строить соответствующие нейтральные кривые с заданной относительной точностью, чего не позволяют делать традиционные подходы.

Разработанный и обоснованный специальный вариант технологии для

исследования устойчивости течений в бесконечных каналах постоянного пря-

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

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

Теоретическая и практическая значимость работы. Теоретическая значимость работы заключается в обосновании предварительных преобразований исходной системы обыкновенных дифференциальных и алгебраических уравнений, полученных после пространственной аппроксимации линеаризованных уравнений вязкой несжимаемой жидкости, обосновании предложенных модификаций процедур Р2ЕГЮ и РМГКГ, обосновании специального варианта метода коллокаций, сохраняющего на дискретном уровне симметричность и отрицательную определенность оператора Лапласа и соотношение (Ну = —дгас1*, теоретическом обосновании зависимости линейного критического числа Рейнольдса течения Пуазейля в бесконечном канале постоянного прямоугольного сечения от отношения длин сторон сечения, и, кроме того, в обосновании методов ньютоновского типа для решения частичных проблем собственных значений с большими разреженными матрицами, возникающих при вычислении характеристик гидродинамической устойчивости.

Практическая значимость работы состоит в создании нового метода вычисления критических чисел Рейнольдса, применимого как в случае плотных, так и в случае разреженных матриц. В отличие от других известных подходов, этот метод позволяет вычислять линейные критические числа Рейнольдса и строить соответствующие нейтральные кривые с заданной относительной точностью, а также является более экономичным. Кроме того, разработаны надежные,

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

Апробация работы. Основные результаты работы докладывались на: международной конференции «Модели и методы аэродинамики» (Евпатория, 2011), международной конференции «Современные проблемы вычислительной математики и математической физики» (Москва, 2014), XX всероссийской конференции «Теоретические основы и конструирование численных алгоритмов и решения задач математической физики» (Новороссийск, пос. Абрау-Дюрсо, 2014), всероссийских конференциях МФТИ (Москва, 2011-2013 г.г.). Результаты диссертации также обсуждались на научных семинарах Института вычислительной математики РАН и Института прикладной математики РАН им. М.В. Келдыша.

Результаты работы были отмечены почетным дипломом на научной конференции МФТИ-56 в 2013-м году.

Личный вклад. Теоретические результаты, представленные в работах [1-5], получены совместно с Ю.М. Нечепуренко. Теоретические результаты, представленные в [6], получены совместно с Ю.М. Нечепуренко и М. Садка-ном (Франция).

Реализация описанных в работах [1-6] алгоритмов и подготовка и проведение соответствующих численных экспериментов была выполнена автором диссертации самостоятельно.

Публикации. Основные результаты диссертации изложены в 9 печатных работах [1-9], из них 3 - в журналах, рекомендованных ВАК [1-3].

Объем и структура работы. Диссертация состоит из введения, трех глав и заключения. Полный объем диссертации составляет 122 страницы с 12 рисунками и 9 таблицами. Список литературы содержит 73 наименования.

Содержание работы

Во введении обосновывается актуальность исследований, проводимых в рамках данной диссертационной работы, формулируется цель, ставятся задачи,

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

Первая глава посвящена описанию, обоснованию и развитию технологии численного анализа систем обыкновенных дифференциальных и алгебраических уравнений, полученных после аппроксимации по пространственным переменным исходных линеаризованных уравнений вязкой несжимаемой жидкости. Технология была предложена в работах A.B. Бойко и Ю.М. Нечепуренко в 2009-м году, однако полного ее обоснования дано не было.

В разделе 1.2 дается постановка задачи. Рассматривается полученная после аппроксимации по пространственным переменным исходных линеаризованных уравнений Навье-Стокса и уравнения неразрывности система обыкновенных дифференциальных и алгебраических уравнений

= jWv+J-jWv + Gp, Fv = 0, (1)

at Ке

где первое уравнение является дискретным аналогом уравнений движения, а второе - уравнения неразрывности, v - покомпонентный вектор дискретных аналогов компонент скорости, р - п;,-компонентный вектор дискретного аналога давления (пр < nv); С и J'1) и J^2' - квадратные матрицы порядка п„, не зависящие от числа Рейнольдса, a G и F - прямоугольные матрицы размеров nv х ПрИПр х nv соответственно.

Предполагается, что (1) наследует известное свойство уравнений вязкой несжимаемой жидкости, а именно, что из уравнений движения можно исключить давление, используя уравнение неразрывности. Можно показать, что (1) обладает этим свойством, если матрица С невырождена и rankG = rank F = rank FC~lG. Действительно, если матрицы F и G не полного ранга и имеют дефект ранга d, то некоторые <1 строк матрицы F являются линейными комбинациями остальных ее строк, а некоторые d столбцов матрицы G являются линейными комбинациями остальных ее столбцов. Поэтому эти строки и столбцы можно сократить, сократив при этом и соответствующие компоненты вектора р. Таким образом из (1) получается система точно такого же вида с пр = n°ld — d и с матрицами F и G, удовлетворяющими условию

rank G = rank F = rank FC_1G = np, (2)

что эквивалентно требованию невырожденности матрицы РС 1С. Тогда, умножив первое уравнение в (1) на матрицу РС-1, получим

ГС^^у + ¿-РС-1^2^ + FC-1Gp = О, 11е

откуда можно найти р.

Предполагая, что система (1) снабжена сеточным аналогом кинетической энергии £(у) = (¿V, V), где Е = Е* > 0, для системы (1) определяются энергетическое и линейное критические числа Рейнольдса соответственно как точные нижние грани таких положительных Яе, при которых существуют решения, кинетическая энергия которых не убывает строго монотонно и просто не убывает при Ь —> оо.

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

Пусть

G = U

Rg О

F* = V

RF О

являются QR-разложениями, т.е. U и V - унитарные матрицы порядка nv, a Rq и Rp - верхние треугольные порядка пр и матрицы U и V разбиты на блоки следующим образом: U = [Щ, С/г], V = [Vi, V-2], где £4 и V/- - матрицы размера Пу х Пк и п\ = пр, П2 = nv — пр. Вводятся следующие обозначения: J\} UfJ^Vj, к = 1,2, Cij = U*CVj и рассматривается система

du ~dt

т(к)

ЯМи+^и, Re

(3)

где Н^ — (C22R~1)~1 Jyz R~l■> а R ~ фактор Холецкого в разложении Холецкого V2EV2 = R*R. Формулируется и доказывается следующая

Теорема 1.3.1 Если выполнено условие (2), то система (1) эквивалентна системе (3). Причем

v(i) = V2R~1u(t), p(t) = R3\C12R-\H^ + - J« - ^J^R-'uit),

а функционал энергии £(v) = ||u|||.

Таким образом, исследование устойчивости нулевого решения системы

(1) сводится к исследованию устойчивости нулевого решения системы (3). При

этом энергетическое Re^ и линейное Re¿ критические числа Рейнольдса, определенные в разделе 1.2, являются точными нижними гранями таких положительных Re, при которых существуют решения системы (3), вторая норма которых не убывает строго монотонно и просто не убывает при t —> оо соответственно. Для обоснования алгоритмов вычисления критических чисел Рейнольдса предполагается, что матрица - квадратная комплексная общего вида, а является дискретным аналогом сужения оператора Лапласа на подпространство соленои-дальных сеточных функций, и, следовательно, при правильной аппроксимации должна быть диссипативной: Н^ + Н^* < 0.

В разделе 1.4 показывается, что вычисление энергетического критического числа Рейнольдса может быть сведено к однократному вычислению максимального собственного значения эрмитового матричного пучка Н^ + Н^* + = 0 со знакоопределенной матрицей при спектральном параметре.

В разделе 1.5 предлагается алгоритм вычисления линейного критического числа Рейнольдса с заданной относительной точностью S с помощью обоснованных в разделе 1.6 модификаций стандартных процедур FZERO и FMIN. В частности, показывается, что Re¿ = 1/ni, где щ- максимальный положительный корень уравнения г (ц) = 0, а г (ß) - максимальная действительная часть собственных значений матрицы Н^ + ßH^2\ Если при Re¿ < Re < Reoo течение неустойчиво, то формулируется и доказывается

Теорема 1.5.1 Если процедуру FZERO( f(n), tol, a, b) применить к функции f(p) = г(р/Reco), с параметрами tol = <5/4, а = 1, b = RCoo/Re^, то Reí, будет вычисляться с заданной относительной точностью 5: |Rc¿, — Re^l/Re^ < 6, где Re^ = 1/Д¿ - вычисленное значение Re¿.

Здесь a, b - концы отрезка, где ищется корень.

Во второй главе рассматривается течение Пуазейля в трехмерном бесконечном в продольном направлении канале постоянного прямоугольного сечения с отношением Z > 1 длин сторон сечения. Исследуется зависимость Re¿(Z).

В разделе 2.2 для численного исследования зависимости Re¿(Z) описывается и обосновывается предложенный в работе [2] вариант обсуждавшейся в первой главе технологии.

В разделе 2.2.1 дается постановка задачи. В трехмерном бесконечном канале {(х,у,г) : —оо < х < оо, — Ь < у < Ь, —а < г < а} постоянного прямоугольного сечения

{(у, г): -Ъ<у <Ъ, -а < г < а}

(4)

с отношением длин сторон = а/Ъ > 1 рассматривается стационарное течение с давлением р = — тх (т > 0), вектором скорости V = (С/, 0,0) и профилем скорости 11 {у, г), удовлетворяющим уравнению Пуассона АН = ~т/(ир) в области (4) с нулевыми граничными условиями, где V и р - соответственно кинематическая вязкость и плотность жидкости, а Д = д2/ду2 + д2/дг2. Такое течение называют течением Пуазейля. Выбирая в качестве нормировки скорости ит ^ = П1 ах и (у, г) и нормировки длины Ь, число Рейнольдса для течения Пуазейля определяется как Не = итлхЬ/и и все уравнения далее записываются в переменных, обезразмеренных обычным образом.

Анализ линейной устойчивости течения Пуазейля сводится к анализу его устойчивости к бесконечно малым возмущениям вида

р'{х,у,г,г)

Р"(у,г)

\{ах—шЬ)

а > 0,

(5)

где а заданная вещественная константа, а амплитуды V , р и комплексная частота ш = сог + йог являются решением следующей проблемы собственных значений:

■ шч = + + Ср, Fv = 0,

Ке

(6)

где

' -\аи -иу -иг ' —1а

0 -'гаи 0 -д/ду

0 0 -юи -д/дг

о2 д2

[ю, д/ду, д/дг], № = сИад(Ь, Ь, Ь), Ь = ^ + — - а2

У~ ду'

~ дг '

и для упрощения обозначений V" и р" обозначены через V и р, соответственно.

Так как возмущение вида (5) нарастает со временем, если ш;>0и убывает, если Шг < 0, то для заданной величины И линейное критическое число

Рейнольдса Rez,(Z) течения Пуазейля, то есть граница устойчивости ко всем возмущениям вида (5), определяется следующим образом:

ReL(Z) = inf ReL(Z,a),

a>0

где ReL(Z,a) = inf{Re : maxImft(Z, a, Re) = 0, Re > 0}, a fi(Z,a,Re) означает спектр проблемы (6).

В разделе 2.2.2 описывается пространственная аппроксимация на гаус-совских сетках типа В: узлы Гаусса для давления и Гаусса-Лобатто для компонент скоростей. Аппроксимация осуществляется методом коллокаций и приводит к проблеме собственных значений вида:

- iwv = J(1)v + J(2)v + Gp, Fv = 0, (7)

где v - nv-компонентный вектор значений компонент скорости, р - пр-компонентный вектор значений давления, пр < nv, jW и J'2' - квадратные матрицы порядка nv, a G и F - прямоугольные матрицы размеров nv х пр и пр х nv, соответственно. Все эти матрицы не зависят от числа Рейнольдса.

Важно аппроксимировать проблему (6) так, чтобы не нарушить соотношение между конвективными и вязкими членами. Прежде всего, для этого необходимо сохранить на дискретном уровне наиболее важные свойства входящих в них операторов: симметричность и отрицательную определенность оператора Лапласа и соотношение div = —grad* между дивергенцией и градиентом. Аппроксимация методом коллокаций эти свойства, вообще говоря, не сохраняет, но выбранные сетки позволяют обосновать в разделе 2.2.3 преобразования проблемы (7), обеспечивающие их сохранение. В частности, формулируется и доказывается следующая теорема:

Теорема 2.2.1 Существуют такие невырожденные диагональные матрицы Kv и Кр, что справедливы следующие равенства:

KpFK~l = {KvGK~l)\ KvJ^K~l = {KvJ^K~l)\

причем

KyJ^K'1 < 0.

После замены переменных V := К„\, р := Крр, Р := КрРКу \ := С := КуСК'1, система (7) принимает вид:

- Ьлг = + 7(2)у + Р*р, Fv = 0, (8)

Ке

где - симметричная и отрицательно определенная матрица.

В разделе 2.2.4 для сокращения вычислительных затрат предлагается преобразование проблемы (8) к проблеме того же вида, но с чисто мнимой матрицей и вещественной матрицей Р. Четность профиля скорости V течения Пуазейля по у и г позволяет свести решение (6) к поиску решений, обладающих одной из четырех указанных в работе симметрий. Обсуждавшиеся выше аппроксимация и преобразования сохраняют это свойство решений. После учета симметрий, из (8) получаются четыре независимые проблемы собственных значений того же вида, но с матрицами примерно вчетверо меньшего порядка и далее рассматриваются решения только одной симметрии, которые являются наиболее неустойчивыми. После этого выполняется алгебраическая редукция, описанная в первой главе. При этом она применяется не к системе вида (7), как в первой главе, а к системе вида (8) с чисто мнимой матрицей симметричной и отрицательно определенной и вещественной F, что по сравнению с общей технологией сокращает число арифметических операций примерно в четыре раза, а объем памяти, необходимой для хранения матриц, - примерно в два раза. Далее описываются алгоритмы вычисления линейного критического числа Рейнольдса и построения соответствующих нейтральных кривых с заданной относительной точностью на основе предложенных и обоснованных в первой главе модификаций процедур Р2ЕЯО и БМШ.

В разделе 2.2.5 обсуждаются результаты численных экспериментов, в ходе которых для различных значений Z вычислялись линейные критические числа Рейнольдса (см. рисунок 1, кривая I) и исследовалась сходимость по шагу сетки.

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

направлению сетки для пространственной аппроксимации. Кроме того, она вы-

числяет линейное критическое число Рейнольдса с заданной относительной погрешностью. Все это позволило получить результаты с гарантированной и более высокой точностью.

В разделе 2.3 дается подробное теоретическое обоснование полученной зависимости линейного критического числа Рейнольдса от величины И.

В разделе 2.3.1 рассматривается плоское течение Пуазейля и обсуждается теорема Сквайра, утверждающая, что минимальное критическое число Рейнольдса для этого течения достигается на возмущениях с равным нулю поперечным волновым числом.

В разделе 2.3.2 рассматривается течение Пуазейля в канале прямоугольного сечения. Объясняется поведение 11сь(^) при больших значениях Z, путем сравнения профиля скорости рассматриваемого течения с соответствующим профилем скорости плоского течения Пуазейля.

В разделе 2.3.3 дается обоснование [1,2] величины такой что при Z < 2С линейное критическое число Рейнольдса бесконечно, а при Z > 2С - конечно. Для этого было рассмотрено некоторое модельное течение: течение Пуазейля в том же самом канале с нулевыми граничными условиями (условия прилипания) в направлении у, но с периодическими граничными условиями на боковых стенках для основного течения и возмущений. Единственным отличием модельного течения от плоского течения Пуазейля является дискретность поперечного волнового числа допустимых возмущений. То есть для него справедлива теорема Сквайра. В отличие от течения Пуазейля с нулевыми граничными условиями на боковых стенках, модельное течение допускает ненулевые возмущения, постоянные в поперечном направлении, которые, как известно, являются наиболее неустойчивыми и на которых достигается критическое число Рейнольдса. Исключив такие возмущения из множества допустимых, для модельного течения получена формула

Ввь{г)= Ы а 2/ Ие?5(а), (9)

а>п/г л/а2 - 7Г2/^2

где Ке£5(а) - критическое число Рейнольдса проблемы Орра-Зоммерфельда для плоского течения Пуазейля, т.е. функция от а, образующая левую ветвь нейтральной кривой плоского течения Пуазейля. Известно, что существует а®^ ~

1.097, такое что Ле^а) < оо при 0 < а < аи 11е£5(а) = со при а > а^ах-Это позволяет для модельного течения из (9) получить оценку 2с « 2.86.

Построенная по формуле (9) зависимость Кизображена на рис. 1 (кривая II). Близость кривых I и II позволяет утверждать, что в канале прямоугольного сечения механизм линейной устойчивости течения Пуазейля с условием прилипания на боковых стенках аналогичен механизму линейной устойчивости течения Пуазейля с условием периодичности на боковых стенках к возмущениям с неравным нулю поперечным волновым числом.

Рисунок 1: Зависимисть линейного критического числа Рейнольдса от 7, для исходного (кивая I) и модельного (кривая II) течений.

Третья глава посвящена решению частичной обычной и обобщенной проблем собственных значений с большими разреженными матрицами. Наиболее популярными методами решения этих задач являются метод Арноль-ди, несимметричный метод Ланцоша и метод Якоби-Дэвидсона. Существенным недостатком этих методов являются высокие требования к объему оперативной памяти, так как приближение к собственному вектору ищется в подпространстве, размерность которого растет на каждом шаге, а для нахождения инвариантного или понижающего подпространства приходится использовать логически более сложные и вычислительно-емкие блочные варианты этих методов.

Поэтому в настоящее время инженеры часто используют метод обратных итераций, который, хотя и обладает лишь линейной скоростью сходимости, требует существенно меньше оперативной памяти. Его недостатком является необходимость точного решения линейных систем с матрицами порядка п на каждой итерации, что может быть слишком затратным при достаточно большом п. В работах М.А. Freitag, А. Spence (2008), М. Robbe, М. Sadkane, А. Spence (2009) и F. Xue, Н.С. Elman (2011) был предложен очень экономичный вариант этого метода - приближенный метод обратных итераций с тюнингом, где упомянутые системы решались приближенно с увеличивающейся на каждой итерации точностью с помощью GMRES с правым предобусловливанием с тюнингом. В таком случае количество итераций GMRES не зависит от номера внешней итерации. Тем не менее, этот вариант метода обратных итераций сходится линейно. В данной главе для решения частичных проблем собственных значений описываются и обосновываются методы ньютоновского типа, являющиеся развитием варианта метода Ньютона для вычисления инвариантного подпространства матрицы, предложенного и обоснованного в работе С.К. Годунов, Ю.М. Нече-пуренко (2002). Начальное приближение для этих методов ищется с помощью описанных и обоснованных в этой главе модификаций приближенного метода обратных итераций с тюнингом.

Раздел 3.2 посвящен вычислению спектрального проектора, отвечающего изолированному подмножеству собственных значений большой разреженной матрицы А G С71ХП, ближайших к заданной точке комплексной плоскости. Без потери общности предполагается, что матрица А невырожденная и требуется найти спектральный проектор, отвечающий ее р <С п (с учетом кратности) минимальным по модулю собственным значениям. В противном случае предлагается рассматривать вместо исходной матрицы А матрицу А — ol, где а -некоторый сдвиг.

Пусть Xi и С С" - соответственно правое и левое инвариантные подпространства, отвечающие р минимальным по модулю собственным значениям матрицы А, а Х\ и Xi € Спхр - матрицы, столбцы которых образуют биорто-гональные базисы в этих подпространствах (так называемые биортогональные матрицы), т.е. = span(X;) и X£Xi = I. Тогда, спектральный проектор, который так же называется проектором Рисса или инвариантным проектором, отве-

чающий р минимальным по модулю собственным значениям матрицы А, может быть представлен в виде

V = (10)

Для вычисления матриц Х\ и Х2 предлагается использовать комбинацию двусторонних методов приближенных обратных итераций с тюнингом и Ньютона. Проектор (10), где Х\ и Хо биортогональные матрицы, является спектральным проектором, отвечающим некоторому подмножеству собственных значений матрицы А, в том, и только том случае, если он коммутирует с матрицей А. Поэтому в критериях останова обсуждаемых методов предлагается следить за величиной нормы коммутатора

Е = AV — РА, (11)

либо за нормами спектральных невязок

Rx = АХг - ХгА, R2 = А*Х2 - Х2Л*, (12)

где Л = XqAXj , которые связаны с Е равенством Е = R\X2 — Х\Щ. Показывается, что если матрицы Х\ и Х2 сбалансированно биортогональные, то есть ЦРЦг = ||Xi= ЦХ2Ц2, то малость нормы коммутатора влечет малость норм спектральных невязок и наоборот. Предложен быстрый алгоритм вычисления ЦЯЦг позволяющий не формировать явно матрицу Е, что привело бы к огромным затратам. При описании алгоритмов в диссертационной работе используются вспомогательные процедуры ort, biort, schür и reord. Первая выполняет орто-нормировку столбцов заданной матрицы, вторая - биортогонализацию столбцов двух заданных матриц, а schür и reord для заданной матрицы вычисляют разложение Шура с диагональными элементами формы Шура, расположенными в заданном порядке.

В разделе 3.2.1 описывается двусторонний метод приближенных обратных итераций с тюнингом. Стартуя с сбалансированно биортогональных матриц Xiß,X2iQ € Спхр, этот метод строит две последовательности матриц и Х2,ь к = 1,2,..., приближенно решая на каждом шаге блочные системы вида АХ\^+1 = Xi k и А*Х-2^+1 = Х2 ь с помощью GMRES с правым предобуслов-ливанием с тюнингом, и, выполняя сбалансированную биортогонализацию найденных матриц Х\ ¡t+i и Х2 jt+i с помощью процедур ort и biort. Перед решением

каждой блочной системы, с помощью schür и reord выполняется преобразование типа Шура-Релея-Ритца, описанного в работе G.W. Stewart (1976).

Алгоритм 1 (Двусторонний метод приближенных обратных итераций с тюнин-

Задать е > 0, р > О, г) > 0 и сбалансированно биортогональные матрицы

1. Вычислить Ак = XTlkAXik и спектральные невязки R\jc = АХitk — Xi,kAk, R.2.1; = A*X2tk — X2,k\%, и проверить сходимость коммутатора Ek =

— XijcB*2 k: если \\Ek\\2 < £, то положить Х:jOTti = X^k, Х2/т1 = X'2]k и остановиться.

2. Вычислить (Qk,Tk) = schür (Afc) и (Qk,Tk) = reord (Qk, Tk, <), положить X\ k := XijcQk, решить приближенно систему AX\ k+i = Xk относительно Xhk+i, так что - AXi,fc+1||2 < где 7i,fc = mm{/C>, r]\\Rhk\\2}

3. Вычислить (Qk,Tk) = reord {Qk,Tk, >), положить X2,k ■= X2¿Qk, решить приближенно систему A*X2]k+i = X2^k относительно X2jc+\, так что \\X2Jt - A*X2Mi||2 < 72,fc. где >y2,k = min{p, r]\\R2tk\\2}

4. Положить (X1)k+1, X2ik+1) = biort(ort(Xuk+1),ort(X2tk+1)).

В разделе 3.2.2 описывается двусторонний метод Ньютона. Пусть Х[ 6 Спхр (I = 1,2) - сбалансированно биортогональные матрицы, столбцы которых образуют базисы в соответственно правом и левом приближенных инвариантных подпространствах матрицы А, достаточно близких к искомым. Тогда V = Х\Х2 - приближенный спектральный проектор и будут малы по норме коммутатор Е и спектральные невязки i?i и R2, определенные в (11) и (12), соответственно. Для матриц

гом).

Xlfi,X2fl € Спхр. Для к = 0,1,...

Х1=Х[- Фь где Ф; £ Cnxp (I = 1,2)- решения системы уравнений

(13)

{I-Т){АФ1-Ф1А~ R1) = 0, РФ 1 = 0,

{I-VУ{A*Ф2-Ф2A*-R2)=0, V* Ф2 = 0,

(14)

и проектора V = X^X^Xi^X^ формулируется и доказывается следующая теорема, показывающая, что формулы (13), (14) описывают шаг ньютоновского уточнения исходных приближенных инвариантных подпространств span(Xi) и span(X2).

Теорема 3.2.1 Если т = Ц-БЦг достаточно мало, то матрица Х2Х\ невырождена и проектор V коммутирует с А с точностью т2, то есть Е = AV — VA = 0(г2).

Затем описывается алгоритм решения систем (14) на основе разложения Шура и GMRES с предобусловливанием. После этого дается формальное описание предложенного метода Ньютона.

Алгоритм 2 (Двусторонний метод Ньютона).

Задать е > 0, 5 > 0 и сбалансированно биортогоналъные матрицы -АГ^у, X2ß G Спхр. Для к = 0,1,...

1. Вычислить Л^ = Х*1кАХ\к и спектральные невязки R\jc = AXi^—

Raj; = A*X'i,k — X2и проверить сходимость коммутатора Ei- = Ri^ - XiikR*l k: если \\Ек\[2 < е. то положить Xit0ut = Х^к, X2:OUt = X-2j; и остановиться.

2. Вычислить (Qk,Tk) = schür (Л^.) и () = reord (Qk, Тк, <), положить Ri,k '■= RijcQb решить приближенно систему (I — "Рк){АФ\^ — Ф1 ,кТк) = (I — Vk)R\,k, *Pk<&\,к = 0 с Vk = Xi^X^k относительно Ф^. и положить

$1,* := ФhkQl

3. Вычислить (Qk,Tk) = reord (Qk,Tk, >), ПОЛОЭЮиШЬ Н-2,к •— R-2,kQk> pCtUUtflb приближенно систему (I-Т>к)*(А*Ф2,к~$2,кТ£) = {I~Vk)*Ri,b ,к = О с Vk = Х\,кХ2к относительно Ф2,(с и положить Фг^ :=

4. Выполнить ньютоновский шаг и сбалансированную биортогонализацию:

(XiMl, X2,k+\) = biort(ort(Xlife - Фг,й),оИ(Х2,ь - Ф2,к)).

В разделе 3.2.3 обсуждаются результаты численных экспериментов с дискретным аналогом неэрмитового эллиптического оператора, которые показали высокую эффективность предложенного подхода.

В разделе 3.3 обсуждается вычисление понижающего подпространства X С С", отвечающего изолированному подмножеству из р конечных собственных значений регулярного матричного пучка вида А — ХВ с большими разреженными матрицами А и В € С"*", причем В может быть вырожденной. Регулярность пучка означает, что существует по крайней мере одно Ао такое, что — Ао В) ф 0. Конечными собственными значениями пучка называют

решения уравнения <Зе<;(Л — ХВ) — 0. Без потери общности, предполагается, что матрица А невырожденная и искомое понижающее подпространство отвечает р конечным минимальным по модулю собственным значениям. В противном случае предлагается заменить А на А —а В, где сг - точка вблизи р рассматриваемых собственных значений такая, что матрица А — аВ невырождена.

Метод, предложенный в [6,9] и описываемый в этом разделе, может рассматриваться как обобщение алгоритма, предложенного в работах С.К. Годунов, Ю.М. Нечепуренко (2002) и О. Е1 КЬоигу, Уи.М. ЫесЬеригепко, М. 8ас1капе (2014), на случай частичной обобщенной проблемы собственных значений. Данный алгоритм является комбинацией приближенного метода обратных итераций с тюнингом и метода Ньютона.

В разделе 3.3.1 описывается обобщение метода приближенных обратных итераций с тюнингом.

Алгоритм 3 (Приближенные обратные итерации с тюнингом). Задать £ > 0, 5г > 0, ё2 > 0 и Х0 е Спхр такую, что Х£Х0 = I, Для к = 0,1,...

1. Вычислить Ук = АХк, IV,. = оП(ВХк) иЬк = \¥£Ук,

2. Положить Як = Ук — \¥кЬк и проверить сходимость: если И-Й^Иг < £, то положить Хои( = Хк и выйти.

3. Решить относительно Хк+\ систему АХк+\ = ВХк приближенно, так что

||ВХк - АХк+11|2 <рк = тшф.йИДкНг}- (15)

4. Положить Хк+\ := ог1(Х^+1).

Полученное приближенное понижающее подпространство используется

как начальное приближение для описанного и обоснованного в разделах 3.3.2

- 3.3.4 приближенного метода Ньютона для вычисления понижающих подпространств. Предложенный метод Ньютона на каждой итерации требует решения обобщенного уравнения Сильвестра, которое решается с помощью алгоритма аналогичного упомянутому выше алгоритму на основе разложения Шура и GMR.ES с предобусловливанием.

Алгоритм 4 (Приближенный метод Ньютона). Задать е > О, 5 > 0 и Х0 е Спхр такую, что ХЦХ0 = I, вычислить Уо = АХ0, Щ = ВХо, \¥о = 011([/о), Для к = 0,1,...

1. Вычислить Як = Ук — \¥к(\У^Ук),

Проверить сходимость: если ЦДьЦг < то выйти.

2. Вычислить (Ук,Бк) = ог^У*;) и положить Хк '■= ХкЯ^1, Ук := Ук.

3. Вычислить {(¿к, Тк) = 5сИиг(У/^..), (С}к, Тк) = геогс) (С}к, Тк, >) и положить Ск = Х^Г1.

4. Положить 11к := икС}к, Ук := УкС}к.

5. Вычислить Д^- = Ук — 11кСк = 01Х(А*Ук).

6. Решить относительно Фк приближенно систему (/ — УкУк)(АФк — ВФкСк) = (I - УкУк*)Ак, 2*кФк = 0.

7. Вычислить Хк+1 - оП(Хк - Фк), Ук+\ = АХк+ъ ик+1 = ВХк+ъ и \Ук+1 = оП{ик+1).

В разделе 3.3.5 кратко описывается обобщенная проблема собственных значений для системы дифференциальных уравнений в частных производных, возникающая при анализе гидродинамической устойчивости течения Пуазейля в канале постоянного прямоугольного сечения. Эта задача подробно обсуждалась во второй главе диссертационной работы. В главе 2 пространственная аппроксимация выполнялась методом коллокаций и приводила к обобщенной проблеме собственных значений с плотными матрицами не очень большого порядка. Для сокращения вычислительных затрат к полученной обобщенной проблеме применялась алгебраическая редукция, сводящая ее к обычной проблеме собственных

значений с плотными матрицами вдвое меньшего порядка. Далее с помощью QR-алгоритма, вычислялись все конечные собственные значения и выбирались те, что были необходимы.

В разделе 3.3.6 рассматривается та же система дифференциальных уравнений, но используется конечно-разностная аппроксимация второго порядка точности на равномерных разнесенных сетках типа С, предложенных В.И. Лебедевым (1964), Эта аппроксимация значительно менее точная, чем метод колло-каций и приводит к частичной обобщенной проблеме собственных значений с большими разреженными матрицами, которая используется в качестве тестовой для предложенного метода Ньютона. Так как этот метод вычисляет только часть собственных значений, для того чтобы продемонстрировать, что найдены именно те собственные значения, что нужны, используются собственные значения, полученные с помощью метода коллокаций и QR-алгоритма. Результаты экспериментов показали высокую эффективность предложенной комбинации алгоритмов (3) и (4).

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

2. Для исследования устойчивости течений в бесконечном канале постоянного прямоугольного сечения предложен, обоснован и реализован специальный вариант технологии, позволяющий сократить требуемый объем оперативной памяти в два раза, а число арифметических операций - примерно в

четыре раза. С помощью этого варианта исследована зависимость линей-

ного критического числа Рейнольдса течения Пуазейля от отношения длин сторон сечения. Выполнены расчеты, уточняющие известные результаты.

3. Для течения Пуазейля в бесконечном канале постоянного прямоугольного сечения впервые дано теоретическое обоснование зависимости линейного критического числа Рейнольдса от величины отношения длин сторон сечения.

4. Предложен и обоснован двусторонний метод Ньютона для вычисления спектрального проектора, отвечающего заданной группе собственных значений большой разреженной матрицы. Для поиска начального приближения разработан двусторонний метод приближенных обратных итераций с тюнингом.

5. Предложен и обоснован метод Ньютона для вычисления понижающего подпространства, отвечающего заданному изолированному подмножеству собственных значений регулярного матричного пучка с большими разреженными матрицами. Для нахождения начального приближения разработан вариант приближенного метода обратных итераций с тюнингом.

Публикации автора по теме диссертации

1. Демьянко К.В., Нечепуренко Ю.М. О зависимости линейной устойчивости течений Пуазейля в прямоугольном канале от отношения длин сторон сечения // Доклады АН. 2011. Т. 440, № 5. С. 618-620.

2. Demyanko K.V., Nechepurenko Yu.M. Linear stability analysis of Poiseuille flow in a rectangular duct // Russ. J. Numer. Anal. Math. Modelling. 2013. T. 28, № 2. C. 125-148.

3. Демьянко K.B., Нечепуренко Ю.М. Двусторонний метод Ньютона для вычисления спектральных проекторов // Вычислительные методы и программирование. 2014. Т. 15. С. 121-129.

4. Нечепуренко Ю.М., Демьянко К.В. О влиянии отношения сторон на устойчивость течений в бесконечных каналах прямоугольного сечения // «Модели

и методы аэродинамики». Материалы 11-й международной школы-семинара. М.: МЦНМО, 2011. С. 149-150.

5. Демьянко К.В., Нечепуренко Ю.М. Двусторонний метод Ньютона для вычисления спектральных проекторов // Тезисы докладов XX Всероссийской конференции и Молодежной школы-конференции «Теоретические основы и конструирование численных алгоритмов решения задач математической физики», посвященной памяти К.И. Бабенко (Дюрсо, 15-20 сентября, 2014 г.). М: Институт прикладной математики им. М.В. Келдыша, 2014. С. 55.

6. Демьянко К.В., Нечепуренко Ю.М., Садкан М. Метод Ньютона для вычисления понижающих подпространств регулярных матричных пучков // Современные проблемы вычислительной математики и математической физики: Международная конференция, Москва, МГУ имени М.В. Ломоносова, 16-17 июня 2014 г.: Тезисы докладов. М.: МАКС Пресс, 2014. С. 42^13.

7. Демьянко К.В. О зависимости линейной устойчивости течений Пуазейля в прямоугольном канале от отношения длин сторон сечения // Труды 54-й научной конференции МФТИ. М.: МФТИ, 2011. С. 72-73.

8. Демьянко К.В. Устойчивость течения Пуазейля в канале прямоугольного сечения // Труды 55-й научной конференции МФТИ. М.: МФТИ, 2012. С. 13-14.

9. Демьянко К.В. Метод Ньютона для решения частичной обобщенной проблемы собственных значений // Труды 56-й научной конференции МФТИ. М.: МФТИ, 2013. С. 136-137.

Заказ № 64-а/12/2014 Подписано в печать 18.12.2014 Тираж 100 экз. Усл. п.л. 1,0

ООО "Цифровичок", тел. (495) 649-83-30 www.cfr.ru; е-таИ:гак@с/г. ги