WWW.DISS.SELUK.RU

БЕСПЛАТНАЯ ЭЛЕКТРОННАЯ БИБЛИОТЕКА
(Авторефераты, диссертации, методички, учебные программы, монографии)

 

Pages:   || 2 | 3 |

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

-- [ Страница 1 ] --

ЦЕНТРАЛЬНЫЙ АЭРОГИДРОДИНАМИЧЕСКИЙ ИНСТИТУТ

имени профессора Н. Е. Жуковского

На правах рукописи

УДК 533.695, 629.7.015.3.036

Кажан Егор Вячеславович

Комбинированный метод численного решения стационарных

уравнений Рейнольдса и его применение к моделированию работы

воздухозаборника вспомогательной силовой установки

в компоновке с фюзеляжем летательного аппарата

Специальность 05.07.01 «Аэродинамика и процессы теплообмена летательных аппаратов»

Диссертация на соискание учной степени кандидата технических наук

Научный руководитель:

кандидат физико-математических наук, доцент Власенко Владимир Викторович

Научный консультант:

доктор технических наук, старший научный сотрудник Босняков Сергей Михайлович г. Жуковский 2013 г.

Оглавление Глава 1. Математическая постановка задачи

Система уравнений Рейнольдса, замкнутая моделью 1. турбулентности (q-)

Математические свойства системы уравнений Рейнольдса... 1. Постановка краевой задачи для уравнений Рейнольдса (на 1. примере задачи о моделировании работы вспомогательной силовой установки в компоновке с фюзеляжем)

Глава 2. Выбор и анализ численного метода для решения поставленной задачи

Анализ некоторых схем на основе модельного уравнения.... 2. 2.2. Базовый явный метод, основанный на схеме ГКР.................. 2.3. Линеаризованная неявная схема на основе схемы ГКР.......... 2.4. Особенности реализации неявной схемы

2.5. Особенности постановки численных граничных условий..... 2.6. Комбинированный метод с явной и неявной частями ............ Глава 3. Тестирование разработанного метода

3.1. Тест 1 – турбулентный пограничный слой на пластине........ 3.2. Тест 2 – профиль NACA0012.

3.3. Тест 3 - компоновка фюзеляж-крыло CRM (Common Research Model). 3.4. Тест 4 – Воздухозаборник и элемент мотогондолы тематической компоновки фюзеляж-крыло-пилон-мотогондола............ Выводы к Главе 3





Глава 4. Применение разработанной методики к моделированию течений в компоновке ВСУ с фюзеляжем ЛА

4.1. Расчтные исследования обтекания хвостовой части фюзеляжа и выбор положения створок ресивера ВСУ

4.2. Оценка влияния формы фюзеляжа

4.3. Анализ физических особенностей течения в воздухозаборном устройстве ВСУ

характеристики воздухозаборного устройства ВСУ

4.5. Полуэмпирическая методика оценки потерь на входе в двигатель ВСУ

Выводы к Главе 4

Выводы

Список использованных источников

Приложение 1. Анализ схем на основе модельного уравнения.................. Приложение 2. Базовый явный численный метод

Приложение 3. Матрицы неявного численного метода

Введение Предлагаемая работа является обобщением опыта проведения расчетных и экспериментальных исследований, накопленного в Центральном Аэрогидродинамическом институте (ЦАГИ) и в Центральном Институте опубликованы юбилейные сборники статей [1][2], в которых демонстрируется высокий уровень расчетных и экспериментальных работ, выполняемых в ЦАГИ. Эти книги являются основой, сформировавшей научные взгляды автора данной работы. В книге [3] представлен опыт разработки пассажирских магистральных самолетов, который был учтен в практической части настоящей работы. Большой задел, который был создан в ЦИАМ в результате многочисленных теоретических и практических исследований, описан в книге [4]. Численный метод, который был создан C.К.Годуновым [5] и в настоящее время имеет всемирное признание, успешно и плодотворно развивался и совершенствовался в обоих институтах. Автор внимательно изучил книгу [6], в которых подробно описаны все особенности реализации метода Годунова. В ЦАГИ впервые был предложен метод Годунова 2-го порядка аппроксимации по пространству [7], который впоследствии был обобщен на многомерные задачи в ЦИАМ [8]. Явная схема Годунова 2-го порядка аппроксимации по времени [9] зарекомендовала себя как весьма надежный инструмент для решения сложных практических задач и в течение нескольких десятилетий успешно эксплуатируется в ЦАГИ [10]. В последнее время в ЦАГИ исследуется возможность применения идей Годунова в конечно-элементных методах высокого порядка аппрксимации [11]. В ЦИАМ был разработан неявный аналог схемы Годунова [12] применительно к уравнениям Эйлера, который был впоследствии обобщен предложена еще одна, весьма эффективная релаксационная неявная схема для решения уравнений Навье-Стокса [14]. Последняя схема была использована при создании различных программ ЦИАМ, включая модификациями) была использована и в настоящей работе для ускорения численной технологии [10], которая использовалась в ЦАГИ к началу данной работы.

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

запуск маршевого двигателя с помощью воздушного стартера на питание сжатым воздухом системы кондиционирования салона и кабины экипажа на земле и в полете;





питание в полете сжатым воздухом системы противообледенения двигателей);

питание бортовой сети самолета электроэнергией переменного тока Актуальность темы. Новый пассажирский самолет МС-21 использует многофункциональную вспомогательную силовую установку (ВСУ), работающую на всех режимах полета. Увеличение энергетической нагрузки на ВСУ и усложнение условий ее функционирования повышает требования по эффективности и устойчивости е работы. Как правило, воздухозаборник установки расположен в хвостовой части фюзеляжа в развитом турбулентном пограничном слое. Поэтому при расчетных исследованиях работы воздухозаборника ВСУ приходится моделировать с высокой точностью обтекание всего планера, что требует использования плотных сеток и значительных ресурсов ЭВМ, включая время расчета.

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

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

В ЦАГИ была разработана и получила широкое применение технология вычислительного эксперимента под названием "Электронной прикладных программ (ППП) EWT – ЦАГИ [17] [18] включает весь набор инструментов для реализации вычислительного эксперимента.

Вычислительный эксперимент позволяет моделировать турбулентные течения вязкого газа при обтекании реальных геометрий. Среди прочих, аэродинамической трубы. Перед автором данной работы стояла задача ускорить получение стационарных решений в рамках описанной технологии и добиться нового качества путем решения класса задач, недоступного ранее. Программный продукт EWT – ЦАГИ включает в себя несколько прикладных программ, в том числе V3Solver [19] и ZEUS [20].

Для получения стационарных решений в этих программах используется метод установления решения по времени. Стационарное течение получается как предел некоторого нестационарного процесса. Возможны разные подходы к моделированию этого процесса. Известно, что численные схемы, которые могут быть использованы в методе установления, можно разделить на два класса - явные и неявные. Они имеют свои достоинства и недостатки. Неявные схемы могут быть абсолютно устойчивы, но требуют больше вычислительных затрат на совершение итерации по сравнению с явными аналогами. И явные и неявные численные схемы широко используются в вычислительной аэродинамике и имеют множество реализаций [21]. В Электронной Аэродинамической Трубе до последнего времени использовались только явные численные схемы. Для ускорения получения стационарных решений применялся метод локального шага по времени. Это означает, что в разных ячейках сетки расчет проводится с разными значениями шага по времени, определяемыми локальными условиями устойчивости.

Этого было недостаточно для эффективного использования в практических расчетах, и перед автором настоящей работы была поставлена задача: в рамках технологии «Электронной аэродинамической трубы» (с использованием уже разработанных функций, алгоритмов и настроек) ускорить получение стационарных решений.

Один возможный подход к решению задачи – это многосеточные методы (multigrid [22] [23]), другой подход – использование неявной схемы. В работе автора был выбран второй подход.

вычислительной аэродинамики, как правило, основаны на принципе линеаризации зависимости от времени уравнений движения газа [21].

Значения параметров на неизвестном временном слое представляются в виде F (u ) F (u ) (u )u, где u u n1 u n, (u ) - некоторая аппроксимация матрицы Якоби, вычисляемая по параметрам на известном временном слое u n. В результате, аппроксимация уравнений движения сводится к системе линейных алгебраических уравнений (СЛАУ) относительно приращения параметров u (неявная схема в “дельта-форме” [24]). При решении трхмерных уравнений Навье-Стокса и Рейнольдса матрица такой СЛАУ имеет блочно ленточную структуру с несколькими ненулевыми блочными диагоналями, разделенными многочисленными нулевыми диагоналями.

Существующие методы решения таких СЛАУ можно подразделить на два больших класса.

