Теплоэнергетика, 2021, № 5, стр. 72-77

Определение истинного объемного паросодержания при гидравлическом расчете трубопроводов пароводяной геотермальной смеси

А. Н. Шулюпин^a, *, Н. Н. Варламова^a

^a Институт горного дела Дальневосточного отделения РАН – обособленное подразделение Хабаровского федерального исследовательского центра Дальневосточного отделения РАН
680000 г. Хабаровск, ул. Тургенева, д. 51, Россия

^* E-mail: ans714@mail.ru

Поступила в редакцию 27.07.2020
После доработки 31.08.2020
Принята к публикации 23.09.2020

DOI: 10.1134/S0040363621050106

Полный текст (PDF)

Аннотация

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

Ключевые слова: геотермальный теплоноситель, пароводяная смесь, трубопровод, истинное объемное паросодержание, модель дрейфа, плотность смеси, параметр распределения, скорость дрейфа

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

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

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

Трасса трубопровода пароводяной смеси от скважины Гео-1 до Верхне-Мутновской ГеоЭС. 1 – Верхне-Мутновская ГеоЭС; 2 – скважина Гео-1; 3 – трубопровод

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

ОПРЕДЕЛЕНИЕ ИСТИННОГО ОБЪЕМНОГО ПАРОСОДЕРЖАНИЯ

Практические задачи, для которых требуется установить истинное объемное паросодержание, решаются с помощью эмпирических формул. Известно множество эмпирических формул для вычисления непосредственно истинного объемного паросодержания (применительно к различным двухфазным средам – газосодержания) или величин, способствующих его дальнейшему теоретическому определению, таких как усредненная по сечению скорость пара или скольжение фаз [4–7]. Вид формул и рекомендуемые значения коэффициентов существенно зависят от условий экспериментов, результаты которых используются для обоснования формул. При этом отсутствуют экспериментальные данные по истинному объемному паросодержанию для трубопроводов пароводяного геотермального теплоносителя. Необходимые данные имеются только для труб диаметром до 0.2 м [7], в таких экспериментах большим считается диаметр 0.15 м [8], в то время как трубопроводы пароводяной смеси на геотермальных месторождениях имеют диаметр от 0.3 до 0.6 м.

Гравитационную составляющую перепада давления (Δp)_grav по направлению потока можно найти по выражению

${{(\Delta p)}_{{gra{v}}}} = - \int\limits_L {\rho g\sin \theta {\text{d}}L} ,$

где g – ускорение свободного падения; θ – угол наклона оси трубы относительно горизонтальной плоскости; L – длина трубы;

$\rho = {{\rho }_{g}}\alpha + {{\rho }_{l}}\left( {1 - \alpha } \right)$

– плотность смеси; ${{\rho }_{g}}$ и ${{\rho }_{l}}$ – плотность пара и воды; α – истинное объемное газосодержание (здесь паросодержание – отношение объема пара к общему объему рассматриваемого элемента).

Данные по трубопроводу от скважины Гео-1 (Верхне-Мутновская ГеоЭС)

Параметр	Опытные и расчетные данные
Параметр	на 16.09.2011 г.	на 11.09.2019 г.
Устьевое давление, МПа	1.130	0.890
Расход смеси, кг/с	65.0	65.3
Энтальпия смеси, кДж/кг	1221	1121
Перепад давления, МПа:
опытный	0.150	0.180
расчетный	0.152	0.176

Ввиду отсутствия экспериментальных данных и обоснованных эмпирических формул для расчета истинного объемного паросодержания в условиях транспортировки геотермального теплоносителя, важно выбрать удобный подход к его определению и провести анализ главных факторов, влияющих на конкретный вид рекомендуемых формул. Существует несколько подходов к определению истинного объемного газосодержания, в частности паросодержания. Наиболее теоретически корректным считается подход на основе модели дрейфа [4–7], согласно которому усредненная по сечению скорость газовой фазы ${{{v}}_{g}}$ вычисляется по формуле

(1)

${{{v}}_{g}} = {{C}_{0}}w + {{{v}}_{d}},$

где C₀ – параметр распределения; w – приведенная скорость смеси, равная скорости смеси по гомогенной модели; ${{{v}}_{d}}$ – скорость дрейфа.

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

(2)

$\alpha = {{{{w}_{g}}} \mathord{\left/ {\vphantom {{{{w}_{g}}} {{{{v}}_{g}}}}} \right. \kern-0em} {{{{v}}_{g}}}},$

