Журнал вычислительной математики и математической физики, 2019, T. 59, № 12, стр. 2046-2059

1. ВВЕДЕНИЕ

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

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

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

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

Метод погруженных границ подразумевает моделирование условия на границе раздела двух сред (твердое тело и окружающая его газовая среда) внутри односвязной области при аппроксимации на сплошной, вообще говоря, неструктурированной сетке. В работе используется разновидность этого подхода, а именно метод Бринкмана штрафных функций [1]. Следуя используемому в данной работе методу Бринкмана штрафных функций, в систему уравнений газовой динамики добавляются штрафные функции, работающие только внутри обтекаемого препятствия. Метод Бринкмана штрафных функций успешно используется при численном моделировании задач внешнего обтекания, как в случае сжимаемых, так и несжимаемых течений. Большинство работ выполнено на криволинейных структурированных сетках. Использование же произвольных неструктурированных сеток привносит дополнительные сложности, связанные с определением поверхности обтекаемого объекта на сетке, особенно при его перемещении. Адаптация метода Бринкмана применительно к моделированию турбулентных течений, в том числе к вихреразрешающим RANS-LES подходам, началась относительно недавно (см., например, [2]– [6]).

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

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

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

2. МАТЕМАТИЧЕСКАЯ МОДЕЛЬ И ЧИСЛЕННЫЙ МЕТОД

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

Вихреразрешающее моделирование турбулентного течения около двухэлементного крылового сегмента проводится с помощью новой модификации гибридного RANS-LES подхода DDES [9], [10], ускоряющего переход к развитой турбулентности в слоях смешения.

Согласно разрабатываемой методике, для моделирования вспомогательного тела с изменяемой геометрией применяется метод погруженных границ, а основного – традиционный подход с использованием согласованной с границей сетки, то есть с определением поверхности обтекаемого твердого тела граничными сеточными узлами. Таким образом, в рассматриваемой задаче решение ищется в области ${{\Omega }_{F}}$, внешней по отношению к крыловому профилю, и во внутренности предкрылка ${{\bar {\Omega }}_{S}}$. Метод погруженных границ, а именно, метод Бринкмана штрафных функций здесь применяется для обеспечения выполнения граничного условия Дирихле для скорости ${\mathbf{u}} = 0$ и эффективной вязкости $\nu = 0$ на поверхности предкрылка $\partial {{\Omega }_{S}}$. Геометрию “погруженного” предкрылка определяет характеристическая функция $\chi $:

Система газодинамических уравнений, включающая DDES модель и записанная в дивергентной форме относительно консервативных переменных ${{\left( {\rho ,\rho {{u}_{i}},E,\rho \nu } \right)}^{{\text{т}}}}$, имеет вид

$\tilde {\mu } = {{\mu }_{{mol}}} + \rho \nu {{f}_{{\nu 1}}}$, ${{\mu }_{{mol}}}$ – коэффициент ламинарной вязкости,

${{\tau }_{{ij}}} = 2\left( {\frac{{\partial {{u}_{i}}}}{{\partial {{x}_{j}}}} + \frac{{\partial {{u}_{j}}}}{{\partial {{x}_{i}}}}} \right) - \frac{2}{3}\frac{{\partial {{u}_{i}}}}{{\partial {{x}_{i}}}}{{\delta }_{{ij}}}$ – компоненты тензора вязких напряжений;

${{q}_{i}} = \left( {\tilde {\mu } + \frac{\gamma }{{\Pr }}{{\mu }_{{mol}}}} \right)\frac{{\partial \varepsilon }}{{\partial {{x}_{i}}}}$ – тепловой поток, ${\text{Pr}} = 0.72$ – число Прандтля;

${{D}_{\nu }}$ и ${{G}_{\nu }}$ – члены, отвечающие за диффузию и генерацию турбулентности;

${{f}_{{\nu 1}}}$ и ${{f}_{w}}$ – демпфирующие функции;

