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

Расчетная модель при численной оптимизации рабочих колес центробежных компрессоров Текст научной статьи по специальности «Механика и машиностроение»

CC BY
294
62
i Надоели баннеры? Вы всегда можете отключить рекламу.
Ключевые слова
РАБОЧЕЕ КОЛЕСО / ЦЕНТРОБЕЖНЫЙ КОМПРЕССОР / РАСЧЕТНАЯ ОБЛАСТЬ / СЕТОЧНАЯ МОДЕЛЬ / ОПТИМИЗАЦИЯ / NUMECA FINE/TURBO / IMPELLER / CENTRIFUGAL COMPRESSOR / COMPUTATIONAL DOMAIN / GRID MODEL / OPTIMIZATION

Аннотация научной статьи по механике и машиностроению, автор научной работы — Неверов Владимир Валерьевич, Кожухов Юрий Владимирович, Яблоков Алексей Михайлович, Лебедев Александр Анатольевич

В статье рассматриваются вопросы выбора параметров сеточной модели и расчетной области при решении задач оптимизации рабочих колес центробежного компрессора методами вычислительной газодинамики. При решении таких задач поиск и применение оптимальных параметров сеточной модели, расчетной области и настроек решателя позволяет обеспечить высокую точность моделирования при наиболее эффективном и производительном использовании машинного времени. Исследование проведено в комплексе Numeca Fine/Turbo для моделей турбулентности Spalart–Allmaras и Shear Stress Transport на примере двух рабочих колес: высоконапорного с ψт = 0,71, Ф = 0,064 и низконапорного ψт = 0,43, Ф = 0,06. Установлено, что выбор оптимальных параметров при постановке задачи значительно сокращает время получения сошедшегося решения и облегчают дальнейшее решение оптимизационных задач.

i Надоели баннеры? Вы всегда можете отключить рекламу.

Похожие темы научных работ по механике и машиностроению , автор научной работы — Неверов Владимир Валерьевич, Кожухов Юрий Владимирович, Яблоков Алексей Михайлович, Лебедев Александр Анатольевич

iНе можете найти то, что вам нужно? Попробуйте сервис подбора литературы.
i Надоели баннеры? Вы всегда можете отключить рекламу.

SIMULATION MODEL IN NUMERICAL OPTIMIZATION OF CENTRIFUGAL COMPRESSOR IMPELLERS

The article deals with choosing the parameters for the grid models and the computational domain while solving problems on optimizing the impellers of centrifugal compressors’s via computational fluid dynamics. In solving such problems, searching and applying the optimal parameters of the grid model, the computational domain, and the solver settings provides a high accuracy of the simulation with the most effective and productive use of machine time. The study was conducted in the Numeca Fine/Turbo program complex for the Spalart-Allmaras and Shear Stress Transport turbulence models on the example of two impellers: high-pressure with ѱт =0,71 and Ф = 0,064, and low-pressure with ѱт =0,43 and Ф = 0,06. It is established that the choice of the optimum parameters can significantly reduce the time required for obtaining a converging solutions and facilitates the further solution of optimization problems.

Текст научной работы на тему «Расчетная модель при численной оптимизации рабочих колес центробежных компрессоров»

DOI 10.5862/JEST.254.6 УДК 621.515

В.В. Неверов, Ю.В. Кожухов, А.М. Яблоков, А.А. Лебедев

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

V.V. Neverov, Yu.V. Kozhukhov, A.M. Yablokov, A.A. Lebedev

simulation model in numerical optimization of centrifugal compressor impellers

В статье рассматриваются вопросы выбора параметров сеточной модели и расчетной области при решении задач оптимизации рабочих колес центробежного компрессора методами вычислительной газодинамики. При решении таких задач поиск и применение оптимальных параметров сеточной модели, расчетной области и настроек решателя позволяет обеспечить высокую точность моделирования при наиболее эффективном и производительном использовании машинного времени. Исследование проведено в комплексе Numeca Fine/Turbo для моделей турбулентности Spalart—Allmaras и Shear Stress Transport на примере двух рабочих колес: высоконапорного с = 0,71, Ф = 0,064 и низконапорного = 0,43, Ф = 0,06. Установлено, что выбор оптимальных параметров при постановке задачи значительно сокращает время получения сошедшегося решения и облегчают дальнейшее решение оптимизационных задач.

РАБОЧЕЕ КОЛЕСО; ЦЕНТРОБЕЖНЫЙ КОМПРЕССОР; РАСЧЕТНАЯ ОБЛАСТЬ; СЕТОЧНАЯ МОДЕЛЬ; ОПТИМИЗАЦИЯ; NUMECA FINE/TURBO.

The article deals with choosing the parameters for the grid models and the computational domain while solving problems on optimizing the impellers of centrifugal compressors's via computational fluid dynamics. In solving such problems, searching and applying the optimal parameters of the grid model, the computational domain, and the solver settings provides a high accuracy of the simulation with the most effective and productive use of machine time. The study was conducted in the Numeca Fine/Turbo program complex for the Spalart-Allmaras and Shear Stress Transport turbulence models on the example of two impellers: high-pressure with =0,71 and Ф = 0,064, and low-pressure with =0,43 and Ф = 0,06. It is established that the choice of the optimum parameters can significantly reduce the time required for obtaining a converging solutions and facilitates the further solution of optimization problems.

IMPELLER; CENTRIFUGAL COMPRESSOR; COMPUTATIONAL DOMAIN; GRID MODEL; OPTIMIZATION; NUMECA FINE/TURBO.

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

жества вариативных расчетов исследуемой проточной части. В отличие от оптимизационных циклов, построенных на алгоритмах решения обратной задачи, математическое моделирование прямой задачи позволяет оценивать вклад всех геометрических параметров в характеристики проточных частей. Опыт кафедры КВиХТ в области оптимизации методами CFD и возможности данной процедуры отражены в работах [1—5].

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

