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

7

МАТЕМАТИЧЕСКОЕ И КОМПЬЮТЕРНОЕ

МОДЕЛИРОВАНИЕ _MODELING AND SIMULATION

УДК 532.529

ПРИМЕНЕНИЕ И РЕАЛИЗАЦИЯ РАЗНОСТНЫХ СХЕМ ВЫСОКОЙ РАЗРЕШАЮЩЕЙ СПОСОБНОСТИ ДЛЯ РЕШЕНИЯ ЗАДАЧ ГАЗОВОЙ ДИНАМИКИ НА НЕСТРУКТУРИРОВАННЫХ СЕТКАХ

К.Н. Волков^

a Университет ИТМО, Санкт-Петербург, 197101, Российская Федерация b Университет Кингстона, Лондон, SW15 3DW, Великобритания, k.volkov@kingston.ac.uk

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

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

APPLICATION AND IMPLEMENTATION OF HIGH-RESOLUTION DIFFERENCE SCHEMES FOR SOLUTION OF GAS DYNAMICS PROBLEMS ON UNSTRUCTURED MESHES K.N. Volkova,b

a ITMO University, Saint Petersburg, 197101, Russian Federation b Kingston University, London, SW15 3DW, United Kingdom, k.volkov@kingston.ac.uk

Abstract. The paper deals with an approach to finite volume discretization of unsteady Navier-Stokes equations on unstructured meshes, and its advantages and development prospects are discussed. Features of inviscid and viscous flux discretization and temporal derivatives are considered. The advantages of the proposed approach include: the ability to operate on both structured and unstructured meshes; usage of high-order finite difference schemes in time and space; selection of median control volume for discretization of governing equations; application of expressions for calculation of the gradient and pseudo-laplasian making it possible to obtain more accurate results on highly stretched meshes in the boundary layer; writing of equations for the calculation of fluxes through the faces of interior and boundary control volumes in the same form, that simplify software implementation. This approach gives the possibility to implement a strategy of mesh adaptation taking into account the features of the certain flow and gives wide opportunities to parallelize computations. Possibilities of the developed approach are demonstrated on the example of the problem solution related to simulation of unsteady flows in the gas turbine engines.

Keywords: fluid dynamics, unstructured mesh, finite-difference scheme, aerofoil.

Введение

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

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

ков. Основной его недостаток состоит в достаточно сложной процедуре сшивки решений, полученных в различных подобластях [2].

Характерной особенностью неструктурированных сеток является произвольное расположение узлов сетки в физической области (отсутствуют выраженные сеточные направления, нет структуры сетки, подобной регулярным сеткам). Число ячеек, содержащих каждый конкретный узел, может изменяться от узла к узлу. Узлы сетки объединяются в многоугольники (двумерный случай) или в многогранники (трехмерный случай). Как правило, на плоскости используются треугольные и четырехугольные ячейки, а в пространстве - тетраэдры и призмы. Основное преимущество неструктурированных сеток перед регулярными состоит в большей гибкости при дискретизации физической области сложной формы, а также в возможности полной автоматизации их построения. Для неструктурированных сеток легче реализуются локальные сгущения и адаптация сетки в зависимости от поведения решения.

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

В настоящей работе разрабатывается подход к дискретизации нестационарных уравнений Навье-Стокса на неструктурированных сетках в рамках метода конечных объемов. Расчетная сетка строится при помощи одного из коммерческих пакетов, таких как Gambit или ICEM CFD. Разработанные программные средства используют трансляцию сетки из формата сеточного генератора в формат общедоступной библиотеки ADF Software Library (Advanced Data Format), которая является частью библиотеки CGNS (CFD General Notation System), разработанной для внутреннего использования в корпорации Boeing и получившей широкое распространение в NASA и компании McDonnel Douglas Aerospace. Возможности разработанного подхода демонстрируются на примере решения задачи, связанной с обтеканием решетки профилей.

Метод конечных объемов

В консервативных переменных уравнение, описывающее нестационарное трехмерное течение вязкого сжимаемого газа, записывается как

^ + V-F (n, Q, VQ) = H (Q, VQ). (1)

dt

Здесь Q(x,t), F(n, Q, VQ), H(Q, VQ) представляют собой вектор консервативных переменных в точке x в момент времени t, вектор потока через поверхность, ориентация которой задается внешней единичной нормалью n, и источниковый член соответственно. При моделировании турбулентных течений уравнение (1) дополняется уравнениями модели турбулентности, а вместо молекулярных коэффициентов переноса используются их эффективные значения.

