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

в потоке сжимаемого газа

Воронич И.В. (voronich@pochta.ru)

Факультет аэромеханики и летательной техники МФТИ

1. Постановка задачи

Шум летательных аппаратов является предметом активного исследования в последние годы. Повышение плотности полетов и влияние шума самолетов и вертолетов на экологию вблизи аэропортов ведут к ужесточению норм, ограничивающих их характеристики. Уменьшение шума летательных аппаратов -проблема комплексная и включает борьбу с шумом от различных источников: винтов, крыльев, планера самолета, двигателей, струй двигателей. Для уменьшения уровня шума необходимо знание исходного акустического поля. В настоящее время разрабатываются как аналитические, так и численные подходы к решению таких задач [1, 2].

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

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

Задача взаимодействия вихревых структур с профилем исследуется около 20 лет как с помощью экспериментальных средств, так и с помощью численного моделирования [3, 4]. Преимущественное внимание уделялось при этом лопастям вертолетов в силу большого практического значения такой задачи. Несмотря на использование экспериментальных методов, полученные данные позволяют судить об акустическом поле только частично. Это связано с объемом и сложностью измерений и проблемами генерации вихревых структур с заданными параметрами. Расчетные исследования вертолетных конфигураций осложняются пространственным характером, необходимостью учета вращения лопастей и динамики вертолета. Представление об основных характеристиках акустического поля и силового воздействия на профиль можно получить на основе решения упрощенных задач.

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

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

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

Представляется целесообразным применять низкодиссипативную базовую схему, основанную на характеристическом подходе, в совокупности с методами повышения точности 2-го порядка, сохраняющими монотонность (TVD) и методами Рунге-Кутта или предиктор-корректор интегрирования по времени. Для вычисления ближнего акустического поля такие методы в совокупности с достаточно подробными расчетными сетками могут дать приемлемые результаты. Для решения задачи распространения волн на большие расстояния можно использовать схемы высокого порядка аппроксимации [7].

Рис. 1. Схема обтекания аэродинамического профиля при наличии вихря в потоке

В работе рассматривается обтекание простого симметричного 12% профиля, состоящего из двух дуг окружности. Параметры задачи следующие: идеальный невязкий сжимаемый газ с показателем адиабаты у = 1.4, числа Маха набегающего потока берутся равными Moo = 0.5, 0.65, 0.8, угол атаки а = 0°, вихревая структура является плоским изоэнтропическим вихрем с диаметром (по полю скорости), близким к хорде профиля, и давлением и плотностью в центре, составляющими около 0.9 от значений в потоке. Направление вращения в вихре против часовой стрелки. Общий вид задачи изображен на рис. 1. Так как число Маха является фактором, который оказывает наибольшее влияние на взаимодействие и характеристики акустического поля [3, 4], представляется необходимым исследовать это влияние в первую очередь. Максимальное число Маха течения, создаваемого покоящимся вихрем, составляет 0.27, так что вихрь создает хотя и небольшое, но ощутимое возмущение.

2. Законы сохранения и численный метод

Рассмотрим законы сохранения массы, импульса и энергии в неподвижном контрольном объеме Q газа в интегральной форме:

J QdQ + J HdS = 0.

здесь

Q = (р, pu, pv, pe)T , E = (pu, pu2 + p, puv, puh)T , F = (pv, puv, pv2 + p, pvh)T , H = Enx + Fny, n = (nx, ny )T .

Здесь Q - вектор консервативных переменных, H - вектор потока, E и F -векторы потока в декартовых координатах, n - единичный вектор нормали к поверхности. Здесь h = e + p/p - полная удельная энтальпия. Будем рассматривать совершенный газ, так что полная удельная эн ергия есть e = (u2 + v2)/2 + p/(p(y - 1)).

В случае дифференцируемости Q, E, F справедливы уравнения Эйлера в дифференциальной дивергентной форме:

д Q д E dF п д t д x д y

В силу свойств векторов потока справедливы соотношения E = AQ, F = BQ , где A и B

- якобианы: A = дEдQ, B = дF/дQ. Поэтому возможна квазилинейная форма уравнений Эйлера:

дQ ЛдQ пд Q

Для вектора потока H имеем: C = дН/дQ = nxA + nyB . Матрица С может быть преобразована к диагональной форме:

.2 i,.2>

Un - c

v Un + c J

где Un = nx u + ny v - нормальная компонента вектора скорости, Ut = -ny u + nx v -касательная компонента, векторы n и т образуют правую систему, ЛС - диагональная

матрица, c = /yp/p - скорость звука, ТС состоит из правых собственных векторов С.

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

Одномерные уравнения в инвариантах Римана можно получить из одномерных уравнений Эйлера для дифференцируемых полей течения:

д s д s л - + u- = 0.

д t д x

+ (u ± а)-- = 0, (5)

д t д x

где R± = u ± 2а/(у - 1), s = p/p1 - инварианты Римана и энтропийная функция соответственно. Инварианты Римана переносятся со скоростями u ± а (акустические характеристики). Энтропия переносится со скоростью u (энтропийная характеристика). Уравнения для инвариантов Римана наглядно отражают характеристические свойства уравнений Эйлера в случае гладкого поля течения. На практике они используются для постановки волновых граничных условий на внешних границах расчетной области, где течение является дозвуковым.

Рассмотрим систему законов сохранения в интегральной форме (1). В качестве контрольного объема выберем четырехугольную ячейку расчетной сетки с площадью AV, полагая ее стороны отрезками длины а/. Зачастую значение интеграла по контуру ячейки представляется суммой произведений аппроксимаций вектора потока в центрах граней Ф на длины граней а/. Таким образом, получается полудискретная задача:

d j QdQ + £ф/°/ = 0. (6)

Введя понятие осредненной по объему ячейки величины AVQn = jQdQ , уравнение (6)

п

можно переписать как дифференциальное уравнение по времени для осредненной величины:

dt AV i

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

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

- \Qt , x < 0

lQR, x > 0

где L и R обозначают состояния слева и справа от грани. Необходимо найти решение при t > 0. Такая задача также называется задачей о распаде разрыва. В методе С.К. Годунова [8] аппроксимация вектора потока через грань вычисляется на основе решения этой задачи. При условии гиперболичности системы (2) автомодельное решение (зависящее от x/t) задачи с начальными данными (8) состоит из m центрированных волн в плоскости (x, t), разделяющих m+1 кусочно-постоянное состояние (включая два начальных).

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

распределение всех параметров: Q(x, t, QL, QR) при -да < x < да, t > 0. Информация, которая необходима для построения аппроксимации вектора потока в методе Годунова, сводится к использованию состояния, образовавшегося после распада разрыва на грани расчетной ячейки (t > 0):

Ф ,+1/2 = E(Q (x+1/2, t, Ql , Qr )). (9)

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

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

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

+A(Ql , Qr ) Q = о. (10)

д t д x

Матрица A (QL, QR ) должна удовлетворять следующим условиям:

1. VQL, Qr AE = A (Ql , Qr )AQ (здесь и далее Af = fR - fL)

2. A (Q, Q) = A(Q)

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

В силу диагонализуемости A(QL, QR) может быть представлена в виде

A(QL, QR ) = TAЛArTA1, где соответствующие матрицы имеют тот же смысл, что и в (4). Учитывая структуру матриц и свойство 1 соответствующее приращение вектора потока можно записать в виде разложения:

m

Здесь а - вектор характеристических приращений: а = TAQ, rk и Xк - правые

собственные векторы и собственные числа A(QL, QR ) .

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

Ф;+1/2 = El + е акХЛ = ER -Е акХЛ ,

Фi+1/2 = 2(El + Er)--еак . (12)

В результате построения матрицы Роу можно найти, что A (QL, QR) = A (u, v, /г), где

и, v, li - средние по Роу величины [9]:

~=4p~lul+Vp~Rur v~=4pLVL+4pRvR r=4pLhL+4pRhR (13)

u/ = I- /- , 1; = I- I- , 1 = I- I- . (13)

VPL wPR VPL wPR VPL wPR

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

1 1 m

Ф = 2(Hl + Hr)-2еак|ч|Гк . (14)

Здесь собственные векторы гк и собственные значения Хк вычисляются на основе

осредненных переменных (13), а = rTC1AQ .

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

Проблема получения нефизичных обобщенных решений связана с условием энтропии [10]. Это условие предполагает выполнение дискретного аналога

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

