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

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

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

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

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

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

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

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

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

— задачи с кусочно-непрерывным спектром для оператора Шредингера.

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. Область локализации на комплексной плоскости всех или части собственных значений (например, собственных значений с ЯеЛ>0).

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

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

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

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

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

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

Т.о. в литературе достаточно широко представлены численные алгоритмы для решения линейных по собственному вектору задач с нелинейностью по собственному значению [79−85], и при этом практически не рассматриваются задачи, нелинейные по вектору. Единого механизма решения таких задач не существует. Большинство работ в этой области посвящено аналитическому исследованию свойств определенных достаточно простых типов задач в некоторых пространствах. Например, в работах [86−87] исследуется задача ф (и)=ХЧ" (и), мбФ" '(а), где ф (и)=— J|Vm|2 dx — — Jm2 dx, vP (m)= определены в пространстве.

2 n 2 n n q (Q). Автор доказывает, что в зависимости от того, является, а положительным или отрицательным, задача имеет одно или два решения.

В конце 70-х годов группа немецких математиков (Kurt, Lackner, Amann и другие) предложили численный метод решения нелинейных задач, основанный на методе Пикара [78,79,88−95]. Рассматривались задачи типа f (x, U (x)) = XL (u), хеП, U (x) = 0, хедП, n где L (u)=- dt{ал (х)дкU (х)) + д (х), QcRN, «к — действительное положительное собственное значение, при условии, что функции/ aik и граница дО. достаточно гладкие, а матрица, А = (aik) симметричная и положительно определенная. Такие задачи возникают в физике плазмы. Авторы предлагают следующий итерационный метод решения дифференциальных задач на собственные значения.

Шаг 1 (Итерационный метод Пикара) < L (un+l)=f (x, un (x)), xeQ, С7л+1(д:)=0, хеП.

Шаг 2 (Нормализация) ^ и -1 ^ |гу л+1 £угл+1.

Таким образом, авторы рассматривали только те задачи, у которых нелинейная часть уравнения является функцией от U, производные U' не входят в нелинейную правую часть уравнения. И при этом в работе не предложен алгоритм для нахождения собственного значения, находится только собственный вектор и не определяется, какой моде соответствует этот вектор. Численные результаты для такого типа задач получаются с достаточной точностью, но четкого доказательства сходимости нет. В работе [114] рассматриваются дифференциальные задачи типа.

— Аи=X/'(«), ueQ, и = 0, uedCl, где и = и (х) — действительная функция, определенная в области QcRw, Xдействительное собственное значение. Авторы получают решение ([и, Х), численно исследуя следующую вариационную задачу минимизации.

