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

( Numerical study of far-field non-local boundary conditions for transonic flow around wing
Preprint, Inst. Appl. Math., the Russian Academy of Science)

Введение

               Работа по созданию алгоритма трехмерных нелокальных граничных условий для задач дозвукового обтекания была начата ранее Софроновым [1], где разработана теоретическая база создания искомых граничных условий и продолжена Софроновым и Нажесткиной [2], где представлены результаты по первой версии программы. В основе метода лежит аналитическое представление течений, описываемых линеаризованными на однородном набегающем потоке уравнениями Эйлера, которые, как предполагается, служат адекватной моделью аэродинамических течений вдали от обтекаемого объекта.  Вопрос численной реализации представляет определенную трудность из-за необходимости аккуратной и эффективной трактовки возникающих поверхностных сингулярных интегралов. Предложенные в [2] алгоритм и компьютерная программа были модифицированы с целью охвата случая произвольной формы выходной границы расчетной области. Новая версия комплекса программ расчета обтекания крыла реализована для многопроцессорной машины MBC-1000 и эффективно распараллеливается в системах MPI и DVM [3]. Этот комплекс включает в себя программу интегрирования уравнений Эйлера внутри расчетной области и программу вычисления нелокальных граничных условий. Первая из указанных программ предоставлена А.Е.Луцким и  использует  разностную схему, аналогичную схеме С.К.Годунова [4]. Стационарное решение находится методом установления. Вторая программа реализована как относительно независимый блок, поставляющий газодинамические функции на внешний слой сетки для очередного шага по времени. 

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

               Отметим, что другой тип нелокальных граничных условий, а именно, полученных на основе метода разностных потенциалов [6], применялся для задач обтекания вязким газом, см., например [7].

Автор благодарен И.Л. Софронову за обсуждение и помощь в исследовании.

Результаты для числа Маха 0.72

               Для проведения расчетов были построены две сетки – мелкая, объемом 276´44´48 ячеек, и грубая, путем разряжения через одну ячейку по каждому направлению. За основу построения мелкой сетки был взят файл из Интернета с сеткой для расчета вязких течений около крыла ONERA-M6. Исходная сетка была переинтерполирована с целью устранения избыточного сгущения к поверхности крыла (поскольку мы рассчитываем невязкое обтекание) и поворота крыла для введения угла атаки 3 градуса. Количество ячеек в разных направлениях было немного уменьшено. Затем эта сетка была разделена на блоки, для  распараллеливания на задаваемое число  процессоров. Расчеты велись на вычислительном комплексе rsc4.kiam.ru.

               Для выявления влияния размера расчетной области на точность решения расчеты проводились на последовательности уменьшающихся вычислительных областей (путем отбрасывания внешних слоев сетки). На Рис. 1 показаны крыло и две сетки, соответствующие самой большой  и самой маленькой  областям.

               Рассматривалась задача обтекания крыла ONERA-M6 при числе Маха 0.72 и угле атаки 3 градуса. Для каждой сетки проводились два расчета: с нелокальными и с характеристическими граничными условиями. Отметим, что для самой большой расчетной области оба решения совпадают с тремя знаками. Среднее арифметическое этих двух решений принималось за эталонное решение.

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

               На Рис. 2 приведены изолинии давления в плоскости симметрии z=0 (т.е. в нижней плоскости из Рис. 1), полученные для самой маленькой расчетной области. Хорошо видно, что решение с нелокальными граничными условиями (слева) гораздо ближе к эталонному решению.  Этот вывод справедлив и для распределения давления на крыле, см. Рис. 3.

               На Рис. 4 приведены значения коэффициента подъемной силы в зависимости от числа отбрасываемых внешних слоев сетки. Коэффициент подъемной силы эталонного решения, а также отклонение от него на ±0.5% изображены горизонтальными линиями. Можно видеть, что нелокальные граничные условия  сохраняют хорошую точность расчета для достаточно малых расчетных областей. Например,  отбрасывание семи линий сетки соответствует сокращению линейного размера области в 3.5 раза.

               Однако полученная точность с нелокальными граничными условиями является, на наш взгляд, еще неудовлетворительной, поскольку алгоритм для двумерных задач демонстрирует более существенную эффективность, см. Приложение. По-видимому, имеется определенное несогласование между методом интегрирования уравнений Эйлера внутри области и методом постановки условий дальнего поля.  В частности, эти методы имеют разный порядок аппроксимации: первый и второй, соответственно.  На несогласованность этих двух методов указывает также и поведение процесса сходимости коэффициента подъемной силы в зависимости от времени. На  Рис. 5 представлены графики Су для расчетов в большой и маленькой областях с характеристическими и нелокальными граничными условиями. Хорошо видна более неравномерная сходимость для нелокальных граничных условий. Это указывает на необходимость дальнейшей модернизации разрабатываемого подхода.