где ${{w}_{g}}$ – приведенная скорость пара, представляющая собой отношение объемного расхода газа к общей площади сечения трубы.

ДЕТАЛИЗАЦИЯ МОДЕЛИ ДРЕЙФА

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

(3)

${{C}_{0}} = 1 + {{k}_{1}}\sin {{\theta }} + {{k}_{2}}{\text{cos}}\,{{\theta }} + {{k}_{3}},$

где k₁, k₂ и k₃ – некие коэффициенты.

На скорость дрейфа влияют два фактора:

гравитация, стремящаяся придать жидкой фазе дополнительную скорость относительно газа по направлению вниз;

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

Если считать второй фактор пропорциональным приведенной скорости смеси, то скорость дрейфа будет равна

(4)

${{{v}}_{d}} = {{k}_{4}}w + {{{v}}_{s}}\sin {{\theta ,}}$

где k₄ – некий коэффициент; v_s – скорость дрейфа, обусловленная гравитацией, для вертикального потока.

Подставляя (3) и (4) в (1), можно получить

(5)

${{{v}}_{g}} = Cw + {{{v}}_{s}}\sin {{\theta ,}}$

где C – параметр структуры (расширенный вследствие гидродинамической части скорости дрейфа параметр распределения), вычисляемый по формуле

(6)

$C = 1 + {{k}_{1}}\sin {{\theta }} + {{k}_{2}}{\text{cos}}\,{{\theta }} + {{k}_{3}} + {{k}_{4}}.$

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

${{k}_{1}} = {{k}_{2}} = {{k}_{3}} + {{k}_{4}}.$

Тогда получается

(7)

$C = 1 + {{k}_{1}}\left( {1 + \sin {{\theta }} + {\text{cos}}\,{{\theta }}} \right).$

Формула (5) предполагает отрицательное значение скорости дрейфа в нисходящих потоках. Таким образом обнаруживается недостаток традиционной модели дрейфа – возможность получения нереальных значений для скорости пара при малой приведенной скорости смеси. Истинная скорость пара не может принимать значения ниже его приведенной скорости, а равенство этих скоростей возможно только при нулевом расходе воды, а также при v_l → ∞. Чтобы исключить такую возможность, для нисходящих потоков предлагается вместо (5) использовать аналогичную формулу, но определяющую усредненную по сечению скорость воды:

(8)

${{{v}}_{l}} = Kw - {{{v}}_{s}}\sin {{\theta ,}}$

где

(9)

$K = 1 + {{k}_{5}}\left( {1 + \sin {{\theta }} + {\text{cos}}\,{{\theta }}} \right)$

– аналог параметра структуры; k₅ – коэффициент структуры нисходящего потока, при этом k₁ приобретает смысл обобщенного коэффициента структуры восходящего потока.

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

(10)

$\alpha = 1 - {{{{w}_{l}}} \mathord{\left/ {\vphantom {{{{w}_{l}}} {{{{v}}_{l}}}}} \right. \kern-0em} {{{{v}}_{l}}}},$

где ${{w}_{l}}$ – приведенная скорость воды, определяемая по ее расходу для полного сечения канала без учета наличия в нем пара.

Расчет коэффициента структуры нисходящего потока осуществляли таким образом, чтобы для горизонтального потока подходы, основанные на формулах (5) и (8), приводили к одинаковому результату. Используя (8) и (9) при условии θ = 0, легко убедиться, что

(11)

${{k}_{5}} = \frac{1}{2}\left( {\frac{{{{{v}}_{{l0}}}}}{w} - 1} \right),$

где

(12)

${{{v}}_{{l0}}} = \frac{{{{w}_{l}}}}{{1 - {{\alpha }_{0}}}}$

– скорость воды для горизонтального течения; ${{\alpha }_{0}}$ – истинное объемное паросодержание для горизонтального течения, которое можно найти по формуле (2) (${{\alpha }_{0}}$ = w_g/v_g0).

С учетом скорости пара согласно (5) и (7)

(13)

${{{v}}_{{g0}}} = w\left( {1 + 2{{k}_{1}}} \right),$

где ${{{v}}_{{g0}}}$ – скорость пара для горизонтального течения.

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

(14)

${{{v}}_{s}} = a{{\left[ {\frac{{g\sigma \left( {{{\rho }_{l}} - {{\rho }_{g}}} \right)}}{{\rho _{l}^{2}}}} \right]}^{{0.25}}},$