В первом классе методов СЛАУ для u не упрощается. Несколько “шагов по времени” с точным решением СЛАУ при бесконечно больших шагах по времени эквивалентны решению нелинейной стационарной системы уравнений методом Ньютона, что обеспечивает квадратичную скорость сходимости при наличии хорошего начального приближения [25]. Однако, как правило, начальное поле сильно отличается от решения, поэтому итерационный процесс приходится начинать с небольших шагов по времени. При точном решении СЛАУ требуемые вычислительные ресурсы неприемлемо велики [21]. Весьма употребительны различные версии метода сопряженных градиентов [26] [27], среди которых особенно эффективным является обобщенный метод минимальных невязок (GMRES [28]) в безматричной реализации [29]. Однако в общем случае для таких методов решающее значение имеет применение различных предобуславливателей [30]. При решении “плохих” практических задач эти методы не обеспечивают достаточную робастность [30].

Более простым в реализации является второй класс методов, основанный на упрощении СЛАУ. В целом, методы этого класса можно подразделить на методы факторизации, связанные с представлением матрицы системы в виде произведения нескольких матриц специального вида, и методы, связанные с представлением матрицы системы в виде суммы нескольких матриц. Высокоэффективная схема первой группы была предложена в работе [31]. Это так называемый метод переменных направлений (ADI), при котором матрица системы факторизуется на три множителя, связанные с численным дифференцированием в каждом из пространственных направлений. Каждый множитель приводится к системе, решаемой матричной прогонкой или даже рядом скалярных прогонок. Однако метод ADI в трехмерном случае обладает условной устойчивостью. Существует ряд других способов факторизации, обладающих абсолютной устойчивостью, например, [32]. Во второй группе схем, основанных на упрощении матрицы системы, используются такие эффективные методы решения СЛАУ, как блочные методы Якоби и Гаусса-Зейделя [33] и методы релаксации, среди которых особенно популярен метод LU-SGS и его модификация LU-SSOR [34]. Как правило, число итераций таких методов на каждом шаге по времени ограничивается [21], так что решение упрощенной СЛАУ не находится, что не мешает сходимости процесса в целом к стационарному решению. В отличие от методов первой группы, методы этого типа легко обобщаются на неструктурированные сетки [35]. Смена направления обхода ячеек структурированной сетки, которая ускоряет сходимость метода ГауссаЗейделя, на неструктурированных сетках может быть успешно заменена произвольной перенумерацией ячеек [36].

При упрощении СЛАУ неявной схемы существенную роль играет сокращение шаблона схемы на неявном слое. Наиболее употребительным способом для достижения этой цели является метод отсроченной коррекции (deferred correction) [37], при котором неявный оператор, применяемый к u, аппроксимируется с первым порядком точности по пространству на компактном шаблоне, и лишь в явном операторе применяется аппроксимация высокого порядка на разврнутом шаблоне.

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

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

Неявный метод [38], предлагаемый в настоящей работе, основан на методе отсроченной коррекции и использует для решения системы линейных уравнений блочный метод Гаусса-Зейделя с перенумерацией ячеек. Этот метод продолжает традиции отечественной школы вычислительной аэродинамики, основанной С.К.Годуновым и его последователями.

произвольного разрыва. Однако неявная схема, предложенная в классической книге С.К.Годунова с соавторами [6], не относится к этому классу. Не относится к этому классу и базовый вариант метода ADI [31].

Однако уже в одной из последующих работ авторов ADI [32] был предложен метод ADI с учетом направления распространения информации, основанный на расщеплении вектора потоков. Родственными к работе [32] являются классические отечественные работы [12], [13]. В этих работах также используется подход ADI, однако матрицы Якоби конвективных потоков аппроксимируются на гранях ячеек, а в явном операторе используется нелинейная схема Годунова второго порядка аппроксимации. Линеаризация Роу признана оптимальным компромиссом между точностью и эффективностью при вычислении параметров на гранях ячеек [21]. В работах [14], [40], [41] описана неявная схема, основанная на решении задачи Римана, которая не использует факторизацию неявного оператора, с решением СЛАУ блочным методом Гаусса-Зейделя. Неоценимые консультации по своему опыту применения подобных схем дали автору разработчики программы КоБра [15].

При построении неявных методов для решения уравнений Рейнольдса, замкнутых дифференциальной моделью турбулентности, возникает проблема выбора аппроксимации источников, связанных с производством и диссипацией турбулентности. Нередко для источников используется неявная аппроксимация; но следует учитывать возможные случаи неустойчивости, описанные в [21]. В работе [42] предложена абсолютно устойчивая схема, в которой для производства турбулентности используется явная аппроксимация, а для диссипации – неявная. В настоящей работе используется другой метод, основанный на анализе собственных чисел матрицы Якоби источников. Этот метод был сформулирован автором самостоятельно как логическое развитие идей, изложенных в [10]. Впоследствии автору стало известно, что аналогичный подход был предложен в работе [43], где сообщается, что он обеспечивает более быструю сходимость, чем метод [42].

Главной оригинальной особенностью численного метода [38], который будет описан ниже, является локальный выбор явной или неявной схемы и способа осуществления шага по времени (глобальный, локальный) в зависимости от соотношения между заданным глобальным шагом по времени и локальным условием устойчивости явной схемы. Далее будет показано, что данный комбинированный подход позволяет существенно ускорить получение стационарного решения по сравнению с методами, использующими одну и ту же схему во всей расчетной области. Следует отметить, что существует обширный класс методов зональной декомпозиции [44], где в разных блоках расчетной сетки используются разные численные методы. В отличие от методов этого класса, в предлагаемом подходе выбор схемы определяется в каждой ячейке сетки на основе локальных параметров численного решения; при этом удается сохранить однородность алгоритма. Локальный выбор схемы обеспечивает большую гибкость метода. Эффективность данного подхода будет продемонстрирована на ряде тестовых задач. Близкий подход, но содержащий ряд существенных отличий, был использован в работах [45-48].

Помимо настоящего введения, диссертация включает 4 главы, заключение, 3 приложения и список использованных источников. Содержание работы изложено на 154 страницах основного текста и 48 страницах приложения.

Список использованных источников содержит 117 наименований. В работе содержится 49 иллюстраций в основном тексте и 13 в приложениях.

В первой главе диссертации описана математическая постановка задачи.

дифференциальные модели турбулентности. Перечислены рассматриваемые граничные условия.

Вторая глава диссертации посвящена выбору и анализу численного метода для решения поставленной задачи. Рассматривается явная численная методология, реализованная к началу данной работы в пакете прикладных программ EWT-ЦАГИ. Для ускорения этой методологии было принято решение разработать ее неявный аналог. Описываются процедура и критерии выбора неявной схемы, включая результаты предварительных тестов. Проводится анализ свойств (аппроксимация, устойчивость и др.) выбранной неявной схемы для модельного уравнения.

Дано детальное описание схемы с линейным неявным сглаживателем для уравнений Рейнольдса. Описываются численные граничные условия.

Вводится комбинированный метод с явной и неявной частью и выбором заданной величины глобального шага по времени.

Третья глава посвящена тестированию разработанного метода. Точность базового явного численного метода, реализованного в пакете прикладных программ EWT-ЦАГИ, была проанализирована в работе [49]. Первая цель тестирования – показать, что стационарное решение, полученное с использованием предлагаемого комбинированного метода, совпадает со стационарным решением, полученным по базовому явному методу.

Вторая и основная цель тестирования – демонстрация ускорения расчета глобально-неявным методом. Для тестирования рассматриваются классические тестовые задачи (пограничный слой на пластине, профиль крыла, модель CRM “фюзеляж-крыло”), а также модельный тестовый пример, основанный на расчете реального ЛА на практической сетке.

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

Работа проводилась по заказу конструкторского бюро и описана в отчетах [50-55]. Эксперимент проводился в аэродинамической трубе СВС-2 ЦАГИ. Дано описание экспериментальной установки и средств измерений. Затем изложено описание расчтных исследований обтекания хвостовой части фюзеляжа самолта МС-21. По результатам этих расчетов проводится выбор пограничного слоя на входе в ресивер. На основании расчетов полного обтекания ЛА оценивается влияние формы фюзеляжа на условия течения в окрестности входа в ресивер ВСУ. Далее проводится анализ физических особенностей течения в воздухозаборном устройстве ВСУ. Делаются рекомендации по выбору угла установки створки ресивера. Для оценки точности численного моделирования течения в ресивере результаты численного моделирования эксперимента в аэродинамической трубе обсуждаются результаты трехмерного расчета реальной компоновки ВСУ с хвостовой частью фюзеляжа и проводится сопоставление этих данных с данными численного моделирования экспериментальной модели.

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

На защиту выносятся:

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

моделировании работы воздухозаборника вспомогательной силовой установки (ВСУ). Анализ физических особенностей течения в ресивере ВСУ для разных режимов ее работы.

