Применение уравнений Эйлера для моделирования вихревых течений· Текст научной статьи по специальности «Физика»

УЧЕНЫЕ ЗАПИСКИ ЦАГИ Том XXX 1999

№1—2

УДК 629.735.33.015.3.025. 1:532.526

ПРИМЕНЕНИЕ УРАВНЕНИЙ ЭЙЛЕРА ДЛЯ МОДЕЛИРОВАНИЯ ВИХРЕВЫХ ТЕЧЕНИЙ*

Н. А. Владимирова, В. С. Сакович

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

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

♦Работа выполнена при финансовой поддержке Международного научно-технического центра (проект № 201-95).

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

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

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

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

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

1. Общая постановка задачи и метод расчета. Дня расчета обтекания трехмерных аэродинамических конфигураций (изолированное крыло или компоновка «крыло — фюзеляж») и структуры вихревых образований в спутном следе используется численное решение системы уравнений Эйлера:

дЁ дРдСт — + — + — + — = 0. (1.1)

д1 дх ду дг

Составляющие векторов <?, Ё, Р и <5 зависят от плотности р, проекций

скорости и, V, м> на оси координат, давления Р и энергии е:

1 т '

q=[ р, ри, ри, pw, е]

Ё = [ри, ри2 + Р, р UV, рии>, и(е + Р)]т, F = [pv, рuv, pv2 + Р, рvw, v(e + .Р)]3", G = [pw, puw, pvw, pw2 + P, w(e + P)]T.

(1.2)

Давление P определяется из уравнения состояния совершенного газа с показателем адиабаты к:

P = (k-l)[e-p/2(u1+v2 + w2)l

Переход к безразмерным переменным в системе уравнений (1.1) производится таким образом, чтобы в безразмерных переменных вид уравнений не менялся: плотность р и давление Р относятся к своим значениям в невозмущенном потоке Рда И Р„, энергия е ----- К Рп, скорость к величине

\Ю 1/9

(Рда / Pqo ) . Единицей времени становится величина L / (/ р^ ) .

На поверхности обтекаемого тела ставится условие непротекания:

ипх + vrty + wn2 = 0, (1.3)

где п = (пх,пу,п2) — нормаль в точке границы. При удалении на бесконечность вверх по потоку р -> рда, е бда, К Кда, вниз по потоку Р-> Рж.

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

Дискретизация уравнений Эйлера (1.1) производится в рамках схемы конечного объема и основывается на представлении искомых составляющих вектора q и потоков Ё, F, G в виде кусочно-постоянных функций и замене системы дифференциальных уравнений (1.1) интегралами по объему каждой ячейки. По теореме Остроградского — Гаусса:

-^jjjqdxdydz + jj^nx +Fny +Gnz)ds = 0, (1.4)

Т S

где т — объем ячейки, s — граница ячейки сетки. Выбор квадратурной формулы для вычисления поверхностного интеграла в (1.4) производится из условия, чтобы равномерный поток Ё(х, у, z) = const, F(x, y,z) = const, G(x, у, z) = const являлся точным решением системы разностных уравнений. Таким требованиям удовлетворяет квадратурная формула прямоугольников с опорными точками в центрах граней ячеек, в которой нормаль определяется как векторное произведение диагоналей грани ячейки. Доопределение величин потоков при вычислении интегралов по границам ячейки производится путем осреднения величин составляющих векторов Ё,Р и G из соседних ячеек, что не приводит к существенному размазыванию скачков уплотнения (так как потоки Ё,Ё и G, в отличие от q, непрерывны в поле течения даже при наличии скачков). На поверхности тела, вследствие выполнения условия непротекания (1.3), ненулевой вклад в потоки через границу дает только давление.

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

^(b-dxdydz = ^Q.-dfi, (1.5)

X S

где Q — вихревой поток:

' 0 -W V

п= W 0 -и

и 0

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

На внешней границе расчетной области ставятся условия излучения, которые должны подавлять отражение сигнала (более детально см. п. 3 данной работы).

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

- а,

др 5р др

дх ду дг

д(р и) д(р и) д(ри)

дх ду дг

д^) д^) а(ру)

дх ду дг

