CC BY
18
Математические структуры и моделирование 2003, вып. 11, с. 67-87
УДК 532.5.031+519.632.4
ОБТЕКАНИЕ РЕШЕТОК ПРОИЗВОЛЬНЫХ ЛОПАСТЕЙ ИДЕАЛЬНОЙ НЕСЖИМАЕМОЙ
ЖИДКОСТЬЮ
А.С. Толстуха
Various aspects of the numerical solution using finite element method for the 3-d potential flow of the ideal fluid through the rotor blades of the turbomachine are discussed. Paper covered the problem stage, wakes influence and secondary flow questions. Also considered discreet solution techniques and its convergence.
1. Введение
Развитие турбостроительной техники выдвигает повышенные требования к знаниям гидродинамических характеристик рабочих колес и направляющих аппаратов, что в свою очередь стимулирует развитие их гидродинамической теории. Современная вычислительная техника позволяет описать сложные гидродинамические процессы в проточной части турбомашин в пространственной постановке соответствующих задач. Однако реализация наиболее полной модели течения, описываемой уравнениями Навье-Стокса, требует слишком больших затрат машинного времени, поэтому в инженерных расчетах они неприемлемы. Результаты многочисленных исследований показывают, что ряд проблем, возникающих при проектировании решеток, с достаточной степенью точности могут быть решены и с помощью более простых моделей. Остается лишь вопрос о пределах их применимости. К числу таких проблем относятся аэроупругие явления в решетках, нестационарные аэродинамические характеристики которых при безотрывном обтекании с достаточной степенью точности определяются в рамках модели идеальной жидкости. Как известно [1], нестационарные аэродинамические характеристики решеток зависят от стационарного поля скоростей, расчету которого и посвящена настоящая работа.
Работы по методам определения стационарного пространственного течения в турбомашинах принято вести от работы [2] - квазитрёхмерная технология, расчёты проводятся последовательно для течений в «тонких» слоях. Осесимметричное в меридианальном слое, приближённо являющемся средней поверх© 2003 А.С. Толстуха
E-mail: ast@univer.omsk.su
Омский государственный университет
Работа выполнена в рамках программы ФЦП «Интеграция»
68
А.С. Толстуха. Обтекание решеток произвольных лопастей ...
ностью тока. Затем в построенном по осесимметричному течению цилиндрическом «тонком» слое, расчёт в котором позволяет учесть форму профилей и, наконец, расчёт в поперечном каналу слое, позволяющий определить вторичные течения. Решение в каждом из трёх слоев уточняет решение в двух других. Применяется в самых различных вариантах и содержательных контекстах, например [3] - для нескольких ступеней со стыковкой условий (склейкой) на границах ступеней; [4] - осесимметричный в меридиональной плоскости МКЭ;
[5] - слой переменной толщины на осесимметрической поверхности.
Много работ посвящено определению течения в поперечных каналу плоскостях - так называемые вторичные течения. Этому вопросу неизменно уделяется внимание в экспериментальных и теоретических исследованиях пространственной структуры потока, [6] - монография, [7] - расчёт в направляющем аппарате, [8] - эксперимент.
Возникновение вторичных течений принято рассматривать в первую очередь как результат вязкого взаимодействия жидкости или газа с твёрдыми поверхностями межлопаточного канала - поверхностями лопастей и цилиндрических обводов. Действительно, к примеру, задача о натекании пограничного слоя с внутреннего обвода на угол в месте сопряжения с поверхностью лопасти в окрестности задней кромки [9] может быть рассмотрена только в рамках модели вязкой среды. Но вклад поверхности лопастей в образование вторичных течений можно учесть в рамках модели идеальной жидкости или газа. Действительно, наблюдаемая картина обтекания лопастных систем со сходом потока с задних кромок, возможная лишь в силу вязких процессов на поверхности лопасти, в отсутствии вязкости реализуется применением условия Жуковекого-Кутта на задней кромке; при этом вокруг лопасти возникает в общем случае переменная по высоте лопасти циркуляция скорости, а с кромки сходят свободные вихри, образующие вихревую пелену. Возникающая картина течения включает и перетекание в поперечных потоку плоскостях.
Таким образом, без учёта влияния свободных вихрей и положения пелены нельзя получить правильную картину пространственного обтекания. В предположении идеальности жидкости вихревая пелена стационарного течения есть поверхность контактного разрыва вектора скорости. Если предположить, что вне поверхности пелены вихри отсутствуют и энтропия постоянна, то есть течение потенциально, тогда в силу ненулевой циркуляции скорости вокруг лопасти потенциал терпит разрыв на пелене. В [10] педложена процедура учёта течения за решёткой лопастей: линейная комбинация решений с единичным разрывом потенциала скоростей на уровнях дискретизации по высоте. Очевидно, что применение его возможно для области с предписанным положением пелены или же требуется итерационная процедура уточнения, на каждом шаге которого требуется решить I + 1 задачу (I - число уровней дискретизации по высоте), что нерационально. Дело осложняется тем, что вихревая пелена неустойчива [11].
В данной работе предложен итерационный процесс, на каждом шаге которого в ограниченной трёхмерной области решается задача относительно потенциала скоростей, коэффициенты и правая часть которой определены предыдущим приближением решения, разрыва потенциала, положения вихревой пелены. Но-
Математические структуры и моделирование. 2003. Вып. 11.
69
вое приближение позволяет скорректировать положение поверхности пелены, исходя из условия контактного разрыва скорости, и скорректировать величины разрыва потенциала вдоль задней кромки, исходя из условия Жуковекого-Купа. Частная задача решается приближённо методом конечных элементов. В свою очередь, использование итерационного метода для решения возникающей при этом системы линейных алгебраических уравнений позволяет эффективно использовать уже полученное приближение.
2. Область и основные уравнения
Рассмотрим нестационарное движение идеальной однородной несжимаемой жидкости. Уравнения сохранения массы и импульса имеют вид
div и = О
d 1
-fU + -VP at р
О,
u(t,x) - вектор скорости, P(t,x) - гидростатическое давление, р - плотность постоянная во всей области течения, х — <пп ”• -ЛТ
д д д
dlV дх + ду + dz,V ^ \дх’ду’дг) ’ dt dt
x,y,z
д д д\т d д
\ , — = —+ u-V
декартова координата, полная производная.
Пусть система вращается с угловой скоростью £ вокруг оси х, £ = (£, О, 0)т -вектор вращения. Обозначим через w относительную скорость. Абсолютная скорость и есть w + £ х х. Абсолютное ускорение есть сумма относительного, переносного (в данном случае центростремительного) и поворотного кориолисова ускорения.
4w + £ х (£ х х) + 2(£ х w) : 1'
dt
—VP.
Р
Уравнение сохранения импульса в форме Громеки-Лэмба
w х rotw + £ х (£ х х) + 2(£ х w)
-VP.
Р
Заметим, что
С X (С X х) = (О, -СУ. -Єр = -V Q (Cr)2) = -V Q'U2) .
где г = \fy2 + z2, ив = £r - окружная скорость вращения системы. Так как
rot w = rot (и — £ х х) = ш — 2£,
получим
dw
ж+у
Ur
W
2
-|—VP = w х ш, Р
70
А. С. Толстуха. Обтекание решеток произвольных лопастей ...
Рис. 1. Область межлопаточного канала.
Рис. 2. Решётка винтовых лопастей. Сплошной линией показаны результаты работы [12].
и - завихренность в неподвижной системе координат.
Предположим, что завихренность в абсолютном движении отсутствует, а в относительном движении постоянна и обусловлена переходом в неинерциальную систему координат. Положим
w(t, х) = УФ(і, х) - ( х х.
Подставим это выражение в предыдущее тождество и получим аналог интеграла Лагранжа
дФ
dt
+
УФ - С х х
— U
в р - + - = !>(().
Уравнение неразрывности примет вид
ДФ = 0,
д2 д2 д2
А = —— + —— + —— - оператор Лапласа. ох2 ду2 oz1
Рассмотрим установившееся течение однородной несжимаемой жидкости через вращающийся осевой венец турбомашины. Перейдём в связанную с венцом декартову систему координат (x,y,z). Осью венца положим ось х. В этой системе рассмотрим область Q (1), ограниченную спинкой Г і и корытцем Г2 двух соседних лопастей венца турбомашины, сходящими с задних кромок соседних лопастей поверхностями вихревой пелены за венцом £1? £2, периодическими границами перед венцом Li, L2, примыкающими к передним кромкам соседних лопастей, а также поверхностями внутреннего (tip) и внешнего (hub) обводов.
Предположим, что течение в неподвижной системе координат является безвихревым. После перехода во вращающуюся с угловой скоростью ( вокруг оси х систему - течение вихревое. Предположим малым вклад составляющей скорости, индуцированный возмущением поля завихренности лопастями, т.е. считаем
Математические структуры и моделирование. 2003. Вып. 11.
71
Рис. 3. Распределение скоростей стационарного течения на поверхности лопасти.
Рис. 4. Зависимость коэффициента Сп от угла атаки.
Рис. 5. Вторичное течение
72
А.С. Толстуха. Обтекание решеток произвольных лопастей ...
Рис. 6. Распределение скоростей на поверхности лопасти кольцевой (N = 8) решётки телесных (Т = 10%) профилей (линии). Точками - скорости на поверхности профиля
плоской решётки.
во всей области вихрь постоянным и равным —2£, В условиях этих предположений потенциал не зависит от времени Ф(хД) = фо(х.) и удовлетворяет однородному уравнению
Аф0(х) = 0, х Є И (1)
Давление р(х) определяется интегралом Бернулли
wg - U20 Р У'І ~ ЦДоо Роо 2 + р 2 + р ’
где wo = Vo(t-(xl Vx - вектор скорости, poo - давление далеко перед венцом в точке, лежащей на той же линии тока. Заметим, что для осевых венцов можно считать отклонения линий тока по радиусу малыми, тогда интеграл Бернулли примет привычный вид
wg р _ У'4 />.х
2 р 2 р
Итак, предполагаем, что далеко перед венцом имеется невозмущенный поток скорости U x = (// х і. u0о2, и0оз)т іакой. что |Uoo| = 1 на среднем радиусе, // х і = const ( ось х направим по оси венца), и002 и и0оз постоянны на каждом радиусе и удовлетворяют условию потенциального вихря. На бесконечности за венцом ставится условие, обеспечивающее сохранение расхода. Например, если радиусы цилиндрических обводов ПОСТОЯННЫ ПО X, МОЖНО ПОЛОЖИТЬ U+O01 = Wool-
Математические структуры и моделирование. 2003. Вып. 11.
73
Условия на других границах следующие: дфо
дп
х
хх • п = О,
х Є Гх U Г2 U ф U hub,
(3)
где п - внешняя нормаль, £ - вектор вращения венца вокруг своей оси,
Фо(х) = A(QX) + Соо, х е А, (4)
дфо,
х
дф0
(Qx), х Є Li, Qx Є Ь2.
(5)
9n v 9n
Константа C00 представляет тангенциальную составляющую скорости перед венцом вида г и = (\/ (aR), где a - угол между соседними лопастями, R -расстояние до оси венца.
0о(х) = фо (Qx) + С0(т}(х)), X Є А, (6)
^(х)- ^(xxj 'П = 0, хеАиА, (7)
г] - естественная координата вдоль задней кромки, Со(г]) - скачок потенциала скорости при переходе через вихревую пелену, являющуюся поверхностью контактного разрыва, при этом величина Со(г]) — С00 равна циркуляции скорости вокруг лопасти по замкнутому контуру, имеющему точку на задней кромке с координатой г]. В свою очередь, г](х) обозначает координату частицы жидкости на задней кромке, затем перенесенную потоком в точку х Є А- Когда х лежит на задней кромке, условия (6), (7) являются одной из форм условия Жуковекого-Кутта, Изменение Со(г}) вдоль задней кромки даёт интенсивность свободных вихрей на пелене. Заметим, что в общем случае функция Co(rf) и расположение поверхностей А, А неизвестны и подлежат определению.
Решение сформулированной краевой задачи осуществляется итерационным методом, В качестве первого приближения используется решение для течения несжимаемой жидкости с (произвольно) заданной функцией Со(г]) и (произвольно) выбранным положением поверхности А- Поверхность А ПРИ этом задаётся поворотом А на угол а вокруг оси венца. На очередном шаге итерационного процесса текущее приближение фо(х) используется ДЛЯ коррекции Со(г]), выполнения условия Жуковекого-Кутта на задней кромке и изменения границ области А, А, таким образом, чтобы они являлись поверхностями контактного разрыва, т.е. на них выполнялось условие (7),
Частные задачи решаются приближённо методом конечных элементов в ограниченной области О' С О, содержащей межлопаточный канал и части О вверх и вниз по потоку, ограниченными плоскостями (! х. G+0о, перпендикулярными оси венца. Условия на этих границах поставим следующие
дфо
дп
х)
^ООІ;
X Є G
—ос
дфо
дп
х
u+ool; X Є G_|
74
А.С. Толстуха. Обтекание решеток произвольных лопастей ...
Затем О' разбивается на элементы - тетраэдры, вершины тетраэдров называются узлами разбиения.
Следуя методу конечных элементов, свяжем е каждым узлом разбиения х/ базисную кусочно-линейную функцию V/(x). На каждом тетраэдре VJ(x) определяется как линейная функция, причём
Щхз)
1,3 = 1
о, j ф I.
Эта функция финитна, носителем её является объединение тетраэдров е вершиной X;.
Приближение потенциала ищем в виде
N
Фо(х) = ^ФоїЩх), (8)
і=і
где фоі - узловые значения приближения на разбиении О'.
В слабой интегральной форме уравнение (1) можно записать в виде
J(Уф)(х) • W/(x)) dn = J (V<Mx) • пЩх) dS. (9)
П' дП>
Подставим (8) в (9), получим
[ (Уф,(х) • n)Ty(x) dS ап>
или
N
^ ^ aij Фої — bi,
3=1
J(VVj(x) • VV/(x)) dQ, bl = f (V<^o(x) • пЩх) dS.
n> an>
С учётом сказанного можно записать:
2J2J ФоіЯцт — k
Тзі зет
(10)
где
gljT = J(VVj(x) • Wi(x)) dT = (VVj\T ■ VVt\t) VolT
- вклад в элемент матрицы оу тетраэдра разбиения Т. Вклад тетраэдра в симметричный относительно диагонали элемент матрицы по величине совпадает с дцт, что обеспечивает симметричность матрицы. VoIt - объём Т.
Математические структуры и моделирование. 2003. Вып. 11.
75
Для удовлетворения условию периодичности полученная система редуцируется отождествлением соответствующих узлов на периодических границах. При этом в некоторых уравнениях редуцированной системы появляется правая часть, зависящая от СД и дискретного представления Со(г]).
Пусть /|. /•_> номера соответственных узлов на периодической границе. Для определенности х/, Є Сі, х;2 Є £•_>. Имеем два уравнения вида (10) (индекс Т в коэффициентах g опущен).
/ \
E Ущ Фоівід + ^ ФокШік = / ((C X x) •n)^i(x) dS, ас
ТЭЦ jeT кет \jeCi HCi J dnl
( \
E УУ Фоз9і23 + Фок9кк = [ ((C X x) •n)^2(x) dS, (12)
ТЭЦ jeT кет \3eC2 HC2 / J dnl
Здесь правые части уравнений определены условием (7), Покажем, что они равны по абсолютной величине и противоположны по знаку.
Во-первых, п|хЄ£2 = — Qn х. с, и ( = QQ (вектор вращения направлен вдоль оси венца); учитывая, что Q - матрица вращения (Q-1 = QT), получим
((С х х) • n)|xe£l = ((С х х) • Q-1Qn)|xe£l =
= ((С х Qx) • Qn)|xe£l = - ((С х х) • п)|хе£2.
Предположение о том, что лопасти имеют заостренные задние кромки, обеспечивает выполнение этих тождеств и для узлов на задних кромках.
Во-вторых, разбиения поверхностей производятся таким образом, что одно переходит в другое при повороте вокруг оси венца, а значения функций ІуДх), Vi2 (х) в любом тетраэдре из носителя определяется только взаимным расположением вершин тетраэдра и не зависит от его расположения относительно оси координат, В частности, значение в точке на грани зависит только от взаимного расположения вершин грани. Из сказанного следует
V)i(x)|xe£i = V)2(Qx)|Qxe£2.
Рассмотренные тождества обеспечивают равенство с точностью до знака интегралов в правых частях (11), (12),
Пусть узлы разбиения занумерованы так, что, если х;- Є С\. то хj+Np Є £2, Xj+Aтр = Qx.j, Np - количество узлов на периодической границе. Тогда, заменяя в уравнении (11) значения потенциала в узлах разбиения Сі по условию (6) на значения в узлах разбиения £•_> плюс разрыв
Ф$з = Фод+мр + CQj Xj Є Сі
76
А.С. Толстуха. Обтекание решеток произвольных лопастей ...
и складывая его с уравнением (12), получим уравнение, не содержащее узловых значений потенциала на С\
(
\ (
\
Е
ТЭЦ
Е °п,/ • \';,.'//и' + Е Фок9кк
зет
\ЗЄС 1
feeT
fe^£i
У
+ Е
Тэт
Е Фоз9і23 + Е ФокЗкк
зет
\jeC-2
кет
к<£С2
У
Е Е с'оійхі- (із)
ТЭЦ je'T jeCt
Итак, два уравнения (11) и (12) можно заменить на одно (13). Уравнения с І ф Сі дают:
/
\
Е
тэг
Е Фо,з+тгР9р + Е
зет \зеС 1
feeT
У
[ ((( х х) ■ n)VJ(x) d5
СІП'
^ЕЕ Го,'//;. (14)
тэг je'T jeCi
Подобным образом удаляем уравнения с номерами 1\ х; Є Li, используя условия (4), (5). Легко видеть, что редуцированная система остаётся симметричной.
Следующий вопрос - о разрешимости полученной системы линейных алгебраических уравнений. Рассмотрим элементы матрицы в строке с номером I. Вклад в диагональный элемент матрицы тетраэдра Т Э I
дат — (VV/1 т
• W/|t) > 0.
С другой стороны, учитывая тождества
Еадіт = і, Evv^ = 0
jeT jeT
которые следуют из определения функций I ;(х). для суммы вкладов в другие элементы строки матрицы тетраэдра Т получим цепочку
/ \
Е дз1т
E(VV)|t • W/|т)
E VVj'|t
• Viz |t —
jeT jeT
зфі зфі
зет
ЧІЛ
У
— — (VI'zIt • VVz|t) — —9ііт■
Из чего следует, что сумма элементов матрицы в строке равна нулю. Матрица вырождена, так как вектор, состоящий из единиц, является собственным вектором, отвечающим нулевому собственному значению.
Математические структуры и моделирование. 2003. Вып. 11.
77
Вырожденноеть матрицы системы без условий закрепления - проявление того, что потенциал скоростей для краевой задачи Неймана определен с точностью до константы. Способ получить решение - найти его в подпространстве, ортогональном собственному вектору системы, который отвечает нулевому собственному числу. Нетривиальное решение существует в силу того, что необходимое условие разрешимости системы с вырожденной матрицей - ортогональность правой части исходной системы решениям однородной системы с сопряжённой матрицей - в рассматриваемой задаче эквивалентно дискретному аналогу условия сохранения массы газа при протекании через область
О = / <livVo0 dx = / (Vo(t • и) OS
П'
dfl'
Для решения линейной системы алгебраических уравнений Ay = b с симметричной матрицей А, которая имеет собственный вектор е = (1,,,,, 1)т, отвечающий нулевому собственному числу, а правая часть b ортогональна е, воспользуемся итерационным методом, например, вида
I(yW^y(fc-1)) + Ay(fe-1) =b.
т
Тогда легко увидеть, что, выбрав вектор начального приближения у(°) ортогональным вектору е, на следующем шаге получим вектор yW с таким же свойством, Если процесс сходящийся, то предел последовательности ортогонален е и является решением системы.
Коррекция Cqi производится по уравнению (11) для узлов на задней кромке X; Є £іПГі_ Как видим, правая часть уравнения пропорциональна перетеканию через границу, именно
[ ((С X х) • п)1Дх) dS ~ ((С X X/) • п).
дП>
Используем (11) как формулу для расчёта правой части по текущему приближению о0(х). Полученная величина характеризует величину перетекания через £i в точках задней кромки. Условие на задней кромке считаем выполненным, если дефект
rf- = EE ФоіЯііт - ((С х х/) • n)Si/S
ТЭ1 ІЄТ
удовлетворяет условию
\di\ < е (15)
для всех узлов на задней кромке, где е - заданная точность, У, - площадь треугольников, образующих разбиение границы и содержащих вершину х/. Внешнюю нормаль п в рассматриваемой точке следует понимать как среднюю нормаль всех примыкающих треугольников
n kSr,
КЭ1
S(=y:sk
КЭ1
78
А.С. Толстуха. Обтекание решеток произвольных лопастей ...
Если (15) не выполнятся - подправим значение Со/. Дефект скорости по нормали к пелене в узле на задней кромке равен
Un = 3^-. (16)
Ь/
Средняя скорость в окружном направлении щ = Coi/(aRi), Щ - расстояние от х/ до оси венца. Проекция скорости на нормаль есть ue(tr п), где t/ - единичный вектор, который, будучи помещен В точку X;, касается окружности с центром на оси венца, окружность лежит в нормальной оси венца плоскости. Учитывая вышесказанное, можно положить
Gif’ = С<рц + hS,
aR[
h > 0 - релаксационный параметр.
Возможность определения величин перетекания в узлах на пелене за кромкой позволяет скорректировать положение пелены. Рассмотрим узел х/ Є £\. Выше по потоку на пелене имеется узел Щ. В Х/ скорость вдоль поверхности пелены обозначим Нс - она вычисляется по значениям потенциала в точках х/. х/,. в предположении, что указанные точки лежат вдоль линии тока на пелене. Нормальная к пелене скорость un вычисляется по формуле (16). Идея заключается в том, чтобы сместить х/ в х( так, чтобы тангенс угла между векторами х/ — х/,. 11 xj — х/,. равнялся її,,/ |іі£і, то есть совместить направление вектора на пелене х( — xfe с направлением вектора скорости + unп.
Построим х( поворотом точки х/ вокруг оси венца. Тогда геометрические рассмотрения приводят к выражению для дуги окружности между этими двумя точками.
_ 'U'n |^Z ^Cfc I
|u?| sin /З + un cos/З’
где COS /3 = (ti ■ n).
Таким образом, при коррекции пелены точку х/ нужно повернуть на угол a/Ri вокруг оси венца.
3. Расчётные аспекты
Итак, для получения численного решения о стационарном течении через венец лопастей в ограниченной по размаху области требуется итерационная процедура. От шага к шагу изменяется расчётная область, что отражается в величинах элементов матрицы и правой части линейной системы в силу изменения носителей базисных функций I /(х). Изменяются величины разрыва потенциала на пелене (правая часть).
Первое замечание касается выбора метода решения линейной системы. Её саму следует решать итерационным методом. Основной аргумент - итераций для решения системы линейных алгебраических уравнений требуется тем меньше,
Математические структуры и моделирование. 2003. Вып. 11.
79
чем точнее задано первое приближение: таковым естественно взять полученный последним результат «внешних» итераций. Если итерационный процесс решения задачи сходится, то количество «внутренних» итераций для частной задачи должно падать до одной итерации на решение линейной системы, что является критерием останова «внешнего» процесса. Стационарные методы решения линейных систем не используют такого рода «априорную» информацию: время их работы фиксировано и зависит от порядка системы и ширины (полу)ленты. Другим аргументом служат меньшие затраты машинной памяти. Итерационный метод требует хранения только ненулевых элементов матрицы; стационарный - (порядок системы х ширина (полу)ленты) элементов. Второй сомножитель зависит от нумерации узлов разбиения области. Применяются различные процедуры перенумерации, предваряющие формирование матрицы, уменьшающие ширину (полу)ленты, Пространственность и наличие периодических границ в расчётной области, когда разнесенные в пространстве узлы взаимодействуют, ведут к большим потерям памяти.
Второе замечание относится к расчёту разрыва и положения пелены. Расчёты показали, что величины разрыва потенциала, определяемые через условие Жуковекого-Кутта на задней кромке, сходятся быстро, характер сходимости монотонный. Вычисленное так распределение циркуляции по размаху лопасти определяет основной вклад в интегральные характеристики лопасти. Следующим по значимости вклад вносится учётом распределения углов схода потока путём подстройки положения пелены. Конечно, в процессе подстройки меняются и величины разрывов потенциала. Влияние на венец дальнейшей эволюции пелены вниз по потоку очень быстро (на порядки) падает.
Как уже отмечалось, пелена неустойчива относительно малых возмущений, следовательно, задача определения её положения и формы является некорректной, Требуется какого-либо рода регуляризация процесса подстройки положения пелены. Нами использован принцип монотонности изменения положения узлов по высоте на каждом шаге. Потеря точности при этом составляет порядок по сравнению с точностью на твёрдых стенках: если на стенках (и внутри области) уравнение вида (11) выполняется с точностью 1(П5, то на пелене она на порядок ниже. Характер сходимости процесса с подстройкой немонотонный и требует в несколько раз большего количества итераций, чем процесс подбора величин разрыва потенциала. Удобна следующая стратегия организации итерационного процесса: вдоль направления потока выбирать ряд узлов по размаху на пелене с максимальным перетеканием и для этого ряда проводить необходимые изменения. Как правило, в начале процесса максимальное перетекание наблюдается на кромке (выстраивается распределение циркуляции), затем в первом ряду за кромкой (выстраивается распределение углов схода потока) и так далее.
Рассмотрим в качестве примера несжимаемое течение через пространственный осевой венец, лопасти которого представляют части винтовых поверхностей таких, ЧТО течение СО скоростью U = 11о + ш х х через венец является однородным, Uo направлен вдоль оси венца. Лопасть представляет собой прямоугольник {г0 < г < Г|. —ф + an < в < ф + an} в координатах (г, в), где в -
80
А.С. Толстуха. Обтекание решеток произвольных лопастей ...
угол в тангенциальном направлении, a - угол между соседними лопастями, п -порядковый номер лопасти. Вдоль оси венца лопасть ограничена координатами 6
х = - (гі — го) etg/З. Мы брали ф = 0,19635, [3 = тг/6. г0 = 1, г і = 2, число
Z __
лопастей - N = 8, 12, 16, При этом густота на среднем радиусе составляла т = 1, 0; 1, 5; 2, 0. На рис.2 приведены результаты расчётов по методу конечных элементов в сравнении е результатами работы [12], в которой применён метод дискретных вихрей. Здесь для расчётов угол атаки задавался равным 7г/6 на среднем радиусе.
На рие.З представлены распределения скоростей вдоль хорды для разных радиусов по высоте (скорости на задней кромке убывают е ростом радиуса).
На рис.4 представлен график изменения коэффициента нормальной силы в зависимости от угла атаки на среднем радиусе. Угол атаки менялся вследствие изменения угловой скорости вращения. Как видно из графика, все точки лежат на прямой Cn = ( ')] a.
Анализируя этот же рисунок, можно заключить, что в данном случае влияние формы вихревой пелены на интегральные характеристики лопасти невелико - коэффициент нагрузки линейно зависит от угла атаки. Тем не менее и здесь можно проследить влияние вихревой пелены (вторичное течение). Рассчитаем положение вихревой пелены и заменим периодические граничные условия за венцом на условие непротекания (стенка). Решив задачу, получим некое течение в системе каналов заданной формы. Затем рассмотрим разность двух решений. Она будет представлять течение, обусловленное стационарной вихревой системой. На рис.5 иллюстрировано вторичное течение в трёх последовательных перпендикулярных оси венца сечениях. Первое содержит задние кромки, последующие расположены ниже по течению на расстоянии 5% и 10% хорды. Слева показаны поля направления скоростей в плоскости сечений, справа -изолинии осевой компоненты скорости. Порядок скорости у задней кромки есть единица (рис.3), поэтому можно заключить, что разница в скоростях двух течений составляет в данном случае рз 0, 5%.
Рассмотрим течение в решётке телесных профилей. На рис.6 приведены распределения скоростей на лопасти при стационарном обтекании неподвижного осевого венца. Количество лопастей N = 8. Относительная толщина профиля Т = 0,1 от величины хорды. Густота решётки на среднем радиусе т = 1. Высота лопастей равна величине хорды. Угол установки /3 = 60°. (1) - скорости на втулке, (2) - на среднем радиусе, (3) - на внешнем обводе. Для сравнения точками показано распределение скоростей в плоской решётке. Как видно из рисунка, имеет место неравномерность распределения поля стационарного течения по высоте лопасти.
4. Сходимость
Здесь рассмотрен вопрос о сходимости приближённого решения трёхмерной задачи Неймана для оператора Лапласа.
Известно (См., например, [13] гл. IV), что (обобщённое) решение данной за-
Математические структуры и моделирование. 2003. Вып. 11.
81
дачи существует и единственно в энергетическом пространстве определённого линейного оператора. Это пространство вкладывается в пространство (локально) интегрируемых функций, поэтому его можно аппроксимировать последовательностью конечномерных подпространств, базисными функциями которых могут быть, в частности, кусочно-линейные. Приближения к решению в каждом из конечномерных подпространств образуют последовательность Ритца, сходящуюся в энергетической норме к решению. Метод Бубнова-Галёркина оказывается методом построения элемента последовательности Ритца,
Представляет интерес вопрос о сходимости в терминах пространств функций, определённых на сетках дискретизации расчётной области. Действовать будем в духе теоремы сходимости [14]. Плюсы: во-первых, требуется изучать свойства разностного оператора, действующих на конечномерных пространствах сеточных функций, во-вторых, имеем в виду относительно простое обобщение на нелинейные задачи (например, течение газа), так как при наличии у разностного оператора аппроксимационных свойств условие устойчивости рассматривается для линеаризованного уравнения.
Рассмотрим разностный оператор относительно фі,..., фм с заданными значениями нормальной производной на границе фп{ = и ,и
У Y, ■ VIj)Volr - У] 6’д =0, і € {1.....N).
ТЭг jeT А Зі
Подставим в разностный оператор значения точного решения в {щ}Дь взятые в узлах разбиения {хД^Д области О, щ = Цх*), (если точка неособая и щ конечно)
,'' = 'E'E".iiVV'-vv.i)Vo1
ТЭг j£T
1
~UTl
Азі
Второе слагаемое есть
• f (Vu(xi) • п)^ dS.
an>
По узловым значениям щ строим кусочно-линейный интерполянт и/(х). Тогда первое слагаемое есть
('Чщ • VVi) d,Q
Am Vi dn + J" (Vu/ • n) Vi dS.
O'
f)T
Интегралы по границам тетраэдров возникают в силу разрывности градиента кусочно-линейной uj. Vuj кусочно-постоянен, Auj = 0,
(Vur • n)Vi dS^ f (Vu(xi) • n)Vi dS.
an>
82
А.С. Толстуха. Обтекание решеток произвольных лопастей ...
Разложив функцию Цх) в ряд в точке х' до второго порядка, учитывая, что V VVj; = 0, получим
ІЄТ
Vmj(x') = V ( ^Цх5ЦЦх') | = ^2u(xj)VVj =
\j£T ) jeT
X fM(x') + (Vm(x0 • (xi - x')) + \(Dj(xj - x') • (xi - x'))) VK
AezT1 \ /
= £(va(x') ■ x,)W, + 1 £(D,(X, - x') ■ (x, - x'))VV
jeT jeT
где Dj - матрица e компонентами {dik}iy=1,2,3
d2u
31
cl
Ik
dx[dxk
(x' + 9j(x.j — x')), 9j Є [0,1].
В точке x' = X,- Є T
1
Vu/fe) = ^(Уи(х4) • Xj)VVj + - ^2(Dj(xj - х^) • (x, - Xi))VV).
jeT jeT
Заметим, что для любого постоянного вектора а Є R3 имеет место тождество
а = Vfa • х) = V а ■
X хз¥з(х)) = Х(а • xi)Wi-
jeT ) jeT
Таким образом,
Vur(xt) = Vu(xt) + - - хі) • (х^ - Xi))VVj.
jeT
Вычислим вклад разрывов производных при переходе через грани тетраэдров, при этом учтём, что, во-первых, Цх) непрерывно дифференцируема (при переходе через грань скачок градиента равен нулю), и, во-вторых, для любого постоянного вектора а Є R3 поток через границу тетраэдра равен нулю
( £ (a - nA)SA = 0).
дет
и/ • п)Т) dS = Ur ■
S,
пд
д
Т QJ,
Тзі Азі
дП'
X X ( ( VM(X0 + \ ^iDAxi ~ xi) • (xi “ X«))V^- ] • Пд j
Тзі Азі V V ^ jeT / /
(Vu(Xi) • n)Vi dS^Yl \\\ J2(DAX3 - Хг) ■ (xi - Хг))^У.
тзі \\ jeT
s,
Д
St
Пі I —,
Математические структуры и моделирование. 2003. Вып. 11.
83
где Si - площадь грани, лежащей напротив вершины е номером г, п* - её внешняя нормаль.
Завершим равномерную оценку
ш <
еМ
Тэг \ j£T
Dj ||2 |xj
x,l2|vijl Т<
1
< - sup II І) 2 x.. о
К
тзі
1
< - sup 2 x.. о
D
hr
E
T'3i
U)iTa,n
dm-Ax - наибольшая длина ребра тетраэдра разбиения области, hm\n - наименьшая высота тетраэдра разбиения области, Шіт - телесный угол, под которым из
вершины е номером і видна противоположная грань тетраэдра Т. ^ wіт Е 4-7г,
тзі
Окончательно
Ф
ОО
max \фі\ < 2л sup *Є{1....JV} x.. о
D
hn
Для ошибки Z{ = фі — щ схема имеет вид
Y,-./(VI} • WjWobr = ~Фі ге{1,...,Щ.
тзі jeT
(17)
В операторном виде
A{N)Z = -Ф.
Линейный оператор A^N\ рассматриваемый в области определения RAM, где е = (1,..., 1)т Є Rv. является обратимым. Покажем сначала, что матрица А(Ю
неотрицательно определена, то есть для любого £ Є R3, £ Ф О выполнено (ЛМ£ • А > 0, Действительно,
N
N
(Т'Щ ■ Є) = Е Е Е ■ Vb)vЫт = у> Y,(Wi ■ £T)voiT
i=1
тзі j£T
і=1
тзі
= ^ ^ £i(VVi • Фг)Уоіт = У^фІт ' &)VolT > 0,
Т і (ф_Т т
Sr = YljcT СіУУіт ~ градиент кусочно-линейной функции, заданной узловыми значениями Е, ' Є {1,,,,, Л }•. Чтобы сделать вывод о невырожденности линейного оператора ЛМ, определённого на подпространстве RJV\{e}, достаточно заметить, что сели у вектора £ не все компоненты одинаковы, то в разбиении области существуют тетраэдры, для которых і, / Ф 0. Следовательно, нестрогое неравенство нужно заменить строгим на RJV\{e}, оператор - положительно определён.
Рассмотрим поведение ошибки в равномерной норме.
5
І
< (A{n)Z-Z)
(Ф -Z) <
І < N Ф і
2 2 ОО
84
А.С. Толстуха. Обтекание решеток произвольных лопастей ...
Или
N
Z < — Ф
оо 5
Оценим S по следующей схеме
(A(N)f()=J2^T ■ £t)VoIt
(18)
Е EE№'«s IтШг ) VolT = ^(РЇ'тМтС • MT£)VolT >
Т \іет j£T
> Y. Or ( - {Мт*.МтЄ)Мте I Voir
^ дт ((МтІ■ МтІ) - -(МтІ■ Мте)2 ) VolT
>
>
min {ртУоїт} у(Мт£ ■ Мт£) — -(Мт£- Мте)
Матрица М/■ = {m i ij )•ц {і 4} есть матрица перевода глобальной нумерации
узлов в локальную нумерацию в элементе Г. символически можно записать так: rriTij = 1, если loqr(j) = г, нуль в остальных случаях. 11 у - матрица размером 4 х 4 с элементами {(VI; у VI; у)}-( Г. 11 у есть вырожденная в силу тождества
ІЄТ
V Г VI, = 0 неотрицательно определённая матрица. //у - минимальное ненулевое собственное значение матрицы Wt- Так как (11 уД//ц • Д/уц ) = (Ст • Ст);
- - / 1
Ду = inf (Ст • ЦТ ) А ( ---------
(MT£-MT£)=l \°тах
(МтуМте)= О
(19)
Пусть /. - матрица 4 х 4 с элементами, равными единице. Заметим, что
1 /і ^ і -\/1
, (Д/,4, • Мте)2 = -LMT£ • -LMT£ = , /. Д/,4, • Д/-К
4 \4 4 J \4
Обозначим до = тт{дтТЫт} и продолжим описание схемы оценки 5:
(И^е-С) > До £ (МтС • Л/7 <Г) - ( jL МГ(■ Мт(
1
До У' ( ( Е - -L ) Мті • МтС
у \\ 4
Дп
ЕЕ^щЕф
Т ЩТ \ j<=T /
Дп
ЕЕ £ («■ -
Т ІЄТ j£T
Математические структуры и моделирование. 2003. Вып. 11.
85
Далее пусть i\...ц - путь, соединяющий узлы е номерами /(t. О,- при дви-
жении по рёбрам тетраэдров разбиения области, при этом
Сго = С,
. ™ах = Cmin = . Д1*11 , Сі-
N
Заметим, что £max и £т;п в силу условия Сг = 0 имеют разные знаки. Тогда
г=1
Е Е Е й - о2 * > 2 Е й, - щ.)2 г t Кт“~ Єті,,)2 >
N
(20)
Итак, следуя неравенствам (19), (20), в (18) можно принять
~ тіиІоД.
4iV • u’i. т
(21)
Оценим снизу объём тетраэдра, \оіт
1 hjT 1 / hjT
2^ ;>-V , = :С-іт
j£T іет
4
>
1 . / hjT
> - mm
_ З ієт V 4
ІЄТ
4тг . о
—— mm hjT.
43.3 ієТ Д
здесь LOjT - телесный угол, под которым грань тетраэдра, лежащая напротив узла с номером j, видна из средней точки ху = - х;-, hjT - высота тетраэдра,
ІЄТ
опущенная из вершины с номером j. В частности,
пип l ot >
(22)
Введём норму у • Цдг (17), используя (22), получим в этой норме
1
j-jjj j] • Ц^, Подстановкой в (18) оценок 6 (21) и
Ф
Z
< 27 • 3 • sup К D
« МЕГ! " AJ*
(23)
что означает сходимость с порядком 2 по шагу сетки в данной норме (в предположении ограниченности || D ||2),
Введённая норма является нормой в пространстве сеточных функций. Её невырожденность (в терминах последовательности сеточных функций это
ZN
означает, что условие lim
ІУ—кх>
= 0 влечёт lim Zдг = 0) обеспечивает
N N^oо
существование согласованной с ней нормы в пространстве решения исходной вариационной задачи.
86
А.С. Толстуха. Обтекание решеток произвольных лопастей ...
Ситуация станет более ясной, когда мы нормируем базисные функции
I;(х) = ,\\;(х). і Є 1....Y. Система нормированных базисных функций
оказывается сильно минимальной [15]. Приближённое решение ищем в виде ф(х) = £ .\'о(1 ((х). В этом случае точное значение щ рз N(f>i, и, как легко
г=1
убедиться, оценка для Ф совпадает с (17). Для ошибки определения коэффици-
7
ЄИТОВ Zs = фл--ПОЛУЧИМ
г Ґг N
N2A(n)Z = -Ф,
где А(Ю
- матрица, рассмотренная ранее. Проведя в точности те же рассу-
ждения, получим для ошибки определения нормированных коэффициентов Z оценку (23) уже в равномерной (|| • || ) норме.
Заключительное замечание относится к актуальному случаю неограниченного решения. Поведение решения и его производных в окрестности особых точек, таких как острая передняя кромка в случае решёток лопастей, известно.
Решение ведёт себя как ~ -, г - расстояние до особой точки, вторые производ-
ные -сетки
как ~
1
г
Принимая во внимание, что можно оценить количество узлов
Vok
< N < 26 • 3 •
л а„
Voln
л/lLv
где Voln - объём всей области, можно усилить нормировку исходных базисных функций, ДОМНОЖИВ ИХ на iV3/2, и получить сходимость коэффициентов приближённого решения в равномерной норме во всей области.
5. Заключение
Рассмотрена задача пространственного обтекания решётки стационарным потоком с учётом влияния вихревой пелены. Течение предполагалось потенциальным. Предложена постановка задачи, в которой вихревая пелена трактуется как поверхность контактного разрыва поля скоростей и реализована итерационная процедура определения распределения разрыва потенциала на задней кромке и положения вихревой пелены. Расчёты подтвердили, что распределение скоростей на поверхности лопастей и их аэродинамические характеристики в первую очередь определяется распределением циркуляции вдоль размаха лопасти, углом выхода потока - во вторую, и в меньшей степени дальнейшей эволюцией пелены.
Использован метод конечных элементов, базисом которых являются кусочно-линейные функции, определённые на тетраэдрах разбиения расчётной области.
Корректность постановки и программной реализации подтверждена сравнением с результатами, полученными ранее другими авторами.
Математические структуры и моделирование. 2003. Вып. 11.
87
Литература
1. Горелов Д.Н., Курзин В.Б., Сарен В.Э. Аэродинамика решёток в нестационарном потоке. - Новосибирск: Наука, 1971. 272 с.
2. Wo С.Н. A general theory of three-dimensional flow in a subsonic and supersonic turbomachines of axial, radial and mixed flow types // Trans. ASM I '.. 1952. V.74. P.1363-1380.
3. Сироткин А.Я., Степанов Г.Ю. Установившееся осесимметричное вихревое течение невязкой жидкости в многоступенчатых турбомашинах // МЖГ. 1981. N.6. С.13-15.
4. Заболотный Ф.Г. Расчёт, установившегося осесимметричного вихревого течения несжимаемой невязкой жидкости в радиально-осевой турбомашине // Изв. АН СССР. Энергетика и транспорт. 1979. N.3. С.147-155.
5. Ершова Л.Н. и др. К расчёту пространственного сжимаемого идеального газа в рабочих колесах турбомашин // Гидродинамика болвших скоростей. Труды 3 Всесоюзной школы-семинара, 8-14 июня 1987г. Красноярск. 1987. С.92-95.
6. Гречаниченко Ю.В., Нестеренко В.А. Вторичные течения в решётках турбомашини. - Харвков: Вища школа, 1983. 119 с.
7. Гнесин В.И., Соколовский Г.А. Расчёт, вторичных течений невязкого газа, в пространственных решётках турбомашин // Проблемы машиностроения (Киев). 1983. N.19. С 91-95.
8. Kawai Т., Adachi Т., Akashita К. Structure and Decay of Secondary Flow / Bull. JSME. 1985. V.28. N.242.
9. Гуревич, Шубин. О взаимодействии пограничных слоев в трёхмерных течениях // МЖГ. 1988. N.3.
10. Worster D.M. The calculation of fully three-dimensional flows in impellers using finite element method / Rep. from Dept, of Mech. Eng. at Heriot Watt Univ. U.K. 1973.
11. Бэтчелор Дж. Введение в динамику жидкости. М.: Мир, 1973. 758 с.
12. Рябченко В.П. Аэродинамические силы, действующие на лопасти, пространственной кольцевой решётки при нестационарном, обтекании // ПМТФ. 1979. N.4.
13. Михайлов В.П. Дифференциальные уравнения в частных производных. М.: Наука, 1983. 424 с.
14. Марчук Г.И. Методы вычислительной математики. М.: Наука, 1980. 534 с.
15. Марчук Г.И., Агошков В.И. Введение в проекционные методы. М.: Наука, 1981. 416 с. *
