книги / Основы проектирования турбин авиадвигаделей

Для замыкания задачи обтекания различными авторами используются раз­ нообразные условия для задания течения в окрестности выходной кромки. В численных методах, применяемых в отечественном турбостроении, чаще всего задается условие, являющееся аналогом условия Кутта —Жуковско­ го для скругленных кромок конечной толщины. Как правило, оно сводит­ ся к положению о равенстве скоростей в точке отрыва потока на спинке

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

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

рi cosj31 = р2 V2cosj32

ивторого из уравнений (5.9). После очевидных преобразований получаем

уравнения

решая которые любым итерационным методом (например, методом Ньютона —Рафсона) получаем интересующие нас значения р2 и V2.

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

5.3. ПРИМЕНЕНИЕ МЕТОДА КОНЕЧНЫХ ЭЛЕМЕНТОВ

ДЛЯ РЕШЕНИЯ ЗАДАЧИ

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

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

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

172

2.Получение матриц элементов.

3.Объединение матриц элементов в общую матрицу и построение век­ тора узловой нагрузки.

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

5.Решение полученной системы уравнений.

6.Расчет любой функции, зависящей от узловых неизвестных.

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

Итак, линеаризация проводится относительно местного значения чис­ ла М. Вводя обычным образом задаваемую на элементе локальную прямо­ угольную систему координат, воспользуемся подходом, аналогичным изло­ женному в уже упоминавшейся работе [18]. Направим ось £э вдоль векто­

ра скорости на элементе Уэ , а ось rj3 —перпендикулярно к оси £э. Обозна­ чая через Мэ и рэ значения числа Маха и плотности на элементе, рассмотрим решение первого из уравнений (5.9) как задачу минимизации функционала

/(Ф) = Б / (э) (Ф),

,5

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

Гэп+ l) = ( V 3 + u ' , v ' ) ,

при обычном условии

173

V

« 1 ;

Ф = ^ э + Ф' « э >ЧЭ);

р = р э + р'-

Подставляя выражения для Ф и р в исходное уравнение, отбрасывая члены второго порядка малости и учитывая, что

ЭКз

= 0, получим

Первый член уравнения (5.10) можно определить из уравнения движения (5.2)

и пренебрегая величинами второго порядка, получим соотношение

подставляя которое в уравнение (5.10) приходим к уравнению Прандтля Глауэрта для возмущения потенциала

построим функционал, соответствующий уравнению (5.11),

174

•Ъп

где 5э и 1э —соответственно площадь и граница элемента в системе коор­ динат | £э' , Т7э } . Представляя уравнение (5.12) в виде Ы = Ы х + 5 / 2, рассмотрим отдельно каждое из слагаемых. Возвращаясь к исходным пере­

Преобразуя выражение для 5/2, обозначая Фэ = Гэ£э и зная, что и' « Уэ, получим

Для общего функционала с учетом (5.13) и (5.14) получим

1 6/ = — 5 2 р эЯ [ ( 1 - М э2)Ф | +

2 1

Сдп

где С —граница области течения, так как сумма криволинейных интегра­ лов, взятых вдоль общих границ соседних элементов, равна нулю. Задача нахождения стационарной точки функционала (5.15) решается итерацион­ ным способом, причем на каждой итерации вычисляются новые значения Уэ, Мэ, р э и направление £э. Дня дискретизации расчетной области в целях более точного описания геометрии ее границ целесообразно использовать криволинейные элементы. Для решения поставленной задачи привлекаются четырехугольные изопараметрические элементы, построение которых осу­ ществляется по восьми узловым точкам. Для каждого элемента в локаль­ ной системе координат (рис. 5.2) имеем

Ф (х ,у) = БФ tNi(E ,Z ), i = 1

где Nf — функции формы элемента; Ф,- —величина потенциала в узловых точках. Функции формы Nj в локальной системе координат { Е, Z^ при нумерации узлов по часовой стрелке, начиная с верхнего правого узла, имеют вид

Рис. 5.2. Изопараметрический конечный элемент

176

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

* = I х (Щ Е , г ) 1

/= 1

у = l yiNt {E9Z).

1 = 1

Далее задача сводится к минимизации функционала (5.15) относительно узловых значений Ф/ (/ = 1 , 2,..., 8 ) на каждом элементе. Рассмотрим пер­ вое слагаемое функционала и выразим все производные, входящие в него, через функции формы

Из условия минимума функционала бД с учетом (5.16) получаем выраже­ ние для элемента матрицы локальной минимизации

и для элемента правой части уравнения

// = - Р э ^ э м э Я ^ p - d ^ d i ,э . 5, Э*э

Производные функций формы в декартовых и локальных координатах свя­ заны соотношением

177

Разрешая записанное матричное уравнение относительно декартовых производных, получим

где у —угол наклона вектора скорости в локальной системе координат. Ве­ личина этого угла определяется для центральной точки каждого элемента на предыдущей итерации по очевидной зависимости

Вычисление коэффициентов Ку для каждого элемента, как это обычно де­ лается в методе конечных элементов, проводится по квадратурной форму­ ле Гаусса на сетке 3X3 и полученные коэффициенты объединяются в гло­ бальную матрицу.

Рассматривая функционал 6 /2, замечаем, что в результате интегрирова­ ния по контуру интегралы вдоль всех линий, кроме ЕТ (см. рис. 5.1), рав­ ны нулю или взаимно уничтожаются. Тогда

где Ьэ —длина границы каждого элемента, принадлежащего ЕТ. Проводя минимизацию относительно Фь Ф2, Ф3, получаем дополнительный вектор правой части полного уравнения минимизации

178

/ 1 4 1

О

О Г *

о

о

о/

Врезультате имеем систему алгебраических уравнений

[К] | ф j = { F }

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

5.4. РАЗБИЕНИЕ РАСЧЕТНОЙ ОБЛАСТИ НА ЭЛЕМЕНТЫ

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

179

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

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

и

(1 + / 2) 3/ 2

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

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

180