Численное исследование устойчивости некоторых физико-технических систем Савенкова Надежда Петровна

Содержание к диссертации

Введение

Глава 1. Исследование устойчивости нестационарных процессов, описываемых дифференциальными операторами единенным вхождением собственной функции 28

1,1. Вычисление граничных точек комплексного спектра линейной незрмитовой задачи на собственные значения 28

1.L1. Некоторые асимптотические свойства системы дифференциальных уравнении с .постоянными коэффициентами 28

1,1*-2. Алгоритм нахождения граничных точек комплексного спектра 30

1.1.j' Модификации алгоритма RBS 33

1.1.4 Численный алгоритм построения кривых устойчивости на плоскости параметров 34

1.2 Нахождение границ спектра обобщенной задачи на собственные значения , 47

1.2.1. Алгоритм нахождения границ спектра обобщенной алгебраической задачи на собственные значения, 47

1.2.2. Численное исследование устойчивости течений Пуазеля и Блазиусав рамках краевой задачи на собственные значення для уравнения Орра-Зоммерфельда 48

1.3 Численное исследование спектра задачи на собственные значения с нелинейным вхождением спектрального параметра 51

1.3.1 Алгоритм нахождения границ спектра симметричной положительно определённой задачи на собственные значения с нелинейным вхождением спектрального параметра 51

1.3.2 Качественное исследование устойчивости в заданном интервале симметричной задачи на собственные значения с нелинейным вхождением спектрального параметра 54

1.3.3 Локализация спектра задач на собственные значения с нелинейным вхождением спектрального параметра 5S

1.3.4 Алгоритм решения нелинейной задачи 61

1.4 Численный метод нахождения кусочно-непрерывного спектра 67

1.4.1. Двумерная задача Шредингера с периодическим потенциалом 67

Ы.2. Модельная задача для уравнения Шредингера с периодическим потенциалом 70

1.4.3. Анализ численных расчётов 73

Приложение к главе 1 76

Глава 2. Устойчивость нестационарных процессов, сводящихся к дифференциаль ным задачам с нелинейным вхождением собственной функции 89

2.1 Построение нейтральных кривых устойчивости на двумерной плоскости параметров 89

2.1 1. Алгоритм нахождения собственных значений и соответствующих им функций нелинейной задачи на собственные значения 89

2.2.Численное исследование нелинейной задачи на собственные значении, имеющей аналитическое решение 92

2.2.[.Модельная задача 92

2,2.2.Построение разностной схемы 93

2.23.Первый способ аппроксимации нелинейного оператора А2 94

2.2.4.Итерационный метод решения разностного уравнения для несимметричной матрицы 96

2.2.5.Результаты тестовых расчетов для несимметричной матрицы , 97

2.2.6.Второй способ аппроксимации нелинейного оператора А: 100

2.2.7.Итерационный метод решения разностного уравнения для симметричной матрицы ...102

2.2.8.Результаты и их обсуждение для симметричной матрицы 103

2.2.9,Общие замечания к расчетам 105

2.3.Математическое моделирование стационарного режима распределении фемтисекундных импульсов в фотонном кристалле 108

2.3,1 .Задача само воздействия светового импульса в фотонном кристалле 108

2.3.2,Построение разностной схемы 110

2.3.3.Итерационный метод решения разностного уравнения 111

2.3,4.Лнализ численных расчетов 111

2.3.5 Лислснное моделирование солитонов сложной временной формы импульсов фемтосекундной длительности 115

2 3.б.Условие действительности собственного значения для решения типа солитона 116

2.3.7.Построение разностной схемы для задачи распространения фемтосекундного импульса в среде с кубичной нелинейностью 117

2.3.8.Итерационный метод решения разностного уравнения 118

23.9.Исследование зависимости формы солитона от параметров у и а на примере первого собственного вектора '. 120

2.3.1О.Числешые расчеты различных собственных значений и соответствующих им собственных векторов 121

2,3 Л1.Численные расчёты для сеток большого размера 125

Приложение к главе 2 127

Глава 3 Исследование устойчивости нестационарных процессов,оппсываемых оператором Навье-Стокса 162

ЗЛ.Система дифференциальных уравнений 162

ЗЛ Л.Постановка задачи 162

ЗЛ.2, Получение двумерной модели, 163

ЗЛ.З. Упрощенная модель 174

3.2.Численный метод решении задачи.., . 176

3.2.1,Метод расщепления по физическим процессам 176

3,2.2-Численное решение отдельных этапов 179

З.З.Анализ численных экспериментов 187

1.3 Л «Модельные расчёты 187

3.3.2Лисленное моделирование аварий с растеканием тяжёлых газов и жидкостей 191

33.3. Растекание жидкости по неровной поверхности 197

3.3.4. Расползание газа по неровной поверхности 204

3.3.5, Моделирование аварии газопровода 209

З.З.б.Мехакикакратерообразования , 214

Глава 4. Математическое моделирование МГД неустойчивостсй в алюминиевом электролизере ,220

4.1 Математическая постановка задачи 220

4.1.1.Введение основных обозначений 220

4.1.2.Вывод уравнения давления в среднем слое 225

4.1.3, Получение уравнения электромагнитной индукции, учитывающего реальную подводку тока к электролизеру 228

4.1 А Моделирование функции помехи , 230

4.1.5* Полная постановка задачи,.,., 233

4.1.6. Начальные и граничные условия 235

4.2. Численный метод решения 237

4.2.1. Расщепление по физическим процессам , .-,.237

4.2.2. Разностный метод решения задачи * 241

4.2.3 Условие устойчивости 253

4.3. Анализ численных расчетов 257

4.3. L Описание возможностей комплекса вычислительных программ 257

4.3.2. Тестовые расчеты: одномерный плоский случай 258

4.3.3. Тестовые расчеты: двумерный случай 265

4.3.4. Расчеты по данным реальной электролизной ванны 297

Заключение 305

Литература

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

Введение к работе

В диссертации изучаются различные методы исследования устойчивости нестационарных нелинейных физических процессов. Рассматриваются три вида физических процессов: 1) процессы, в которых возможно проведение линеаризации относительно известного невозмущенного решения ; 2) процессы, в которых линеаризация недопустима, т.к. ведёт к искажению сущности физического явления; 3) существенно нелинейные физические процессы, математическая постановка которых базируется на уравнениях адвекции и магнитной гидродинамики.

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

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

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

В диссертации предлагаются методы исследования устойчивости перечисленных выше трёх типов физических процессов- К методам предъявляются следующие требования:

1. Предсказать динамику протекания процесса на основе исследования спектра соответствующей стационарной задачи в случаях возникновения : -линейной задачи на собственные значения типа Штурма-Лиувилля

7 -задачи с нелинейным вхождением собственного числа и линейным вхождением собственной функции

-задачи с линейным вхождением собственного числа и нелинейным вхождением собственной функции -задачи с кусочно-непрерывным спектром для оператора Шредингера.