д(рмО а(рм') д(р'И')

дх дг

д(е) д(е) д(е)

дх ду дг

дър д3р а3р

дх3 ду3 а*3

д\ри) а3(р и) а3(р и)

дх3 ду3 дг3

а3(р и) а3(рг>) а3(ро)

дх3 ду3 а*3

д3(рм>) д3(ру^) а3(Рм')

дх3 ду3 дг3

д\е) д3(е) д\е)

дх3 ду3 дг3

Ах

Ау

■ Аг

(1.6)

Ч4) Ах3

Ау3

■Аг3

где в^~§гас1.Р, =шах(0,1/32-е^).

Впервые введение диссипации таким способом было предложено в [2]. «Борьба» с разрывами в решении (в отсутствие вязкостной диффузии на разностной сетке нельзя воспроизвести волны, длина которых меньше шага сетки) приводит к краевой задаче с уравнениями, содержащими дисси-

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

Авторы данной работы ограничиваются рассмотрением стационарных задач. Стационарное решение достигается методом установления. Для ускорения сходимости используется многосеточный алгоритм. Применение многосеточного алгоритма требует особенно сильного затухания в коротковолновом спектре, необходимого для сглаживания осцилляций, возникающих при переходах с грубой сетки на мелкую, поэтому кроме линейной диссипации в метод расчета вводится нелинейный механизм подавления осцилляций в областях больших градиентов параметров [4]. В данной работе используется трехсеточный алгоритм, который заключается в последовательном переходе на более грубые сетки с одним временным шагом на каждой из них и интерполяцией результатов от самой грубой до самой мелкой. Такая методика использована в работах [5]—[8] и дает приблизительно шестикратное уменьшение времени счета.

Для расчета течения около изолированного крыла используется сетка типа С — Н, считающаяся наиболее удачной для этого класса задач. В каждой плоскости Z= const вдоль размаха строится двумерная сетка типа С вокруг сечения крыла с отображением внешности сечения на область, близкую к полуплоскости (конформное преобразование типа квадратного корня). Затем сдвиговое преобразование отображает эту область на полуплоскость.

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

Более подробно детали численного метода, а также вопросы тестирования программы рассматриваются в работе [6]. Примеры его успешного использования для решения ряда задач прикладной аэродинамики, связанных с развитием вихревых систем, могут быть найдены в работах [9]—[13].

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

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

Математическая постановка рассматриваемой задачи предполагает дополнительно выполнение граничных условий непротекания (1.3) на твердой бесконечной плоскости, моделирующей поверхность земли, и генерацию расчетной сетки в конечной области около конфигурации «изолированное крыло — твердая граница».

2.1. Алгоритм построения расчетной сетки. Трехмерная расчетная сетка около конфигурации «крыло — плоскость земли» строится как суперпозиция двумерных сеток в плоскостях Z= const (ось OZ направлена вдоль размаха крыла). Декартовы координаты узлов сетки в этих плоскостях вычисляются по аналитическим формулам расчета сетки типа С, основанным на суперпозиции простых конформных преобразований и преобразования сдвига. Для этого с помощью преобразования

хх =x(a + pxI£),

У\=У( а + рх),

где а = ■

71 ХТЕ XLE в- 71 Г 1 її