Е (и):=— [[Vm|2 dx -" min на Jf (u)dx = y 2 n n для различных ограниченных значений у. При этом авторы получают только собственный вектор, соответствующий минимальному собственному значению. С помощью рассмотренного метода нельзя найти остальные собственные значения и векторы.

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

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

Для разрешения многих вопросов относительно механизма процессов переноса вещества и энергии в тяжелых газах были поставлены специальные эксперименты. Наиболее полная серия опытов по рассеянию тяжелого газа была выполнена в Англии в 1982;1984 гг. на полигоне острова Торни-Айленд [30]. Исследовалось рассеяние больших объемов предварительно подогретой до температуры окружающего воздуха смеси фреона-12 и азота. Контейнер с газом объемом 2000 м³ представлял собой разборный пластиковый цилиндр высотой 13 м. Его конструкция предусматривала возможность быстрого раскрытия стенок и практически мгновенного истечения смеси. Эксперименты финансировались группой промышленных компаний топливно-энергетического направления различных стран мира.

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

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

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

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

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

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

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

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

1.Получение основных полей в стационарном состоянии.

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

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

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

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

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

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

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

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

О.

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

В § 1.2 исследуется обобщённая задача на собственные значения с линейным вхождением собственных функций <�р (у) и собственного числа X. Задача возникает при исследовании течений Пуазейля и Блазиуса в рамках краевой задачи на собственные значения для уравнения Орра-Зоммерфельда и в исходной постановке описывается дифференциальным уравнением четвёртого порядка (см. [ 8 ]) р{4) =d (y, A).

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

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

РУ-± 1 =H,-±i =0.

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

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

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

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

0</С, < + К (х, Л)<�К2, Л> О, К (х, Л) е С' [0,1] по х, г (х, Х) е. С1 по Л, < Кг.

Эта задача является модельной для исследования устойчивости электронного течения в коаксиальной магнитной линии при условии самоизоляции и для задачи с квадратичным вхождением собственного числа. Подобные уравнения возникают при решении задачи флаттера, математическая постановка которой сводится к решению матричного дифференциального уравнения {ЛгЕ + ЛО + В) й = 0, Х = 5 + ico.

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

В этом же параграфе 1.3 преложенный метод применяется к решению задачи с запаздыванием, в которой собственное число входит экспоненциально: D (A)u =(А + ЛВ + Се'*)и = 0.

Где В = С = Е (Еединичная матрица), Атрёхдиагональная .

Предложенный в § 1.3 численный метод исследования устойчивости нестационарных процессов был отработан на тестовых примерах, для которых известны аналитические решения, и применён к конкретным прикладным задачам с различным типом нелинейного вхождения собственного числа. м (0) = м (1) = 0, К (х, Л) = -{Щ^.

Л + <�р (х при выполнении условий.

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

Г" «(Х'У) + У1у (х>У) + (л~ V{x, y))if/{x, y) = О.

— 00 < +00, —со<�у< +00.

И*оО|±в1 <�">

V (x, y) = V (x + naXiy + ma2) п, т = 0,±1,+2,.

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

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

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

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

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

В § 2.3 демонстрируется применение предлагаемой методики к задаче исследования распределения фемтосекундных импульсов в фотонном кристалле. Процесс распространения фемтосекундного импульса в слоистой среде с учетом кубичной нелинейности описывается следующим уравнением [28]: т л2г г z{z)-— + i + + fa (z)|f/|2f/ = 0, 0 0,.

U (0,t)=U{L, t)=0, t>0, U (z, 0) = U0(z), 0.

Здесь введены обозначения: ч fe,. (dl+d2Xj-).

4Z) |e2, dx +djj-l).

4пП p=-7iQСО 2тtC 27IC fi=——, где со =- — размерная частота света, ^ 'а — частота периодическои Р структуры, ^ + d2, следовательно, Q =.

Я, = Jl? i.

Представив решение задачи в специальном виде U (z, t)=A{z)eiX', получим задачу на собственные значения с нелинейным вхождением собственной функции в оператор

D д2Л.

— Р А—А2 А = ХА, 0.

8 dz a (o)=a (l)=o.

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

1аиги + уаЦиги)=0, 0О dz dt1 11 dt4 1 ' U (z, 0)=U (z, L)=0, u (o, t)=uM.

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

U{z, i)=A{t)eiX Соответствующая задача на собственные значения имеет вид: aA2A'iya-^fA2А)=-А, О.

И<11,=М<�Р<=>> о где |л (-с)|2 — интенсивность импульса, а и у — действительные параметры, «к — собственное значение (декремент (инкремент) затухания), t — время, L — размер исследуемой области. Для заданных значений параметров, а и у необходимо найти действительные собственные значения X и отвечающие им собственные векторы. В этом случае собственный вектор представляет собой солитон, не изменяющийся вдоль координаты распространения светового импульса. На двумерной плоскости параметров (а, у) проводится вычисление собственных функций высокого порядка, отвечающих устойчивому режиму протекания процесса..

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

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

• + dt дх ду дрд дрйЗ dpvS FDw + FwH +w|ieSoQ, Q, dt дх ду д~рсд дрис5 dpvcd dt дх ду FDc + FcH +FC2q +czeSo.

Q, dph5 dpuhd dpvhS | at ax ду о ^ где д (^дйЛ д ^дй.

Du дх J ду ^ ду.

F9H =-^p0U.

F =—{pv5— ckv дх) ду «dw pv5 ду аналогично определяются FDc, FDh. Начальное распределение скоростей, размеры облака и функции источника считаются заданными..

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

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

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

В § 4.1 рассматривается математическая модель процесса электролиза алюминия. Система дифференциальных уравнений в частных производных описывает гидродинамические процессы и электромагнитные процессы как в среде металла, так и в среде электролита. В силу слабого изменения решения по высоте возможно провести осреднение исходной системы уравнений Навье — Стокса и тем самым получить двумерную систему гидродинамических уравнений, что существенно облегчает её численное решение. Для вычисления давления в начальный момент времени решается уравнение Лапласа в двух средах (соответственно с условиями Дирихле и Неймана). Выписывается полная осреднённая система дифференциальных уравнений: Слой эюидкого металла: дрА dpxuxhx | dpxvxhx dt дх ду dpxuxhx^dpxu2hx ^ dpiu]v^hl д dt дх ду дх дрх у, hx др, ц, уД dpxvxhx д.

Н—:—г dt дх.

W)*(Pig-A,)-Pi (t о) ду ду |2 fd*1 +Гум1 +FyZo +fXyhx dPlCvJA, dPl" lCv, ТА dPiVxcVi Txh dt дх by d dt.

11,1 ГнЛ fH, J) +(V, V) 1 = —, V.