Личный вклад автора:

1. Модификация вычислительной методологии (программы) путем применения комбинированной вычислительной схемы.

2. Демонстрация эффективности (ускорения) разработанной программы путем проведения верификационных тестовых расчетов.

3. Подготовка и проведение расчетных исследований ВСУ как в условиях аэродинамической трубы СВС-2, так и в компоновке с фюзеляжем самолета МС-21.

4. Участие в эксперименте, подготовленном и осуществленном в АДТ СВС-2 ЦАГИ, с обработкой, анализом и сопоставлением расчетных и экспериментальных данных.

Научная новизна работы состоит в следующем:

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

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

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

Достоверность полученных результатов обосновывается:

Сопоставлением результатов расчетов с классическими решениями (раздел 3.1).

Сопоставлением результатов расчетов с экспериментальными данными (разделы 3.2 и 4.4).

Сопоставлением результатов расчетов с расчетными данными других авторов (раздел 3.2).

Сопоставлением результатов расчетов, полученных с использованием различных численных подходов (разделы 3.1, 3.2, 3.3, 3.4) и разных физических моделей (раздел 3.2).

Практическая ценность работы состоит в следующем. Полученные в результате работы технические решения приняты к исполнению при прохождении макетной комиссии Инженерного центра им Яковлева ОАО «Иркут». Разработанная численная методология внедрена в практику расчетных работ в ЦАГИ. Разработанная автором программа COMGLEI включена в пакете прикладных программ EWT-ЦАГИ. Получено свидетельство о государственной регистрации №2013610173 [56].

Апробация работы:

Материалы работы докладывались и обсуждались на 14 отраслевых и международных конференциях [57-75].

Публикации. Результаты диссертации изложены в 3 печатных работах [76] [38] [77].

В заключении автор хотел бы выразить глубокую благодарность к.т.н.

С.В. Михайлову и к.т.н. С.В. Матяшу за неоценимую помощь в разработке программы «COMGLEI», а также с.н.с. В.Ф. Третьякову и в.и.к.

ОАО «Корпорация Иркут» Е.П. Быкову за помощь в проведении исследований вспомогательной силовой установки самолта МС-21.

Глава 1. Математическая постановка задачи 1.1 Система уравнений Рейнольдса, замкнутая моделью Прогресс в развитии авиации в последние годы в значительной степени определяется оптимизацией аэродинамических обводов и применением новых технических решений с использованием существенно нелинейных режимов обтекания, требующих обязательного учта вязкости и состояния пограничного слоя. К таким задачам можно отнести использование больших углов атаки с управляемым отрывом потока, разработку высоконесущих скоростных профилей с развитыми сверхзвуковыми зонами с замыкающим скачком уплотнения, исследование технических шумоглушения в соплах ТРДД, охлаждения элементов конструкции, струйного управления вектором тяги, применение интенсивных щелевых струй для управления обтеканием внешним потоком) и т.д.

Численное решение системы уравнений Рейнольдса (осредннных по времени уравнений Навье-Стокса) [10][21][78] позволяет детально исследовать физическую картину течения, получить количественную оценку параметров течения и аэродинамических характеристик, исследовать влияние различных факторов и параметров. На основе этого можно сформировать рекомендации по повышению эффективности исследуемых технических решений и сэкономить время и объмы необходимых экспериментов, а также облегчить понимание и трактовку их результатов.

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

Считается [18], что в настоящее время научная разработка методов расчта с использованием осредннной по Рейнольдсу (RANS) системы уравнений Навье-Стокса (с замыканием в виде дифференциальной модели турбулентности) в основном закончена. Это не означает, что работы в данном направлении полностью остановлены. Тем не менее, основные усилия научных школ направлены на развитие методов LES, с которыми связывается будущее вычислительных работ Действительно, опубликованные дифференциальные модели турбулентности [42], как правило, настроены на решение узкого класса задач и их использование правомерно только в рамках зонального подхода. Например, известная модель (k-) [80] дат приемлемые результаты в пристеночных областях, где имеется развитый турбулентный пограничный слой. Другая использования в слоях смешения, и е применение за пределами этих слов проблематично. Удачным решением является комбинированный подход, использованный в модели SST [84], в котором происходит плавный переход от модели (k-) к модели (k-) в зависимости от функции удалнности тврдых стенок. Такой подход может быть применен и к другим моделям турбулентности – например, к модели (q-) [85].

Рассмотрим произвольное течение газа. Вырежем в этом течении фиксированный объем V : поверхность S этого объема неподвижна во времени, и сквозь нее течет газ. Запишем для этого объема законы сохранения массы, импульса и энергии.

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

В левой части этого уравнения стоит скорость изменения полного количества величины a в объеме V. Правая часть описывает причины, вызывающие это изменение:

потоки величины a сквозь поверхность S ( F a - поток величины a сквозь единицу площади в единицу времени, ds - вектор элемента площади, перпендикулярный поверхности);

локальные источники и стоки величины a ( W a - это источниковый член, который описывает скорость возникновения или расходования величины a за счет локальных источников и стоков).

Применив к интегралу по поверхности теорему Гаусса-Остроградского, внеся производную по времени под знак интеграла (поскольку V const ) и устремив объем к нулю, получим дифференциальное уравнение, выражающее этот закон сохранения:

Здесь Fxa, Fya, Fza - потоки величины a сквозь единичную площадку, ориентированную перпендикулярно осям x, y, z. Дифференциальное интегральный закон (1), так как (2) неприменимо к разрывным течениям.

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

Система уравнений Рейнольдса – это система уравнений для средних параметров газа – температуры T, трех компонент скорости u, v, w и давления p ; она получается осреднением уравнений Навье-Стокса по времени. В результате возникают дополнительные члены, связанные с пульсациями. Для замыкания системы уравнений среднего течения призваны так или иначе связать турбулентные потоки с параметрами рассматривается только одна модель турбулентности - q модель Коукли [85], [87], [88]. Эта модель включает два дополнительных турбулентности:

кинетическая энергия турбулентности. Параметр q измеряется в м/с и кинетической энергии турбулентности ( - плотность газа, ij - тензор вязких напряжений). Этот параметр измеряется в Гц и пропорционален характерной частоте турбулентных пульсаций.

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

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

Система уравнений Рейнольдса, замкнутая моделью турбулентности ( q ), представляет собой систему законов сохранения (2), записанных для величин консервативных параметров U, u, v, w, E, q,.

Объединяя уравнения (2) для каждого компонента вектора U, получим систему уравнений следующего вида:

Потоки вектора u в системе Рейнольдса целесообразно разделить на конвективные и диффузионные потоки: Fi Fi конв Fi дифф, где Здесь были введены следующие обозначения:

I ij - суммарный диффузионный поток i -й компоненты импульса в направлении оси x j, т.е. сумма молекулярных и турбулентных потоков импульса – вязких и турбулентных напряжений:

I jiui I jxu I jyv I jz w - работа вязких и турбулентных напряжений при конвективном переносе газа;

j - суммарный диффузионный поток тепловой энергии вдоль x j :

где C p - удельная теплоемкость газа при постоянном давлении, Pr 0.72, PrT 0.9 ;

T jk - суммарный диффузионный поток кинетической энергии турбулентности вдоль x j ; модель для этого потока описана в [10] [91];

T jq - суммарный диффузионный поток параметра q вдоль x j :

T j - суммарный диффузионный поток параметра вдоль x j.

вычисляется по формуле Сатерленда [92]. Коэффициент турбулентной вязкости в модели ( q ) вычисляется по формуле ( dW – расстояние от данной точки до ближайшей твердой поверхности), сжимаемости, которая основана на результатах работы [93] ( M T 2q / c – турбулентное число Маха, c – скорость звука).

Вектор источниковых членов W 0, 0, 0, 0, 0, S q, S в правой части системы (3) содержит ненулевые элементы только в уравнениях модели турбулентности ( q ). Источниковые члены S q и S турбулентности; их удобно представить в следующем виде:

В отличие от оригинальной версии модели турбулентности ( q ), в данной работе используется коррекция к производству турбулентности, предложенная в [94]. Эта коррекция позволяет существенно снизить нетонкослойных течениях (в частности, на скачках уплотнения и в окрестности точек растекания при натекании потока на твердую поверхность). Для этой цели производство кинетической энергии турбулентности вычисляется по формуле Pk T | G |2, где G – “тонкослойная” часть градиента проецированием градиента модуля скорости g |V | на направление e, перпендикулярное преимущественному направлению векторов приращений скорости в окрестности данной точки: G g e e, g.