В качестве практического решения указанной проблемы был предложен такой способ ограничения снизу модулей собственных значений в (14), чтобы не допускать обращения их в нуль, особенно для волн разрежения, при сохранении гладкости аппроксимации вектора потока [10]. Соответствующее выражение для скорректированных модулей собственных значений имеет вид:

v 5k J

(15)

Здесь вводятся величины 8k, вообще разные для каждого характеристического поля. По виду зависимостей A,/c(u/c) можно судить о схемной диссипации разностного метода. По сравнению со схемами расщепления вектора потока метод расщепления разности вектора потока обладает заметно меньшей схемной диссипацией, что обеспечивает лучшее разрешение контактных разрывов и сдвиговых слоев.

Не любые значения 5k обеспечивают необходимую коррекцию. Задачи с сильными волнами разрежения предъявляют особые требования к схемной вязкости при наличии звуковой точки. Как показывает практика, для контактного разрыва и волны сжатия можно обойтись без внесения дополнительной диссипации, тогда как для волны разрежения со звуковой точкой такая мера необходима. Нужное количество диссипации определяется исходя из перепада значений соответствующей характеристической скорости: 5k = 2 max(AXk, 0). В случае растяжения характеристического поля (AXk > 0) вблизи звуковой точки минимальное значение Xk определяется по перепаду характеристической скорости AXk, тогда как при сжатии

(AXk < 0) минимальное значение не отграничивается от нуля. Такая коррекция обеспечивает хорошие характеристики схемного решения в случае волн разрежения и позволяет получать физичные обобщенные решения задач газовой динамики [11].

Повышение точности метода по пространственным переменным и по времени в формулировке схем, сохраняющих монотонность (или TVD), предполагает вычисление уточненных значений переменных слева и справа от грани расчетной ячейки и применение методов Рунге-Кутта или предиктор-корректор интегрирования по времени [10, 11, 12].

В настоящей работе использовались следующие выражения для уточненных значений примитивных переменных (p, u, v,p) на грани расчетной ячейки [12, 13], вычисленные по принципу наименьших приращений:

gL = gi + min m0d(PL (gi - gi-1 X PR (gi+1 - gi )) ,

PL = di-1/2,i /(di-1,i-1/2 + di-1/2,i X PR = di,i+1/2 /(di,i+1/2 + di+1/2,i+1)

gR = gi+1 - min m0d(PL (gi+1 - gi X PR (gi+2 - gi+1 )) , (16)

PL = di+1/2,i+1 /(di,i+1/2 + di+1/2,i+1) , PR = di 3/2 /(di+1,i+3/2 + di+3/2,i+2) ,

здесь di, k обозначает расстояние между соответствующими точками расчетной сетки. Такой же подход применяется для величин в индексном направлении j. На функцию minmod(x, y) налагается ряд требований, обеспечивающих повышение порядка

аппроксимации до второго для гладких решений и сохранение монотонности решения (при монотонных начальных данных для скалярного закона сохранения) [14]. Типичные функции minmod(x, y) должны быть близки к базовой функции

minmod(x, y) = (sgn(x) + sgn(y))min( x , y )/2, (17)

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

minmod(x,y) = (sgn(x) + sgn(y))min(a x , a y , ( x + y )/2)/2. (18)

Такая функция обеспечивает гладкость переключения производных при x = y, в то же время не слишком отличаясь от базовой функции (непрерывно переходя в нее при a - 1), надежные значения a можно брать в интервале 1.05-1.25.

Для интегрирования полудискретной задачи (7) по времени для явных схем используются методы Рунге-Кутта или методы типа предиктор-корректор [15]. В настоящей работе использовался метод второго порядка точности

Q = Q + 2 AtLA ({Q}),

= Qj +AtLA ({Q*}). (19)

Здесь LA - соответствующий конечно-разностный оператор. Такой метод обладает небольшой схемной диссипацией, что в сумме с описанным методом вычисления потоков дает численный алгоритм, пригодный для моделирования задач со сдвиговыми слоями и контактными разрывами (для задач с сильными скачками предпочтительнее использовать многошаговые методы Рунге-Кутта [15]). В частности, расчет используемой для оценки диссипативных и дисперсионных ошибок численных методов задачи о переносе вихря показывает, что при переносе вихря на расстояние десяти диаметров численные ошибки незначительны, форма вихря сохраняется очень хорошо, метод уверенно показывает второй порядок сходимости к точному решению, в том числе на хаотически возмущенных сетках.

Разностные граничные условия являются реализацией условия непротекания на твердых границах

Un\w = 0 (20)

с учетом кривизны стенки (при использовании уравнения (3.2), условий постоянства энтропии и полной энтальпии) [16] и разностную реализацию уравнений (5) для нахождения параметров на внешних дозвуковых границах либо экстраполяцию параметров на внешних сверхзвуковых границах.

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

3. Условия расчета и результаты

Результаты были получены для режимов с числом Маха M = 0.5, 0.65, 0.8 при угле атаки а = 0° и показателе адиабаты у = 1.4. При этом использовалась вихревая структура со следующим распределением окружной скорости в системе центра вихря:

V (r) = 2-R exp ((1 - (r / Я)2)/2). (1)

Здесь Г - интенсивность вихря, R - радиус ядра вихря. Для расчетов брались величины Г = 2, R = 0.25 (при этом давление и плотность в центре составляют около 0.9 от невозмущенных значений, а диаметр вихря близок к единице). Распределения плотности и давления получаются при использовании (1) и уравнения импульса в полярной системе координат

£ = Р (2)

r r

при дополнительном энергетическом условии постоянства энтропии в вихре: p / pY = s0 Таким образом, распределения плотности и скорости даются выражениями:

2 v(y-1)

1 --8 exp (1 - (r / R)2)

(& 8v) = (-y, x). (4)

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

f±(x) = ±(b + ylR2 -(x-0.5)2), 0 < x <1 (5)

R = d/2 + 1/(8d), b = d/2- 1/(8d)

где d - полутолщина профиля.

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

Первый и второй невозмущенные режимы течения являются дозвуковыми с ощутимым влиянием сжимаемости. Третий режим трансзвуковой с замкнутыми изолированными сверхзвуковыми зонами на верхней и нижней поверхностях профиля. Для каждого режима приводятся по четыре картины изолиний давления на последовательных этапах взаимодействия (рис. 2, 4, 6). Также для всех режимов приведены зависимости от времени (в единицах L/cx) аэродинамических коэффициентов Cx, Cy, и Mz, позволяющие получить представление о силовом воздействии на профиль (рис. 3, 5, 7).

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

На основе анализа результатов (рис. 2, 4, 6) можно отметить, что при возрастании числа Маха увеличивается интенсивность взаимодействия и уменьшается длина волны акустических возмущений. Появляются инициированные структуры в ближнем поле течения в виде волн сжатия и разрежения, которые существуют еще ощутимое время (~L/UX.) после того, как акустические возмущения покинут зону взаимодействия. Характерным является различие структуры акустических волн, идущих от разных сторон профиля: разрежение-сжатие-разрежение (волны I типа) и сжатие-разрежение-сжатие (волны II типа), при этом волны II типа имеют более отчетливую форму.

Анализ полученных данных позволяет определить базовые частотные и амплитудные характеристики волн. Значения амплитуд AI, AII, длин волн AI, AII и частот f ближнего акустического поля представлены в таблице 1. По этим данным видно, что, несмотря на различие форм, волны I и II типов отличаются амплитудами, но мало отличаются по длине. Длина волны варьируется вдоль ее фронта, несколько возрастая в направлении от носка профиля к задней кромке, и уменьшается с ростом числа Маха. Амплитуда волн растет с ростом числа Маха, при этом амплитуда волн I типа больше, чем у волн II типа. Направление распространения имеет большое значение: максимальную амплитуду имеет часть волны, идущая в приблизительно поперечном к течению направлении.

Диапазон частот волн можно определить исходя из скоростей распространения разных частей волн и их длин. Для оценки принималось, что в системе координат, связанной с потоком, эти скорости близки к скорости звука 340 м/с при 18°С. Хорда профиля бралась равной 0.2 м. Как видно из таблицы, имеет место интенсивный рост частоты волн при приближении числа Маха к единице. Это обусловлено проявлением нелинейных эффектов сжатия с ростом скорости потока.

Таблица 1

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