Px) I л J IA) FDTx+FHi-FZ9+f hx I.

J, ca.

V,—rof (-4 eP.

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

Мт.

F =JL ?El)+JL.

D" 1 дху т дх) ду.

F =А (lLA йх^ т дх J ду.

Мт.

DT> дх дТ, дх ay! д] ду хнх' ^уП," компоненты силы трения на границе раздела жидкого металла хн, ~ G я, и электролита:.

Р +Р2 /77 м, -M2)V (Щи2У +(v, -v2)2.

20 > Fyz" ' компоненты силы трения на нижней границе жидкого металла:.

Руг,.

Fllt — теплообмен между металлом и электролитом на границе раздела: рн, =<�хх (Т2 -Тх).

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

F2q — теплообмен между металлом и дном ванны, температуру которой будем считать постоянной:.

Fz0 -«offi-Ji).

J — ток в слое металла..

Слой электролита: dp2h2 [ dp2u2h2 | dp2v2h2 = dt дх ду dp2u2h2 | dp2u22h2 | dp2u2v2h2 8 1 a, & - * dp2u2y2h2 | dp2v2h2 д dt dx dy dy 12 1 h2-h2(t0))*(p2g-f22)-p2(t0) Foh + ^ ~ Fy, h — gp2h2 ^ + f2yh2 dp2cv T2h2 dp2u2cv T2h2 dp2v2cv T2h2 dt.

8 dt рг J (V2,V) dx X ч Рг.

Рг.

S2p2.

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

F = д (дм, дл: ду ди2 Л ду г S Г ^ V дх J ду.

3 f dv, Л гМт.

М д ду к2Т + к2) — 2 ду.

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

У" 2 анодами.

Fxh2 =<�н1Ргиг^иг рун2 = ~Си2 РгЪ, а в остальной области эти силы равны нулю..

FH — характеризует теплообмен с воздухом в каналах: Рщ =сс2{Т2-Тв) ..

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

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

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

-н,| =0 dt 1|г.

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

5* 0 г.

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

Sn ! Г.

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

-н2| =0 dt 2lr.

Для плотности электромагнитных сил в электролите: f2| =0 zlr.

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

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

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

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

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

Основные выводы диссертации сводятся к следующему:.

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

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

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

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

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

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

Заключение

..