Цель работы и объект исследования

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

Объекты исследования — два центробежных радиальных рабочих колеса с цилиндрическими лопатками:

РК1* — рабочее колесо со средней линией лопаток по дуге окружности. Условный расчетный коэффициент расхода Фр = 0,06; расчетный коэффициент теоретического напора ф^ = 0,43; диаметр рабочего колеса D2 = 0,82 м; относительный диаметр втулки = 0,23; угол выхода лопаток вл2 = 32°; количество лопаток z = 13;

РК2 — рабочее колесо с средней линией лопаток, образованной по двум сопряженным па-

* Список принятых обозначений и сокращений — см. Приложение в конце статьи

раболам. Условный расчётный коэффициент расхода Фр = 0,064; расчётный коэффициент теоретического напора ^ = 0,715; диаметр рабочего колеса D2 = 0,402 м; относительный диаметр втулки D = 0,3; угол выхода лопаток вл2 = = 75,7°; количество лопаток z = 24.

Построение сеточных моделей и проведение расчетов выполнены с помощью программного комплекса Numeca Fine/Turbo, преимущества которого описаны в работе [8], а основные аспекты работы приведены в [9]. Для всех расчетов задавались одинаковые условия: стационарная постановка; рабочее тело — идеальный воздух; граничное условие входа — полное давление и температура; граничное условие выхода — массовый расход. Моделирование произведено для двух моделей турбулентности — Spalart-Allmaras и Shear Stress Transport. В расчете не учитывались области уплотнений по покрывающему и основному дискам. Таким образом, приведенные значения политропного КПД — гидравлические. Все расчеты проведены для одного межлопаточного канала. Такой подход возможен, так как в других каналах течение периодически повторяется.

Результаты расчета отражены в значениях эффективности, которые рассчитаны по сечениям 0—0 (вход в РК) и 2'— 2' (выход из РК) согласно рис. 1; в количестве итераций, необходимых для получения условно сошедшегося решения (критерии — падение среднеквадратичных невязок до 10-4, выход на неизменное значение параметров эффективности и напора в процессе решения); в обезразмеренном времени, приходящемся на одну итерацию (значения отнесены к минимальному времени на одну итерацию в каждой серии расчетов); в обезразмеренном машинном времени, необходимом для получения условно сошедшегося решения (значения отнесены к минимальному машинному времени в каждой серии расчетов). Время, приходящееся на одну итерацию, в общем случае зависит от размерности сетки и, в случае применения, от разбития решения на процессорные ядра. Машинное время в совокупности зависит от количества итераций для получения решения и от времени на одну итерацию.

Исследование влияния положения выходной границы рабочего колеса

Отрывные течения и так называемый след после рабочего колеса могут оказать существенное влияние на результаты расчета, если выходная область находится близко к выходу из рабочего колеса. Исследование проводилось при фиксированной входной границе путем удлинения участка безлопаточного диффузора с Ь3 = Ь2 после рабочего колеса. Размер сеточных моделей увеличивался соответственно увеличению элементов на участке безлопаточного диффузора (БЛД) при его удлинении. Эскиз вариаций выходных границ расчетной области приведен на рис. 1. Съем параметров для определения эффективности произведен в сечениях 0—0 и 2'—2' (В2 = 1,05А2). Результаты расчета для обоих рабочих колес представлены на рис. 2 и 3.

Для низконапорного колеса разница полученных значений КПД при отдалении выходной границы расчетной области невелика и при значении Авых = 1,15В2 уже выходит на постоянную величину. Количество необходимых итераций для сведения решения имеет степенную зависимость и значительно снижается при отдалении выходной границы на 1,3А2 и далее. Однако при отдалении границы возрастает и общее количество элементов, что ведет к увеличению времени, приходящегося на одну итерацию. В совокуп-

ности наиболее оптимальным вариантом с точки зрения использования машинного времени оказался вариант с выходной границей, расположенной на 1,3А2. Дальнейшее увеличение не дает значительного уменьшения необходимых для сходимости итераций, но ведет к увеличению затрачиваемого времени на одну итерацию. В случае высоконапорного колеса положение выходной границы сильнее влияет на получаемую в результате решения эффективность рабочего колеса из-за большей неоднородности потока на выходе с лопаток. Даже по этой причине в высоконапорных рабочих колесах следует отодвигать выходную границу на 1,25А2 и более. С другой стороны, по результатам расчета сходимость решения резко ухудшается после 1,2А2. Количество необходимых итераций возрастает в 2,3 раза на модели турбулентности SA при удалении границы с 1,2А2 до 1,4А2. Это связано с развитыми отрывными течениями в диффузоре на удлиненных участках и, следовательно, с пульсациями и неоднородностями массового расхода на выходной границе, что замедляет схождение решения. На приведенных графиках ширина диффузора равна ширине канала на выходе из рабочего колеса, т.е. Ь/Ь2 = 1. Большая ширина диффузора способствует возникновению в нем отрывного течения, поэтому, кроме исходного варианта, был рассмотрен вариант

1,4 А

Вариации выходной границы

1,35А2 1,25А2 1,15 А2 1,05А.

2-2

Вариации входной границы 0-0

Рис. 1. Схема контрольных сечений и вариации расчетной области

а)

Пп 0,948

0.946

0,944

0,942

0,940

0,938

0,936

0,934

0,932

0,930

Относительная итерационная сходимость

б) Относительное время на Относительное машинное 1 итерацию время сходимости

1 1

-■»■■PK1-SA -■А- PTC1-SST -■-PK1-SA -H-PK1-SST ■юлитропный КПД политропный КПД

dTepat^ итерац онная сходимость ионная сходимость

\ ... ....... ■ - - ■

\