Автором модели ( q ), Т.Коукли, разработаны две версии этой модели турбулентности, одна из которых, опубликованная в работе [87], лучше описывает свободную турбулентность, а другая, опубликованная в работе [88], лучше описывает турбулентные пограничные слои. Эти версии отличаются наборами коэффициентов.

В работе [85] по аналогии с известной моделью SST Ментера [84] было предложено использовать переходную функцию, которая близка к единице на твердой поверхности и обращается в нуль в свободной турбулентности. В качестве такой функции был использован упрощенный формуле C (1 F )С1 FС2, где C1 – любой коэффициент из 1-го варианта модели ( q ), а C2 – аналогичный коэффициент из 2-го варианта этой модели турбулентности. Именно этот вариант модели ( q ) используется в настоящей работе.

1.2 Математические свойства системы уравнений Рейнольдса турбулентности ( q ), можно выделить три группы членов, которые очень сильно отличаются по математическим свойствам: конвективные потоки, диффузионные потоки и источниковые члены.

Для анализа математических свойств конвективных потоков обычно рассматривают систему уравнений вида (3), из которой выброшены диффузионные и источниковые члены. В случае одномерного течения без разрывов, когда все параметры изменяются только вдоль оси x, а вдоль осей y и z, перпендикулярных оси x, течение равномерно:

система уравнений запишется в виде:

Далее проводится характеристический анализ матрицы Якоби системы (12). Известно, (см. например [78]), что эта матрица имеет следующие собственные числа 1,2,3,4,5 u (пятикратное), 6 u c, 7 u c. Поскольку все собственные числа этой матрицы действительны, а все собственные векторы линейно независимы, система уравнений (12) является гиперболической [10], т.е. описывает направленный перенос возмущений вдоль характеристик. Система (12) имеет три семейства (траектории малых звуковых возмущений). Вдоль характеристик первого инвариант), z 2 v, z3 w, z4 q и z5 ; вдоль характеристик второго и третьего семейств – инварианты z6,7 u соответственно [10] [78].

Если течение изэнтропично ( Диффузионные члены присутствуют во всех уравнениях системы Рейнольдса, кроме уравнения неразрывности. Старшие производные в этих уравнениях (производные второго порядка) связаны именно с присутствием диффузионных членов.

Грубую оценку математических свойств уравнений Рейнольдса можно провести, если рассматривать каждое из уравнений системы как уравнение для какого-нибудь одного газодинамического параметра, а все остальные параметры газа, входящие в это уравнение, - как заданные переменные коэффициенты. Такой подход используется, например, в [10].

Например, если рассматривать уравнение неразрывности, как уравнение для, то это уравнение в одномерном нестационарном случае можно представить в виде стоит гиперболический оператор, то это уравнение, как уравнение для, гиперболично. Уравнение для x -компоненты импульса в одномерном как уравнение для u, то оно параболично.

Результаты подобного грубого анализа уравнений Рейнольдса, проведенного в [10], можно сформулировать следущим образом.

В стационарном случае система уравнений Рейнольдса может быть отнесена к смешанному гиперболически-эллиптическому типу: уравнение неразрывности гиперболично, остальные – эллиптичны (если рассматривать их по отдельности, как уравнения для u, v, w, E, q, ).

распространяется направленно (вдоль линий тока), а информация об остальных шести перечисленных параметрах передается во все стороны.

В нестационарном случае система уравнений Рейнольдса может быть отнесена к смешанному гиперболически-парабоическому типу: уравнение неразрывности гиперболично, остальные – параболичны. Информация об одном параметре (возмущения ) распространяется направленно (вдоль траекторий объемов газа). Информация об остальных шести параметрах передается по пространству - во все стороны, а по времени – направленно (из прошлого в будущее).

Для анализа математических свойств системы уравнений Рейнольдса, связанных с присутствием источниковых членов модели турбулентности, в работе [10] рассматривается система уравнений, в которой оставлены только источниковые члены и лагранжевы производные вдоль траекторий объемов газа Линеаризуем источниковые члены в окрестности некоторой точки пространства-времени ( r0, t0 ) :

где R - матрица Якоби источниковых членов. Тогда решение системы (13) имеет вид [10]:

где ek - собственный вектор матрицы R0, соответствующий собственному числу k этой матрицы, а z0 P 1 R0 1 W0, где P - матрица, столбцами которой являются собственные векторы ek. В формуле (14) используется лагранжево представление: под значением параметра в момент времени t понимается значение параметра в этот момент времени в том объеме газа, который в момент времени t0 находился в точке r0.

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

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

Стационарное решение получается как предел нестационарной адаптации течения к заданным граничным условиям. При таком подходе расчет ведется на базе нестационарных уравнений Рейнольдса маршем по времени (т.е. путем последовательного перехода от одного слоя t const к следующему). Считается, что стационарное решение достигнуто, если в течение контрольного периода все интегральные характеристики обтекаемых тел (включая расход в контрольном сечении тракта силовой установки) не выходят за границы заданного узкого “коридора”.

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

эта информация влияет на течение внутри.

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

Пусть Vn - проекция скорости на внешнюю нормаль к границе расчетной области, а M n Vn / c - число Маха по нормали к границе.

Рассмотрим самые распространнные варианты границ:

инварианты передаются изнутри расчтной области. На границе не нужно задавать ни одного параметра.

2. Граница дозвукового вытекания потока ( 0 M n 1 ). Внутрь по который вместе с инвариантами, приходящими из расчтной области, позволяет однозначно определить входящий от границы инвариант. Часто задают давление.

3. Граница дозвукового втекания потока ( 1 M n 0 ). Внутрь по скорости, касательные к границе расчетной области), z4 q, параметров газа. Например, на входе в сопло часто задают турбулентности и накладывают условие V 1 V 2 0.

4. Граница сверхзвукового втекания потока ( M n 1 ). Внутрь по характеристикам передаются все инварианты. На границе нужно задавать значения всех параметров потока.

5. Тврдая поверхность с непротеканием газа (Vn 0 ). Внутрь по z5 Vn. На границе нужно задавать один параметр газа. Для тврдой границы это условие Vn 0.

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

теоретически требуется задать на границах различных типов в случае системы уравнений Эйлера и в случае системы уравнений Навье-Стокса.

Таблица 1. Количество параметров, которые нужно задать на границе Рейнольдса Однако в работе [10] показано, что из-за постепенного затухания возмущений при диффузионном переносе область влияния границы, обусловленного диффузией, является конечной. В [10] дана оценка размера этой области - Lдифф. Если внешняя граница удалена от интересующей области течения на расстояние большее, чем Lдифф, то на такой границе краевые условия можно ставить так же, как и для невязкого (теплоизолированная твердая стенка).

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

Например, на теплоизолированной стенке необходимо, чтобы все три компоненты скорости были равны нулю (прилипание) и чтобы нулю был равен тепловой поток через стенку ( T / n 0 ). Из-за прилипания должны обнулиться и пульсации скорости на стенке, поэтому там следует задать q 0. В работе [95] для модели ( q ) рекомендовано задавать на стенке Рассмотрим в качестве примера постановку краевой задачи для моделирования работы вспомогательной силовой установки (ВСУ), которая исследована в Главе 4.

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

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

Рисунок 1. Схема компоновки ВСУ в корпусе ЛА Для моделирования работы ВСУ необходимо воспроизвести обтекание всего ЛА в целом, чтобы получить правильную структуру пограничного слоя перед входом в ресивер ВСУ. Соответственно, расчетная область должна включать поверхность всего ЛА (рис.2). На всей поверхности ЛА ставится описанное выше граничное условие теплоизолированной стенки с прилипанием потока.

Рисунок 2. Общий вид расчетной области для моделирования ВСУ При моделировании обтекания ЛА на крейсерском режиме полета ( M 0.8) внешняя граница расчетной области удаляется от поверхности ЛА на расстояние нескольких длин фюзеляжа. Это расстояние намного превосходит характерный размер области влияния границы за счет диффузии Lдифф, поэтому на внешней границе ставятся “невязкие” граничные условия, основанные на анализе числа Маха по нормали к границе (см. выше). Инварианты Римана, входящие через границу внутрь расчетной области, вычисляются по параметрам невозмущенного потока на бесконечности.

В настоящей работе участок тракта ВСУ, содержащий компрессор, камеру сгорания и турбину, исключается из расчетной области, а в выходном сечении воздухозаборника ВСУ и во входном сечении выпускного устройства ВСУ ставятся специальные граничные условия – “активные диски” [96] – см. рис.1.

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

Выходное сечение воздухозаборника является границей дозвукового вытекания потока из расчетной области. Поэтому на этой границе должен быть задан только один параметр потока. В качестве такого параметра задается расход газа через выходное сечение (в предположении, что параметры газа во входном сечении равномерны). Этот расход принимается равным расходу через выпускное устройство ВСУ - Gexit.

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

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

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