Рис. 1: Большая и малая расчетные области для крыла ONERA-M6

Рис. 2: Изолинии давления на плоскости симметрии z=0. Слева – расчет с нелокальными ГУ (NBC), справа – с характеристическими ГУ (CBC). Сплошные линии соответствуют эталонному решению в большой области. Пунктирные линии и заливка соответствуют расчетам  в маленькой области. Граница расчетной области изображена на нижних картинках. Верхние картинки дают более подробное представление о решении вблизи крыла. Все изолинии приведены с одинаковым шагом.

               Аналогичная серия расчетов на мелкой сетке не дала существенного улучшения результата. Практически произошел небольшой, почти одинаковый сдвиг для значений аэродинамических коэффициентов в решениях с нелокальными и характеристическими граничными условиями. Для сравнения решений мы приводим Рис. 6 и  Рис. 7, которые аналогичны Рис. 2 и Рис. 3, соответственно, но на мелкой сетке.

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

Рис. 3: Изолинии давления на верхней поверхности крыла. Слева – расчет с нелокальными ГУ (NBC), справа – с характеристическими ГУ (CBC). Сплошные линии соответствуют эталонному решению в большой области. Пунктирные линии и заливка соответствуют расчетам  в маленькой области.

Рис. 4: Коэффициент подъемной силы Cy в зависимости от числа отбрасываемых внешних поверхностей сетки (от размера области). Расчетам с нелокальными ГУ соответствует  линия NBC.  Расчетам с характеристическими ГУ– линия CBC.  Горизонтальные линии показывают Cy и его 0.5%-ое отклонение для эталонного решения.

Рис. 5: Установление коэффициента подъемной силы Cy по времени для большой и маленькой областей.  Расчетам с нелокальными ГУ (NBC) соответствуют  пунктирные линии.  Расчетам с характеристическими ГУ (CBC) – сплошные линии.

Рис. 6:  Расчет на мелкой сетке. Изолинии давления на плоскости симметрии z=0. Слева – расчет с нелокальными ГУ (NBC), справа – с характеристическими ГУ (CBC). Сплошные линии соответствуют эталонному решению в большой области. Пунктирные линии и  заливка соответствуют расчетам  в маленькой области.

Рис. 7: Расчет на мелкой сетке. Изолинии давления на верхней поверхности крыла. Слева – расчет с нелокальными ГУ (NBC), справа – с характеристическими ГУ (CBC). Сплошные линии соответствуют эталонному решению в большой области. Пунктирные линии и заливка соответствуют расчетам  в маленькой области.

Рис. 8: Изолинии давления, вид сверху на расчетную область. Внизу показана поверхность на предпоследнем слое сетки при расчете в маленькой области. Вверху – та же поверхность, вырезанная из расчета в большой области. Слева – расчет с нелокальными ГУ (NBC), справа – с характеристическими ГУ (CBC). Видно, что оба расчета в большой области совпадают (верхний ряд). Очень близок к ним и расчет в маленькой области с нелокальными ГУ  (слева внизу). В то же время заметны различия для расчета в маленькой области с характеристическими ГУ (внизу справа). Расчет на мелкой сетке. (Белая полоска непрорисована из-за разрыва топологии графического представления сетки)

Рис. 9: Изолинии давления, вид сзади на расчетную область (подветренная сторона). Внизу показана поверхность на предпоследнем слое сетки при расчете в маленькой области. Вверху – та же поверхность, вырезанная из расчета в большой области. Слева – расчет с нелокальными ГУ (NBC), справа – с характеристическими ГУ (CBC). Видно, что оба расчета в большой области совпадают (верхний ряд). Близок к ним и расчет в маленькой области с нелокальными ГУ  (слева внизу). В то же время заметны различия для расчета в маленькой области с характеристическими ГУ (внизу справа). Расчет на мелкой сетке. (Белая полоска непрорисована из-за разрыва топологии графического представления сетки)

Результаты для других чисел Маха

Здесь описаны расчеты обтекания крыла ONERA-M6 с использованием нелокальных и  характеристических граничных условий. Были рассмотрены числа Маха 0.80 и 0.85 (угол атаки 3 градуса), крупная сетка 138´22´24 и мелкая сетка 276´44´48 ячеек.

А. Расчеты с числом Маха 0.80

1. Зависимость величины коэффициента подъемной силы Су от размера расчетной области (усечения) для нелокальных искусственных граничных условий (НИГУ) и характеристических граничных условий (ХГУ) показана в таблице 1, расчеты на крупной сетке. За эталонное решение принималось среднее арифметическое двух решений для неусеченной области, т.е. Су=(0.06111 +  0.06097)/2 = 0.06104.