Вектор потока расщепляется на невязкую (индекс I) и вязкую (индекс V) составляющие:

F (n, Q, VQ) = F1 (n, Q) + FV (n, Q, VQ), (2)

Введем вектор невязки

R(Q) = V-F(n, Q, VQ) -H(Q, VQ).

Уравнение (2) примет следующий вид:

ddQ+R(Q) = о. (3)

dt

Для дискретизации уравнений, записанных в виде (1), (2) или (3), используется метод конечных объемов на гибридной сетке.

Интегрируя уравнение (1) по контрольному объему V, с границей dV,, ориентация которой задается внешней единичной нормалью n = {nx,ny,nz}, и применяя теорему Гаусса-Остроградского, получим:

dt J QdО + <j) [F(n, Q, VQ) -(v4 • n) Q] dS = J H (Q, VQ)dО. (4)

V dV V

Преобразуем уравнение (4) к виду dQ

idr+R (Q) = о. (5)

dt

Вектор невязки в уравнении (5) находится из соотношения

R(Q) = V[F(n, Q, VQ) -(Vb n)Q] dS - JH(Q, VQ)dО |. (6)

V [dV V J

Под vb = {ub,vb,wb} понимается скорость перемещения границы dV, контрольного объема V,-.

Среднемедианный контрольный объем V, связанный с узлом /=1,..., N гибридной сетки, где N - число узлов, строится таким образом, что геометрические центры ячеек сетки с вершиной в узле / соединяются друг с другом через середины разделяющих их граней. Пример контрольного объема показан на рис. 1.

- сетка

..... контрольный объем, полученный

соединением центров соседних ячеек

_ среднемедианный контрольный

объем

Рис. 1. Среднемедианный контрольный объем

Весовые множители (площади граней) внутренних граней контрольного объема являются антисимметричными, As. = -As. для V/ е E¡, а весовые множители его граничных граней - симметричными, Asik = Aski для Vk е B¡ [1]. При этом имеет место следующее соотношение:

£ n-ASp. + £ n¡k Asft = 0. (7)

jeE¡ keB¡

Здесь E¡ - множество внутренних граней, связанных с узлом i; B¡ - множество граничных граней, связанных с узлом i; n. - внешняя единичная нормаль, задающая ориентацию грани (i, j); As. - площадь грани, соединяющей узлы i и j; nik - внешняя единичная нормаль к граничной грани (i, k); Asik - площадь граничной грани, соединяющей узлы i и k.

Расчет потоков

Интеграл от потока в соотношении (6) разделяется на два слагаемых, связанных с внутренними и граничными гранями контрольного объема:

<fF(n,Q,VQ)dS = £F(n.,Q,VQ) | ) n.As. + £F(na,Q,VQ) nikAs* . (8)

V je-E, x'+x/) keB,

С учетом (8) соотношение (6) примет следующий дискретный вид:

i' ^

s¡=VT

(9)

£ Рпу+ £ Рлпл- Н V

V кеБ1

Здесь Р'у - поток через внутреннюю грань (/ , у), ориентация которой задается внешней единичной нормалью п у Р' к - поток через граничную грань (, к), ориентация которой задается внешней единичной нормалью пк; Д5у и Алк - площади внутренней ((, у) и граничной ((, к) граней контрольного объема (рис. 2).

i j

1-1 -1

i

J

а б

Рис. 2. Внутренний (а) и граничный (б) контрольные объемы

Внутренние грани. Поток через внутреннюю грань ( , у) контрольного объема вычисляется в серединной точке грани,

x у =1 ( + x у ),

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

Fj = 2 (F + Fj ) n j ч.

Суммируя по всем внутренним граням (для Vj e E¡ ), получим

F=2 Х(+ ^ )п« ч •

2 уеД

С учетом соотношения (7) для замкнутого контрольного объема можно записать:

F = 2 Х( + ) ^ Ч - F X Пу ^ •

2 уеЕ, уеЕ,

Тогда

F = 1 X(F/ -F) Пу•

2 ]еЕ,

Учитывая, что Длу = -А^;,- для V/ е Е, вклад каждой грани (/, у) контрольного объема представляется в виде

F = 2 (( - F ) n j Ч- •

Граничные грани. Поток через граничную грань (, , к) контрольного объема рассчитывается в точке

L_ JT [i + {D + 2)8]к ]

2Б + 2

где ^ - размерность задачи; Ок - множество узлов в граничной ячейке к, 5кк - символ Кронекера. Вклад грани, лежащей на границе области, записывается отдельно. Учитывая соотношение (7), имеем

1 Г л F = 2 X ( + ^ ) П; Ч + FikП к Ал,к - Fiк П к Ал,к + X П у Ч

2 уеЕ, V уеЕ,

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

F = 2 Х(( - F ) n ] Ч •

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

Невязкие потоки

Из соотношений (6) и (9) имеем

1 1 ( ^

R!(Q) = -fF1 (n,Q)dS = V XFZП]AS] + XFÎnaAs

