Содержание к диссертации
Введение
Глава 1 Моделирование в газовой динамике 14-55
1.1. Значение математических моделей в технике и физике 15-23
1.1.1. Постановка задачи моделирования 15-16
1.1.2. Схема процесса моделирования 16-20
1.1.3. Классификация математических моделей 20-21
1.1.4. Основы математического моделирования 21-23
1.2. Математические модели в газовой динамике 23-54
1.2.1.1. Двумерные стационарные течения: плоские и осесимметричные 23-26
1.2.1.2. Теорема Крокко о вихрях 26-28
1.2.1.3. Потенциал скорости; уравнение для потенциала скорости в сжимаемом газе 29-31
1.2.1.4. Переменные годографа. Уравнение Чаплыгина 31-32
1.2.1.5. Особенности решения задач в переменных годографа 32-34
1.2.1.6. Уравнение Эйлера-Трикоми - трансзвуковой аналог уравнения Чаплыгина 34-37
1.2.1.7. Уравнение для потенциала плоского почти однородного трансзвукового потока газа 37-38
1.2.2. Невязкие сжимаемые течения 38-45
1.2.3. Трансзуковые невязкие течения 45-54
1.2.3.1. Общие замечания 45-48
1.2.3.2. Трансзвуковое уравнение малых возмущений 48-51
1.2.3.3. Полное уравнение для потенциала 51-54
1.3. Метод фиктивного газа 54-55
1.4. Выводы 55
Глава 2 Численные методы решения уравнений газовой динамики 56-79
2.1. Численные методы решения уравнений математической физики 57
2.1.1. Теория аппроксимации, устойчивости и сходимости разностных схем 58-62
2.1.2. Методы численного решения задач математической физики 62-69
2.1.3. Вычислительные методы в линейной алгебре 69-76
2.1.4. Вопросы оптимизации численных методов 76-79
2.2. Метод простых итераций 79-95
2.2.1. Основные разновидности итерационных процессов 79-82
2.2.2, Метод простой итерации 82-95
2.3. Метод сопряженных градиентов 95-109 2.3.1. Вариант метода сопряженных градиентов 104-109
2.4. Метод характеристик 110-129
2.4.1. Характеристики дифференциального уравнения
в частных производных второго порядка 110-113
2.4.2. Приближенное построение сетки характеристик 119-123
2.4.3. Нахождение характеристических полосок 124-129
2.5. Выводы 129
Глава 3 Алгоритмы моделей, численные схемы и результаты 130-152
3.1. Алгоритмическая схема поставленной задачи
нахождения безударного аэродинамического профиля 131-133
3.2. Численные схемы реализации алгоритма модели нахождения безударного профиля 134-139
3.3. Результаты численного эксперимента 139-150
3.4. Выводы 150
Заключение 151-153
Список литературы
- Схема процесса моделирования
- Особенности решения задач в переменных годографа
- Методы численного решения задач математической физики
- Численные схемы реализации алгоритма модели нахождения безударного профиля
Введение к работе
При обтекании аэродинамического профиля потоком воздуха с большой до звуковой скоростью, характерной для пассажирской и транспортной авиации, около профиля образуется зона сверхзвуковой скорости течения, которая обычно заканчивается ударной волной. Наличие ударной волны может вызвать более ранний переход ламинарного течения в турбулентное и связанное с этим резкое снижение аэродинамического качества крыла. Кроме того, возникновение ударных волн средней и большой интенсивности приводит к существенным волновым потерям и, что более важно, может инициировать отрыв пограничного слоя и возникновение нестационарных явлений типа бафтинга.
Поэтому одной из важнейших задач аэродинамики летательных аппаратов является проектирование профилей и крыльев безударной формы, для которых ускорение и торможение потока в местной сверхзвуковой зоне происходит без образования ударных волн, и реализуется гладкое обтекание.
Методы проектирования делятся на две категории. К первой категории относятся комбинированные численные методы, объединяющие решение прямой задачи аэродинамики и решение задачи оптимизации, позволяющей в ходе итерационного расчета изменять форму профиля, чтобы минимизировать целевую функцию, например, сопротивление или разность между полученным на данной итерации и желаемым распределениями давления [67], [81]. Недостатком таких методов являются большие затраты времени расчетов, вследствие чего их применение в практике аэродинамического проектирования связано со значительными неудобствами и затратами средств. Ко второй категории относятся обратные методы, предполагающие задание распределения давления или скорости на профиле и определение формы профиля в результате решения соответствующей краевой задачи. В исследованиях, связанных с решением обратных задач для трансзвуковых течений, сформировалось несколько подходов. В методе Вольпа и Мельника [10] проводится решение полного уравнения для потенциала скорости с пространственными координатами в качестве независимых переменных; при этом в заданное распределение давления вводится свободный параметр, который в процессе расчета подбирается так, чтобы обеспечить существование решения обратной задачи. В работах [63], [96] применено преобразование годографа, позволяющее получить линейные уравнения для описания изоэнтропических течений. В методе фиктивного газа [97] используется степенное выражение для плотности газа, обеспечивающее эллиптический тип уравнения для потенциала скорости во всей расчетной области. Были предложены также численные методы, позволяющие задавать как форму профиля, так и распределение давления на профиле или исследовать смешанную задачу, в которой часть профиля задана, а на другой его части задается распределение давления или скорости [77].
Преимущество обратных методов аэродинамического
проектирования перед указанными выше методами оптимизации -более высокая скорость счета и экономичность, позволяющие осуществлять проектирование профилей в интерактивном режиме работы на ПК или рабочей станции. Недостатком указанных обратных методов является возникновение некоторых принципиальных трудностей, связанных с необходимостью выполнения условия замкнутости контура и условия совпадения скорости набегающего потока с заданным значением. Как показывают предварительные исследования, задание распределения давления может приводить к наличию вогнутых участков спроектированного аэродинамического профиля и к возникновению внутренних ударных волн в местной сверхзвуковой зоне, поскольку оно представляет собой жесткое граничное условие, не адаптированное к структуре течения. Результаты некоторых расчетов [98], действительно, демонстрируют наличие внутренних ударных волн.
Схема процесса моделирования
Решение прикладной проблемы с применением математики неизбежно требует математического моделирования. Общая схема организации такого моделирования показана на рис. 1.1. Выделенные штриховым контуром элементы являются компонентами программного обеспечения ЭВМ или этапами машинной обработки.
Адекватность модели объекту всегда ограничена и зависит от цели моделирования. Рекомендуется выбирать модель минимальной сложности при заданной точности либо максимальной точности - при заданной сложности. Понятие сложности обычно оценивается длиной программы. Принцип баланса точности требует соизмеримости погрешностей, вызываемых различными причинами: неполнотой модели, неточным заданием ее параметров, погрешностью выбранного численного метода, счетом в ограниченной разрядной сетке.
При исследовании сложных систем может потребоваться разработка набора моделей, соответствующих различным иерархическим уровням рассмотрения и функциональным разрезам деятельности системы. Такое описание на каждом уровне использует свой набор понятий и терминов.
В сравнении с натурным экспериментом математическое моделирование имеет следующие преимущества: Экономичность (в частности, сбережение ресурсов реальной системы). Возможность моделирования гипотетических, т.е. не реализованных в натуре объектов (прежде всего на разных этапах проектирования). Возможность реализации режимов, опасных или трудновоспроизводимых физически. Возможность изменения масштаба времени. Легкость многоаспектного анализа. Большая прогностическая сила вследствие возможности выявления общих закономерностей. Универсальность технического и программного обеспечения проводимой работы (ЭЦВМ, системы программирования и пакеты прикладных программ широкого назначения).
Математическое моделирование на определенных этапах исследования может сочетаться с натурным. Классическим примером такого подхода является исследование динамики летательного аппарата на комплексе из математической модели самого аппарата, воспроизводимой на аналоговой или цифровой ЭВМ, и макета реальной аппаратуры управления.
Практическое использование модели возможно лишь после тщательного ее исследования и настройки, в процессе которых необходимо решить задачи проверки адекватности, идентификации параметров, оценки значимости параметров и структурного преобразования.
Адекватность модели устанавливается проверкой для нее основных законов предметной области типа законов сохранения и сопоставлением результатов моделирования частных вариантов с известными для этих вариантов решениями.
Задачей идентификации является определение значений рабочих параметров модели по набору исходных данных, получаемому в результате наблюдения над реальной системой. При этом тип модели предполагается известным. Эта задача обычно ставится в форме минимизации функционала отклонения траектории модели от траектории исследуемой системы. На практике для ее решения традиционно применяются методы наименьших квадратов и наибольшего правдоподобия.
Выявление значимых параметров (и пренебрежение остальными) позволяет уменьшить размерность пространства параметров. Оно базируется на определении коэффициентов чувствительности (регрессии) выходных показателей по отношению к входным и представляет особую ценность при решении задач оптимизации.
Задачи преобразования модели обычно решаются после накопления некоторого опыта работы с ней. К ним относятся упрощение, усложнение, структурные изменения - в зависимости от результатов анализа этого опыта. Важным способом уменьшения размерности модели является рекомендуемый в теории подобия переход к обобщенным параметрам.
Особенности решения задач в переменных годографа
Оказывается, что в переменных годографа плоское потенциальное движение газа описывается линейным уравнением. Под переменными годографа обычно понимают переменные, связанные так или иначе со скоростью движения газа. Это могут быть либо компоненты и, v, либо модуль и угол скорости w, в, либо их некоторые комбинации. Ниже приводится уравнение для описания плоского потенциального движения газа, известное в литературе как уравнение С. А. Чаплыгина.
Введем новую переменную Ф = - р + хи + уи (25) Здесь р = р{х, у) потенциал скорости. Его полный дифференциал запишется так: dtp - udx + vdy = d(ux) - xdu + d(uy) - ydv (26) d $ = d(- p + их + uy) = xdu + ydv Переход от независимых переменных х, у к переменным u, v совершается с помощью преобразования Лежандра, даваемого функцией (25) Ф = Ф{и,и) и формулами ЭФ дФ .... ди до вытекающими из второй формулы (26). Удобнее от переменных и,и перейти к переменным w,$. Получим из (27) Уравнения газовой динамики приводятся к одному уравнению С. А. Чаплыгина: Э2Ф w2 д2Ф дФ л + - TTTZrT + w r = 0- (29) дв2 \-w2/a2 d2 dw Здесь a=a{w) из уравнения Бернулли. Уравнение (29) широко использовалось для исследования околозвуковых течений.
Решение задачи в переменных годографа состоит из двух этапов. 1. Определение функции Ф = Ф(и-, в), удовлетворяющей гра ничным условиям задачи. 2. Переход в физическую плоскость (w, в) = (х, у).
Следовательно, для получения решения исходной задачи нужно не только найти единственное решение Ф(№, в), но и осуществить однозначный переход в физическую плоскость (х, у). Необходимым условием такого перехода является следующее требование: якобиан Д = (1) не должен менять знак, проходя через нуль[43]. Введем определение: линия w=w(6) (или у = у(х)), в каждой точке которой Д = 0, называется предельной линией.
Покажем, что если при прохождении предельной линии определитель Д меняет знак, то решение в плоскости (х,у) становится комплексным.
Пусть Д = 0 на линии w=w( 6) а производная ду/д&\ Ф О т. е. поток в направлении в этом месте переменный. Тогда имеет место дв d{x,y)d{w,6) д{х,у) dz = 0 (2) ду W d(w,0) d(w,y) d(w,y) &w так как Д = 0. Следовательно, функцию х — x(yv,y), которая получится изх-x(w,в),у-y(w,9), путем исключения в, в окрестности линии w=wd0), npny=const можно представить так: 1 д2х х-х0 = (3) 2dw: (w-w0)" где вторая производная определяется, очевидно, на самой линии. Она отлична от нуля, когда якобиан и вместе с ним производная от д: в (2) изменяют свой знак при пересечении предельной линии. Из (3) видно, что по одну сторону от линии w=w0{&) функция w=w(x), при j =const будет комплексной.
Формула для Д получится, если производные подставить из формул замены переменных. После достаточно громоздких вычислений получим
Отсюда видно, что в дозвуковом потоке (w а) всегда Д 0. Предельные линии возможны только в сверхзвуковом потоке. Появление предельной линии означает, что в данном месте невозможно непрерывное течение. В потоке должны возникнуть ударные волны. Отметим, что в общем случае положение предельной линии не совпадает с положением ударной волны.
Остановимся на одном важном факте, вытекающем из условия положительности якобиана в дозвуковом потоке. Имеем dw = 0 (5) 1 d(6,w) d(0,w)d(x,w) дв ду X А д(х,у) d(x,w) д(х,у) дх При (w а) всегда Д 0, т.е. производные дв/дх\п и дв/ду\ имеют одинаковые знаки. Иначе говоря, если мы будем в физической плоскости двигаться вдоль линии w=wo=const так, что область w w0 остается справа от линии, то вектор скорости будет монотонно поворачиваться против часовой стрелки, так как угол в должен монотонно возрастать.
Этот результат принадлежит А. А. Никольскому и Г. И. Таганову. Он используется иногда при анализе трансзвуковых течений, так как остается справедливым и на звуковой линии -линии перехода от дозвукового течения к сверхзвуковому, на которой w=a (критическая скорость).
Методы численного решения задач математической физики
Понятия аппроксимации, устойчивости и сходимости создали необходимую базу широкого поиска эффективных разностных схем для решения задач математической физики. Алгоритмы решения задач с помощью конечно-разностных методов, как правило, представляют собой сочетание методов построения разностных аналогов задач и методов их решения. Поэтому прогресс в конструктивной теории конечно-разностных методов обязан взаимосогласованному развитию указанных двух направлений исследований.
Если попытаться просуммировать богатейший опыт в развитии конечно-разностных методов последних лет, то условно можно выделить следующие важные направления.
Построение разностных схем. Одно из таких направлений связано с разработкой методов построения консервативных разностных схем, основанных на законах сохранения, свойственных большинству физических процессов. Для конструирования консервативных разностных схем исходят из уравнений балансов, записанных для отдельной ячейки сеточной области, с последующим использованием квадратурных и интерполяционных формул. Построенные разностные уравнения после необходимых преобразований и суммирования по всем точкам сеточной области удовлетворяют дискретным аналогам интегральных законов сохранения.
Такие подходы рассмотрены О. А. Ладыженской [29], в работах которой построены разностные операторы для уравнений с разрывными коэффициентами, имеющие единый вид для любой внутренней точки области. Для обоснования алгоритма использовано понятие обобщенного решения и доказано, что решение разностной задачи образует некоторый функционал, переходящий при h - 0 в функционал дифференциальной задачи.
Консервативные разностные схемы сквозного счета в гидродинамике разработаны С. К. Годуновым [13,15], Лаксом и Вендрофом [84] на основе явных разностных аппроксимаций. Большое значение для решения задач гидродинамики имеет метод интегральных соотношений, предложенный А. А. Дородницыным [17,18] и развитый О. М. Белоцерковским, П. И. Пушкиным [4] и др., в котором использована частичная разностная аппроксимация уравнений, записанных в дивергентной форме, на основе метода прямых. Эти методы сыграли существенную роль в формировании общего взгляда на конструкцию разностных схем для квазилинейных уравнений. Интересные общие подходы к интегрированию уравнений гидродинамики также предложены в работах К. И. Бабенко, В. В. Русанова, Фромма [75], Кроули [70], В. Ф. Куропатенко.
В последнее время большое внимание уделяется построению решений задач математической физики высокого порядка точности.
Здесь в основном определились два направления. Первое связано с повышенным порядком аппроксимации разностных уравнений. Такие идеи рассмотрены в исследованиях А. А. Самарского, Н. Н. Яненко [56] и А. Н. Валиуллина [6], Митчелла [88] и др. Второе направление связано с построением решений на основе разностных уравнений сравнительно невысокого порядка на последовательности сеток с изменяющимися шагами. Эти методы получили название экстраполяции по Ричардсону и нашли отражение в исследованиях Е. А. Волкова [8,9], Фокса [74], В. В. Шайдурова [54] и др.
Методы построения разностных схем для уравнений эллиптического и параболического типов в классе разрывных коэффициентов разработаны на основе интегро-интерполяционного метода А. Н. Тихоновым, А. А. Самарским [44,45] и др.
Вариационные методы. В последние годы заметно повысился интерес к вариационным методам решения задач математической физики. Вариационные методы Ритца, Галёркина, Трефтца и др. давно заняли в вычислительной математике важное место. Особенно эффективны эти методы в тех задачах, где искомыми являются функционалы от решения. Оказалось, что уже при сравнительно невысоких приближениях функционалы получаются с большой точностью.
Численные схемы реализации алгоритма модели нахождения безударного профиля
При построении первой модели ставилась задача получить результат с минимальными затратами компьютерного времени поэтому использовались простейшие схемы и методы. Применялась разностная схема «простого креста» для уравнения и схемы численного интегрирования второго порядка для граничных условий, причем граничное условие на контуре реализовывалось в предположении достаточной тонкости профиля и близкому к параллельному Ох течению при помощи «сноса» граничного условия на профиле на ось х (нормаль п »(0,1)). Дня решения задачи Коши использовался метод характеристик на основе метода Эйлера численного интегрирования обыкновенных дифференциальных уравнений (смотри Глава 2 4.3), причем в качестве начальных условий взяты координаты положения точки на V профиле и угол наклона вектора скорости V к оси Ox or = arctg—, " х результате при достаточно малых числах Маха (порядка 0.675) был получен профиль, демонстрировавший практически безударное течение. Однако при увеличении числа Маха данная численная схема не давала результата, что выражалось в отсутствии сходимости метода характеристик. Программа расчета течения газа в канале методом фиктивного газа и расчета течения реального газа методом характеристик была написана в среде Digital Visual Fortran, который был выбран для решения поставленной задачи, как наиболее адекватное средство для проведения больших объемов вычислений с крупными массивами исходных данных. Однако на ряде этапов, (например, при нахождении исходной точки безударного профиля) оказалось необходимым подключить математический пакет (в данном случае Mathematica).
Графическое представление потока газа в канале, формы звуковой линии, и собственно, безударного профиля выполнено также в среде Mathematica.
Была разработана более гибкая и точная схема численной реализации модели. Введены разностные схемы более высокого порядка, вплоть до четвертого для точек внутри области, для вычисления скорости VVL для граничных условий на стенках канала. Удалось уйти от предположения о тонкости профиля, которое влечет соответствующую ошибку, перейдя к разностным эи уравнениям для производной по нормали _ дп
Предусмотрена возможность использования контура, задаваемого любой гладкой кривой. Вместо метода простых итераций был применен метод сопряженных градиентов. Он применялся для решения линейной системы, возникающей при линеаризации исходной нелинейной системы Ь(р(у)Уи = В, где V — скорость с потенциалом U за счет фиксации скорости, т. е. постоянной плотности. Так как метод сопряженных градиентов может быть применен для системы с положительно определенной матрицей, была разработана соответствующая подпрограмма умножения матрицы системы на транспонированную ей, причем эта подпрограмма сохраняет все преимущества ленточного характера матрицы системы. Разумеется прямое умножение на транспонированную матрицу неприемлемо в силу ее чрезвычайно больной размерности (п2хт2\ а также двойной нумерации элементов искомого вектора U. Используя ленточный характер операция применения такой матрицы сводится к скалярному умножению 5 - 16 (в зависимости от точности разностной схемы) векторов порядка яхт. Выбирается начальное приближение Е7„, находится матрица скорости V0 и полученная система L[p{V0)yU = В методом сопряженных градиентов решается с точностью до є по невязке системы в некоторой метрике (разумеется с приведением матрицы системы к симметричной положительно определенной). После достижения промежуточного результата Ux скорость вновь пересчитывалась на Vx и система L\p(V\)) U = В решалась с новой фиксированной плотностью и так далее до сходимости исходной нелинейной системы с точностью f, то есть получения Un, для которого невязка л/Ш/?(Уп)]/„-Я . Таким образом, метод сопряженных градиентов применяется не как метод точного решения линейной системы, что нереально вследствие ошибок округления и большой размерности системы, а как итерационный метод соответствующий по скорости сходимости геометрической прогрессии.
