Solving of the inverse problem for a multielement airfoil in a compressible viscous gas flow

Full Text

1. Введение

Многозвенные или разрезные профили обладают рядом особенностей обтекания, отсутствующих у изолированных профилей. Скорость на задней кромке впереди расположенных элементов может быть заметно выше, чем скорость набегающего потока (Cp < 0) (рис. 1, [1]), в связи с чем перепад давления, который должен преодолевать каждый пограничный слой, снижается. Весь перепад давления от передней до задней кромки разрезного профиля преодолевается не одним пограничным слоем, как в случае изолированного профиля, а «эстафетой» свежих пограничных слоев на каждом элементе.

 

Рис. 1. Расчетное и экспериментальное [1] распределения давления на многозвенном профиле, М = 0.2, Re = 9×10⁶, α = 16°

 

Этот благоприятный для повышения несущих свойств разрезных профилей факт был отмечен еще в довоенных работах В.В. Голубева [2, 3]. Дальнейшие объяснения физики обтекания механизированных профилей и крыльев при малых скоростях были даны значительно позже [4, 5].

При больших скоростях в сжимаемом газе на верхней поверхности изолированного профиля развивается местная сверхзвуковая зона (М > 1), заканчивающаяся скачком уплотнения, и далее дозвуковой поток тормозится к задней кромке. На многозвенных профилях, ввиду меньших перепадов давления, можно обеспечить практически бесскачковое обтекание (рис. 2, [6]), получив тем самым выигрыш в подъемной силе или в продвижении по числу Маха.

 

Рис. 2. Концепция щелевого крыла NASA [6]

 

Для построения формы разрезного профиля в сжимаемом газе по заданному распределению давления необходим метод решения обратной задачи. Ранее авторами был разработан соответствующий метод для изолированного профиля в трансзвуковом потоке [7] и для многозвенного профиля при малых скоростях [8]. В данной работе дано описание итерационного метода решения обратной/смешанной краевой задачи аэродинамики для многозвенного профиля при больших скоростях в рамках осредненных уравнений Рейнольдса (RANS-методы). Он является развитием предыдущих методов.

Дано краткое описание алгоритма и используемых методов. Приведены примеры построения геометрии элементов двухзвенного профиля по заданному распределению давления, в том числе и при наличии скачков уплотнения.

2. Проектирование многозвенного профиля

Надежные прямые методы CFD, особенно с использованием RANS-уравнений, позволили существенно снизить объем экспериментальных исследований и глубже проникнуть в физику явлений, сопровождающих обтекание разрезных крыльев. Однако для эффективности процедуры аэродинамического проектирования одних лишь прямых методов недостаточно, необходимо также иметь в арсенале оптимизационные и обратные методы [8–12]. Обратные методы позволяют построить геометрию многоэлементного профиля по заданному распределению давления Cp^* на поверхности. Кроме обратных задач, на практике часто решаются и смешанные задачи, в которых часть элементов фиксирована, и на них необходимо найти распределение давления, а геометрия других должна быть определена по заданному целевому Cp^*.

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

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

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

 

Рис. 3. Принцип построения методов остаточной коррекции

 

В данной работе метод остаточной коррекции используется для построения геометрии многоэлементного профиля по заданному распределению давления в сжимаемом потоке вязкого газа при больших скоростях. В качестве прямого метода выбран расчет обтекания многозвенного профиля в рамках уравнений RANS. Блок коррекции геометрии использует панельные методы решения прямой/обратной задачи для механизированного профиля в идеальной несжимаемой жидкости (рис. 4). Ранее авторы использовали метод остаточной коррекции при построении изолированных профилей в трансзвуковом потоке вязкого газа [7], при построении механизированных профилей при малых скоростях [8], а также применительно к решению пространственных обратных задач для крыльев магистральных самолетов при трансзвуковых скоростях [12].

 

Рис. 4. Структура используемого метода

 