(10)

V уеЕ кеВ,

Для дискретизации невязких потоков в (10) используется модифицированный вариант схемы Ми8СЬ, которая представляет собой комбинацию центрированных конечных разностей 2-го и 4-го порядка, для переключения между которыми служит сглаживатель потока, построенный на основе характеристических переменных.

Разностная схема. Рассмотрим уравнение

dQ + AVQ = 0.

dt

Якобиан находится из соотношения

(11)

A =

dF1 dF' dF' dFl

dQ dQ dQ dQ Схема MUSCL для уравнения (11) записывается в виде

2 ¡V (Q] )+F (Q. )-| A |(Q +- Q-)].

F] = ] 2

(12)

Соотношение (12) представляет собой комбинацию центральной разностной производной 2-го порядка и диссипативного члена. Диссипативное слагаемое представляется в виде [5]

w (Q+-q -)=2 (1 -к) w [( 1 Q] +- Q] +1 q H 2 Q] - Q'+2 Q-)

(13)

где ке [0,1], а Qj■+, Qj, Q,-, Q^ соответствуют точкам х,+, ху, х,, х,_, лежащим на равном расстоянии друг от друга. Соотношения (12) и (13) дают следующую разностную схему для расчета потока:

К = 2 {К (Qj)+^ ©)]-2 (-к) Ы (((Q) - А (Q))}, (14)

+

+

где L(Q) - псевдолапласиан. Для сохранения монотонности решения в схему (14) вводится ограничитель потока [6]. После линеаризации потоков получим

fi =1 \[AA + aQ MAI 3 (-ф)(( (Q) - L*(Q) - Q')

(15)

где L"(Q) - модифицированный псевдолапласиан. Полагается, что к=1/3. Соотношение (15) представляет собой комбинацию конечных разностей 2-го и 4-го порядка точности. Для переключения между ними служит функция (сглаживатель потока) [6, 7].

На твердой стенке, вследствие условия непротекания, = 0 для У к е Б1. На входной границе расчетной области полагается

К = 1 [ Рк (Qk)+Рк (Q») -| А| ( -Q»)],

где Qoo - консервативные переменные на бесконечности.

Расчет псевдолапласиана. Псевдолапласиан представляет собой обобщение центрированной разностной производной 2-го порядка на неструктурированную сетку:

Ц (Q)=щ Z( - Q■ )

(1б)

jeE.

где \Е,\ - число элементов множества Е. С использованием (16) оценки показывают, что

Ь (Q) ~ 0(к 2)v2Q\x=х_,

поэтому диссипативный член в (14) имеет порядок 0(НЪ). После подстановки в соотношение (15) и деления на объем в (9) получается погрешность порядка 0(Н2).

Разложение в ряд Тейлора в окрестности точки x ■ дает ь ш)=Ь (x)VQ \x

где Ь) = {Ь ¡х,Ь ¡у,Ь, х}'. Следовательно, на неравномерной сетке схема не имеет 2-го порядка точности. Если Q является линейной функцией пространственных координат, то Ь(0) Ф 0. Указанные обстоятельства приводят к потере точности решения, в связи с чем псевдолапласиан переопределяется следующим образом [9]:

L"(Q) = Ь (Q) (x). (17)

Если Q = ax+c, то Ь(@) = aLI■(x), VQ = а , поэтому Ц^) = 0 . Вместе с тем, представление псевдолапласиана в виде (17) допускает потерю устойчивости численного решения и дает неточные результаты на сильно растянутых сетках, используемых для расчета течения в пограничном слое. Для обеспечения устойчивости решения вводится оператор масштабирования:

( V1 ^ ^ ( .у1

L,(Q) =

Z

1

VjeE |xj -xi\

Z Qj-Qt

jeE,

xj -x¡

L¡ (x) =

Z-

jeE,

1

xj -x¡

Z

jeE,

xj -x¡

xj -x¡

После этого модифицированный псевдолапласиан находится из соотношения L* (Q) = Ц (Q)- VQ,L, (x). (18)

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

Расчет градиента. Рассмотрим контрольный объем abcde, связанный с узлом i неструктурированной сетки (рис. 3). С учетом того, что Qi = const и VQ = 0 , используя формулы Грина и теорему Стокса, получим

Л 1 ( Л

VQ = 1I\VQdю -\VQdю = 11 (Qdl-(Qtdl

^ ¿ \ r^ I ^ ¿ \ r T

(19)

уо о У " V ь Ь У

Здесь ^ - расширенный контрольный объем 12345, участвующий в вычислении градиента; Ь - площадь граней расширенного контрольного объема; V - исходный контрольный объем аЬсёв; Б - площадь граней контрольного объема. Интегралы, входящие в (19), представляются в виде

( Qd> = Z а + Q

n,.

( Qidl = QL .

L 1,2,3,4,5 ^ l

При этом имеют место следующие соотношения:

S = Sab + Sic + Scrf + SÄ + Sa;

L — Lj2 + L23 + L34 + L45 + L51 — 3 х S . Q —3х-.

4

Рис. 3. Вычисление градиента Нормали к граням контрольного объема находятся из соотношений

Используя формулы Грина и учитывая (7), получим

(VQ ), —

X 2 (Q] - q ) n]Ч

]eE. 2

(20)

Вязкие потоки

Из соотношений (6) и (9) имеем

R-(Q) — - f F-(n, Ö, VQ)dS — - (X F.4 AS] +X F- nit As.

—' — l jeE,

л

_ , (21)

1 д¥, ' 1 V уеЕ кеВ

Учитывая, что вследствие условий прилипания и непротекания на стенке = 0 для Ук е В1, и пренебрегая вязкими силами на входной границе расчетной области, соотношение (21) можно переписать в следующем виде:

< = V X Е; п , Ал,.

У1 /ед

Соотношение (22) используется для всех контрольных объемов, включая граничные. После линеаризации соотношения (22) получим

R- — 1

' -

X

dF" -v— n As

d(VQj ) ] j

dF-

d(VQ )

■ — BM

-1 5Q dl

(22)

(23)

где М - матрица перехода от консервативных переменных к примитивным, В - вязкий якобиан. Соотношение (23) приобретает вид

r- — 1XQl-Q BM 1 n..As ..

]eE |X] - Xi

Для расчета градиента VQ в серединной точке каждой грани х, = (х ,+х;)/2 в соотношении (23) используется полусумма соответствующих узловых значений. Для расчета градиентов (VQ), и (VQ)j в

узлах сетки применяется соотношение (20). Однако среднее арифметическое центральных разностей не демпфирует высокочастотных гармоник решения [8, 9]. Хотя выражение для расчета невязких потоков и включает диссипативные слагаемые, демпфирующие высокочастотные осцилляции решения, этого недостаточно в пограничном слое, где вязкие члены становятся доминирующими. Исходя из этого, составляющая градиента в направлении наиболее короткой грани заменяется простыми разностями:

VQ] — (VQj )* -

(VQ] )* -Ч -

Q. - Q

X . - X .

] '

8s .. —,

X f - X..

(24)

X . - X .

] .

Дискретизация по времени

Перепишем уравнение (3) в виде

dQ

dt

— L(Q) :

(25)

3

5

где L(Q) - дифференциальный оператор. После линеаризации уравнение (25) примет вид

dQ = CQ, dt

где C - квадратная матрица.

Для дискретизации уравнения (25) используется k-шаговый метод Рунге-Кутта

Q("+1) = A(kC)Q(n\ Л(z) = iamzm ,

m=0

где a0 = a1 = 1 и ap Ф 0. Область устойчивости S = {z = x + iy: |Л(z)| < 1} имеет вид окружности с радиусом rc. На комплексной плоскости область устойчивости представляется в виде круга [7]: z = r exp (i 6).

Радиус области устойчивости находится из соотношения

л ^ 3л г = min r (6), — <6<—. c 6 2 2

Приведем расчетные соотношения пятишагового метода Рунге-Кутты (rc = 2,7):

Q,(m) = Q,(0) -«mAtR(m-1) (m = 1,...,5),

где

R(m-1) = C. (Q(m-1)) — B(m-1) •

B(m-1) = ßmD (Q}m-1) ) + (1 -ßm ) B(m-2).

При этом C (Q(m-1)) представляет собой вклад конвективных слагаемых, Ц (Q,<m-1)) учитывает