2 .Исследовать динамику протекания нестационарного процесса , предварительно проведя математическое моделирование процесса на основе использования многомерного оператора Навье-Стокса для :

-задачи растекания тяжёлых жидкостей и газов по орографически неоднородной поверхности -задачи электролиза алюминия.

При исследовании устойчивости физических процессов в механике жидкости и газа [46-50], технических конструкций [51-52], обычно приходят к решению дифференциальных задач на собственные значения. При этом требуется найти одно или несколько собственных чисел и собственных функций, характеризующих устойчивость или неустойчивость исследуемого процесса. Если дифференциальная задача на собственные значения решается разностным методом [53], то возникает необходимость исследования спектра алгебраической задачи на собственные значения для разреженных матриц большого порядка. Причём часто требуется найти собственные числа, имеющие максимальную действительную или мнимую часть, и собственные векторы, отвечающие этим собственным числам, поскольку такие собственные числа определяют инкремент роста или декремент затухания процесса. Аналогичные собственные числа требуются и при рассмотрении чисто математических проблем ( например^ в теории разностных схем и при построении вычислительных алгоритмов) [53-54].

Хорошо известный класс степенных методов [55-57] в случае комплексных собственных чисел и отсутствия начальных приближений к искомым собственным числам позволяет найти только экстремальные по модулю собственные числа. Поэтому в случаях, когда спектр задачи комплексный, для определения инкремента (декремента) неустойчивости (устойчивости) приходится решать полную проблему собственных значений [58],что является не всегда осуществимо в силу ограниченности памяти ЭВМ и длительности времени расчёта. В работах [58-60] предложен метод, позволяющий вычислить несколько собственных чисел с максимальной или минимальной действительными частями большой разреженной несимметричной матрицы, по при этом метод требует наличия априорной информации о расположении спектра. Для численного исследования свойств решений некоторых задач, например, при построении нейтральных

кривых устойчивости, также неэффективно искать весь спектр матрицы, часто достаточно найти только одно собственное число и соответствующий ему собственный вектор.

Наиболее используемыми методами решения частичной проблемы собственных значений являются следующие методы: метод бисекции, степенной метод, метод обратных итераций и его различные модификации [3,53,55,79]» Однако при отсутствии начального приближения эти методы также не позволяют вычислить инкремент (декремент) процесса.

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

В диссертации демонстрируется применение предложенного алгоритма к

численному исследованию устойчивости течений Пуазейля и Блазиуса в рамках краевой задачи на собственные значения для уравнения Орра-Зоммерфельда. Известны несколько подходов к решению задач такого рода. Один из них изложен в [62], где дифференциальная задача Орра- Зоммерфельда с помощью конечно-разностного метода сводилась к алгебраической проблеме собственных значений, для решения которой использовался один из вариантов прогонки, либо решалась полная проблема собственных значений [9],

Другой подход в решении задач гидродинамической устойчивости заключается в непосредственном численном интегрировании обыкновенных днфференциальных уравнений возмущённого движения [63-65], Обычно при этом используется метод стрельбы для решения краевой задачи. Основной недостаток метода состоит в необходимости наличия априорной информации о распределении спектра, а эффективность метода определяется выбором начальных приближений искомых параметров, что само по себе является самостоятельной проблемой.

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

достаточно хорошо изучено [55-57,70] и соответствующая литература представлена

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

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

свести к линейной задаче на собственные значения, как правило требуется определить:

1.Область локализации на комплексной плоскости всех или части собственных значений (

например, собственных значений с ReX>0 ).

2.Количество собственных значений в заданной части комплексной плоскости( например,

в правой полуплоскости).

З.Величины всех или части собственных значений с максимальной действительной

частью,

С математической точки зрения в задачах устойчивости наиболее трудным является

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

надо знать, является ли оно искомым, т.е, тем, которое отвечает за устойчивость или

неустойчивость процесса и должно быть крайним в области локализации ( например,

самым правым или самым верхним). Поэтому чаще всего спектральную задачу решают в

уже заданной области локализации спектра, которая, определяется из физических

соображений. Однако это не всегда удаётся сделать,

В диссертации предлагается и обосновывается подход, позволяющий локализовать спектр для широкого класса алгебраических задач на собственные значения с нелинейным вхождением спектрального параметра ( к нему, например, относятся задачи, которые получаются из систем обыкновенных дифференциальных уравнений с запаздыванием[66-67], дробно-линейные задачи на собственные значения [68]), Заметим, что в линейном случае локализовать спектр можно с помощью теоремы Гершгорина [71],

Определить число собственных значений нелинейной спектральной задачи, находящихся в заданной области комплексной плоскости, можно с помощью известных алгоритмов, основанных иа теоремах о нулях аналитических функций [72-73]. Зная, где на комплексной плоскости находятся все или часть собственных значений и сколько их, можно приступить к их нахождению. Тут в основном используются итерационные метод, которые в большинстве случаев сводятся к вычислению нулей некоторой функции f(X). Например, в работах [18,74] в качестве /(Л.) берётся младший диагональный элемент треугольного разложения матрицы В(Л)-А(Л)-Я Е в работе [75] f(X)-detB(X) и f(A.)—P_n(A), в работах [72-73], где Р_п(Л) — полином степени п от А.

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

обратных итераций на нелинейную проблему собственных значений [76], итерационный метод последовательных линейных задач [77], а также метод сеток [78].

Т.о. в литературе достаточно широко представлены численные алгоритмы для решения линейных по собственному вектору задач с нелинейностью по собственному значению [79-85], и при этом практически не рассматриваются задачи, нелинейные по вектору. Единого механизма решения таких задач не существует. Большинство работ в этой области посвящено аналитическому исследованию свойств определенных достаточно простых типов задач в некоторых пространствах. Например, в работах [86-87] исследуется задача

ф(и)=ЛЧ"(и)^иєФ^ч(а),

где ф(ы)=— J|Va|²<& - — \u²dx, Ч*(и)= Jg(x,m(*))g!a: определены в пространстве
2 q 2 _а п

H_Q[l). Автор доказывает, что в зависимости от того, является се положительным или

отрицательным, задача имеет одно или два решения.

В конце 70-х годов группа немецких математиков (Kurt, Lackner, Amann и другие) предложили численный метод решения нелинейных задач, основанный на методе Пикара [78,79,88-95]. Рассматривались задачи типа

f(x_fU(x))=XL{U), х<=П,

U{x)=0, хедП,

где L(u) = - ^da_a(x)d_kU(xy} + a{x)u(x)_t ПсЛ", X -действительное положительное

/.АН

собственное значение* при условии, что функции/ <х_1к и граница dd достаточно гладкие, а матрица А = [а^.) симметричная и положительно определенная. Такие задачи

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