Заметим, что любой обратный метод, так или иначе, связан с конкретным прямым методом. Расчетные распределения давления, полученные различными прямыми методами, могут немного отличаться как друг от друга, так и от экспериментальных значений, и поэтому вопрос о “точности” обратного метода некорректен, так как она непосредственно зависит от точности используемого прямого метода.

3. Постановка задачи и алгоритм расчета

Решается обратная задача нахождения геометрии элементов многозвенного профиля по заданному распределению давления в сжимаемом потоке вязкого газа. Возможна и смешанная постановка задачи, когда ряд элементов выстраивается по заданному распределению давления, а остальные элементы не изменяются. Набор элементов, геометрия которых будет изменяться, задается в исходных данных вместе с координатами исходных узловых точек контура каждого j-го элемента многозвенного профиля (x_ij,y_ij) и желаемым распределением давления на варьируемых элементах ($C{p}_{ij}^{*}$). Для корректной постановки задачи необходимо также задать точки элементов, которые должны оставаться неподвижными при изменении геометрии (например, передняя или задняя кромка элемента).

Обратная задача решается методом итераций. На каждой k-й итерации выполняются следующие действия.

1. Строится сетка вокруг многозвенного профиля, геометрия которого получена на предыдущей (k–1)-й итерации; рассчитывается его обтекание при помощи RANS-метода; определяется невязка распределения давления $\mathrm{\delta }C{p}^k$ как разница между текущим и заданным распределениями давления для вязкого обтекания.
2. Панельным методом рассчитывается обтекание того же профиля несжимаемой жидкостью. При этом отрывные «пузыри» в вырезах предкрылка и основного профиля, полученные по RANS-методу, необходимо “добавить” к исходной геометрии для моделирования обтекания «жидкого» контура.
3. а. В случае докритического обтекания целевое распределение давления для профиля в несжимаемой жидкости формируется путем прибавления разницы ${\mathrm{\delta }}_{0}C{p}^k=r\mathrm{\delta }C{p}^k$ (где коэффициент демпфирования r ≤ 1) к распределению давления, найденному на предыдущем шаге.
4. б. При наличии в области течения сверхзвуковых зон разность между рассчитанным и заданным распределением давления на каждом варьируемом элементе расщепляется на дозвуковую и сверхзвуковую части (рис. 5):

$\mathrm{\delta }C{p}^k={\mathrm{\delta }}_{-}C{p}^k+{\mathrm{\delta }}_{+}C{p}^k$

 

 

Рис. 5. Расщепление невязки давления на дозвуковую и сверхзвуковую части.

 

Сверхзвуковая часть ${\mathrm{\delta }}_{+}C{p}^k$ используется для получения модификации поверхностных наклонов в сверхзвуковой зоне по формуле Аккерета