ж'' ь...... tiVj к1"""* ........ *.......

2,2 2 1,8 1,6 1,4 1,2 1

1,05 1,10 1,15 1,20 1,25 1,30 1,35 1,40

1,24 1,21 1,18 1,15 1,12 1,09 1,06 1,03 1

А- ■РК1-5А время на 1 итерацию - -А- 'РКХ-ББТ время на 1 итерацию —РКХ-БА машинное время сходимости и РК1-55Т машинное время сходимости

/ с'

/

> 1

» к

** \ с'

Л* . J ...... X

л-—

1,8 1,7 1,6 1,5 1,4 1,3 1,2 1,1 1

Относительный диаметр границы выхода

1,05 1,10 1,15 1,20 1,25 1,30 1,35 1,40

Относительный диаметр границы выхода

Рис. 2. Результаты оценки влияния выходной границы для низконапорного колеса РК1 по сечениям

0-0 (а) и 2'-2' (б)

а)

Относительная итерационная сходимость

б) Относительное время на 1 итерацию

Относительное машинное время сходимости

0,954 0,952 0,950 0,948 0,946 0,944 0,942 0,940 0,93В 0,936 0,934 0,932 0,930 0,928 0,926

-А- •■а.» PK2-SA политропный КПД PK2-SST политропный КПД FK2-SA итерационная сходим PK2-SST итерационная сходи

А. ость

**1 чисть

4 /

*

; ;

к \

; .......1 1г......А- ■■■/ч

; :

; Г : у

L. д

1,05 1Д0 1Д5 1,20 1,25

3,8 3,6 3,4

3.2 3

2.3 2,6

2.4 2,2 2

1.3 1,6

1.4 1,2 1

1,30 1,35 1,40 Относительный диаметр границы выхода

1,28 1,26 1,24 1,22 1,2 1,18 1,16 1,14 1,12 1.1 1,08 1,06 1,04 1,02 1

•■А --A--F —я— JK2-SA время на 1 итерацию >K2-SST время на 1 итерацию JK2-SA машинное время сходимости 'K2-SST машинное время сходимости

L

M yi.-'i

* У

.J V

....."

я'' / г fl j я

/

й

iНе можете найти то, что вам нужно? Попробуйте сервис подбора литературы.

1,05 1,10 1,15 1,20

3,8 3,6 3,4 3,2 3 2,8 2,6 2,4 2,2 2 1,8 1,6 1,4 1.2 1

1,25 1,30 1,35 1,40

Относительный диаметр границы выхода

Рис. 3. Результаты оценки влияния выходной границы для высоконапорного колеса РК2 по сечениям

0-0 (а) и 2—2' (б)

с Ь3/Ь2 = 0,8. Сужение диффузора по результатам расчета не дало эффекта — после удаления выходной границы от 1,2А2 не только не улучшилась сходимость, но и общее время сходимости увеличилось по сравнению с исходной постановкой задачи. В данном случае наилучший вариант — это выбор выходной границы 1,25А2, где в результате решения значение эффективности уже можно считать приемлемым, но и скорость сходимости еще не критично возросла.

Суммируя данные по влиянию выходной границы, можно утверждать: для высоконапорных колес положение выходной границы более значительно влияет на значение КПД численного эксперимента; из-за большей неоднородности потока по сравнению с низконапорными колесами на сведение задачи требуется большее количество итераций и, следовательно, машинного времени. На участке 1,05—1,2 А2 количество необходимых итераций уменьшилось в 1,4 раза против 2,1 у низконапорного варианта; при увеличении Авых свыше 1,25А2 сходимость задачи резко ухудшается из-за отрывных течений и нестационарных явлений на участке БЛД. Для низконапорных колес влияние положения выходной границы на эффективность сравнительно мало; при удлинении безлопаточного участка сходимость значительно улучшается, что позволило добиться в рассматриваемом случае снижения машинного времени расчета более чем в 1,7 раза. Логически удлинение участка оказывается ограничено оптимальным соотношением числа итераций, необходимого для сведения задачи, и времени, приходящегося на одну итерацию.

Во всех сериях расчетов отмечено следующее: модель турбулентности SA, будучи более простой, чем SST, дает на 10—20 % более быструю итерационную сходимость решений и на 5—15 % меньшее время, приходящееся на итерацию. А модель SST, в свою очередь, завышает КПД в абсолютных значениях на 0,5—1 % по сравнению с моделью SA, что было установлено не только в данной работе, но и в других исследованиях, проводимых на кафедре КВиХТ СПбПУ в пакете Numeca. В целом кривые для обеих моделей турбулентности эквидистантны.

Исследование влияния положения входной границы рабочего колеса

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

В контексте оценки оптимальной длины входного участка перед рабочим колесом рассматривались значения длин 0—0,75 А2. Результаты расчета приведены на рис. 4 и 5.

Для низконапорного колеса уменьшение длины входного участка привело к линейному увеличению эффективности на модели турбулентности SA (для рассмотренной задачи — 0,3 %) и к постоянному значению КПД для модели SST. При нулевой длине произошло резкое падение эффективности на 0,2 % для обеих используемых моделей турбулентности; оно вызвано тем, что профиль скорости на данной координате уже перестраивается для входа на лопатки и отличается от задаваемого однородным полем в случае расположения входной границы на этой координате, что ведет к изменению обтекания лопаток и разнице в результатах. Поэтому, несмотря на продолжающееся ускорение сходимости, не следует выбирать слишком короткий входной участок.

В общем, положение координаты, где профиль скорости значительно перестроится, зависит от комплекса факторов, например от ширины

а)

Пп

0,946

Относительная итерационная сходимость

б) Относительное время на 1 итерацию

Относительное машинное время сходимости

д «РК1-5А политропный КПД ■А«-РКг^Т политропный КПД —и^РК1-5А итерационная сходимость ш РК1-55Т итерационная сходимость

к.......ч

.......*....... *....... ч

**

..Л Г"

• -1

0,75

0,625

0,5

0,375 0,25 0,125

2,6 2,4 2,2 2 1,8 1,6 1,4 1,2

Длина входного участка, отнесенная к диаметру А

2,6

2.5

2.4

2.3 2,2 2,1

2 1,9 1,8 1,7

1.6

1.5

1.4 1,3 1,2 1,1

1

-■А- РК1-5А время на 1 итерацию - -А- -РК1-55Т время на 1 итерацию —а— РК1-5БТ машинное время сходимости ■ РК1-5А машинное время сходимости

д

\ ь.. ч

'-1 \

5 \

Хг

.....

/ *

0,75 0,625 0,5 0,375 0,25 0,125 0

Длина входного участка, отнесенная к диаметру А2

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

0-0 (а) и 2—2' (б)

а)

Пп

0,948

0,946 0,944 0,942 0,940 0,938 0,936 0,934 0,932

Относительная итерационная сходимость

б) Относительное время на Относительное машинное 1 итерацию время сходимости

Д РК2-5А политропный КПД -чА"РК2-Я5Т политропный КПД п РК2-5А итерационная сходимость РК2-55Т итерационная сходимость

0,75 0,625

0,5

1,6 1,4 1,2 1

0,375 0,25 0,125 О

Длина входного участка, отнесенная к диаметру А2

I*

**

г*

■-, ^ д 1* ^^

к***

2,6 2,4

1,8

1,65 1,6 1,55 1,5 1,45 1,4 1,35 1,3 1,25 1,2 1,15 1,1 1,05 1

л \ - Д РК2-5А время на 1 итерацию - -а - -РК2-55Т время на 1 итерацию —В— РК2-5А машинное время сходимости ■нн-РК2-55Т машинное время сходимости

\

1\ ■ ( у

>2

0,75 0,625 0,5 0,375

3,6 3,4 3,2 3 2,8 2,6 2,4 2,2 2 1,8 1,6 1,4 1,2 1

0,25 0,125 О Длина входного участка, отнесенная к диаметру А2

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

0-0 (а) и 2—2' (б)

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

Исследование влияния выбора сеточных параметров

Вопрос выбора размерности сетки при верификации численных экспериментов решается путем стандартного исследования на сеточную независимость. Такое исследование позволяет определить, начиная с какого количества узлов получаемое решение практически перестает зависеть от густоты сеточной модели. В оптимизационных задачах необходим более точный выбор сеточных параметров, обеспечивающих приемлемые результаты при минимальном затрачиваемом времени на решение. Сеточный генератор AutoGrid 5 программного комплекса №теса позволяет автоматически генерировать блочно-структурированные параметризированные сетки. В автоматическом режиме для радиальных колес обычно можно использовать два вида топологии: НОН и Н&1. Принципиальная разница между сетками с разными топологиями показана на рис. 6 и заключается в том, что для сетки типа Н&1, в отличие от типа НОН, характерно межсекторное разделение вблизи поверхности ло-

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

В общем случае выбор варианта топологии обосновывается качеством получаемой сеточной модели: минимальной скошенностью элементов, максимальным коэффициентом расширения соседних элементов и т.д. Поэтому не для всех конфигураций рабочих колес можно применить оба варианта топологии. В случае, когда для рабочего колеса можно построить равноценные по качеству сетки для о беих топологий, следует про -изводить выбор топологии с точки зрения скорости получаемого решения. Рассматривались сетки, построенные с помощью двух топологий для рабочих колес. Количество сеточных элементов сеток РК1 и РК2 — 1,15 млн, минимальный угол скошенности для РК1 — 45°, для РК2 — 36°. Результаты расчета для режима расчетного расхода представлены в табл. 1.

Таким образом, выбор оптимальной сеточной топологии позволяет по результатам исследования экономить до 20 % машинного времени, что ощутимо при проведении большого количества расчетов.

Стандартная процедура при проведении разнообразных численных экспериментов — исследование на сеточную независимость. Выводы по таким исследованиям приведены во многих работах [10—12].

а) б)

Рис. 6. Сеточная топология НОН (а) и Н&1 (б)

Таблица 1

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

РК1 РК2

Модель турбулентности Сеточная топология Машинное время, сек. Число необходимых итераций Машинное время, сек. Число необходимых итераций

SA НОН 1538 696 1611 1070

Н&1 1581 835 1633 1305

SST НОН 1630 830 1890 1390

Н&1 1652 1022 1900 1606

С уменьшением размеров элементов сетки происходит увеличение их количества, в связи с чем при расчетах изменение диссипации энергии происходит более плавно, а расчетные потери КПД в элементах проточной части уменьшаются. Для рассматриваемых рабочих колес произведена оценка и выбор оптимального количества элементов для дальнейших оптимизационных расчетов. Результаты для низконапорного РК1 представлены на рис. 7. Результаты для РК2 качественно повторяют результаты для РК1.

iНе можете найти то, что вам нужно? Попробуйте сервис подбора литературы.

С ростом числа элементов сетки КПД низконапорного РК1 практически перестает меняться,

начиная с 0,8 млн элементов. Время, приходящееся на одну итерацию, растет линейно. Количество итераций для сведения решения имеет зависимость от числа элементов, близкую к линейной. В итоге машинное время получения сошедшегося решения нелинейно зависит от размеров сеточных моделей и значительно возрастает при размерах, превышающих 1 млн элементов. Поэтому для данного рабочего колеса рационален выбор сеточной модели размером 0,8—1 млн элементов для одного лопаточного сектора — это практически не оказывает влияния на получаемые результаты при экономии вычислительных и временных ресурсов.

а) п

•и

0,945

Относительное время итерационной сходимости

б) Относительное время Относительное машинное на 1 итерацию время сходимости

д-р ■■А--Р —НН*Р 'К1-БА политропный КПД 'К1-55Т политропный КПД 'К1-БА итерационная сходимос 'К1-55Т итерационная сходимо

сть ¿А

..... ....А .....

д д

4

»* д

д

2,6

2,4

2,2

1,6 1,4 1,2

О 0,2 0,4 0,6 0,8 1 1,2 1,4 1,6 1,8 2 2,2

Количество сеточных элементов, млн

А

Д РК1-БА время на 1 итерацию ■ -а - ■РК1-55Т время на 1 итерацию =И1— РК1-БА машинное время сходимос ■ РК1-55Т машинное время сходимо

сти / А

/

4 /у

у д

/ /

£ д

А д >

Л Ал

25

22

19

16

13

10

0 0,2 0,4 0,6 0,8 1 1,2 1,4 1,6 1,8 2 2,2

Количество сеточных элементов, млн

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

по сечениям 0—0 (а) и 2—2' (б)

Сеточный параметр у+ — безразмерное расстояние от первого узла сетки до стенки. Данный критерий важен при выборе модели турбулентности (низко- или высокорейнольдсовая), которая определяет требование к сетке. Например, при использовании стандартных пристенных функций (высокоренольдсовые модели) необходимо избегать значений у+ меньших, чем приблизительно 30, и больших, чем 300, для задач расчёта потока газа. С масштабируемыми пристенными функциями (низкорейнольдосвые модели) значения этого параметра дают информацию о разрешающей способности сетки вблизи стенки. Использование масштабируемой пристенной функции обязывает иметь достаточно тонкий слой элементов, чтобы разложить пограничный слой. Рекомендуемый в литературе диапазон у+ для таких моделей находится в интервале от 0,001 до 1. Для определения данного параметра и его дальнейшей корректировки обычно проводится пробный расчет. Определяется значение критерия у+ во всей расчетной области, после чего корректируется глобально или локально размер первой пристенной ячейки для удовлетворения критерия. Влияние параметра рассмотрено в работе [13].

Считается, что низкорейнольдсовые модели точнее в предсказании точек отрывов, но соот-

ветственно требовательнее к размерности сетки. Рассматриваемые в работе модели — низкоре-нойльдсовые. Рассмотрено влияние параметра у+ в рекомендуемом диапазоне от 0,001 до 1 и при выходе за него (от 1 до 48) на решение и сходимость задачи.

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

Как видно из графиков, попадание у+ в рекомендуемый диапазон от 0,001 до 1 оказывает минимальное воздействие на сходимость решения. То есть скорость сходимости в данном диапазоне остается постоянной. Однако получаемые результаты эффективности достаточно сильно различаются: при увеличении среднего значения параметра у+ с 0,17 до 0,82 КПД низконапорного рабочего колеса увеличивается на 0,5 %, а затем при у+ = 1,75 падает на 0,7 %. Для высоконапорного колеса результаты качественно аналогичны, за исключением некоторого ускорения сходимости для модели SST при стремлении у+ к 1,17.

Таблица 2

Вариации размера ячеек и значения параметра у+

Рабочее колесо и модель Номер Размер ячейки, Средний Максимальный

турбулентности варианта мм у+ у+

1 0,001 0,17 0,59

2 0,0015 0,25 0,86

3 0,005 0,82 3,03

РК1 8А&88Т 4 5 0,01 0,03 1,75 5,5 6,23 14,6

6 0,08 12,6 25,2

7 0,15 19,7 34,5

8 0,30 30,4 47,7

1 0,0005 0,12 0,165 0,276 1,17 2,57,3 16,0 24,4 0,38

2 0,0007 0,53

3 0,0012 0,92

РК2 8А&88Т 4 5 0,005 0,01 4,08 7,8

6 0,03 16,1

7 0,08 28,1

8 0,15 39,9

а)

Пп 0,99

0,98

0,97

0,96

0,95

0,94

0,93

0,92

0,91

0,90

0,89

0,88

Относительное время итерационной сходимости

д РК1-БА политропный КПД ■ -а - -РК^БТ политропный КПД РК1-5А итерационная сходимость и РК1-58Т итеоапионная сходимость

V

1*7

• ■ ч /

3,2 3 2,8 2,6 2,4 2,2 2 1,8 1,6 1,4 1,2 1

5 6 7 8

Варианты высот первой пристенной ячейки

б) Относительное время на 1 итерацию

1,24 1,22 1,2 1,18 1,16 1,14 1,12 1.1

1,06 1,04 1,02 1

А РК1-5А время на 1 итерацию ■ ■А,,РК1-55Т время на 1 итерацию —РК1-5А машинное время сходимости и РК1-55Т машинное время сходимости

........ . ■. • • ^ г*" к

...... г* т

г*"*"

Относительное машинное время сходимости

3,4 3,2 3 2,8 2,6 2,4 2,2 2 1,8 1,6 1,4 1.2

5 6 7 8

Варианты высот первой пристенной ячейки

Рис. 8. Оценки влияния параметра у+ сеточной модели для низконапорного колеса

Увеличение у+ за рекомендуемый интервал ведет к некорректному росту эффективности колеса вплоть до 99 % КПД абс. Это связано со значительным уменьшением потерь, в том числе и в результате нивелирования отрывных течений из-за недостаточной размерности сетки для разрешения пограничного слоя. Вместе с некорректными результатами значительно возрастает и количество итераций, необходимых для получения сошедшегося решения. При примерно одинаковом времени, затрачиваемом на одну итерацию, машинное время получения сведенного решения возрастает только пропорционально необходимому количеству итераций. В случае высоконапорного колеса РК2 результаты и время сходимости ведут себя аналогично варианту низконапорного РК1. Однако решение для модели ББТ показало более жесткие требования к параметру у+: при у+>2,5 решение прерывается. Таким образом, для модели ББТ при расчетах высоконапорных рабочих колес прерывание процесса решения может быть связано не только с общим качеством сетки, но и с неприемлемо большим размером первого пристенного элемента. Модель 8Л в этом случае пока-

зала возможность расчета даже при значительном отклонении параметра у+ от рекомендуемого.

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

iНе можете найти то, что вам нужно? Попробуйте сервис подбора литературы.

при выборе случайных параметров входная граница определена на удалении 0,5D2, выходная граница — 1,1-^; выбрана сеточная топология Н&1; высота пристенной ячейки 0,01мм (соответствующее среднее значение у+ = 1,75); общее количество элементов 1,15 млн;

при выборе оптимальных параметров входная граница определена на удалении 0,125—2, выходная граница — 1,3—2; выбрана сеточная топология НОН; высота пристенной ячейки 0,001 мм (соответствующее среднее значение у+ = 0,17); общее количество элементов 0,92 млн.

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

данные цифры соответственно равны 2,4 секунды на итерацию, 227 итераций для сошедшегося решения и общее время расчета — 535 секунд, что в 4,25 раза быстрее, чем при неоптимальных параметрах расчетной модели. Таким образом, 100 оптимизационных шагов заняли бы 64 часа для случайных и 15 часов для оптимальных параметров. Это существенно, учитывая тот факт, что визуально данные расчетные области и сеточные модели схожи. Также в комплексе Nu-meca реализованы специальные экспертные настройки CPU Booster, которые позволяют производить моделирование с числом Куранта, равным 1000. Число Куранта (называемое также число Куранта—Фридрихса—Леви) характеризует физический масштаб времени расчёта и определяется по формуле

uAt

CFL =-,

Ax

где u — скорость переноса, At — временной шаг, Ах — пространственный шаг.

Увеличение числа Куранта вручную позволяет ускорить процесс расчета задачи при большой серии вариантных оптимизационных расчётов. Данная опция показывает свою работоспособность не для всех видов задач. Однако, если удается произвести запуск с такими настройками, то получение сошедшегося решение значительно ускоряется. Для оптимальных модельных параметров с функцией решателя CPU Booster 100 оптимизационных шагов по результатам расчета заняли бы 11,5 часов.

Корректно подходя к задаче оптимизации проточной части, будь то рабочее колесо центробежного компрессора или ступень турбины, можно экономить вычислительные ресурсы и трудозатраты расчетчика. Это особенно актуально в настоящее время, так как перечень решаемых задач постоянно расширяется, а объемы вычислительных процессов неуклонно растут. Проведение не усложненных предварительных исследований позволяет выполнить в разы больше оптимизационных шагов или даже задач при тех же затрачиваемых вычислительных машинных ресурсах. А наработки по валидации расчетных моделей [14—16] позволяют более точно оценивать итоговые характеристики оптимизированной проточной части.

Заключение

В результате исследования выявлены наиболее значимые при проведении больших серий расчетов параметры расчетной модели рабочих колес центробежного компрессора — положение выходной границы и общая размерность сетки. Эти параметры вне зависимости от типа рабочего колеса при неоптимальном выборе увеличивают время получения сошедшегося решения в 2 раза и более, а также оказывают существенное влияние (до 1,7 % КПД для выходной границы и 1,5 % КПД для сеточной модели) на результаты решения. Это неприемлемо при оптимизации геометрической формы рабочих колес как для получения аргументированного результата, так и с точки зрения рационального использования вычислительных ресурсов.

Кроме того, для адекватности получаемых результатов следует соблюдать сеточный параметр у+ для выбранной модели турбулентности и не допускать влияния на общие результаты расчётов и картину обтекания лопаток слишком близко расположенной к рабочему колесу входной границы. Слишком большие значения критерия расчетной сетки у+ оказывают наибольшее влияние на результаты решения. Из-за низкой разрешающей способности сетки в пограничном слое нивелируются отрывные течения и уменьшаются потери трения. Поэтому при неправильном выборе размера пристенной ячейки политропный КПД рабочего колеса может увеличиваться до 99 % и более.

По результатам предварительной оценки определены оптимальные параметры расчетной модели при постановке задачи оптимизации. Для низконапорного колеса РК1 выходная граница определена на расстоянии 1,3А2; входная граница отодвинута от входа в рабочее колеса на 0,125А2; сеточная топология — НОН; размер первой пристенной ячейки — 0,002 мм; общий размер сеточной модели лопаточного сектора — 0,9 млн элементов. Для высоконапорного рабочего колеса РК2 выходная граница определена на 1,25А2; входная граница отодвинута от входа в рабочее колеса на расстояние 0,125А2; сеточная топология—НОН; размер первой пристенной ячейки — 0,001мм; общий размер сеточной модели лопаточного сектора 0,7 млн элементов. Данные наборы параметров модели обеспечивают наименьшее машинное время сходимости одного оптимизационного приближения.

Приложение

Принятые обозначения: "фт — коэффициент теоретического напора; Ф — условный коэффициент расхода;

D2 — наружный диаметр рабочего колеса центробежного компрессора;

b2 — ширина канала на выходе из рабочего колеса;

b3 — ширина канала на входе в диффузор;

y+ — сеточный параметр;

РК — рабочее колесо;

БЛД — безлопаточный диффузор;

КПД — коэффициент полезного действия, эффективность;

CFD — Computational Fluid Dynamics (вычислительная гидрогазодинамика);

SA — модель турбулентности Spalart—Allmaras;

SST — модель турбулентности Shear Stress Transport.

список

1. Kozhukhov y.v., Danilishm A.M., Yun v.K. Multi-objective optimization for impeller shroud contour, width of vane diffuser and number of blades of the centrifugal compressor stage based on the CFD calculation // IOP Conference Series Materials Science and Engineering 08/2015. Vol. 90(1):012046. DOI:10.1088/1757-899X/90/1/012046

2. Frese F., Einzinger J., Will J. Design optimization of an impeller with CFD and Meta-Model of optimal Prognosis (MoP) //10th International conference of tur-bocharges and turbocharging. 2012, London. P. 121—135.

3. Bonaiuti D., Zangeneh M. On the Coupling of Inverse Design and Optimization Technigues for the Multi-objective, Multipoint Desing of Turbomachinery Blades // Journal of turbomachinery-transaction of the ASME. Vol. 131(2). Article 021014. DOI: 10.1115/1.2950065.

4. demeulenaere a., bonaccorsi J.-c., Gutzwiller D., Hu L., sun H. Multi-Disciplinary Multi-Point Optimization of a Turbocharger Compressor Wheel. // ASME Turbo Expo 2015: Turbine Technical Conference and Exposition (Montreal, Canada). Paper № GT2015-43631. P. V02CT45A020. DOI: 10.1115/GT2015-43631.

5. dickmann H.P. Shroud contour optimization for a turbocharger centrifugal compressor trim family // 10th European Turbomachinery Conference. 04/2013. Lap-peenranta, Finland. Paper № ETC2013-015. 12 p.

6. boldyrev у., Rubtsov A., Kozhukhov у., lebedev A., cheglakov I., danilishin A. Simulation of unsteady processes in turbomachines based on nonlinear harmonic NLH-method with the use of supercomputers // CEUR Workshop Proceedings. 2015. Vol. 1482. P. 273-279.

7. Блинов В.Л., Бродов Ю.М., Седунин В.А., Комаров О.В. Выбор параметров расчетной модели при решении задач многокритериальной оптимизации плоских компрессорных решеток // Компрессорная техника и пневматика. 2015. № 1. С. 32-36.

8. Чеглаков И.В. Разработка методики проведения численного эксперимента и оптимизация не-

литературы

подвижных элементов малорасходных ступеней в среде моделирования Fine/Turbo на примере ступени СВД-22: Дис. ... магистра / СПбПУ, Санкт-Петербург. 2015.

9. Батурин О.В. Методика цифрового моделирования осевых многоступенчатых турбин низкого давления с учётом неравномерного поля параметров на входе в турбину, трёхмерной структуры потока в лопаточных венцах и утечек через радиальные зазоры лабиринтных уплотнений. [Электронный ресурс]: Электрон. учебное пособие / Батурин О.В., Колма-кова Д.А., Матвеев В.Н., Попов Г.М., Шаблий Л.С. Самарский гос. аэрокосм. университет им. С.П. Королева. Самара, 2012. 122 с. (1 CD-ROM).

10. Свобода Д.А., Жарковский А.А. Экспериментальные и расчетные исследования осевого насоса с быстроходностью ns = 560 // Известия Самарского научного центра Российской академии наук. 2013. Т. 15 № 4(2). С. 573-578.

11. Le sausse P., Fabrie P., Arnou D., clunet F. CFD comparison with Centrifugal Compressor Measurements on a wide Operating Range // Johnson Controls Industries, In EPJ Web of Conferences. 2013. Vol. 45. P. 01059.

12. v. v. N. K. satish Koyyalamudi, Quamber H. Nagpurwala. Stall Margin Improvement in a Centrifugal Compressor through Inducer Casing Treatment // International Journal of Rotating Machinery. 2016. Article ID 2371524, 19. D0I:10.1155/2016/2371524

13. Aghaeitog R., Tousi a. M., Tourani a. Comparison of turbulence methods in CFD analysis of compressible flows in radial turbomachines. // Aircraft Engineering and Aerospace Technology. 2008. Vol. 80, Iss. 6. P. 657-665. ISSN 1748-8842.

14. Яблоков А.М., Кожухов Ю.В., Лебедев А.А. Исследование течения в малорасходной ступени центробежного компрессора методами вычислительной газодинамики // Научно-технические ведомости

enemy. 2015. № 4 (231). C. 59-69. D01:10.5862/ JEST.231.7. (rus.)

15. Kozhukhov Y.v., Yun v.K., Reshetnikova L.v., Prokopovich M.v. Numerical Investigation of Different Radial Inlet Forms for Centrifugal Compressor and Influence of the Deflectors Number by Means of Computational Fluid Dynamics Methods with Computational

Model Validation // IOP Conference Series Materials Science and Engineering. August 2015. Vol. 90(1). Iss. 012047. DOI: 10.1088/1757-899X/90/1/012047

16. Kang shun, Lm Qiang, Ming-Xu Qi.CFD validation of a high speed centrifugal compressor impeller // Journal of Engineering Thermophysics. 2005. Vol. 26(3). P. 400-404.

references

1. Kozhukhov Y.v., Danilishin A.M., yun v.K. Multi-objective optimization for impeller shroud contour, width of vane diffuser and number of blades of the centrifugal compressor stage based on the CFD calculation. IOP Conference Series Materials Science and Engineering. 08/2015. Vol. 90(1):012046. DOI:10.1088/1757-899X/90/1/012046

2. Frese F., Einzinger J., Will J. Design optimization of an impeller with CFD and Meta-Model of optimal Prognosis (MoP). 10th International conference of turbocharges and turbocharging. 2012, London. P. 121-135.

3. Bonaiuti D., Zangeneh M. On the Coupling of Inverse Design and Optimization Technigues for the Multi-objective, Multipoint Desing of Turbomachinery Blades // Journal of turbomachinery-transaction of the ASME. Vol. 131(2). Article 021014. DOI: 10.1115/1.2950065.

4. demeulenaere a., bonaccorsi J.-c., Gutzwiller D., Hu L., sun H. Multi-Disciplinary Multi-Point Optimization of a Turbocharger Compressor Wheel. ASME Turbo Expo, 2015. Turbine Technical Conference and Exposition. Montreal, Canada. Paper № GT2015-43631. P. V02CT45A020. DOI: 10.1115/GT2015-43631.

5. Dickmann H. P. Shroud contour optimization for a turbocharger centrifugal compressor trim family. 10th European Turbomachinery Conference. Lappeenranta, Finland. 04/2013. Paper № ETC2013-015. 12 p.

6. boldyrev Y., RuMsov A., Kozhukhov Y., lebedev A., cheglakov i., danilishin A. Simulation of unsteady processes in turbomachines based on nonlinear harmonic NLH-method with the use of supercomputers. CEUR Workshop Proceedings. 2015. Vol. 1482. P. 273-279.

7. blinov v.L., brodov Yu.M., sedunin v.A., Koma-rov o.v. Vybor parametrov raschetnoy modeli pri resh-enii zadach mnogokriterialnoy optimizatsii ploskikh kompressornykh reshetok. Kompressornaya tekhnika i pnevmatika. 2015. № 1. S. 32-36 (rus.)

8. cheglakov i.v. Razrabotka metodiki provedeniya chislennogo eksperimenta i optimizatsiya nepodvizhnykh elementov maloraskhodnykh stupeney v srede mode-lirovaniya Fine/Turbo na primere stupeni SVD-22: Dis. ... magistra / SPbPU. Sankt-Peterburg, 2015. (rus.)

9. baturin o.v. Metodika tsifrovogo modelirovaniya osevykh mnogostupenchatykh turbin nizkogo davleniya s uchetom neravnomernogo polya parametrov na vkhode v

turbinu, trekhmernoy struktury potoka v lopatochnykh ventsakh i utechek cherez radialnyye zazory labirintnykh uplotneniy. [Elektronnyy resurs]. : Elektron ytshebnoe posobie / Baturin O.V., Kolmakova D.A., Matveev V.N., Popov G.M., Shabliyy L.S. Samarskiy gos. aerokosm. universitet im. S. P. Koroleva. Samara, 2012. 122 s. (rus.)

10. svoboda D.A., Zharkovskiy A.A. Eksperimental-nyye i raschetnyye issledovaniya osevogo nasosa s bys-trokhodnostyu ns = 560. Izvestiya Samarskogo nauchnogo tsentra Rossiyskoy akademii nauk. 2013. T. 15, № 4(2). S. 573-578. (rus.)

11. Le sausse P., Fabrie P., Arnou d., clunet F. CFD comparison with Centrifugal Compressor Measurements on a wide Operating Range, Johnson Controls Industries. EPJ Web of Conferences. 2013. Vol. 45. P. 01059.

12. v. v. N. K. satish Koyyalamudi, Quamber H. Nagpurwala. Stall Margin Improvement in a Centrifugal Compressor through Inducer Casing Treatment, International Journal of Rotating Machinery. 2016. Vol. 2016, Article ID 2371524, 19 p. D0I:10.1155/2016/2371524

13. Aghaeitog R., Tousi a. M., Tourani a. Comparison of turbulence methods in CFD analysis of compressible flows in radial turbomachines. In: Aircraft Engineering and Aerospace Technology. 2008. Vol. 80, Iss. 6. P. 657665. ISSN 1748-8842.

14. Yablokov A.M., Kozhukhov Yu.v., lebedev A.A. Issledovaniye techeniya v maloraskhodnoy stupeni tsen-trobezhnogo kompressora metodami vychislitelnoy gazo-dinamiki. Nauchno-tekhnicheskiye vedomosti SPbGPU. 2015. № 4(231). S. 59-69. D0I:10.5862/JEST.231.7 (rus.)

15. Kozhukhov Y.v., yun v.K., reshetnikova L.v., Prokopovich M.v. Numerical Investigation of Different Radial Inlet Forms for Centrifugal Compressor and Influence of the Deflectors Number by Means of Computational Fluid Dynamics Methods with Computational Model Validation. IOP Conference Series Materials Science and Engineering. August 2015. Vol. 90(1):012047. DOI: 10.1088/1757-899X/90/1/012047

16. Kang shun, Lm Qiang, Ming-Xu QI. CFD validation of a high speed centrifugal compressor impeller. Journal of Engineering Thermophysics. 2005. Vol. 26(3). P. 400-404.

сведения об abtopax/authors

НЕВЕРОВ Владимир Валерьевич — епециалист Санкт-Петербургского политехнического университета Петра Великого.

195251, Россия, г. Санкт-Петербург, Политехническая ул., 29. E-mail: [email protected]

NEVEROV Vladimir V. — Peter the Great St. Petersburg Polytechnic University. 29 Politechnicheskaya St., St. Petersburg, 195251, Russia. E-mail: [email protected]

КОЖУХОВ Юрий Владимирович — кандидат технических наук заведующий кафедрой Санкт-Петербургского политехнического университета Петра Великого. 195251, Россия, г. Санкт-Петербург, Политехническая ул., 29. E-mail: [email protected]

KOZHUKHOV Yurii V. - Peter the Great St. Petersburg Polytechnic University. 29 Politechnicheskaya St., St. Petersburg, 195251, Russia. E-mail: [email protected]

ЯБЛОКОВ Алексей Михайлович — ассистент Санкт-Петербургского политехнического университета Петра Великого.

195251, Россия, г. Санкт-Петербург, Политехническая ул., 29. E-mail: [email protected]

YABLOKOV Aleksei M. — Peter the Great St. Petersburg Polytechnic University. 29 Politechnicheskaya St., St. Petersburg, 195251, Russia. E-mail: [email protected]

ЛЕБЕДЕВ Александр Анатольевич — кандидат технических наук доцент Санкт-Петербургского

политехнического университета Петра Великого.

195251, Россия, г. Санкт-Петербург, Политехническая ул., 29.

E-mail: [email protected]

LEBEDEV Aleksandr A. — Peter the Great St. Petersburg Polytechnic University. 29 Politechnicheskaya St., St. Petersburg, 195251, Russia. E-mail: [email protected]

© Санкт-Петербургский политехнический университет Петра Великого, 2016

i Надоели баннеры? Вы всегда можете отключить рекламу.