Таблица 1: Величина Су и разность с эталонным решением; крупная сетка.

               Усечение «6» соответствует уменьшению линейных размеров расчетной области примерно в 3 раза. Мы наблюдаем, что точность расчета с НИГУ при таком усечении превосходит в 4 раза точность расчета с ХГУ и составляет 0.4%.

2. Расчеты на мелкой сетке представлены в таблице 2. Коэффициент подъемной силы принимался равным Су=(0.06172 +  0.06160)/2 = 0.06166.  Усечения расчетной области проводятся с удвоенным числом поверхностей сетки, что соответствует тем же линейным размерам, что и для крупной сетки. 

Таблица 2: Величина Су и разность с эталонным решением; мелкая сетка.

               Здесь выигрыш в точности не столь существенен (по сравнению с крупной сеткой) и составляет около 2-х раз.

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

      Су=0.06233      для НИГУ, и

Су=0.06223      для ХГУ.

               Таблица 3 содержит погрешности расчетов для этих значений Су.

Таблица 3: Величина погрешности для асимптотических значений Су (по Ричардсону)

               Мы отчетливо наблюдаем улучшение точности расчета при уменьшении размера расчетной области для случая НИГУ. Это можно объяснить тем, что диссипация разностной схемы на крупных ячейках вблизи открытой границы довольно велика и, тем самым, ухудшает точность при использовании большой расчетной области (усечение 0); в то же время использование НИГУ для малых расчетных областей позволяет более аккуратно моделировать течение Эйлера в отброшенных крупных ячейках. 

               В то же время мы наблюдаем в таблице 3 ухудшение точности Су на малых расчетных областях при использовании традиционных ХГУ.

               Таким образом, данный асимптотический анализ показывает, что точность Су на малых областях примерно в 2 раза выше при использовании НИГУ по сравнению с ХГУ.

B. Расчеты с числом Маха 0.85.

Мы провели расчеты для большой области (усечение 0) и самой маленькой (усечение 6).

               Таблицы 4, 5, 6 содержат данные, аналогичные данным из таблиц 1, 2, 3.

Таблица 4: Величина Су и разность с эталонным решением; крупная сетка, М=0.85;

Эталонный Су=0.07721

Таблица 5: Величина Су и разность с эталонным решением; мелкая сетка, М=0.85;

Эталонный Су=0.07823

               Асимптотические значения коэффициента подъемной силы при М=0.85 составляют

      Су=0.07937      для НИГУ, и

Су=0.07916      для ХГУ.

Таблица 6: Величина погрешности  для асимптотических значений Су (по Ричардсону)

               Здесь асимптотический анализ показывает, что точность Су на малых областях примерно в 3 раза выше при использовании НИГУ по сравнению с ХГУ.

C. Асимптотический анализ для расчетов с числом Маха 0.72

Асимптотические значения коэффициента подъемной силы составляют

      Су=0.04565      для НИГУ, и

Су=0.04551      для ХГУ.

Таблица 7: Величина погрешности с для асимптотических значений Су (по Ричардсону), М=0.72

Как и для случая М=0.80, точность Су на малых областях примерно в 2 раза выше при использовании НИГУ по сравнению с ХГУ.

D. Комбинация НИГУ и ХГУ

Из таблиц 1, 2, 4, 5 видно, что погрешности  для  НИГУ и ХГУ имеют противоположные знаки. Например, полусумма соответствующих значений Су  находится ближе к эталонному значению. Кроме того, сходимость метода установления к решению с ХГУ быстрее, чем с НИГУ. Поэтому возникла идея попробовать комбинированный подход, а именно, использовать ХГУ, в которых значения величин «на бесконечности» берутся не из набегающего потока, а из значений посчитанных в алгоритме НИГУ. Была проведена модификация алгоритма и выполнены необходимые тестовые расчеты. Результат оказался практически тот же, что и при использовании исходных ХГУ. Как по точности, так и по скорости.

Выводы

Можно сделать следующие выводы по результатам моделирования обтекания:

·         Применение нелокальных граничных условий дальнего поля (НИГУ) в рассмотренной задаче обтекания крыла ONERA-M6 при Махах от 0.72 до 0.85 (угол атаки 3 градуса) позволяет использовать расчетные области сравнительно малого размера (средний радиус области составляет примерно 2 размаха крыла); при этом точность Су находится, как правило, в пределах 0.5%, что в два-три раза лучше, чем для характеристических граничных условий;