вклад источниковых членов, а также физической и численной диссипации. Коэффициенты cm и ßm имеют следующие значения:

113 1,

а1 =—, а2 = —, а3 = —, а, = —, а5 = 1;

1 4 2 6 3 8 4 2 5

14 11

ß = 1, ß2 = 0, ß3 = — ß, = 0, ß5 = — .

Поскольку ß2 = 0 и ß4 = 0, то Ц (Q(2) ) и Ц (Qf} )

не рассчитываются.

Ускорение сходимости

Для ускорения сходимости используется многосеточный метод решения системы разностных уравнений. Для построения последовательности вложенных неструктурированных сеток используется метод схлопывающихся граней [9]. Два узла 1 и; неструктурированной сетки, связанных гранью, заменяются одним узлом, расположенным посередине между ними. Схлопывание ячейки производится в направлении наиболее короткой грани. Градиент находится из соотношения (24).

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

Результаты расчетов

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

Профиль лопатки представляет собой модифицированный профиль МЛСЛ0006 [12], который строится при помощи суперпозиции распределений кривизны и толщины профиля вдоль координаты 0 < х / Ь < 1, где Ь - хорда профиля. Распределения кривизны и толщины профиля даются соотношения-

C ( x) — Hc - R +[ R2 - ( x - 0,5)2 J"2; T(x) — Ht (2,969x1/2 -1,26x - 3,51x2 - 2,843x3 -1,036x4 ).

ми

Здесь Я = {н] + 0,25) / 2Ис. Величины И1 и Ис представляют собой максимальную толщину профиля и радиус кривизны в серединной точке хорды профиля соответственно. В расчетах принимается, что Н = 0,06, Нс = 0,05 (рис. 4).

Рис. 4. Профиль лопатки турбины

Геометрия расчетной области показана на рис. 5. Угол расположения лопатки относительно оси х равняется в = 45°. Поток поступает через входное сечение под углом а к горизонтальной оси (угол атаки профиля равняется а - в). В качестве рабочей среды принимается воздух.

L/2 L/2

Рис. 5. Геометрия расчетной области

Рассматриваются два варианта течения в межлопаточном канале, соответствующие различным граничным условиям во входном сечении. Во входном сечении задается направление течения и число Маха (M = 0,7, а = 55° в случае 1 и M = 0,8, а = 58° в случае 2). Приведенные значения чисел Маха соответствуют перепадам давлений pi/p0 = 0,8716 в случае 1 и p1/p0 = 0,8740 в случае 2. Перепад давления в выходном сечении расчетной области определяется при помощи решения соответствующей стационарной задачи (течение в канале с неподвижными лопатками) при нескольких заданных значениях числа Маха на входе. Такая процедура требуется потому, что в [11] давление в выходном сечении не указывается (приводится только число Маха на входе). Полученный перепад давлений задается в качестве граничного условия в выходном сечении расчетной области. На направлении оси z выставляются периодические граничные условия. На поверхности профиля используется граничное условие непротекания для нормальной составляющей скорости.

Неструктурированная сетка в серединном сечении канала содержит 5443 треугольных элемента (5 вложенных сеток, сетка наилучшей разрешающей способности содержит 115 узлов на поверхности профиля, V-цикл). Сгущение узлов сетки производится у поверхности профиля. Шаг интегрирования по времени полагается равным At = 3,25-10-5 с.

Вибрация лопаток вопроизводится в виде двух последовательных циклов - восходящего и нисходящего движения профиля по гармоническому закону. Восходящее движение профиля происходит с амплитудой 0,01 м по нормали к его хорде. Амплитуда нисходящего движения составляет 2° относительно серединной линии профиля. Рассматривается 3 частоты (ю = 0,25, 0,5, 0,75 рад/с) и 2 фазовых угла (ф = 0°, 90°). Производится расчет одного периода колебаний T = 2п/ю. В качестве начального приближения используется стационарное решение задачи.

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

x(t) = x0 + sin (cot)(a+8x+ + a~8x~).

Амплитуды возмущений вычисляются следующим образом:

a+ = 211 + sgn[sin(at)]}, a~ = 2{1-sgn[sin(at)]} .

Здесь x0 - начальное невозмущенное положение узлов; 5x+ и 5x- - максимальные отклонения узлов от невозмущенного положения.

Новые координаты граничных узлов находятся из соотношения

L

x"+' (x, t) = x0 (xn) + Real [ [dxr (xn) + idxt (xn)] exp (iat)J, причем

