CC BY
212
Том ХЫУ
УЧЕНЫЕ ЗАПИСКИ ЦАГИ
2013
№ 3
УДК 533.6.011
РЕШЕНИЕ ОБРАТНОЙ ЗАДАЧИ ДЛЯ ПРОФИЛЯ В РАМКАХ УРАВНЕНИЙ НАВЬЕ — СТОКСА, ОСРЕДНЕННЫХ ПО РЕЙНОЛЬДСУ
А. Л. БОЛСУНОВСКИЙ, Н. П. БУЗОВЕРЯ, И. А. ГУБАНОВА, М. А. ГУБАНОВА
Разработан итерационный метод решения обратной задачи для профиля в сжимаемом потоке газа в рамках осредненных по Рейнольдсу уравнений Навье — Стокса, относящийся к классу методов остаточной коррекции. Деформация контура профиля, уменьшающая невязку между рассчитанным и целевым распределением давления, определяется с помощью решения обратной задачи для профиля в потенциальном потоке сжимаемого невязкого газа. Приведены примеры построения формы контура по заданному распределению давления (со скачками уплотнения и без них), демонстрирующие широкие возможности метода.
Ключевые слова: обратная задача, уравнения Навье — Стокса, ЛЫБУБ СЕХ, профиль, распределение давления
ВВЕДЕНИЕ
Обратные методы аэродинамики, определяющие геометрию элементов летательного аппарата (ЛА), в первую очередь крыла, по заданному (целевому) распределению давления, являются мощным инструментом аэродинамического проектирования. Они позволяют устранить или ослабить скачки уплотнения, снизить уровень возмущений потока в заданном месте, реализовать распределение давления, благоприятное для развития ламинарного или турбулентного пограничного слоя. Конечно, не всякое заданное распределение давления можно реализовать физически, так как в общем случае обратная задача является некорректной. Условия разрешимости точно сформулированы только для двумерных потенциальных течений [1—3], а для реальных течений — вязких, трансзвуковых и, в особенности, трехмерных — неизбежно использование инженерных подходов.
БОЛСУНОВСКИИ Анатолий Лонгинович
кандидат технических наук, начальник отдела ЦАГИ
БУЗОВЕРЯ Николай Петрович
кандидат технических наук, ведущий научный сотрудник ЦАГИ
ГУБАНОВА Ирина Анатольевна
инженер ЦАГИ
ГУБАНОВА Мария Анатольевна
младший научный сотрудник ЦАГИ
Один из наиболее распространенных подходов к построению обратных методов основывается на принципе остаточной коррекции, согласно которому деформация поверхности профиля определяется тем или иным способом через невязку между рассчитанным и заданным
с*р ) распределениями давления [4—13]. Задача решается итерационно путем попеременных вызовов прямого метода анализа и блока коррекции (рис. 1). Время работы блока коррекции геометрии обычно составляет лишь незначительную долю от времени прямого расчета. Вследствие итерационного характера решения задачи правило деформации поверхности может быть достаточно грубым, например, в [4, 5, 10] используются аналитические до- и сверхзвуковые решения по распределению давления на волнистой стенке, поэтому число необходимых итераций исчисляется десятками и даже сотнями. Естественно, чем более точно определяется коррекция геометрии, тем меньше потребуется общее число итераций.
В методах остаточной коррекции используют хорошо развитый аппарат решения прямых задач обтекания. Они не связаны с конкретным методом решения прямой задачи, а используют его как «черный ящик», так что последний легко может быть заменен впоследствии на более совершенный. Развитие прямых методов, как правило, опережает развитие чисто обратных методов
(основанных на решении задачи обтекания с граничными условиями типа ер = е*р), в результате
чего c помощью принципа остаточной коррекции можно решать обратные задачи для более широкого класса течений и компоновок.
Бурное развитие вычислительной аэродинамики в последнее время привело к широкому внедрению в практику аэродинамического проектирования прямых методов, основанных на решении осредненных уравнений Навье — Стокса (RANS-методы), замкнутых той или иной моделью турбулентности [14—19]. Такой подход представляет собой разумный компромисс между сложностью решаемых уравнений и потребными для этого вычислительными и временными ресурсами. Использование более простых уравнений хотя и ведет к существенному ускорению вычислений, не носит столь универсальный характер и требует тщательного обоснования ввиду опасности упустить какое-либо важное физическое свойство течения, влияющее, в конечном счете, на распределенные и суммарные аэродинамические характеристики ЛА.
Предлагаемый в данной работе итерационный метод решения обратной задачи для профиля в рамках уравнений RANS также принадлежит к классу методов остаточной коррекции. Основная идея метода заключается в том, что в качестве корректора, генерирующего поправки к текущим значениям координат контура профиля, используется метод решения обратной задачи для профиля в потенциальном потоке сжимаемого невязкого газа. Данный метод — TRAINV — был разработан в ЦАГИ много лет тому назад [12] и широко применяется как для проектирования изолированных скоростных профилей, так и в качестве корректора в 3-мерном методе решения обратной задачи для трансзвуковых крыльев [13].
1. ПОСТАНОВКА ЗАДАЧИ И ОПИСАНИЕ АЛГОРИТМА РАСЧЕТА
В данной работе решается обратная задача нахождения координат контура профиля по заданному распределению давления. Метод решения обратной задачи относится к классу методов остаточной коррекции. Для создания эффективного алгоритма используется принцип вложенных уровней, согласно которому в качестве блока расчета коррекции геометрии целесообразно применять метод обратной задачи более низкого уровня, где под уровнем подразумевается сложность решаемых уравнений или рассматриваемых форм. Типичными уровнями геометрии в аэродинамических задачах являются профиль, крыло, крыло + фюзеляж, полная компоновка. Типичными уровнями сложности уравнений являются уравнение Лапласа, уравнение полного потенциала, уравнения Эйлера, уравнения Навье — Стокса. Как правило, переход к следующему уровню увеличивает вычислительные затраты на один-два порядка.
Рис. 2. Принцип вложенных уровней
Рис. 3. Применение принципа вложенных уровней в данной задаче
Использование принципа вложенных уровней позволяет постепенно наращивать уровень решаемых обратных задач при сохранении ясного понимания физики явлений, а также рационально распределять расчеты на разных уровнях, с тем чтобы минимизировать общий объем необходимых вычислений. Применяя последовательно данный принцип, можно создать «-уровне-вые комплексы решения обратных задач для сложных компоновок и сложных уравнений, используя фактически только один простейший обратный метод на самом низшем уровне. Например, разработанные в ЦАГИ методы решения обратной задачи для крыла в трансзвуковом потоке газа — TRAWDES [13] состоят из трех уровней (рис. 2).
Предлагаемый в данной работе итерационный метод также состоит из трех уровней (рис. 3). Верхний уровень включает в себя метод прямого расчета обтекания профиля вязким сжимаемым потоком — программный пакет ANSYS CFX [20]. На втором уровне осуществляется прямой расчет обтекания профиля в сжимаемом невязком потоке методом полного потенциала. Первый (низший) уровень представляет собой метод решения прямой/обратной задачи для профиля в несжимаемой жидкости. Первый и второй уровни в совокупности составляют метод решения обратной задачи для профиля в сжимаемом невязком потоке— программу TRAINV [12].
Обратная задача для профиля решается в процессе итераций. На каждой к-й итерации выполняются следующие действия.
1. Строится сетка вокруг профиля, геометрия которого получена на предыдущей (к - 1)-й итерации, рассчитывается его обтекание при помощи CFX, определяется невязка распределения
давления 5ср, как разница между текущим и заданным распределениями давления.
2. Рассчитывается обтекание того же профиля в сжимаемом невязком потоке методом полного потенциала.
3. Формируется целевое распределение давления для профиля в сжимаемом невязком потоке путем прибавления разницы г5ср (где г < 1) к распределению давления с шага 2.
4. Решается обратная задача в сжимаемом невязком потоке по целевому распределению давления с шага 3 и получается к-я геометрия профиля.
Для полученной геометрии необходимо снова решить прямую задачу RANS, найти распределение давления и т. д. На каждой итерации вычисляется невязка давления
N
гер
N Т(СР, СРг
определяющая отклонение расчетного распределения давления от за-
данного. Процесс повторяется несколько раз до сходимости, либо до выхода на минимально достижимую невязку. Достаточная для практики сходимость соответствует уровню sCp < 0.01.
Таким образом, для решения обратной задачи по предложенному алгоритму требуется последовательный вызов различных пакетов программ. Для этой цели авторами был написан скрипт — командный файл на языке программирования Python.
2. КРАТКОЕ ОПИСАНИЕ ИСПОЛЬЗУЕМЫХ РАСЧЕТНЫХ МЕТОДОВ
2.1. МЕТОД РЕШЕНИЯ ПРЯМОЙ ЗАДАЧИ В СЖИМАЕМОМ ВЯЗКОМ ПОТОКЕ
Программный комплекс ЛК8У8 СБХ предназначается для решения широкого круга задач гидрогазодинамки, химической кинетики, горения, динамики многофазных потоков и др. [20]. Для конечно-разностной аппроксимации уравнений Навье — Стокса, осредненных по Рейнольд-су, используется метод конечного объема. Имеется возможность проводить расчеты как на многоблочных структурированных, так и на неструктурированных и гибридных сетках. В данной работе расчеты проводились на структурированных многоблочных сетках, позволяющих осуществить подробную дискретизацию области пограничного слоя и вязкого следа профиля.
Основные уравнения сохранения массы, импульса и энергии для каждого контрольного объема записываются в интегральной форме законов сохранения. Дискретизация конвективных слагаемых по пространству проведена с помощью противопотокового метода второго порядка аппроксимации. Для аппроксимации производных по времени в уравнениях используется неявная схема Эйлера первого или второго порядка.
Расчеты проводились на сетке «С» с размером ~100 тыс. ячеек. Сетка строилась при помощи сеткопостроителя 1СЕМ СББ. При изменении контура профиля в процессе итераций созданная начальная сетка автоматически перестраивалась для соответствующей геометрии. Высота
первой ячейки равна Ау ^ 10-6 (значения параметра у + вблизи стенки находятся в пределах от 0.3 до 0.6). Границы расчетной области удалены от профиля на 50 хорд. В качестве модели турбулентности выбрана двухпараметрическая дифференциальная модель 88Т. Течение полагалось полностью турбулентным.
2.2. РЕШЕНИЕ ПРЯМОЙ ЗАДАЧИ ОБТЕКАНИЯ ПРОФИЛЯ МЕТОДОМ ПОЛНОГО ПОТЕНЦИАЛА
Для расчета трансзвукового обтекания профиля потоком идеального газа применяется программа О. В. Карася, в основу которой положен многосеточный смешанный метод решения полного уравнения для потенциала скорости. Он включает в себя быстрый метод решения уравнения Пуассона для дозвуковой части течения и метод последовательной линейной сверхрелаксации для сверхзвуковой части. Решение ищется на сетке «О» в преобразованной плоскости переменных, получаемых путем последовательных комфорных преобразований с внешности профиля на внутренность единичного круга. Расчет ведется последовательно на трех сетках: 32 х 8, 64 х 15 и 128 х 15 ячеек. Расчет обтекания профиля на одном угле атаки требует приблизительно 0.1 секунды работы современного ПК. В связи с тем, что при решении используется неследящая разностная схема, применение программы ограничено умеренно закритическими режимами обтекания (местное число М < 1.25).
2.3. МЕТОД РЕШЕНИЯ ПРЯМОЙ/ОБРАТНОЙ ЗАДАЧИ В НЕСЖИМАЕМОЙ ЖИДКОСТИ
Для расчета обтекания профиля в несжимаемой жидкости применяется панельный метод с кусочно-линейным распределением плотности вихревого слоя на поверхности (рис. 4). Прямая и обратная задача решаются с идентичными граничными условиями и распределениями вихревых особенностей.
Условие постоянства функции тока в узлах контура профиля
= y. cos а - x. sin а - y(s)ín^[x. - x(5)]2 +[y. - y (5)]2ds
i = 1, 2, ..., N +1
= const,
приводит в прямой задаче к системе линейных алгебраических уравнений относительно неизвестных интенсивностей вихревого слоя у ^:
N
X ( - aN+1 з ) У] = ( - %+1) ^ а - (У - УN+1) С08 а 3=1
I = 1, 2, ..., N.
При решении обратной задачи для профиля в несжимаемой жидкости интегральное уравнение разрешается непосредственно относительно координат контура:
1
y =-
cos а
* с + Y(s )W (Xi - x (s ))2 + (y - y (s ))2 ds
- xi tg а
(1)
или в панельном представлении
1
У =-
cos а
N
* с -Z ^ Y J
J =1
+ X. tg а, i = 1,2, ..., N +1.
Это выгодно отличает подход с применением метода функции тока от других известных методов, в которых, как правило, определяются величины наклонов.
Уравнение (1) нелинейное, так как неизвестные значения у{ входят под знак интеграла. Решение ищется в процессе итераций, в результате определяется новая форма контура. На каждой итерации осуществляется преобразование координат в канонический вид и подправляется значение угла атаки. В случае самопересечения контура или чрезмерного раскрытия задней кромки осуществляется клиновидная добавка к ординатам профиля, выводящая затупление задней кромки на заданную величину. Обычно для сходимости процесса требуется 10—30 итераций. Поскольку решение обратной задачи для несжимаемой жидкости является лишь подэтапом решения общей задачи, добиваться полной сходимости нецелесообразно.
2.4. РЕШЕНИЕ ОБРАТНОЙ ЗАДАЧИ В СЖИМАЕМОМ ПОТОКЕ ПОТЕНЦИАЛЬНОГО ГАЗА
Решение обратной задачи для профиля в сжимаемом потоке газа осуществляется по принципу остаточной коррекции, где в качестве прямого метода выступает метод, описанный в разделе 2.2, а в качестве корректора — метод, описанный в разделе 2.3.
В случае докритического обтекания невязка давления с верхнего уровня передается непосредственно на нижний, возможно с некоторой релаксацией. После решения обратной задачи в несжимаемой жидкости геометрия контура сглаживается в районе критической точки и передается на верхний уровень для прямого расчета.
При наличии в области течения сверхзвуковых зон итерационный процесс необходимо видоизменить в связи с различным типом уравнений, описывающих до- и сверхзвуковое течение. Определение поправок к координатам поверхности профиля в области сверхзвукового течения должно отражать гиперболический характер уравнения.
В целом разность между рассчитанным и заданным распределением давления расщепляется на дозвуковую и сверхзвуковую части [6] (рис. 5):
дск = 5 ск + 5 ск
Рис. 5. Расщепление разности между рассчитанным и заданным распределением давления на дозвуковую и сверхзвуковую части
Сверхзвуковая часть 5сзвср используется для получения модификации поверхностных наклонов в сверхзвуковой зоне по формуле Аккеретта:
5y' = Ü.WM2 - 15^.
Затем поправки поверхностных наклонов в сверхзвуковой зоне преобразуют в эквивалентен экв тт
ное дозвуковое распределение поправок давления одзвср . Для этого проводят два последовательных расчета дозвукового обтекания профиля панельным методом: один раз с исходными поверхностными наклонами, второй раз — с модифицированными. Общая поправка давления в несжимаемой жидкости определяется как сумма двух составляющих
Ъскр = Г§дзвСр + ^дзв^Г, Г ^ 1 Ч ^ 1
Далее процедура аналогична предыдущей, но требует дополнительного усиленного сглаживания геометрии профиля в районе скачка уплотнения.
3. ПРИМЕРЫ РАСЧЕТОВ
Для проверки обратных методов часто применяют тесты по восстановлению геометрии. Они заключаются в том, что заданное распределение давления берется с прямого расчета существующего профиля, а исходная геометрия выбирается произвольно. Такие тесты позволяют контролировать сходимость метода не только по распределению давления, но и по геометрии. Заметим, что в общем случае решения обратной задачи искомая геометрия неизвестна, и, более того, в общем случае не существует доказательства ее существования, т. е. обратная задача часто является некорректной.
На рис. 6 представлен пример решения обратной задачи для заданного распределения значений коэффициента давления, полученного при расчете дозвукового обтекания профиля NACA 4412 на режиме М = Ü.2, а = 6°. Рассматривается полностью турбулентное течение при числе Рей-
0 12 3 4 5
Рис. 6. Решение обратной задачи для дозвукового обтекания M ш = 0.2, а = 6°
нольдса Re = 4 -1Ü6. В качестве начального приближения был взят симметричный профиль NACA ÜÜ1Ü под углом атаки а = 3°. На рисунке представлены результаты, полученные после пяти глобальных итераций (каждая итерация включает прямой расчет и один шаг решения обратной задачи). Построенный контур профиля, угол атаки и рассчитанное по ним распределение давления практически идентичны задаваемым. Так как корректор хорошо описывает физику течения, то достаточно хороший результат получается уже на первой итерации. На рисунке также приведен график уменьшения невязки по итерациям и график сходимости по углу атаки.
В следующем примере (рис. 7) обратная задача решалась для заданного распределения давления, полученного при расчете обтекания профиля NACA 4412 на режиме М = 0.67, а = 2°. На этом режиме у данного профиля наблюдается скачок уплотнения и отрыв пограничного слоя на верхней поверхности. Начальное распределение давления соответствовало профилю NACA 0010, обтекаемому под углом атаки а =1°. Для решения задачи проектирования потребовалось двенадцать итераций. Достаточно хорошая сходимость распределения давления получена на седьмой итерации, но было проведено еще пять итераций, чтобы получить лучший результат для геометрии. На рисунке также приведен график уменьшения невязки по итерациям и график сходимости значения угла атаки.
Рис.
0 1 234 56 789 10 11 12
7. Решение обратной задачи при трансзвуковом обтекании Мш = 0.67, а = 2°
Рис. 8. Решение обратной задачи для целевого распределения давления, заданного произвольно
В третьем примере целевое распределение давления задано произвольно (рис. 8) с целью реализации бесскачкового обтекания на верхней поверхности профиля на режиме М = 0.67. Распределение давления на нижней поверхности не менялось. В качестве начального приближения выбран профиль NACA 4412. Решение задачи с приемлемой точностью получено за десять итераций.
ЗАКЛЮЧЕНИЕ
Разработан алгоритм и реализована программа решения обратной задачи для профиля при помощи решения уравнений Навье — Стокса программным пакетом АК8У8 СБХ. Решение осуществляется итерационно по методу остаточной коррекции. В качестве корректора использован разработанный ранее в ЦАГИ метод решения обратной задачи для профиля в невязком сжимаемом потенциальном потоке — ТИАТКУ.
Приведены примеры решения обратных задач для профиля при дозвуковом и трансзвуковом обтекании, демонстрирующие высокую эффективность метода. В дальнейшем он может быть развит для решения обратных задач более сложного вида, в том числе двумерных задач для многозвенного профиля, и трехмерных задач.
ЛИТЕРАТУРА
1. Елизаров А. М., Ильинский Н. Б., Поташев А. В. Обратные краевые задачи аэро-гидродинамики. — М.: Физматлит, 1994.
2. V o l p e G. Inverse design of airfoil contours: constraints, numerical method and applications // AGARD-CP-463, 1989, p. 4.1—4.18.
3. N i k o l s k i А. А. Some aspects of helicopter airfoil design // Twenty-first European Ro-torcraft Forum, Aug. 31. — Sept. 1, 1995, St-Petersburg, paper N 17.
4. Davis W. H. Technique for developing design tools from the analysis methods of computational aerodynamics // AIAA Paper 79-1529, 1979.
5. Jimenes-Varona J. An inverse method for transonic wing design // Intern. J. of Numerical Methods in Engineering. 1999. V. 44, N 2, p. 249 —264.
6. Hirose N., Takanashi S., Kawai N. Transonic airfoil design procedure utilizing a Navier — Stokes analysis code // AIAA J. 1987. V. 25, N 3, p. 353—359.
7. Fray J. M. J., S l o o f J. W., Boerstol J. W., Kassies A. Inverse method with geometric constraints for transonic aerofoil design // Intern. J. for Numerical Methods in Engineering. 1986. V. 22, N 2, p. 327—339.
8. Labrujere Th. E. Residual-correction type and related computational method for aerodynamic design / VKI Lecture Notes on Optimum Design Methods in Aerodynamics. — Brussels, April, 1994.
9. Malone J. B., Narramore J. C., Sankar L. N. An efficient airfoil design method using the Navier — Stokes equations // AGARD CP-463, 1989, p. 5.1—5.18.
10. Campbell R. L. Efficient viscous design of realistic aircraft configurations // AIAA Paper 98-2539, 1998.
11. Dulikravich G. S. Aerodynamic shape design and optimization: status and trends // J. of Aircraft. 1992. V. 29, N 6, p. 1020—1026.
12. Болсуновский А. Л., Бузоверя Н. П. Метод решения обратной задачи для профиля в околозвуковом потоке газа // Труды ЦАГИ. 1993, вып. 2497, с. 17—27.
13. Болсуновский А. Л., Бузоверя Н. П., Карась О. В., Ковалев В. Е. Развитие методов аэродинамического проектирования крейсерской компоновки дозвуковых самолетов // Труды ЦАГИ. 2002, вып. 2655, с. 133 — 145.
14. V o s J. B., R i z z i A., D a r r a c q D., H i r s c h e l E. H. Navier — Stokes solvers in European aircraft design // Progress in Aerospace Sciences. 2002. V. 38, p. 601 —697.
15. Волков А. В. Методы решения сеточных уравнений конечно-элементной аппроксимации пространственных течений // Ученые записки ЦАГИ. 2010. Т. XLI, № 3, с. 52—68.
16. Вышинский В. В., Судаков Г. Г. Применение численных методов в задачах аэродинамического проектирования // Труды ЦАГИ. 2007, вып. 2673, 141 с.
17. Босняков С. М., Акинфиев В. О., Власенко В. В., Глазков С. А., Горбушин А. Р., Кажан Е. В., Курсаков И. А., Лысенков А. В., Матяш С. В., Михайлов С. В. Использование методов вычислительной аэродинамики в экспериментальных работах ЦАГИ // Матем. моделирование, 23:11(2011), с. 65—98.
18. Матяш С. В. Новый метод постановки граничных условий на удаленной границе в конечно-объемных методах численного решения аэродинамических задач // Ученые записки ЦАГИ. 2009. Т. XL, № 1, с. 42—51.
19. Головина Н. В. Сравнение результатов численных расчетов методом, основанным на разностной схеме Годунова — Колгана — Родионова, с экспериментальными данными для случая трансзвукового обтекания профиля RAE 2822 // Ученые записки ЦАГИ. 2009. Т. XL, № 5, с. 41—52.
20. www.ansys.com. Лицензия ANSYS Academic Research № 607044.
Рукопись поступила 16/V 2012 г.