ХТЕ ~XLE {HLE НТЕ, 9 Р ХТЕ ~XLE <НтЕ HLE ,

в каждом

продольном сечении крыла Z~ const точки передней и задней кромок с координатами (*/,£, Hie) и (хте> Нте) соответственно (в исходной декартовой системе координат х, у, z) переводим на высоту я. Затем применяем конформное преобразование комплексной плоскости toj =xj +iy\ на комплексную плоскость а>2 = х2 + *У2 следующего вида: coj = а>2 + е-£°2, которое переводит внешнюю область, ограниченную контуром профиля и плоскостью земли, на область, близкую к полосе шириной л. Для нахождения обратного преобразования a>2 =/((0i) и образа контура профиля ус(х2) в плоскости со2 используется приближенный метод Ньютона. Далее простым сдвиговым преобразованием

х3 = х2,

Уз=У2 -Ус(х2)(п~У2)

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

На рис. 1 показана форма сеточных линий в плоскости симметрии (Z = 0) стреловидного крыла, расположенного вблизи поверхности земли.

Рис. 1. Фрагмент сетки в вертикальной плоскости симметрии стреловидного крыла, расположенного вблизи твердой границы

2.2. Результаты расчетов. Компьютерное моделирование интерференции крыла и вихревого следа за ним с твердой поверхностью проводилось на сетевой рабочей станции с быстродействием 20—25 МРЬОРв. Для исследований было взято стреловидное крыло самолета В 747. В расчетной области вокруг крыла и бесконечной плоскости, моделирующей поверхность земли, с помощью алгоритма, приведенного выше, была построена расчетная сетка с числом узлов 128 х 24 х 32 (в физическом полупространстве). Протяженность расчетной области составляла величину порядка 6—8 хорд во всех координатных направлениях от его поверхности (см. рис. 1). Для достижения сходимости осуществлялось около 1000 итераций, что приводило К получению стационарного решения методом установления с точностью 10'7—10'8 (для максимальной невязки во всем поле течения). Результаты расчета интегральных аэродинамических характеристик подъемной силы (су) и индуктивного сопротивления (сх.) крыла В 747 для дозвукового режима обтекания

Мда =0,2 (число Маха набегающего потока) и а = 7° (угол атаки) в зависимости от расстояния плоскости крыла от поверхности земли Н (в долях корневой хорды) приведены в табл. 1. Расчетные значения коэффициентов су и сх показывают увеличение несущих свойств крыла и уменьшение его

индуктивного сопротивления при приближении к земле на расстояния порядка 0,5 + 2 хорд, что находится в полном соответствии с давно известными результатами по учету влияния близости земли.

Результаты расчета абсолютной величины завихренности |й| в поле течения вниз по потоку от задней кромки крыла для 4 различных положений крыла над твердой

Таблица 1

М„ = 0,2, а = 7°

Я Су СХ1

2,30 0,761 0,0292

1,53 0,771 0,0276

0,77 0,800 0,0240

0,46 0,842 0,0213

К)

Ov

В 747

М=0,2

а=7°

Рис. 2. Изолинии завихренности (изостеры |<в|= const) в горизонтальной плоскости течения (вид сверху) для 4 значений высоты крыла над

поверхностью земли

поверхностью представлены на рис. 2. Распределение завихренности показано в срединной горизонтальной плоскости (плоскость вихревого следа) в форме изолиний (изостер) в диапазоне от нулевого минимального значения до максимального (0,45) с шагом 0,05. Видно, как при приближении крыла к границе земли (два верхних фрагмента на рис. 2, Н = 0,46 и 0,77 ) интенсивность вихревой пелены, сходящей с задней кромки, и зона ее формирования увеличиваются. В окрестности конца крыла из такой вихревой структуры формируется концевой вихрь. На малых расстояниях от границы такое вихревое образование является достаточно мощным и устойчивым (левый верхний фрагмент, Н = 0,46). С увеличением расстояния между крылом и твердой границей интенсивность сходящего концевого вихря заметно уменьшается (правый нижний фрагмент, Н = 2,30).

3. Дифракция вихрей. Численно моделируется воздействие вихревых образований, рассматриваемых в качестве дальнего вихревого следа, сформированного «самолетом-генератором», на структуру ближнего поля течения и аэродинамические характеристики стреловидного несущего крыла. В рамках полных трехмерных уравнений Эйлера решается нелинейная стационарная задача для случая продольного входа крыла в вихревое образование (продольная дифракция вихря на препятствии). Для инициализации поля скоростей в дальнем следе за летательным аппаратом используется эмпирическая модель [14]—[16], краткое описание которой будет дано ниже.

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

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

Продольные оси двух бесконечных вихрей одинаковой интенсивности разного знака предполагаются направленными вдоль вектора скорости набегающего потока и параллельными вертикальной плоскости симметрии задачи. На части внешней границы, расположенной на расстояниях порядка 3—5 хорд от обтекаемого крыла по всем направлениям, где дозвуковой поток втекает в область расчета, граничные условия записываются в следующей форме. Задаются величина энтропии £, тангенциальная скорость иТ и инвариант Римана Яех( = + 2а / (к -1), где ьп — нормальная к границе ско-

рость, а — скорость звука, к — показатель адиабаты. Второй инвариант Римана И.;п1 =ип-2а / (к -\) экстраполируется на границу из расчетной области. На части границы, где поток вытекает из области, задается только Яех1, а К;п|:, энтропия и тангенциальная скорость экстраполируются из области расчета. Рассмотрим ячейку сетки, одна из граней которой АВСИ является частью внешней границы. Тогда нормаль к внешней границе может быть вычислена при помощи векторного произведения диагоналей грани:

- (гА^гс)х(гв-^)

I (гА-гс)х(?в-?0)\’

Компоненты скорости и, V, м> вне расчетной области вычисляются в системе координат, связанной с набегающим потоком, по следующим формулам:

и = Ко;

V = V (у1(у-уо)2+(г-г0)2) -Т....—

у/(у~уо)2 +(2~2о)

ивихр(у1(У-У0)2+(2 + 20)2) I '

у/(У~Уо)2 +(2 + го)

т = -увихр ( ^(У-Уо)2 +(2~2о)2) -I---

(у-уо)2+(2~2о)2

+увихр[а/ (У -Уо)2 +(2 + го)2 ] —

У-Уо Уо)2+(2 + 2о)2

(3.1)

где — скорость набегающего потока, ивихр — характерная скорость

вихря в заданной точке пространства (х, у, г), (уо. г0) — координаты вихря в этой системе координат. Нормальная составляющая скорости, необходимая для определения инварианта Римана, вычисляется по формуле ип = (V, Я), тангенциальная — по формуле = V - • п. Величина скоро-

сти звука вычисляется из условия постоянства энтальпии по формуле

о = ^(*-1)(£-(и2+г/2 + м'2)/2),

где Е —энтальпия на бесконечности перед крылом, которая считается постоянной в набегающем потоке.

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

Остановимся несколько подробнее на используемой в данной работе эмпирической модели следа [14]—[16].

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

— В центре каждого из них наблюдается турбулентное вихревое ядро, которое начинает формироваться сразу за задней кромкой крыла. Это ядро является центром, на который в дальнейшем наматывается вихревая поверхность. В работе [14] приведен эмпирический закон роста его размеров Яс по мере удаления от крыла на расстояние х:

с nL

і

0,421 + 0,02 LGx

Poo^S2

где S' и I — площадь и размах крыла, G — вес самолета, Vx — скорость полета, Рад — плотность атмосферы.

— Поведение окружной скорости в следе таково, что сначала при увеличении расстояния от центра она растет от нуля до некоторого максимального значения, а затем медленно падает (приблизительно как 1/Л). Положение максимума определяет размер турбулентного ядра Я* и его циркуляцию Ге. Согласно экспериментальным данным величина Гс составляет 0,2 — 0,7 от полной циркуляции.

— Наличие в атмосфере турбулентной диффузии приводит к падению полной циркуляции Г, которое описывается эмпирическим соотношением [15]:

Г = Г0ехр

где Г0 — полная циркуляция, сошедшая с половины крыла, q — уровень турбулентности атмосферы, Ъ — расстояние между центрами вихревых образований.

В данной работе для вычисления поля скоростей и положения вихрей в дальнем следе за летательным аппаратом используется эмпирическая модель [16], разработанная группой под руководством Г. Г. Судакова (ЦАГИ). Алгоритм [16] базируется на экспериментальных данных, аккумулирует в себе результаты [14], [15] и учитывает все основные факторы: вид распределенной нагрузки на крыле, рост турбулентного ядра вихря, влияние турбулентности и стратификации атмосферы, влияние поверхности земли и т. п. Для описания внешнего невязкого поля течения использована модель Бетца [17], скорректированная для учета вытесняющего действия турбулентного ядра. Для оценки размера и циркуляции турбулентного ядра вихря с учетом турбулентной диффузии в атмосфере, положения вихревых ядер с учетом стратификации атмосферы и времени разрушения следа (до момента касания двух вихревых нитей) используются результаты работы [18].

3.2. Результаты расчета. Эмпирическая методика [16] была использована для расчета структуры свернувшейся вихревой пелены за самолетом класса Боинг 747 в посадочной конфигурации (скорость посадки 70 м/с, высота полета 100 м, угол атаки 10°, угол отклонения закрылков 30°, с,, = 1,44, степень турбулентности атмосферы <7=1,0 м/с, стратификация атмосферы, определяемая частотой Бранта — Вайсала,

рф>

= 0,01 1/с, где g — ускорение свободного падения, у

вертикальная координата, расстояние между центрами вихрей Ъ = 42 м, начальная циркуляция Г0 = 593 м2/с). В рассчитанное при таких значениях параметров поле скоростей дальнего следа помещалось изолированное стреловидное крыло со следующими геометрическими параметрами: размах X = 16,2 м, удлинение X = 3,81, сужение г| = 1,33, стреловидность % = = 30°, в сечениях крыла установлен сверхкритический профиль с максимальной относительной толщиной с = 9%.

Ниже приводятся результаты численного моделирования аэродинамического взаимодействия вихревых образований за самолетом Боинг 747 (рассмотрены несимметричный случай одного вихревого жгута и симметричная задача для вихревой пары) с указанным стреловидным крылом (число Маха набегающего потока Мда = 0,2, угол атаки а = 10°). Рассмотрены два значения времени существования внешнего вихревого поля: / = 10 с («молодой» след ) и Г = 40 с («старый» след, полное время жизни следа до первого касания вихрей составляет Т = 65 с). Рассчитанное по методике [16] поле окружной скорости в следе подставляется в качестве значений увихр в граничные условия (3.1). Следует отметить, что «молодой» и «старый» следы существенно отличаются по значению максимальной окружной скорости на границе турбулентного ядра вихря: для времен существования следа 10 и 40 с это отличие составляет приблизительно 3,5 раза.

Все расчеты проводились на сгенерированной в расчетной области вокруг крыла расчетной сетке типа С — Н с числом узлов 128 х 24 х 48 (в полупространстве). Границы расчетной области удалены на расстояние 1,5 — 2 размаха крыла по всем направлениям от его поверхности. Для получения устойчивого стационарного решения (при этом достигается уменьшение невязки на 7 — 9 порядков) осуществлялось в среднем 800 — 1000 итераций в схеме установления.

Таблица 2

Эмпирическая модель следа (16], В 747, * = 10 с, I = 16,2 м, Ь = 42 м

Мю = 0,2, а =10° Эйлер Линейная теория

СУ £у/СуИЗОЛ СУ Су /СуИЗОл

Изолированное крыло 0,806 1,000 0,760 1,000

Высота вихрей К= 0 0,344 0,427 0,332 0,437

Высота вихрей К= 5 м 0,374 0,464 0,364 0,479

Высота вихрей Г= -5 м 0,372 0,462 0,362 0,476

В табл. 2 представлены полученные в расчетах значения коэффициента подъемной силы су рассматриваемого стреловидного крыла с размахом Ь= 16,2 м, помещенного в симметричное вихревое поле двух вихрей, расположенных на расстоянии Ъ = 42 м и моделирующих «молодой» (/ = 10 с) след за самолетом В 747. Рассмотрены три взаимных положения плоско-

стей крыла и следа по высоте: 7= 0; 5 м; -5 м. Как видно, воздействие смоделированного внешнего вихревого поля приводит к существенному (в 2 — 2,5 раза ) снижению несущих свойств крыла. В правой части табл. 2 для сравнения приведены результаты линейной теории, которая использует гипотезу «замороженности» вихрей. Отклик несущих свойств крыла на воздействие поля от двух вихрей практически (с точностью до 2 — 5% ) одинаков по результатам линейной теории (метод особенностей) и решения уравнений Эйлера, в рамках которого гипотеза «замороженности» вихрей a priori не закладывается. Эти и другие результаты авторов приводят к выводу, что в задачах дифракции, когда размер «препятствия» заметно (в 2—3 раза и более) меньше расстояния между вихрями, вихри в области ближнего поля около препятствия можно считать «замороженными», т. е. не учитывать обратное влияние препятствия на форму и структуру вихревых образований.

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

На рис. 3—5 и в табл. 3 представлены результаты расчетов для 5 случаев расположения вихревого жгута в плоскости крыла (7=0): достаточно далеко от крыла на расстоянии Z=21 м, вблизи концов крыла (Z = ± 8 м), посередине левого полукрыла (Z = 4 м) и в корневом сечении крыла (Z = = 0). Ось Z направлена вдоль размаха крыла. Значения коэффициентов су,

сх. и момента крена тх получены в результате воздействия «молодого»

(t =10 с) и «старого» (/ = 40 с) вихря, моделирующего след.

Таблица 3

Размах крыла L = 16,2 м, сутол = 0,806

Положение вихря, м с с с. тх

< = Юс < = 40 с < = Юс < = 40 с <= 10 с II О а

21 0,571 0,063 - 0,009

8 0,331 0,537 0,051 0,064 -0,016 -0,011

4 0,477 0,051 0,052

0 0,762 0,748 0,048 0,056 0,086 0,053

-8 1,223 1,068 0,070 0,045 - 0,002 - 0,007

Заметим, что знак вращения одиночного вихря, моделирующего в расчетах вихревой след, совпадает со знаком левого вихря, сходящего с левого конца крыла (т. е. вращение по часовой стрелке).

Анализ приведенных в табл. 3 результатов численного моделирования воздействия вихревого следа в виде одиночного вихря на обтекание несущего крыла позволил выявить два наиболее критичных с точки зрения безопасности полета режима резкого изменения несущих и моментных характеристик крыла: ЭТО резкое уменьшение (или увеличение) подъемной СИЛЫ крыла Су (в сравнении со значением с^30л для изолированного крыла, помещенного в равномерный безвихревой набегающий поток) в ситуации, когда ядро вихревого шнура «задевает» конец крыла, и появление значительного момента крена, если корневое сечение крыла «натыкается» на вихрь.

10 с ,

/ = 40 с

Изостеры в вертикальной плоскости на «входе» в расчетную область

Изостеры в вертикальной плоскости К2Г, пересекающей заднюю кромку крыла

Изостеры в вертикальной плоскости YZ на «выходе» из расчетной области

и>

-й-

Координаты вихря: У=0 (плоскость крыла), '/=0 (середина крыла) / = 10 с / = 40 с

Изостеры в вертикальной плоскости У2 на «входе» в расчетную область

Изостеры в вертикальной плоскости пересекающей заднюю кромку крыла

Изостеры в вертикальной плоскости У2 на «выходе» из расчетной области

Координаты вихря: У= 0 (плоскость крыла), Х= 8 м (левый конец крыла)

< = 10 с

/ = 40 с

Изостеры в вертикальной плоскости на «входе» в расчетную область

Изостеры в вертикальной плоскости У2, пересекающей заднюю кромку крыла

Изостеры в вертикальной плоскости К?на «выходе» из расчетной области Рис. 5

и>

-4

Проанализировать структуру возникающего течения и картину взаимной интерференции продольного одиночного вихря, моделирующего след, с вихревой системой, формируемой несущим крылом, можно на рис. 3—5. Информация о течении, полученная в результате численного моделирования, представлена на рис. 3—5 в форме изостер (линий постоянства абсолютной величины завихренности |co|=const)B вертикальных плоскостях х = const, соответствующих: сечению вблизи входной границы расчетной области, расположенному на расстоянии порядка размаха перед крылом (верхние фрагменты); сечению вблизи задней кромки крыла (средние фрагменты) и сечению на расстоянии размаха вниз по потоку от крыла (нижние фрагменты). Подобная численная информация о полях завихренности также может быть представлена в форме черно-белых или цветных изображений с непрерывным спектром.

Вихревые образования близкой к цилиндрической формы (верхние фрагменты на рис. 3—5) соответствуют одиночному вихрю, моделирующему дальний след за самолетом В 747. Изображенная на фрагментах конфигурация отражает структуру этого вихря на входе в расчетную область, полученную в результате расчета взаимной интерференции стреловидного крыла и дальнего следа путем решения уравнений Эйлера. Следует отметить, что в рассматриваемых примерах максимальное значение завихренности для изостер на рис. 3—5 равно 0,05 условных единиц и соответствует максимальной величине завихренности в центре одиночного вихря, в то время как вихревая система, генерируемая крылом, имеет на порядок большие значения завихренности (0,6—0,8). На рис. 3 вихрь, моделирующий след, вращается по часовой стрелке и «натыкается» на правый конец крыла, на котором формируется собственный вихрь с противоположным вращением. Их взаимодействие приводит к образованию достаточно сложной многоядерной структуры, дроблению вихрей и возникновению вторичных вихрей (левый нижний фрагмент на рис. 3). Указанные эффекты проявляются слабее в случае дифракции «старого» вихря (/ = 40 с) в сравнении с «молодым» вихрем (Г = 10 с, левые и правые фрагменты на рис. 3). По мере продвижения внешнего вихря к левому концу крыла, на котором формируется вихрь того же знака, т. е. вращение по часовой стрелке (рис. 4, 5), можно наблюдать процессы сближения и слияния вихревых ядер (рис. 5) и формирования единой вихревой системы с двумя близко расположенными ядрами.

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

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

2. Разработаны алгоритм и программа расчета аэродинамического взаимодействия одиночного вихревого образования или симметричной

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

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

ЛИТЕРАТУРА

1. Jameson A., Schmidt W., Turkel Е. Numerical solution of the Euler equation by finite volume methods using Runge — Kutta time-stepping schemes//AIAA Paper.— 1981, N 1259.

2. Jameson A. Transonic airfoil calculations using the Euler equations//Numer. Meth. in Aeronautical Fluid Dynam.— 1982.

3. Belotserkovsky О. М., Davydov Yu. M. Numerical approach for investigating some transonic flow//Lect Notes in Phys., Springer-Verlag.—1973. Vol. 19.

4. P u 11 i a m Т. H. Artificial dissipation models for the Euler equation//AIAA Paper.— 1985, N438.

5.Jameson A. Solution of the Euler equation for two-dimensional transonic flow by a multigrid method//Applied Math, and Comput—1983. Vol. 13, N 3—4.

6. С а к о в и ч B.C. Многосеточный алгоритм интегрирования уравнений Эйлера для расчета обтекания профиля трансзвуковым потоком/ЛГруды ЦАГИ.— 1991, №2517.

7. Jameson A. Vertex based multigrid algorithm for three-dimensional compressible flow calculation//Numerical methods for compressible flows.

New York.— 1986.

8. Sakovich V. S. Multiple-grid solution of Euler’s equations on unstructured grids//Comp. Math.& Math. Phys.— 1994. Vol. 34, N 12.

9. Vyshinsky V. V., Kravchenko S. A. On the possibility of mathematical modelling of separated flows by solving the Euler equations//Proceedings of the 19-th Congress of the International Council of the Aeronautical Sciences.— 1994. Vol. 3. 18—23 Sept., Anaheim, California, US.

10. Vyshinsky V. V., Kravchenko S. A. To calculation of threedimensional flows around the prolate bodies//RJCM.— 1994, N 2.

11. Владимирова H. А., Вышинский В. В., Кравченко С. А., Сакович B.C. Исследование вихревых структур в рамках краевой задачи для уравнений Эйлера//Труды ЦАГИ.— 1996, № 2622.

12. V1 a d i m i г о v а N. A., Sakovich V. S. Vortex flow simulation in the framework of a boundary-value problem for the Euler equations/ZTrudy TsAGI.— 1997. Vol. 2627.

13. Vladimirova N. A., Sakovich V. S. Vortex structures simulation j using the Euler equations//Proceedings of the Fifth Russian — Chinese Symposium

on Aerodynamics and Flight Dynamics.— 1997. Vol. 1.

14. McCormic B. W., Tangier J. L. and Sherrieb H. E. Structure of trailing vortices//Joumal of Aircraft.— 1969. Vol. 5, N 3.

15. Donaldson C. A brief review of the aircraft trailing vortex problem//ARAP.— 1971. Rep. N 155.

16. Воеводин А. В., Гайфуллин А. М., Захаров С. Б., Судаков Г. Г. Зональный метод расчета следа за летательным аппаратом//Труды ЦАГИ,— 1996, № 2622.

17. Betz A. Behaviour of vortex systems//Zeitschift fllr angewandte Mathemathik und Mechanik.— 1932. Vol. 12, N 3; also NASA TM 713.

18. Stuever R. A., Greene G. C. An analysis of relative wake-vortex hazards for typical transport aircraft//AIAA Paper.— 1994, N 810.

Рукопись поступила 20/11998 г.