Шаг\ (Итерационный метод Пикара) *Lip**^l)=f(x_BU"{x% хєП,

U^n+[(x)=O_t хєП.

Шаг2 (Нормализация)
тгл-4-1 11/7^л*чІ /7^л+|

Таким образом, авторы рассматривали только те задачи, у которых нелинейная часть уравнения является функцией от t/, производные U* не входят в нелинейную правую часть уравнения, И при этом в работе не предложен алгоритм для нахождения

собственного значения, находится только собственный вектор и не определяется, какой моде соответствует этот вектор. Численные результаты для такого типа задач получаются с достаточной точностью, но четкого доказательства сходимости нет. В работе [114] рассматриваются дифференциальные задачи тина -Au=\f'(u\ иєП_г

ы=о, и еда,

где и = и(х) — действительная функция, определенная в области QciR , X —

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

Е(и):=— jjVw| dx -> тіп на [f(u)dx-y

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

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

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

Для разрешения многих вопросов относительно механизма процессов переноса вещества и энергии в тяжелых газах были поставлены специальные эксперименты. Наиболее полная серия опытов по рассеянию тяжелого газа была выполнена в Англии в 1982-1984 гг. на полигоне острова Торни-Айленд [30], Исследовалось рассеяние больших объемов

12 предварительно подогретой до температуры окружающего воздуха смеси фреона-12 и азота. Контейнер с газом объемом 2000 м представлял собой разборный пластиковый цилиндр высотой 13м, Его конструкция предусматривала возможность быстрого раскрытия стенок и практически мгновенного истечения смеси. Эксперименты финансировались группой промышленных компаний топливно-энергетического направления различных стран мира.

Так как проведение натурных экспериментов в условиях реальных предприятий в большинстве случаев является невозможным, то основным инструментом исследования процесса рассеяния облаков тяжелых газов, испускаемых промышленным объектом, является метод математического моделирования. Он включает в себя анализ физики всех фаз явления, построение на этой основе математической модели общего процесса, разработку численных методов и алгоритма решения задачи, написание компьютерной программы и проведение вычислительных экспериментов. На сегодняшний день задача описания рассеяния облака тяжелого газа в условиях термической и орографической неодиородностей, как отмечается в обзорах [96], [97], является одной из наиболее актуальных в промышленной безопасности.

При численной реализации трехмерных моделей рассеяния тяжелых газов используются конечно-разностные методы, реже - метод конечных элементов [98], [99]. Для численного моделирования пространственных течении горячих масс газа (в качестве математической модели выбрана полная система уравнений Навье - Стокса для сжимаемого теплопроводного газа) использовалась конечно-разностная методика, основанная на явпой 3-х шаговой схеме расщепления по физическим процессам [100]. Иногда используется описание диффузии в переменных Лагранжа наряду с описанием конвекции в переменных Эйлера [99]. В работе [101] на основе явного метода Годунова построена схема для численного моделирования 3-х мерных нестационарных уравнений Навье - Стокса, Конвективный член дискретизируется явно, а вязкий - неявно.

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

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

13 ограниченность методов решения системы связана с невозможностью развития их па случай исследования пространственных потоков вязкой жидкости, а также течений сжимаемого газа, объясняет возросший в последнее время интерес к численному решению уравнений Навье -Стокса, записанных в естественных переменных.

По данным работы [96] к 1987 г. в мире существовало более сорока специализированных программ расчета аварийного рассеяния в трехмерном пространстве. Основной трудностью при численном моделировании процесса распространения примеси в атмосфере является ограниченная емкость памяти и скорость современных вычислительных машин.

Отметим еще один важный момент. Одной из главных задач вычислительных методов расчета выброса в окружающую среду вредных химических отходов предприятий является расчет полей температур и концентраций, вызываемых такими выбросами в атмосфере и поверхностных водах [103], Распространение загрязняющих веществ в жидких средах определяется двумя процессами: конвективным переносом вследствие осредненного движения среды и диффузией за счот турбулентности. Поэтому математическая модель должна правильно описывать как поле средних скоростей, так и характеристики турбулентной диффузии. Система точных уравнений, описывающих во времени все детали эволюции поля скорости и скалярных полей в практической задаче ис может быть решена с помощью современных вычислительных средств. В настоящее время существует единственный экономически оправданный выход: решать уравнения осредненного движения, которыми определяется распределение осреднепных по Бремени величин [103]. При этом время осреднения должно быть много больше временного масштаба турбулентности, но много меньше временного масштаба осредненного течения (например, суточного цикла в пограничном слое атмосферы). Обычно только средние величины и имеют практический смысл. Уравнения осредненного движения содержат члены турбулентного переноса. Для замыкания системы уравнений эти члены должны быть аппроксимированы с помощью определенной модели турбулентности. В некоторых случаях вполне достаточно лишь приближенного описания характера турбулентности. Так, в задачах о больших водных массах, значения турбулентной вязкости и коэффициент турбулентной диффузии обычно принимают постоянными. Более сложные модели в таких задачах себя не оправдывают из-за значительной неопределенности в задании граничных условий и погрешности в численном решении. Для других задач, в частности, при расчетах полей вблизи источников, требуются уже более точные модели. Сложность модели турбулентности зависит от условий конкретной задачи.

14 Во многих случаях * в потоках со свободными поверхностями средние значения характеристик потока слабо меняются в вертикальном направлении, так что изменение этих величин в горизонтальном направлении можно определить, решая двумерные уравнения среднего течения для осредыенных по глубине значений [103], Уравнения среднего течения для осредненных по глубине величин можно вывести, интегрируя 3-х мерные уравнения по глубине в предположении, что распределение давления является гидростатическим. Такая модель была построена для жидкости, однако газ автором модели не рассматривался.

К недостаткам известных численных моделей рассеяния тяжелых примесей следует отнести излишнюю идеализацию граничных условий на поверхности земли [96]. Обычно она полагалась гладкой, что не соответствует наиболее важным случаям реальной действительности. Это говорит лишь о начальной стадии исследования явления с помощью современных математических и вычислительных средств.

Ещё более сложной является задача изучения магнитно-гидродииамических неустойчивостей в алюминиевом электролизёре. Математическое моделирование проводится здесь сразу в двух средах, металле и электролите, и имеет две основные цели: I .Получение основных полей в стационарном состоянии. 2.Изучение механизмов, вызывающих волнообразование на поверхности металла.

МГД-явления изучаются комплексно, включая натурные измерения, физическое и математическое моделирование [104-113] и др. Измерение циркуляции электролита и алюминия на промышленных электролизёрах методом радиоактивных изотопов впервые было предпринято в работе [104], позднее, также методом радиоактивных изотопов в [105] были получены схемы движения расплавов, подобные [104,106], Траектории движения металла, скорость в металле были измерены методом растворения железных стержней [107],

Преимущества математического моделирования заключаются в быстроте, многовариантности и экономичности проведения исследований, В настоящее время можно выделить три вида моделей для описания физических полей в алюминиевом электролизёре. Первый из них, например, [108] основан на двумерных уравнениях Навье -Стокса и к - є модели турбулентности; учитывается только движение в горизонтальной плоскости, не учитывается вертикальный перенос импульса, т.е. трение относительно узких слоев жидкости о дно ванны, нижнюю поверхность анода и между собой. Второй подход - это рассмотрение движения металла, электролита, распределения температуры в вертикальном разрезе, представляющий собой систему уравнений тепло-электропереноса* Максвелла и Навье - Стокса [109] Этот подход не позволяет получить распределение

температур и линии течения в пленарной плоскости, тем не менее, конвекция учитывается с помощью эффективного коэффициента теплопроводности. Третья, наиболее распространённая в настоящее время модель Моро-Эванса [110], Одна из главных идей в формулировке этой модели состоит в том, что основной интерес представляет описание течения средних слоев жидкостей. Как следствие этого в уравнении движения в горизонтальной плоскости учитывается трение слоев друг о друга. Следующий шаг - это оценка значений вклада каждого члена в уравнение для моментов движения. В результате разность между электромагнитным силами и градиентом давления уравновешивается силами трения, которые предполагаются пропорциональными скорости. Для определения граничных условий в электролите включаются большие каналы. Позднее эта модель получила своё развитие в работе [111], где предложен метод для расчета стационарного электромагнитного поля в предположении, что распределение электромагнитных величин при продольном расположении электролизеров в серии двухмерно и не зависит от координаты х вдоль длинной стороны ванны. Для расчета магнитного поля от соседнего ряда электролизеров последний представляется как линия с током, равным по модулю току серии. Магнитное поле от соседнего ряда электролизеров и от ошиновки считается по закону Био-Савара-Лапласа для отрезка с током. Для нахождения магнитного поля от токов, текущих внутри электролизера из всех возможных решений уравнения V В=0 выбираются те, которые соответствую простейшему распределению плотности тока в жидкости, как линейной функции координат. Это позволяет получить решение гидродинамических уравнений в виде аналитических выражений. В [112] авторы описывают методику расчёта распределения тока и магнитных полей для конкретной конструкции электролизёра. Далее в ряде работ, например [45,113] также используется модель Моро - Эваиса. Для расчёта магнитного поля используются замеры магнитного поля по периметру ванны. Наблюдается совпадение рассчитанных траекторий движения металла с экспериментальными, полученными в [107].

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

Приведём краткий обзор содержания диссертации.

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

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

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

В 1.1 рассматривается задача па собственные значения с линейным вхождением собственной функции V(x) и собственного числа h которая возникает при исследовании устойчивости бриллиюэпового течения плазмы. Задача формулируется в виде дифференциального уравнения четвёртого порядка:

ch х dx \ dx at x )

К(фГ()=0_?

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

17 спектре существует хотя бы одно положительное собственное значение). В основе метода находится разработанный в [1-2 ] алгоритм нахождения границ спектра на комплексной плоскости (RBS-метод) для стандартной алгебраической задачи на собственные значения. При этом нет необходимости решать полную проблему собственных значений, что существенно экономит время расчётов. Кроме того, экономия оперативной памяти ЭВМ, присущая данному методу , позволяет брать достаточно мелкие сетки расчета для достаточно больших значений длины отрезка расчёта L .

В 1.2 исследуется обобщённая задача на собственные значения с линейным вхождением собственных функций <р(у) и собственного числа L Задача возникает при исследовании течений Пуазейля и Блазиуса в рамках краевой задачи на собственные значения длл уравнения Орра-Зоммерфельда и в исходной постановке описывается дифференциальным уравнением четвёртого порядка ( см. [ 8 ])

d(y,A) = 2a²+AR + iaRU(y)

д(у_гЯ) = a⁴ + iaRU\y) + Ra²X + ia*RU(y)

Здесь Х= -іас, і -мнимая единица, а -волновое число, R -число Рейнольдса невозмущенного течения, с -фазовая скорость, U(y) -профиль безразмерной скорости. Если 1т с < 0» то происходит затухание, если 1т с > 0, то имеет место нарастание колебаний.

Граничные условия для течения в канале записываются в виде

Граничные условия для течения в пограничном слое имеют вид

В работе предлагается обобщённый метод RBS для построения в плоскости параметров (a , R) кривой устойчивости. Метод эффекгивен и экономичен, поскольку не требует вычисления в каждой точке плоскости параметров всего спектра обобщённой матричной задачи на собственные значения вида Аи^ХВи. Благодаря использованию обобщённого метода RBS в данном случае достаточно исследовать знак ReX на плоскости параметров (a , R) .

IS В 1.3 рассматривается задача с линейным вхождением собственной функции и нелинейным вхождением собственного числа. Предлагается численный метод, позволяющий ответить на следующие вопросы: -Где на комплексной плоскости локализован спектр задачи? -Сколько собственных чисел имеет задача заданной области плоскости? -Позволяет найти отдельные (или все) собственные числа. -Позволяет найти соответствующие собственным числам собственные функции.

Исследуется сходимость метода» Демонстрируется применение метода для задачи с дробно-линейным вхождением собственного значения вида (см. [12-17 ]):

dht d [ „. -,. АЛ - _Л _

dx² dx\ dx)

м(0) = (/(1) = Q, Я(х,Л) = .^г(*'^А)

X + <р(х)

при выполнении условий

0<К_г1 + К{х>Л)Къ* Л>0, К(х>Л)єС^}[ОДпох, г(х,Д) є С¹ поЯ, \к'_х(х,А)\< К₂.

Эта задача является модельной для исследования устойчивости электронного течения в коаксиальной магнитной линии при условии самоизоляции и для задачи с квадратичным вхождением собственного числа. Подобные уравнения возникают при решении задачи флаттера, математическая постановка которой сводится к решению матричного дифференциального уравнения {X^lE + AD + B)u^Q, X = 5 + i&

Квадратные матрицы В и D заданы и зависят от параметра v (скорости полёта самолёта). Требуется с фиксированной точностью є определить критическую скорость v*p, равную минимальному значению скорости, для которого в спектре задачи появляется хотя бы одно собственное значение Акр, не принадлежащее заданной области G.

В этом же параграфе L3 преложенный метод применяется к решению задачи с запаздыванием, в которой собственное число входит экспоненциально : D(X)u = (А + АВ + GT^r> = О Где В - С = Е (Е- единичная матрица), А -трЕхдиагональпая .

Предложенный в 1.3 численный метод исследования устойчивости

нестационарных процессов был отработан на тестовых примерах, для которых известны аналитические решения, и применён к конкретным прикладным задачам с различным типом нелинейного вхождения собственного числа,

В 1А исследуется задача на собственные значения, имеющая кусочно-непрерывный спектр и возникающая при изучении процесса каналирования электронов (см.[23]). Физический процесс описывается двумерным дифференциальным уравнением Шредингсра с условиями периодичности на концах отрезка:

pL(x>y) + vsl_y(x*y) + №-V(x,y)№x,y) = o

-CQ<у <+CQ

V{x_yy) = V{x±na_uy + ma₂)\ п,т = 0_S±1,±2,—

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

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

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

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

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

В 2.2 численно исследуется нелинейная задача, имеющая аналитическое решение. Приводятся практические рекомендации по использованию разработанной в 2.1 методики решения.

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

распространения фемтосекундного импульса в слоистой среде с учетом кубичной нелинейности описывается следующим уравнением [28]:

z{z)~+iD^^+i$t{z)U + ia{z)\U\²U = 0, 00,

U(0,t)=U(L,t)=0, *>0, U(z,0)=U_o(z), 0

Здесь введены обозначения:

**)-

ґє„ (d, +djj-\)<₁z^d_i + {d, + djj-l),

E,.

<*,+(&+<*_аХ/-і)^(<*.+^)Л

_t и e₂ -диэлектрические проницаемости слоев;

Ч. (^+^./-1)^^,+(^+^./-1),

*{*)=

[а₂₁ ^,+(^+^-1)^(^+с/,)Л

коэффициенты a_{ и a₂ характеризуют кубичную нелинейность отдельного слоя; j= 1,...,Л^ где ЛГ„, - число слоев;

Д=—

4лП р=-п1;

-СО _ 2тСС 2тгг

Q = , где w= - размерная частота света, -q =Z—— — частота периодической

структуры, p_sfr -d_x ^-\-d₂ ^/є7, следовательно, п ₌ і;

rfl =

^ь'^г ,^ - P^Jfr- -толщины слоев.

^{4^7 4V^T}

Представив решение задачи в специальном виде

U(z,t)=A{z)e^iX',

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

Dd²A

—$А-—\А\^гА = \А₉ 0

г &'

a(o)=a(l)=q.

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

21 распространяется без пространственных искажений. Безразмерные уравнения, описывающие рассматриваемый процесс, впервые введены в [29]. Уравнения имеют вид

W+i^ + iapfU+ya-bfu^O. 00 U{z,0)=U(z,L) = 0,

v(oj)=uM

Решение типа солитопа может быть записано в специальном виде

U{z,t)=A{t)e^. Соответствующая задача на собственные значения имеет вид: d'A

о]Л|³Л-/уа— А\²а)=-\А, 0

dt² ""'""' "" ''"dt [Л(0)=Л()=0,

с условием нормировки

где \л(хТ - интенсивность импульса, а и у — действительные параметры, X — собственное

значение (декремент (инкремент) затухания), і - время, L - размер исследуемой области. Для заданных значений параметров а и у необходимо найти действительные собственные значения X и отвечающие им собственные векторы. В этом случае собственный ветстор представляет собой солитон, не изменяющийся вдоль координаты распространения светового импульса. На двумерной плоскости параметров (а ,у) проводится вычисление собственных функций высокого порядка, отвечающих устойчивому режиму протекания процесса,

В главе 3 исследуется динамика развития сложных нелинейных нестационарных процессов, которые описываются системами существенно нелинейных уравнений типа уравнений Навье-Стокса (т,е. уравнениями с сильной обратной связью), по физическому смыслу не допускающими линеаризацию. Это затрудняет и фактически делает невозможным проведение анализа развития процесса на основе исследования спектра и соответствующих собственных функций, как это делалось в главах 1-2 для нестационарных процессов, описываемых операторами типа Штурма-Лиувилля и Шрёдингера и их нелинейными модификациями.

В 3.1 приводится математическая модель растекания тяжёлых жидкостей и газов по орографически неоднородной поверхности. Для процесса характерно слабое изменение

22 решения по одной из трёх пространственных компонент , что позволяет провести осреднение исходной трёхмерной системы дифференциальных уравнений Навье-Стокса по этой компоненте и затем исследовать полученную двумерную модель в среднем слое. Соответствующая система дифференциальных уравнений имеет вид :

дриб Ври²5 dpuvS _ яд{р-рь)3² _л^^&о,р ,? ,їг

cpuwo apvwo - _ і _

dpvd dpuvS dpv²S __ g d(p- p₀)S² _

Ы дх dy 2 By

dpwS dpuwS dpvwS

p₍₎)5^!L+F_Dr+K_t+F_vH,

dt dx dy

- + -

Bp5 дрйд dpvS

- +

_=a

5* dx By

. + ,

дред dpucS dpved

- + -

= ^+^+^+4-,,

dt dx dy

-+-

dph 8 dpuh б dpvh д

+

= ^*+^+^.+А

_F ^±{-_у3дй_\±i-_vSEL]
3x\ дх)ду{ By J

Ff*i = ^<рРъР\ф>^гдеФ = («,v, w,c_f A),

firV tk/ By

^=-/^-1 + -

pvS—

dy)

аналогично определяются F^^F^. Начальное распределение скоростей, размеры облака и функции источника считаются заданными.

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

23 отдельных этапах решения применяются как разностный метод решения, так и получение аналитического решения.

В 3.3 демонстрируется применение разработанного метода математического моделирования для модельных задач и для моделирования последствий промышленных аварий на примере аварии газопровода под Уфой [39], а также для моделирования ударного лунного кратера [40]. Приводятся результаты сравнения численного расчёта с известными данными МЧС.

В главе 4 исследуется возможность применения разработанного в главе 3 метода математического моделирования для моделирования МГД неустойчивости в двухслойной среде.

В 4.1 рассматривается математическая модель процесса электролиза алюминия. Система дифференциальных уравнений в частных производных описывает гидродинамические процессы и электромагнитные процессы как в среде металла, так и в среде электролита. В силу слабого изменения решения по высоте возможно провести осреднение исходной системы уравнений Навье - Стокса и тем самым получить двумерную систему гидродинамических уравнений, что существенно облегчает её численное решение. Для вычисления давления в начальный момент времени решается уравнение Лапласа в двух средах (соответственно с условиями Дирихле и Неймана). Выписывается полная осреднённая система дифференциальных уравнений: Слой жидкого металла:

ФА _| $Р\^иА , др^А -^
8t дх By

ду дх\_2

«і 2,

dpyU^ dp_lu_l Л, dp_xu_xv_xh_{ __ д
dt дх

т№-^А.С'о))*(лг-Л)-А(0

+

ду ду12

ЭауД , Эр,ы,уА j Эр, у,'ft, _ в

_{і(й,-м>о))*(л-/;г)--лС'о)}

+

-+-

= ^+^-^₀+/_? К

'іО

\.Рх)

dt

dt ' ⁺ ' Эх + (V„V)

Рх) \Р\ ) lA ^а1