x = x0 + dxr cos (at)-dx¿ sin (at), = -ro[dxr sin (at) + dxi cos (at)].

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

Распределения коэффициента давления по поверхности профиля показаны на рис. 6 и рис. 7 для случаев 1 и 2 при двух положениях профиля. Коэффициент давления определяется следующим образом:

C = Pd - Ри

Р тт2\ I ,

pU aro

где pU и pD представляют собой давления на верхней и нижней поверхности профиля.

а б

Рис. 6. Распределения коэффициента давления по поверхности профиля при восходящем (а) и нисходящем (б) движении в случае 1. Значки • соответствуют данным [13]

0 0,2 0,4 0,6 0

х

0 0,2 0,4 0,6 0,:

х

б

Рис. 7. Распределения коэффициента давления по поверхности профиля при восходящем (а) и нисходящем (б) движении в случае 2. Значки • соответствуют данным [13]

Заключение

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

а

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

Литература

1. Волков К.Н., Емельянов В.Н. Вычислительные технологии в задачах механики жидкости и газа. М.: ФИЗМАТЛИТ, 2013. 468 с.

2. Cagnone J.S., Sermeus K., Nadarajah S.K., Laurendeau E. Implicit multigrid schemes for challenging aerodynamic simulations on block-structured grids // Computers and Fluids. 2011. V. 44. N 1. P. 314-327.

3. Deng X., Mao M., Tu G., Zhang H., Zhang Y. High-order and high accurate CFD methods and their applications for complex grid problems // Communications in Computational Physics. 2012. V. 11. N 4. P. 10811102.

4. Wang Z.J., Fidkowski K., Abgrall R., Bassi F., Caraeni D., Cary A., Deconinck H., Hartmann R., Hillewaert K., Huynh H.T., Kroll N., May G., Persson P.-O., van Leer B., Visbal M. High-order CFD methods: current status and perspective // International Journal for Numerical Methods in Fluids. 2012. V. 72. N 8. P. 811845.

5. Luo H., Baum J.D., Lohner R. Edge-based finite element scheme for the Euler equations // AIAA Journal. 1994. V. 32. N 6. P. 1183-1190.

6. Crumpton P.I., Moinier P., Giles M.B. An unstructured algorithm for high Reynolds number flows on highly stretched grids // Proc. 10th International Conference on Numerical Methods in Laminar and Turbulent Flow. Swansea, UK, 1997. P. 561-572.

7. Brognies Z., Rajasekharan A., Farhat C. Provably stable and time-accurate extensions of Runge-Kutta schemes for CFD computations on moving grids // International Journal for Numerical Methods in Fluids. 2012. V. 69. N 7. P. 1249-1270.

8. Moinier P., Giles M.B. Stability analysis of preconditioned approximations of the Euler equations on unstructured meshes // Journal of Computational Physics. 2002. V. 178. N 2. P. 498-519.

9. Crumpton P.I., Giles M.B. Implicit time accurate solutions on unstructured dynamic grids // International Journal for Numerical Methods in Fluids. 1997. V. 25. N 11. P. 1285-1300.

10. Brandt A. Multi-level adaptive solutions to boundary-value problems // Mathematics of Computation. 1977. V. 31. N 138. P. 333-390.

11. Fransson T.H., Verdon J.M. Standard configurations for unsteady flow through vibrating axial-flow turbo-machine-cascades // Unsteady Aerodynamics, Aeroacoustics and Aeroelasticity of Turbomachines and Propellers. Ed. H.M. Atassi. NY: Springer-Verlay, 1993. P. 859-889.

12. Verdon J.M. Linearized unsteady aerodynamic theory // United Technologies Research Center Report R85-151774-1. 1987.

13. Lawrence C., Spyropoulos E., Reddy T.S.R. Unsteady cascade aerodynamic response using a multiphysics simulation code // NASA Report TM-2000-209635. 2000.

Волков Константин Николаевич - доктор физико-математических наук, научный сотрудник, Универ-

ситет ИТМО, Санкт-Петербург, 197101, Российская Федерация; старший лектор, Университет Кингстона, Лондон, SW15 3DW, Великобритания, k.volkov@kingston.ac.uk

Konstantin N. Volkov - D.Sc., Researcher, ITMO University, Saint Petersburg, 197101, Russian

Federation; Senior Lecturer, Kingston University, London, SW15 3DW, United Kingdom, k.volkov@kingston.ac.uk

Принято к печати 14.05.14 Accepted 14.05.14