Вначале в Приложении 1 проведено построение неявной численной схемы для решения стационарных задач методом установления. Схема строится методом конечного объма для модельного уравнения где источниковый член S (u) и вязкость (u) - произвольные нелинейные функции.

Уравнение (15) допускает непрерывное нетривиальное стационарное решение в невязком случае (при v 0, S (u) 0, (u) 0 ) – веер разрежения, ограниченный характеристиками и.

допускает нетривиальное стационарное решение типа пограничного слоя.

установления строится на основе следующих принципов:

аппроксимация должна быть абсолютно устойчивой;

2) поскольку в стационарном пределе зависимость от времени исчезает, достаточно ограничиться первым порядком аппроксимации по времени;

3) в стационарном пределе, когда совпасть с базовой явной схемой.

В базовой явной схеме [10] все потоки вычисляются на явном слое, а для источниковых членов записывается полунеявная (второго порядка) аппроксимация. Аппроксимации конвективных потоков в явной схеме основаны на точном решении задачи Римана о распаде разрыва.

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

Исходя из принципа 1, мы вычисляем все потоки на неявном слое, т.к.

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

Введен оператор приращения параметров за шаг по времени:

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

Неявная часть все равно должна в конце концов обратиться в нуль.

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

Производится линеаризация приращений потоков, основанная на принципе Роу. Матрица Роу фиксируется на явном слое. Аналогично, значение вязкости в диффузионном члене также фиксируется на явном слое.

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

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

Также в приложении 1 приведн анализ явной и неявной схем с использованием двух стандартных методов вычислительной аэродинамики – метода дифференциального приближения и метода Неймана.

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

При помощи этого метода в Приложении 1 даны ответы на вопросы:

почему явная схема теряет устойчивость при числе Куранта CFL 1 и почему построенная неявная схема устойчива при описании нестационарного развития течения с большими числами Куранта CFL 1.

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

Метод Неймана [98] (известный также как метод рядов Фурье) является наиболее универсальным способом точного определения устойчивости схемы. Условие устойчивости, которое получается методом Неймана, является не только необходимым, но и достаточным [98]. Кроме того, метод Неймана позволяет количественно исследовать амплитудные и фазовые ошибки схем. Но, к сожалению, его удается применить только к линейным задачам.

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

взаимодействии пограничного слоя и движущегося скачка уплотнения.

Используется модельное уравнение (15) при v 0, S (u) 0, (u) const.

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

Скачок уплотнения был задан таким образом, чтобы одна из границ типа “joint” рассекала его пополам. Вторую границу он пересекал при движении всем фронтом.

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

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

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

В следующем разделе построенная неявная схема будет применена к полной системе уравнений Рейнольдса.

Базовый явный метод, основанный на схеме ГКР К началу данной работы в пакете прикладных программ EWT-ЦАГИ был реализован явный численный метод для решения системы уравнений Рейнольдса, замкнутой моделью турбулентности q- [85], [87], [88]– см. (3)-(5). Подробное описание этого метода дано в Приложении 2; оно основано на работах [10] и [16] [91]. В настоящем разделе будут лишь кратко перечислены ключевые характеристики этой базовой явной схемы.

Численная схема строится в рамках конечно-объемного подхода.

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