${{l}_{T}}$ – пространственный масштаб турбулентности;

$\eta $ – коэффициент проницаемости.

Выражения для ${{D}_{\nu }}$ и ${{G}_{\nu }}$, а также демпфирующих функций ${{f}_{{\nu 1}}}$ и ${{f}_{w}}$ приведены в [11]. Масштаб турбулентности ${{l}_{T}}$ определяется в соответствии с DDES методом [12].

Источниковые члены в уравнениях системы (1) обеспечивают быструю релаксацию (порядка $\exp ( - {t \mathord{\left/ {\vphantom {t \eta }} \right. \kern-0em} \eta })$, $\eta \ll 1$) решения для скорости и эффективной вязкости к нулю в точках сетки, находящихся в области ${{\bar {\Omega }}_{S}}$, т.е. там, где характеристическая функция $\chi \ne 0$.

Для достижения приемлемой точности метода погруженных границ для моделирования граничного условия на поверхности предкрылка необходимо обеспечить достаточное сгущение расчетной сетки в ее окрестности. С этой целью используется стационарная адаптация расчетной сетки согласно вариационному методу c контролем качества сетки [7], [13]. Заметим, что используемая формулировка метода удобна для эффективного распараллеливания.

Как уже было сказано, описание поверхности обтекаемого тела, имитируемого методом погруженных границ, дается характеристической функцией $\chi $, заданной в узлах расчетной сетки. Вместе с тем в гибридных вихреразрешающих моделях, описывающих пристенные турбулентные течения, необходимым параметром является расстояние до стенки, которое определяет масштаб турбулентности ${{l}_{T}}$. Расстояние до стенки, а также его градиент используются также при построении неявной обобщенной функции, согласно которой проводится сеточная адаптация. Чтобы избежать трудоемкой процедуры вычисления расстояния до поверхности на каждом временном слое, что необходимо при перемещении тела, предлагается задавать функцию линий уровня на специальной интерополяционной решетке, движущейся вместе с телом по расчетной сетке [14]. Для удобства хранения такая решетка в работе строится на сетке типа восьмеричного дерева. Пример интерополяционной решетки для предкрылка приводится на фиг. 1.

Фиг. 1.

Интерполяционная решетка, определяющая геометрию предкрылка: а – общий вид, б – участок около верхней кромки предкрылка.

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

Для проведения численных расчетов применяется вершинно-центрированная конечно-объемная EBR (Edge-Based Reconstruction) схема [15]. С теоретической точки зрения, на произвольной неструктурированной сетке данная схема обладает вторым порядком точности, однако обеспечивает точность в смысле величины ошибки заметно выше по сравнению с традиционными схемами годуновского типа второго порядка. Это достигается благодаря квазиодномерной реконструкции переменных на расширенных реберно-ориентированных шаблонах. Реконструкции строятся таким образом, что на трансляционно-симметрических сетках (сетках, остающихся неизменными при перемещении на вектор их любого ребра) EBR схема совпадает с конечно-разностным методом высокого порядка (пятого или шестого).

Описанная математическая модель, все описанные алгоритмы и численный метод реализованы в программном комплексе NOISEtte.

3. АДАПТАЦИЯ И ДЕФОРМАЦИЯ СЕТКИ ДЛЯ ЗАДАНИЯ ГЕОМЕТРИИ ПРЕДКРЫЛКА И ЕЕ МОДИФИКАЦИИ

3.1. Расчетные сетки и их построение

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

Первая сетка (BF или Body-Fitted) представляет собой неструктурированную сетку, состоящую из треугольных призм и гексаэдральных слоев вблизи поверхности обтекаемых тел. Поверхности предкрылка и крыла описываются граничными узлами расчетной сетки. Фрагмент сетки в поперечном сечении вблизи предкрылка приведен на фиг. 2. Полная трехмерная сетка состоит из 14.7 млн. узлов и 17.7 млн. элементов.

Фиг. 2.

Фрагмент расчетной сетки BF, согласованной с поверхностью предкрылка.