где a – коэффициент для скорости дрейфа; σ – коэффициент поверхностного натяжения.

Согласно литературным данным в различных формулах для определения значения a варьируют от 1.18 до 2.9 [4–7]. Следует заметить, что обычно (14) применяется для потоков с малым газосодержанием, когда дрейф вызван процессом всплытия газа в жидкости. Рассматриваемый случай больше соответствует движению жидкости в газовой среде, поэтому скорость дрейфа должна зависеть от плотности пара в большей степени, чем предписывает соотношение (14). Например, для падения капель используется аналогичная формула, но в знаменателе дроби стоит не плотность воды, а плотность пара [9], соответственно, скорость дрейфа характеризуется бо́льшими значениями. Поэтому в формуле (14) для рассматриваемого случая предлагается взять значение близкое к максимальному, но лежащее внутри указанного диапазона: a = 2.8.

Результаты опытов по критическому истечению пароводяных смесей в условиях геотермальных месторождений, в том числе значения динамического давления, указывают на наличие полной гомогенности смеси [10]. Для гомогенного потока C = 1, v_s = 0. При скоростях, приближающихся к критическим, эти величины должны стремиться к таким же значениям.

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

(15)

${\text{M}} = {{w}_{g}}{{\left( {\frac{{1.1xp}}{{{{\rho }_{g}}}}} \right)}^{{ - 0.5}}},$

где х – массовое расходное паросодержание; p – давление.

Учитывая в простейшем виде указанные факторы, для скорости дрейфа, обусловленной гравитацией, окончательно получим

(16)

${{{v}}_{s}} = 2.8\left( {1 - {\text{M}}} \right){{\left[ {\frac{{g\sigma \left( {{{\rho }_{l}} - {{\rho }_{g}}} \right)}}{{\rho _{l}^{2}}}} \right]}^{{0.25}}}.$

Часто в модели дрейфа для параметра распределения используют значение 1.2 [4–7], рекомендуемое для снарядного режима течения [9]. В рассматриваемом случае паросодержание обычно превышает значения, характерные для снарядного течения. В пароводяных скважинах, находящихся в условиях, близких к описываемым, для лучшего согласования с экспериментальными данными рекомендовано значение 1.1 [11]. Если ориентироваться в (7) на значение 0.05 для ${{k}_{1}},$ то применительно к вертикальному восходящему потоку (5) будет соответствовать обычно используемому виду классической модели дрейфа с рекомендованным для скважины коэффициентом распределения.

Среди величин, влияющих на значение параметра распределения, следует также отметить массовое расходное паросодержание. В рассматриваемых условиях стремление к единице расходного паросодержания должно приводить к такому же стремлению параметра распределения. Также параметр распределения должен зависеть от числа Маха. Как отмечалось, в критическом потоке смесь соответствует гомогенной модели (C₀ = 1). Гомогенная модель также будет характерна при близких плотностях фаз.

С учетом всех отмеченных факторов для обобщенного коэффициента структуры восходящего потока в (7) получим

(17)

${{k}_{1}} = 0.05\left( {1 - x} \right)\left( {1 - {\text{M}}} \right)\left( {1 - \frac{{{{\rho }_{g}}}}{{{{\rho }_{l}}}}} \right).$

Таким образом, для определения истинного объемного паросодержания в трубопроводах, транспортирующих пароводяной геотермальный теплоноситель, в восходящих потоках рекомендуется применять формулы (2), (5), (7) и (15)–(17), а в нисходящих – (8)–(13) и (15)–(17). Данные рекомендации были использованы при создании программы по гидравлическому расчету трубопроводов пароводяной смеси на геотермальных месторождениях. В программе наряду с гравитационной составляющей перепада давления проводится расчет составляющих перепада давления на трение и местные сопротивления по формулам, которые дают результаты, хорошо согласующиеся с опытными данными [12].

СРАВНЕНИЕ РАСЧЕТНЫХ И ЭСПЕРИМЕНТАЛЬНЫХ ДАННЫХ

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

Для верификации программы в части учета гравитационной составляющей было проведено сравнение расчетных и опытных данных по перепаду давления в трубопроводе от скважины Гео-1 (Верхне-Мутновская ГеоЭС), имеющем максимальный перепад высот из всех действующих отечественных трубопроводов пароводяной смеси, для двух моментов времени: 16.09.2011 и 11.09.2019. Схема трассы трубопровода представлена на рисунке.