Здесь n t n 1 t n - шаг по времени (индексом n нумеруются времени потоков вектора u сквозь грань ячейки Ei ( ds - элемент площади, n (nx ; n y ; nz ) - единичный вектор внешней нормали к грани).

W n WdVdt - аппроксимация осредненных за шаг по времени и по объему ячейки источниковых членов.

Для описания конвективных потоков используется явная монотонная схема Годунова-Колгана-Родионова (ГКР) [5][6][7][9], для описания диффузионных членов – явная центрально-разностная аппроксимация с модификацией для работы на неравномерных сетках, а для источниковых членов – локально-неявная аппроксимация. Схема номинально имеет второй порядок аппроксимации по пространству и времени.

Для совершения шага по времени в базовом явном методе используется двухшаговая процедура предиктор-корректор [9]. На шагепредикторе потоки F аппроксимируются по параметрам на явном (известном) временном слое t n, а на шаге-корректоре – по параметрам на слое t n 1/ 2 t n n / 2. (Эти параметры определяются на шаге-предикторе.) Источниковые члены на корректоре вычисляются на слое t n 1 / 2 ; на предикторе они не вычисляются, а используется их значение на слое t n 1 / 2, которое было найдено при совершении прелдыдущего шага по времени.

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

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

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

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

Линеаризованная неявная схема на основе схемы ГКР В качестве основы для ускорения базового явного численного метода нужна неявная схема, которую можно было бы рассматривать как дополнение к базовому явному методу и которая позволила бы сохранить точность стационарных решений, обеспечиваемую базовым явным методом. Желательно также, чтобы эта схема обеспечивала не слишком сильный рост затрат на одну итерацию по сравнению с базовым явным методом. После анализа литературы автор остановился на неявной реалаксационной схеме для уравнений Навье-Стокса, предложенной в используется для получения стационарных решений системы уравнений Рейнольдса в различных программах ЦИАМ, в т.ч. – в известной программе КоБра [15].

Поскольку интерес представляет только стационарное решение уравнений Рейнольдса, данная неявная схема имеет первый порядок по аппроксимируются по параметрам на неявном (неизвестном) временном слое t n 1 :

Линеаризация потоков позволяет представить их в виде:

Здесь и далее оператор обозначает приращение параметров за один шаг по времени: a a n 1 a n. uside - приращение вектора консервативных параметров на грани ячейки, ui - в текущей ячейке, а uNBR - в соседней ячейке, расположенной по другую сторону от данной грани текущей ячейки. Ai и ANBR - некоторые матрицы, вычисляемые по параметрам на явном слое. Проведя также линеаризацию источниковых членов в текущей ячейке i, представим неявную схему (19) в виде:

турбулентности ( q ), матрицы Ri и RNBR имеют размерность 77. Из-за присутствия uNBR уравнение (20) для данной ячейки связано с аналогичными уравнениями для соседних ячеек. Таким образом, (20) представляет собой набор из 7 строчек системы линейных алгебраических уравнений (СЛАУ), имеющей размерность 7 N ( N - число ячеек в текущем блоке расчетной сетки). Многоблочный случай рассматривается далее при обсуждении граничных условий (раздел 2.5). Каждой ячейке текущего блока расчетной области соответствует набор из 7 строчек вида (7). В случае трхмерных уравнений Рейнольдса, замкнутых моделью турбулентности ( q ), получающаяся система линейных алгебраических уравнений (20) имеет ленточную матрицу, с количеством ненулевых диагоналей, соответствующим шаблону схемы на неявном слое. Если потоки на неявном слое удается представить в виде (19), то шаблон схемы содержит 7 ячеек на неявном слое. Таким образом, матрица СЛАУ содержит 7 ненулевых диагоналей, элементами которых являются заполненные блоки размерности 77. Ненулевые диагонали располагаются в матрице не подряд, а с разрывами, определяемыми сдвигом индекса i между текущей ячейкой и ее соседями.

Для сравнения рассмотрим аналогичную запись базовой явной схемы:

В неявной схеме (20) можно выделить слагаемые, отличающие е от явной схемы (21):

Здесь I - единичная матрица.

Дополнительные члены, содержащие ui и uNBR, обеспечивают устойчивость схемы, сглаживая возмущения, порождаемые неустойчивым при больших шагах по времени явным оператором (см. Приложение 1).

сглаживателем.

При сходимости к стационарному решению неявный сглаживатель исчезает, т.к. ui 0. Если использовать в явной части схемы ту же аппроксимацию потоков и источников, что и в базовом явном методе, то стационарное решение, полученное по схеме (3), совпадет с решением, удовлетворяются все требования к неявной схеме, которые были поставлены в пункте 2.1.

Для решения СЛАУ (20) в схеме [14][15] используется поблочный итерационный метод Гаусса-Зейделя [33]. Это просто означает, что для приращений параметров в соседних ячейках ( uNBR ) используется последние найденные в ходе итераций значения. Поэтому для нахождения нового итерационного значения приращения параметров в текущей ячейке ( ui ) достаточно обратить матрицу Ri размером 77. Для этого применяется метод LU-разложения с неполным выбором ведущего элемента [100].

Использование поблочного метода Гаусса-Зейделя обеспечивает минимальный рост вычислительных затрат на одну итерацию по сравнению с явной схемой.

Особенности реализации неявной схемы В ходе практической реализации неявного численного метода использовались лишь основные идеи схемы [14][15], описанные в предыдущем разделе. Построенный неявный метод в своих деталях неизбежно отличается от метода, реализованного, например, в программе КоБра [15], - по следующим причинам:

При построении явного оператора схемы (22) используются те же способы аппроксимации потоков и источниковых членов, что и в базовом явном численном методе, который реализован в пакете прикладных программ EWT-ЦАГИ (см. раздел 2.2) и который содержит отличия от методов, используемых в ЦИАМ.

Численный метод расписывается применительно к системе уравнений Рейнольдса, замкнутой моделью турбулентности ( q ), которая не используется в ЦИАМ.

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

Перечислим важнейшие особенности реализации неявной схемы.

Особенности аппроксимации конвективных потоков Вектор диффузионных потоков через грань ячейки Ei имеет следующий вид (см. (17) и (4)):

где s ( sx ; s y ; sz ) SEi n - вектор внешней нормали к грани ячейки, конвективный поток массы через грань ячейки, H E p / - полная энтальпия единицы массы газа.

При построении неявного сглаживателя используется линеаризация потоков – см. (19). Для линеаризации конвективных потоков использован известный метод Роу [39].

Конвективные потоки через грань Ei, согласно [14][15], запишем в линеаризованном виде:

где матрица Роу AEi и вектор CEi считаются постоянными на данной грания ячейки в течение всего шага по времени. Матрицы Роу для конвективных потоков системы уравнений Рейнольдса, замкнутой моделью турбулентности ( q ), выписаны в Приложении 3.

Тогда приращение вектора конвективных потоков будет равно Величины на гранях uEi, uEi1 и их приращения uEi uEi1 uEi могут быть найдены из решения задачи Римана о распаде произвольного разрыва методом Роу [39]:

В этих формулах u L и u R - векторы консервативных параметров слева и справа от грани ячейки (направление L R совпадает с направлением увеличения индексов структурированной сетки);

sign A P sign P, где PP AEi ; P - матрица, столбцами которой являются собственные векторы матрицы Роу AEi ; - диагональная матрица собственных чисел матрицы AEi. Тогда приращение вектора u представится в виде:

Чтобы получить неявную схему первого порядка, будем брать u L и u R из центров ячеек, прилегающих к данной грани. Для левых граней (внешняя нормаль к которым направлена в сторону уменьшения индексов ячеек структурированной сетки) возьмем uL Ei uNBR, uR Ei ui, а для правых граней (внешняя нормаль к которым направлена в сторону увеличения индексов ячеек структурированной сетки) - uL Ei ui, uR Ei uNBR. (Напомним, что индекс NBR обозначает соседнюю ячейку, имеющую общую грань Ei с ячейкой i.) AEi P | | P 1, | | - диагональная матрица, у которой на диагонали стоят модули собственных чисел матрицы Роу AEi.

Пусть Ei - грань ячейки, внешняя нормаль к которой направлена в сторону увеличения индексов ячеек структурированной сетки. Для таких граней получим следующую аппроксимацию конвективных потоков:

Аналогично, пусть Ei - грань ячейки, внешняя нормаль к которой направлена в сторону уменьшения индексов ячеек структурированной конвективных потоков:

Введем параметр sig E, который равен (+1) для граней Ei и (-1) для граней Ei.Тогда аппроксимации (26), (27) можно записать одной формулой:

Следует подчеркнуть, что аппроксимация (28) имеет вид (19).

Матрица Роу AEi по структуре совпадает с матрицей Якоби конвективных потоков. Все ее элементы могут быть вычислены через параметры из набора f u; v; w; H ; q;. Эти параметры вычисляются путем специального осреденения [39] параметров на явном (известном) временном слое n в ячейках, расположенных слева и справа от грани Ei :

где f - любой параметр из набора f u; v; w; H ; q; T.

По матрицам AEi однозначно вычисляются матрицы AEi P | | P 1.

Удобно вычислять не матрицу Роу AEi, а сразу собственные векторы и собственные числа. Все компоненты этих матриц вычисляются на явном слое. Матрицы выписаны в Приложении 3.

аппроксимируются при помощи нелинейной схемы Колгана 2-го порядка по пространству [7]. Вместо функции-лимитера Колгана (minmod) используется менее диссипативный лимитер Ван-Лира [101] в векторной версии [10] [102]. Эта аппроксимация конвективных потоков совпадает по структуре с аппроксимацией конвективных потоков на шаге-корректоре базовой явной схемы [10]. Отличие между ними заключается в том, что в базовой явной схеме эта аппроксимация применяется на слое (n 1 / 2), а здесь – на явном слое n. При сходимости к стационарному решению, аппроксимации становятся одинаковыми.

Особенности аппроксимации диффузионных потоков Вектор диффузионных потоков через грань ячейки Ei имеет следующий вид (см. (17) и (5)):

где s ( sx ; s y ; sz ) SEi n - вектор площади-нормали к грани ячейки, Tnk Txk sx Tyk s y Tzk sz.

(неизвестном) временном слое (n 1) в виде:

аппроксимация 2-го порядка по пространству, которая по структуре совпадает с аппроксимацией этоих потоков в базовой явной схеме [10], только все параметры берутся на явном (известном) слое n.

В величины I in, n, Tnq и Tn входят коэффициенты молекулярной вязкости и турбулентной вязкости T, “турбулентное давление” компоненты усеченного вектора примитивных переменных f (u; v; w; T ; q; )T. При построении аппроксимации для приращений за шаг по времени F дифф производится линеаризация величин I in, n, Tnq и молекулярной и турбулентной вязкости, а также “турбулентное давление” замораживаются на явном слое n и вычисляются так же, как в базовом явном методе [10]. Диффузионные потоки кинетической энергии турбулентности Tnk, которые обычно малы по сравнению с остальными потоками энергии, замораживаются на слое n целиком, так что их приращения полагаются равнми нулю.

диффузионных потоков импульса) при конвективном переносе газа со скоростью u; v; w. Линеаризация этого члена позволяет представить его приращение в виде Поскольку величины u E, vE, wE связаны с конвективным переносом газа, их следует вычислять так же, как и при вычислении приращений конвективных потоков, то есть по методу Роу.

целесообразно представить в виде где диффузионных потоков, входящих в F1дифф, на примере приращения диффузионного потока x -компоненты импульса I xn. При вычислении молекулярной и турбулентной вязкости, а также “турбулентное давление” замораживаются на явном слое n :

Рассмотрим криволинейную систему координат (,, ), связанную с индексными направлениями структурированной сетки. В этой системе координат ось пересекает грань ячейки, а оси и касательны к грани. Тогда производные можно представить в виде:

криволинейным координатам, касательным к граням. Это позволяет представить диффузионные потоки в виде (19) и тем самым сокращает шаблон схемы на неявном слое до 7 ячеек. Подобный подход использован, например, в программе EDGE (Швеция, FOI) [103]. Для запишем простейшую центрально-разностную аппроксимацию 1-го порядка точности. Получим:

где - расстояние между центром ячейки и центром грани Ei.

Параметр sig Ei равен (+1) для граней Ei, внешняя нормаль к которым направлена в сторону увеличения индексов ячеек структурированной сетки, и равен (-1) для граней Ei, внешняя нормаль к которым направлена в сторону уменьшения индексов ячеек структурированной сетки. Подставив (33) в (32), получим для I xn I xn (t n 1 ) I xn (t n ) следующую аппроксимацию:

Аналогично конструируются аппроксимации для I yn, I zn. Таким же способом получаются и приращения потоков тепла и параметров турбулентности:

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

Здесь - матрица преобразования от консервативных переменных u ; u; v; w; E; q; к примитивным переменным f (u; v; w; T ; q; )T, а A - матрица коэффициентов, зависящих от криволинейных координат. Заметим, что матрица преобразования не квадратная - так же, как и матрица A. Квадратной матрицей будет их произведение. Матрицы приведены в Приложении 3.

Теперь рассмотрим устройство F2дифф. Величины ( I xn )Ei, ( I yn )Ei и ( I zn )Ei берутся из аппроксимации диффузионных потоков на явном слое.

Как было сказано выше, величины u E, vE, wE следует вычислять по методу Роу. Поэтому следует извлечь эти величины из вектора приращений консервативных переменных uEi, вычисленного по формуле (25). Тогда Особенности аппроксимации источниковых членов источниковых членов имеет вид W 0, 0, 0, 0, 0, S q, S. Выражения для S q и S выписаны в разделе 1.1 – см. (11). Будем вычислять коэффициенты Aq, Bq, Cq, A, B, C по параметрам на известном временном слое n. Тогда W оказывается линейной функцией от двух последних компонент вектора консервативных переменных u - от q и от (см. (11)). Эту функцию для ячейки i можно представить следующим образом:

где матрица Якоби источниковых членов имеет вид:

числа этой матрицы.

Схема (18) аппроксимирует систему уравнений (13) в ячейке i на источниковых членов на слое (n 1), которое используется в схеме (18), можно предложить две аппроксимации: явную аппроксимацию и неявную аппроксимацию:

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