Вторая сетка (IB или Immersed-Boundary) строится с применением стационарной сеточной адаптации следующим образом. Область предкрылка и окружающих его особенностей течения сначала покрывается квазиравномерной сеткой, как показано на фиг. 3 слева. Далее она адаптируется к поверхности предкрылка согласно вариационному методу [7]. Результат адаптации демонстрируется на фиг. 3 справа. Трехмерная сетка IB так же, как и сетка BF, состоит из треугольных призм и гексаэдральных пограничных слоев и содержит 16.8 млн. узлов и 19.4 млн. элементов.

Фиг. 3.

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

3.2. Деформация сетки при модификации внутренней поверхности предкрылка

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

Фиг. 4.

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

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

1. Определяется дискретная область деформации расчетной сетки ${{\Omega }_{h}}$ с криволинейными неподвижными границами ${{{\mathbf{R}}}_{L}}$и ${{{\mathbf{R}}}_{R}}$, соединяющими концевые точки закрылка A и B. Область деформации окружает внутреннюю поверхность предкрылка ${{{\mathbf{R}}}_{C}}$ (см. фиг. 4б).

2. Внутри области деформации строятся базовые узлы ij решетки. Узлы решетки задаются радиус вектором ${\mathbf{r}}$ и вычисляются следующим образом.

2.1. Границы ${{{\mathbf{R}}}_{L}}$и ${{{\mathbf{R}}}_{R}}$ равномерно разбиваются на N отрезков, обеспечивая взаимно-однозначное соответствие между узлами левой и правой границ и определяя, таким образом, сеточные линии индекса $i$ искомой решетки ($i = 0,\; \ldots ,\;N$).

2.2. Граница ${{{\mathbf{R}}}_{L}}$ определяет левую сеточную линию индекса $j = 0$. Внутренняя граница ${{{\mathbf{R}}}_{C}}$ предкрылка определяет центральную сеточную линию решетки индекса $j = MC$. $MC$ – число сеточных линий индекса $j$ слева от ${{{\mathbf{R}}}_{C}}$. Граница ${{{\mathbf{R}}}_{R}}$ задает правую сеточную линию индекса $j = M$, где $M = MC + MR$, $MR$ – число сеточных линий индекса $j$ справа от ${{{\mathbf{R}}}_{C}}$. Пересечение сеточных линий индексов $i$ и $j$ определяет базовые узлы ${{{\mathbf{r}}}_{{i,j}}}$ для $j = 0$, $MC$, $M$и $0 \leqslant i \leqslant N$. Заметим, что координаты узлов ${{{\mathbf{r}}}_{{0,j}}} = {{{\mathbf{r}}}_{A}}$ и ${{{\mathbf{r}}}_{{N,j}}} = {{{\mathbf{r}}}_{B}}$ не зависят от индекса $j$.

3. Для $0 \leqslant i \leqslant N$ вычисляются расстояния между левой и центральной сеточными линиями ${{D}_{i}} = \left| {{{{\mathbf{r}}}_{{i,0}}} - {{{\mathbf{r}}}_{{i,MC}}}} \right|$ и определяется единичный вектор направления сеточных линий индекса $i$: ${{{\mathbf{d}}}_{i}} = {{\left( {{{{\mathbf{r}}}_{{i,0}}} - {{{\mathbf{r}}}_{{i,MC}}}} \right)} \mathord{\left/ {\vphantom {{\left( {{{{\mathbf{r}}}_{{i,0}}} - {{{\mathbf{r}}}_{{i,MC}}}} \right)} {{{D}_{i}}}}} \right. \kern-0em} {{{D}_{i}}}}$.