$\mathrm{\delta }{y}^{\text{'}}=\frac{1}{2}\sqrt{M^2-1}{\mathrm{\delta }}_{+}Cp$

Затем поправки поверхностных наклонов в сверхзвуковой зоне преобразуют в эквивалентное дозвуковое распределение поправок давления ${\mathrm{\delta }}_{-}C{p}^{\mathrm{equi}}$. Для этого проводят два последовательных расчета дозвукового обтекания профиля панельным методом: один раз с исходными поверхностными наклонами, второй раз – с модифицированными. Общая поправка давления в несжимаемой жидкости определяется как сумма двух составляющих

${\mathrm{\delta }}_{0}C{p}^k=r({\mathrm{\delta }}_{-}Cp+{\mathrm{\delta }}_{-}C{p}^{\mathrm{equi}})$; $r\le 1$

4. По целевому распределению давления с шага 3 решается обратная задача в несжимаемой жидкости и получается (k+1)-я геометрия профиля. Для большей устойчивости алгоритма предусмотрено сглаживание полученных координат вблизи носка контура и в районе скачка уплотнения.

Для полученной геометрии необходимо снова решить прямую задачу RANS, найти распределение давления и т.д. На каждой итерации контролируется невязка давления по элементам ${\epsilon }_{C{p}_j}=\sqrt{\frac{\text{1}}{N_j}\sum _i^{N_j}{\left(C{p}_{ij}-C{p}_{ij}^{{}_{*}}\right)}^{\text{2}}}$, определяющая отклонение текущего распределения давления от заданного. Процесс повторяется несколько раз до выхода на некоторую минимально достижимую невязку. Достаточная для практики сходимость соответствует уровню ε_Ср ≤ 0.02.

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

4. Используемые расчетные методы

В данной работе расчеты вязкого течения сжимаемого газа проводились на структурированных многоблочных сетках, позволяющих осуществить подробную дискретизацию области пограничного слоя и вязкого следа профиля. Число ячеек равно 150–300 тыс. в зависимости от количества элементов. Границы расчетной области удалены от профиля на 50 хорд. Пример структурированной расчетной сетки, использованной для решения представленных в работе задач, приведен на рис. 6. При изменении контуров элементов профиля в процессе итераций созданная начальная сетка автоматически перестраивается для соответствующей геометрии.

 

Рис. 6. Структурированная многоблочная сетка

 

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

Для прямого/обратного расчета обтекания механизированного профиля несжимаемой жидкостью применяется панельный метод с кусочно-линейным распределением плотности вихревого слоя на плоских панелях (рис. 7). Для замкнутого контура при выполнении условия непротекания касательная скорость равна интенсивности вихревого слоя. Для выполнения условия Жуковского, интенсивности вихревого слоя на задней кромке равны и противоположны по знаку. В случае разомкнутой задней кромки, на ней размещается панель с источниками, интенсивность которых равна интенсивностям вихревого слоя [9].

 

Рис. 7. Вихревые особенности на панелях многозвенного контура

 

В основу расчета положено постоянство функции тока на контуре каждого элемента, из которого следует интегральное уравнение, связывающее координаты узловых точек с плотностью вихревого слоя (для упрощения понимания представлено уравнение для одного элемента):

${\mathrm{\psi }}_i=y_i\text{cos}\alpha -x_i\text{sin}\alpha -\frac{\text{1}}{\text{2}\mathrm{\pi }}\oint \gamma \left(s\right)\text{ln}\sqrt{{(x_i-x(s\left)\right)}^{\text{2}}+{(y_i-y(s\left)\right)}^{\text{2}}}ds=\mathrm{const}$

$i=\mathrm{1,}\mathrm{2,}\dots ,N+1$

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

$y_i=\frac{\text{1}}{\text{cos}\alpha }\left[{\mathrm{\psi }}_c+\frac{\text{1}}{\text{2}\pi }\underset{C}{\oint }\gamma \text{(s)ln}\sqrt{{(x_i-x(s\left)\right)}^{\text{2}}+{(y_i-y(s\left)\right)}^{\text{2}}}ds\right]+x_i\mathrm{tg}\alpha$

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

Последнее уравнение нелинейное, так как неизвестные значения у_i входят под знак интеграла. Решение ищется в процессе итераций, в результате определяется новая форма контура. Значение ${\mathrm{\psi }}_c$ выбирается в зависимости от того, какая точка на контуре является фиксированной. Это может быть, например, передняя или задняя кромка. Распределение интенсивностей вихревого слоя является заданным для варьируемых элементов, а для фиксированных элементов его необходимо периодически пересчитывать при помощи прямого метода. В данной работе пересчет циркуляций проводился после каждого шага решения обратной задачи.

В целом, расчетные методы аналогичны тем, которые используются при решении обратной задачи при малых скоростях [8]. Единственным дополнительным методом, который применяется для расчета эквивалентной дозвуковой поправки давления ${\mathrm{\delta }}_{\mathrm{ДЗВ}}C{p}_{}^{\mathrm{ЭКВ}}$, является панельный метод с тем же распределением гидродинамических особенностей, но с явным выполнением условия непротекания V_N=0 в центрах панелей.

5. Примеры расчетов

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

В качестве начальной геометрии был выбран разрезной профиль с относительной толщиной $\overline{c} = 22.3\%$, разработанный авторами ранее [8] для малых скоростей и испытанный в нескольких аэродинамических трубах. Расчеты проводились при числе Рейнольдса Re = 5×10⁶ и полностью турбулентном обтекании. При числе Маха М = 0.5 и угле атаки α = 5° на исходном профиле появляется скачок уплотнения (рис. 8). Целевое распределение давления задавалось так, чтобы убрать этот скачок, сохраняя при этом подъемную силу и относительную толщину основного звена. Форма закрылка не менялась. Итерационный процесс сошелся за 12 итераций. Расчетное сопротивление упало на 37 каунтов (один каунт для C_x равен 0.0001) за счет снижения интенсивности скачка.

 

Рис. 8. Пример решения обратной задачи при М = 0.5, ∆С_у, ∆С_x – изменения интегральных характеристик по сравнению с характеристиками исходного профиля

 

Полученный при М = 0.5 профиль был выбран во втором примере в качестве начального, но уже при М = 0.6 и α = 2°. При этих условиях на верхней поверхности возникает сильный скачок уплотнения, вызывающий отрыв пограничного слоя (рис. 9). Заданное распределение давления носило полочный характер, часть подъемной силы была дополнительно реализована за счет модификации нижней поверхности. Отметим, что заданное распределение давления несколько раз модифицировалось в процессе проектирования на основании получаемых промежуточных результатов. Закрылок снова был фиксирован.

 

Рис. 9. Пример решения обратной задачи при М = 0.6

 

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

 

Рис. 10. Сходимость решения обратной задачи при М = 0.6 

 

Итерационный процесс сошелся за 14 итераций, относительная толщина снизилась до $\overline{c} = 20.8\%$. Расчетное сопротивление упало на 260 каунтов за счет снижения интенсивности скачка и ликвидации отрыва. Получено безотрывное обтекание при C_y = 1.31 – такой уровень несущих свойств недостижим на бесщелевом профиле.

При помощи разработанного метода был спроектирован разрезной профиль, предназначенный для более высоких скоростей (М = 0.768), представленный в третьем примере. На рис. 11 приведены геометрия профиля и распределение числа Маха в пространстве, а на рис. 12 – распределение давления по поверхности профиля на режиме М = 0.768, Re = 10×10⁶, C_y = 0.75. На рис. 13 приведена поляра С_x(C_y) для полученного профиля в сравнении с полярой для бесщелевого профиля, также спроектированного под данный режим. При малых C_y сопротивление щелевого профиля больше на 10–20 каунтов за счет наличия щели. Его преимущество реализуется при более высоких C_y. Он позволяет достигать C_y = 0.75, в то время, как на изолированном профиле резкий рост C_x начинается уже с C_y = 0.55.

 

Рис. 11. Поле числа Маха на режиме М = 0.768, С_у = 0.75

 

Рис. 12. Распределение давления на режиме М = 0.768, С_у = 0.75

 

Рис. 13. Расчетные поляры разрезного и изолированного профилей

 

Заключение

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

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

About the authors

A. L. Bolsunovsky

Central Aerohydrodynamic Institute

Author for correspondence.
Email: bolsmail@mail.ru
Russian Federation, Zhukovsky

N. P. Busoverya

Central Aerohydrodynamic Institute

Email: bolsmail@mail.ru
Russian Federation, Zhukovsky

S. V. Gerasimov

Central Aerohydrodynamic Institute

Email: bolsmail@mail.ru
Russian Federation, Zhukovsky

M. A. Gubanova

Central Aerohydrodynamic Institute

Email: bolsmail@mail.ru
Russian Federation, Zhukovsky

References

Supplementary files