В Приложении 1 применительно к модельному уравнению показано, что неявная аппроксимация источниковых членов (38) не гарантирует абсолютной устойчивости схемы. Источники порождают экспоненциально растущие моды решения вида e k t, где k - собственные числа матрицы (см. раздел 1.2, формула (14)). Неявная аппроксимация (38) абсолютно устойчива только в случае, если все моды являются затухающими, т.е.

если q 0 и 0. Если есть растущие моды ( k 0 ), то неявная схема устойчива лишь при условии, что шаг по времени 1 / k. Наоборот, явная аппроксимация (37) абсолютно устойчива при описании растущих мод, но условно устойчива при описании затухающих мод.

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

PP 1, где - диагональная матрица собственных чисел матрицы W / u i. Тогда источниковые члены аппроксимируются следующим образом:

где S P P - матрица, соответствующая только отрицательным собственным числам матрицы W / u.

Смысл аппроксимации (39) можно пояснить на примере линейной системы уравнений без потоковых членов:

Для данной системы уравнений схема (18)+(39) имеет вид:

Умножив (40) слева на P 1 и введя замену переменных z P 1u, представим данную численную схему в виде:

Но эта система уравнений представляет собой набор независимых друг от друга скалярных уравнений для каждой компоненты вектора z :

устойчива). Получаем в линейном случае абсолютно устойчивую схему.

впоследствии автор обнаружил, что такой же способ был использован в работе [43].

Матрица S выписана в Приложении 3.

Особенности решения системы линейных уравнений диффузионных потоков (31), (34), (35) и источниковых членов (39) в (18), представим неявную схему в виде (20), где Как было сказано в разделе 2.3, система линейных алгебраических уравнений (20) решается при помощи поблочного итерационного метода Гаусса-Зейделя.

В общем случае, матрица системы (20) не удовлетворяет известным достаточным условиям сходимости метода Гаусса-Зейделя [104]. Однако на каждом шаге по времени решение системы линейных уравнений можно не доводить до сходимости. При этом аппроксимация нестационарных уравнений Рейнольдса теряется, но в конце концов все равно удается стационарному решению. Стационарное решение не зависит от точности выполнения шага по времени и будет аппроксимировать стационарное решение уравнений Рейнольдса со 2-м порядком. Подобный подход использован и в [14], [15].

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

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

Расчеты ведутся с небольшими значениями числа Куранта. На рис.3,а представлено решение, которое получается, если выполнять на каждом шаге по времени несколько однонаправленных итераций (с одинаковым направлением обхода ячеек). Видно, что решение симметричной задачи стало асимметричным. Если же менять направление обхода ячеек на каждой итерации, то симметрия практически восстанавливается.

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

Рисунок 3. Зависимость решения от направления обхода ячеек Рисунок 4. Зависимость решения от числа Куранта Рис.4 показывает, что качество описания нестационарного процесса быстро ухудшается при увеличении числа Куранта. Но этого качества и не требуется от схемы, которая предназначена для решения стационарных задач.

В пространственных задачах для получения того же результата приходится использовать шесть итераций метода Гаусса-Зейделя на каждом шаге по времени:

итерация с последовательным обходом сечений j const, k const в положительном направлении индексной оси i (здесь i, j, k - индексы, структурированной сетке);

итерация с последовательным обходом сечений j const, k const в отрицательном направлении индексной оси i ;

итерация с последовательным обходом сечений i const, k const в положительном направлении индексной оси j ;

итерация с последовательным обходом сечений i const, k const в отрицательном направлении индексной оси j ;

итерация с последовательным обходом сечений i const, j const в положительном направлении индексной оси k ;

итерация с последовательным обходом сечений i const, j const в отрицательном направлении индексной оси k.

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

Вычислить поток на явном временном слое Fn i через грань ячейки, расположенную на границе блока расчетной области. Для определения этого потока используются те же численные граничные условия, что и в базовом явном численном методе. Эти граничные условия подробно описаны в работе [10] и здесь рассматриваться не Определить приращение параметров в отсутствующей соседней ячейке. Пусть это ячейка с номером NBR NBRb. Тогда в формулу (20) подставляется uNBRb, рассчитанное особым образом, в зависимости от типа граничного условия.

Опишем специфические детали реализации каждого типа граничного условия, которые предусмотрены в пакете прикладных программ EWT-ЦАГИ.

1. Граница «joint» - граница регулярной стыковки двух соседних блоков расчетной области (с непрерывными на границе сеточными линиями). На каждой итерации метода Гаусса-Зейделя последовательно обходятся все блоки расчетной области, и при работе в приграничной ячейке i в качестве uNBRb берется последнее найденное значение u из соответствующей приграничной ячейки соседнего блока. Таким образом, метод Гаусса-Зейделя работает так же, как если бы совместно решалась система линейных уравнений для всех ячеек всех блоков расчетной области. Влияние многоблочной структуры сетки проявляется лишь в специфическом порядке обхода ячеек (с последовательным обходом блоков сетки).