Трубопровод имеет длину 1050 м, внутренний диаметр 0.406 м, суммарный коэффициент местных сопротивлений 8, транспортировка осуществляется вниз по рельефу с перепадом высот 110 м. Данный трубопровод характеризуется незначительным отложением солей в процессе эксплуатации, что позволяет принять единое значение коэффициента шероховатости (0.2 мм) для обоих моментов времени.

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

Погрешность опытного определения перепада давления оценивается как ±0.2 × 10⁵ Па. Расхождение расчетного и опытного перепадов не выходит за указанный диапазон, что косвенно подтверждает адекватность рекомендованной процедуры определения истинного объемного паросодержания.

ВЫВОДЫ

1. Истинное объемное паросодержание при гидравлическом расчете трубопроводов пароводяной геотермальной смеси рекомендуется определять на основе модели дрейфа. При этом для усредненной скорости пара целесообразно использовать формулу (5), включающую в себя расширенный благодаря гидродинамической части скорости дрейфа параметр распределения и обусловленную исключительно гравитацией скорость дрейфа.

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

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

Список литературы

Lee K.C., Jenks D.G. Ohaaki geothermal steam transmission pipelines // Proc. of the 11th New Zealand Geothermal Workshop. 1989. P. 25–30. https://www.geothermal-energy.org/pdf/IGAstandard/ NZGW/1989/Lee.pdf
Rizaldy, Zarrouk S.J. Pressure drop in large diameter geothermal two-phase pipelines // Proc. of the 38th New Zealand Geothermal Workshop. 23–25 Nov. 2016. Auckland, New Zealand, 2016. P. 1–5. https://www.geothermal-energy.org/pdf/IGAstandard/NZGW/2016/ 065_Rizaldy_Final.pdf
Шулюпин А.Н., Чермошенцева А.А., Варламова Н.Н. Новые вызовы при освоении месторождений парогидротерм с транспортировкой пароводяной смеси // Горный информационно-аналитический бюллетень. 2019. № 2. С. 43–49. https://doi.org/10.25018/0236-1493-2019-02-0-43-49
Woldesemayat M.A., Ghajar A.J. Comparison of void fraction correlations for different flow patterns in horizontal and upward inclined pipes // Int. J. Multiphase Flow. 2007. V. 33. Is. 4. P. 347–370. https://doi.org/10.1016/j.ijmultiphaseflow.2006.09.004
Bhagwat S.M., Ghajar A.J. Similarities and differences in the flow patterns and void fraction in vertical upward and downward two phase flow // Exp. Therm. Fluid Sci. 2012. V. 39. P. 213–227. https://doi.org/10.1016/j.expthermflusci.2012.01.026
Xu Y., Fang X. Correlations of void fraction for two-phase refrigerant flow in pipes // Appl. Therm. Eng. 2014. V. 64. Is. 1–2. P. 242–251. https://doi.org/10.1016/j.applthermaleng.2013.12.032
Bhagwat S.M., Ghajar A.J. A flow pattern independent drift flux model based void fraction correlation for a wide range of gas-liquid two phase flow // Int. J. Multiphase Flow. 2014. V. 59. P. 186–205. https://doi.org/10.1016/j.ijmultiphaseflow.2013.11.001
Experimental study on void fraction, pressure drop and flow regime analysis in a large ID piping system / Z. Dang, Z. Yang, X. Yang, M. Ishii // Int. J. Multiphase Flow. 2019. V. 111. P. 31–41. https://doi.org/10.1016/j.ijmultiphaseflow.2018.10.006
Теоретические основы теплотехники. Теплотехнический эксперимент: справочник / под ред. Григорьева В.А. и Зорина В.М. М.: Энергоатомиздат, 1988.
Шулюпин А.Н. Вопросы гидравлики пароводяной смеси при освоении геотермальных месторождений. Владивосток: Дальнаука, 2011.
Дрознин В.А. Физическая модель вулканического процесса. М.: Наука, 1980.
Shulyupin A.N., Chermoshentseva A.A., Varlamova N.N. Numerical study of the stability of the steam-water flow in pipelines of geothermal gathering system // CEUR Workshop Proc. (Information Technologies and High-Performance Computing). Khabarovsk, Russia, 16–19 Sept. 2019. P. 103–109. http://ceur-ws.org/Vol-2426/paper15.pdf

Дополнительные материалы отсутствуют.

Инструменты

следующая статья выпуска предыдущая статья выпуска содержание выпуска

Теплоэнергетика

Архивы выпусков Информация о журнале Отправить рукопись в журнал