·         Причина более низкой эффективности НИГУ по сравнению с двумерным случаем кроется, на наш взгляд, в слишком высокой диссипативности разностной схемы, применяемой для интегрирования уравнений Эйлера  (схема первого порядка). Дальнейшее уточнение модели и более высокая эффективность НИГУ (как, например, в ранее исследованном двумерном случае) возможны, прежде всего, при повышении порядка точности разностной схемы интегрирования уравнений Эйлера внутри области с первого на второй.

Приложение. Результаты двумерных расчетов из [5]

Было рассмотрено три режима обтекания:

I. Дозвуковой режим, число Маха  M=0.3 , с малым углом атаки alpha=5 градусов;

II. Дозвуковой режим, число Маха М=0.3, с  углом атаки alpha=15 градусов;

III. Трансзвуковой режим, число Маха М=0.825, с углом атаки alpha=2 градуса.

               В качестве профиля был выбран руль Жуковского. На Рис.10 представлена сетка объемом 128*83 ячеек. Одна и та же базовая сетка 128*128 применялась во всех расчетах. Каждая сетка 128*K получалась из базовой путем отбрасывания внешних линий. Тем самым варьировался размер расчетной области D. Угол атаки задавался путем поворота сетки относительно вектора скорости набегающего потока, который всегда совпадал с направлением оси  x.

Рис

Рис.10: Расчетная область с размером сетки 128х83 (руль Жуковского)

               Для интегрирования уравнений Эйлера внутри D использовался сеточный метод второго порядка точности, описанный в [7].

               Примем за «эталонное» решение то решение, которое соответствует максимальному размеру расчетной области, 128*128. Известно, что величины коэффициентов подъемной силы  Cy и волнового сопротивления Cx очень чувствительны к точности алгоритма. Поэтому возьмем за меру отклонения решения от эталонного величину

               В Табл. 1.1 приведены результаты расчетов режима I. Здесь R/L – отношение среднего радиуса расчетной области к хорде профиля; N – число итераций, необходимое для установления четырех значащих цифр в Cx и Cy; T/T* – отношение времен расчета вариантов, где T* – время счета в такой области, что Cy совпадает с эталонным по четырем значащим цифрам. Мы видим, что НИГУ  обеспечивают точность в пределах 0.3% для областей с размером R/L вплоть до 0.68. Экономия времени счета при использовании малых областей составляет 11 раз. Сетка и график решения изображены на Рис.11.

               В Табл. 1.2 и на Рис.12 даны результаты для режима II. В этом примере решение сильно отклоняется от набегающего потока. Поэтому приемлемая погрешность линейной модели течения в дальнем поле достигается при бoльших размерах расчетной области по сравнению с режимом I. Здесь оптимальный размер R/L расчетной области составляет 1.1. Экономия времени счета по сравнению с максимальной областью составляет 5.3 раз.

               Табл. 1.3 содержит результаты расчетов для режима III. Мы видим, что области размером R/L=2.1 обеспечивают необходимую точность при использовании НИГУ. Для таких областей сетка типа O еще достаточно мелкая в области следа, что открывает возможность обоснованного применения таких сеток при расчетах трансзвуковых режимов обтекания. Экономия времени счета достигает 6  крат.

                                                                                  Рис.11                                                        

                                                                                  Рис.12                                                       

                                                                                  Рис.13                                                       

Список литературы

[1] И.Л.Софронов, Нелокальные искусственные граничные условия для задач трехмерного стационарного обтекания. // Матем. Модел. Т.10. №9, 1998, 64–86.

[2] Э.И.Нажесткина, И.Л.Софронов, Численная реализация граничных условий дальнего поля для задач трансзвукового аэродинамического обтекания. // Препринт ИПМ им. М.В.Келдыша РАН, №93, 2001.

[3] Н.Коновалов, В.Крюков, Параллельные программы для вычислительных кластеров и сетей. // Открытые системы,№3,2002.

[4] Численное решение многомерных задач газовой динамики // Под ред. Годунова С.К. “Наука”,М.,1976.

[5] И.Л.Софронов, Точные искусственные граничные условия для некоторых задач аэродинамики и дифракции // Дисс. д.ф.-м.н, Москва, ИПМ им. М.В.Келдыша РАН, 2000.

[6] В.С.Рябенький,  Метод  разностных  потенциалов  для  некоторых задач механики сплошных сред // М.: Наука, 1987.

[7] S.V.Tsynkov and V.N.Vatsa,  An improved treatment of external boundary for three-dimensional flow computations // AIAA Paper No. 97-2074, in: Proceedings of the 13th AIAA  Computational Fluid Dynamics Conference, Part 2, Snowmass Village, Colorado, 1997, pp. 1139–1149.

[8] И.Л.Софронов, Быстросходящийся метод решения уравнений Эйлера //  ЖВМиМФ, Т. 31, № 4, 1991, 575 – 591.