Диффузионные члены определяются по формулам:

ди. V ду

+—

с- д[дй,

„ д( dvA д( dv_x

ф\

^щ дх

+ ^К*.Г+*,)

F_xjj ,F - - компоненты силы трения на границе раздела жидкого металла

и электролита:

Ъ.=-<*

^ц|^а<«. -«_a)V(Si -щ)² +0 -%У ^F& = -<* ЦЧ*. -^(51-3-^+(7,-)¹

/^ , F_;2 - компоненты силы трения на нижней границе жидкого металла:

^v₀ =-^л^?1л/"1²⁺^²

і⁷,, - теплообмен между металлом и электролитом на границе раздела:

Так как удельная электрическая проводимость электролита много меньше проводимости алюминия, то Т₂ > Г,.

j%_o - теплообмен между металлом и дпом ванны, температуру которой будем считать

постоянной:

F_Za=a,(f_]-f_a)

J\ — ток в слое металла. Слой электролита;

.+-

+

dp₂h_{2 i} dp₇u₂h₂ ^dp₇v₂h₂ _

dp₂u₂h_{2 |} dp₂u₂²h_{2 |} dp₂u₂v₂h₂ _ S

+

+ ^ +^, -^, -gp₂h₂^ + f_2ih₂

^v_a^i , др_гіі_гУ_гЬ_{г |} dp₂v₂h₂ _ 3

By ду

^Vh-h₂(.i_a))*(p₂g-f_2l)-p₂(t₀)

