Содержание к диссертации
Введение
1. Методика численного интегрирования системы уравнений Навье - Стокса 15
1.1. Основные уравнения. Схема расщепления . 15
1.2. Аппроксимация конвективных членов 18
1.3. Учет эффектов вязкости 21
1.4. Граничные условия на поверхности обтекаемого тела 23
1.5. Граничные условия на внешних границах расчетной области 25
1.6. Приведение полного потока жидкости через внешние границы расчетной области к нулю 28
1.7. Решение уравнения Пуассона. Граничные условия для давления 30
1.8. Исследование устойчивости и аппроксимации . 34
2. Результаты численного моделирования 37
2.1. Краткая характеристика созданного комплекса программ 37
2.2. Течение в кубической каверне 38
2.3. Пространственное обтекание тонкой прямоугольной пластины под углом атаки 48
2.4. Перспективность разработанного алгоритма . 63
3. Моделирование нестационарных режимов обтекания неравномерно движущихся тел 65
3.1. Предварительные замечания 65
3.2. Основные уравнения. Схема расщепления . 67
3.3. Модификация методики постановки граничных условий 69
3.4. Колебания пластины в покоящейся жидкости . 72
3.5. Колебания пластины в потоке жидкости . 77
3.6. Поворот пластины в потоке жидкости. Плоская и пространственная задачи 84
Заключение 97
Литература 100
- Основные уравнения. Схема расщепления .
- Краткая характеристика созданного комплекса программ
- Предварительные замечания
Введение к работе
Бурное развитие вычислительной техники и расширение возможностей современных ЭВМ с одной стороны, совершенствование существующих и разработка новых более экономичных алгоритмов с другой привели к тому, что появился новый метод исследования - численный эксперимент [I]. Численный эксперимент должен базироваться на гармоничном сочетании основных аспектов научной проблемы. К таковым относятся: четкое представление о физических особенностях явления; рациональная математическая постановка задачи, позволяющая реализовать экономный численный метод решения; приспособленность всего математического описания явления и численного метода к возможностям ЭВМ [ 2 ]. Численный эксперимент в ряде случаев может заменить физический или натурный эксперименты, проведение которых является либо практически невозможным, либо дорогостоящим.
Все элементы технологической цепочки численного моделирования физических процессов, состоящей из обоснования физико-математической модели, выбора или разработки численного метода и создания работающего комплекса программ, существенно определяют экономичность вычислительного процесса в целом. Несмотря на существующий прогресс в этой области, задача численного моделирования сложных физических явлений еще очень далека от своего завершения.
Моделирование на ЭШ пространственных задач динамики вязкой жидкости и газа на основе уравнений Навье-Стокса является сравнительно молодым, но интенсивно развивающимся направлением вычислительной механики. По оценкам советских и американ- ских специалистов для полного.решения указанных задач требуются мощные суперкомпьютеры с быстродействием до 10 млрд. операций в секунду и практически неограниченной памятью [3]. Стало общепризнанным, что такая производительность может быть достигнута только на ЭВМ новой архитектуры с использованием параллельных процессоров. Однако, ожидать, пока такие машины войдут в практику, было бы непростительной потерей времени. Отработка методики, подбор наиболее перспективных разностных алгоритмов, получение первых приближенных результатов возможны и даже необходимы уже сейчас. Необходимы потому, что с точки зрения современной вычислительной механики, разностная схема является не только аппроксимацией наперед заданного дифференциального, интегро-дифференциального или интегрального уравнения, а представляет собой самостоятельную математическую модель [ 3].
Этой тематике и посвящена настоящая работа. Разработан и реализован на практике алгоритм численного моделирования пространственных течений несжимаемой вязкой жидкости. Особое внимание уделено нестационарным задачам, в которых рассматривается неравномерное движение тел в среде. Моделирование проводится путем прямого численного интегрирования системы уравнений Навье-Стокса, записанных в естественных переменных (скорость V и давление р ) по методу расщепления. В связи с моделированием пространственных течений, в классическую схему метода [4-7] были внесены следующие изменения [ 8-Ю ] : для аппроксимации конвективных членов использованы ориентированные разности, разработана новая эффективная методика постановки граничных условий на внешних границах расчетной области, для решения уравнения Пуассона использован более современный попеременно-треугольный метод ГЦ].
Проведено численное моделирование и получены основные характеристики следующих гидродинамических задач: пространственного обтекания тонкой прямоугольной пластины под углом атаки; течения вблизи колеблющейся пластины в плоском потоке и покоящейся жидкости; плоского и пространственного обтекания поворачивающейся пластины.
В настоящее время существует большое количество методов численного интегрирования системы уравнений Навье-Сток-са для несжимаемой вязкой жидкости. Однако, большая часть этих методов была разработана применительно к решению двумерных задач для уравнений относительно функции тока УР и вихря и) [12-15]. Это связано с тем, что при моделировании плоских или осесимметричных течений (",«*>)-система оказывается, как правило, более экономичной, чем (V,p)~cncTe-ма [12]. Для пространственных течений положение существенно изменяется. Более экономичной оказывается уже (V,p)--система, где на каждом временном слое нужно решать только одно уравнение Пуассона, а не три, как это имеет место для трехмерного аналога (У,С*))-системы [16-17].
Основная трудность при численном интегрировании системы уравнений Навье-Стокса, записанных в естественных переменных, связана с расчетом поля давления. Один из путей преодоления указанной трудности был разработан на основе введения "искусственной сжимаемости" [18-19], когда вместо уравнения неразрывности использовалось одно из уравнений ^ї/ЬЬ + «Ltf Va С , ИЛИ
Такой подход применяется и в настоящее время при моделировании пространственных течений [20-211, хотя его ограниченность очевидна и связана с тем, что уравнение неразрывности аппроксимируется точно только в режиме установления.
Значительный успех был достигнут благодаря разработке и совершенствованию методов, типа MAC и его различных модификаций [ 22-26 ]. Методам этого типа присуща специфическая конечно-разностная схема и специфическая структура ячейки: давление определяется в центре ячейки, а компоненты скорости - в центрах соответствующих сторон. В этих методах в рассмотрение вводятся частицы-маркеры, которые не обладают массой и переносятся со скоростью конвекции. Непосредственно в вычислениях частицы-маркеры не участвуют и во внутренних точках вообще не оказывают обратного влияния на поток. Однако, имея ввиду, что метод MAC интенсивно применялся для моделирования течений несжимаемой вязкой жидкости со свободной поверхностью, отметим, что форма свободной поверхности, не известная априори, определялась в процессе счета по положению маркеров.
В методе расщепления, разработанном на основе метода MAC, была сохранена специфическая структура ячейки, использована схема расщепленияj лежащая в основе методов РІС [22] и Чорина [18-19].
Первая глава диссертации посвящена описанию- алгоритма численного интегрирования системы уравнений Навье-Стокса. В этой главе приведены основные уравнения и схема расщепления, изложены способы конечно-разностной аппроксимации конвективных и диссипативных членов. Особое внимание уде- лено постановке граничных условий на внешних границах расчетной области. Предложен универсальный трехшаговыи алгоритм, пригодный для описания границ любого типа (входной, боковой и выходной). На первом шаге нормальные составляющие скорости определяются по методике "открытых граничных условий" (0ІУ) [27 ], суть которых заключается в разностной аппроксимации на границе условия lCbj/г^.) + с (Э5/Эп.яг =о, где J - исследуемая функция, С - характерная скорость переноса возмущений по сетке, п. - нормаль к границе Г .
На втором шаге учитывается невозмущенная скорость набегающего потока. На третьем - приводится к нулю полный поток жидкости через внешние границы области интегрирования. Необходимость включения в схему этого шага связана с тем, что граничные условия для давления приводятся к разностному аналогу однородных условий вида
Сдр/ЭЮ г = О.
Полученная задача Неймана для уравнения Пуассона, как известно, имеет решение только в том случае, если
5І5 die у ль = & ^. af = о. я. $ где Q - область интегрирования, S - поверхность, ограничивающая данную область.
В [ 28 ], чтобы избежать необходимости приведения к нулю полного потока жидкости через внешние границы, предлагается использовать для давления неявный аналог "открытых граничных условий". Получаемая при этом третья краевая задача сводится к задаче Дирихле на расширенной области. Последняя, как известно, имеет решение всегда и решается обычно проще и быст- рее, чем задача Неймана. Однако, несмотря на указанные преимущества, возможность использования 0ІУ для давления представляется весьма проблематичной. Связано это с тем, что при аппроксимации градиента давления на верхнем временном слое полностью теряется смысл скорости С как характерной скорости переноса возмущений. Возникают определенные сложности с вычислением указанной скорости, и, в результате, существенно снижается эффективность предлагаемого подхода. Кроме того, ОГУ для давления в совокупности с 017 для скорости приводят к одновременному неявному заданию на внешней границе значений скорости и давления, а при этом может нарушаться соотношение Вернули.
Как отмечалось выше, давление в методе расщепления определяется путем решения уравнения Пуассона. В разработанном алгоритме для решения этого уравнения был выбран попеременно--треугольный метод (ЇЇТМ) [III. Основное преимущество ПТМ заключается в том, что итерационный оператор расщепляется на произведение двух операторов вне зависимости от размерности задачи. Использование ПТМ в трехслойной модификации с выбором итерационных параметров по методу сопряженных градиентов позволяет значительно снизить необходимое число итераций при незначительном увеличении объема вычислений, по сравнению с двуслойной схемой. Например, на поле до 10 ячеек для умень- шения квадрата невязки в 10 раз максимальное число итераций по давлению для трехслойной схемы не превышает 45-50. Дяя двуслойной при тех же начальных условиях - примерно на порядок больше, а метод верхней релаксации требует до 10 тысяч итераций.
Проведенное исследование использованной разностной схе- мы показало, что условие устойчивости первого этапа выполняется, если шаг по времени t не превышает некоторого максимального значения X мах » а второй и третий этапы безусловно устойчивы. Схема аппроксимирует исходные уравнения Навье-Стокса с первым по X и вторым по 'iv порядками точности ( 4с-максимальный из шагов по пространственным направлениям).
Предложена методика, с помощью которой можно построить схему второго порядка точности по X .
Вторая глава диссертации посвящена описанию результатов численных расчетов. Для проверки разработанного алгоритма была выбрана задача о течении несжимаемой вязкой жидкости в кубической каверне с движущейся верхней крышкой. Получено хорошее соответствие результатов расчетов с данными других авторов [29-321 как для квадратной, так и для кубической каверн.
Большой научный и практический интерес представляет задача пространственного обтекания тонкой прямоугольной пластины под углом атаки. Одно из возможных приложений связано с рассмотрением пластины как упрощенной модели крыла конечного размаха.
Теория крыла конечного размаха была заложена в работах Прантля и получила дальнейшее широкое развитие как в нашей стране, так и за рубежом [33-36]. Теоретические исследования, однако, ограничены из-за большой сложности возникающей задачи. Значительный прогресс в данной области был достигнут за счет использования ЭВМ. Особенно следует выделить широкий цикл работ, связанных с моделированием по методу дискретных вихрей [2,37-39]. С помощью указанного метода удается получить динамические характеристики не только тонких крыльев простой геометрической формы, но и сложных крыльев с изменяю- - II - щейся геометрией и даже летательных аппаратов в целом. Он дает возможность исследовать процессы сворачивания вихревой пелены, ее разрушения и формирования ближнего следа.
Математическое обоснование метода дискретных вихрей базируется на приближении идеальной несжимаемой жидкости. Задача моделирования обтекания тонкой несущей поверхности (системы поверхностей) сводится к решению уравнения Лапласа с соответствующими граничными условиями. Формирование вихревого следа во времени описывается дифференциальными уравнениями движения свободных вихрей по траекториям жидких частиц. Напряженность вихрей не изменяется со временем, а меняется лишь их положение в пространстве.
При реализации на ЭВМ непрерывные вихревые слои, которыми моделируются несущие поверхности и их следы, заменяются системами дискретных вихрей. Замыкание задачи проводится с помощью граничных условий, которые выполняются в конечном числе контрольных точек на несущих поверхностях. Неизвестные циркуляции дискретных вихрей находятся из решения системы линейных алгебраических уравнений.
Метод дискретных вихрей экономичен, поскольку в рассмотрение вводится только область с ненулевой завихренностью. Он дает возможность моделировать пространственное обтекание тел сложной формы как в стационарном, так и нестационарном режимах. Моделирование, однако, проводится в рамках приближения идеальной жидкости на основе линейных уравнений, а нелинейные эффекты учитываются только при постановке, граничных условий. Положение точек отрыва должно быть известно заранее.
Описанный в диссертации метод моделирования пространственных течений несжимаемой вязкой жидкости основан на прямом численном интегрировании системы уравнений Навье-Стокса. Такой подход представляется наиболее общим и в настоящее время интенсивно развивается. В работах [5,20,21,29,40] предлагаются различные варианты разностных схем, разработанных на основе методов МАС,' РІС и Чорина. Широкое распространение для моделирования плоских и осесимметричных течений получили методы конечных и граничных элементов. В настоящее время изучаются возможности использования этих методов для пространственных течений [41-45]. Изучаются также возможности построения гибридных разностных схем, учитывающих априорную информацию о поведении решения [46-48].
Возвращаясь к задаче о пластине, отметим, что мощности доступных ЭВМ весьма ограничены. Максимально возможное реальное количество узлов сетки составляет пока несколько десятков тысяч. Для адекватного моделирования локальной структуры те-, чения во всей области интегрирования, включая зоны пограничного слоя и вязкого следа за пластиной (даже при не очень больших числах Рейнольдса порядка 1000) общее количество узлов сетки следовало бы увеличить на 3 - 4 порядка. Такое требование пока невыполнимо. Поэтому, представленные во второй главе диссертации результаты -расчетов дают лишь общую картину течения, которая характеризуется сходом вихревой пелены с острых кромок пластины с образованием вихревых жгутов. Полученные результаты соответствуют физической сущности задачи и согласуются с аналитическими расчетами по асимптотической теории [49].
В третьей главе диссертации рассматриваются нестационарные режимы обтекания неравномерно движущихся тел несжимаемой вязкой жидкостью. Задачи, возникающие в связи с изучением та- - ІЗ - ких режимов привлекает большое внимание специалистов аэро- и гидродинамики. Подробное изложение этих задач дано в [503. Это перегрузки крыльев и хвостового оперения самолетов, совершающих быстрые маневры, флаттер крыла или гребного винта и т.д.
Особый интерес к нестационарному движению профиля или пластины в жидкости (или газе) связан с изучением машущего полета птиц и насекомых [ 51-53 ]. Проведенные экспериментальные исследования [ 54-57 ] показали неприменимость квазистационарного приближения и необходимость учета нестационарных эффектов. Полученные экспериментальные значения подъемной силы при быстром изменении угла атаки значительно отличались от рассчитанных в рамках квазистационарного приближения. Экспериментальные данные были подтверждены и дополнены расчетами на ЭВМ [58-61].
В настоящей работе продемонстрирована возможность моделирования указанных задач с помощью метода расщепления. Учет эффектов, связанных с нестационарностью, привел только к необходимости некоторой перестройки алгоритма и модификации методики постановки граничных условий на внешних границах расчетной области. В диссертации представлены результаты расчетов следующих гидродинамических задач: колебаний пластины в покоящейся жидкости; колебаний пластины в потоке жидкости; плоского и пространственного обтеканий поворачивающейся пластины. Для второй и третьей задач . проведено исследование нестационарных значений коэффициента подъемной силы.
Отметим, что значения коэффициента кинематической вязкости "V , заданные при моделировании указанных задач, соответствовали числам Рейнольдса порядка 10 . Однако, это соответст- вие чисто формальное. Здесь, так же как и в задаче пространственного обтекания прямоугольной пластины под постоянным углом атаки, полученные результаты дают лишь общую картину течений. Локальная структура этих течений определяется схемной вязкостью.
Основные результаты диссертации изложены в работах [8-101, докладывались и обсуждались: - на научных семинарах лаборатории вычислительной физики ВЦ АН СССР под руководством академика О.М.Белоцерковского и доктора физ.-мат. наук, профессора А.И.Толстых; '- на IX Всесоюзной школе-семинаре по численным методам механики вязкой жидкости (Ольгино, 1982); на ХХУП конференции МФТИ (Долгопрудный, 1982); на Всесоюзной школе-семинаре "Методы гидрофизических исследований" (Солнечногорск, 1983); на X Всесоюзной школе-семинаре по численным методам механики вязкой жидкости (Новосибирск, 1984); на научных семинарах ЦАГИ, МФТИ, ИШЛех. АН СССР.
Основные уравнения. Схема расщепления .
В формулах (I.3-I.5) использованы следующие обозначения: t - шаг по времени, верхний индекс К (или К + І .) -номер временного слоя, {} - разностная аппроксимация соответствующего дифференциального оператора. Аппроксимация осуществляется на разнесенных сетках, когда давление определяется в центре ячейки, а компоненты скорости - в центрах соответствующих граней (рис. I).
В связи с переходом к моделированию пространственных течений в классическую разностную схему метода расщепления были внесены следующие изменения:
- Для аппроксимации конвективных членов использованы ориентированные, а не центральные разности.
- Постановка граничных условий реализована после первого, а не после третьего этапа.
- Разработана эффективная методика постановки граничных условий на внешних границах расчетной области.
- Для решения уравнения Пуассона использован попеременно-треугольный метод.
- Переработаны методики учета эффектов вязкости и постановки граничных условий на теле.
Введенные изменения и дополнения заимствованы из различных источников и принципиально новыми не являются. Однако, использование их в едином алгоритме осуществлено, по-видимому, впервые.
Краткая характеристика созданного комплекса программ
Краткая характеристика созданного комплекса программ Все проведенные расчеты были выполнены на импортной электронной вычислительной машине средней мощности, которая по своим параметрам примерно соответствует ЕС - 1060.
Комплекс счетных программ, реализованный на языке Фортран, содержал около 2 тыс. операторов, без учета блока решения уравнения Пуассона, который для ускорения счета был реализован на автокоде. Весьма существенным моментом является выбор формы представления результатов расчетов пространственных течений. Для оперативной обработки получаемых результатов был создан сервисный комплекс, включающий в себя блоки вывода информации на АЦПУ и на графопостроитель. Реализованный на языке Фортран, сервисный комплекс содержал около 3.5 тыс. операторов.
Предварительные замечания
Изучение нестационарных режимов обтекания различных тел потоком жидкости или газа является одной из наиболее актуальных проблем аэро- и гидродинамики. Внимание исследователей привлекает трудность и новизна возникающих задач, необходимость решения которых диктуется развитием современной техники - это движение самолетов и судов на подводных крыльях, поиски новых типов движителей, изучение и использование эффектов машущего полета птиц и насекомых и т.д.
Наряду с интенсивными теоретическими и экспериментальными исследованиями [51-57], все большее значение уделяется моделированию указанных задач на ЭВМ. Моделирование проводится с помощью различных методов и подходов: метода дискретных вихрей [2,58,59], прямого численного интегрирования уравнений Навье-Стокса, записанных относительно завихренности и функции тока [60,61,77,78]. Изучаются механизмы образования подъемной силы при взмахе крыльев, эффекты возникающие при обтекании вращающегося или колеблющегося глиптического цилиндра или профиля потоком несжимаемой вязкой жидкости. Однако, во всех рассмотренных работах по данной тематике вычислители ограничивались моделированием только плоских течений. В настоящей работе наряду с двумерными задачами, связанными с моделированием поворота тонкой пластины в потоке жидкости, рассматривается и аналогичная трехмерная. Численное интегрирование системы згравнений Навье-Стокса проводится по методу расщепления в системе отсчета, связанной с вращающимся телом. Учет эффектов, связанных с неинерциаль-ностью принятой системы отсчета, привел к необходимости перестройки описанного выше алгоритма и модификации методики постановки граничных условий на внешних границах расчетной области.