2. Граница “connect” - граница нерегулярной стыковки двух соседних блоков расчетной области (с непрерывными на границе интерполированными параметрами из соответствующей ячейки соседнего блока. В качестве uNBRb берутся интерполированные приращения параметров из ближайшей приграничной ячейки соседнего блока.

3. Граница “riemann” – внешняя граница расчетной области, удаленная от обтекаемых тел. Опыт показал, что на такой границе можно ставить явное граничное условие, т.е. полагать uNBRb 0. Это не влияет ни на устойчивость расчета, ни на скорость сходимости к стационарному решению.

Для определения uNBRb на границах остальных типов деалется предположение, что значения параметров в заграничной ячейке линейно связаны с параметрами в приграничной ячейке. Эта линейная связь может быть представлена в виде uNBRb Eb ui C, где Eb - некоторая постоянная матрица, а C - некоторый постоянный вектор. Отсюда следует uNBRb Eb ui. Подставив это соотношение в (20), получим следующее уравнение для приграничной ячейки i :

4. Граница “symmetry” – плоскость симметрии. В этом случае все параметры, кроме нормальной к границе компоненты скорости, должны совпадать с параметрами из приграничной ячейки, а нормальная к границе компонента скорости отличается знаком. Этому соответствует матрица где n (nx ; n y ; nz ) - единичный вектор внешней нормали к границе.

5. Граница “solid_insulated” – теплоизолированная твердая стенка с прилипанием потока. В этом случае предполагается, что отсутствующая соседняя ячейка NBRb имеет нулевой объем, а ее центр расположен на твердой стенке. Тогда все параметры в этой ячейке, кроме скорости и характерной величины турбулентных пульсаций скорости (параметра турбулентности q ), должны совпадать с параметрами из приграничной ячейки, а скорость и q должны быть равны нулю. Этому соответствует матрица, которая в безразмерных переменных записывается следующим образом:

(Используется система обезразмеривания, в которой компоненты скорости и величина q относятся к заданной характерной скорости потока газовая постоянная для воздуха, а характерная частота турбулентности к V* / L*, где L* 1 м.) Комбинированный метод с явной и неявной частями При проведении расчета с заданным глобальным шагом по времени t 0 неявная схема (20) оказывается эффективной не во всей расчетной области. Для примера на рис.5 показана сетка для расчета обтекания модели CRM (Common Research Model, которая рассматривалась на семинаре 4th AIAA CFD Drag Prediction Workshop в 2009 г.). При такой структуре расчетной сетки (которая типична для практических расчетов) максимальное обобщенное число Куранта (см. раздел 2.2), определяемое мельчайшими ячейками, прилегающими к поверхности модели, составляет CFLmax 105. При этом почти во всей расчетной области, кроме тонкого слоя вокруг модели, локальное число Куранта оказывается меньшим единицы. В этих областях более эффективным с вычислительной точки зрения является расчет по явной схеме, т.к. она не требует вычисления компонент матриц Ri, RNBR и обращения в каждой ячейке системы из линейных уравнений. Эффективность явной схемы с точки зрения скорости сходимости к стационарному решению можно еще повысить, применяя локальный шаг по времени [10] [99]. Напротив, при проведении расчетов с использованием неявной схемы с большими значениями числа Куранта в окрестности обтекаемых тел применение локального шага оказывается в общем случае ненадежным, т.к. нарушения законов сохранения, обусловленные локальным шагом по времени, при больших значениях числа Куранта оказываются слишком большими. Это может приводить к плохой сходимости или даже к неустойчивости расчета.

Рисунок 5. Расчетная сетка для моделирования течения Эти соображения приводят к идее комбинированного метода, который и выдвигается на защиту в данной работе.

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

В тех ячейках, где локальное число Куранта CFL больше единицы, применяется неявная схема с заданным глобальным шагом t 0. Это обычно области сгущения сетки, где расчет по явной схеме шел бы с очень малым шагом по времени iуст и был бы неэффективен.

В тех ячейках, где выполняется условие 0.1 CFL 1, применяется базовая явная схема с локальным шагом по времени, равным t iуст.

Это заметно сокращает объем вычислений в ячейке. К тому же обеспечивается максимально возможная величина шага по времени в ячейке, а значит, более быстрая сходимость, чем при использовании неявной схемы с t 0.

Наконец, в буферных зонах, удаленных от обтекаемых тел, где нужно как можно быстрее прогнать и погасить все возмущения, применяется неявная схема с большим и локальным шагом по времени, равным t 30 iуст. Эта группа ячеек выделяется условием CFL 0.1.

На рис.6 для расчета обтекания модели CRM показаны области использования разных типов схем.

Рисунок 6. Области использования разных типов схем комбинированный метод может быть реализован без существенных изменений распараллеленного численного алгоритма.

Чтобы превратить неявную схему (20) в базовую явную схему, достаточно выполнить следующие действия:

Положить в неявном операторе Ri I ( I - единичная матрица), использованием параметров на явном (известном) временном слое n, определить параметры на временном слое (n 1 / 2) и вычислить явную часть схемы (сумму потоков и источниковых членов) по этим параметрам.

Порядок вычислений при использовании комбинированного метода показан на рис.7.

Рисунок 7. Порядок вычислений в комбинированном методе Поскольку совершение шага по времени по явной схеме производится с использованием только параметров на известном временном слое n, а ui находится за одну итерацию, то вычисления по явной схеме проводятся лишь в ходе первой из шести итераций метода Гаусса-Зейделя.

На остальных итерациях “явные” ячейки пропускаются.

Рассмотрим работу комбинированного метода на границах стыковки ячеек, в которых используются явная и неявная схемы. Пусть на данном шаге по времени вычисления в ячейке с номером i ведутся по явной схеме, а в соседней ячейке с номером NBR - по неявной. Тогда поток через общую грань этих двух ячеек вычисляется в ячейке i на временном слое (n 1 / 2), а в ячейке NBR - на временном слое n. Если теперь учесть, что расчет в “явной” ячейке ведется с локальным шагом по времени, равным t iуст, в во всех “неявных” ячейках первого типа – с глобальным шагом по времени, равным t 0, то станет ясно, что поток через общую грань “явной” и “неявной” ячеек вычисляется в ячейке i в момент времени t t n iуст, а в ячейке NBR - в момент времени t t n 0. Эти моменты времени, вообще говоря, не совпадают. Поэтому значения потока через эту грань в “явной” и “неявной” ячейках будут различными и, следовательно, расчет будет проводиться с нарушением консервативности (количества физических величин, выносимые через данную грань из ячейки i, будут отличаться от количеств этих величин, вносимых через данную грань в ячейку NBR ). Однако следует напомнить, что в расчете по неявной схеме итерации метода Гаусса-Зейделя не доводятся до сходимости. Поэтому аппроксимация нестационарного процесса все равно теряется. Теряется и она в расчете по явной схеме с локальным шагом по времени (т.е. с разными шагами по времени в разных ячейках). Поэтому нельзя даже считать, что к моменту начала шага по времени обе ячейки находились на одном временном слое t n. В этих условиях требование консервативности на этапе сходимости становится несущественным, а основное значение имеют устойчивость алгоритма и его сходимость к стационарному решению. При достижении же стационарного решения зависимость от “времени” исчезает, и потому формулы для потока через данную грань, используемые в “явной” и “неявной” схеме, совпадают. Таким образом, если стационарное решение достигается, то численное решение становится консервативным.

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

Только после того, как комбинированный метод был полностью реализован и успешно применен к решению ряда практических задач, автору стало известно, что идея перехода от неявной схемы к явной, в зависимости от локального условия устойчивости, использовалась также в серии работ Ю.Б.Радвогина и Н.А.Зайцева из ИПМ им. М.В.Келдыша РАН [45-48]. Нужно подчеркнуть, что, несмотря на использование сходной идеи, комбинированный метод, предлагаемый в настоящей работе, содержит ряд существенных особенностей, отличающих его от метода Радвогина и Зайцева:

Метод Радвогина и Зайцева был разработан и исследован только применительно к гиперболическим дифференциальным уравнениям в частных производных (в частности, к уравнениям Эйлера).

В работах Радвогина и Зайцева рассматривается совершенно другой тип неявной схемы, связанный с расщеплением по пространственным направлениям и сведением задачи к последовательному решению одномерных задач. В работе [46] проведен анализ устойчивости схемы и показано, что для достижения абсолютной устойчивости такой схемы необходимо осреднять решение, полученное за одну итерацию по явной и неявной схемам, с весами (для неявной схемы) и. По опыту использования комбинированного метода, который предложен в настоящей работе, устойчивость решения сохраняется, по крайней мере, вплоть до числа Куранта CFL~106 без использования дополнительного сглаживания.

Метод Радвогина и Зайцева предназначен для эффективного решения нестационарных задач. Поэтому в работах Радвогина и Зайцева при CFL1 используется неявная схема с глобальным шагом по времени, а при CFL1 – явная схема также с глобальным шагом. В настоящей работе предложен метод для решения стационарных задач. Поэтому предлагается использовать явную схему с локальным шагом по времени и, кроме того, при CFL0.1 переходить на неявную схему с локальным шагом по времени, соответствующим CFL=30.

Глава 3. Тестирование разработанного метода Описанный выше метод был реализован в виде программы-солвера COMGLEI (Combination of Global and Local tau type with Explicit and Implicit schemes). Эта программа является модификацией неявного блока программы ZEUS, разработанного совместно автором настоящей статьи и С.В. Михайловым. После модификаций удалось добиться совпадения результатов программы COMGLEI с результатами промышленного солвера V3Solver, в котором реализован базовый явный численный метод.

И ZEUS, и V3Solver входят в пакет прикладных программ EWT-ЦАГИ.



Pages:   || 2 | 3 |
 
Похожие работы:

«Чигиринский Юлий Львович ОБЕСПЕЧЕНИЕ ТОЧНОСТИ И КАЧЕСТВА ПОВЕРХНОСТЕЙ ПРИ МНОГОПЕРЕХОДНОЙ МЕХАНИЧЕСКОЙ ОБРАБОТКЕ НА ОСНОВЕ СОВЕРШЕНСТВОВАНИЯ ИНФОРМАЦИОННЫХ И МАТЕМАТИЧЕСКИХ СРЕДСТВ ПРОЕКТИРУЮЩЕЙ ПОДСИСТЕМЫ САПР ТП 05.02.08 – Технология машиностроения 05.13.06 – Автоматизация и управление технологическими процессами и производствами (в машиностроении) диссертация на...»

«ГОРЕЛКИН Иван Михайлович РАЗРАБОТКА И ОБОСНОВАНИЕ СПОСОБОВ ПОВЫШЕНИЯ ЭНЕРГОЭФФЕКТИВНОСТИ НАСОСНОГО ОБОРУДОВАНИЯ КОМПЛЕКСОВ ШАХТНОГО ВОДООТЛИВА Специальность 05.05.06 – Горные машины Диссертация на соискание ученой степени кандидата технических наук Научный руководитель...»

«Карапузова Марина Владимировна УДК 621.65 ГИДРОДИНАМИЧЕСКИЕ ОСОБЕННОСТИ КОНСТРУИРОВАНИЯ КОМБИНИРОВАННОГО ПОДВОДА ЦЕНТРОБЕЖНОГО НАСОСА Специальность 05.05.17 – гидравлические машины и гидропневмоагрегаты Диссертация на соискание научной степени кандидата технических наук Научный руководитель Евтушенко Анатолий Александрович канд. техн. наук, профессор Сумы – СОДЕРЖАНИЕ ПЕРЕЧЕНЬ...»

«УДК 622.673.4:621.625 Васильев Владимир Иванович ОБОСНОВАНИЕ РАЦИОНАЛЬНЫХ ДИНАМИЧЕСКИХ ПАРАМЕТРОВ ПРЕДОХРАНИТЕЛЬНОГО ТОРМОЖЕНИЯ ШАХТНЫХ ПОДЪЕМНЫХ УСТАНОВОК Специальность 05.02.09 – динамика и прочность машин Диссертация на соискание научной степени кандидата технических наук Научный руководитель – доктор технических наук, профессор В. М. Чермалых Киев - СОДЕРЖАНИЕ...»






 
© 2013 www.diss.seluk.ru - «Бесплатная электронная библиотека - Авторефераты, Диссертации, Монографии, Методички, учебные программы»

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