др_гс Т₂Ь₂ др_гй₂с f₂h₂ dp₂v₂cf₂h₂

. + _^ + 1 _{= F}^ +^ __F//| +/^+2^

+

= 1—-V V₂ roti-Л.)

І є_гр_г а₂

Диффузионные члены определяются формулами:

F_n!r =

ди₂} д Г ди₂ дх J ду{ ду

э=т л

_{F =}_а_[ э^_

+

dv.

а*

ад:

^««Ы(*_аг+*_а>-^

+

а^

ду j

F _п , F тт -компоненты силы трения на верхней границе электролита, при этом под

анодами:

F_lMl =-С,_1гР₂^иг4^и2^г+у.

а в остальной области эти силы равны нулю.

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

Qx* - джоулево тепло, h — ток в электролите. Заданы следующие граничные условия:

Для скоростей в слое жидкого металла:

-V,

= о

Для магнитного поля в металле;

-H.I =о

Для плотности электромагнитных сил в металле:

= 0 г

Для скоростей в слое электролита: = 0

&

Для магнитного поля в электролите:

-й₂| =0 dt ^г|г

Для плотности электромагнитных сил в электролите:

^?*1г=

Силы на границе в электролите заданы равными нулю, потому что плотность тока на границе в электролите есть 0, так как гарнисаж не проводящий.

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

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

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

27 Автор выражает искреннюю признательность заведующему кафедрой вычислительных методов факультета ВМиК МГУ им. М В. Ломоносова академику РАН, профессору Александру Андреевичу Самарскому за внимание и интерес к моей деятельности в области вычислителыюй математики; автор горячо благодарит академика РАЕН, профессора физического факультета МГУ Рунара Николаевича Кузьмина за помощь в постановке задач и обсуждения результатов; автор приносит глубокую благодарность профессору факультета ВМиК МГУ Алексею Владимировичу Гулину за постоянное внимание к работе и ценные замечания; автор благодарит заместителя директора ИММ РАН профессора Владимира Фёдоровича Тишкина и к.ф.-м,н, , с.н.с. ИММ РАН Андрея Александровича Кулешова за полезные обсуждения результатов работы; автор благодарит заведующего лабораторией вычислительных методов в физике факультета ВМиК МГУ профессора Вячеслава Анатольевича Трофимова за совместную работу и постановку задач.

Некоторые асимптотические свойства системы дифференциальных уравнении с .постоянными коэффициентами

Очевидно, что доказательство соотношений (1.1,5) и (1.1.6) не изменятся в случае, когда для кратных корней Х\ и Х имеется столько линейно-независимых собственных векторов, какова их кратность. Не ограничивая общности, предположим теперь для простоты, что собственному числу Х\ кратности 2 соответствует один собственный вектор Vi и один присоединенный вектор V2 (все остальные собственные числа считаются однократными). Тогда, разложив УГ по каноническому базису, получим в этом случае = (1+ +0 1 + + x -xf

Рассмотрев далее отношение при s- co, нетрудно убедиться в справедливости (1,1.5), Однако в этом случае скорость сходимости замедляется. Аналогично доказывается (используя правило Лопиталя) справедливость (1 Л.5) и (1.1.6) в случае более сложной структуры жордановой формы матрицы А. В случае наличия в спектре матрицы А комплексно-сопряжённых собственных чисел Х\ и Xi доказательство Теоремы 3 проводится для матрицы A-icE, где с -константа, і - мнимая единица.

Выпишем алгоритм нахождения Х\ и Х , в основе которого находится явный метод Эйлера (1.1.4) (метод RBS). Будем считать, что шаг интегрирования удовлетворяет условиям Теоремы 3. Сходимость метода RBS вытекает из сформулированных выше теорем. Алгоритм метода RBS: 1. задаём ;Л И = Г, где .("-I) _ U + I //-(Hi) =?( +1К&5 2. выбираем т 0 (т 0) и вычисляем х{ } = x[stiW(x 0+0 х{ } вычисляется по формуле (1.1.4) 3. если Ідг(н +1 — JC10 , ,то XN = (Sj +Ij »xs,)/n:i0s иначе идём к пункту 2.

Здесь ег 0 - заданный параметр, х1 }&0. При т 0 по пунктам 1,2 3 вычисляется а.1.Для вычисления границ мнимой части спектра в алгоритме RBS достаточно заменить матрицу А в (1.1.4) на ІА.

Рассмотрим теперь матрицу А, у которой на прямой х = ReA/ находятся два собственных значения X } и X 2- Теорема 3 говорит о том, что пределов (1.1.5) и (1.1.6) в этом случае не существует. Однако найти границы спектра матрицы А можно найти при помощи метода RBS, если предварительно осуществить поворот комплексной плоскости на некоторый угол ф (острый). При этом возникает задача Ах = fix (1 Л.9) где A = Aei9 9 fi = № При этом на прямой х = Re/if, где /Jj - собственное значение с минимальной действительной частью задачи (1Л.9)? уже находится одно собственное значение, и по Теореме 3 для его нахождения можно применить метод RBS.

В некоторых случаях выбранное согласно Теореме 3 значение т является слишком мелким. Выход из этой ситуации заключается в применении неявных разностных методов, т.е, в замене (1Л.4) на неявную схему Эйлера. Подробно такая модификация метода RBS рассмотрена в работе [2].

Алгоритм RBS применялся к расчету границ спектра различных матриц, спектр которых был известен заранее. Результаты расчётов позволили сделать вывод о том, что алгоритм RBS может быть использован для быстрого нахождения граничных точек комплексного спектра и соответствующих им собственных векторов в случае отсутствия начальных приближений к исходным собственным числам (см. [1]). Особенно эффективен метод RBS в том случае, когда исследуется знак вещественной или мнимой части собственного числа, отвечающий за устойчивость исследуемого процесса.

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

Итак, рассмотрим следующую алгебраическую задачу, элементы матрицы которой зависят от двух действительных параметров є, а : А(ъйа)у= %у, (1.1.10) где А - квадратная матрица порядка (ЛМ) с комплексными элементами, А,-собственное значение. Требуется па плоскости параметров найти область, в которой все собственные значения задачи (1.1.10) удовлетворяют заданному свойству (действительные, отрицательные). Ниже предлагается численный алгоритм для построения нейтральных кривых устойчивости физических процессов, зависящих от двух действительных параметров.

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

Предположим, для устойчивости процесса достаточно, чтобы собственное значение задачи (1.1Л0), имеющее максимальную действительную часть, было отрицательным. В [3] показано, что в случае действительного спектра погрешность между собственными значениями дифференциальной и л точное "і \kif « разностной задач Кк \ \ минимальна для малых номеров к и сильно возрастает с ростом к. Таким образом, для действительного спектра численное исследование знака минимального собственного значения даст гораздо более точный результат, чем численное исследование знака максимального значения.

Умножив правую н левую части уравнения (1.1.10) на (-1) и обозначив В —— А, \ — — Х, приходим к задаче, устойчивость решения которой зависит от знака собственного значения с минимальной действительной частью

Алгоритм нахождения собственных значений и соответствующих им функций нелинейной задачи на собственные значения

Численные исследования показали, что итерационные процессы (2,1.4) сходятся к к-ому собственному значению и соответствующему собственному вектору задачи (2Л .1) при выполнении определенных условий:

I. Начальное приближение для итерационных процессов (2.1.4) должно иметь структуру, соответствующую точному решению, то есть количество нулей в точном решении J ив начальном приближении у должно совпадать. При этом не важна форма начального приближения. При начальном приближении различной формы (треугольники, ступеньки или синус) метод сходится к к-ому собственному вектору н соответствующему собственному значению. Различные начальные приближения для первого собственного вектора показаны на рисунках 2ЛЛ -2,1.3, для второго - на рисунках 2.1,4 -2Л.6, для третьего - на рисунках 2Л.7-2Л.9. Начальные приближения для остальных собственных векторов берутся аналогично. Этот вывод хорошо согласуется с теоремами А.А,Самарского [26] для линейной задачи на собственные значения: ТеоремаЬ При кеази-закреплении основной области Т по множеству емкости нуль собственные значения уравнения Аи + Aw = 0 при граничном условии на Г и=0 остаются неизменными.

Теорема2,Для того, чтобы наименьшее собственное значение уравнения Дн +Aw = 0 при граничном условии на Ги=0 оставалось неизменным при квазизакреплении основной облети Т по множеству Е, необходимо и достаточно, чтобы ёмкость множества Е была равна нулю.

Параметр т должен быть меньше некоторого т0? определяемого экспериментально. Значение Т0 зависит от вида нелинейности элементов матрицы А и от начального приближения у . Чем больше особенностей имеют функции, входящие в элементы матрицы А", тем меньше должен быть параметр

V Численно показано, что для сходимости метода (2,1,4) требуется очень маленький итерационный параметр т, следовательно, метод сходится медленно. Если для данного т процесс расходится, надо уменьшить шаг. Также могут возникнуть случаи, когда при некоторых х} и т2 метод сходится, но к различным решениям.

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

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

IV. Численные исследования показали, что сходимость метода немонотонная, т.е. величина max 1=1 лм у" —у"\ может увеличиваться на некотором этапе вычислений, но во всех сделанных расчетах при приближении к точному решению метод всегда начинает сходиться монотонно. V. Численные расчеты также показали, что абсолютная погрешность вычисления собственного значения Хк минимальна для номеров Л, соответствующих собственным числам с минимальной по модулю действительной частью, и сильно возрастает с ростом модуля действительной части. Относительная погрешность также возрастает, по не так быстро. Таким образом, для собственных чисел с большой по модулю действительной частью молото говорить только о порядке и знаке собственного значения Хк, не оценивая его величину. При этом погрешность вычисления собственного вектора значительно меньше погрешности вычисления соответствующего ему собственного числа.

Рассмотрим применение численного метода нахождения комплексного спектра к нелинейной задаче на собственные значения, имеющей аналитическое решение. Рассмотрим задачу d2V = Л V, где 0 x L, dr (2.2.1) іф)=г(фо Собственные значения и собственные функции задачи (2.2.1) определяются по формулам: А-=—I Vk(x) sin-—t А: = 1Д.... Для численного решения задачи (2.2.1) введем равномерную сетку с шагом А (й=\х = jht j = 0.„N,h=—к Определим сеточную функцию v= voiv, v2 " v.v-i v,v) Рассмотрим разностную задачу, аппроксимирующую дифференциальную задачу (2.2.1) со вторым порядком точности v = vi i=lX...,N-l, vo=vw=0, h-N = L (2.2.2)

Собственные значения и собственные векторы задачи (2.2.2) известны и находятся по формулам ([3] стр.40): 4 Л. = —5Ш h2 2L ук(х,) т , k=l,2,...,N-l, i=0,\,...,N . Также в [3] показано, что собственные векторы задачи (2.2.2) образуют ортогональный базис. Пусть V=lf , тогда задача (2.2.1) переписывается в виде dU - = Яи\ 0 x L, (2.2.3) dx2 U2{O)=U2(L)=O , x = i-h, і = 0,1,...Л k=\,2f...,N-l, h-N = L. С помощью тождественных преобразований, обозначив // = —, и так как U тождественно не равно 0, получим новую задачу с нелинейным вхождением собственных функций в оператор d2U \(dU\ Т! /(0)=/(L 0, Задача (2.2.4) имеет решение (2.2.4) А =2" ) к=и іп ,АИД... .

При этом pk(x)f =\Vk(x)\ , =1,2,,... 2,2.2,Построение разностной схемы. На равномерной сетке о определим разностную функцию (Уо р - Л-м ) и построим для дифференциальной задачи (2.2.4) разностную задачу со вторым порядком аппроксимации. Оператор A{U) представляется в виде сумми двух. операторов А(и)=А{(и)+Аг{и), d2U где At{u)= - линейная частії . ,__х і (duY h, Л2\у } = —т I - нелинейная часть dx2 " /v U\dx оператора А. Линейный оператор Аг аппроксимируется со вторым порядком аппроксимации по известным формулам А}у=г — у х. В диссертации исследованы различные способы аппроксимации нелинейного оператора A2\U)= , В первом случае матрица А2 U V dxJ несимметричной, во втором случае - самосопряженной.

Получение двумерной модели,

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

Основная концепция самой математической модели была заложена в [30], где исходной является трехмерная система для осредненных по времени параметров течения :

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

Поэтому мы можем проинтегрировать исходные трёхмерные уравнения газовой динамики по высоте от ZQ\x y) до rL[x9y9t) и получить двумерную систему уравнений для определения изменения средних по высоте характеристик потока в горизонтальном направлении. Для адекватного учета явления разбавления облака газа окружающим воздухом будем считать, что плотность р концентрация с, вертикальная скорость w и энтальпия h имеют неоднородное распределение по высоте облака р = p(x,y,ztt)r c=c(x,yrz,t)r w = w(x,yfz,t)9 h=h(xry,z,t); z0(xry) z H(x,yrt)r где ZQ- рельеф подстилающей поверхности, H -верхняя граница облака.

Будем предполагать, что все частицы облака газа или жидкости по высоте при любых фиксированных (xtyj) движутся с одинаковыми горизонтальными скоростями, за исключением лишь окрестности нижней границы облака, то есть it(xty,zj) ъй (x,yj), v(x,y,z4t) nv(x,ytt).

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

В окрестности нижней границы течения у поверхности земли скорости резко замедляются и на нижней границе выполнены условия u=v J=I -0, uA = uy = uz — vx = vy = v = wx = Wy _ 0 при любых (x,yrt). Нас будут интересовать осредиекные по высоте облака параметры течения p{x9y,t) = — \p{x4y t)dz, (3.1.2) где b(x,y\t) = H(x,y,t) - zc(x,y) - толщина облака газа или жидкости. Аналогично определяются с (xy9t)t h(xfytt).

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

Проинтегрируем уравнения системы (3.1.2) по высоте от нижней границы z0(x,y) до верхней границы H{x,y,t) и получим двумерную систему уравнений для определения изменения в горизонтальном направлении средних по высоте характеристик течения. При этом толщину 5(х, у, t} будем также рассматривать как одну из искомых функций двумерной модели. Рассмотрим третье уравнение исходной системы: дон д/л/w Bpvw dpw2 dp -—+--— + - -—+ ——+— = Fw-pg, сі дх ду dz dz где ,, д , dw. д , 9и\ д , dwx дх дх ду ду dz dz

Так как. газ (жидкость) - молекулярно слабо связанная система, можем расщепить уравнение для w на два уравнения, одно из которых описывает процесс переноса вертикальной скорости,, а второе - изменение давления за счет гравитации: dpw , dptiw , dpw , dpw2 _ +— + . . +—_ (3.1.3) o( ex су oz %—( (3.1.4) dz

Рассмотрим уравнение (3.1.3). Пусть po , Po - известные фоновые значения давления и плотности на верхней границе, которые предполагаются постоянными. Тогда р = р- po , р—р - ро - избыточная плотность и избыточное давление в облаке.

Получение уравнения электромагнитной индукции, учитывающего реальную подводку тока к электролизеру

Рассматриваемая модель описывает в динамике все термодинамические параметры течения, электромагнитное поле в электролите и жидком металле, колебание поверхности раздела: металл-электролит. Однако эта модель не учитывает изменение магнитного поля, обусловленное токами проводимости, подведёнными к электролизной ванне, поля возмущающих электромагнитных сил, давление на настыль (гарнисаж) в средах металла и электролита, а также конкретную конфигурацию настыли и анодов. Т.о. в модель (4.1.1)-(4Л.З) необходимо внести соответствующие изменения, которые неизбежно повысят адекватность самой модели.

В первую очередь необходимо ввести удобным образом декартову систему координат. Пусть ось х будет направлена вдоль длины ванны, ось у от глухой стороны к лицевой, а ось z направлена вертикально вверх. Начало координат совмещено с углом ванны? к которому ближе первый анод, на уровне подошвы слоя металла. Введенная система координат изображена на рис. 4Л. 1. 224 Будем использовать следующие обозначения высоты слоев: zQ(xty) - высота настыли, отсчитываемая от нулевого уровня (уровня подины); Aj( ,y,f) -толщина слоя металла; Ht(x,ytt) = zQ(x,y) + h{(x,y,t) -высота границы раздела; h2(xrytt) -толщина слоя электролита; H2(x,y,t) = Н\(x,y9t) + h2(x,y,t) - высота поверхности электролита от нулевого уровня.

Величины, относящиеся к металлу, будем снабжать индексом "Iй, к электролиту -индексом "2". Значение величины в среднем слое при необходимости будем обозначать чертой сверху, например: р. Впрочем, во многих случаях значение величины в среднем слое будет подразумеваться (за очевидностью) и без черты. Мараліетішзісидкого_штадп а: VY =(«j,vMW!) - вектор скоростей, рх - плотность жидкого металла, рг - давление в жидком металле, Т} -температура жидкого металла, Н,= (HixfHi уУНи) -вектор напряженности магнитного поля в слое жидкого металла, f\ =(fwf}y fiz) —плотность электромагнитных сил в металле, fT - нагревание металла за счет действия электромагнитных сил, Е( - вектор напряженности электрического поля в жидком металле, , = 1 - диэлектрическая постоянная жидкого металла, vim - магнитная вязкость среды, 7 - электрическая проводимость жидкого металла. Параметры электролита: V2 = (u2,v2,w2) -векторскоростей, р2 -плотность электролита, р2 - давление в электролите, Т2 -температураэлектролита, H2 = {II2х, II2у,Н1г) - вектор напряженности магнитного поля в слое электролита, 1г "(Лл Л Лг) -плотность электромагнитных сил в электролите, fT - нагревание электролита за счет действия электромагнитных сил, у2гя - магнитная вязкость среды,

Е: - вектор напряженности электрического поля в электролите, сг - диэлектрическая постоянная электролита, т2 - электрическая проводимость электролита. J„p - ток проводимости. Q & 2 -10"7 м I сек — скорость осаждения металла из электролита. о(х У) " функция помехи (образованной анодами, погруженными в электролит) в уравнении движения. Рассмотрим среду жидкого металла. Уравнение для давления в металле [42] имеет вид: "Л + /.. (4-1.3) проинтегрировав по z от zo до z» получаем: А(2)"А(2о)3 Лг(г-го) + /и(2-2о) или, что то же самое: Л(г) = Яі( о)"Аг(г- ) + /іД2-го) (4Л-4) Далее получим осредненное значение, проинтегрировав (4.1.4) по z от z0 до Нь учитывая что А, =НХ -z0: +/?,(X, , Z," )

В данном выражении неизвестно значение давления pi(x,y,zo,t). Поскольку это давление значительно не изменяется с течением времени, делаем приближение, что Pi(x,y,zo,t) равно своему значению в начальный момент времени. Теперь встала задача о поиске давления на настьшь в начальный момент времени.

Похожие диссертации на Численное исследование устойчивости некоторых физико-технических систем

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

Итерационные алгоритмы анализа стохастической устойчивости колебательных системГубкин Андрей Анатольевич

Математическое моделирование и анализ статической устойчивости автономных электроэнергетических системМадонов Анатолий Николаевич

Методы качественного анализа устойчивости математических моделей динамических системКарпухин Владимир Борисович

Устойчивость математических моделей систем фазовой синхронизацииПеркин, Алексей Александрович

Устойчивость и предельное поведение открытых репликаторных системПавлович, Екатерина Николаевна

Резонансные свойства детерминированных математических моделей и устойчивость функционирования технических систем с гистерезисными нелинейностямиТолоконников Павел Вячеславович

Оценки устойчивости для нестационарных марковских моделей в системах массового обслуживанияКоротышева Анна Владимировна

Математическое моделирование напряженно-деформированного состояния и устойчивости пластин и пологих оболочек с построением систем аппроксимирующих функцийАбросимов Алексей Анатольевич

Компьютерные модели и программные средства для исследования биомедицинских системНасер Нихад Махмуд