Содержание к диссертации
Введение
Глава 1. Квазигазодинамические уравнения с источниками тепла 8
Вывод КГД уравнений при наличии внешних сил и тепловых источников . 9
Уравнение баланса энтропии 13
Выводы 17
Глава 2. Алгоритм решения КГД уравнений на треугольной сетке 19
Аппроксимация КГД уравнений 20
Аппроксимация производной по времени 20
Контрольный объем 21
Аппроксимация пространственных производных 23
Аппроксимация граничных условиях 23
Построение вычислительного алгоритма 24
Выбор шага по времени 27
Глава 3. Дозвуковое обтекание цилиндра. Дорожка Кармана . 30
Постановка задачи 31
Выбор параметра г 34
Моделирование течения 36
Выводы 41
Глава 4. Задача о сверхзвуковом обтекании пластинки 42
Постановка задачи 43
Сравнение методов на примере однородной стационарной задачи 47
Влияние правой границы на результаты моделирования 53
Влияние внешних сил 56
Влияние источников тепла 61
Заключение 63
Глава 5. Задача о структуре одномерной ударной волны 67
Математическая модель 68
Алгоритм решения 71
Результаты численного решения для аргона и азота 73
Заключение 80
Направления дальнейшего развития 81
Приложения 83
- Вывод КГД уравнений при наличии внешних сил и тепловых источников
- Аппроксимация производной по времени
- Моделирование течения
- Влияние правой границы на результаты моделирования
Введение к работе
Актуальными задачами численного моделирования газодинамических течений являются задачи с внешними источники энергии. К таким проблемам относятся, например, расчеты течений излучающего газа, исследование возможностей управления потоками с помощью энерговложения, расчеты активных сред в резонаторах газовых лазеров, задачи горения и многие другие практически важные вопросы.
Обычно такие задачи решаются с использованием системы уравнений Навье-Стокса, однако на этом пути существует целый ряд сложностей, например, необходимость введения дополнительных слагаем ых-регулярнзаторов для улучшения сходимости или использования специальных схем.
В последнее время на основе квазигазодинамических (КГД) уравнений были построены эффективные численные алгоритмы решения задач газовой динамики. КГД уравнения отличаются от уравнений Навье-Стокса дополнительными слагаемыми с малым параметром т размерности времени в качестве коэффициента. Эти дополнительные слагаемые имеют дис-сипативный характер, что демонстрируется справедливостью теоремы о неубывании полной термодинамической энтропии в замкнутом объеме, доказанной для этих уравнений. Диссипативные слагаемые выполняют роль регулярнзаторов и позволяют строить эффективные численные алгоритмы расчета нестационарных течений газа в дозвуковых и сверхзвуковых режимах.
Целью настоящей работы является обобщение КГД уравнений на течения в присутсвии внешних источников (или стоков) тепла, создание алгоритма на базе этих уравнений для расчета течений на неструктуриро-
ванных сетках, а также тестирование этого алгоритма и решение с его помощью практически полезной задачи об управлении аэродинамическими свойствами летательного аппарата.
Первая глава посвящена выводу КГД уравнений для идеального газа в присутствии источников тепла и внешних сил. При этом уравнения без источников уже были получены ранее другими авторами (см. [31]). Следуя тому же формализму, удается не только добавить слагаемые с Q в КГД систему, но и показать диссипативпость этих добавок при малых т. При этом в полученных выражениях легко прослеживается связь как с КГД системой в отсутствие энергоподвода, так и с уравнениями Навье-Стокса.
Во второй главе формулируется вычислительный алгоритм на базе КГД уравнений. При этом приоритетами являются простота, универсальность и эффективность в максимально широком классе задач. Для этого выбираются треугольные нерегулярные сетки и схема относительно невысокого порядка, однако даже при невысоком порядке апроксимации удается получать хороший результат за приемлемое время. Кроме того, описанный алгоритм позволяет моделировать как дозвуковые, так и сверхзвуковые течения при минимальном вмешательстве со стороны человека.
В третьей главе проводится тестирование алгоритма в отсутствие источников тепла на задаче об обтекании прямого кругового цилиндра. Эта задача хорошо исследована экспериментально и в настоящее время является довольно общепринятым тестом численного алгоритма на состоятельность. В определенном диапазоне чисел Рейиольдса за цилиндром образуется устойчивый периодический режим (так называемая "дорожка Кармана"), для периода колебаний в котором имеется эмпирически полученное выражение. Критерием состоятельности алгоритма служит совпадение наблюдаемого на модели периода колебаний с вычисленным по формуле.
В четвертой главе рассматривается модельная задача о воздействии электрического разряда на характер обтекания пластинки. При этом результаты имеют важное практическое значение для оценки параметров экс-
периментальной установки. Работа проводилась совместно с французской научной группой, в которой, помимо эксперимента, проводилось математическое моделирование методом Монте-Карло. Это дало возможность сравнить результаты алгоритма не только с экспериментальными данными, но и с результатами другого численного метода, полученными независимо.
КГД уравнения позволяют достаточно просто строить эффективные численные алгоритмы, однако для построения аналитического решения удобнее использовать исходные уравнения Навье-Стокса. В пятой главе рассмотрена одна из таких задач — задача о структуре фронта одномерной ударной волны. При помощи интегрирования системы Навье-Стокса и сведения к двум обыкновенным дифференциальным уравнениям задачу о фронте одномерной ударной волны становится возможным решать методом стрельбы по параметру. Традиционно считается, что получающаяся в уравнениях Навье-Стокса ширина волны отличается от экспериментально наблюдаемой в несколько раз. Однако, используя уточненное выражение для коэффициента вязкости и вводя вторую вязкость для многоатомного газа, результат удается значительно улучшить так, что отклонение от эксперимента не превосходит 30%. В пятой главе проводится расчет ширины фронта для ударной волны в аргоне и азоте, на примере последнего демонстрируется важность учета второй вязкости в уравнениях Навье-Стокса.
В заключении формулируются основные результаты работы и возможные направления дальнейших исследований.
В приложении 1 приведен краткий обзор существующих коммерческих и свободных программ для построения сеток и работы с ними. Результаты настоящей работы получены с использованием пакетов Comsol Multiphysics и MATLAB, в приложении даны некоторые пояснения относительно построения и экспорта сеток в этих пакетах.
Приложение 2 содержит наиболее важные фрагменты программного кода, использовавшегося при решении КГД системы на треугольной сетке. Программа написана на языке Java версии 5.0.
Обобщение КГД уравнений на случай неидеальных газов важно для практического применения. Некоторые полученные в этом направлении результаты вынесены в приложение 3. В общем случае удается построить КГД систему, однако диссипативпость добавок в общем случае доказать пока не удается. Есть уверенность, что для газа Ван-дер-Ваальса вдали от критической точки слагаемые имеют диссипативный характер, однако точные границы применимости КГД уравнений определены не были.
Автор выражает благодарность Татьяне Геннадьевне Елизаровой, под руководством которой была выполнена настоящая работа, и Юрию Владимировичу Шеретову за помощь в интерпретации полученных результатов.
Вывод КГД уравнений при наличии внешних сил и тепловых источников
К квазигазодинамическим уравнениям можно относиться двояко. С одной стороны, они представляют собой способ построения регуляризаторов специального вида для эффективного численного решения задач динамики вязкого газа. И с этого они и начинались. С другой стороны, в работах [10, 11, 28, 9, 7] вывод соответствующих регуляризаторов в отсутствие внешних сил и источников тепла был произведен на основе уравнения Больцмана для некоторой кинетической модели. С этой точки зрения квазигазодинамические уравнения можно рассматривать как естественное обобщение уравнений Навье-Стокса, а появляющиеся дополнительные слагаемые — как добавки, имеющие вполне определенный физический смысл. Данный подход становится еще более привлекательным, если заметить, что для КГД уравнений доказана теорема о балансе энтропии, а сами они естественным образом записываются в виде законов сохранения.
В книге [31] изложен альтернативный "эмпирический" вывод КГД уравнений, менее строгий, но применимый и в случае наличия внешних сил. Однако выведенные в указанной работе уравнения не охватывают задачи с внешними источниками тепла. В данной главе мы воспользуемся этой же методикой для обобщения КГД уравнений на системы с внешними источниками тепла, а также выведем соответствующее уравнение баланса энтропии, чтобы показать сохранение диссипативного характера полученных добавок.
Вывод КГД уравнений при наличии внешних сил и тепловых источников
Запишем интегральные законы сохранения для малого фиксированного объема VB конечно-разностном виде, заменяя производные по времени разностным отношением для моментов времени t и t + At, где At - малый промежуток времени. Тогда законы сохранения массы, импульса и полной энергии можно представить в виде, соответственно: / d3 + j/7/ t , =0, (1.1) Хтл+ll/ # +lis = Jp F\ Px + Jjllbsijdoj, (1.2) р(й /2 + е)-р(и /2 + є) , ff [(и У t tp _ d3x + ff p u i[ r- + e + )dai+ Jv A Jh 4 2 p + JJ lNSidi = J P Fid x + JJ Щду« і + J Q d x- (L3)
Здесь и далее использованы обычные обозначения: р - удельная плотность, щ - компонента скорости, р - давление, F{ - компонента внешней силы, є - удельная внутренняя энергия, Q - мощность тепловых источников. В формулах подразумевается суммирование по повторяющимся индексам. Вектор теплового потока и тензор вязких напряжений вычисляются как: qmi = -x — T, (1.4) ( д д 2г д \
Пад = 7]\-дх- +ЩЩ- 3Sijdx kUk) (L5) ц\\ к- коэффициенты вязкости и теплопроводности. Величинами со звездочками в (1.1)- (1.3) обозначены значения газодинамических параметров в момент времени t t t + At. Обозначим At/2 = т, и определим параметры газа в средней точке t = t , ограничиваясь первым порядком малости по т: Р = Р + тр (1.6) u i=Ui + TW (L7) Р = Р + т% (L8)
Подставим эти выражения в формулы (1.1)-(1.3). Величины Fi и Q считаем мало меняющимися за время т, поэтому звездочки возле них просто отбросим. Заметим также, что фактически использованное нами соотношение ф-ф д ± Atd2 ± верно с точностью до членов второго порядка по At, поэтому слагаемыми порядка 0(т2) естественно пренебречь. Возвращаясь вновь к дифференциальному виду производной по времени в интегралах по объему, получим / —р(Рх+ (риі + т — риЛ(Іаі = 0, (1.10) / — рщ(і3х + II (pui + т — рщ j Ujdaj + // три{ — щ Ь}+ +JL{P+Titp) =Xip+TTtp)3 +п» a-») l JiTpUi {u4tt + Fte+pFt-p + -pWtp)dai + tiq,iSidai = = / l рщ + T — рщ J Fid3x + // ITNS ІjUidaj + I QcPx. (1.12)
Следуя методике, изложенной в [31], рассмотрим систему уравнений Эйлера (при наличии внешних сил и источников тепла) №+ _№tl.+ _p = , (L14) д fu2 \ д fu2 p\ г, Л /, 1Гч 91Р{-2+е) + Г1{-1 + г + р)=Ри Р1 + С3 (U5) заметив, что для идеального политропного газа с уравнением состояния p = pRT, є = -Т (1.16) 7-1 из нее следуют тождества д 1 5 115 п ,., _ — - + «,-- -w/ = 0, 1.17 at р axi р р ах{ _Ui + u.__Uj + __p_Fj = 0, (1.18) У axi р axi р p + B;.p + 7p.Mi_(7_l)(? = 0. (1.20)
Справедливость тождеств (1.17)-(1.20) проверяется непосредственной подстановкой в них выражений для производных по времени из уравнений Эйлера (1.13)-(1.15) с последующим приведением подобных слагаемых.
Будем считать, что в нулевом по г приближении для нашего газа справедливы уравнения Эйлера. Из физических соображений это допустимо, если г — характерный временной параметр газа, в качестве которого обычно выбирается время свободного пробега частиц или пропорциональное ему характерное время максвеллизации системы. При этом вслед за отброшенными слагаемыми порядка 0(т2) мы пренебрегаем и величинами 0{тг)).
Аппроксимация производной по времени
В качестве неизвестных выберем р, р и щ. Заметим, что каждое из уравнений КГД системы имеет вид dt дхі При этом величины ВІ и С выражаются через значения функций на предыдущем слое и производные от этих значений. Воспользуемся явной по вре мени схемой, то есть будем строить решение на следующем временном слое по формуле
Порядок решения понятен из вида КГД системы: 1. сначала из первого уравнения получаем значение р на новом слое, 2. затем, используя полученное р, из второго уравнения находим компоненты Ui 3. и на последнем этапе, подставляя р и щ, находим новое значение р, вычислив є и воспользовавшись уравнением состояния.
Эта процедура применима, конечно, только к внутренним узлам сетки. Для определения значений функций на границе необходимо воспользоваться граничными условиями.
Контрольный объем
Для аппроксимации правой части (2.1) воспользуемся методом конечного объема, подробно описанным в [5]. Заметим, что выражения, обозначенные условно буквами В{ и С, содержат только первые пространственные производные. Вокруг каждого внутреннего узла М сетки построим контрольный объем, как это показано на рисунке 2.1. Точки Mk — соседние узлы сетки (не обязательно внутренние). Точки Р — центры пересечения медиан1 треугольников, составляющих сетку. Точки Nk — середины ребер, соединяющих соседние узлы сетки, а заштрихованный многоугольник Рк образует контрольный объем. Точки iVJt для краткости в программе называются полуцелыми точками (из-за аналогии с одномерным случаем, когда значения пространственных производных вычисляются в полуцелых точках). J3TO ОДИН ИЗ ВОЗМОЖНЫХ способов построения контрольного объема. Он обеспечивает высокую точность на достаточно регулярной сетке [1], [22]. g-Bi й 4 В,п, fl а ] B,(Nt)n,(Nt) Д(р,р»_,. (2.2) Здесь S - площадь контрольного объёма-многоугольника, щ — координаты внешней нормали, Діркд._, — расстояние между центрами соседних треугольников.
На рисунке 2.1 обозначена точка N!, - середина отрезка Р-\р2- Ее отличие от точки JVa хорошо заметно. В общем случае она не лежит на отрезке ММ%, Мы отождествляем точки N% (к ней относится среднее арифметическое параметров газа в точках М и М2). Щ (используется в интеграле 2.2) и точку пересения отрезков Р1Р2 и MM-L (в ней вычисляются пространственные производные).
Аппроксимация пространственных производных
Рассмотрим произвольный четырехугольник, изображенный на рисунке 2.2. На рисунке 2.1 ему соответствует, например, четырехугольник MP\MiPi-В работе [5] выводится следующее выражение, аппроксимирующее с первым порядком значения производных / в точке К: (2.3) (2.4) df (/2 - /0) (уз - У\) + (/і - /з) (У2 - уо) дх 2S df (/2 - /0) (si - 2) + (/1 - /3) (so - х2) ду 25
При этом подразумевается, что нумерация вершин производится в положительном направлении обхода внутренности (против часовой стрелки). Xfa У к координаты вершин четырехугольника.
Аппроксимация граничных условиях
Как уже говорилось, схема с контрольным объёмом не охватывает граничные точки области. Соответствующие конечно-разностные выражения для граничных условий имеют следующий вид [5]. Для условий Дирихле /г = /о = f(Mt) = /о, для условия, содержащего производную, :««: = /0 = /(мо = -Е /(Лй) + /о Efc Ымм N + PZkAlMM
Здесь суммирование ведется по всем соседним внутренним узлам сетки, Ahh-Mi — расстояние между точками М& и МІ, п — внешняя нормаль к границе области. Для /? = О формула аппроксимирует условие второго рода, для (3 0 — третьего.
Построение вычислительного алгоритма
Мы не будем выписывать окончательное выражение для разностной схемы. Во-первых, оно уже выписано в работе [5]. Во-вторых, при построении программы окончательное выражение не использовалось. Вместо этого вычисления строились поэтапно, а промежуточные значения сохранялись и использовались повторно для улучшения читаемости программы и уменьшения количества арифметических операций.
Сначала кратко опишем необходимые структуры данных: 1. массив координат узлов сетки п —+ (х, у), 2. массив треугольников п — (прі,пр2,прз), 3. массив координат середин треугольников п — (хс,ус), 4. массив иолу целых точек (JVjt на рисунке 2.1), нумерующий пары соседних треугольников, п — (па,Щ2) и таблица обратного соответствия {гкі,па)- п, 5. массив списков треугольников, имеющих данную вершину, пр\ — {щк}, причем треугольники в списках должны быть отсортированы так, чтобы их центры обходились против часовой стрелки.
Моделирование течения
Итак, рассмотрим результаты работы вычислительного алгоритма. Стадии развития течения изображены па рисунке 3.8. Сначала за цилиндром образуется симметричная отрывная область, состоящая из двух вихрей, вращающихся в разных направлениях. Далее область отрыва удлиняется, при этом увеличиваются и вихри. При малых числах Рейнольдса на этом этапе решение выходит па стационар, в пашем же случае вихри продолжают расти, пока один из них не отрывается. При этом область отрыт теряет симметрию. Дорожка продолжает удлиняться, на месте оторвавшегося вихря зарождается новый меньшего размера. Постепенно устанавливается периодическое течение: за цилиндром образуется изгибающаяся дорожка из отрывающихся вихрей. Расчитаннып по формуле (3.3) период колебаний дорожки равен Тсхр = 0.053с. (3.17)
Чтобы определить период в нашем случае, построим одномерный график зависимости одного из параметров задачи (мы выбрали вертикальную компоненту скорости) в фиксированной точке от времени. Получившиеся графики в разных точках течения представлены на рисунке 3.9.
В правой части графика колебания установились и можно определять период. Мы рассмотрели расстояние между тремя последними максимумами, получившаяся в результате величина периода колебаний равна (в независимости от координаты точки наблюдения)
Tmodel = 0-059c- (ЗЛ8)
Результат превосходит предсказанное теоретически значение. Это отклонение вполне понятно. В наших теоретических расчетах в формулу для периода мы подставили в качестве D диаметр цилиндра. В то же время, даже в работах самого Кармана под D подразумевается некоторая эффективная толщина препятствия. Возможность увеличения "эффективного" диаметра цилиндра в наших расчетах изображена на рисунке 3.10. Штриховой красной линией проведена окружность большего диаметра, которая, однако, не содержит внутри себя ни одного внутреннего узла сетки. В пределах порядка аппроксимации граничных условий, а он у нас первый, диаметр изображенной окружности и исходный диаметр цилиндра равнозначны2.
Проведем моделирование с использованием более подробной сетки (изображена на рисунке 3.11). Количество узлов и треугольников указано на рисунке, рекомендуемый размер треугольника при генерации сетки был уменьшен на 20%. Качественно характер колебаний не изменился (визуально он почти не отличается от изображенного на рисунке 3.9), однако получившееся значение периода 7 12 = 0.057 с (3-19) уже ближе к экспериментальному. Заметим, что несмотря на увеличение общего количества треугольников размер треугольников вблизи границы цилиндра сильно не изменился. Однако увеличилась плотность узлов сетки в области образования и развития вихрей и этого оказалось достаточным для улучшения результата.
Выводы
Была рассмотрена модельная нестационарная задача, результат которой хорошо известен и служит хорошим критерием состоятельности алгоритма. На ее примере было изучено влияние величины параметра г на характеристики получающегося решения, продемонстрирован способ введения регуляризации на основании КГД подхода.
Надо отметить, что КГД уравнения позволяют совершенно естественным образом вводить регуляризацию на сильно неравномерных треугольных сетках, на которых для уравнений Навье-Стокса требуется гораздо более сложное исследование.
Полученная величина периода колебаний согласуется с теоретической, хотя и отличается от нее в большую сторону. Причину расхождения нужно искать в конечном размере треугольника и вытекающей из этого неизбежной ошибке аппроксимации. При этом моделирование на более подробной сетке позволило улучшить результат. Подобный эффект нужно учитывать при построении сеток для промышленных и инженерных задач, когда 10%-ное расхождение результата с теоретической величиной имеет серьезные последствия.
Идея искусственной ионизации газа с целью изменения обтекательиых свойств поверхности была выдвинута больше 10 лет назад, но все еще исследуется. Многообещающей областью исследования является уменьшение силы трения, действующей на летательный объект с помощью электрогидродинамического соленоида, представляющего собой набор электродов, создающих электрический разряд вокруг тела.
Научной группой из Орлеана3 ставился эксперимент в аэродинамической трубе по воздействию постоянного электрического поля на характеристики обтекания пластинки с прикрепленными к ней электродами. Схематически экспериментальная установка изображена на рисунке 4.12.
К неподвижной пластинке прикреплены два электрода, к которым подведено постоянное напряжение. Вся конструкция обтекается сверхзвуковым потоком разреженного воздуха, измеряемая величина — сила трения, действующая на пластинку. Задача математического моделирования заключалась в сначала качественном, а потом, возможно, и количественном предсказании поведения силы трения в зависимости от величины электрического поля. Параллельно с КГД подходом независимо производилось моделирование методом Монте-Карло, что позволило сравнить полученные результаты.
Влияние правой границы на результаты моделирования
Моделируемый эксперимент проводится в лаборатории большого размера, по сравнению с которым пластинка мала. Естественной моделью было бы обтекание в бесконечном пространстве, однако численный эксперимент требует введения границ и постановки граничных условий.
Рассматриваемая геометрия области была выбрана изначально для решения методом Монте-Карло, а потом без изменений использована и при решении рассматриваемым алгоритмом. Результаты разных методов хорошо согласуются между собой, но возникает естественный вопрос, насколько выбранный размер области достаточен для моделирования происходящих явлений. При этом наибольшие опасения возникают но поводу правой границы, так как условия на ней носят наиболее искусственный характер.
Отодвинем правую границу и построим новую сетку так, чтобы ни одно из ее ребер не пересекало прежнюю правую границу х — 0.12. Сетка изображена на рисунке 4.19. Преимущество такой сетки понятно: мы можем оценить влияние правой границы, напрямую сравнивая результаты без необходимости интерполировать значения. Получившиеся результаты визуально не отличимы на пересечении областей, так что наложим на распределения линии уровней разностей значений.
На рисунке 4.20 изображены получившиеся двумерные графики. Видно, что влияние границы ограничивается, в основном, несколькими прилегающими к ней треугольниками. Максимальные значения отклонений соста ВИЛИ Д« шах = 190 м/с, Дртах = 4.1 Па, Артах = 0.9 10"4 кг/м3, (4.17) относительные отклонения тах = 0 31j Айвах = Q Дртах = Q 3? (4Д8) итах Ртах Ртах
Линии уровня проведены через равные промежутки, вне области, ограниченной линиями относительная ошибка составляет менее 3%. При этом, как видно из графиков, максимальных и близких к ним значений отклонение достигает лишь в нескольких прилегающих к границе узлах сетки.
Влияние внешних сил
Как уже говорилось, мы рассмотрим случай, когда электрическое поле однородно и постоянно внутри эффективного объема (обозначенного прямо — угольником). Вектор Е направлен параллельно оси х в положительном направлении (вправо). Так как в формуле для внешней силы, напряженность поля умножается на коэффициент ионизации а, в качестве параметра задачи выберем это произведение.
Было сделано несколько пробных вычислений, после которых стало понятно, при каких значениях параметров течение практически не изменяется, а при каких не остается и следов от первоначального (с выключенным полем) режима. Наибольший интерес представляют промежуточные значения, которые и были выбраны для расчетов. В таблице 4.2 представлены выбранные значения.
Электрическое поле ускоряет молекулы воздуха, их скорость возрастает. Вблизи левой границы поля образуется дефицит частиц: из-за большей скорости исходящий поток справа больше входящего потока слева, в этой области давление и плотность падает. Так как молекулы движутся быстрее, относительно далекие слои воздуха не успевают прогреться, их температура меньше температуры в отсутствие поля, а из-за разрежения, даже меньше температуры невозмущенного потока. Распределение скорости, температуры и плотности изображено на рисунке 4.21.
Из-за простоты рассматриваемой модели результат имеет простую и красивую физическую интерпретацию, что особенно важно, так как исследованиями этой же задачи занималась группа физиков-экспериментаторов.
Профили скоростей для различных нанряженностей поля изображены на рисунке 4.22, а на рисунке 4.23 изображена подсчитанная зависимость силы трения вдоль пластинки от приложенного поля.
На графике силы трения хорошо виден дискретный характер алгоритма: при дифференцировании скорости их график теряет свою гладкость. Тем не менее, этого достаточно для вычисления суммарной силы трения и изучения общей закономерности: сила трения возрастает при увеличении ускоряющего поля.
