CC BY
43
УДК 621.452
ИДЕНТИФИКАЦИЯ ВЫЧИСЛИТЕЛЬНОГО АЛГОРИТМА ДЛЯ РАСЧЕТА ОБТЕКАНИЯ ТРАНСЗВУКОВОЙ КОМПРЕССОРНОЙ РЕШЕТКИ
Д.Л. БУТРИМОВ, К.С.ФЕДЕЧКИН Статья представлена доктором технических наук, профессором Умушкиным Б.П.
В статье описывается вычислительная методика решения уравнений Эйлера, рассматриваются процедуры, позволяющие ускорить получение установившегося решения.
Ключевые слова: вычислительная методика, решения уравнений Эйлера, обтекание трансзвуковой решетки, задача оптимизации.
Введение
Математическое моделирование в настоящее время широко применяется для описания процессов и явлений в различных технических системах и их элементах. Для моделирования течений жидкостей и газов наиболее универсальными являются модели, базирующиеся на уравнениях Навье-Стокса.
Существующие методики решения уравнений газовой динамики в ряде случаев позволяют успешно решать практические задачи применительно к определенным объектам исследования. Однако для получения достоверных результатов необходимо идентифицировать математическую модель применительно к конкретному объекту исследования. В настоящей статье рассматривается возможность идентификации вычислительного метода на основе выбора оптимальных значений коэффициентов применительно к моделированию течений в компрессорной решетке профилей.
Математическая модель течения жидкости и газа
Для моделирования течения газа используется система дифференциальных уравнений Эйлера
Щ-+уР=°, (1)
где вектор консервативных переменных и конвективных потоков равен
' Р г ( Г? Л Р^
и = ри Р = рuV + рі
рv 5 рvV + р]
V РЕ , V рн^? )
(2)
где і - время; р, р, и, V - давление, плотность, декартовы компоненты скорости; Е - полная энергия; Н - полная энтальпия; і и ] - единичные вектора, связанные с декартовой системой
координат; V = иі + V] - вектор скорости.
Объект
исследования
Ф
Г еометричес-кая модель
с>
Расчетная
сетка
Эксперимент
ж
YЯI=
А® °
(к1,.., кг),
К '
Математическая
модель
^+vF=0,
F=F(V,p, p, к, к,,..., к
кі, k2,.., К
Итерационная процедура
Рис. 1
Для замыкания системы уравнений используется уравнение состояния для идеального невязкого нетеплопроводного газа
Р = (к - 1)Р
Е-
2 2 Л
и + V '
2
(3)
где к - показатель адиабаты.
Для получения численного решения исходную систему уравнения преобразуют в дискретную форму. Дискретизация системы осуществляется по пространству и по времени. После проведения процедуры дискретизации получается алгебраическая система уравнений, для решения которой существуют эффективные вычислительные алгоритмы и специальные процедуры ускорения сходимости (рис. 1).
Пространственная дискретизация
Интегральная форма уравнений Эйлера дискретизируется методом конечных объемов
(т т 4
\^ЛП +1 [Р •ЙА5]к =°, (4)
у к = 1
где - контрольный объем ячейки; Й - вектор нормали к грани контрольного объема; Д5 -площадь (в двумерном случае - длина) к-й грани; к - грань контрольного объема.
Рассчитываемые параметры течения относятся к центрам контрольных объемов (рис. 2). Конечно-объемные методы подразумевают преобразования конвективных потоков и допускают их постоянство вдоль границ.
Конвективный поток рассчитывается по центрально-разностной схеме
' и у + и ц+1Л 2
і >] +,
(5)
г
У
х
• 1+І0 М 4+1/2 / ,
V 1
£ ?и 4 § !<Ы , Ь ?и+1
&> |
к
Рис. 2
Для подавления нелинейной неустойчивости, которая может проявляться на разрывах и в областях с большим градиентом параметров потока, при расчете конвективных потоков в дискретную систему уравнений добавляется дополнительный член, называемый искусственной вязкостью, в форме предложенной Джеймсоном [1]
где
а
и+у +Уг5и<-1+/2 + +уг8іи>-і+1 ’
зи,, і+% = и,і і+, -
е3и,, + = и
І, і +
і , і + 2
3иг,і+! + 3ии -ииі-!.
, і
Є,і + X = 0,5 к 1 І, і + 1/2таХ(и, і-1,иі, і ,иі, і +, і + 2 )
£1 + К = -тах(°;°,5к(4)Лщи+1/2 -еі*. +у)
(6)
(7)
Л І,і + 1/2 = (
+ с А 5) І, і
і +1/2
(8)
Здесь переменная Ц у. = ру.+1 — 2ру. + рг.,у—1 /( рг.,у+1 + 2ру. + рг.,у—1) является сенсором,
активирующим второго порядка разностную диссипацию в области сильных градиентов давления, и деактивируется в остальных областях, к2-'1 и - константы, значения которых
л *
варьируются соответственно от 0,25 до 1,0 и от 0,015625 до 0,1, 1 i j+^2 - спектральный радиус
матрицы [ЭЕ/Эи ] А8,,] +^2 конвективных потоков (с - локальная скорость звука).
В большинстве современных методик [2],[3],[4] предлагается учитывать для ячеек, у которых соотношение смежных сторон значительно отличается от единицы, величину спектрального радиуса по всем координатным направлениям, например, через следующее соотношение
Л,,і+1,2 = тах(Л і,і+1,2,(Л\у+1/2)(1-")Л,Г).
(9)
2
л *
Спектральный радиус 1,+^2 ^ определяется таким же образом, что и в соотношении (8). В
работе [5] отмечается повышение устойчивости вычислительного алгоритма и скорости сходимости к установившемуся решению применительно к обтеканию крыльевого профиля для а = 0,5. Параметр а позволяет изменять величину искусственной диссипации по координатным направлениям и варьируется в диапазоне от 0 до 1.
Стандартные значения к2, к4 и а принимаются
к(2) = 1,0, к(4) = 0,1, а = 0,5. (10)
Дискретизация по времени
Для дискретизации по времени система уравнений (4) с учетом соотношения (5) преобразуется к виду
(и пі - и п,1)
Гп+1
’} ~ І’і ' — А I о/цпч ПЛ
—[ я (ип) - я: ] (11)
А* пІ, Л ' ' :
где п — индекс временного слоя; Я(Оп) — среднеобъемная невязка конвективных потоков
через грани контрольного объема; Яа — численная диссипация в ячейке (І,і).
Полученная дискретная система уравнений решается методом Рунге-Кутта по гибридной схеме [6]
и ° = Оп
и' = и"-“1‘АГ[Я^°) - Я:]
О2 = Оп-аг А- [ Я (і;1) - Я: ]
и3 = Оп-а, А-[Я(1}2)-Я'їт] (12)
о 4 = Оп-а4 — [ Я (О3) - Я'2т ]
>(2,°)
- I Я (О ) — Я
п
А г
~п
А г
~п
ип+1 = о5,
О! = Оп-а, п [ Я (О4) - ЯГ> ]
Я'/’ = X ¿к- (и “) (13)
к
Я <2-0> =Ьз ЯГ + (1 _ Ьз) ЯГ я = Ь Я Г + (1 _ Ь) Я Г0).
В таблице приведены значения коэффициентов а и в, значения которых оптимизированы для применения данной схемы совместно с многосеточным методом.
Таблица 1
№ этапа а в
1 0.2500 1.00
2 0.1667 0.00
3 0.3750 0.56
4 0.5000 0.00
5 1.0000 0.44
Как видно из таблицы, численная диссипация может не рассчитываться на втором и четвертом этапе гибридной схемы.
Процедуры ускорения сходимости решения
1. Локальный шаг интегрирования по времени
При численном интегрировании дискретной системы уравнений по времени используется максимально возможный по условию устойчивости шаг интегрирования для каждой ячейки
СРЬ _
D , = —------------- W , ,
г’7 Л' + Л ] ' ,7
(14)
где С¥Ь - число Куранта-Фридрихса-Леви; X и X - определяются подобно выражению (8).
Следует заметить, что строго математического обоснования данного подхода нет, однако на практике вычислительный алгоритм получается более устойчивым и повышается скорость сходимости решения. Очевидно, что для нестационарных задач локальный шаг интегрирования по времени неприменим.
2. Процедура неявного сглаживания невязок
Как известно, явные схемы имеют жесткие ограничения на максимально допустимый шаг интегрирования по времени. Данные ограничения обобщаются комплексным критерием, называемым числом Куранта-Фридрихса-Леви (СБЬ). Для описанной выше гибридной схемы максимальное значение СБЬ, при котором возможен устойчивый счет, равно 3,6. Для смягчения значения СБЬ, а следовательно, увеличения шага интегрирования по времени, разработана процедура неявного сглаживания невязок, суть которой заключается в замене невязок в каждой ячейке комбинацией осредненных невязок по нескольким ячейкам. Существуют различные варианты определения сглаживающих коэффициентов для комбинирования невязок [2],[7],[8]. В данной работе используется методика, предложенная в работе [2]
є 1Я *І_х,3 + (1 + 2 є 1 ) Я * є 3 Я **І,3 _ 1 + (1 + 2 є 3 ) Я
єІЯ і += ЯІ,
(15)
** *
І , 3 + 1 “ Я І , 3
Здесь Я ,,, Я
,,, я означают сглаженные невязки в направлениях і и , соответственно.
к,,=[ Я(0") _ Я: ] — несглаженная невязка.
Система уравнений (15) можно представить в виде матрицы, которая имеет трехдиагональную структуру. Уравнения в этом виде эффективно решаются методом прогонки [9].
Сглаживающие коэффициенты определяются следующим выражением
є
є
тах
тах
С¥Ь*
Л1
С¥Ь Л1 + вЛ
СП*
ЛЯ
С¥Ь вЛ + Л^
_ 1
_ 1
0
(16)
2
2
0
в = 0,125. СБЬ и СБЬ сглаженной и несглаженной схемы такое, что
С^<^747. <17>
СРЬ
Из опыта практического применения вычислительного алгоритма установлено, что оптимальное соотношение чисел СБЬ равно 2. Таким образом, «эффективное» число СБЬ увеличивается в два раза.
3. Многосеточный метод
Существенно повысить скорость сходимости решения позволяет геометрический многосеточный метод. Он основывается на решении дискретных уравнений на различных сеточных уровнях. В процессе решения задачи на исходной подробной сетке эффективно подавляются высокочастотные компоненты ошибки, низкочастотные демпфируются слабо. С увеличением размера ячеек низкочастотные компоненты становятся высокочастотными, их можно быстро подавить стандартными вычислительными процедурами, которые применяются на подробной сетке. Детально многосеточная технология описана в работах [2],[6],[10],[11].
Отметим основные этапы применяемого метода.
Для простоты рассмотрим две сетки. Индексом И обозначается подробная сетка, 2И - грубая, получаемая удалением каждой второй линии сетки в каждом координатном направлении.
А. Расчет решения на подробной сетке.
После одного шага интегрирования по времени дискретных уравнений (11) по схеме (12) получаем решение на следующем временном слое
и 5
• (18)
1.}
При этом для дальнейших вычислений сохраняются невязки. полученные на 5-ом этапе схемы Рунге-Кутта
Я*(0Ґ) = [Я(И4) - Я'4-2']* (19)
Б. Перенос значений переменных на грубую сетку.
Значения переменных на грубой сетке получается весовым осреднением по контрольному объему
£ [п *и» ”*■"
и " ' '
* * J к
£ [п .и
0 _ к _ 1_______________
2 * _ £[п* ] ■ (20)
к _ 1
Здесь производится суммирование по четырем ячейкам подробной сетки. которые составляют одну ячейку грубой сетки. После переноса решения на грубую сетку осуществляется реализация первого этапа гибридной схемы. на котором находится невязка
Я 2 * (и 2*) _ Я (И0) - Я 0.
В. Ограничение невязок на грубой сетке.
Для сглаживания низкочастотных компонент ошибки осуществляется перенос невязок с подробной сетки на грубую. Оператор переноса невязок такой же. как и при переносе значений переменных на грубую сетку
£ [п*й „(ИГ1)'
і;пк- * (и; *‘) _
£ [п * і
0, к _1
т2п
где - оператор ограничения.
Г. Получение решения на грубой сетке.
Решение на грубой сетке получается по схеме (12) с добавлением дополнительного члена, называемого силовой функцией
г~- — - - (22)
(й, )2* = [^ *<°Г1) - *2*(02*)
В итоге схему (12) можно записать в виде
А і
2 *
ик = и - а
2 * 2 * мк
0п+1 = и5
2* 2*
2*
о
2*
К (и1-1) - +(о, )2 *
к = 1...5
(23)
Как видно из соотношения (23), на первом этапе схемы решение на грубой сетке полностью формируется из невязок, определенных на подробной сетке. На последующих этапах силовая функция также «управляет» решением, получаемым на грубой сетке. Это гарантирует соответствие «точного» решения, получаемого на подробной сетке, и «приближенного», получаемого на грубой. Расчет параметров потока на грубой сетке производился с добавлением численной диссипации второго порядка с постоянными коэффициентами.
Для следующих уровней грубой сетки, например, (4Ь), силовая функция формируется следующим образом
( й> )4 *
14^2 * (иг* ‘) - К 4 * (и 4* )
4*
4Ь
2Ь
[(К(6«) - К>к X* + (й, )
К4* (и4*)
(24)
4*
'2к /2к .
Процесс повторяется до достижения заданного уровня грубой сетки.
Д. Расчет коррекции решения на грубой сетке.
После одной или нескольких итераций, осуществляемых по схеме (23), вычисляется коррекция решения
• п + 1 т— 0
(25)
Поправка к решению А и 2 к рассчитывается в каждой ячейке и далее ее значения
интерполируются на подробную сетку.
робной сетке.
Е. Расчет решения на под Окончательное решение
и,
і, і
на подробной сетке получается на основе решения (18) и
поправки (25), переносимой с грубой сетки
"и* " = - ип+1 - +
1, і 1 , і
2 *
(26)
где 12* - билинейный оператор, называемый оператором продолжения.
Алгоритм оптимизации
Как отмечалось выше, целью проведенных исследований являлось изучение возможности повышения адекватности и точности вычислительного метода на основе выбора оптимальных значений коэффициентов применительно к моделированию течений в компрессорной решетке профилей. Для определения значений констант вычислительной методики решается задача оптимизации.
Для решения оптимизационной задачи используется метод непрямой оптимизации на основе самоорганизации, реализованный в программном комплексе ІОБО [12]. Данный метод основывается на процедурах аппроксимации функций многих переменных. На начальной
I
стадии процесса оптимизации точность функций аппроксимации может быть весьма незначительной из-за небольшого количества точек в плане эксперимента и сравнительно большой текущей области поиска. Однако по мере решения задачи количество точек в плане эксперимента, находящихся в окрестности искомых точек экстремумов, возрастает. В то же время уменьшается размер текущей области поиска. Эти тенденции приводят к увеличению точности функции аппроксимации и, следовательно, к повышению эффективности работы алгоритма оптимизации. Фактически в процессе оптимизации происходит постоянное накопление информации об исследуемом объекте и использование этой информации для определения направления дальнейшего поиска. Важнейшим достоинством используемого алгоритма оптимизации является его способность решать практические задачи для случая невыпуклых, недифференцируемых и стохастических целевых функций и ограничений.
Следует особо отметить, что в процессе решения задачи оптимизации количество обращений к математической модели, для которой решается задача оптимизации, может достигать нескольких сотен или даже нескольких тысяч раз, поэтому применение процедур, позволяющих повысить скорость сходимости решений, является одним из ключевых элементов вычислительной методики. Так, например, включение всех описанных выше методов, ускоряющих сходимость решения, позволило сократить время вычислений в 2.. .3 раза.
Численные исследования
В качестве объекта исследования была взята решетка профилей, соответствующая цилиндрическому сечению лопаток рабочего колеса трансзвукового компрессора. Детальное описание геометрии и результаты экспериментальных исследований представлены в статье [13]. В представленной решетке профилей рассматривалось течение идеального невязкого нетеплопроводного газа.
Расчеты проводились на структурированной сетке типа Н размерности 149х49 (75 ячеек по хорде) (рис. 3).
Граничные условия
На входной границе задавались полное давление, полная температура и угол атаки, на выходной - статическое давление. На профиле лопатки принималось условие непротекания для нормальной компоненты скорости. Периодические граничные условия выставлялись на
противоположных границах расчетной области, обладающих свойством повторяемости.
Рис. 3
Для расчетного угла атаки были получены зависимости коэффициента потерь полного давления в скачке уплотнения (£ск) и отношения статического давления на выходе из решетки к статическому давлению на входе от числа Маха на входе (р2/р1) (рис. 4) и (рис. 5) — соответственно кривая «Расчет ("стандартные" коэффициенты)». При этом для расчета были приняты значения коэффициентов, соответствующие «стандартным» (подобные значения используются в большинстве аналогичных методик расчета).
Видно, что результаты расчета отношения статического давления на выходе из решетки к статическому давлению на входе довольно хорошо согласуются с экспериментом (погрешность расчета не более 3%), а коэффициент потерь полного давления в скачке уплотнения существенно отличается от экспериментальных данных (в некоторых случаях более, чем в два раза). Отмечена следующая корреляция: в точках, где достигается минимальная погрешность для отношения статических давлений, наблюдается максимальная погрешность в определении потерь в скачке уплотнения и наоборот.
Задача оптимизации
Для повышения эффективности трансзвуковых лопаточных машин необходимо детально изучить природу потерь в них. Существенную часть потерь полного давления в трансзвуковых и сверхзвуковых решетках составляют потери, которые вызваны ударными волнами на входе в межлопаточный канал. Для рассматриваемой решетки потери в скачках уплотнения составляют 40.60% общих потерь в зависимости от угла атаки и числа Маха на входе. Точное определение потерь в скачках уплотнения является важным фактором для создания высокоэффективных турбомашин.
Для повышения точности расчета была решена задача оптимизации значений коэффициентов «искусственной» вязкости соотношения (7), т.к. именно этот член в разностных уравнениях отвечает за величину потерь.
В качестве варьируемых переменных были выбраны коэффициенты к2, к4) из соотношения (7) и а из соотношения (9). Коэффициент к2 варьировался в диапазоне от 0.25 до 5.0, к4 — от 0,015625 до 0,3, а — от 0 до 1. Критерием оптимизации выступала минимизация квадрата разности рассчитанных и экспериментально определенных потерь полного давления для заданного числа Маха на входе.
Для каждого заданного числа Маха потребовалось 100.150 обращений к математической модели, чтобы получить достаточное число решений для определения оптимальных соотношений параметров и величины целевой функции.
В результате анализа данных, полученных в процессе решения задачи оптимизации, установлено определяющее влияние параметра а на коэффициент потерь полного давления. С уменьшением а величина потерь, как правило, снижается. Рекомендуемое значение а = 0. Остальные параметры оказывали существенно меньшее влияние на величину потерь полного давления в рассматриваемом диапазоне варьирования, однако наилучшее совпадение с экспериментом было получено при следующих значениях коэффициентов: к2 = 2,0, к4 = 0,12. На рис. 4 и 5 кривые «Расчет ("оптимальные" коэффициенты)» представлены результаты расчета коэффициента потерь полного давления в скачке уплотнения и отношения статического давления на выходе из решетки к статическому давлению на входе с оптимизированным набором коэффициентов.
Заключение
Предложен подход, позволяющий идентифицировать вычислительный алгоритм для расчета обтекания трансзвуковой компрессорной решетки по результатам эксперимента путем решения задачи оптимизации. Анализ решения оптимизационной задачи позволил выявить параметры алгоритма и их значения, оказывающие определяющее влияние на целевую
функцию (критерий оптимизации) - коэффициент потерь полного давления. Показана возможность адаптации вычислительного алгоритма для повышения достоверности результатов расчета. Данный подход предполагается использовать для расчета вязких течений, основываясь на решении уравнений Навье-Стокса, осредненных по Рейнольдсу.
0.08
ск
0.07
0.06
0.04
0.03
0.8
0.85
0.9
0.95 Рис. 4
• 3 а счет ("с 3 а счет ("о Эксперт! гандартнь птимальн ент 1е" коэффициенты) ые" коэффициенты /1
1 ^
А •
V 1
1.05
1.1 М 1'1->
Рис. 5
ЛИТЕРАТУРА
1. Jameson A., Schmidt W., Turkel E. Numerical Solutions of the Euler Equations by Finite Volume Methods Using Runge-Kutta Time-Stepping Schemes // AIAA Paper 81-1259, 1981
2. Swanson R.C., Turkel E. Multistage Schemes with Multigrid for Euler and Navier-Stokes Equations // NASA Technical Paper N3631, 1997.
3. Kok J.C., Boerstoel J.W., Kassies A., Spekreijse S.P. A robust multi-block Navier-Stokes flow solver for industrial applications // Proceedings of the ECCOMAS Computational Fluid Dynamics Conference, September 1996. France, Paris, 1996.
4. Anand Kumar Role of Flow-consistent grid in CFD // Project Document CM 9805. India, Bangalore, 1998.
5. Swanson R.C. Turkel E. Artificial dissipation and central difference schemes for the euler and Navier-Stokes equations // NASA Contractor Report N178296, 1987.
6. Blazek J. Computational Fluid Dynamics: Principles and Applications. Amsterdam, London, New York, Oxford, Paris, Shannon, Tokyo: Elsevier, 2001.
7. Jameson A., Baker T.J. Solution of the Euler Equations for Complex Configurations // AIAA Paper 83-1929, 1983.
8. User Manual. FINE™/Turbo v8 (including Euranus) Flow Integrated Environment. - Brussels: Numeca Int., 2007.
9. Ворожцов Е.В. Разностные методы решения задач механики сплошных сред: учеб. пособие. -Новосибирск: Изд-во НГТУ, 1998.
10. Jameson A. Solution of the Euler Equations for Two-Dimensional Transonic Flow by a Multigrid Method // MAE Report No. 1613, Dept. of Mechanical and Aerospace Engineering, Princeton University, 1983.
11. Mavriplis D.J. Three-Dimensional Multigrid for the Euler Equations // AIAA Journal, 1992.
12. Егоров И.Н. Методы непрямой статистической оптимизации на основе самоорганизации и их использование в оптимизационных задачах авиационных ГТД. - М: ВИНИТИ, 1988.
13. Шрейбер Г.А., Штаркен Г. Исследование течения в элементарном венце трансзвукового компрессора методом испытаний решетки // Труды американского общества инженеров-механиков, серия Энергетические машины и установки, 1984. - № 2.
IDENTIFY A COMPUTING ALGORITHM FOR CALCULATION OF TRANSONIC FLOW AROUND COMPRESSOR LATTICE
Butrimov D.L., Fedechkin K.S.
This article describes a computational method for solving Euler equations are considered procedures to expedite obtaining the steady solutions.
Key words: computational method, Euler equations.
Сведения об авторах
Бутримов Дмитрий Леонидович, 1981 г.р., окончил ВВИА (2003), адъюнкт кафедры теории авиационных двигателей ВВИА им. проф. Н.Е. Жуковского, автор 7 научных работ, область научных интересов - математическое моделирование течения в осевых компрессорах газотурбинных двигателей.
Федечкин Константин Сергеевич, 1975 г.р., окончил ВВИА (1998),старший научный сотрудник, кандидат технических наук, начальник отделения научно - исследовательского отдела ВВИА им. проф. Н.Е. Жуковского, автор 30 научных работ, область научных интересов - математическое моделирование течения в осевых компрессорах газотурбинных двигателей, оптимизация элементов авиационных газотурбинных двигателей.