4. Определение модуля скорости движения центральной сеточной линии ($j = MC$) вдоль сеточных линий индекса $i$: ${{V}_{i}} = {{\alpha \,{{D}_{i}}} \mathord{\left/ {\vphantom {{\alpha \,{{D}_{i}}} {\left( {{{T}_{2}} - {{T}_{1}}} \right)}}} \right. \kern-0em} {\left( {{{T}_{2}} - {{T}_{1}}} \right)}}$, где ${{T}_{1}}$ и ${{T}_{2}}$ – время начала и конца перемещения внутренней поверхности предкрылка. Задаваемый параметр $0 < \alpha < 1$ определяет относительную величину смещения внутренней поверхности предкрылка в области между кривой ${{{\mathbf{R}}}_{L}}$ и начальным положением кривой ${{{\mathbf{R}}}_{С}}$.

На втором этапе алгоритма строится двумерная решетка и определяются новые координаты ее базовых узлов на момент времени $t + \Delta \,t$, где $t \geqslant {{T}_{1}}$, а $\Delta \,t$ – шаг интегрирования по времени системы уравнений. Этот этап включает в себя следующие шаги.

5. Определение новых координат базовых узлов решетки. Сеточные линии с узлами ${{{\mathbf{r}}}_{{i,0}}}$ и ${{{\mathbf{r}}}_{{i,M}}}$ остаются неподвижными, а центральная сеточная линия ($j = MC$) смещается на величину $\Delta \,{{{\mathbf{r}}}_{{i,C}}} = {{V}_{i}}{{{\mathbf{d}}}_{i}}\Delta \,t$ вдоль сеточных линий индекса $i$. Таким образом, новые координаты внутренней поверхности предкрылка переопределяются по формуле ${{{\mathbf{r}}}_{{i,MC}}}: = {{{\mathbf{r}}}_{{i,MC}}} + \Delta {{{\mathbf{r}}}_{{i,C}}}$

6. На основе базовых узлов, определенных на шаге 2.2 и шаге 5, промежуточные узлы решетки вычисляются по следующим параметрическим формулам:

$\begin{gathered} {{{\mathbf{r}}}_{{i,j}}} = (1 - {{s}_{j}}){{{\mathbf{r}}}_{{i,L}}} + {{s}_{j}}{{{\mathbf{r}}}_{{i,MC}}},\quad 1 \leqslant j \leqslant M - 1,\quad 1 \leqslant i \leqslant N - 1, \\ {{{\mathbf{r}}}_{{i,L}}} = T\left\{ \begin{gathered} {{{\mathbf{r}}}_{{i,0}}},\quad j < MC, \hfill \\ {{{\mathbf{r}}}_{{i,M}}},\quad j \geqslant MC, \hfill \\ \end{gathered} \right. \\ \end{gathered} $(**)

где параметром является сеточная функция растяжения ${{s}_{j}}$ [16, с. 123], меняющееся в диапазоне от 0 до 1, и определенная как

${{s}_{j}} = P{{\eta }_{j}} + \left( {1 - P} \right)\left( {1 - \frac{{\operatorname{th} \left[ {Q\left( {1 - {{\eta }_{j}}} \right)} \right]}}{{\operatorname{th} Q}}} \right),\quad {{\eta }_{j}} = \left\{ \begin{gathered} \frac{j}{{MC}},\quad j = 0,\; \ldots ,\;MC - 1, \hfill \\ \frac{{M - j}}{{MR}},\quad j = M,\; \ldots ,\;MC. \hfill \\ \end{gathered} \right.$

Здесь $P$ и $Q$ – числовые константы, обеспечивающие контроль распределения узлов решетки вдоль отрезка $[{{{\mathbf{r}}}_{{i,L}}},{{{\mathbf{r}}}_{{i,MC}}}]$. Так, при $P = 1$ получаем равномерную решетку слева и справа от центральной сеточной линии (см. фиг. 4б). При $P = 1.8$ и $Q = 2$ имеем неравномерную решетку, сгущающуюся к центральной сеточной линии ($j = MC$) с левой и правой сторон (см. фиг. 5а).

Фиг. 5.

а – начальная неравномерная решетка на момент времени $t = {{T}_{1}}$, б – финальная решетка на момент времени $t = {{T}_{2}}$.

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

7. Для каждого узла исходной расчетной сетки ${\mathbf{r}}_{k}^{{node}} \in {{\Omega }_{h}}$ определяется принадлежность ячейкам решетки на основе двоичного дерева поиска.

8. Используя кусочно-линейную интерполяцию по скоростям вершин ячейки решетки, вычисляются скорости узлов сетки ${\mathbf{V}}_{k}^{{node}}$, попавших в эту ячейку. Далее определяются координаты узлов сетки на новом шаге по времени, как ${\mathbf{\hat {r}}}_{k}^{{node}} = {\mathbf{r}}_{k}^{{node}} + {\mathbf{V}}_{k}^{{node}}\Delta t.$

На последнем шаге алгоритма происходит проверка условия $t + \Delta t \leqslant {{T}_{2}}$ и переход к шагу 5.

Пример работы алгоритма построения дополнительной решетки, приведенной на фиг. 4 и фиг. 5, использовал следующие значения входных параметров: $N = 20$, $MC = 8,$ $MR = 10,$ ${{T}_{1}} = 0,$ ${{T}_{2}} = 0.5,$ $\alpha = {1 \mathord{\left/ {\vphantom {1 3}} \right. \kern-0em} 3},$ $\Delta \,t = 0.01$.

Таким образом, построенная на момент времени ${{T}_{1}} < t \leqslant {{T}_{2}}$ новая деформированная сетка может быть использована для расчета обтекания предкрылка новой геометрической формы. На фиг. 6 изображена начальная и финальная стадия деформации и движения, используемых в приводимых ниже расчетах, сетки и решетки ($N = 101$, $MC = 16,$ $MR = 17$), на фиг. 6 изображены не все сеточные линии индекса i).

Фиг. 6.

Слева – сетка и решетка перед началом деформации, справа – сетка и решетка после деформации.

4. ЧИСЛЕННЫЕ РЕЗУЛЬТАТЫ

4.2. Анализ численных результатов

Вихреразрешающему моделированию турбулентного течения около крылового сегмента в рамках гибридного метода DDES предшествовали расчеты по модели RANS. Как видно из фиг. 7, где представлено поверхностное распределение коэффициента давления, расчеты на всех сетках показали очень близкие результаты. Отметим, что в расчетах использовалась еще одна сетка – сетка IB-BF, которая строилась путем заполнения внутренности предкрылка в сетке BF.

Фиг. 7.

Поверхностное распределение коэффициента давления для крыла и предкрылка.

Полученное подходом RANS стационарное поле течения взято в качестве начального распределения для дальнейших расчетов с использованием модели DDES. На фиг. 8 приведены результаты вихреразрешающего моделирования в виде полей завихренности и производной давления по времени: для BF сетки (фиг. 8а), а также для IB сетки при исходной и деформированной форме внутренней поверхности предкрылка (фиг. 8б, 8в соответственно). В целом поля течения схожи между собой. Вместе с тем при использовании метода погруженных границ (фиг. 8б, 8в) заметны высокочастотные осцилляции вблизи поверхности предкрылка. Хотя их возможное влияние находится в области высоких частот и, как показывает опыт, существенного влияния на течение в целом они не оказывают, избежать их помогает более аккуратная настройка параметров метода погруженных границ и, в частности, снижение проницаемости обтекаемого тела или добавление искусственной вязкости внутрь обтекаемого тела, где результаты не имеют физического смысла.

Фиг. 8.

Мгновенные поля завихренности и производной давления по времени: а – BF сетка, б – IB сетка с исходной формой предкрылка, в – IB сетка с деформированной формой предкрылка.

Более корректное сопоставление результатов вихреразрешающего моделирования, полученных традиционным подходом на согласованной с поверхностью BF сетке и при использовании метода погруженных границ на IB сетке, дает анализ спектрального состава пульсаций давления в контрольных точках области, расположение которых дано на фиг. 9. Приведем спектральную плотность мощности пульсаций давления в наиболее чувствительных к методу точках S1 и S6 (см. фиг. 10). Данные точки лежат на внутренней поверхности предкрылка, которые в случае применения метода погруженных границ попадают на неявно заданную границу раздела сред. Как видно из приведенных графиков, полученные результаты хорошо соотносятся друг с другом, а также с результатами вычислительного эксперимента, проведенного авторами данной постановки [18].

Фиг. 9.

Расположение контрольных точек.

Фиг. 10.

Спектральная плотность мощности пульсаций давления в точках S1 и S6 на внутренней поверхности предкрылка.

Аэроакустические свойства двухэлементного крылового сегмента характеризует, в первую очередь, генерируемый им шум в дальнем поле. В данной работе шум в дальнем поле вычисляется с помощью интегрального метода Фокса–Уильяма Хокинга путем пространственно-временного интегрирования данных, накопленных в процессе расчета на поверхности, окружающей конфигурацию и расположенной на небольшом расстоянии от нее. На фиг. 11 показаны спектры уровня звукового давления в контрольных точках на удалении 2 метра, или примерно 7 длин хорды крыла, от исследуемой конфигурации. Можно видеть, что численные результаты, полученные традиционным подходом и методом погруженных границ (зеленая и красная линии соответственно) для исходной формы предкрылка весьма хорошо согласуются друг с другом, а также с экспериментальными данными, отмеченными треугольниками и квадратами. Также заметно повышение амплитуды акустического излучения на пиках частот в диапазоне до 5 КГц в случае деформированной внутренней поверхности предкрылка. Особенно отчетливо это видно в точке измерений, находящейся над крылом. Это является ожидаемым результатом, так как с демонстрационной целью проводилась достаточно произвольная деформация относительно хорошо выверенной конструкции. Также в спектре, соответствующем деформированному предкрылку, заметно сильное завышение акустического давления на высоких частотах в районе 10 КГц. Этот эффект, вероятнее всего, имеет численную природу, связанную с повышением искусственной проницаемости предкрылка после его деформации, уменьшающей его толщину. Следует отметить, что коэффициент проницаемости, изначально подобранный для исходной формы предкрылка, в процессе счета не изменялся, в то время как толщина предкрылка заметно уменьшилась.

Фиг. 11.

Спектральный состав уровня звукового давления в контрольных точках дальнего поля, расположенных на расстоянии 2 м от точки начала крыла в направлениях на 100° по часовой стрелке (а) и 90° против часовой стрелки относительно крыловой хорды.

Проведенные расчеты демонстрируют работоспособность созданной вычислительной методики и показывают ее эффективность, которая особенно наглядно проявляется в выигрыше во времени расчета. Так, при использовании гибридного вихреразрешающего подхода DDES для выхода на статистически стационарное турбулентное течение, стартуя с решения RANS, требуется порядка 10–15 безразмерных единиц времени, нормированных на длину хорды крыла и скорость набегающего потока. Для накопления корректных данных после этого необходимо еще, как минимум, 30 единиц времени. На деформацию поверхности предкрылка было отведено 0.5 единиц времени. К окончанию плавного процесса деформации турбулентное течение сразу выходило на квазистационарный режим. Таким образом, существенное, порядка 25% от счета задачи при фиксированной конфигурации, время на установление тратится только один раз на старте с RANS решения для начальной геометрии. Затем переход от одной конфигурации к другой происходит практически без временных потерь на установление. На фиг. 12 проиллюстрирован переход от одной формы предкрылка к другой на примере графиков пульсаций давления по времени в двух контрольных точках, расположенных над и под крылом (см. фиг. 9).

Фиг. 12.

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

6. ЗАКЛЮЧЕНИЕ

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

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

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

Работа была выполнена с использованием оборудования центра коллективного пользования “Комплекс моделирования и обработки данных исследовательских установок мега-класса” НИЦ “Курчатовский институт”, http://ckp.nrcki.ru/.